CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                  C
C     This program is Monte Carlo simulation of dynamically        C
C     triangularized random surfaces.                              C
C     Writen by Enzo Marinari, June 1992                           C
C     Optimized by Leping Han, Aug. 1992                           C
C     Syracuse University, Aug. 1992                               C
C                                                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      PROGRAM pippo

      INCLUDE 'string.h'
      
      CALL init

      IF(start_mode.ne.-1)CALL init_field

      DO i_sweep=1, n_discard
        CALL flip_many_links
        CALL update_many_sites
      ENDDO

      DO i_sweep = 1, n_sweeps
        DO i_correlated = 1, n_correlated
          DO i_one_sweep =1, time_scale
            CALL flip_many_links
            CALL update_many_sites
          ENDDO
          CALL measure
          CALL average_1
        ENDDO
        CALL average_2
      ENDDO

      CALL check_mesh

      CALL ciao

      STOP
  123 FORMAT(a1)
      END

C______________________________________________________________________

      SUBROUTINE ciao

      INCLUDE 'string.h'
      INTEGER id, is, im, io, i_two
      REAL * 8 diff

      PRINT *,'Flip acceptance rate _ ',tot_lacc
      PRINT *,'Site acceptance rate _ ',tot_sacc

      WRITE(f_out,*)'Flip acceptance rate _ ',tot_lacc
      WRITE(f_out,*)'Site acceptance rate _ ',tot_sacc
      
      DO io = 1, n_obs
        DO i_two=1,2
          diff =  (tot_obs(i_two,2,io)-tot_obs(i_two,1,io)**2)
     $          * inv_n_sweeps_1
          IF(diff.GT.0)THEN
            tot_obs(i_two,3,io) = sqrt(diff)
          ELSE
            tot_obs(i_two,3,io) = -999999.
          ENDIF
        ENDDO
      ENDDO

      DO io = 1, n_obs
        DO i_two=1,2
          PRINT*,i_two,phrase(io),tot_obs(i_two,1,io),' +/- '
     $            ,tot_obs(i_two,3,io)
          WRITE(f_out,*)i_two,phrase(io),tot_obs(i_two,1,io),' +/- '
     $            ,tot_obs(i_two,3,io)
        ENDDO
      ENDDO

      WRITE(f_oco,*)nx,ny,n,i_seed1,i_seed2
     $     ,time_scale,n_sweeps,n_correlated
     $     ,n_discard,lambda,trial,write_all
        write(f_oco,387)(rand_u(id),id=1,97)
        write(f_oco,387)rand_c, rand_cd, rand_cm
        write(f_oco,*)rand_i97, rand_j97, rand_ivec, rand_test
        write(f_oco,387)((field           (id,is),id=0,dm),is=0, nm)
        write(f_oco,387)(square_m        (   im)        ,im=0,n6m)
        write(f_oco,387)((lno_m           (id,im),id=0,dm),im=0,n6m)
        write(f_oco,387)((rno_m           (id,im),id=0,dm),im=0,n6m)
        write(f_oco,*)(m_start_s       (   is)        ,is=0, nm)
        write(f_oco,*)(n_connectivity_s(   is)        ,is=0, nm)
        write(f_oco,*)(edge_m          (   im)        ,im=0,n6m)
        write(f_oco,*)(s_site_m        (   im)        ,im=0,n6m)
        write(f_oco,*)(s_neighbor_m    (   im)        ,im=0,n6m)
        write(f_oco,*)(m_next_m        (   im)        ,im=0,n6m)
        write(f_oco,*)(m_former_m      (   im)        ,im=0,n6m)
        write(f_oco,*)(m_back_m        (   im)        ,im=0,n6m)
 387    format(1x,4g16.7)
      IF(write_all .EQ. 1)THEN
         CLOSE(f_all)
         CLOSE(f_all2)
      ENDIF
      CLOSE(f_oco)
      CLOSE(f_out)

      RETURN
      END

C_______________________________________________________________________

      SUBROUTINE init

      INCLUDE 'string.h'

      INTEGER i, io, i_two, i_three, m0
      INTEGER id, is, im
      INTEGER rnx, rny, rn, ri_seed1, ri_seed2, rtime_scale, rn_sweeps
      INTEGER rn_correlated, rn_discard,rwrite_all
      REAL*4 a, rlambda, rtrial
      real ranmar

CC    please, keep the phrases 'one word' ('xysfsdf', no blanks)
      phrase( 1) = '|q-6|     '
      phrase( 2) = 'log(q)    '
      phrase( 3) = 'GaussS    '
      phrase( 4) = 'q_max     '
      phrase( 5) = 'Ra        '
      phrase( 6) = 'Se        '
      phrase( 7) = 'I(|K|)    '
      phrase( 8) = 'Area      '
      phrase( 9) = 'I(|K-R|)  '
      phrase(10) = 'I(R**2)   '
      phrase(11) = 'I(K**2)   '
      phrase(12) = 'I(K*R)    '

      OPEN(unit=f_out,status='new',file='aaa.out')

      WRITE(f_out,*)' nx = ', nx, ' ny = ', ny, ' n = ',n
      PRINT*,' nx = ', nx, ' ny = ', ny, ' n = ',n
      PRINT *, 'Input first i_seed (0 <= seed <= 31328) _ '
      PRINT *, 'NOT USED IF START_MODE = -1'
      READ(5,*)i_seed1
      PRINT *,' i_seed1 = ', i_seed1
      WRITE(f_out,*)' i_seed1 = ', i_seed1
      PRINT *, 'Input second i_seed (0 <= seed <= 30081) _ '
      PRINT *, 'NOT USED IF START_MODE = -1'
      READ(5,*)i_seed2
      PRINT *,' i_seed2 = ', i_seed2
      WRITE(f_out,*)' i_seed2 = ', i_seed2
      PRINT *, 'Input time_scale (how many multiples of 3n per measure'
      READ(5,*)time_scale
      PRINT *,' time_scale = ', time_scale
      WRITE(f_out,*)' time_scale = ', time_scale
      PRINT *, 'Input n_sweeps (independent measurements)_ '
      READ(5,*)n_sweeps
      PRINT *,' n_sweeps = ', n_sweeps
      WRITE(f_out,*)' n_sweeps = ', n_sweeps
      PRINT *, 'Input n_correlated (internal loop)_ '
      READ(5,*)n_correlated
      PRINT *,' n_correlated = ', n_correlated
      WRITE(f_out,*)' n_correlated = ', n_correlated
      PRINT *, 'Input n_discard (internal loop)_ '
      READ(5,*)n_discard
      PRINT *,' n_discard = ', n_discard
      WRITE(f_out,*)' n_discard = ', n_discard
      PRINT *, 'Input lambda_ '
      READ(5,*)lambda
      PRINT *,' lambda = ', lambda
      WRITE(f_out,*)' lambda = ', lambda
      PRINT *, 'Input trial_ '
      READ(5,*)trial
      PRINT *,' trial = ', trial
      WRITE(f_out,*)' trial = ', trial
      PRINT *,'Input 1 if write all _'
      READ(5,*)write_all
      PRINT *,'write_all = ', write_all
      WRITE(f_out,*)'write_all = ', write_all
      PRINT *,'Input start_mode (-1 for continuing from aaa.ico)'
      READ(5,*)start_mode
      PRINT *,'start_mode = ', start_mode
      WRITE(f_out,*)'start_mode = ', start_mode

