$DEBUG
C      ANALYSIS PROGRAM FOR ONE DIMENSIONAL SOLUTE MOVEMENT AT SPEED u  
C                              IN SOIL
C
C                 MAIN PROGRAM
C
      PARAMETER(NNX1=101)
      DIMENSION TH0(NNX1),THI(NNX1),THH(NNX1),DDD(NNX1),UUU(NNX1),
    +TH0D(NNX1,10000),THHD(NNX1,10000)
    DOUBLE PRECISION DD0,TH0,THA,THB,AL,T0,AT0,DDX,DDT,TH0D,THHD,
    +R,TI,TI0,XJ,TH1,TH2,THI,AI,AKI,BETA1,BETA2,BETB1,BETB2,
     +SFCB1,SFCB2,UU0,SFC1,SFC2,BBU,BBS,DDD,UUU,Q3,QD
C
      CHARACTER AAAB*14,AAAC*14
C
      DATA MDLNO,THA,THB,AL/1,0.20D0,1.0D0,3.5D0/
      DATA DD0,BBS,UU0,BBU/0.03D0,0.4D0,1.D0,0.4D0/
    DATA T0,DDT,AT0/4.D0,1.D-3,0.6D0/
C
     NNX=NNX1-1
C
    WRITE(*,*)'Input Out-Data Filename 1 Please!'
      READ (*,*)AAAB
      WRITE(*,15)'Inputed Out-Data Filename 1 Is ',AAAB
      WRITE(*,*)'Input Out-Data Filename 2 Please!'
      READ (*,*)AAAC
      WRITE(*,15)'Inputed Out-Data Filename 2 Is ',AAAC
    WRITE(*,*)
15    FORMAT (1X,14A)
C
      OPEN(6,FILE=AAAB)
    OPEN(7,FILE=AAAC)
      WRITE(6,*)'     **************************************************
     +****'
      WRITE(6,*)'     *
     +   *'
      WRITE(6,75)
C
      WRITE(6,*)'     *
     +   *'
      WRITE(6,*)'     **************************************************
     +****'

75    FORMAT(6X,'* ONE DIMENSION FLOW ANALYSIS PROGRAM OF LOESS SOIL  *'
     +)
18    FORMAT (6X,14A/)
C
C     Input Main Parameters of One Dimensional Water Transplation
C
      DDX=AL/NNX
    NNT=T0/DDT
C
      WRITE(6,35)NNT,NNX,DDT,DDX,THA,THB,AL,T0,AT0,MDLNO,DD0,BBS,UU0,BBU
C
35    FORMAT(/1X,'THE TOTAL NUMBER FOR TIME STEPS',I7/,
     +1X,'THE TOTAL NUMBER FOR SPACE STEPS',I7/,
     +1X,'THE TIME STEP LENGTH: DDT=',E12.4/,
     +1X,'THE SPACE STEP LENGTH: DDX=',E12.4/,
     +1X,'THE INITIAL WATER CONTENT IN SOIL AREA: THA=',E12.4/,
     +1X,'THE WATER CONTENT AT X=0.0: THB=',E12.4/,
     +1X,'THE LENGTH OF ONE DIMENSIONAL SOIL AREA: AL=',E12.4/,
    +1X,'THE TOTAL TIME OF COMPUTATION: T0=',E12.4/,
    +1X,'THE TIME PERIOD OF PENETRATION: AT0=',E12.4/,
     +1X,'THE MODEL NO:=1(Lin)=2(Exp)=3(Pow): MDLNO=',I7/, 
     +1X,'THE SOLUTE DISPERSION COEFFICENT D0: DD0=',E12.4/
    +1X,'THE SOLUTE DISPERSION MODEL COEFFICENT b: BBS=',E12.4/
    +1X,'THE SOLUTE MOVEMENT SPEED u: UU0=',E12.4/
     +1X,'THE SOLUTE MOVEMENT MODEL COEFFICIENT bu: BBU=',E12.4/)
C    
C
C     Output Control Parameters of Structure
C
      WRITE(6,130)
    TI=0.0D0
C
    DO 33 J=1,NNX1
    XXJ=DFLOAT(J-1)/DFLOAT(NNX)
