C***********************************************************************
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C***********************************************************************
C
C
C
      PROGRAM DIFF_HX
! 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
! 
      REAL*8 Ki,Kd,KL,kk
      INTEGER xtyp
! 
      PARAMETER (PI=3.14159265358979324D0)
      PARAMETER (NM=25,NMT=2*NM,NME=(NM+1)*(NM+1),NMQ=NM*(NM+1))
!     
      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 Psi(0:NMT),Theta(0:NMT)
! 
      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)
! 
      DIMENSION LEN(0:NMT),LENE(-1:NM),LENQ(-1:NM)
      DIMENSION XX(0:500)
! 
      CHARACTER*20 DATANM
      CHARACTER*25 OOUUT1
! 
! 
! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>  Main body of DIFF_HX
! 
! 
! ********************************************************************
! *                           DATA INPUT                             *
! ********************************************************************
! 
! 
      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)'
! 
      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,*) M,FACTOR,ERRMAX
      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
! 
! ********************************************************************
! *              BASIC PARAMETERS & ARRAY DIMENSIONING               *
! ********************************************************************
! 
      LENE(-1) = 0
      LENQ(-1) = 0
      DO 5 I=0,NMT
         LEN(I)  = 0
    5 CONTINUE
      DO 6 I=0,NM
         LENE(I) = 0
         LENQ(I) = 0
    6 CONTINUE
! 
      M2      = 2*M
      LEN(0)  = M2
      LEN(1)  = M2
      DO 8 I=2,M2
         LEN(I) = LEN(I-1)-1
    8 CONTINUE
!     
      NSUME = 0
      NSUMQ = 0
      DO 9 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
    9 CONTINUE
C 
      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)
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 
      DO 1000 NUU=0,nxx
! 
         TPRIOD = FACTOR*TIME
         GAMMA = -DLOG(ERRMAX)/(2.0D0*TPRIOD)
!              
         Sfact = (1.0D0/TPRIOD)*DEXP(GAMMA*TIME)
! 
         DO 10 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 
            Psi(I) = beta*(gama*exp(rm*(XX(nuu)-2.0d0*X))
     &                    +exp(-rm*XX(nuu)))
C 
            IF(itype.EQ.2) Theta(I) = Psi(I)*u/w0
            IF(itype.EQ.3) Theta(I) = Psi(I)*v/w0
C
            IF(I.EQ.0) THEN
               Psi(0)   = 5.0d-1*Psi(0)
               Theta(0) = 5.0d-1*Theta(0)
            END IF
!     
   10    CONTINUE
! 
! ********************************************************************
! *                EVALUATION OF EH, QH, DH, AAH, BBH                *
! ********************************************************************
! 
         DO 15 I=1,M2
            EH1(I) = (0.0D0,0.0D0)
            IF(itype.GT.1) EH2(I) = (0.0D0,0.0D0)
   15    CONTINUE
!     
         DO 20 I=1,M2
            I1    = I-1
            QH1(I) = Psi(I)/(Psi(I1)+1.0d-35)
            IF(itype.GT.1) QH2(I) = Theta(I)/(Theta(I1)+1.0d-35)
   20    CONTINUE
!     
         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
!              
               EH1(NN) = QH1(NN1+1)-QH1(NN1)+EH1(NN2+1)
               IF(itype.GT.1) EH2(NN) = QH2(NN1+1)
     &                                 -QH2(NN1)+EH2(NN2+1)
   25       CONTINUE
!              
            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
!              
               QH1(NN) = QH1(NN2+1)*EH1(NN1+1)/(EH1(NN1)+1.0d-35)
               IF(itype.GT.1) 
     &            QH2(NN) = QH2(NN2+1)
     &                     *EH2(NN1+1)/(EH2(NN1)+1.0d-35)
   30       CONTINUE
   40    CONTINUE
!              
         DH1(0) = Psi(0)
         IF(itype.GT.1) DH2(0) = Theta(0)
!              
         DO 50 I=1,M
            I1  = 2*I-1
            I2  = 2*I
            NN1 = LENQ(I-1)+1
            NN2 = LENE(I-1)+1
!                
            DH1(I1) = -QH1(NN1)
            DH1(I2) = -EH1(NN2)
!                
            IF(itype.GT.1) THEN
               DH2(I1) = -QH2(NN1)
               DH2(I2) = -EH2(NN2)
            END IF
   50    CONTINUE
!              
!              
!              
         AAH1(-1) = (0.0D0,0.0D0)
         BBH1(-1) = (1.0D0,0.0D0)
         AAH1(0)  = DH1(0)
         BBH1(0)  = (1.0D0,0.0D0)
!              
         IF(itype.GT.1) THEN 
            AAH2(-1) = (0.0D0,0.0D0)
            BBH2(-1) = (1.0D0,0.0D0)
            AAH2(0)  = DH2(0)
            BBH2(0)  = (1.0D0,0.0D0)
         END IF
!              
         ARGG  = PI*TIME/TPRIOD
         AR1   = DCOS(ARGG)
         AR2   = DSIN(ARGG)
         ZZ    = DCMPLX(AR1,AR2)
!              
         DO 60 I=1,M2-1
            IM1   = I-1
            IM2   = I-2
!              
            ZZP     = DH1(I)*ZZ
            AAH1(I) = AAH1(IM1)+ZZP*AAH1(IM2)
            BBH1(I) = BBH1(IM1)+ZZP*BBH1(IM2)
!              
            IF(itype.GT.1) THEN 
               ZZP     = DH2(I)*ZZ
               AAH2(I) = AAH2(IM1)+ZZP*AAH2(IM2)
               BBH2(I) = BBH2(IM1)+ZZP*BBH2(IM2)
            END IF
   60    CONTINUE
!              
         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))
!              
         AAH1(M2) = AAH1(M2-1)+RHH*AAH1(M2-2)
         BBH1(M2) = BBH1(M2-1)+RHH*BBH1(M2-2)
         RATIO    = AAH1(M2)/BBH1(M2)
!              
         Cx = Sfact*REAL(RATIO) 
!              
         IF(itype.EQ.1) THEN 
            Fx = Kd*Ki*Cx
         ELSE 
            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))
!              
            AAH2(M2) = AAH2(M2-1)+RHH*AAH2(M2-2)
            BBH2(M2) = BBH2(M2-1)+RHH*BBH2(M2-2)
            RATIO    = AAH2(M2)/BBH2(M2)
!              
            Fx = Sfact*REAL(RATIO) 
		 END IF
!              
         IF(DABS(Cx).LT.1.0d-20) Cx = 0.0d0
         IF(DABS(Fx).LT.1.0d-20) Fx = 0.0d0
!              
         WRITE(2,6005) TIME,xx(nuu),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_XH
! 
! 
      STOP
      END
C
C
C
C***********************************************************************
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C***********************************************************************