CC    calculate all inverse number 

      inv_n = 1.0/float(n)
      inv_3n = 1.0/float(3*n)
      inv_n6 = 1.0/float(n6)
      inv_n_sweeps = 1.0/float(n_sweeps)
      inv_n_sweeps_1 = 1.0/float(n_sweeps-1)
      inv_n_correlated = 1.0/float(n_correlated)
      third_two_pi = two_pi/3.0

      IF(write_all .EQ. 1)THEN
        OPEN(unit=f_all ,status='new',file='aaa.all')
        OPEN(unit=f_all2,status='new',file='aaa.all2')
        write(f_all,*)n,time_scale,n_correlated, n_sweeps
        write(f_all2,*)n,time_scale,n_correlated, n_sweeps
      ENDIF
      IF(start_mode .EQ. -1)THEN
        OPEN(unit=f_ico,status='old',file='aaa.ico')
        read(f_ico,*)rnx,rny,rn,ri_seed1,ri_seed2,rtime_scale,rn_sweeps
     $       ,rn_correlated,rn_discard,rlambda,rtrial,rwrite_all
        if(rnx.ne.nx)then
           print*,rnx,nx
           stop'rnx .ne. nx'
        endif
        if(rny.ne.ny)then
           print*,rny,ny
           stop'rny .ne. ny'
        endif
        if(rn.ne.n)then
           print*,rn,n
           stop'rn .ne. n'
        endif
        read(f_ico,387)(rand_u(id),id=1,97)
        read(f_ico,387)rand_c, rand_cd, rand_cm
        read(f_ico,*)rand_i97, rand_j97, rand_ivec, rand_test
        read(f_ico,387)((field           (id,is),id=0,dm),is=0, nm)
        read(f_ico,387)(square_m        (   im)        ,im=0,n6m)
        read(f_ico,387)((lno_m           (id,im),id=0,dm),im=0,n6m)
        read(f_ico,387)((rno_m           (id,im),id=0,dm),im=0,n6m)
        read(f_ico,*)(m_start_s       (   is)        ,is=0, nm)
        read(f_ico,*)(n_connectivity_s(   is)        ,is=0, nm)
        read(f_ico,*)(edge_m          (   im)        ,im=0,n6m)
        read(f_ico,*)(s_site_m        (   im)        ,im=0,n6m)
        read(f_ico,*)(s_neighbor_m    (   im)        ,im=0,n6m)
        read(f_ico,*)(m_next_m        (   im)        ,im=0,n6m)
        read(f_ico,*)(m_former_m      (   im)        ,im=0,n6m)
        read(f_ico,*)(m_back_m        (   im)        ,im=0,n6m)
 387    format(1x,4g16.7)
        close(f_ico)
      ELSE
        CALL rmarin(i_seed1,i_seed2)
      ENDIF
      OPEN(unit=f_oco,status='new',file='aaa.oco')
      scale = 1.0
      md_half = - float(dim) / 2.0
      anx2 = float(nx) / 2.0 - .5
      any2 = float(ny) / 2.0 - .5

      if(start_mode .ne. -1)then
         DO i=1,200
            a = ranmar()
         ENDDO
         CALL set_up_mesh
      endif

      DO m0=1,n6m
        tlog(m0)=md_half*log(float(m0))
      ENDDO
      tlog(0)=-1.e28

      ind_lacc  = 0.0
      ind_sacc  = 0.0
      DO io=1,n_obs
        DO i_two =1,2
          ind_obs(i_two,io)=0.0
          DO i_three=1,3
            tot_obs(i_two,i_three,io)=0.0
          ENDDO
        ENDDO
      ENDDO

      tot_lacc   = 0.0
      tot_sacc   = 0.0

      RETURN
      END

C_______________________________________________________________________

      SUBROUTINE init_field

      INCLUDE 'string.h'

      INTEGER d0, s0, m0, s1
      INTEGER random_start, old_conf
      INTEGER flat_start_xy, flat_start_xz, flat_start_yz
      INTEGER bulk_start_xy, bulk_start_xz, bulk_start_yz

      REAL scratch
      real ranmar

      old_conf = -1
      random_start = 0
      flat_start_xy = 1
      flat_start_xz = 2
      flat_start_yz = 3
      bulk_start_xy = 4
      bulk_start_xz = 5
      bulk_start_yz = 6

      IF(start_mode .EQ. random_start)THEN
        DO d0=0,dm
          DO s0=0,nm
            field(d0,s0) = scale * (.5-ranmar())
          ENDDO
        ENDDO
      ELSEIF(start_mode .EQ. flat_start_xy)THEN
        DO s0=0,nm
          field(0,s0) = scale * (.5-ranmar())
          field(1,s0) = scale * (.5-ranmar())
          field(2,s0) = 0.0
        ENDDO
      ELSEIF(start_mode .EQ. flat_start_xz)THEN
        DO s0=0,nm
          field(0,s0) = scale * (.5-ranmar())
          field(1,s0) = 0.0
          field(2,s0) = scale * (.5-ranmar())
        ENDDO
      ELSEIF(start_mode .EQ. flat_start_yz)THEN
        DO s0=0,nm
          field(0,s0) = 0.0
          field(1,s0) = scale * (.5-ranmar())
          field(2,s0) = scale * (.5-ranmar())
        ENDDO
      ELSEIF(start_mode .EQ. bulk_start_xy)THEN
        DO s0=0,nm
          field(0,s0) = float(mod(s0,nx))
          field(1,s0) = float((s0-mod(s0,nx))/nx)
          field(2,s0) = 0.0
        ENDDO
      ELSEIF(start_mode .EQ. bulk_start_xz)THEN
        DO s0=0,nm
          field(0,s0) = float(mod(s0,nx))
          field(1,s0) = 0.0
          field(2,s0) = float((s0-mod(s0,nx))/nx)
        ENDDO
      ELSEIF(start_mode .EQ. bulk_start_yz)THEN
        DO s0=0,nm
          field(0,s0) = 0.0
          field(1,s0) = float(mod(s0,nx))
          field(2,s0) = float((s0-mod(s0,nx))/nx)
        ENDDO
      ELSE
         stop'which start_mode????'
      ENDIF

      DO m0 = 0, n6m
        square_m(m0) = 0.0
        s0 = s_site_m(m0)
        s1 = s_neighbor_m(m0)
        DO d0 = 0, dm
          scratch = field(d0,s0)-field(d0,s1)
          square_m(m0) = square_m(m0) + scratch * scratch
        ENDDO
      ENDDO

      CALL init_curvature

      RETURN
  123 FORMAT(a1)
      END

