C
C      ********************************************************
C      *                                                      *
C      *                   **** SEMINF ****                   *
C      *                                                      *
C      * ONE-DIMENSIONAL GROUND-WATER SOLUTE TRANSPORT MODEL  *
C      *                                                      *
C      *   FOR A SEMI-INFINITE SYSTEM WITH A FIRST-TYPE OR    *
C      *                                                      *
C      *         THIRD-TYPE BOUNDARY CONDITION AT X=0         *
C      *                                                      *
C      *            VERSION CURRENT AS OF 11/30/88            *
C      *                                                      *
C      ********************************************************
C
C       THE FOLLOWING CARD MUST BE CHANGED IF PROBLEM DIMENSIONS ARE
C       GREATER THAN THOSE GIVEN HERE.
C         MAXX = MAXIMUM NUMBER OF X-VALUES
C         MAXT = MAXIMUM NUMBER OF TIME VALUES
      PARAMETER (MAXX=100,MAXT=20)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL XP(MAXX),CP(MAXX),TP,XSCLP
      CHARACTER*10 CUNITS,VUNITS,DUNITS,KUNITS,LUNITS,TUNITS
      DIMENSION CXT(MAXX,MAXT),X(MAXX),T(MAXT)
      COMMON /IOUNIT/ IN,IO
C
C     PROGRAM VARIABLES
C
C           NOTE:  ANY CONSISTANT SET OF UNITS MAY BE USED IN THE
C           MODEL.  NO FORMAT STATEMENTS NEED TO BE CHANGED AS
C           LABELS FOR ALL VARIABLES ARE SPECIFIED IN MODEL INPUT.
C
C     C0       SOLUTE CONCENTRATION AT THE INFLOW BOUNDARY [M/L**3]
C     DX       LONGITUDINAL DISPERSION COEFFICIENT [L**2/T]
C     VX       GROUND-WATER VELOCITY IN X-DIRECTION [L/T]
C     DK       FIRST-ORDER SOLUTE DECAY CONSTANT [1/T]
C     X        X-POSITION AT WHICH CONCENTRATION IS EVALUATED [L]
C     T        TIME AT WHICH CONCENTRATION IS EVALUATED [T]
C     CN       NORMALIZED CONCENTRATION C/C0 [DIMENSIONLESS]
C     CXT      SOLUTE CONCENTRATION C(X,T) [M/L**3]
C
C     NBC      SOURCE BOUNDARY CONDITION TYPE (1 OR 3)
C     NX       NUMBER OF X-POSITIONS AT WHICH SOLUTION IS EVALUATED
C     NT       NUMBER OF TIME VALUES AT WHICH SOLUTION IS EVALUATED
C     IPLT     PLOT CONTROL. IF IPLT>0, CONCENTRATION PROFILES ARE PLOTTED
C
C       CHARACTER VARIABLES USED TO SPECIFY UNITS FOR MODEL PARAMETERS
C     CUNITS   UNITS OF CONCENTRATION (M/L**3)
C     VUNITS   UNITS OF GROUND-WATER VELOCITY (L/T)
C     DUNITS   UNITS OF DISPERSION COEFFICIENT (L**2/T)
C     KUNITS   UNITS OF SOLUTE DECAY CONSTANT  (1/T)
C     LUNITS   UNITS OF LENGTH (L)
C     TUNITS   UNITS OF TIME (T)
C
C       DEFINE INPUT/OUTPUT FILES AND PRINT TITLE PAGE
      CALL OFILE
      CALL TITLE4
      WRITE(IO,201)
C
C       READ IN MODEL PARAMETERS
      READ(IN,101) NBC,NX,NT,IPLT
      IF(NBC.EQ.1) WRITE(IO,202)
      IF(NBC.EQ.3) WRITE(IO,203)
      WRITE(IO,205) NX,NT
      READ(IN,105) CUNITS,VUNITS,DUNITS,KUNITS,LUNITS,TUNITS
      READ(IN,110) C0,VX,DX,DK,XSCLP
      WRITE(IO,210) C0,CUNITS,VX,VUNITS,DX,DUNITS,DK,KUNITS,XSCLP
      READ(IN,110) (X(I),I=1,NX)
      WRITE(IO,215) LUNITS
      WRITE(IO,220) (X(I),I=1,NX)
      READ(IN,110) (T(I),I=1,NT)
      WRITE(IO,225) TUNITS
      WRITE(IO,220) (T(I),I=1,NT)
C
C       BEGIN TIME LOOP
      DO 40 IT=1,NT
C
C       BEGIN X-COORDINATE LOOP
      DO 50 IX=1,NX
