C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C PROGRAM DIFF_TH C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C REAL*8 Ki,Kd,KL,kk INTEGER timtyp C PARAMETER (PI=3.14159265358979324D0) PARAMETER (NM=25,NMT=2*NM,NME=(NM+1)*(NM+1),NMQ=NM*(NM+1)) C COMPLEX*16 SST,ZZ,ZZP,HHH,FHH,RHH,RATIO COMPLEX*16 beta,gama,rm,pp1,pp2,pp3 COMPLEX*16 v,u,DD,R,gg1,gg2,gg3,eta1,eta2 COMPLEX*16 Psi0(0:NMT),PsiX(0:NMT) COMPLEX*16 rms(0:NMT),rmw(0:NMT) C COMPLEX*16 EH1(NME),QH1(NMQ) COMPLEX*16 DH1(0:NMT),AAH1(-1:NMT),BBH1(-1:NMT) COMPLEX*16 EH2(NME),QH2(NMQ) COMPLEX*16 DH2(0:NMT),AAH2(-1:NMT),BBH2(-1:NMT) COMPLEX*16 EH3(NME),QH3(NMQ) COMPLEX*16 DH3(0:NMT),AAH3(-1:NMT),BBH3(-1:NMT) COMPLEX*16 EH4(NME),QH4(NMQ) COMPLEX*16 DH4(0:NMT),AAH4(-1:NMT),BBH4(-1:NMT) C DIMENSION LEN(0:NMT),LENE(-1:NM),LENQ(-1:NM) DIMENSION TIME(0:500) DIMENSION Ctt(0:1001) C CHARACTER*20 DATANM CHARACTER*25 OOUUT1,OOUUT2,OOUUT3 C C C =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of DIFF_TH C C C ********************************************************************** C * DATA INPUTS * C ********************************************************************** C C WRITE (6,6000) READ (5,5000) DATANM DO 1 I=1,20 IF(DATANM(I:I).EQ.' ') GO TO 2 1 CONTINUE 2 IF(I.NE.20) THEN IM = I-1 ELSE IM = I END IF C OOUUT1 = DATANM(1:IM) //'.C(t)' OOUUT2 = DATANM(1:IM) //'.M(t)' OOUUT3 = DATANM(1:IM) //'.Hist' C OPEN(UNIT=1,FILE=DATANM,STATUS='OLD') OPEN(UNIT=2,FILE=OOUUT1) OPEN(UNIT=3,FILE=OOUUT2) OPEN(UNIT=4,FILE=OOUUT3) C READ(1,*) itype,nequ IF(itype.NE.1.AND.itype.NE.2.AND.itype.NE.3) THEN WRITE(2,6001) STOP END IF IF(nequ.NE.1.AND.nequ.NE.2) THEN WRITE(2,6002) STOP END IF C READ(1,*) phi,rho,Sr READ(1,*) D0,Ds,decay READ(1,*) tortp,torti,torts READ(1,*) Kd,Ki,KL,kk C IF(Sr.EQ.0.0d0) Ki = 1.0d0 IF(Ki.EQ.1.0d0) Sr = 0.0d0 C READ(1,*) Vup,Vdn,A,X READ(1,*) M,FACTOR,ERRMAX READ(1,*) NPERIOD,timtyp C IF(timtyp.NE.1.AND.timtyp.NE.2) THEN WRITE(2,6003) STOP END IF C ntime = 0 TIME(0) = 0.0d0 IF(timtyp.EQ.1) THEN DO n8=1,NPERIOD ntime = n8 READ(1,*) TIME(n8) END DO ELSE DO 3 n=1,NPERIOD READ(1,*) NPOINT,TINCR C DO 3 I=1,NPOINT ntime = ntime+1 TIME(ntime) = TIME(ntime-1)+TINCR 3 CONTINUE END IF C C ********************************************************************** C * BASIC PARAMETERS & ARRAY DIMENSIONING * C ********************************************************************** C LENE(-1) = 0 LENQ(-1) = 0 DO 4 I=0,NMT LEN(I) = 0 4 CONTINUE DO 5 I=0,NM LENE(I) = 0 LENQ(I) = 0 5 CONTINUE C M2 = 2*M LEN(0) = M2 LEN(1) = M2 DO 6 I=2,M2 LEN(I) = LEN(I-1)-1 6 CONTINUE C NSUME = 0 NSUMQ = 0 DO 7 I=0,M2 IF(MOD(I,2).EQ.0) THEN IE = I/2 NSUME = NSUME+LEN(I) LENE(IE) = NSUME ELSE IQ = I/2+1 NSUMQ = NSUMQ+LEN(I) LENQ(IQ) = NSUMQ END IF 7 CONTINUE C h = 1.0d0-Sr+Sr*Ki w0 = ((1.0d0-phi)/phi)*rho*Ki w = w0*Kd DT = D0*(tortp*(1.0d0-Sr)+torti*Sr*Ki) C Cup = 1.0d0 Cdn = 0.0d0 Clp = 0.0d0 Cratio = 0.0d0 WRITE(2,6005) TIME(0),Cup,Cdn,Clp,Cratio C Ctt(0) = 0.0d0 Ctt(2*ntime+1) = 0.0d0 C rmu = Cup*Vup rmww = 0.0d0 rmss = 0.0d0 rmd = 0.0d0 rmt = rmu C WRITE(3,6005) TIME(0),rmu,rmww,rmss,rmd,rmt C C C C ********************************************************************** C * DIFFUSION WITH LINEAR EQUILIBRIUM SORPTION * C ********************************************************************** C IF(itype.EQ.1) THEN KL = 0.0d0 kk = 0.0d0 C C ********************************************************************** C * DIFFUSION WITH KINETIC LINEAR SORPTION * C ********************************************************************** C ELSE IF(itype.EQ.2) THEN KL = 0.0d0 u0 = w*kk C C ********************************************************************** C * DIFFUSION WITH IRREVERSIBLE LINEAR SORPTION * C ********************************************************************** C ELSE kk = 0.0d0 v0 = w0*KL END IF C C C C C ********************************************************************** C ********************************************************************** C ********************************************************************** C * THE TIME LOOP * C ********************************************************************** C ********************************************************************** C ********************************************************************** C C C C DO 1000 NUU=1,NTIME C TPRIOD = FACTOR*TIME(nuu) GAMMA = -DLOG(ERRMAX)/(2.0D0*TPRIOD) C Sfact = (1.0D0/TPRIOD)*DEXP(GAMMA*TIME(nuu)) C DO 100 I=0,M2 ARGG = I*PI/TPRIOD SST = DCMPLX(GAMMA,ARGG) C C ********************************************************************** C * DIFFUSION WITH LINEAR EQUILIBRIUM SORPTION * C ********************************************************************** C IF(itype.EQ.1) THEN R = DCMPLX(h+w,0.0d0) DD = DCMPLX((DT+Ds*torts*w),0.0d0) C C ********************************************************************** C * DIFFUSION WITH KINETIC LINEAR SORPTION * C ********************************************************************** C ELSE IF(itype.EQ.2) THEN u = u0/(sst+kk+decay) R = h+u DD = DT+Ds*torts*u C C ********************************************************************** C * DIFFUSION WITH IRREVERSIBLE LINEAR SORPTION * C ********************************************************************** C ELSE v = v0/(sst+decay) R = h+v DD = DT+Ds*torts*v END IF C rm = cdsqrt(R*(sst+decay)/DD) C gg1 = phi*DD*A*rm gg2 = Vdn*(sst+decay) IF(Vdn.NE.0.0d0) THEN gama = (gg1-gg2)/(gg1+gg2) ELSE gama = DCMPLX(1.0d0,0.0D0) END IF C IF(nequ.EQ.1) THEN pp1 = Vup*(1.0d0+gama*exp(-2.0d0*rm*X)) pp2 = (phi*A*R/rm)*(1.0d0-gama*exp(-2.0d0*rm*X)) pp3 = (Vdn*(1.0d0+gama)-(phi*A*R/rm)*(1.0d0-gama)) & *exp(-rm*X) beta = Vup/((sst+decay)*(pp1+pp2+pp3)) ELSE gg3 = Vup*(sst+decay) eta1 = gg3-gg1 eta2 = gg3+gg1 beta = Vup/(eta2+gama*eta1*exp(-2.0d0*rm*X)) END IF C Psi0(I) = beta*(1.0d0+gama*exp(-2.0d0*rm*X)) PsiX(I) = beta*(1.0d0+gama)*exp(-rm*X) rmw(I) = (beta/rm)*(1.0d0-exp(-rm*X)) & *(1.0d0+gama*exp(-rm*X)) C IF(itype.EQ.1) THEN rms(I) = rmw(I)*w ELSE IF(itype.EQ.2) THEN rms(I) = rmw(I)*u ELSE rms(I) = rmw(I)*v END IF C IF(I.EQ.0) THEN Psi0(0) = 5.0d-1*Psi0(0) PsiX(0) = 5.0d-1*PsiX(0) rmw(0) = 5.0d-1*rmw(0) rms(0) = 5.0d-1*rms(0) END IF C 100 CONTINUE C C ********************************************************************** C * EVALUATION OF EH, QH, DH, AAH, BBH * C ********************************************************************** C DO 15 I=1,M2 EH1(I) = (0.0D0,0.0D0) EH2(I) = (0.0D0,0.0D0) EH3(I) = (0.0D0,0.0D0) EH4(I) = (0.0D0,0.0D0) 15 CONTINUE C DO 20 I=1,M2 I1 = I-1 QH1(I) = Psi0(I)/Psi0(I1) QH2(I) = PsiX(I)/(PsiX(I1)+1.0d-35) QH3(I) = rmw(I)/(rmw(I1)+1.0d-35) QH4(I) = rms(I)/(rms(I1)+1.0d-35) 20 CONTINUE C DO 40 J=1,M DO 25 I1=0,M2-2*J I = I1+1 NN = LENE(J-1)+I NN1 = LENQ(J-1)+I NN2 = LENE(J-2)+I C EH1(NN) = QH1(NN1+1)-QH1(NN1)+EH1(NN2+1) EH2(NN) = QH2(NN1+1)-QH2(NN1)+EH2(NN2+1) EH3(NN) = QH3(NN1+1)-QH3(NN1)+EH3(NN2+1) EH4(NN) = QH4(NN1+1)-QH4(NN1)+EH4(NN2+1) 25 CONTINUE C DO 30 I1=0,M2-2*J-1 I = I1+1 NN = LENQ(J)+I NN1 = LENE(J-1)+I NN2 = LENQ(J-1)+I C QH1(NN) = QH1(NN2+1)*EH1(NN1+1)/EH1(NN1) QH2(NN) = QH2(NN2+1)*EH2(NN1+1)/(EH2(NN1)+1.0d-35) QH3(NN) = QH3(NN2+1)*EH3(NN1+1)/(EH3(NN1)+1.0d-35) QH4(NN) = QH4(NN2+1)*EH4(NN1+1)/(EH4(NN1)+1.0d-35) 30 CONTINUE 40 CONTINUE C DH1(0) = Psi0(0) DH2(0) = PsiX(0) DH3(0) = rmw(0) DH4(0) = rms(0) C DO 50 I=1,M I1 = 2*I-1 I2 = 2*I NN1 = LENQ(I-1)+1 NN2 = LENE(I-1)+1 C DH1(I1) = -QH1(NN1) DH1(I2) = -EH1(NN2) C DH2(I1) = -QH2(NN1) DH2(I2) = -EH2(NN2) C DH3(I1) = -QH3(NN1) DH3(I2) = -EH3(NN2) C DH4(I1) = -QH4(NN1) DH4(I2) = -EH4(NN2) C 50 CONTINUE C C C AAH1(-1) = (0.0D0,0.0D0) BBH1(-1) = (1.0D0,0.0D0) AAH1(0) = DH1(0) BBH1(0) = (1.0D0,0.0D0) C AAH2(-1) = (0.0D0,0.0D0) BBH2(-1) = (1.0D0,0.0D0) AAH2(0) = DH2(0) BBH2(0) = (1.0D0,0.0D0) C AAH3(-1) = (0.0D0,0.0D0) BBH3(-1) = (1.0D0,0.0D0) AAH3(0) = DH3(0) BBH3(0) = (1.0D0,0.0D0) C AAH4(-1) = (0.0D0,0.0D0) BBH4(-1) = (1.0D0,0.0D0) AAH4(0) = DH4(0) BBH4(0) = (1.0D0,0.0D0) C ARGG = PI*TIME(nuu)/TPRIOD AR1 = DCOS(ARGG) AR2 = DSIN(ARGG) ZZ = DCMPLX(AR1,AR2) C DO 60 I=1,M2-1 IM1 = I-1 IM2 = I-2 C ZZP = DH1(I)*ZZ AAH1(I) = AAH1(IM1)+ZZP*AAH1(IM2) BBH1(I) = BBH1(IM1)+ZZP*BBH1(IM2) C ZZP = DH2(I)*ZZ AAH2(I) = AAH2(IM1)+ZZP*AAH2(IM2) BBH2(I) = BBH2(IM1)+ZZP*BBH2(IM2) C ZZP = DH3(I)*ZZ AAH3(I) = AAH3(IM1)+ZZP*AAH3(IM2) BBH3(I) = BBH3(IM1)+ZZP*BBH3(IM2) C ZZP = DH4(I)*ZZ AAH4(I) = AAH4(IM1)+ZZP*AAH4(IM2) BBH4(I) = BBH4(IM1)+ZZP*BBH4(IM2) C 60 CONTINUE C C C HHH = 5.0D-1*(1.0D0+ZZ*(DH1(M2-1)-DH1(M2))) FHH = 1.0D0+(ZZ*DH1(M2)/(HHH*HHH)) RHH =-HHH*(1.0D0-CDSQRT(FHH)) C AAH1(M2) = AAH1(M2-1)+RHH*AAH1(M2-2) BBH1(M2) = BBH1(M2-1)+RHH*BBH1(M2-2) RATIO = AAH1(M2)/BBH1(M2) C Cup = Sfact*REAL(RATIO) C C HHH = 5.0D-1*(1.0D0+ZZ*(DH2(M2-1)-DH2(M2))) FHH = 1.0D0+(ZZ*DH2(M2)/(HHH*HHH)) RHH =-HHH*(1.0D0-CDSQRT(FHH)) C AAH2(M2) = AAH2(M2-1)+RHH*AAH2(M2-2) BBH2(M2) = BBH2(M2-1)+RHH*BBH2(M2-2) RATIO = AAH2(M2)/BBH2(M2) C Cdn = Sfact*REAL(RATIO) C C HHH = 5.0D-1*(1.0D0+ZZ*(DH3(M2-1)-DH3(M2))) FHH = 1.0D0+(ZZ*DH3(M2)/(HHH*HHH)) RHH =-HHH*(1.0D0-CDSQRT(FHH)) AAH3(M2) = AAH3(M2-1)+RHH*AAH3(M2-2) BBH3(M2) = BBH3(M2-1)+RHH*BBH3(M2-2) RATIO = AAH3(M2)/BBH3(M2) C clp = Sfact*REAL(RATIO)*h/X rmww = Sfact*REAL(RATIO)*h*A*phi C HHH = 5.0D-1*(1.0D0+ZZ*(DH4(M2-1)-DH4(M2))) FHH = 1.0D0+(ZZ*DH4(M2)/(HHH*HHH)) RHH =-HHH*(1.0D0-CDSQRT(FHH)) AAH4(M2) = AAH4(M2-1)+RHH*AAH4(M2-2) BBH4(M2) = BBH4(M2-1)+RHH*BBH4(M2-2) RATIO = AAH4(M2)/BBH4(M2) C rmss = Sfact*REAL(RATIO)*A*phi C rmu = Cup*Vup rmd = Cdn*Vdn rmt = rmu+rmww+rmss+rmd C IF(dabs(rmww).GT.1.0d-20) THEN Cratio = rmss/rmww ELSE Cratio = 0.0d0 END IF IF(dabs(Cratio).LT.1.0d-20) Cratio = 0.0d0 C IF(dabs(Cup).LT.1.0d-20) Cup = 0.0d0 IF(dabs(Cdn).LT.1.0d-20) Cdn = 0.0d0 IF(dabs(Clp).LT.1.0d-20) Clp = 0.0d0 IF(dabs(rmd).LT.1.0d-20) rmd = 0.0d0 C IF(dabs(rmww).LT.1.0d-20) rmww = 0.0d0 IF(dabs(rmss).LT.1.0d-20) rmss = 0.0d0 C WRITE(2,6005) TIME(nuu),Cup,Cdn,Clp,Cratio WRITE(3,6005) TIME(nuu),rmu,rmww,rmss,rmd,rmt C Ctt(nuu) = Cdn C 1000 CONTINUE C WRITE(4,6006) (Ctt(n), n=0,ntime) C C C ************************************************************************ C * F O R M A T S * C ************************************************************************ C C 5000 FORMAT(3(A20)) C 6000 FORMAT(T5,'ENTER THE NAME OF THE DATA FILE') 6001 FORMAT(T5,'!!!!! UNAVAILABLE TYPE OF DIFFUSION !!!!!',/, & T5,'THE RUN WAS ABORTED - CORRECT AND TRY AGAIN',) 6002 FORMAT(T5,'!!!!! UNAVAILABLE SET OF EQUATIONS !!!!!',/, & T5,'THE RUN WAS ABORTED - CORRECT AND TRY AGAIN',) 6003 FORMAT(T5,'!!!!! UNAVAILABLE TIME OPTION !!!!!',/, & T5,'THE RUN WAS ABORTED - CORRECT AND TRY AGAIN',) 6005 FORMAT(T5,6(1PE15.8,5x)) 6006 FORMAT(4(1PE15.8,5x)) 6008 FORMAT(T5,//) C C C =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> End of DIFF_TH C C STOP END C C C C********************************************************************* C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************