C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C PROGRAM DIFF_XS ! IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER xtyp REAL*8 Ki,Kd,KL,kk ! DIMENSION FCTRAL(0:30),VV(30) DIMENSION XX(0:100) ! CHARACTER*20 DATANM CHARACTER*25 OOUUT1 ! ! ! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of DIFF_XS ! ! ! ! ********************************************************************** ! * DATA INPUTS * ! ********************************************************************** ! 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 ! OOUUT1 = DATANM(1:IM) //'.C(x)S' ! OPEN(UNIT=1,FILE=DATANM,STATUS='OLD') OPEN(UNIT=2,FILE=OOUUT1) ! 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 ! READ(1,*) phi,rho,Sr READ(1,*) D0,Ds,decay READ(1,*) tortp,torti,torts READ(1,*) Kd,Ki,KL,kk ! IF(Sr.EQ.0.0d0) Ki = 1.0d0 IF(Ki.EQ.1.0d0) Sr = 0.0d0 ! READ(1,*) Vup,Vdn,A,X READ(1,*) NITR READ(1,*) TIME READ(1,*) NSUBD,xtyp IF(xtyp.NE.1.AND.xtyp.NE.2) THEN WRITE(2,6003) STOP END IF ! nxx = 0 XX(0) = 0.0d0 IF(xtyp.EQ.1) THEN xincr = X/REAL(NSUBD) DO 3 n=1,NSUBD nxx = n XX(n) = n*xincr 3 CONTINUE ELSE DO 4 n=1,NSUBD READ(1,*) NPOINT,xxincr C DO 4 I=1,NPOINT nxx = nxx+1 xx(nxx) = xx(nxx-1)+xxincr IF(xx(nxx).GT.X) THEN WRITE(2,6004) nxx,xx(nxx) STOP END IF 4 CONTINUE END IF ! ! ********************************************************************** ! * EVALUATION OF THE FACTORIALS * ! ********************************************************************** ! FCTRAL(0)=1.0D0 DO 10 I=1,NITR RI=1.0D0*I FCTRAL(I)=FCTRAL(I-1)*RI 10 CONTINUE ! ! ********************************************************************** ! * EVALUATION OF THE V'S FOR THE STEHFEST ALGORITHM * ! ********************************************************************** ! NHAF =NITR/2 DO 30 J=1,NITR ! MUPL=MIN0(J,NHAF) MLOL=(J+1)/2 SUMV=0.0D0 ! 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 ! ! h = 1.0d0-Sr+Sr*Ki w0 = ((1.0d0-phi)/phi)*rho w = w0*Kd*Ki DT = D0*(tortp*(1.0d0-Sr)+torti*Sr*Ki) DL2 = DLOG(2.0D0) 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*Ki END IF C C C C C ********************************************************************** C ********************************************************************** C ********************************************************************** C * THE TIME LOOP * C ********************************************************************** C ********************************************************************** C ********************************************************************** C C C C TINV = 1.0D0/TIME DO 1000 I=0,nxx SUM1 = 0.0D0 SUM2 = 0.0D0 SUM3 = 0.0D0 ! DO 100 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 ! Psi = beta*(gama*dexp(rm*(XX(i)-2.0d0*X))+dexp(-rm*XX(i))) ! IF(itype.EQ.2) Theta = Psi*u/w0 IF(itype.EQ.3) Theta = Psi*v/w0 ! SUM1 = SUM1+VV(J)*Psi IF(itype.GT.1) SUM2 = SUM2+VV(J)*Theta 100 CONTINUE ! Cx = DL2*TINV*SUM1 ! IF(itype.EQ.1) THEN Fx = Kd*Ki*Cx ELSE Fx = DL2*TINV*SUM2 END IF ! WRITE(2,6005) Time,XX(i),Cx,Fx 1000 CONTINUE ! ! ! ********************************************************************** ! * F O R M A T S * ! ********************************************************************** ! ! 5000 FORMAT(3(A20)) ! 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 SPACE OPTION !!!!!',/, & T5,'THE RUN WAS ABORTED - CORRECT AND TRY AGAIN',) 6004 FORMAT(T5,'!!!!! XX(',i3,') = ',1pe11.4, & ' > SAMPLE WIDTH X !!!!!',/, & T5,'THE RUN WAS ABORTED - CORRECT AND TRY AGAIN',) 6005 FORMAT(T5,5(1PE15.8,5x)) ! ! ! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> End of DIFF_XS ! ! STOP END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C***********************************************************************