!L=12;b=0.2;for n in 1 2 3 4;do echo "$L ${b}$n 1000 1" > if00$n;done
!caf IsingEmbarrassing.f90 -o is; cafrun -n 4 ./is
!
!If you need more processes than the cores of your CPU, use:
!
!cafrun -n 20 -oversubscribe ./is
!
!CAF=/usr/lib/x86_64-linux-gnu/open-coarrays/openmpi/lib
!mpifort -fcoarray=lib IsingEmbarrassing.f90 -L $CAF -lcaf_mpi -o is ; mpirun -n 4 ./is
!
!================== IsingEmbarrassing.f90 =======
!--------------------------------------------------
module          IsingEmbarrassing
 use, intrinsic          :: iso_fortran_env
 implicit none
 SAVE
 integer, parameter      :: dp = real64
 integer                 :: L
 integer                 :: N
 integer                 :: XNN, YNN
 integer ,allocatable    :: s(:)
 real(dp),dimension(0:4) :: prob
 real(dp)                :: beta
 integer                 :: nsweep,start
 real(dp)                :: r
 integer                 :: fin, fout, status
 integer                 :: nimg, myimg !, myi, myj
 character(1000)         :: fname_in,fname_out
!--------------------------------------------------
contains
!--------------------------------------------------
 subroutine     init !inherits implicit none!
  integer :: i
!----------------------
! Images:
  nimg  = num_images()
  myimg = this_image()
  write(fname_in ,'(A,I0.3)') 'if',myimg
  write(fname_out,'(A,I0.3)') 'of',myimg
  open(newunit=fin ,file=fname_in ,status='old',iostat=status)
  if(status /= 0 ) then
   write(error_unit,*) 'Error opening file ', trim(fname_in)
   error stop
  end if
  open(newunit=fout,file=fname_out)
!----------------------
!Ask user for input:
  !#               Enter L, beta, nsweep, start (0 cold/1 hot):'
  read (fin , *)         L, beta, nsweep, start
  write(fout, *)'#       L  beta  nsweep  start  myimg  nimg     (start: 0 cold/1 hot)'
  write(fout, *)'#     ',L, beta, nsweep, start, myimg, nimg
  N = L*L; XNN = 1; YNN = L
  ALLOCATE(s(N))

  call random_init(.false.,.true.) ! random_init(REPEATABLE, IMAGE_DISTINCT)

!initialize probabilities for E_\nu > E_mu
  prob=0.0_dp
! Sum_<ik> s_i^\mu = -4, -2, 0, +2, +4
! Need:     exp(-2*beta*4), exp(- 2*beta*2) 
! and not:  exp( 2*beta*0), exp(  2*beta*2),exp(2*beta*4)
  do i=2,4,2 !i = dE/2 = (E_nu-E_mu)/2=2,4
   prob(i) = exp(-2.0_dp*beta*i)
  enddo
  !initial configuration:
  select case(start)
  case(0)!cold:
   s = 1 !all s(i) = 1
  case(1)!hot:
   do i=1,N
    call random_number(r)
    if(r .lt. 0.5_dp)then
     s(i) =  1
    else
     s(i) = -1
    endif
   enddo
  case default
   print *,'init: start= ',start,' not valid. Exiting...',myimg
   error stop
  end select
  
 end subroutine init
!--------------------------------------------------
subroutine      met()
 integer :: i,k
 integer :: nn,snn,dE

 do k=1,N
!pick a random site:
  call random_number(r)
  i = INT(N*r)+1
!snn=sum of neighboring spins:
  snn = 0
  nn=i+XNN;if(nn.gt.N)nn=nn-N;snn = snn + s(nn)
  nn=i-XNN;if(nn.lt.1)nn=nn+N;snn = snn + s(nn)
  nn=i+YNN;if(nn.gt.N)nn=nn-N;snn = snn + s(nn)
  nn=i-YNN;if(nn.lt.1)nn=nn+N;snn = snn + s(nn)
!dE=change in energy/2:
  dE=snn*s(i)
!flip:
  if(dE.le.0)then
   s(i) = -s(i) !accept
  else 
   call random_number(r)
   if(r < prob(dE)) s(i) = -s(i) !accept
  end if
 enddo !do k=1,N: end sweep
end subroutine  met
!--------------------------------------------------
subroutine      measure()
 write(fout,*) E(),M()
end subroutine  measure
!--------------------------------------------------
function        E()
 integer :: E
 integer :: en,sum,i,nn
 en = 0
 do i=1,N
!Sum of neighboring spins: only forward nn necessary in the sum
  sum = 0
  nn=i+XNN;if(nn.gt.N)nn=nn-N;sum = sum + s(nn)
  nn=i+YNN;if(nn.gt.N)nn=nn-N;sum = sum + s(nn)
  en=en+sum*s(i)
 enddo
 e = -en
end function    E
!--------------------------------------------------
function        M()
 integer :: M
 M=SUM(s)
end function    M
!--------------------------------------------------
end module      IsingEmbarrassing
!==================-------------===================
program         ising2dProgram
 use IsingEmbarrassing
 implicit none
 integer isweep
 call init
 do isweep = 1, nsweep
  call met
  call measure
 end do
end program ising2dProgram
!=================================================
!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2022)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------