C
C       CALL ROUTINE TO CALCULATE NORMALIZED CONCENTRATION
C       BASED ON TYPE OF BOUNDARY CONDITION SPECIFIED
      IF(NBC.EQ.1) CALL CNRML1(DK,T(IT),X(IX),DX,VX,CN)
      IF(NBC.EQ.3) CALL CNRML3(DK,T(IT),X(IX),DX,VX,CN)
      CXT(IX,IT)=CN*C0
50    CONTINUE
C
C       CONVERT X AND C TO SINGLE PRECISION AND DIVIDE BY C0 TO
C       PLOT NORMALIZED CONCENTRATION PROFILE FOR EACH TIME VALUE.
      IF(IPLT.LT.1) GO TO 40
      DO 60 I=1,NX
      XP(I)=SNGL(X(I))
60    CP(I)=SNGL(CXT(I,IT)/C0)
      TP=SNGL(T(IT))
      CALL PLOT1D(XP,CP,NX,TP,IT,NT,TUNITS,LUNITS,XSCLP)
40    CONTINUE
C
C       PRINT OUT TABLES OF CONCENTRATION VALUES
      NPAGE=1+(NT-1)/9
      DO 80 NP=1,NPAGE
      IF(NP.EQ.1) WRITE(IO,230) TUNITS
      IF(NP.NE.1) WRITE(IO,231) TUNITS
      NP1=(NP-1)*9
      NP2=9
      IF((NP1+NP2).GT.NT) NP2=NT-NP1
      WRITE(IO,235) (T(NP1+J),J=1,NP2)
      WRITE(IO,236) CUNITS,LUNITS
      DO 70 IX=1,NX
      WRITE(IO,240) X(IX),(CXT(IX,NP1+J),J=1,NP2)
      IF(MOD(IX,45).NE.0) GO TO 70
      WRITE(IO,231) TUNITS
      WRITE(IO,235) (T(NP1+J),J=1,NP2)
      WRITE(IO,236) CUNITS,LUNITS
70    IF(MOD(IX,5).EQ.0 .AND. MOD(IX,45).NE.0) WRITE(IO,241)
80    CONTINUE
C
      CLOSE (IN)
      CLOSE (IO)
      STOP
C
C        FORMAT STATEMENTS
101   FORMAT(20I4)
105   FORMAT(8A10)
110   FORMAT(8F10.0)
201   FORMAT(/////1H ,30X,'ANALYTICAL SOLUTION TO THE ONE-DIMENSIONAL'/
     1 1H ,28X,'ADVECTIVE-DISPERSIVE SOLUTE TRANSPORT EQUATION'/
     2 1H ,38X,'FOR A SEMI-INFINITE SYSTEM'//
     3 1H ,35X,'PROGRAM WRITTEN BY E.J. WEXLER'/
     4 1H ,39X,'U.S. GEOLOGICAL SURVEY'/
     5 1H ,38X,'WATER RESOURCES DIVISION'///
     6 1H0,40X,'INPUT DATA'/1H ,40X,10(1H-))
202   FORMAT(1H0,25X,'FIRST-TYPE BOUNDARY CONDITION AT X = 0.0')
203   FORMAT(1H0,25X,'THIRD-TYPE BOUNDARY CONDITION AT X = 0.0')
205   FORMAT(1H0,25X,'NUMBER OF X-COORDINATES (NX) = ',I4/1H ,25X,
     1 'NUMBER OF TIME VALUES (NT) = ',I4)
210   FORMAT(1H0,25X,'SOLUTE CONCENTRATION ON MODEL BOUNDARY (C0) =',
     1 1P1E13.6,1X,A10/1H ,25X,
     2 'GROUND-WATER VELOCITY IN X-DIRECTION (VX) =',1P1E13.6,1X,A10/
     3 1H ,25X,'DISPERSION IN THE X-DIRECTION (DX) =',1P1E13.6,1X,A10/
     4 1H ,25X,'FIRST-ORDER SOLUTE DECAY RATE (DK) =',1P1E13.6,1X,A10/
     5 1H ,25X,'PLOT SCALING FACTOR (XSCLP) =',1P1E13.6)
215   FORMAT(1H0,25X,'X-COORDINATES AT WHICH SOLUTE CONCENTRATIONS ',
     1 'WILL BE CALCULATED, IN ',A10/1H ,25X,78(1H-)/)
220   FORMAT(1H ,5X,8F12.4)
225   FORMAT(1H0,25X,'TIMES AT WHICH SOLUTE CONCENTRATIONS '
     1 'WILL BE CALCULATED, IN ',A10/1H ,25X,70(1H-)/)
230   FORMAT(1H1/1H0,15X,'SOLUTE CONCENTRATION AS A FUNCTION OF TIME'/
     2 1H0,25X,'TIME VALUES, IN ',A10)