C
    IF(MDLNO.EQ.1)THEN
      DDD(J)=DD0*(1.D0+BBS*XXJ)
    UUU(J)=UU0*(1.D0+BBU*XXJ)
      ENDIF
C
    IF(MDLNO.EQ.2)THEN
      DDD(J)=DD0*(1.D0+XXJ**BBS)
    UUU(J)=UU0*(1.D0+XXJ**BBU)
      ENDIF
C
    IF(MDLNO.EQ.3)THEN
      DDD(J)=DD0*(1.D0+BBS**XXJ)
    UUU(J)=UU0*(1.D0+BBU**XXJ)
      ENDIF
33    CONTINUE
C
    IF(MDLNO.EQ.1)THEN
        UU0=UU0*(1.D0+0.5D0*BBU)
    DD0=DD0*(1.D0+0.5D0*BBS)
      ENDIF
C
    IF(MDLNO.EQ.2)THEN
        UU0=UU0*(1.D0+1.0D0/(1.0D0+BBU))
    DD0=DD0*(1.D0+1.0D0/(1.0D0+BBS))
      ENDIF
C
    IF(MDLNO.EQ.3)THEN
      UU0=UU0*(1.D0+(BBU-1.0D0)/DLOG(BBU))
    DD0=DD0*(1.D0+(BBS-1.0D0)/DLOG(BBS))
      ENDIF
C
      DO 5 J=2,NNX
     TH0(J)=THA
    THH(J)=THA
5     CONTINUE
      THH(1)=THB
    THH(NNX1)=THA
C
     KI=NNT/10
    KKI=NNT/100
      R=DDT/(DDX*DDX)
C
    LL=0
      DO 215 I=1,NNT
    TI=DFLOAT(I)*DDT
      TH0(1)=THB
    IF(TI.GT.AT0)TH0(1)=THA
    IF(TI.GT.AT0)THH(1)=THA
C
    TH0(NNX1)=THA
C
    DO 6 J=1,NNX1
    XJ=(J-1)*DDX
      IF(I.EQ.1)WRITE(7,95)TI,XJ,TH0(J),THH(J),DDD(J),UUU(J)
6     CONTINUE 
C
    AI=DFLOAT(I)
    AKI=DFLOAT(KI)
    IF(AI/AKI.EQ.DFLOAT(INT(AI/AKI)))WRITE(*,3)I,DFLOAT(I)/NNT*100.
3     FORMAT(1X,' K=',I10,' (',F5.1,'%)')
C
    IK=DINT(DFLOAT(I)/DFLOAT(KKI))
C
    DO 20 J=2,NNX
    XJ=(DFLOAT(J)-1.D0)*DDX
C    
     TH1=TH0(J-1)
     TH2=TH0(J+1)
    DD2=0.5D0*(DDD(J+1)+DDD(J))
    DD1=0.5D0*(DDD(J-1)+DDD(J))
C
    THI(J)=R*DD1*TH1+(1.D0-R*(DD1+DD2))*TH0(J)+R*DD2*TH2
    +-UUU(J)*0.5D0*DDT/DDX*(TH2-TH1)
C     
C    COMPUTING ANALYTICAL SOLUTION
C
      IF(DFLOAT(I)/DFLOAT(KKI).EQ.DFLOAT(IK))THEN
C
    XJU=XJ-UU0*TI
    BETA1=ABS(XJU)/(2.D0*DSQRT(DD0*TI))
     CALL EFC(BETA1,SFC1)
     IF(XJU.LT.0.0D0)SFC1=2.0D0-SFC1
    BETA2=(XJ+UU0*TI)/(2.D0*DSQRT(DD0*TI))
      CALL EFC(BETA2,SFC2)
C
    Q3=UU0*XJ/DD0/3.D0
C
    IF(Q3.LT.700.D0)THEN 
    QD=DEXP(Q3)
    ELSE
    WRITE(*,*)'U*L/D IS TOO LARGE'
    STOP
    ENDIF
C
    SSA=0.5D0*(THB-THA)*(SFC1+SFC2*QD*QD*QD)+THA
C    
    THH(J)=SSA