C_______________________________________________________________________

      SUBROUTINE set_up_mesh

      INCLUDE 'string.h'

      INTEGER site, site1, s0, s2, m_start, m_start2
      integer nc,nc1,i,j
      DO site = 0 , nm
        m_start_s(site) = 6 * site
        n_connectivity_s(site) = 6
      ENDDO

      DO site = 0 , nm
        s_neighbor_m(6*site+0) = site + 1
        s_neighbor_m(6*site+1) = site + nx + 1
        s_neighbor_m(6*site+2) = site + nx
        s_neighbor_m(6*site+3) = site - 1
        s_neighbor_m(6*site+4) = site - nx - 1
        s_neighbor_m(6*site+5) = site - nx

        DO site1 = 0, 5
          m_next_m  (6*site+site1) = 6*site + site1 + 1
          m_former_m(6*site+site1) = 6*site + site1 - 1
          s_site_m  (6*site+site1) = site
        ENDDO
        m_next_m   ( 6*site + 5 ) = 6*site
        m_former_m ( 6*site     ) = 6*site + 5

        IF(site .LT. nx)THEN
          s_neighbor_m(6*site+5) = site + nx*(ny-1)
          s_neighbor_m(6*site+4) = site + nx*(ny-1) - 1
        ENDIF
        IF(nx*ny-1-site .lt. nx)then
          s_neighbor_m(6*site+2) = site - nx*(ny-1)
          s_neighbor_m(6*site+1) = site - nx*(ny-1) + 1
        ENDIF
        IF(mod(site,nx) .eq. nx-1)then
          s_neighbor_m(6*site+0) = site - nx + 1
          s_neighbor_m(6*site+1) = site + 1
        ENDIF
        IF(mod(site,nx) .eq. 0)then
          s_neighbor_m(6*site+3) = site + nx - 1
          s_neighbor_m(6*site+4) = site - 1
        ENDIF
      ENDDO

      s_neighbor_m(0+4) = nx * ny - 1
      s_neighbor_m(6*(nx*ny-1)+1) = 0

      DO s0=0, nm
        m_start = m_start_s(s0)
        nc = n_connectivity_s(s0)
        do i=1,nc
           s2 = s_neighbor_m(m_start)
           m_start2 = m_start_s(s2)
           nc1 = n_connectivity_s(s2)
           do j=1,nc1
              IF(s0.EQ.s_neighbor_m(m_start2))GOTO 52
              m_start2 = m_next_m(m_start2)
           enddo
 52        CONTINUE

CC    here m_start and m_start2 are the two opposite links.
           m_back_m(m_start ) = m_start2
           m_back_m(m_start2) = m_start
           m_start = m_next_m(m_start)
        enddo
      ENDDO

      RETURN
      END
C_______________________________________________________________________

      SUBROUTINE average_1

      INCLUDE 'string.h'
      INTEGER io

      ind_lacc   = ind_lacc  + sw_lacc * inv_n_correlated
      ind_sacc   = ind_sacc  + sw_sacc * inv_n_correlated

      DO io=1,n_obs
        ind_obs(1,io) = ind_obs(1,io)
     $        +sw_obs(io) * inv_n_correlated
        ind_obs(2,io) = ind_obs(2,io)
     $        +sw_obs(io)*sw_obs(io) * inv_n_correlated
      ENDDO

      RETURN
      END

C_______________________________________________________________________

      SUBROUTINE average_2

      INCLUDE 'string.h'
      INTEGER i, io

      PRINT *,i_sweep
      PRINT *,'lacc,sacc ',ind_lacc, ind_sacc
      DO io=1,n_obs
        PRINT *, phrase(io),ind_obs(1,io), ind_obs(2,io)
      ENDDO

      tot_lacc  = tot_lacc + ind_lacc * inv_n_sweeps
      tot_sacc  = tot_sacc + ind_sacc * inv_n_sweeps

      DO io=1,n_obs
        DO i=1,2
          tot_obs(i,1,io) = tot_obs(i,1,io)
     $           + ind_obs(i,io) * inv_n_sweeps
          tot_obs(i,2,io) = tot_obs(i,2,io)
     $           + ind_obs(i,io)*ind_obs(i,io) * inv_n_sweeps
        ENDDO
      ENDDO

      ind_lacc = 0.0
      ind_sacc = 0.0

      DO io=1,n_obs
        DO i=1,2
          ind_obs(i,io) = 0.0
        ENDDO
      ENDDO

      RETURN
      END

