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     Express version by Leping Han, July 1992                     C
C     NPAC                                                         C
C     Syracuse University, Aug. 1992                               C
C                                                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      PROGRAM pippo

      INCLUDE 'string.h'
      real data_buffer(1000,n_obs)
      common /data_buffer/ data_buffer
      integer counter,row,col
      common /counter/counter

      CALL KXINIT
      CALL KXPARA(NDDATA)

      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
          counter = i_correlated
          CALL measure
          CALL average_1
        ENDDO
        IF(write_all .EQ. 1)THEN
           call kmulti(f_all)
           call kmulti(f_all2)
           if(nddata(1) .lt. 10)then
              write(f_all,150)nddata(1)
              write(f_all2,150)nddata(1)
           else
              write(f_all,151)nddata(1)
              write(f_all2,151)nddata(1)
           endif
           WRITE(f_all ,181)((data_buffer(row,col),
     $          col=1,6),row=1,n_correlated)
           WRITE(f_all2 ,181)((data_buffer(row,col),
     $          col=7,12),row=1,n_correlated)
           call kflush(f_all)
           call kflush(f_all2)
        ENDIF
        CALL average_2
      ENDDO

      CALL check_mesh

C      call g_compute_bounds
C      call g_set_up
C      call g_draw_mesh
C      call g_draw_links_with_triangles
C      call g_wait
      call ksingl(f_all)
      call ksingl(f_all2)

      CALL ciao

      STOP
  123 FORMAT(a1)
  150 format('n',i1)
  151 format('n',i2)
  181 FORMAT(1x,6g13.6)
      END