C
    IF(TI.GT.AT0)THEN
    TI0=TI-AT0
    XJU=XJ-UU0*TI0
    BETB1=ABS(XJU)/(2.D0*DSQRT(DD0*TI0))
C
      CALL EFC(BETB1,SFCB1)
    IF(XJU.LT.0.0D0)SFCB1=2.0D0-SFCB1
    BETB2=(XJ+UU0*TI0)/(2.D0*DSQRT(DD0*TI0))
      CALL EFC(BETB2,SFCB2)
C
    SSB=0.5D0*(THB-THA)*(SFCB1+SFCB2*QD*QD*QD)
    THH(J)=SSA-SSB
    ENDIF
    ENDIF
C
20    CONTINUE
C
    DO 30 J=2,NNX
    TH0(J)=THI(J)
30    CONTINUE
C
      IF(DFLOAT(I)/DFLOAT(KKI).EQ.DFLOAT(IK))THEN
    LL=LL+1
C    
      DO 22 J=1,NNX1
    XJ=(J-1)*DDX
    TH0D(J,LL)=TH0(J)
    THHD(J,LL)=THH(J)
     WRITE(7,95)TI,XJ,TH0(J),THH(J),DDD(J),UUU(J)
      WRITE(6,95)TI,XJ,TH0(J),THH(J),DDD(J),UUU(J)
22    CONTINUE
    ENDIF
C
215   CONTINUE
C
    WRITE(7,*)
    WRITE(6,*)
C
      DO 222 J=1,NNX1
    XJ=(J-1)*DDX
    DO 223 LA=1,LL
    TI=DFLOAT(LA*KKI)*DDT
     WRITE(7,95)XJ,TI,TH0D(J,LA),THHD(J,LA),DDD(J),UUU(J)
      WRITE(6,95)XJ,TI,TH0D(J,LA),THHD(J,LA),DDD(J),UUU(J)
223   CONTINUE
222   CONTINUE
C
C     Output Data of Elements, Subdomains & Domains
C
95    FORMAT(1X,1E10.5,6E13.5)
C
130   FORMAT(/1X,'TIME      COORDINATE   WAT CNT NUM  WAT CNT ANA  DISPE
     +RSION D VELOCITY V'/)
C
    CLOSE(6)
    CLOSE(7)
    STOP
    END
C
    SUBROUTINE EFC(BETA,SFC)
    DOUBLE PRECISION BETA,S,TERM,SFC,X1,X2,A,P2,XX1,XX2,TNNB,Y1,Y2,
    +DIS,XO,DISA,DISAP
C    
    P2=DSQRT(3.14159265358979323846D0)
    IF(BETA.LT.1.D0)THEN
    S=BETA
    TERM=BETA
    DO 1 J=1,30
    TERM=TERM*(-BETA*BETA)/DFLOAT(J)
    S=S+TERM/(2.D0*J+1.D0)
1    CONTINUE
    S=S*2.D0/P2
    SFC=1.D0-S
C
    ELSE IF(BETA.LT.5.D0)THEN
    SFC=1.0D0
    NNB=BETA*100
    A=1.D0/DFLOAT(NNB)
    DO 2 I=1,NNB
    X1=(I-1.D0)*BETA*A
    X2=I*BETA*A
    SFC=SFC-A*(DEXP(-X1*X1)+DEXP(-X2*X2))*BETA/P2
2     CONTINUE
C
    ELSE IF(BETA.LT.27.D0)THEN
    SFC=0.0D0
     TNNB=1.D4
    NNB=10000
    A=1.D0/DFLOAT(NNB)
    XO=1.D0/(BETA+TNNB)
    DIS=TNNB*XO/BETA
    DISA=DIS*A
    DISAP=DISA/P2
    DO 3 I=1,NNB
    XX1=XO+DISA*DFLOAT(I-1)
    XX2=XX1+DISA
    Y1=1.D0/(XX1*XX1)
    Y2=1.D0/(XX2*XX2)
    SFC=SFC+(Y1*DEXP(-Y1)+Y2*DEXP(-Y2))*DISAP
3     CONTINUE
      ELSE
    SFC=0.0D0
C
    ENDIF
      RETURN
      END