C_______________________________________________________________________

      SUBROUTINE init_curvature

      INCLUDE 'string.h'
      INTEGER m_s, m_s_l,m_s_r
      INTEGER i,nc,idim
      INTEGER s0, s1, s3, s4
      REAL dv(0:dm), f(0:dm), g(0:dm)
      REAL d31(0:dm), d41(0:dm), norm_do
      REAL dot_31, factor31, h_ms_r(0:dm)
      REAL norm_old_right, norm_old_left
      REAL dot_41, factor41, h_ms_l(0:dm)
      REAL inorm_31, inorm_41

      DO s0=0, nm
        m_s = m_start_s(s0)
        nc = n_connectivity_s(s0)
        do i=1,nc
           s1 = s_neighbor_m(m_s)
           m_s_r = m_former_m(m_s)
           s3 = s_neighbor_m(m_s_r)
           m_s_l = m_next_m(m_s)
           s4 = s_neighbor_m(m_s_l)
           dot_31 = 0
           dot_41 = 0
           do idim=0,dm
              f(idim) = field(idim,s0)
              g(idim) = field(idim,s1)
              dv(idim) = g(idim)  - f(idim)
              d31(idim) = field(idim,s3) - f(idim)
              d41(idim) = field(idim,s4) - f(idim)
              dot_31    = dot_31 + dv(idim)*d31(idim)
              dot_41    = dot_41 + dv(idim)*d41(idim)
           enddo

           norm_do = square_m(m_s)
           factor31 = dot_31/norm_do
           factor41 = dot_41/norm_do

           norm_old_right = 0
           norm_old_left = 0
           do idim=0,dm
              h_ms_r(idim) = factor31 * dv(idim) - d31(idim)
              h_ms_l(idim) = factor41 * dv(idim) - d41(idim)
              norm_old_right = norm_old_right + 
     $             h_ms_r(idim)*h_ms_r(idim)
              norm_old_left = norm_old_left + 
     $             h_ms_l(idim)*h_ms_l(idim)
           enddo

           inorm_31 = 1.0 / sqrt ( norm_old_right )
           inorm_41 = 1.0 / sqrt ( norm_old_left  )

           do idim=0,dm
              rno_m(idim,m_s) = h_ms_r(idim) * inorm_31
              lno_m(idim,m_s) = h_ms_l(idim) * inorm_41
           enddo

           edge_m(m_s) = rno_m(0,m_s) * lno_m(0,m_s)
     $          + rno_m(1,m_s) * lno_m(1,m_s) 
     $          + rno_m(2,m_s) * lno_m(2,m_s)

           m_s = m_next_m(m_s)

        enddo
      ENDDO

      RETURN
      END
C_______________________________________________________________________

      SUBROUTINE measure

      INCLUDE 'string.h'

      INTEGER io, s0, m0, m_start, nc, sn, sr,link
      REAL f0, f1, f2, qmax, arg_scratch, arg_scratch1
      REAL fn0, fn1, fn2, fr0, fr1, fr2
      REAL site_d_angle, scalar_p, square_r, sq1, sq2
      REAL area, site_d_area
      REAL bari(0:dm), scratch(0:dm)
      REAL tmp
      INTEGER i,tmp_link

CC    unroll DO loop here gives a little improvement
      DO io=1, n_obs
        sw_obs(io) = 0.0
      ENDDO

      qmax = 0.0
      bari(0) = 0.0
      bari(2) = 0.0
      bari(3) = 0.0

      DO s0 = 0, nm
        nc = n_connectivity_s(s0)
        IF(float(nc) .GT. qmax)qmax=float(nc)
        sw_obs(1)  = sw_obs(1)  + abs(nc-6)
        sw_obs(2) = sw_obs(2) + log(float(nc))
        bari(0) = bari(0) + field(0,s0)
        bari(1) = bari(1) + field(1,s0)
        bari(2) = bari(2) + field(2,s0)
      ENDDO


      bari(0) = bari(0) * inv_n
      bari(1) = bari(1) * inv_n
      bari(2) = bari(2) * inv_n

      DO s0 = 0, nm
         scratch(0) = field(0,s0)-bari(0)
         scratch(1) = field(1,s0)-bari(1)
         scratch(2) = field(2,s0)-bari(2)
         sw_obs(5) = sw_obs(5) + scratch(0)*scratch(0) +
     $        scratch(1)*scratch(1) + scratch(2)*scratch(2)
      ENDDO
      sw_obs(1) =  third_two_pi * sw_obs(1) * inv_n
      sw_obs(2) = sw_obs(2) * inv_n
      sw_obs(4) = qmax
      sw_obs(5) = sw_obs(5) * inv_n

      DO m0 = 0, n6m
        sw_obs(3) = sw_obs(3) + square_m(m0)
      ENDDO

      sw_obs(3) = 0.5 * sw_obs(3)

      do link=0,n6m
         sw_obs(6) = sw_obs(6) + edge_m(link) + 1.0
      enddo

      sw_obs(6) = 0.5 * sw_obs(6) * inv_n

      DO s0=0, nm
        nc = n_connectivity_s(s0)
        site_d_angle = 0.0
        site_d_area  = 0.0

        m_start = m_start_s(s0)
        do i=1,nc
           sn = s_neighbor_m(m_start)
           tmp_link = m_former_m(m_start)
           sr = s_neighbor_m(tmp_link)
           f0 = field(0,s0)
           f1 = field(1,s0)
           f2 = field(2,s0)
           fn0 = field(0,sn)
           fn1 = field(1,sn)
           fn2 = field(2,sn)
           fr0 = field(0,sr)
           fr1 = field(1,sr)
           fr2 = field(2,sr)

           scalar_p = (fn0-f0)*(fr0-f0)
     $          + (fn1-f1)*(fr1-f1)+(fn2-f2)*(fr2-f2)
           square_r = (fr0-f0)**2+(fr1-f1)**2+(fr2-f2)**2

           sq2 = square_r
           sq1 = square_m(m_start)

CC    scalar_p here is square of cos(\theta), arg_scratch
CC    is square of sin(\theta)
             
           arg_scratch1 = sq1 * sq2
           if(arg_scratch1 .gt. 0)then
              scalar_p = scalar_p*scalar_p / arg_scratch1
           else
              print *, 'sq1,sq2,sq1*sq2 ',sq1, sq2, arg_scratch1
              print *,'MEASURE SOS, SQ1*SQ2=0, SCALAR_P=0 SET'
              scalar_p = 0
           endif

           arg_scratch = 1.-scalar_p
           if(arg_scratch .ge. 0)then
              area = 0.5 * sqrt(arg_scratch1 * arg_scratch)
           else
              print *,'scalar_p, 1-sp**2 ', scalar_p, arg_scratch
              print *,f0,f1,f2
              print *,fn0,fn1,fn2
              print *,fr0,fr1,fr2
              print *,'MEASURE SOS, 1-sca**2 < 0, AREA=0 set'
              area = 0.0
           endif

