C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C PROGRAM DIFF_TS C IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL*8 Ki,Kd,KL,kk INTEGER timtyp C DIMENSION FCTRAL(0:30),VV(30) DIMENSION TIME(0:1000) DIMENSION Ctt(0:1001) C CHARACTER*20 DATANM CHARACTER*26 OOUUT1,OOUUT2,OOUUT3 C C C =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of DIFF_TS C C C C ********************************************************************** C * DATA INPUTS * 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)S' OOUUT2 = DATANM(1:IM) //'.M(t)S' OOUUT3 = DATANM(1:IM) //'.HistS' 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,*) NITR READ(1,*) NPERIOD,timtyp 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 * EVALUATION OF THE FACTORIALS * C ********************************************************************** C FCTRAL(0)=1.0D0 DO 10 I=1,NITR RI=1.0D0*I FCTRAL(I)=FCTRAL(I-1)*RI 10 CONTINUE C C ********************************************************************** C * EVALUATION OF THE V'S FOR THE STEHFEST ALGORITHM * C ********************************************************************** C NHAF =NITR/2 DO 30 J=1,NITR C MUPL=MIN0(J,NHAF) MLOL=(J+1)/2 SUMV=0.0D0 C DO 20 M=MLOL,MUPL RM = M*1.0D0 M1 = NHAF-M M2 = M-1 M3 = J-M M4 = 2*M-J M5 = 2*M U0 = FCTRAL(M5)/FCTRAL(M4) U1 =(RM**NHAF)/FCTRAL(M1) U2 = FCTRAL(M)*FCTRAL(M2)*FCTRAL(M3) VU = U0*U1/U2 SUMV = SUMV+VU 20 CONTINUE MEXP=NHAF+J VV(J)=((-1.0D0)**MEXP)*SUMV 30 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) DL2 = DLOG(2.0D0) 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 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 I=1,ntime TINV = 1.0D0/TIME(I) SUM1 = 0.0D0 SUM2 = 0.0D0 SUM3 = 0.0D0 SUM4 = 0.0D0 C DO 40 J=1,NITR RJ = 1.0D0*J SI = DLOG(2.0D0)*RJ*TINV SINV = 1.0D0/SI C C ********************************************************************** C * DIFFUSION WITH LINEAR EQUILIBRIUM SORPTION * C ********************************************************************** C IF(itype.EQ.1) THEN R = h+w DD = DT+Ds*torts*w C C ********************************************************************** C * DIFFUSION WITH KINETIC LINEAR SORPTION * C ********************************************************************** C ELSE IF(itype.EQ.2) THEN u = u0/(si+kk+decay) R = h+u DD = DT+Ds*torts*u C C ********************************************************************** C * DIFFUSION WITH IRREVERSIBLE LINEAR SORPTION * C ********************************************************************** C ELSE v = v0/(si+decay) R = h+v DD = DT+Ds*torts*v END IF C rm = dsqrt(R*(si+decay)/DD) C gg1 = phi*DD*A*rm gg2 = Vdn*(si+decay) IF(Vdn.NE.0.0d0) THEN gama = (gg1-gg2)/(gg1+gg2) ELSE gama = 1.0d0 END IF C IF(nequ.EQ.1) THEN pp1 = Vup*(1.0d0+gama*dexp(-2.0d0*rm*X)) pp2 = (phi*A*R/rm)*(1.0d0-gama*dexp(-2.0d0*rm*X)) pp3 = (Vdn*(1.0d0+gama)-(phi*A*R/rm)*(1.0d0-gama)) & *dexp(-rm*X) beta = Vup/((si+decay)*(pp1+pp2+pp3)) ELSE gg3 = Vup*(si+decay) eta1 = gg3-gg1 eta2 = gg3+gg1 beta = Vup/(eta2+gama*eta1*dexp(-2.0d0*rm*X)) END IF C Psi0 = beta*(1.0d0+gama*dexp(-2.0d0*rm*X)) PsiX = beta*(1.0d0+gama)*dexp(-rm*X) wmmm = (beta/rm)*(1.0d0-dexp(-rm*X)) & *(1.0d0+gama*dexp(-rm*X)) C SUM1 = SUM1+VV(J)*Psi0 SUM2 = SUM2+VV(J)*PsiX SUM3 = SUM3+VV(J)*wmmm C IF(itype.EQ.1) THEN SUM4 = SUM4+VV(J)*wmmm*w ELSE IF(itype.EQ.2) THEN SUM4 = SUM4+VV(J)*wmmm*u ELSE SUM4 = SUM4+VV(J)*wmmm*v END IF C 40 CONTINUE C Cup = DL2*TINV*SUM1 Cdn = DL2*TINV*SUM2 Clp = DL2*TINV*SUM3*h/X C rmu = Cup*Vup rmw = DL2*TINV*SUM3*h*A*phi rms = DL2*TINV*SUM4*A*phi rmd = Cdn*Vdn rmt = rmu+rmw+rms+rmd C IF(dabs(rmw).GT.1.0d-10) THEN Cratio = rms/rmw ELSE Cratio = 0.0d0 END IF IF(dabs(Cratio).LT.1.0d-10) Cratio = 0.0d0 C IF(dabs(Cup).LT.1.0d-10) Cup = 0.0d0 IF(dabs(Cdn).LT.1.0d-10) Cdn = 0.0d0 IF(dabs(Clp).LT.1.0d-10) Clp = 0.0d0 IF(dabs(rmw).LT.1.0d-10) rmw = 0.0d0 IF(dabs(rms).LT.1.0d-10) rms = 0.0d0 IF(dabs(rmd).LT.1.0d-10) rmd = 0.0d0 C WRITE(2,6005) TIME(I),Cup,Cdn,Clp,Cratio WRITE(3,6005) TIME(I),rmu,rmw,rms,rmd,rmt C Ctt(i) = 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_TS C C STOP END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C***********************************************************************