program FortranFigure4_1 implicit none save integer, parameter :: long = selected_real_kind(15,50) integer, parameter :: numbermesh=40 integer, parameter :: nahbon=100000 real(kind=long), dimension(nahbon) :: RanNum real(kind=long), dimension(numbermesh) :: A,k2S,DhS integer, dimension(numbermesh,2) :: Asave real(kind=long), parameter :: DiffConst=0.0001_long ! diffusion constant real(kind=long), parameter :: k1=0.001_long ! degradation rate real(kind=long), parameter :: k2=0.3_long ! creation rate real(kind=long) :: time,Dh,h,ss,a0,savingtime,finaltime,ate,tau integer :: ii,ir,ns logical :: savingdata integer size, seed(2) data seed /10,10/ size = 2 h=0.025_long ! h=1/real(numbermesh) Dh=DiffConst/(h*h) time=0_long ir=0 a0=0.0_long savingtime=10.0_long*60.0_long finaltime=30.0_long*60.0_long savingdata=.TRUE. call random_seed(SIZE=size) call random_seed(PUT=seed(1:size)) call random_number(RanNum) do ii=1,numbermesh A(ii)=0.0_long end do do ii=2,numbermesh-1 DhS(ii)=2.0_long*Dh end do ii=1 DhS(ii)=Dh ii=numbermesh DhS(ii)=Dh do ii=1,numbermesh/5 k2S(ii)=k2 end do do ii=numbermesh/5+1,numbermesh k2S(ii)=0.0_long end do do ii=1,numbermesh a0=a0+DhS(ii)*A(ii)+k1*A(ii)+k2S(ii) end do do while (timenahbon) then call random_number(RanNum) ir=1 end if tau=(1/a0)*log(1/RanNum(ir)) ir=ir+1 if (ir>nahbon) then call random_number(RanNum) ir=1 end if ss=0_long ii=0 do while ((ssRanNum(ir)*a0) then A(ii)=A(ii)-1_long a0=a0-k1-DhS(ii) else ii=0 do while ((ssRanNum(ir)*a0) then A(ii)=A(ii)+1_long a0=a0+k1+DhS(ii) else ii=0 do while ((ssRanNum(ir)*a0) then A(ii)=A(ii)-1_long A(ii+1)=A(ii+1)+1_long a0=a0+DhS(ii+1)-DhS(ii) else ii=1 do while ((ssRanNum(ir)*a0) then A(ii)=A(ii)-1_long A(ii-1)=A(ii-1)+1_long a0=a0+DhS(ii-1)-DhS(ii) else write(*,*) 'Error.' STOP end if end if end if end if time=time+tau if (int(0.01_long*time)-int(0.01_long*(time-tau))>0.5_long) then ate=0_long do ii=1,numbermesh ate=ate+DhS(ii)*A(ii)+k1*A(ii)+k2S(ii) end do write(*,*) 'Time computed:',nint(time),'seconds (final time 1800 sec)' end if if ((time>savingtime).AND.savingdata) then ns=1 savingdata=.FALSE. do ii=1,numbermesh Asave(ii,ns)=nint(A(ii)) end do end if end do ns=2 do ii=1,numbermesh Asave(ii,ns)=nint(A(ii)) end do open(unit=1, file='dataFigure4_1.dat', status='replace') 100 format(i6,X,i6) do ii=1,numbermesh write(1,100) Asave(ii,1),Asave(ii,2) end do close(unit=1) end program FortranFigure4_1