CC    square of cos(2\theta)=2*(square of cos(\theta)) - 1

           scalar_p = 2*scalar_p - 1.0
           if(abs(scalar_p).le.1.0)then
              site_d_angle = site_d_angle + 0.5*acos(scalar_p)
           elseif(scalar_p .gt.0)then
              print *,'scalar_p',scalar_p
              print *,'MEASURE SOS, scalar_p > 1.0, angle++0 set'
              site_d_angle = site_d_angle + 0.0
           else
              print *,'scalar_p',scalar_p
              print *,'MEASURE SOS, scalar_p < 1.0, angle++pi set'
              site_d_angle = site_d_angle + two_pi * 0.25
           endif
           site_d_area  = site_d_area  + area

           m_start = m_next_m(m_start)
        enddo
        tmp = float(nc)
        sw_obs(7) = sw_obs(7) + abs(two_pi - site_d_angle)
        sw_obs(8) = sw_obs(8) + site_d_area
        sw_obs(9) = sw_obs(9) + abs(site_d_angle - 0.5*third_two_pi*tmp)
        sw_obs(10) = sw_obs(10) + tmp * 0.33333333
     $       * ( two_pi*two_pi * ((6-tmp)/tmp)**2)
        sw_obs(11) = sw_obs(11) + tmp * 3.
     $       * ( (two_pi-site_d_angle) / site_d_area ) ** 2
        sw_obs(12) = sw_obs(12) + tmp
     $       * ( two_pi * ((6.-tmp)/tmp)
     $       * ( two_pi - site_d_angle)/site_d_area)
      ENDDO
      sw_obs(12) = sw_obs(12)/sqrt(sw_obs(10)*sw_obs(11))
      sw_obs(7) = sw_obs(7) * inv_n
      sw_obs(8) = sw_obs(8) * inv_n6
      sw_obs(9) = sw_obs(9) * inv_n
      sw_obs(10) = sw_obs(10) * inv_n
      sw_obs(11) = sw_obs(11) * inv_n

      IF(write_all .EQ. 1)THEN
        WRITE(f_all ,181)(sw_obs(io),io=1,6)
        WRITE(f_all2,181)(sw_obs(io),io=7,n_obs)
      ENDIF

      RETURN
  123 FORMAT(a1)
  181 FORMAT(1x,6g13.6)
      END

C______________________________________________________________________

      SUBROUTINE flip_many_links

      INCLUDE 'string.h'

      INTEGER link, s1, s2, s3, s4, i
      INTEGER c1, c2, c3, c4
      INTEGER i_flip
      INTEGER m_start, site, ma ,mb 
      INTEGER m_f1, m_f2, m_f3, m_f4
      INTEGER m_b1, m_b2, m_b3, m_b4
      
      REAL len2_31, len2_23, len2_42, len2_14
      REAL dot_13, factor13
      REAL dot_23, factor23
      REAL factor51
      REAL factor63
      REAL dot_72,factor72
      REAL dot_84,factor84
      REAL e51t,e63t,e72t,e84t

      INTEGER idim
      REAL d31(0:dm),d41(0:dm),d42(0:dm),d23(0:dm),dn(0:dm)
      REAL h_new_r(0:dm),h_new_l(0:dm),h_13_l(0:dm)
      REAL h_23_l(0:dm),h_24_l(0:dm),h_14_r(0:dm)
      REAL norm_scratch23,norm_scratch13,norm_scratch51
      REAL norm_scratch63,norm_scratch72,norm_scratch84
      REAL edget1
      REAL len2_o, len2_n, edge, edget,inv_len2_n
      real ranmar

      sw_lacc = 0

      DO i_flip = 1, 3 * n

        link = n6 * (1.0-ranmar())

        IF(link.EQ.n6)THEN
          link = n6m
        ENDIF

        s1 = s_site_m(link)
        c1 = n_connectivity_s(s1)
        IF(c1 .LT. 4)GOTO 10
        s2 = s_neighbor_m(link)
        c2 = n_connectivity_s(s2)
        IF(c2 .LT. 4)GOTO 10

        m_f1 = m_former_m(link)
        m_b4 = m_next_m(link) 
        m_b1 = m_back_m(m_f1)
        m_f2 = m_former_m(m_b1)
        m_f4 = m_back_m(m_b4)
        m_b3 = m_next_m(m_f4)
        m_f3 = m_back_m(m_b3)
        m_b2 = m_back_m(m_f2)

        s3 = s_neighbor_m(m_f1)
        s4 = s_neighbor_m(m_b4)

        c3 = n_connectivity_s(s3)
        c4 = n_connectivity_s(s4)

        m_start = m_start_s(s3)
        do i=1,c3
           site = s_neighbor_m(m_start)
           IF(site .EQ. s4)GOTO 10
           m_start = m_next_m(m_start)
        enddo

