! APPENDIX: Code for Simple Adsoprtion-Desorption Dynamic Monte Carlo
! created by Nikolai (23 May 2001) for D.D. Johnson Lectures
! see: K.A.Fichthorn, W.H.Weinberg, J. Chem. Phys. 95 (2) p.1090 (1991).
!!!!!^!!!1!!!!^!!!!2!!!!^!!!!3!!!!^!!!!4!!!!^!!!!5!!!!^!!!!6!!!!^!!!!7!!
! Required Structure of INPUT DATA FILE:
!!  128  128   ! dimensions of simulation box (<Lx=256)
!!  2.5        ! max real time (secs) **use 2.5 to start**
!!  1.0 2.0    ! rate parameters (rA,rD)
!!!!!^!!!1!!!!^!!!!2!!!!^!!!!3!!!!^!!!!4!!!!^!!!!5!!!!^!!!!6!!!!^!!!!7!!
      program dMCs   
      implicit none

!glsc:
      include '../glsc-3.5/libs/src/glsc_ftn.h' !!!check_the_path!!!
      
      integer Lx, Ly, LT, N
      parameter( Lx=512, Ly=512, LT= 100000, N=10 )
      integer isite(Lx,Ly)
      real*8  WA, WD, rA, rD   
      real*8  r !random_number
      real*8  t,dt,tl,tt !real_time
      !theta = fractional surface covarage: 0<theta<1 !
      real*8  rc,tc,theta,theta0,xc, exptt,expt1
      integer iT
      integer mx,my
      integer i,j,ic,nc 
      integer Input, NoTheory
      real*8 pt !length of plot in x-direction = max real time!
      
      
  ! initial default parameters:
      NoTheory=0 !do not draw theoretical curve!
      Input=0    !do not read Input!
      mx=Lx
      my=Ly
      pt=LT
      rA=1.0
      rD=2.0
  ! initial TIME and Fraction Coverage:
      iT=0    !integer time 
      t=0     !real time
      tl=t
      dt=0.001  !real time step
      theta=0 !fractional surface coverage, 0 <= theta <= 1 .   
      
      

  !here read parameters...: mx,my,pt,rA,rD:
      !.........................input............................!
      Input=1
      NoTheory=1
      if(Input.ne.0) then
        read(5,*) mx,my !dimensions of simulation box
        read(5,*) pt    !max real time
        read(5,*) rA,rD !parameters
      endif
      if(mx>Lx) then
        mx=Lx 
        print*, "Input Error: reducing Lx!"
      endif
      if(my>Ly) then
        my=Ly 
        print*, "Input Error: reducing Ly!"
      endif
      
  !parameters mx,my,pt,rA,rD remain unchanged after this line!
      rc=theta
      nc=mx*my
      xc=1/dfloat(nc)
      if(rA > rD) then 
        WA=1.0
        WD=rD/rA
      else
        WA=rA/rD
        WD=1.0
      endif
      if(rc<0.or.rc>1) stop "Check: should be 0 < rc < 1 !!!"
      theta0=rA/(rA+rD)
      
  !initialize lattice at random:
      ic=0
      call RANDOM_SEED
      do i=1,mx
       do j=1,my
        call RANDOM_NUMBER(r)
        if(r<rc) then
          isite(i,j)=1
          ic=ic+1
        else 
          isite(i,j)=0
        endif
       enddo
      enddo

!glsc:   !Initialization and definitions
	call g_init ('Meta', 4, 250.0d0, 160.0d0)
	call g_device (G_BOTH)
          call g_def_scale (5, 0.0d0, pt, 0.0d0, theta0+0.05d0,
     5                         10.0d0, 10.0d0, 220.0d0, 120.0d0)
!glsc:    g_def_scale(#number, ?,?,?,?, x0,y0,x_width,y_width)
	call g_def_line (0, G_BLACK, 0, G_LINE_SOLID)
	call g_def_line (1, G_RED, 0, G_LINE_SOLID)
	call g_def_line (2, G_GREEN, 0, G_LINE_DASHED)
	call g_def_line (3, G_BLUE, 0, G_LINE_SOLID)
	call g_def_text (0, G_BLACK, G_FONT_TIMES_12)
	call g_def_text (1, G_RED, G_FONT_TIMES_18)
	call g_def_text (2, G_GREEN, G_FONT_TIMES_18)
	call g_def_text (3, G_BLUE, G_FONT_TIMES_24)	

!glsc:Plot graphics:
      call g_cls () !clear_screen!
      call g_sel_scale (5)

      !draw box:
      call g_sel_line (0)
      call g_box (0.0d0, pt, 0.0d0, theta0+0.05d0, G_YES, G_NO) 
      
  !draw text:
      call g_sel_text (1)
      call g_text (110.0d0, 20.0d0,'theta(t)',8)
      call g_sel_text (0)
      call g_text (110.0d0, 140.0d0,'real time t',14)


  !draw line y=theta0:
      call g_sel_line (2)
      call g_move (0.0d0, theta0)
      call g_plot (pt,   theta0) !draw line y=theta0!
      !call g_text (0.0d0, theta0,'theta0',7)
      
      if(NoTheory.ne.0) then
      !draw theoretical line theta(t)
      call g_sel_line (3)
      call g_move (0.0d0, theta)
      expt1=exp(-(rA+rD)*dt) !/dfloat(mx*my))
      exptt=1
      tt=0
      do while(tt<pt)
        tc=theta0*(1-exptt)
        tt=tt+dt
        exptt=exptt*expt1
        call g_plot (tt,tc)
        call g_move (tt,tc)
      enddo
      endif
      
      !draw line theta(t)
      call g_sel_line (1)
      
  !!!run Monte Carlo: Poisson process:
      do while(t<pt)
        !Select a random site (i,j):
        call RANDOM_NUMBER(r)
        i=int(r*mx)+1
        call RANDOM_NUMBER(r)
        j=int(r*my)+1
        
   !Generate a random number 0<r<1:
        call RANDOM_NUMBER(r)
        
        !Occupied?
        if(isite(i,j).ne.0) then 
        !YES
          if(r.le.WD) then
            !Remove species from lattice
            isite(i,j)=0
            ic=ic-1
      ! TIME from Poisson Distribution
            call RANDOM_NUMBER(r)
            dt=-log(r)/(ic*rD+(nc-ic)*rA)
            t=t+dt
          endif
        else
        !NOT Occupied
          if(r.le.WA) then
            !Add species to lattice
            isite(i,j)=1
            ic=ic+1
      ! TIME from Poisson Distribution
            call RANDOM_NUMBER(r)
            dt=-log(r)/(ic*rD+(nc-ic)*rA)
            t=t+dt
          endif
        endif
        iT=iT+1
        rc=ic*xc

!glsc: !Plot_graphic: theta(t):
        call g_move (tl,theta)
        theta=rc
        tl=t
        call g_plot (tl,theta) 
                
        !print*, t,rc
        
      enddo !while(t<pt)

!glsc:Termination:
      pause
      !call g_sleep (3.0d0)
      call g_term ()
      
      end
!end_of_file!