C______________________________________________________________________

      SUBROUTINE ciao

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

      call kmulti(6)
      write(6,*)'n=',nddata(1)
      write(6,*)'Flip acceptance rate _ ',tot_lacc
      write(6,*)'Site acceptance rate _ ',tot_sacc

      call kmulti(f_out)

      write(f_out,*)'n=',nddata(1)
      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)
     $              /float(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
          write(6,*)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
      call kflush(6)
      call kflush(f_out)
      call ksingl(6)
      call ksingl(f_out)
      
      call kmulti(f_oco)
      WRITE(f_oco,*)nx,ny,n,i_seed1,i_seed2
     $     ,time_scale,n_sweeps,n_correlated
     $     ,n_discard,lambda,trial,write_all
      call kflush(f_oco)
        write(f_oco,*)(rand_u(id),id=1,97)
      call kflush(f_oco)
        write(f_oco,*)rand_c, rand_cd, rand_cm
      call kflush(f_oco)
        write(f_oco,*)rand_i97, rand_j97, rand_ivec, rand_test
      call kflush(f_oco)
        write(f_oco,*)((field           (id,is),id=0,dm),is=0, nm)
      call kflush(f_oco)
        write(f_oco,*)(square_m        (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)((lno_m           (id,im),id=0,dm),im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)((rno_m           (id,im),id=0,dm),im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(m_start_s       (   is)        ,is=0, nm)
      call kflush(f_oco)
        write(f_oco,*)(m_end_s         (   is)        ,is=0, nm)
      call kflush(f_oco)
        write(f_oco,*)(n_connectivity_s(   is)        ,is=0, nm)
      call kflush(f_oco)
        write(f_oco,*)(s_site_m        (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(s_neighbor_m    (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(s_right_m       (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(s_left_m        (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(m_next_m        (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(m_former_m      (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(m_back_m        (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(c_copyx_m       (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(c_copyy_m       (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(c_coprx_m       (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(c_copry_m       (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(c_coplx_m       (   im)        ,im=0,n6m)
      call kflush(f_oco)
        write(f_oco,*)(c_coply_m       (   im)        ,im=0,n6m)
      call kflush(f_oco)
 387    format(1x,4g16.7)
      IF(write_all .EQ. 1)THEN
         CLOSE(f_all)
         CLOSE(f_all2)
      ENDIF
      call ksingl(f_oco)
      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
      integer rand_set1(32),rand_set2(32),j

C     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)    '

      color(1) = default
      color(2) = white
      color(3) = red
      color(4) = green
      color(5) = blue
      color(6) = cyan
      color(7) = magenta
      color(8) = yellow
      color(9) = black

      g_i_n_copy = 3
      g_i_n_add = g_i_n_copy/2
      g_n_copy = g_i_n_copy
      g_n_add = g_i_n_add
      offx = g_n_add * nx
      offy = g_n_add * ny

      IF(mod(g_i_n_copy,2).ne.1)then
        PRINT *, mod(g_i_n_copy, 2)
        stop'Please use odd number of copies for graphics'
      ENDIF

      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    produce the set of random number

c      i_seed1 = i_seed1 + nddata(1)
c      i_seed2 = i_seed2 + nddata(1)
     
      call RAND(i_seed1,rand_set1)
C      write(6,*)(rand_set1(i),i=1,32)
      call RAND(i_seed2,rand_set2)
C      write(6,*)(rand_set2(i),i=1,32)

      do i=1,nddata(2)
         do j=i+1,nddata(2)
            if((rand_set1(i) .eq. rand_set1(j)) .or.
     $           (rand_set2(i) .eq. rand_set2(j)))then
               stop'have equal random number seed, re_type seed'
            endif
         enddo
      enddo
               
      i_seed1 = rand_set1(nddata(1)+1)/2
      i_seed2 = rand_set2(nddata(1)+1)/2
      
      write(6,*)'== The set of random number for each processors =='
      call kmulti(6)
      write(6,*)'processor n=',nddata(1)
      write(6,*)i_seed1,i_seed2
      call kflush(6)
      call ksingl(6)


CC    assign lambda to each processor, so the following 5 line maybe 
CC    varied depending on the requirements

      if(nddata(1) .le. 8)then
         lambda = 1.40
      else 
         lambda = 0.8
      endif

      write(6,*)'== The set of lambda for each processor =='
      call kmulti(6)
      write(6,*)'processor n=',nddata(1)
      write(6,*)lambda
      call kflush(6)
      call ksingl(6)

      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')
        call kmulti(f_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,*)(rand_u(id),id=1,97)
        read(f_ico,*)rand_c, rand_cd, rand_cm
        read(f_ico,*)rand_i97, rand_j97, rand_ivec, rand_test
        read(f_ico,*)((field           (id,is),id=0,dm),is=0, nm)
        read(f_ico,*)(square_m        (   im)        ,im=0,n6m)
        read(f_ico,*)((lno_m           (id,im),id=0,dm),im=0,n6m)
        read(f_ico,*)((rno_m           (id,im),id=0,dm),im=0,n6m)
        read(f_ico,*)(m_start_s       (   is)        ,is=0, nm)
        read(f_ico,*)(m_end_s         (   is)        ,is=0, nm)
        read(f_ico,*)(n_connectivity_s(   is)        ,is=0, nm)
        read(f_ico,*)(s_site_m        (   im)        ,im=0,n6m)
        read(f_ico,*)(s_neighbor_m    (   im)        ,im=0,n6m)
        read(f_ico,*)(s_right_m       (   im)        ,im=0,n6m)
        read(f_ico,*)(s_left_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)
        read(f_ico,*)(c_copyx_m       (   im)        ,im=0,n6m)
        read(f_ico,*)(c_copyy_m       (   im)        ,im=0,n6m)
        read(f_ico,*)(c_coprx_m       (   im)        ,im=0,n6m)
        read(f_ico,*)(c_copry_m       (   im)        ,im=0,n6m)
        read(f_ico,*)(c_coplx_m       (   im)        ,im=0,n6m)
        read(f_ico,*)(c_coply_m       (   im)        ,im=0,n6m)
 387    format(1x,4g16.7)
        call ksingl(f_ico)
        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
        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
        a = rand_uni
      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

      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
            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
            field(d0,s0) = scale * (.5-rand_uni)
          ENDDO
        ENDDO
      ELSEIF(start_mode .EQ. flat_start_xy)THEN
        DO s0=0,nm
          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
          field(0,s0) = scale * (.5-rand_uni)
          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
          field(1,s0) = scale * (.5-rand_uni)
          field(2,s0) = 0.0
        ENDDO
      ELSEIF(start_mode .EQ. flat_start_xz)THEN
        DO s0=0,nm
          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
          field(0,s0) = scale * (.5-rand_uni)
          field(1,s0) = 0.0
          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
          field(2,s0) = scale * (.5-rand_uni)
        ENDDO
      ELSEIF(start_mode .EQ. flat_start_yz)THEN
        DO s0=0,nm
          field(0,s0) = 0.0
          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
          field(1,s0) = scale * (.5-rand_uni)
          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
          field(2,s0) = scale * (.5-rand_uni)
        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

      DO site = 0 , nm
        m_start_s(site) = 6 * site
        m_end_s (site) = 6 * site + 5
        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

        s_right_m(6*site+1) = site + 1
        s_right_m(6*site+2) = site + nx + 1
        s_right_m(6*site+3) = site + nx
        s_right_m(6*site+4) = site - 1
        s_right_m(6*site+5) = site - nx - 1
        s_right_m(6*site+0) = site - nx

        s_left_m(6*site+5) = site + 1
        s_left_m(6*site+0) = site + nx + 1
        s_left_m(6*site+1) = site + nx
        s_left_m(6*site+2) = site - 1
        s_left_m(6*site+3) = site - nx - 1
        s_left_m(6*site+4) = site - nx

        DO site1 = 0, 5
          c_copyx_m(6*site+site1) = 0
          c_copyy_m(6*site+site1) = 0
          c_coprx_m(6*site+site1) = 0
          c_copry_m(6*site+site1) = 0
          c_coplx_m(6*site+site1) = 0
          c_coply_m(6*site+site1) = 0
          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 ) = - 1
        m_former_m ( 6*site     ) = - 1

        IF(site .LT. nx)THEN
          s_neighbor_m(6*site+5) = site + nx*(ny-1)
          c_copyy_m(6*site+5) = -1
          s_neighbor_m(6*site+4) = site + nx*(ny-1) - 1
          c_copyy_m(6*site+4) = -1
          s_right_m(6*site+0) = site + nx*(ny-1)
          s_left_m(6*site+4) = site + nx*(ny-1)
          s_right_m(6*site+5) = site + nx*(ny-1) - 1
          s_left_m (6*site+3) = site + nx*(ny-1) - 1
          c_copry_m(6*site+0) = - 1
          c_coply_m(6*site+4) = - 1
          c_copry_m(6*site+5) = - 1
          c_coply_m(6*site+3) = - 1
        ENDIF
        IF(nx*ny-1-site .lt. nx)then
          s_neighbor_m(6*site+2) = site - nx*(ny-1)
          c_copyy_m(6*site+2) = +1
          s_neighbor_m(6*site+1) = site - nx*(ny-1) + 1
          c_copyy_m(6*site+1) = +1
          s_right_m(6*site+3) = site - nx*(ny-1)
          s_left_m(6*site+1)  = site - nx*(ny-1)
          s_right_m(6*site+2) = site - nx*(ny-1) + 1
          s_left_m(6*site+0) = site - nx*(ny-1) + 1
          c_copry_m(6*site+3) = + 1
          c_coply_m(6*site+1) = + 1
          c_copry_m(6*site+2) = + 1
          c_coply_m(6*site+0) = + 1
        ENDIF
        IF(mod(site,nx) .eq. nx-1)then
          s_neighbor_m(6*site+0) = site - nx + 1
          c_copyx_m(6*site+0) = +1
          s_neighbor_m(6*site+1) = site + 1
          c_copyx_m(6*site+1) = +1
          s_right_m(6*site+1) = site - nx + 1
          s_left_m(6*site+5) = site - nx + 1
          s_right_m(6*site+2) = site + 1
          s_left_m(6*site+0) = site + 1
          c_coprx_m(6*site+1) = + 1
          c_coplx_m(6*site+5) = + 1
          c_coprx_m(6*site+2) = + 1
          c_coplx_m(6*site+0) = + 1
        ENDIF
        IF(mod(site,nx) .eq. 0)then
          s_neighbor_m(6*site+3) = site + nx - 1
          c_copyx_m(6*site+3) = -1
          s_neighbor_m(6*site+4) = site - 1
          c_copyx_m(6*site+4) = -1
          s_right_m(6*site+4) = site + nx - 1
          s_left_m(6*site+2) = site + nx - 1
          s_right_m(6*site+5) = site - 1
          s_left_m(6*site+3) = site - 1
          c_coprx_m(6*site+4) = - 1
          c_coplx_m(6*site+2) = - 1
          c_coprx_m(6*site+5) = - 1
          c_coplx_m(6*site+3) = - 1
        ENDIF
      ENDDO
      s_neighbor_m(0+4) = nx * ny - 1
      c_copyx_m(0+4) = -1
      c_copyy_m(0+4) = -1
      s_neighbor_m(6*(nx*ny-1)+1) = 0
      c_copyx_m(6*(nx*ny-1)+1) = +1
      c_copyy_m(6*(nx*ny-1)+1) = +1

      s_right_m(0+5) = nx * ny - 1
      s_left_m(0+3) = nx * ny - 1
      c_coprx_m(0+5) = -1
      c_copry_m(0+5) = -1
      c_coplx_m(0+3) = -1
      c_coply_m(0+3) = -1
      s_right_m(6*(nx*ny-1)+2) = 0
      s_left_m(6*(nx*ny-1)+0) = 0
      c_coprx_m(6*(nx*ny-1)+2) = +1
      c_copry_m(6*(nx*ny-1)+2) = +1
      c_coplx_m(6*(nx*ny-1)+0) = +1
      c_coply_m(6*(nx*ny-1)+0) = +1

      DO s0=0, nm
        m_start = m_start_s(s0)
   41   CONTINUE
        IF(m_start .GE. 0)THEN
          s2 = s_neighbor_m(m_start)
          m_start2 = m_start_s(s2)
   51     CONTINUE
          IF(m_start2 .GE. 0)THEN
            IF(s0.EQ.s_neighbor_m(m_start2))GOTO 52
            m_start2 = m_next_m(m_start2)
            GOTO 51
          ENDIF
   52     CONTINUE

C         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)
          GOTO 41
        ENDIF
      ENDDO

      RETURN
      END
C_______________________________________________________________________

      SUBROUTINE average_1

      INCLUDE 'string.h'
      INTEGER io

      ind_lacc   = ind_lacc  + sw_lacc/float(n_correlated)
      ind_sacc   = ind_sacc  + sw_sacc/float(n_correlated)

      DO io=1,n_obs
        ind_obs(1,io) = ind_obs(1,io)
     $        +sw_obs(io)/float(n_correlated)
        ind_obs(2,io) = ind_obs(2,io)
     $        +sw_obs(io)*sw_obs(io)/float(n_correlated)
      ENDDO

      RETURN
      END

C_______________________________________________________________________

      SUBROUTINE average_2

      INCLUDE 'string.h'
      INTEGER i, io

      
      write(6,*) 'i_sweep= ', i_sweep
      call kmulti(6)
      write(6,*)'processor n=',nddata(1)
      write(6,*)'lacc,sacc ',ind_lacc, ind_sacc
      DO io=1,n_obs
        write(6,*) phrase(io),ind_obs(1,io), ind_obs(2,io)
      ENDDO
      call kflush(6)
      call ksingl(6)

      tot_lacc  = tot_lacc + ind_lacc/float(n_sweeps)
      tot_sacc  = tot_sacc + ind_sacc/float(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)/float(n_sweeps)
          tot_obs(i,2,io) = tot_obs(i,2,io)
     $           + ind_obs(i,io)*ind_obs(i,io)/float(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 s00, s2, m_s

      INTEGER s0, s1, s3, s4, m_start
      REAL dv0, dv1, dv2, f0, f1, f2, g0, g1, g2
      REAL d31_0, d31_1, d31_2, d41_0, d41_1, d41_2, norm_do
      REAL ab_31, factor31, norm31_0, norm31_1, norm31_2
      REAL norm_old_right, norm_old_left
      REAL ab_41, factor41, norm41_0, norm41_1, norm41_2
      REAL inorm_31, inorm_41

      DO s00=0, nm
        m_start = m_start_s(s00)
   41   CONTINUE
        IF(m_start .GE. 0)THEN
          s2 = s_neighbor_m(m_start)
          m_s = m_start
          s0 = s_site_m(m_s)
          f0 = field(0,s0)
          f1 = field(1,s0)
          f2 = field(2,s0)
          s1 = s_neighbor_m(m_s)
          s3 = s_right_m(m_s)
          s4 = s_left_m(m_s)
          g0 = field(0,s1)
          g1 = field(1,s1)
          g2 = field(2,s1)
          dv0 = g0  - f0
          dv1 = g1  - f1
          dv2 = g2  - f2
          d31_0 = field(0,s3) - f0
          d31_1 = field(1,s3) - f1
          d31_2 = field(2,s3) - f2
          d41_0 = field(0,s4) - f0
          d41_1 = field(1,s4) - f1
          d41_2 = field(2,s4) - f2
          norm_do  = dv0*dv0   +  dv1*dv1   + dv2*dv2
          ab_31    = dv0*d31_0 +  dv1*d31_1 + dv2*d31_2
          factor31 = ab_31/norm_do
          norm31_0 = factor31 * dv0 - d31_0
          norm31_1 = factor31 * dv1 - d31_1
          norm31_2 = factor31 * dv2 - d31_2
          ab_41    = dv0*d41_0 +  dv1*d41_1 + dv2*d41_2
          factor41 = ab_41/norm_do
          norm41_0 = factor41 * dv0 - d41_0
          norm41_1 = factor41 * dv1 - d41_1
          norm41_2 = factor41 * dv2 - d41_2

          norm_old_right = norm31_0*norm31_0
     $         + norm31_1*norm31_1 + norm31_2*norm31_2
          norm_old_left = norm41_0*norm41_0
     $         + norm41_1*norm41_1 + norm41_2*norm41_2

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

          rno_m(0,m_s) = norm31_0 * inorm_31
          rno_m(1,m_s) = norm31_1 * inorm_31
          rno_m(2,m_s) = norm31_2 * inorm_31
          lno_m(0,m_s) = norm41_0 * inorm_41
          lno_m(1,m_s) = norm41_1 * inorm_41
          lno_m(2,m_s) = norm41_2 * inorm_41

          m_start = m_next_m(m_start)

          GOTO 41
        ENDIF
      ENDDO

      RETURN
      END
C_______________________________________________________________________

      SUBROUTINE measure

      INCLUDE 'string.h'

      INTEGER io, id, s0, m0, m_start, nc, sn, sr
      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, curvature
      REAL bari(0:dm), scratch(0:dm)
      real data_buffer(1000,n_obs)
      common /data_buffer/ data_buffer
      integer counter
      common /counter/counter

      DO io=1, n_obs
        sw_obs(io) = 0.0
      ENDDO

      qmax = 0.0
      DO id=0,dm
        bari(id) = 0.0
      ENDDO
      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))
        DO id=0,dm
          bari(id) = bari(id) + field(id,s0)
        ENDDO
      ENDDO
      DO id=0,dm
        bari(id) = bari(id)/float(n)
      ENDDO
      DO s0 = 0, nm
        DO id=0,dm
          scratch(id) = field(id,s0)-bari(id)
C            field(id,s0) = scratch(id)
          sw_obs(5) = sw_obs(5) + scratch(id)**2
        ENDDO
      ENDDO
      sw_obs(1) = ( two_pi / 3. ) * sw_obs(1) / float(n)
      sw_obs(2) = sw_obs(2) / float(n)
      sw_obs(4) = qmax
      sw_obs(5) = sw_obs(5) / float(n)

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

      sw_obs(3) = sw_obs(3) / 2.0

      DO s0=0, nm
        m_start = m_start_s(s0)
   12   CONTINUE
        IF(m_start .GE. 0)THEN
          curvature = 1.
     $      + lno_m(0,m_start)*rno_m(0,m_start)
     $      + lno_m(1,m_start)*rno_m(1,m_start)
     $      + lno_m(2,m_start)*rno_m(2,m_start)
          sw_obs(6) = sw_obs(6) + curvature
          m_start = m_next_m(m_start)
          GOTO 12
        ENDIF
      ENDDO
      sw_obs(6) = sw_obs(6) / float(n) / 2.

      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)
   22   CONTINUE
        IF(m_start .GE. 0)THEN
          sn = s_neighbor_m(m_start)
          sr = s_right_m(m_start)
          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

          arg_scratch = square_m(m_start)
          if(arg_scratch .ge. 0)then
             sq1 = sqrt(arg_scratch)
          else
             print *, arg_scratch
             print *,'MEASURE SOS, SQ(M_START)<0, SQ1=0 SET'
             sq1=0.0
          endif

          arg_scratch = square_r
          if(arg_scratch .ge. 0)then
             sq2 = sqrt(arg_scratch)
          else
             print *, arg_scratch
             print *,'MEASURE SOS, SQUARE_R<0, SQ2=0 SET'
             sq2=0.0
          endif
             
          arg_scratch1 = sq1 * sq2
          if(arg_scratch1 .gt. 0)then
             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**2
          if(arg_scratch .ge. 0)then
             area = 0.5 * arg_scratch1 * sqrt(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

          if(abs(scalar_p).le.1.0)then
             site_d_angle = site_d_angle + 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 / 4.0
          endif
          site_d_area  = site_d_area  + area

          m_start = m_next_m(m_start)
          GOTO 22
        ENDIF
        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 - two_pi/6.*float(nc))
        sw_obs(10) = sw_obs(10) + float(nc) / 3. 
     $       * ( two_pi**2 * ((6-float(nc))/float(nc))**2)
        sw_obs(11) = sw_obs(11) + float(nc) * 3.
     $       * ( (two_pi-site_d_angle) / site_d_area ) ** 2
        sw_obs(12) = sw_obs(12) + float(nc)
     $       * ( two_pi * ((6.-float(nc))/float(nc))
     $       * ( 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) / float(n)
      sw_obs(8) = sw_obs(8) / float(n6)
      sw_obs(9) = sw_obs(9) / float(n)
      sw_obs(10) = sw_obs(10) / float(n)
      sw_obs(11) = sw_obs(11) / float(n)

      do io=1,12
         data_buffer(counter,io) = sw_obs(io)
      enddo


      RETURN
  123 FORMAT(a1)
      END

c______________________________________________________________________

      SUBROUTINE flip_many_links

      INCLUDE 'string.h'

      INTEGER link, s1, s2, s3, s4, s0, s6
      INTEGER c1, c2, c3, c4, c0, c6, c3x, c3y, c4x, c4y
      INTEGER nc1, nc2, nc3, nc4
      INTEGER i_flip
      INTEGER m_start, site, ma ,mb, me3, me4, mnma, mnmb
      INTEGER mfma, mfmb
      INTEGER m_start31, m_start32, m_start41, m_start42
      INTEGER m_start13, m_start23, m_start14, m_start24
      INTEGER m_f1, m_f2, m_f3, m_f4
      INTEGER m_b1, m_b2, m_b3, m_b4

      REAL dn0, dn1, dn2, norm_scratch, d42_0, d42_1, d42_2
      REAL norm_d31, norm_d23, norm_d42, norm_d14
      REAL d31_0, d31_1, d31_2
      REAL norm31_0, norm31_1, norm31_2
      REAL ab_13, factor13
      REAL norm13_0, norm13_1, norm13_2
      REAL d41_0, d41_1, d41_2
      REAL norm41_0, norm41_1, norm41_2
      REAL d23_0, d23_1, d23_2, ab_23, factor23
      REAL norm23_0, norm23_1, norm23_2
      REAL ab_51,factor51, norm51_0, norm51_1, norm51_2
      REAL ab_63,factor63, norm63_0, norm63_1, norm63_2
      REAL ab_72,factor72, norm72_0, norm72_1, norm72_2
      REAL ab_84,factor84, norm84_0, norm84_1, norm84_2
      REAL norm51t_0, norm51t_1, norm51t_2
      REAL norm63t_0, norm63t_1, norm63t_2
      REAL norm72t_0, norm72t_1, norm72t_2
      REAL norm84t_0, norm84t_1, norm84t_2
      REAL norm51o_0, norm51o_1, norm51o_2
      REAL norm63o_0, norm63o_1, norm63o_2
      REAL norm72o_0, norm72o_1, norm72o_2
      REAL norm84o_0, norm84o_1, norm84o_2
      REAL e51o, e51t
      REAL e63o, e63t
      REAL e72o, e72t
      REAL e84o, e84t

      REAL norm_do, norm_dn, edge, edget

      sw_lacc = 0

      DO i_flip = 1, 3 * n

        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
        link = n6 * (1.0-rand_uni)

        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
        s3 = s_right_m(link)
        s4 = s_left_m(link)

        c3 = n_connectivity_s(s3)
        c4 = n_connectivity_s(s4)
        c3x = c_coprx_m(link)
        c3y = c_copry_m(link)
        c4x = c_coplx_m(link)
        c4y = c_coply_m(link)
C     now checking that 3 and 4 are not already linked
        IF(c4 .LT. c3) THEN
          s0 = s4
          c0 = c4
          s6 = s3
          c6 = c3
        ELSE
          s0 = s3
          c0 = c3
          s6 = s4
          c6 = c4
        ENDIF
        m_start = m_start_s(s0)
   11   CONTINUE
        IF(m_start .GE. 0)THEN
          site = s_neighbor_m(m_start)
          IF(site .EQ. s6)GOTO 10
          m_start = m_next_m(m_start)
          GOTO 11
        ENDIF
C Here we know that topology is ok: we can accept.
C Now we check the action.

        m_f1 = m_start_s(s1)
  212   CONTINUE
        IF(m_f1 .GE. 0)THEN
          IF(s_neighbor_m(m_f1) .EQ. s3)GOTO 213
          m_f1 = m_next_m(m_f1)
          GOTO 212
        ENDIF
        PRINT *,'LINK UPDATE ABORT: forgot and not found 1'
        stop'TERROR ABORT !!!'
  213   CONTINUE

        m_f2 = m_start_s(s3)
  312   CONTINUE
        IF(m_f2 .GE. 0)THEN
          IF(s_neighbor_m(m_f2) .EQ. s2)GOTO 313
          m_f2 = m_next_m(m_f2)
          GOTO 312
        ENDIF
        PRINT *,'LINK UPDATE ABORT: forgot and not found 2'
        stop'TERROR ABORT       !!!'
  313   CONTINUE

        m_f3 = m_start_s(s2)
  412   CONTINUE
        IF(m_f3 .GE. 0)THEN
          IF(s_neighbor_m(m_f3) .EQ. s4)GOTO 413
          m_f3 = m_next_m(m_f3)
          GOTO 412
        ENDIF
        PRINT *,'LINK UPDATE ABORT: forgot and not found 3'
        stop'TERROR ABORT !!!'
  413   CONTINUE

        m_f4 = m_start_s(s4)
  512   CONTINUE
        IF(m_f4 .GE. 0)THEN
          IF(s_neighbor_m(m_f4) .EQ. s1)GOTO 513
          m_f4 = m_next_m(m_f4)
          GOTO 512
        ENDIF
        PRINT *,'LINK UPDATE ABORT: forgot and not found 4'
        stop'TERROR ABORT !!!'
  513   CONTINUE

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

        d31_0 = field(0,s3) - field(0,s1)
        d31_1 = field(1,s3) - field(1,s1)
        d31_2 = field(2,s3) - field(2,s1)

        d41_0 = field(0,s4) - field(0,s1)
        d41_1 = field(1,s4) - field(1,s1)
        d41_2 = field(2,s4) - field(2,s1)

        dn0 = field(0,s4) - field(0,s3)
        dn1 = field(1,s4) - field(1,s3)
        dn2 = field(2,s4) - field(2,s3)

        d23_0 = field(0,s2) - field(0,s3)
        d23_1 = field(1,s2) - field(1,s3)
        d23_2 = field(2,s2) - field(2,s3)

        d42_0 = field(0,s4) - field(0,s2)
        d42_1 = field(1,s4) - field(1,s2)
        d42_2 = field(2,s4) - field(2,s2)

        norm_do  = square_m(link)
        norm_dn  = dn0*dn0   +  dn1*dn1   + dn2*dn2

        norm_d31  = square_m(m_f1)
        norm_d23  = square_m(m_f2)
        norm_d42  = square_m(m_f3)
        norm_d14  = square_m(m_f4)

        norm31_0 = rno_m(0,link)
        norm31_1 = rno_m(1,link)
        norm31_2 = rno_m(2,link)

        norm41_0 = lno_m(0,link)
        norm41_1 = lno_m(1,link)
        norm41_2 = lno_m(2,link)

        ab_23    = dn0*d23_0 +  dn1*d23_1 + dn2*d23_2
        factor23 = ab_23/norm_dn
        norm23_0 = factor23 * dn0 - d23_0
        norm23_1 = factor23 * dn1 - d23_1
        norm23_2 = factor23 * dn2 - d23_2
        norm_scratch = norm23_0*norm23_0
     $       + norm23_1*norm23_1 + norm23_2*norm23_2
        norm_scratch = 1.0 / sqrt ( norm_scratch )
        norm23_0 = norm23_0 * norm_scratch
        norm23_1 = norm23_1 * norm_scratch
        norm23_2 = norm23_2 * norm_scratch

        ab_13    = dn0*d31_0 +  dn1*d31_1 + dn2*d31_2
        factor13 = - ab_13/norm_dn
        norm13_0 = factor13 * dn0 + d31_0
        norm13_1 = factor13 * dn1 + d31_1
        norm13_2 = factor13 * dn2 + d31_2
        norm_scratch = norm13_0*norm13_0
     $       + norm13_1*norm13_1 + norm13_2*norm13_2
        norm_scratch = 1.0 / sqrt ( norm_scratch )
        norm13_0 = norm13_0 * norm_scratch
        norm13_1 = norm13_1 * norm_scratch
        norm13_2 = norm13_2 * norm_scratch
C        print *,'norm13', norm13_0, norm13_1, norm13_2

        norm51_0 = rno_m(0,m_f1)
        norm51_1 = rno_m(1,m_f1)
        norm51_2 = rno_m(2,m_f1)

        norm51o_0 = lno_m(0,m_f1)
        norm51o_1 = lno_m(1,m_f1)
        norm51o_2 = lno_m(2,m_f1)

        ab_51    = d31_0*d41_0 +  d31_1*d41_1 + d31_2*d41_2
        factor51 = ab_51/norm_d31
        norm51t_0 = factor51 * d31_0 - d41_0
        norm51t_1 = factor51 * d31_1 - d41_1
        norm51t_2 = factor51 * d31_2 - d41_2
        norm_scratch = norm51t_0*norm51t_0
     $       + norm51t_1*norm51t_1 + norm51t_2*norm51t_2
        norm_scratch = 1.0 / sqrt ( norm_scratch )
        norm51t_0 = norm51t_0 * norm_scratch
        norm51t_1 = norm51t_1 * norm_scratch
        norm51t_2 = norm51t_2 * norm_scratch

        norm63_0 = rno_m(0,m_f2)
        norm63_1 = rno_m(1,m_f2)
        norm63_2 = rno_m(2,m_f2)

        norm63o_0 = lno_m(0,m_f2)
        norm63o_1 = lno_m(1,m_f2)
        norm63o_2 = lno_m(2,m_f2)

        ab_63    = d23_0*dn0 +  d23_1*dn1 + d23_2*dn2
        factor63 = ab_63/norm_d23
        norm63t_0 = factor63 * d23_0 - dn0
        norm63t_1 = factor63 * d23_1 - dn1
        norm63t_2 = factor63 * d23_2 - dn2
        norm_scratch = norm63t_0*norm63t_0
     $       + norm63t_1*norm63t_1 + norm63t_2*norm63t_2
        norm_scratch = 1.0 / sqrt ( norm_scratch )
        norm63t_0 = norm63t_0 * norm_scratch
        norm63t_1 = norm63t_1 * norm_scratch
        norm63t_2 = norm63t_2 * norm_scratch

        norm72_0 = rno_m(0,m_f3)
        norm72_1 = rno_m(1,m_f3)
        norm72_2 = rno_m(2,m_f3)

        norm72o_0 = lno_m(0,m_f3)
        norm72o_1 = lno_m(1,m_f3)
        norm72o_2 = lno_m(2,m_f3)

        ab_72    = d42_0*d23_0 +  d42_1*d23_1 + d42_2*d23_2
        factor72 = - ab_72/norm_d42
        norm72t_0 = factor72 * d42_0 + d23_0
        norm72t_1 = factor72 * d42_1 + d23_1
        norm72t_2 = factor72 * d42_2 + d23_2
        norm_scratch = norm72t_0*norm72t_0
     $       + norm72t_1*norm72t_1 + norm72t_2*norm72t_2
        norm_scratch = 1.0 / sqrt ( norm_scratch )
        norm72t_0 = norm72t_0 * norm_scratch
        norm72t_1 = norm72t_1 * norm_scratch
        norm72t_2 = norm72t_2 * norm_scratch

        norm84_0 = rno_m(0,m_f4)
        norm84_1 = rno_m(1,m_f4)
        norm84_2 = rno_m(2,m_f4)

        norm84o_0 = lno_m(0,m_f4)
        norm84o_1 = lno_m(1,m_f4)
        norm84o_2 = lno_m(2,m_f4)

        ab_84    = d41_0*dn0 +  d41_1*dn1 + d41_2*dn2
        factor84 = - ab_84/norm_d14
        norm84t_0 = factor84 * d41_0 + dn0
        norm84t_1 = factor84 * d41_1 + dn1
        norm84t_2 = factor84 * d41_2 + dn2
        norm_scratch = norm84t_0*norm84t_0
     $       + norm84t_1*norm84t_1 + norm84t_2*norm84t_2
        norm_scratch = 1.0 / sqrt ( norm_scratch )
        norm84t_0 = norm84t_0 * norm_scratch
        norm84t_1 = norm84t_1 * norm_scratch
        norm84t_2 = norm84t_2 * norm_scratch
C        print *,'norm84t', norm84t_0, norm84t_1, norm84t_2

        edge = norm31_0*norm41_0 +norm31_1*norm41_1 +norm31_2*norm41_2
        e51o = norm51_0*norm51o_0+norm51_1*norm51o_1+norm51_2*norm51o_2
        e63o = norm63_0*norm63o_0+norm63_1*norm63o_1+norm63_2*norm63o_2
        e72o = norm72_0*norm72o_0+norm72_1*norm72o_1+norm72_2*norm72o_2
        e84o = norm84_0*norm84o_0+norm84_1*norm84o_1+norm84_2*norm84o_2

        edget = norm23_0*norm13_0 +norm23_1*norm13_1 +norm23_2*norm13_2
        e51t  = norm51_0*norm51t_0+norm51_1*norm51t_1+norm51_2*norm51t_2
        e63t  = norm63_0*norm63t_0+norm63_1*norm63t_1+norm63_2*norm63t_2
        e72t  = norm72_0*norm72t_0+norm72_1*norm72t_1+norm72_2*norm72t_2
        e84t  = norm84_0*norm84t_0+norm84_1*norm84t_1+norm84_2*norm84t_2

        edge  = edge  + e51o + e63o + e72o + e84o
        edget = edget + e51t + e63t + e72t + e84t

        old_action = lambda * edge
        new_action = lambda * edget

        old_action = old_action + norm_do
     $    + tlog(c1)+tlog(c2)+tlog(c3)+tlog(c4)
        new_action = new_action + norm_dn
     $    + tlog(c1-1)+tlog(c2-1)+tlog(c3+1)+tlog(c4+1)
        IF(new_action .GT. old_action)THEN
          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
          IF(exp(-(new_action-old_action)).lt.rand_uni)goto10
        ENDIF

        sw_lacc = sw_lacc + 1
        n_connectivity_s(s1) =  n_connectivity_s(s1) - 1
        n_connectivity_s(s2) =  n_connectivity_s(s2) - 1
        n_connectivity_s(s3) =  n_connectivity_s(s3) + 1
        n_connectivity_s(s4) =  n_connectivity_s(s4) + 1
        nc1 = n_connectivity_s(s1)
        nc2 = n_connectivity_s(s2)
        nc3 = n_connectivity_s(s3)
        nc4 = n_connectivity_s(s4)
        ma = link
        mb = m_back_m(link)

        square_m(ma) = norm_dn
        square_m(mb) = norm_dn

        rno_m(0,ma) = norm23_0
        rno_m(0,mb) = norm13_0
        lno_m(0,ma) = norm13_0
        lno_m(0,mb) = norm23_0
        lno_m(0,m_f1) = norm51t_0
        lno_m(0,m_f2) = norm63t_0
        lno_m(0,m_f3) = norm72t_0
        lno_m(0,m_f4) = norm84t_0
        rno_m(0,m_b1) = norm51t_0
        rno_m(0,m_b2) = norm63t_0
        rno_m(0,m_b3) = norm72t_0
        rno_m(0,m_b4) = norm84t_0

        rno_m(1,ma) = norm23_1
        rno_m(1,mb) = norm13_1
        lno_m(1,ma) = norm13_1
        lno_m(1,mb) = norm23_1
        lno_m(1,m_f1) = norm51t_1
        lno_m(1,m_f2) = norm63t_1
        lno_m(1,m_f3) = norm72t_1
        lno_m(1,m_f4) = norm84t_1
        rno_m(1,m_b1) = norm51t_1
        rno_m(1,m_b2) = norm63t_1
        rno_m(1,m_b3) = norm72t_1
        rno_m(1,m_b4) = norm84t_1

        rno_m(2,ma) = norm23_2
        rno_m(2,mb) = norm13_2
        lno_m(2,ma) = norm13_2
        lno_m(2,mb) = norm23_2
        lno_m(2,m_f1) = norm51t_2
        lno_m(2,m_f2) = norm63t_2
        lno_m(2,m_f3) = norm72t_2
        lno_m(2,m_f4) = norm84t_2
        rno_m(2,m_b1) = norm51t_2
        rno_m(2,m_b2) = norm63t_2
        rno_m(2,m_b3) = norm72t_2
        rno_m(2,m_b4) = norm84t_2

        s_site_m(ma) = s3
        s_site_m(mb) = s4
        s_neighbor_m(ma) = s4
        s_neighbor_m(mb) = s3
        c_copyx_m(ma) = c4x - c3x
        c_copyy_m(ma) = c4y - c3y
        c_copyx_m(mb) = c3x - c4x
        c_copyy_m(mb) = c3y - c4y

        me3 = m_end_s(s3)
        me4 = m_end_s(s4)
        mnma = m_next_m(ma)
        mnmb = m_next_m(mb)
        mfma = m_former_m(ma)
        mfmb = m_former_m(mb)

        m_end_s(s3) = ma
        m_end_s(s4) = mb
        m_former_m(ma) = me3
        m_former_m(mb) = me4
        m_next_m(ma) = -1
        m_next_m(mb) = -1
        m_next_m(me3) = ma
        m_next_m(me4) = mb

        IF(mfma .NE. -1)THEN
          m_next_m(mfma) = mnma
        ELSE
          m_start_s(s1) = mnma
        ENDIF

        IF(mfmb .NE. -1)THEN
          m_next_m(mfmb) = mnmb
        ELSE
          m_start_s(s2) = mnmb
        ENDIF

        IF(mnma .NE. -1)THEN
          m_former_m(mnma) = mfma
        ELSE
          m_end_s(s1) = mfma
        ENDIF

        IF(mnmb .NE. -1)THEN
          m_former_m(mnmb) = mfmb
        ELSE
          m_end_s(s2) = mfmb
        ENDIF

        m_start13 = m_f1
        m_start31 = m_b1
        m_start32 = m_f2
        m_start23 = m_b2
        m_start24 = m_f3
        m_start42 = m_b3
        m_start41 = m_f4
        m_start14 = m_b4

        s_right_m(ma) = s2
        s_left_m(ma)  = s1
        s_right_m(mb) = s1
        s_left_m(mb)  = s2
        s_right_m(m_start14) = s3
        s_left_m(m_start41) = s3
        s_left_m(m_start13) = s4
        s_right_m(m_start31) = s4
        s_right_m(m_start42) = s3
        s_left_m(m_start24) = s3
        s_left_m(m_start32) = s4
        s_right_m(m_start23) = s4

        c_coprx_m(ma) = c_copyx_m(m_start32)
        c_copry_m(ma) = c_copyy_m(m_start32)
        c_coplx_m(ma) = c_copyx_m(m_start31)
        c_coply_m(ma) = c_copyy_m(m_start31)
        c_coprx_m(mb) = c_copyx_m(m_start41)
        c_copry_m(mb) = c_copyy_m(m_start41)
        c_coplx_m(mb) = c_copyx_m(m_start42)
        c_coply_m(mb) = c_copyy_m(m_start42)

        c_coprx_m(m_start14) = c_copyx_m(m_start13)
        c_copry_m(m_start14) = c_copyy_m(m_start13)
        c_coplx_m(m_start41) = c_copyx_m(mb)
        c_coply_m(m_start41) = c_copyy_m(mb)
        c_coplx_m(m_start13) = c_copyx_m(m_start14)
        c_coply_m(m_start13) = c_copyy_m(m_start14)
        c_coprx_m(m_start31) = c_copyx_m(ma)
        c_copry_m(m_start31) = c_copyy_m(ma)
        c_coprx_m(m_start42) = c_copyx_m(mb)
        c_copry_m(m_start42) = c_copyy_m(mb)
        c_coplx_m(m_start24) = c_copyx_m(m_start23)
        c_coply_m(m_start24) = c_copyy_m(m_start23)
        c_coplx_m(m_start32) = c_copyx_m(ma)
        c_coply_m(m_start32) = c_copyy_m(ma)
        c_coprx_m(m_start23) = c_copyx_m(m_start24)
        c_copry_m(m_start23) = c_copyy_m(m_start24)

   10   CONTINUE

      ENDDO

      sw_lacc = sw_lacc / (3.*n)

      RETURN
  123 FORMAT(a)
      END
C_______________________________________________________________________

      SUBROUTINE update_many_sites

      INCLUDE 'string.h'

      INTEGER i_flip, m_start, s0, s1, s3, s4, s5
      INTEGER m_forgot, i, icount, m_start1, m_start2, mem

      REAL ft0, ft1, ft2
      REAL f0, f1, f2
      REAL g0, g1, g2
      REAL h0, h1, h2
      REAL norm_do, norm_dt, abt_31, abt_41
      REAL norm_d23
      REAL factort31, factort41
      REAL norm31_0, norm31_1, norm31_2
      REAL norm41_0, norm41_1, norm41_2
      REAL normt31_0, normt31_1, normt31_2
      REAL normt41_0, normt41_1, normt41_2
      REAL curvature, curvature_t
      REAL gauss, gausst, edge, edget
      REAL dt0, dt1, dt2
      REAL dt31_0, dt31_1, dt31_2, dt41_0, dt41_1, dt41_2
      REAL d23_0, d23_1, d23_2
      REAL scratch
      REAL ab_01, factor01, curvature1, curvature_t1
      REAL norm01_0, norm01_1, norm01_2
      REAL normt01_0, normt01_1, normt01_2
      REAL norm51_0, norm51_1, norm51_2

      INTEGER iv_scra(0:n6m), m_scra(0:n6m)
      REAL v_scra(0:n6m)
      REAL rno_scra(0:dm,0:n6m)
      REAL lno_scra(0:dm,0:n6m)
      REAL fno_scra(0:dm,0:n6m)

      sw_sacc = 0.0

      DO i_flip = 0, 3*n-1

        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
        s0 = n * (1.0-rand_uni)

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

c        s0 = mod(i_flip,n)

        f0 = field(0,s0)
        f1 = field(1,s0)
        f2 = field(2,s0)

        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
        ft0 = f0 + trial * (rand_uni-.5)
        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
        ft1 = f1 + trial * (rand_uni-.5)
        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
        ft2 = f2 + trial * (rand_uni-.5)
        gauss = 0.0
        gausst = 0.0
        edge = 0.0
        edget = 0.0
        icount = 0

        m_start = m_start_s(s0)
   12   CONTINUE

        IF(m_start .GE. 0)THEN

          s1 = s_neighbor_m(m_start)
          s3 = s_right_m   (m_start)
   11     s4 = s_left_m    (m_start)

          m_forgot = m_start_s(s3)
  212     CONTINUE
          IF(m_forgot .GE. 0)THEN ! 1686
            IF(s_neighbor_m(m_forgot) .EQ. s1)GOTO 213
            m_forgot = m_next_m(m_forgot)
            GOTO 212
          ENDIF
          PRINT *,'UPDATE ABORT: forgot and not found'
          stop'TERROR ABORT !!!'
  213     CONTINUE

          s5  = s_right_m(m_forgot)

          g0 = field(0,s1)
          g1 = field(1,s1)
          g2 = field(2,s1)
          h0 = field(0,s3)
          h1 = field(1,s3)
          h2 = field(2,s3)
          d23_0 = g0 - h0
          d23_1 = g1 - h1
          d23_2 = g2 - h2
          dt0 = g0  - ft0
          dt1 = g1  - ft1
          dt2 = g2  - ft2
          dt31_0 = h0 - ft0
          dt31_1 = h1 - ft1
          dt31_2 = h2 - ft2
          dt41_0 = field(0,s4) - ft0
          dt41_1 = field(1,s4) - ft1
          dt41_2 = field(2,s4) - ft2

          norm_do = square_m(m_start)
          norm_dt  = dt0*dt0   +  dt1*dt1   + dt2*dt2
          iv_scra(icount) = m_start
          v_scra(icount)  = norm_dt
          m_scra(icount)  = m_forgot

          norm_d23  = d23_0*d23_0   +  d23_1*d23_1   + d23_2*d23_2

          norm31_0 = rno_m(0,m_start)
          norm31_1 = rno_m(1,m_start)
          norm31_2 = rno_m(2,m_start)

          norm41_0 = lno_m(0,m_start)
          norm41_1 = lno_m(1,m_start)
          norm41_2 = lno_m(2,m_start)

          abt_31    = dt0*dt31_0 +  dt1*dt31_1 + dt2*dt31_2
          factort31 = abt_31/norm_dt
          normt31_0 = factort31 * dt0 - dt31_0
          normt31_1 = factort31 * dt1 - dt31_1
          normt31_2 = factort31 * dt2 - dt31_2
          scratch   = normt31_0 * normt31_0
     $              + normt31_1 * normt31_1
     $              + normt31_2 * normt31_2
          scratch   = sqrt(scratch)
          normt31_0 = normt31_0 / scratch
          normt31_1 = normt31_1 / scratch
          normt31_2 = normt31_2 / scratch
          rno_scra(0,icount)=normt31_0
          rno_scra(1,icount)=normt31_1
          rno_scra(2,icount)=normt31_2

          abt_41    = dt0*dt41_0 +  dt1*dt41_1 + dt2*dt41_2
          factort41 = abt_41/norm_dt
          normt41_0 = factort41 * dt0 - dt41_0
          normt41_1 = factort41 * dt1 - dt41_1
          normt41_2 = factort41 * dt2 - dt41_2
          scratch   = normt41_0 * normt41_0
     $              + normt41_1 * normt41_1
     $              + normt41_2 * normt41_2
          scratch   = sqrt(scratch)
          normt41_0 = normt41_0 / scratch
          normt41_1 = normt41_1 / scratch
          normt41_2 = normt41_2 / scratch
          lno_scra(0,icount)=normt41_0
          lno_scra(1,icount)=normt41_1
          lno_scra(2,icount)=normt41_2

          curvature =
     $      norm31_0*norm41_0+norm31_1*norm41_1+norm31_2*norm41_2
          curvature_t =
     $      normt31_0*normt41_0+normt31_1*normt41_1+normt31_2*normt41_2

          norm51_0 = rno_m(0,m_forgot)
          norm51_1 = rno_m(1,m_forgot)
          norm51_2 = rno_m(2,m_forgot)

          norm01_0 = lno_m(0,m_forgot)
          norm01_1 = lno_m(1,m_forgot)
          norm01_2 = lno_m(2,m_forgot)

          ab_01    = dt31_0*d23_0 +  dt31_1*d23_1 + dt31_2*d23_2
          factor01 = - ab_01/norm_d23
          normt01_0 = factor01 * d23_0 + dt31_0
          normt01_1 = factor01 * d23_1 + dt31_1
          normt01_2 = factor01 * d23_2 + dt31_2
          scratch   = normt01_0 * normt01_0
     $              + normt01_1 * normt01_1
     $              + normt01_2 * normt01_2
          scratch   = sqrt(scratch)
          normt01_0 = normt01_0 / scratch
          normt01_1 = normt01_1 / scratch
          normt01_2 = normt01_2 / scratch

          fno_scra(0,icount)=normt01_0
          fno_scra(1,icount)=normt01_1
          fno_scra(2,icount)=normt01_2

          curvature1 =
     $      norm51_0*norm01_0+norm51_1*norm01_1+norm51_2*norm01_2
          curvature_t1 =
     $      norm51_0*normt01_0+norm51_1*normt01_1+norm51_2*normt01_2

          gauss  = gauss  + norm_do
          gausst = gausst + norm_dt
          edge  = edge  + curvature + curvature1
          edget = edget + curvature_t + curvature_t1

          m_start = m_next_m(m_start)
          icount = icount + 1

          GOTO 12
        ENDIF

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

        IF(new_action .GE. old_action)THEN
          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

          IF(exp(-(new_action-old_action)).lt.rand_uni)goto22
        ENDIF

        sw_sacc = sw_sacc + 1.0
        field(0,s0) = ft0
        field(1,s0) = ft1
        field(2,s0) = ft2

        DO i = 0, icount - 1
          mem = iv_scra(i)
          s1 = s_neighbor_m(mem)
          s3 = s_right_m(mem)

          square_m(mem) = v_scra(i)

          rno_m(0,mem)=rno_scra(0,i)
          rno_m(1,mem)=rno_scra(1,i)
          rno_m(2,mem)=rno_scra(2,i)

          lno_m(0,mem)=lno_scra(0,i)
          lno_m(1,mem)=lno_scra(1,i)
          lno_m(2,mem)=lno_scra(2,i)

          m_start = m_back_m(mem)

          square_m(m_start)=v_scra(i)

          rno_m(0,m_start)=lno_scra(0,i)
          rno_m(1,m_start)=lno_scra(1,i)
          rno_m(2,m_start)=lno_scra(2,i)

          lno_m(0,m_start)=rno_scra(0,i)
          lno_m(1,m_start)=rno_scra(1,i)
          lno_m(2,m_start)=rno_scra(2,i)

          m_start1 = m_scra(i)

          lno_m(0,m_start1) = fno_scra(0,i)
          lno_m(1,m_start1) = fno_scra(1,i)
          lno_m(2,m_start1) = fno_scra(2,i)

          m_start2 = m_back_m(m_start1)

          rno_m(0,m_start2) = fno_scra(0,i)
          rno_m(1,m_start2) = fno_scra(1,i)
          rno_m(2,m_start2) = fno_scra(2,i)

        ENDDO

   22   CONTINUE
      ENDDO
      sw_sacc = sw_sacc / (3.*n)
      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_end_tot, m_start2, s1, s2, counter

      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

      m_end_tot = 0
      DO s0 = 0, nm
        m_start = m_start_s(s0)
        m_end   = 0
   11   CONTINUE
        IF(m_start .GE. 0)THEN
          m_end = m_end+1
          m_start = m_next_m(m_start)
          GOTO 11
        ENDIF
        m_end_tot = m_end_tot + m_end
        IF(n_connectivity_s(s0).NE.m_end)THEN
          PRINT*,'TERROR!!!! NEXT, s0, conn,m_end = '
     $           ,s0,n_connectivity_s(s0),m_end
          stop'CHECK_MESH TERROR ABORT: n_conn .ne. m_end'
        ENDIF
      ENDDO
      IF(m_end_tot.NE.n6)THEN
        PRINT*,'TERROR!!!! NEXT, m_end_tot = ',m_end_tot
        stop'CHECK_MESH TERROR ABORT: m_end_tot .ne. n6'
      ENDIF

      m_end_tot = 0
      DO s0 = 0, nm
        m_start = m_end_s(s0)
        m_end   = 0
   12   CONTINUE
        IF(m_start .GE. 0)THEN
          m_end = m_end+1
          m_start = m_former_m(m_start)
          GOTO 12
        ENDIF
        m_end_tot = m_end_tot + m_end
        IF(n_connectivity_s(s0).NE.m_end)THEN
          PRINT*,'TERROR!!!! FORMER, s0, conn,m_end = '
     $           ,s0, n_connectivity_s(s0),m_end
          stop'CHECK_MESH TERROR ABORT: n_conn .ne. m_end'
        ENDIF
      ENDDO
      IF(m_end_tot.NE.n6)THEN
        PRINT*,'TERROR!!!! FORMER, m_end_tot = ',m_end_tot
        stop'CHECK_MESH TERROR ABORT: m_end_tot .ne. n6'
      ENDIF
      DO s0 = 0, nm
        IF(n_connectivity_s(s0).LT.3)THEN
          PRINT*,'CHECK_MESH TERROR!!!! n_connectivity_s(s0).lt.3'
          PRINT*,'s0=',s0,' n_connectivity_s(s0)=',
     $           n_connectivity_s(s0)
          stop'TERROR ABORT'
        ENDIF
      ENDDO

      DO s0=0, nm
        m_start = m_start_s(s0)
        m_end   = 0
   41   CONTINUE
        IF(m_start .GE. 0)THEN
          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
C
          m_start2 = m_start_s(s2)
   51     CONTINUE
          IF(m_start2 .GE. 0)THEN
            IF(s1.EQ.s_neighbor_m(m_start2))GOTO 52
            m_start2 = m_next_m(m_start2)
            GOTO 51
          ENDIF
          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
C
          counter = 0
          m_start1 = m_start_s(s1)
   71     CONTINUE
          IF(m_start1 .GE. 0)THEN
            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)
            GOTO 71
          ENDIF
          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
C
          m_end = m_end+1
          m_start = m_next_m(m_start)
          GOTO 41
        ENDIF
      ENDDO
C      print*, 'CHECK MESH: OK _ ', i_sweep

      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 **************************************************************************