CC    Here we know that topology is ok: we can accept.
CC    Now we find what the links are around the pologon

        len2_n = 0.0
        dot_23 = 0.0
        dot_13 = 0.0
        dot_72 = 0.0
        dot_84 = 0.0
        do idim=0,dm
           d31(idim) = field(idim,s3) - field(idim,s1)
           d41(idim) = field(idim,s4) - field(idim,s1)
           dn(idim) = field(idim,s4) - field(idim,s3)
           d23(idim) = field(idim,s2) - field(idim,s3)
           d42(idim) = field(idim,s4) - field(idim,s2)
           len2_n  = len2_n + dn(idim)*dn(idim)
           dot_23    = dot_23 + dn(idim)*d23(idim)
           dot_13    = dot_13 + dn(idim)*d31(idim)
           dot_72    = dot_72 + dn(idim)*d42(idim)
           dot_84    = dot_84 + dn(idim)*d41(idim)
        enddo

        len2_o  = square_m(link)
        len2_31  = square_m(m_f1)
        len2_23  = square_m(m_f2)
        len2_42  = square_m(m_f3)
        len2_14  = square_m(m_f4)

        inv_len2_n = 1.0/len2_n
        factor23 = dot_23 * inv_len2_n
        factor13 = - dot_13 * inv_len2_n
        factor51 =  dot_13/len2_31
        factor63 = dot_23/len2_23
        factor72 = - dot_72/len2_42
        factor84 = - dot_84/len2_14

        norm_scratch23 = 0.0
        norm_scratch13 = 0.0
        norm_scratch51 = 0.0
        norm_scratch63 = 0.0
        norm_scratch72 = 0.0
        norm_scratch84 = 0.0
        
        do idim=0,dm
           h_new_r(idim) = factor23 * dn(idim) - d23(idim)
           norm_scratch23 = norm_scratch23 + h_new_r(idim)*h_new_r(idim)
           h_new_l(idim) = factor13 * dn(idim) + d31(idim)
           norm_scratch13 = norm_scratch13 + h_new_l(idim)*h_new_l(idim)
           h_13_l(idim) = factor51 * d31(idim) - dn(idim)
           norm_scratch51 = norm_scratch51 + h_13_l(idim)*h_13_l(idim)
           h_23_l(idim) = factor63 * d23(idim) - dn(idim)
           norm_scratch63 = norm_scratch63 + h_23_l(idim)*h_23_l(idim)
           h_24_l(idim) = factor72 * d42(idim) + dn(idim)
           norm_scratch72 = norm_scratch72 + h_24_l(idim)*h_24_l(idim)
           h_14_r(idim) = factor84 * d41(idim) + dn(idim)
           norm_scratch84 = norm_scratch84 + h_14_r(idim)*h_14_r(idim)
        enddo

        norm_scratch23 = 1.0 / sqrt ( norm_scratch23 )
        norm_scratch13 = 1.0 / sqrt ( norm_scratch13 )
        norm_scratch51 = 1.0 / sqrt ( norm_scratch51 )
        norm_scratch63 = 1.0 / sqrt ( norm_scratch63 )
        norm_scratch84 = 1.0 / sqrt ( norm_scratch84 )
        norm_scratch72 = 1.0 / sqrt ( norm_scratch72 )

        edget1 = 0.0
        e51t = 0.0
        e63t = 0.0
        e72t = 0.0
        e84t = 0.0

        do idim=0,dm
           h_new_r(idim) = h_new_r(idim) * norm_scratch23
           h_new_l(idim) = h_new_l(idim) * norm_scratch13
           edget1 = edget1 + h_new_r(idim)*h_new_l(idim)
           h_13_l(idim) = h_13_l(idim) * norm_scratch51
           e51t  = e51t + rno_m(idim,m_f1)*h_13_l(idim)
           h_23_l(idim) = h_23_l(idim) * norm_scratch63
           e63t  = e63t + rno_m(idim,m_f2)*h_23_l(idim)
           h_24_l(idim) = h_24_l(idim) * norm_scratch72
           e72t  = e72t + rno_m(idim,m_f3)*h_24_l(idim)
           h_14_r(idim) = h_14_r(idim) * norm_scratch84
           e84t  = e84t + rno_m(idim,m_f4)*h_14_r(idim)
        enddo

        edge  = edge_m(link) + edge_m(m_f1) + edge_m(m_f2) +
     $       edge_m(m_f3) + edge_m(m_f4)
        edget = edget1 + e51t + e63t + e72t + e84t

        old_action = lambda * edge
        new_action = lambda * edget

        old_action = old_action + len2_o
     $    + tlog(c1)+tlog(c2)+tlog(c3)+tlog(c4)
        new_action = new_action + len2_n
     $    + tlog(c1-1)+tlog(c2-1)+tlog(c3+1)+tlog(c4+1)

        IF(new_action .GT. old_action)THEN
          IF(exp(-(new_action-old_action)).lt.ranmar())goto10
        ENDIF

        sw_lacc = sw_lacc + 1

        ma = link
        mb = m_back_m(link)

        m_former_m(m_b4) = m_f1
        m_next_m(m_f1) = m_b4
        if(link .eq. m_start_s(s1))m_start_s(s1) = m_f1

        m_former_m(m_b2) = m_f3
        m_next_m(m_f3) = m_b2
        if(mb .eq. m_start_s(s2))m_start_s(s2) = m_f3

        m_former_m(link) = m_f2
        m_next_m(link) = m_b1
        m_next_m(m_f2) = link
        m_former_m(m_b1) = link

        m_former_m(mb) = m_f4
        m_next_m(mb) = m_b3
        m_next_m(m_f4) = mb
        m_former_m(m_b3) = mb

        n_connectivity_s(s1) =  c1 - 1
        n_connectivity_s(s2) =  c2 - 1
        n_connectivity_s(s3) =  c3 + 1
        n_connectivity_s(s4) =  c4 + 1

        square_m(ma) = len2_n
        square_m(mb) = len2_n

        do idim=0,dm
           rno_m(idim,ma) = h_new_r(idim)
           rno_m(idim,mb) = h_new_l(idim)
           lno_m(idim,ma) = h_new_l(idim)
           lno_m(idim,mb) = h_new_r(idim)
           lno_m(idim,m_f1) = h_13_l(idim)
           lno_m(idim,m_f2) = h_23_l(idim)
           lno_m(idim,m_f3) = h_24_l(idim)
           lno_m(idim,m_f4) = h_14_r(idim)
           rno_m(idim,m_b1) = h_13_l(idim)
           rno_m(idim,m_b2) = h_23_l(idim)
           rno_m(idim,m_b3) = h_24_l(idim)
           rno_m(idim,m_b4) = h_14_r(idim)
        enddo

        edge_m(ma) = edget1
        edge_m(mb) = edget1
        edge_m(m_f1) = e51t
        edge_m(m_b1) = e51t
        edge_m(m_f2) = e63t
        edge_m(m_b2) = e63t
        edge_m(m_f3) = e72t
        edge_m(m_b3) = e72t
        edge_m(m_f4) = e84t
        edge_m(m_b4) = e84t

        s_site_m(ma) = s3
        s_site_m(mb) = s4
        s_neighbor_m(ma) = s4
        s_neighbor_m(mb) = s3

   10   CONTINUE

      ENDDO

      sw_lacc = inv_3n * sw_lacc 

      RETURN
  123 FORMAT(a)
      END
