[tsscds] Diff of /trunk/src/kmc.f90
Diff of /trunk/src/kmc.f90
 Parent Directory
|  Revision Log
|  Patch
| 8 |
! ndd species that dissapear |
! ndd species that dissapear |
| 9 |
character*80 :: title |
character*80 :: title |
| 10 |
integer, parameter :: dp = selected_real_kind(15,307) |
integer, parameter :: dp = selected_real_kind(15,307) |
| 11 |
integer :: m,nran,nr,nesp,inran,i,j,k,l,mu,ndisp,pd,kk,ijk,tcont,tint,rep |
integer :: m,nran,nr,nesp,inran,i,j,k,l,mu,ndisp,pd,kk,ijk,tcont,tint,lexit |
| 12 |
real(dp) :: t,a0,r2a0,suma,ptot,diff |
real(dp) :: t,a0,r2a0,suma,ptot |
| 13 |
real(dp) :: rnd(2) |
real(dp) :: rnd(2) |
| 14 |
integer,dimension(:),allocatable :: p,p0,re,pr,n,ndd,cont,pp,ppold |
integer,dimension(:),allocatable :: p,p0,re,pr,n,ndd,cont,pold,sig,sigold,rep |
| 15 |
real (dp) ,dimension(:),allocatable :: a,rate |
real (dp) ,dimension(:),allocatable :: a,rate |
| 16 |
read(*,"(a80)") title |
read(*,"(a80)") title |
| 17 |
print "(t3,a80)",title |
print "(t3,a80)",title |
| 18 |
read(*,*) m,nesp,nran |
read(*,*) m,nesp,nran |
| 19 |
allocate(re(m),pr(m),a(m),rate(m),p0(nesp),p(nesp),n(nesp),cont(m),pp(nesp),ppold(nesp)) |
allocate(re(m),pr(m),a(m),rate(m),p0(nesp),p(nesp),n(nesp),cont(m),pold(nesp),sig(nesp),sigold(nesp),rep(nesp)) |
| 20 |
n=(/ (l,l=1,nesp) /) |
n=(/ (l,l=1,nesp) /) |
| 21 |
do i=1,m |
do i=1,m |
| 22 |
read(*,*) rate(I),re(i),pr(i) |
read(*,*) rate(I),re(i),pr(i) |
| 43 |
big: do inran=1,nran |
big: do inran=1,nran |
| 44 |
print "(/,t3,a,i4,/,t3,a,9999(i7))","Calculation number",inran,"Time(ps)",n |
print "(/,t3,a,i4,/,t3,a,9999(i7))","Calculation number",inran,"Time(ps)",n |
| 45 |
p=p0 |
p=p0 |
| 46 |
pp=p0 |
pold=p0 |
| 47 |
|
sig=0 |
| 48 |
t=0.d0 |
t=0.d0 |
| 49 |
tcont=0 |
tcont=0 |
| 50 |
|
lexit=0 |
| 51 |
rep=0 |
rep=0 |
| 52 |
print "(e10.4,9999(i7))",t,p |
print "(e10.4,9999(i7))",t,p |
| 53 |
do j=1,m |
do j=1,m |
| 59 |
t=t-log(rnd(1))/a0 |
t=t-log(rnd(1))/a0 |
| 60 |
r2a0=rnd(2)*a0 |
r2a0=rnd(2)*a0 |
| 61 |
suma=0.d0 |
suma=0.d0 |
| 62 |
|
pold=p |
| 63 |
s1: do mu=1,m |
s1: do mu=1,m |
| 64 |
suma=suma+a(mu) |
suma=suma+a(mu) |
| 65 |
if(suma>=r2a0) exit s1 |
if(suma>=r2a0) exit s1 |
| 66 |
enddo s1 |
enddo s1 |
| 67 |
p(re(mu))=p(re(mu))-1 |
p(re(mu))=p(re(mu))-1 |
| 68 |
p(pr(mu))=p(pr(mu))+1 |
p(pr(mu))=p(pr(mu))+1 |
| 69 |
|
|
| 70 |
|
do j=1,m |
| 71 |
|
sigold(j)=sig(j) |
| 72 |
|
if(p(j)>pold(j)) sig(j)=0 |
| 73 |
|
if(p(j)<pold(j)) sig(j)=1 |
| 74 |
|
if(sig(j).ne.sigold(j).and.p(j)>0.and.pold(j)>0) then |
| 75 |
|
rep(j)=rep(j)+1 |
| 76 |
|
if(rep(j)==1000) lexit=1 |
| 77 |
|
endif |
| 78 |
|
enddo |
| 79 |
|
if(lexit==1) exit |
| 80 |
cont(mu)=cont(mu)+1 |
cont(mu)=cont(mu)+1 |
| 81 |
tcont=tcont+1 |
tcont=tcont+1 |
| 82 |
do j=1,m |
do j=1,m |
| 83 |
a(j)=rate(j)*p(re(j)) |
a(j)=rate(j)*p(re(j)) |
| 84 |
enddo |
enddo |
| 85 |
a0=sum(a) |
a0=sum(a) |
| 86 |
if(mod(tcont,tint)==0.or.a0==0) then |
if(mod(tcont,tint)==0.or.a0==0) print "(e10.4,9999(i7))",t,p |
|
print "(e10.4,9999(i7))",t,p |
|
|
ppold=pp |
|
|
pp=p |
|
|
diff=0 |
|
|
do j=1,m |
|
|
diff=diff+(ppold(j)-pp(j))**2 |
|
|
enddo |
|
|
if(diff==0) rep=rep+1 |
|
|
if(rep==1000) exit |
|
|
endif |
|
| 87 |
enddo |
enddo |
| 88 |
enddo big |
enddo big |
| 89 |
print*,"Population of every species" |
print*,"Population of every species" |
|
|
Legend:
| Removed from v.287 |
|
| changed lines |
| |
Added in v.288 |
|
|
|