231   FORMAT(1H1/1H0,15X,'SOLUTE CONCENTRATION AS A FUNCTION OF TIME =',
     1 5X,'(CONTINUED)'/
     2 1H0,25X,'TIME VALUES, IN ',A10)
235   FORMAT(1H ,20X,9F12.4)
236   FORMAT(1H ,19X,'*',108(1H-)/
     1 1H ,4X,'X-COORDINATE,',2X,'!',44X,'SOLUTE CONCENTRATION, IN '
     2 A10/1H ,4X,'IN ',A10,2X,1H!/1H ,19X,'!')
240   FORMAT(1H ,5X,F12.4,2X,'! ',9F12.5)
241   FORMAT(1H ,19X,'!')
      END
      SUBROUTINE CNRML1(DK,T,X,D,V,CN)
      IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
C
C       SOLUTION FOR THE ONE-DIMENSIONAL SOLUTE TRANSPORT EQUATION
C       FOR A SEMI-INFINITE SYSTEM WITH A FIRST-TYPE SOURCE
C       BOUNDARY CONDITION. VALUE RETURNED IS THE NORMALIZED SOLUTE
C       CONCENTRATION AT A GIVEN X-COORDINATE AND TIME.
C       FOR NO SOLUTE DECAY, A SIMPLIFIED SOLUTION IS USED.
C
      ALPHA=2.0D0*DSQRT(D*T)
      U=DSQRT(V*V+4.0D0*D*DK)
      X2D=X/(2.0D0*D)
C
C       SOLUTION WITH SOLUTE DECAY
      IF(DK.EQ.0.0) GO TO 10
C
C       TERM 1
      X1=X2D*(V-U)
      Y1=(X-U*T)/ALPHA
      CALL EXERFC(X1,Y1,Z1)
C
C       TERM 2
      X2=X2D*(V+U)
      Y2=(X+U*T)/ALPHA
      CALL EXERFC(X2,Y2,Z2)
      CN=(Z1+Z2)/(2.0D0)
      RETURN
C
C       SOLUTION WITH NO SOLUTE DECAY
C       TERM 1
10    Y1=(X-V*T)/ALPHA
      CALL EXERFC(0.0D0,Y1,Z1)
      X2=X*V/D
C
C       TERM 2
      Y2=(X+V*T)/ALPHA
      CALL EXERFC(X2,Y2,Z2)
      CN=(Z1+Z2)/(2.0D0)
      RETURN
      END
      SUBROUTINE CNRML3(DK,T,X,D,V,CN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C       SOLUTION FOR THE ONE-DIMENSIONAL SOLUTE TRANSPORT EQUATION
C       FOR A SEMI-INFINITE SYSTEM WITH A THIRD-TYPE SOURCE
C       BOUNDARY CONDITION. VALUE RETURNED IS THE NORMALIZED SOLUTE
C       CONCENTRATION AT A GIVEN X-COORDINATE AND TIME.
C       FOR NO SOLUTE DECAY, A SIMPLIFIED SOLUTION IS USED.
C
      ALPHA=2.0D0*DSQRT(D*T)
      U=DSQRT(V*V+4.0D0*D*DK)
      X2D=X/(2.0D0*D)
      VXD=V*X/D
C
C       SOLUTION WITH SOLUTE DECAY
      IF(DK.EQ.0.0) GO TO 10
C
C       TERM 1
      X1=VXD-DK*T
      Y1=(X+V*T)/ALPHA
      CALL EXERFC(X1,Y1,Z1)
      Z1=Z1*2.0D0
C
C       TERM 2
      X2=X2D*(V-U)
      Y2=(X-U*T)/ALPHA
      CALL EXERFC(X2,Y2,Z2)
      Z2=Z2*(U/V-1.0D0)
C
C       TERM 3
      X3=X2D*(V+U)
      Y3=(X+U*T)/ALPHA
      CALL EXERFC(X3,Y3,Z3)
      Z3=Z3*(U/V+1.0D0)
      CN=V*V*(Z1+Z2-Z3)/(4.0D0*D*DK)
      RETURN
C
C       SOLUTION FOR NO SOLUTE DECAY
10    PI=3.14159265358979D0
C
C       TERM 1
      Y1=(X-V*T)/ALPHA
      CALL EXERFC(0.0D0,Y1,Z1)
      Z1=0.50D0*Z1
C
C       TERM 2
      X2=X*V/D
      Y2=(X+V*T)/ALPHA
      CALL EXERFC(X2,Y2,Z2)
      Z2=Z2*0.50D0*(1.0D0+V*(X+V*T)/D)
C
C       TERM 3
      Z3=DEXP(-1.0D0*(X-V*T)*(X-V*T)/(4.0D0*D*T))
      Z3=Z3*V*DSQRT(T/(PI*D))
      CN=Z1-Z2+Z3
      RETURN
      END