C_______________________________________________________________________

      SUBROUTINE update_many_sites

      INCLUDE 'string.h'

      INTEGER i_flip, m_start, s0, s1, s3, s4
      INTEGER m_forgot, i, icount, m_start1, m_start2, mem
      integer nc,r_link,l_link

      REAL ft(0:dm)
      REAL r_s1(0:dm),r_s3(0:dm)
      REAL norm_dt, dott_31, dott_41
      REAL factort31, factort41
      REAL h_dt_r(0:dm),h_dt_l(0:dm)
      REAL curvature_t
      REAL gauss, gausst, edge, edget
      REAL dt(0:dm)
      REAL dt31(0:dm),dt41(0:dm)
      REAL d23(0:dm)
      REAL scratch31, scratch41, scratch01
      REAL dot_01, factor01, curvature_t1
      REAL h_forgot_l(0:dm)

      INTEGER iv_scra(0:max_num), m_scra(0:max_num)
      REAL v_scra(0:max_num)
      REAL rno_scra(0:dm,0:max_num)
      REAL lno_scra(0:dm,0:max_num)
      REAL fno_scra(0:dm,0:max_num)
      INTEGER idim
      REAL edget_scra(0:max_num)
      REAL edget_m_scra(0:max_num)
      REAL inv_norm_dt
      real ranmar

      sw_sacc = 0.0

      DO i_flip = 0, 3*n-1

        s0 = n * (1.0-ranmar())

        IF(s0.EQ.n)THEN
          s0=nm
        ENDIF

        ft(0) = field(0,s0) + trial * (ranmar()-.5)
        ft(1) = field(1,s0) + trial * (ranmar()-.5)
        ft(2) = field(2,s0) + trial * (ranmar()-.5)

        gauss = 0.0
        gausst = 0.0
        edge = 0.0
        edget = 0.0
        icount = 0

        nc = n_connectivity_s(s0)
        m_start = m_start_s(s0)

        do i=1,nc
          iv_scra(icount) = m_start
          s1 = s_neighbor_m(m_start)
          r_link = m_former_m(m_start)
          s3 = s_neighbor_m(r_link)
          l_link = m_next_m(m_start)
          s4 = s_neighbor_m(l_link)

          m_forgot = m_former_m(m_back_m(r_link))
          m_scra(icount)  = m_forgot

          norm_dt = 0.0
          dott_31 = 0.0
          dott_41 = 0.0
          dot_01  = 0.0

          do idim=0,dm
             r_s1(idim) = field(idim,s1)
             r_s3(idim) = field(idim,s3)
             dt(idim) = r_s1(idim)  - ft(idim)
             d23(idim) = r_s1(idim) - r_s3(idim)
             dt31(idim) = r_s3(idim) - ft(idim)
             dt41(idim) = field(idim,s4) - ft(idim)
             norm_dt  = norm_dt + dt(idim)*dt(idim)
             dott_31    = dott_31 + dt(idim)*dt31(idim)
             dott_41    = dott_41 + dt(idim)*dt41(idim)
             dot_01     = dot_01 + dt(idim)*d23(idim)
          enddo

          v_scra(icount)  = norm_dt

          inv_norm_dt = 1.0/norm_dt
          factort31 = dott_31 * inv_norm_dt
          factort41 = dott_41 * inv_norm_dt
          factor01 = - dot_01/square_m(m_forgot)
          scratch31   = 0.0
          scratch41   = 0.0
          scratch01   = 0.0

          do idim=0,dm
             h_dt_r(idim) = factort31 * dt(idim) - dt31(idim)
             h_dt_l(idim) = factort41 * dt(idim) - dt41(idim)
             h_forgot_l(idim) = factor01 * d23(idim) + dt(idim)
             scratch31   = scratch31 + h_dt_r(idim) * h_dt_r(idim)
             scratch41   = scratch41 + h_dt_l(idim) * h_dt_l(idim)
             scratch01   = scratch01 + h_forgot_l(idim) 
     $            * h_forgot_l(idim)
          enddo

          scratch31   = 1.0/sqrt(scratch31)
          scratch41   = 1.0/sqrt(scratch41)
          scratch01   = 1.0/sqrt(scratch01)

          curvature_t = 0.0
          curvature_t1 = 0.0
          
          do idim=0,dm
             h_dt_r(idim) = h_dt_r(idim) * scratch31
             rno_scra(idim,icount)=h_dt_r(idim)
             h_dt_l(idim) = h_dt_l(idim) * scratch41
             lno_scra(idim,icount)=h_dt_l(idim)
             h_forgot_l(idim) = h_forgot_l(idim) * scratch01
             fno_scra(idim,icount)=h_forgot_l(idim)
             curvature_t = curvature_t + h_dt_r(idim)*h_dt_l(idim)
             curvature_t1 = curvature_t1 
     $            + rno_m(idim,m_forgot)*h_forgot_l(idim)
          enddo

          edget_scra(icount) = curvature_t
          edget_m_scra(icount) = curvature_t1

          gauss  = gauss  + square_m(m_start)
          gausst = gausst + norm_dt
          edge  = edge  + edge_m(m_start) + edge_m(m_forgot)
          edget = edget + curvature_t + curvature_t1

          m_start = m_next_m(m_start)
          icount = icount + 1
          if(icount .ge. max_num)then
             print *,'icount=',icount
             stop'too many links around'
          endif

        enddo  ! loop over the current sites

        old_action = lambda * edge +  gauss
        new_action = lambda * edget + gausst

        IF(new_action .GE. old_action)THEN
          IF(exp(-(new_action-old_action)).lt.ranmar())goto22
        ENDIF

        sw_sacc = sw_sacc + 1.0

        field(0,s0) = ft(0)
        field(1,s0) = ft(1)
        field(2,s0) = ft(2)

        DO i = 0, icount - 1
          mem = iv_scra(i)
          m_start = m_back_m(mem)
          square_m(m_start)=v_scra(i)
          square_m(mem) = v_scra(i)
          m_start1 = m_scra(i)
          m_start2 = m_back_m(m_start1)
          edge_m(mem) =  edget_scra(i)
          edge_m(m_start) =  edget_scra(i)
          edge_m(m_start1) = edget_m_scra(i)
          edge_m(m_start2) = edget_m_scra(i)

          do idim=0,dm
             rno_m(idim,mem) = rno_scra(idim,i)
             lno_m(idim,mem) = lno_scra(idim,i)
             rno_m(idim,m_start) = lno_scra(idim,i)
             lno_m(idim,m_start) = rno_scra(idim,i)
             lno_m(idim,m_start1) = fno_scra(idim,i)
             rno_m(idim,m_start2) = fno_scra(idim,i)
          enddo
        ENDDO

   22   CONTINUE
      ENDDO
      sw_sacc = sw_sacc * inv_3n
      RETURN
  123 FORMAT(a)
      END

C_______________________________________________________________________

      SUBROUTINE check_mesh

      INCLUDE 'string.h'

      INTEGER local_connectivity, total_connectivity
      INTEGER site, s0, m_start, m_end, m_start1
      INTEGER m_start2, s1, s2, counter
      integer i,j,k,nc,nc1,nc2

      total_connectivity = 0
      DO site = 0, nm
        local_connectivity = n_connectivity_s(site)
        total_connectivity = total_connectivity
     $        + local_connectivity
        IF (local_connectivity .LT. 3)THEN
          PRINT *,'TERROR!!!! site = ',site,' conn = ',
     $           local_connectivity
          stop'CHECK_MESH TERROR ABORT: small connectivity'
        ENDIF
      ENDDO

      IF(total_connectivity .NE. (n6m+1))then
        PRINT *,'TERROR!!!! total connectivity = ',
     $        total_connectivity
        stop'CHECK_MESH TERROR ABORT: wrong total connectivity'
      ENDIF

      DO s0=0, nm
        m_start = m_start_s(s0)
        nc = n_connectivity_s(s0)
        m_end   = 0
        do i=1,nc
           s1 = s_site_m(m_start)
           s2 = s_neighbor_m(m_start)
           IF(s1.NE.s0)THEN
              PRINT *,'CHECK_MESH TERROR !!!! s1.ne.s0'
              PRINT *,'s0=',s0,'s1=',s1
              stop'TERROR ABORT !!!!'
           ENDIF
           IF(s2.EQ.s1)THEN
              PRINT *,'CHECK_MESH TERROR !!!! s1.eq.s2'
              PRINT *,'s0=',s0,'s1=',s1,'s2=',s2
              stop'TERROR ABORT !!!!'
           ENDIF

          m_start2 = m_start_s(s2)
          nc1 = n_connectivity_s(s2)
          do j=1,nc1
             IF(s1.EQ.s_neighbor_m(m_start2))GOTO 52
             m_start2 = m_next_m(m_start2)
          enddo
          PRINT *,'CHECK_MESH TERROR!!!! s1     neighbor of s2,'
          PRINT *,'CHECK_MESH TERROR!!!! s2 not neighbor of s1'
          PRINT *,s1,s2
          stop'TERROR ABORT!!!!'
 52       CONTINUE

          counter = 0
          m_start1 = m_start_s(s1)
          nc2 = n_connectivity_s(s1)
          do k=1,nc2
             IF(s1.EQ.s_neighbor_m(m_start1).OR.
     $            s2.EQ.s_neighbor_m(m_start1))counter = counter + 1
             m_start1 = m_next_m(m_start1)
          enddo
          IF(counter.NE.1)THEN
             PRINT *,'CHECK_MESH TERROR !!!!',
     $            'Too many s1 and s2 in s1 list,'
             PRINT *,'s1,s2,counter',s1,s2,counter
             stop'TERROR ABORT  !!!!'
          ENDIF

          m_end = m_end+1
          m_start = m_next_m(m_start)
       enddo
      ENDDO

      RETURN
      END
C **************************************************************************
C **************************************************************************
      SUBROUTINE rmarin(ij,kl)
C This is the initialization routine for the random number generator
C     RANMAR()
C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
C                                                      0 <= KL <= 30081
      REAL u(97), c, cd, cm
      INTEGER i97, j97, ivec
      LOGICAL test
      COMMON /raset1/ u, c, cd, cm, i97, j97, ivec, test
      DATA test /.false./
      IF( ij .LT. 0  .OR.  ij .GT. 31328  .OR.
     *    kl .LT. 0  .OR.  kl .GT. 30081 ) THEN
        PRINT '(A)', ' The first random number seed must have a value
     *between 0 and 31328'
        PRINT '(A)',' The second seed must have a value between 0 and
     *30081'
        STOP
      ENDIF
      i = mod(ij/177, 177) + 2
      j = mod(ij    , 177) + 2
      k = mod(kl/169, 178) + 1
      l = mod(kl,     169)
      DO 2 ii = 1, 97
        s = 0.0
        t = 0.5
        DO 3 jj = 1, 24
          m = mod(mod(i*j, 179)*k, 179)
          i = j
          j = k
          k = m
          l = mod(53*l+1, 169)
          IF (mod(l*m, 64) .ge. 32) then
            s = s + t
          ENDIF
          t = 0.5 * t
    3   CONTINUE
        u(ii) = s
    2 CONTINUE
      c = 362436.0 / 16777216.0
      cd = 7654321.0 / 16777216.0
      cm = 16777213.0 /16777216.0
      i97 = 97
      j97 = 33
      test = .true.
      RETURN
      END
C **************************************************************************
C **************************************************************************
      FUNCTION ranmar()
      IMPLICIT NONE
      REAL rand_u(97), rand_c, rand_cd, rand_cm, rand_uni, ranmar
      INTEGER rand_i97, rand_j97, rand_ivec
      LOGICAL rand_test
      COMMON /raset1/ rand_u, rand_c, rand_cd, rand_cm, rand_i97
     $     , rand_j97, rand_ivec, rand_test

      IF( .NOT. rand_test ) THEN
        PRINT '(A)',' Call the init routine (RMARIN) before calling RAN
     *MAR'
        STOP
      ENDIF
      rand_uni = rand_u(rand_i97) - rand_u(rand_j97)
      IF( rand_uni .LT. 0.0 ) rand_uni = rand_uni + 1.0
      rand_u(rand_i97) = rand_uni
      rand_i97 = rand_i97 - 1
      IF(rand_i97 .EQ. 0) rand_i97 = 97
      rand_j97 = rand_j97 - 1
      IF(rand_j97 .EQ. 0) rand_j97 = 97
      rand_c = rand_c - rand_cd
      IF( rand_c .LT. 0.0 ) rand_c = rand_c + rand_cm
      rand_uni = rand_uni - rand_c
      IF( rand_uni .LT. 0.0 ) rand_uni = rand_uni + 1.0
      ranmar = rand_uni
      RETURN
      END
C **************************************************************************

