TYPCRV
C     ******************************************************************00000010
C                                                                       00000020
C     PURPOSE                                                           00000030
C        TO COMPUTE TYPE CURVE FUNCTION VALUES FOR PARTIAL PENETRATION  00000040
C        IN A NONLEAKY AQUIFER USING EQUATIONS 1 AND 9A OF HANTUSH,M.S.,00000050
C        1961,DRAWDOWN AROUND A PARTIALLY PENETRATING WELL: HYDRAULIC   00000060
C        DIV. JOUR., AM. SOC. CIVIL ENGINEERS PROC., P. 83-98.          00000070
C     INPUT DATA                                                        00000080
C        1 CARD  -  FORMAT (3F5.1,I5,2E10.4)                            00000090
C             B - AQUIFER THICKNESS                                     00000100
C             L - DEPTH, BELOW TOP OF AQUIFER, TO BOTTOM OF PUMPING     00000110
C             WELL SCREEN                                               00000120
C             D - DEPTH, BELOW TOP OF AQUIFER, TO TOP OF PUMPING WELL   00000130
C             SCREEN                                                    00000140
C             NUM - NUMBER OF OBSERVATION WELLS OR PIEZOMETERS TIMES    00000150
C             NUMBER OF VALUES OF KZ/KR.                                00000160
C             SMALL - SMALLEST VALUE OF 1/U FOR WHICH COMPUTATION IS    00000170
C             DESIRED                                                   00000180
C             LARGE - LARGEST  VALUE OF 1/U FOR WHICH COMPUTATION IS    00000190
C             DESIRED                                                   00000200
C        NUM CARDS (ONE FOR EACH OBS. WELL OR PIEZOMETER AND FOR EACH  100000210
C             VALUE OF R*SQRT(KZ/KR). - FORMAT (3F5.1)                  00000220
C             R - RADIAL DISTANCE FROM PUMPED WELL TIMES SQRT(KZ/KR).   00000230
C             LPRIME - DEPTH, BELOW TOP OF AQUIFER, TO BOTTOM OF OBS.   00000240
C             WELL SCREEN (ZERO FOR PIEZOMETER)                         00000250
C             DPRIME - DEPTH, BELOW TOP OF AQUIFER, TO TOP OF OBS. WELL 00000260
C             SCREEN (TOTAL DEPTH FOR PIEZOMETER)                       00000270
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00000280
C        DQL12,SERIES,BESK,FCT,L,F,EXPI                                 00000290
C                                                                       00000300
C     ******************************************************************00000310
      REAL*8 U                                                          00000320
      REAL*4 L,LB,LPB,LPRIME,LARGE                                      00000330
      DIMENSION ARRAY(13,12), IARG(12), ARG(13), A(12), C(12)           00000340
      DATA ARG/1.,1.2,1.5,2.,2.5,3.,3.5,4.,5.,6.,7.,8.,9./              00000350
      DATA A/12*'  N*'/,C/12*'10**'/                                    00000360
      IRD=5                                                             00000370
      IPT=6                                                             00000380
      READ (IRD,6) B,L,D,NUM,SMALL,LARGE                                00000390
      LB=L/B                                                            00000400
      DB=D/B                                                            00000410
      IBEGIN=ALOG10(SMALL)                                              00000420
      IEND=ALOG10(LARGE)+1.                                             00000430
      JLIMIT=IEND-IBEGIN                                                00000440
      IF (JLIMIT.GT.12) JLIMIT=12                                       00000450
      DO 5 K=1,NUM                                                      00000460
      READ (IRD,6) R,LPRIME,DPRIME                                      00000470
      RB=R/B                                                            00000480
      LPB=LPRIME/B                                                      00000490
      DPB=DPRIME/B                                                      00000500
      DO 1 I=1,13                                                       00000510
      ARGI=ARG(I)                                                       00000520
      DO 1 J=1,JLIMIT                                                   00000530
      IARG(J)=IBEGIN+J-1                                                00000540
      Y=ARGI*10.**(IBEGIN+J-1)                                          00000550
      U=1./Y                                                            00000560
      X=U                                                               00000570
      CALL EXPI(X,WU,DUMMY)                                             00000580
    1 ARRAY(I,J)=WU+F(U,RB,LB,DB,LPB,DPB)                               00000590
      IF (LPB-0.) 2,2,3                                                 00000600
    2 WRITE (IPT,7) DPB,RB,LB,DB                                        00000610
      GO TO 4                                                           00000620
    3 WRITE (IPT,8) LPB,DPB,RB,LB,DB                                    00000630
    4 WRITE (IPT,9) (A(I),C(I),IARG(I),I=1,JLIMIT)                      00000640
      DO 5 I=1,13                                                       00000650
      WRITE (IPT,10) ARG(I),(ARRAY(I,J),J=1,JLIMIT)                     00000660
    5 CONTINUE                                                          00000670
      STOP                                                              00000680
C                                                                       00000690
C                                                                       00000700
    6 FORMAT (3F5.1,I5,2E10.4)                                          00000710
    7 FORMAT ('1','W(U)+F(U,R/B,L/B,D/B,Z/B), Z/B=',F5.2,', SQRT(KZ/KR)*00000720
     1R/B=',F5.2,', L/B=',F5.2,', D/B=',F5.2,', U=1/N')                 00000730
    8 FORMAT ('1','W(U)+F(U,R/B,L/B,D/B,L''/B,D''/B), L''/B=',F5.2,', D'00000740
     1'/B=',F5.2,', SQRT(KZ/KR)*R/B=',F5.2,', L/B=',F5.2,', D/B=',F5.2,'00000750
     2, U=1/N')                                                         00000760
    9 FORMAT ('0',2X,'N',1X,12(2A4,I2))                                 00000770
   10 FORMAT ((' ',F4.1,12(F9.4,1X)))                                   00000780
      END                                                               00000790
      REAL FUNCTION F*4(U,RB,LB,DB,LPB,DPB)                             00000800
C     ******************************************************************00000810
C                                                                       00000820
C     FUNCTION F                                                        00000830
C                                                                       00000840
C     PURPOSE                                                           00000850
C        TO COMPUTE DEPARTURES FROM THEIS CURVE CAUSED BY PARTIAL       00000860
C        PENETRATION OF PUMPED WELL.                                    00000870
C     USAGE                                                             00000880
C        F(U,RB,LB,DB,LPB,DPB)                                          00000890
C     DESCRIPTION OF PARAMETERS                                         00000900
C        ALL REAL, U DOUBLE PRECISION                                   00000910
C        U - R**2*S/4*T*TIME (RADIAL DISTANCE SQUARED * STORAGE         00000920
C            COEFFICIENT / 4*TRANSMISSIVITY * TIME                      00000930
C        RB - R/B ( RADIAL DISTANCE / AQUIFER THICKNESS )               00000940
C        LB - L/B ( FRACTION OF AQUIFER PENETRATED BY PUMPED WELL)      00000950
C        DB - D/B ( FRACTION OF AQUIFER ABOVE PUMPED WELL SCREEN)       00000960
C        LPB - L'/B (FRACTION OF AQUIFER PENETRATED BY OBS. WELL, ZERO  00000970
C              FOR PIEZOMETER)                                          00000980
C        DPB - D'/B (FRACTION OF AQUIFER ABOVE OBS. WELL SCREEN, TOTAL  00000990
C              DEPTH FOR PIEZOMETER)                                    00001000
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00001010
C        DQL12,SERIES,BESK,FCT,L                                        00001020
C     METHOD                                                            00001030
C        SUMS THE SERIES THROUGH N*PI*R/B EQ 20                         00001040
C                                                                       00001050
C     ******************************************************************00001060
      REAL*8 U,V                                                        00001070
      REAL*4 L,N,LB,LPB                                                 00001080
      SUM=0.                                                            00001090
      N=0.                                                              00001100
      PIRB=3.141593*RB                                                  00001110
      PILB=3.141593*LB                                                  00001120
      PIDB=3.141593*DB                                                  00001130
      IF (LPB-0.) 1,1,4                                                 00001140
C     CHECKS FOR WELL OR PIEZOMETER                                     00001150
    1 PIZB=3.141593*DPB                                                 00001160
    2 N=N+1.                                                            00001170
      V=N*PIRB/2.                                                       00001180
      IF (V.GT.10.) GO TO 3                                             00001190
C     TRUNCATES SERIES WHEN V>10                                        00001200
      X=L(U,V)/N                                                        00001210
      SUM=SUM+(SIN(N*PILB)-SIN(N*PIDB))*COS(N*PIZB)*X                   00001220
      GO TO 2                                                           00001230
    3 F=.6366198*SUM/(LB-DB)                                            00001240
      GO TO 7                                                           00001250
    4 PILPB=3.141593*LPB                                                00001260
      PIDPB=3.141593*DPB                                                00001270
    5 N=N+1                                                             00001280
      V=N*PIRB/2.                                                       00001290
      IF (V.GT.10.) GO TO 6                                             00001300
C     TRUNCATES SERIES WHEN V>10                                        00001310
      X=L(U,V)/N                                                        00001320
      SUM=SUM+(SIN(N*PILB)-SIN(N*PIDB))*(SIN(N*PILPB)-SIN(N*PIDPB))*X/N 00001330
      GO TO 5                                                           00001340
    6 F=.2026424*SUM/((LB-DB)*(LPB-DPB))                                00001350
    7 RETURN                                                            00001360
      END                                                               00001370
      SUBROUTINE EXPI(X,RES,AUX)                                        00001380
C                                                                       00001390
C     ..................................................................00001400
C                                                                       00001410
C        SUBROUTINE EXPI                                                00001420
C                                                                       00001430
C        PURPOSE                                                        00001440
C           COMPUTES THE EXPONENTIAL INTEGRAL -EI(-X)                   00001450
C                                                                       00001460
C        USAGE                                                          00001470
C           CALL EXPI(X,RES)                                            00001480
C                                                                       00001490
C        DESCRIPTION OF PARAMETERS                                      00001500
C           X     - ARGUMENT OF EXPONENTIAL INTEGRAL                    00001510
C           RES   - RESULT VALUE                                        00001520
C           AUX   - RESULTANT AUXILIARY VALUE                           00001530
C                                                                       00001540
C        REMARKS                                                        00001550
C           X GT 170 (X LT -174) MAY CAUSE UNDERFLOW (OVERFLOW)         00001560
C           WITH THE EXPONENTIAL FUNCTION                               00001570
C           FOR X = 0 THE RESULT VALUE IS SET TO -1.E75                 00001580
C                                                                       00001590
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  00001600
C           NONE                                                        00001610
C                                                                       00001620
C        METHOD                                                         00001630
C           DEFINITION                                                  00001640
C           RES=INTEGRAL(EXP(-T)/T, SUMMED OVER T FROM X TO INFINITY).  00001650
C           EVALUATION                                                  00001660
C           THREE DIFFERENT RATIONAL APPROXIMATIONS ARE USED IN THE     00001670
C           RANGES 1 LE X, X LE -9 AND -9 LT X LE -3 RESPECTIVELY,      00001680
C           A POLYNOMIAL APPROXIMATION IS USED IN -3 LT X LT 1.         00001690
C                                                                       00001700
C     ..................................................................00001710
C                                                                       00001720
      IF(X-1.)2,1,1                                                     00001730
    1 Y=1./X                                                            00001740
      AUX=1.-Y*(((Y+3.377358E0)*Y+2.052156E0)*Y+2.709479E-1)/((((Y*     00001750
     11.072553E0+5.716943E0)*Y+6.945239E0)*Y+2.593888E0)*Y+2.709496E-1) 00001760
      RES=AUX*Y*EXP(-X)                                                 00001770
      RETURN                                                            00001780
    2 IF(X+3.)6,6,3                                                     00001790
    3 AUX=(((((((7.122452E-7*X-1.766345E-6)*X+2.928433E-5)*X-2.335379E-400001800
     1)*X+1.664156E-3)*X-1.041576E-2)*X+5.555682E-2)*X-2.500001E-1)*X   00001810
     2+9.999999E-1                                                      00001820
      RES=-1.E75                                                        00001830
      IF(X)4,5,4                                                        00001840
    4 RES=X*AUX-ALOG(ABS(X))-5.772157E-1                                00001850
    5 RETURN                                                            00001860
    6 IF(X+9.)8,8,7                                                     00001870
    7 AUX=1.-((((5.176245E-2*X+3.061037E0)*X+3.243665E1)*X+2.244234E2)*X00001880
     1+2.486697E2)/((((X+3.995161E0)*X+3.893944E1)*X+2.263818E1)*X      00001890
     2+1.807837E2)                                                      00001900
      GOTO 9                                                            00001910
    8 Y=9./X                                                            00001920
      AUX=1.-Y*(((Y+7.659824E-1)*Y-7.271015E-1)*Y-1.080693E0)/((((Y     00001930
     1*2.518750E0+1.122927E1)*Y+5.921405E0)*Y-8.666702E0)*Y-9.724216E0) 00001940
    9 RES=AUX*EXP(-X)/X                                                 00001950
      RETURN                                                            00001960
      END                                                               00001970
C     ******************************************************************00001980
C                                                                       00001990
C     PURPOSE                                                           00002000
C        TO COMPUTE A TABLE OF VALUES OF THE LEAKY AQUIFER WELL         00002010
C        FUNCTION - W(U,R/B) - HANTUSH,M.S., AND JACOB,C.E., 1955,      00002020
C        NON-STEADY RADIAL FLOW IN AN INFINITE LEAKY AQUIFER: AM.       00002030
C        GEOPHYS. UNION TRANS., V. 36, NO. 1, P. 95-100.                00002040
C     INPUT DATA                                                        00002050
C        1 CARD - FORMAT(2E10.5)                                        00002060
C             USMALL - SMALLEST VALUE OF 1/U FOR WHICH COMPUTATION IS   00002070
C             DESIRED.                                                  00002080
C             ULARGE - LARGEST  VALUE OF 1/U FOR WHICH COMPUTATION IS   00002090
C             DESIRED.                                                  00002100
C        2 CARDS - FORMAT(8E10.5)                                       00002110
C             BDAT - 12 VALUES OF R/B FOR TABLE.                        00002120
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00002130
C        L,SERIES,FCT,BESK,DQL12                                        00002140
C                                                                       00002150
C     ******************************************************************00002160
      REAL*4 L                                                          00002170
      REAL*8 U,V                                                        00002180
      DIMENSION ARRAY(73,12), Y(73), BDAT(12), YNUM(6)                  00002190
      DATA YNUM/1.,1.5,2.,3.,5.,7./                                     00002200
      IRD=5                                                             00002210
      IPT=6                                                             00002220
      READ (IRD,6) USMALL,ULARGE                                        00002230
      READ (IRD,6) BDAT                                                 00002240
      IBEGIN=ALOG10(USMALL)                                             00002250
      IEND=ALOG10(ULARGE)+.99999                                        00002260
      ILIMIT=(IEND-IBEGIN)*6+1                                          00002270
      IF (ILIMIT.GT.73) ILIMIT=73                                       00002280
      DO 1 I=1,12                                                       00002290
      IF (BDAT(I).EQ.0.) GO TO 2                                        00002300
    1 CONTINUE                                                          00002310
      NB=12                                                             00002320
      GO TO 3                                                           00002330
    2 NB=I-1                                                            00002340
    3 II=0                                                              00002350
      DO 4 I=1,ILIMIT                                                   00002360
      II=II+1                                                           00002370
      IF (II.GT.6) II=1                                                 00002380
      IEXP=IBEGIN+(I-1)/6                                               00002390
      Y(I)=YNUM(II)*10.**IEXP                                           00002400
      U=1./Y(I)                                                         00002410
      DO 4 J=1,NB                                                       00002420
      V=BDAT(J)/2.                                                      00002430
    4 ARRAY(I,J)=L(U,V)                                                 00002440
      WRITE (IPT,7) (BDAT(I),I=1,NB)                                    00002450
      DO 5 I=1,ILIMIT                                                   00002460
    5 WRITE (IPT,8) Y(I),(ARRAY(I,J),J=1,NB)                            00002470
      STOP                                                              00002480
C                                                                       00002490
C                                                                       00002500
    6 FORMAT (8E10.5)                                                   00002510
    7 FORMAT ('1','W(U,R/B)'/'0',10X,'| R/B'/' ',6X,'1/U |',12E10.2)    00002520
    8 FORMAT (' ',E10.3,12F10.4)                                        00002530
      END                                                               00002540
C     ******************************************************************00002550
C                                                                       00002560
C     PURPOSE                                                           00002570
C        TO COMPUTE TYPE CURVE FUNCTION VALUES FOR H(U,BETA) --         00002580
C        HANTUSH,M.S.,1960, MODIFICATION OF THE THEORY OF LEAKY         00002590
C        AQUIFERS: JOUR. GEOPHYS. RES., V. 65, NO. 11,  P. 3713-3725.   00002600
C        THE COMPUTATIONAL ALGORITHM WAS DEVISED AND PROGRAMMED BY      00002610
C        S.S.PAPADOPULOS.                                               00002620
C     INPUT DATA                                                        00002630
C        1 CARD - FORMAT(2E10.5)                                        00002640
C             USMALL - SMALLEST(BEGINNING) VALUE OF 1/U.                00002650
C             ULARGE - LARGEST(ENDING) VALUE OF 1/U.                    00002660
C        2 CARDS - FORMAT(8E10.5)                                       00002670
C             BDAT - 12 VALUES OF BETA (ZERO OR BLANK VALUES ARE        00002680
C             PERMISSIBLE IF LESS THAN 12 DESIRED, WILL TERMINATE       00002690
C             AT FIRST ZERO OR BLANK VALUE).                            00002700
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00002710
C        H,DQG32,HUB,W - MUST BE INCLUDED IN DECK.                      00002720
C        DSQRT,DEXP,DERFC,DLOG - MUST BE IN COMPUTER LIBRARY.           00002730
C                                                                       00002740
C     ******************************************************************00002750
      REAL*8 U,BETA,H                                                   00002760
      DIMENSION ARRAY(73,12), Y(73), BDAT(12), YNUM(6)                  00002770
      DATA YNUM/1.,1.5,2.,3.,5.,7./                                     00002780
      IRD=5                                                             00002790
      IPT=6                                                             00002800
      READ (IRD,6) USMALL,ULARGE                                        00002810
      READ (IRD,6) BDAT                                                 00002820
      IBEGIN=ALOG10(USMALL)                                             00002830
      IEND=ALOG10(ULARGE)+.99999                                        00002840
      ILIMIT=(IEND-IBEGIN)*6+1                                          00002850
      IF (ILIMIT.GT.73) ILIMIT=73                                       00002860
      DO 1 I=1,12                                                       00002870
      IF (BDAT(I).EQ.0.) GO TO 2                                        00002880
    1 CONTINUE                                                          00002890
      NB=12                                                             00002900
      GO TO 3                                                           00002910
    2 NB=I-1                                                            00002920
      II=0                                                              00002930
    3 DO 4 I=1,ILIMIT                                                   00002940
      IEXP=IBEGIN+(I-1)/6                                               00002950
      II=II+1                                                           00002960
      IF (II.GT.6) II=1                                                 00002970
      Y(I)=YNUM(II)*10.**IEXP                                           00002980
      U=1./Y(I)                                                         00002990
      DO 4 J=1,NB                                                       00003000
      BETA=BDAT(J)                                                      00003010
    4 ARRAY(I,J)=H(U,BETA)                                              00003020
      WRITE (IPT,7) (BDAT(I),I=1,NB)                                    00003030
      DO 5 I=1,ILIMIT                                                   00003040
    5 WRITE (IPT,8) Y(I),(ARRAY(I,J),J=1,NB)                            00003050
      STOP                                                              00003060
C                                                                       00003070
    6 FORMAT (8E10.5)                                                   00003080
    7 FORMAT ('1','H(U,BETA)'/'0',10X,'| BETA'/' ',6X,'1/U |',12E10.2)  00003090
    8 FORMAT (' ',E10.3,12F10.4)                                        00003100
      END                                                               00003110
      DOUBLE PRECISION FUNCTION H(U,B)                                  00003120
C     ******************************************************************00003130
C                                                                       00003140
C     FUNCTION H                                                        00003150
C     PURPOSE                                                           00003160
C        TO COMPUTE THE INTEGRAL OF                                     00003170
C        EXP(-Y)*ERFC(B*SQRT(U)/SQRT(Y*(Y-U)))/Y SUMMED OVER Y          00003180
C        FROM U TO INFINITY (FUNCTION H(U,BETA) OF HANTUSH).            00003190
C     DESCRIPTION OF PARAMETERS                                         00003200
C        BOTH DOUBLE PRECISION                                          00003210
C        U - R**2*S/(4*T*TIME), (RADIAL DISTANCE SQUARED * STORAGE      00003220
C        COEFFICIENT / (4 * TRANSMISSIVITY * TIME). U MUST BE > 1.D-60. 00003230
C        B - (R/4)*(SQRT(K'*S'/(B'*T*S)+K''*S''/(B''*T*S)).             00003240
C             K',S',B' - HYD. COND., STORAGE COEFF., THICKNESS OF       00003250
C                  UPPER CONFINING BED.                                 00003260
C             K'',S'',B'' - HYD. COND., STORAGE COEFF., THICKNESS OF    00003270
C                  LOWER CONFINING BED.                                 00003280
C     METHOD                                                            00003290
C        I. FOR U < 1.D-60, NO COMPUTATION IS MADE.                     00003300
C        II. FOR B=0, H(U,0)=W(U) (THEIS WELL FUNCTION).                00003310
C        III. H(U,B)=0 IF                                               00003320
C             1. U > 10.                                                00003330
C             2. B > 1 AND B**2*U > 300.                                00003340
C        IV. ERFC(ARG)=0 FOR ARG > 40 AND H(U,B) = H(UB,B)              00003350
C             FOR U < Y < UB WHERE UB IS THE U CORRESPONDING TO ARG = 4000003360
C             SINCE H(UB,B) < W(UB) THEN FOR UB > 10, H(U,B) = 0.       00003370
C             ERFC(ARG) = 1 FOR ARG < 2.E-10 AND H(UUB,B) = W(UUB)      00003380
C             WHERE UUB IS THE U CORRESPONDING TO ARG = 2.E-10.         00003390
C             IF UUB > 10, H(U,B) = INTEGRAL FROM UB TO 10.             00003400
C             IF UUB < 10, H(U,B) = INTEGRAL FROM UB TO UUB + W(UUB)    00003410
C                                                                       00003420
C     ******************************************************************00003430
      IMPLICIT REAL*8(A-H,O-Z)                                          00003440
      COMMON UUU,BBB                                                    00003450
      EXTERNAL HUB                                                      00003460
      UUU=U                                                             00003470
      BBB=B                                                             00003480
      IF (U.GT.1.D-60) GO TO 1                                          00003490
      WRITE (6,7)                                                       00003500
      STOP                                                              00003510
    1 IF (B.EQ.0.0) GO TO 5                                             00003520
      IF (U.GT.10.0) GO TO 6                                            00003530
      BU=B*B*U                                                          00003540
      IF (B.GT.1.0.AND.BU.GE.300.0) GO TO 6                             00003550
      H1=0.0                                                            00003560
      UP=10.0                                                           00003570
      UB=0.5*U*(1.0+DSQRT(1.0+0.025*B*B/U))                             00003580
      IF (UB.GT.UP) GO TO 6                                             00003590
      UUB=0.5*U*(1.0+DSQRT(1.0+1.D20*B*B/U))                            00003600
      IF (UUB.GT.UP) GO TO 2                                            00003610
      H1=W(UUB)                                                         00003620
      UP=UUB                                                            00003630
    2 H2=0.0                                                            00003640
      XL=UB                                                             00003650
    3 XU=10.*XL                                                         00003660
      IF (XU.GE.UP) XU=UP                                               00003670
      CALL DQG32(XL,XU,HUB,AREA)                                        00003680
      H2=H2+AREA                                                        00003690
      XL=XU                                                             00003700
      IF (XL.EQ.UP) GO TO 4                                             00003710
      GO TO 3                                                           00003720
    4 H=H1+H2                                                           00003730
      RETURN                                                            00003740
    5 H=W(U)                                                            00003750
      RETURN                                                            00003760
    6 H=0.0                                                             00003770
      RETURN                                                            00003780
C                                                                       00003790
    7 FORMAT ('0','U TOO SMALL FOR COMPUTATION')                        00003800
      END                                                               00003810
      SUBROUTINE DQG32(XL,XU,FCT,Y)                                     00003820
C                                                                       00003830
C     ..................................................................00003840
C                                                                       00003850
C        SUBROUTINE DQG32                                               00003860
C                                                                       00003870
C        PURPOSE                                                        00003880
C           TO COMPUTE INTEGRAL(FCT(X), SUMMED OVER X FROM XL TO XU)    00003890
C                                                                       00003900
C        USAGE                                                          00003910
C           CALL DQG32 (XL,XU,FCT,Y)                                    00003920
C           PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT                00003930
C                                                                       00003940
C        DESCRIPTION OF PARAMETERS                                      00003950
C           XL     - DOUBLE PRECISION LOWER BOUND OF THE INTERVAL.      00003960
C           XU     - DOUBLE PRECISION UPPER BOUND OF THE INTERVAL.      00003970
C           FCT    - THE NAME OF AN EXTERNAL DOUBLE PRECISION FUNCTION  00003980
C                    SUBPROGRAM USED.                                   00003990
C           Y      - THE RESULTING DOUBLE PRECISION INTEGRAL VALUE.     00004000
C                                                                       00004010
C        REMARKS                                                        00004020
C           NONE                                                        00004030
C                                                                       00004040
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  00004050
C           THE EXTERNAL DOUBLE PRECISION FUNCTION SUBPROGRAM FCT(X)    00004060
C           MUST BE FURNISHED BY THE USER.                              00004070
C                                                                       00004080
C        METHOD                                                         00004090
C           EVALUATION IS DONE BY MEANS OF 32-POINT GAUSS QUADRATURE    00004100
C           FORMULA, WHICH INTEGRATES POLYNOMIALS UP TO DEGREE 63       00004110
C           EXACTLY. FOR REFERENCE, SEE                                 00004120
C           V.I.KRYLOV, APPROXIMATE CALCULATION OF INTEGRALS,           00004130
C           MACMILLAN, NEW YORK/LONDON, 1962, PP.100-111 AND 337-340.   00004140
C                                                                       00004150
C     ..................................................................00004160
C                                                                       00004170
      DOUBLE PRECISION XL,XU,Y,A,B,C,FCT                                00004180
      A=.5D0*(XU+XL)                                                    00004190
      B=XU-XL                                                           00004200
      C=.4986319309247408D0*B                                           00004210
      Y=.3509305004735048D-2*(FCT(A+C)+FCT(A-C))                        00004220
      C=.4928057557726342D0*B                                           00004230
      Y=Y+.813719736545284D-2*(FCT(A+C)+FCT(A-C))                       00004240
      C=.4823811277937532D0*B                                           00004250
      Y=Y+.1269603265463103D-1*(FCT(A+C)+FCT(A-C))                      00004260
      C=.4674530379688698D0*B                                           00004270
      Y=Y+.1713693145651072D-1*(FCT(A+C)+FCT(A-C))                      00004280
      C=.4481605778830261D0*B                                           00004290
      Y=Y+.2141794901111334D-1*(FCT(A+C)+FCT(A-C))                      00004300
      C=.4246838068662850D0*B                                           00004310
      Y=Y+.2549902963118809D-1*(FCT(A+C)+FCT(A-C))                      00004320
      C=.3972418979839712D0*B                                           00004330
      Y=Y+.2934204673926777D-1*(FCT(A+C)+FCT(A-C))                      00004340
      C=.3660910593701448D0*B                                           00004350
      Y=Y+.3291111138818092D-1*(FCT(A+C)+FCT(A-C))                      00004360
      C=.3315221334651076D0*B                                           00004370
      Y=Y+.3617289705442425D-1*(FCT(A+C)+FCT(A-C))                      00004380
      C=.2938578786203812D0*B                                           00004390
      Y=Y+.3909694789353515D-1*(FCT(A+C)+FCT(A-C))                      00004400
      C=.2534499544661147D0*B                                           00004410
      Y=Y+.4165596211347338D-1*(FCT(A+C)+FCT(A-C))                      00004420
      C=.2106756380653177D0*B                                           00004430
      Y=Y+.4382604650220191D-1*(FCT(A+C)+FCT(A-C))                      00004440
      C=.1659343011410638D0*B                                           00004450
      Y=Y+.4558693934788194D-1*(FCT(A+C)+FCT(A-C))                      00004460
      C=.1196436811260685D0*B                                           00004470
      Y=Y+.4692219954040228D-1*(FCT(A+C)+FCT(A-C))                      00004480
      C=.722359807913982D-1*B                                           00004490
      Y=Y+.4781936003963743D-1*(FCT(A+C)+FCT(A-C))                      00004500
      C=.2415383284386916D-1*B                                          00004510
      Y=B*(Y+.4827004425736390D-1*(FCT(A+C)+FCT(A-C)))                  00004520
      RETURN                                                            00004530
      END                                                               00004540
      DOUBLE PRECISION FUNCTION HUB(X)                                  00004550
C     ******************************************************************00004560
C                                                                       00004570
C     FUNCTION HUB                                                      00004580
C     PURPOSE                                                           00004590
C        TO COMPUTE VALUES OF THE INTEGRAND OF H(U,B)                   00004600
C     DESCRIPTION OF PARAMETER                                          00004610
C        X - DOUBLE PRECISION, POINT AT WHICH INTEGRAND IS EVALUATED.   00004620
C     METHOD                                                            00004630
C        FORTRAN EVALUATION OF FUNCTION.                                00004640
C                                                                       00004650
C     ******************************************************************00004660
      IMPLICIT REAL*8(A-H,O-Z)                                          00004670
      COMMON UUU,BBB                                                    00004680
      ARG=DSQRT((BBB*BBB*UUU)/(X*X-X*UUU))                              00004690
      HUB=DEXP(-X)*DERFC(ARG)/X                                         00004700
      RETURN                                                            00004710
      END                                                               00004720
      DOUBLE PRECISION FUNCTION W(U)                                    00004730
C     ******************************************************************00004740
C                                                                       00004750
C     FUNCTION W                                                        00004760
C     PURPOSE                                                           00004770
C        TO EVALUATE THE WELL FUNCTION OF THEIS.                        00004780
C     DESCRIPTION OF PARAMETER                                          00004790
C        U - DOUBLE PRECISION, ARGUMENT FOR WELL FUNCTION.              00004800
C                                                                       00004810
C     ******************************************************************00004820
      IMPLICIT REAL*8 (A-H,O-Z)                                         00004830
      IF (U.LE.0.0) GO TO 2                                             00004840
      IF (U.GT.100.) GO TO 3                                            00004850
      IF (U.GE.1.0) GO TO 1                                             00004860
      W=-.57721566+U*(.99999193+U*(-.24991055+U*(.05519968+U*(-.0097600400004870
     1+.00107857*U))))-DLOG(U)                                          00004880
      GO TO 4                                                           00004890
    1 ENUM=DEXP(-U)*(.2677737343+U*(8.6347608925+U*(18.0590169730+U*(8.500004900
     1733287401+U))))                                                   00004910
      DEN=U*(3.9584969228+U*(21.0996530827+U*(25.6329561486+U*(9.573322300004920
     1454+U))))                                                         00004930
      W=ENUM/DEN                                                        00004940
      GO TO 4                                                           00004950
    2 WRITE (6,5) U                                                     00004960
      STOP                                                              00004970
    3 W=0.0                                                             00004980
    4 RETURN                                                            00004990
C                                                                       00005000
    5 FORMAT ('0',5X,'W(U) NOT DEFINED FOR U=',1PD15.8)                 00005010
      END                                                               00005020
C     ******************************************************************00005030
C                                                                       00005040
C     PURPOSE                                                           00005050
C        TO COMPUTE TYPE CURVE FUNCTION VALUES FOR PARTIAL PENETRATION  00005060
C        IN A LEAKY AQUIFER USING EQ. 73 OF HANTUSH,M.S., 1964,         00005070
C        HYDRAULICS OF WELLS IN CHOW, VEN TE, ADVANCES IN HYDROSCIENCE, 00005080
C        VOL. 1: ACADEMIC PRESS, NEW YORK, P. 281-442.                  00005090
C     INPUT DATA                                                        00005100
C        1 CARD  -  FORMAT (3F5.1,I5,2E10.4)                            00005110
C             B - AQUIFER THICKNESS                                     00005120
C             E - DEPTH, BELOW TOP OF AQUIFER, TO BOTTOM OF PUMPING     00005130
C             WELL SCREEN                                               00005140
C             D - DEPTH, BELOW TOP OF AQUIFER, TO TOP OF PUMPING WELL   00005150
C             SCREEN                                                    00005160
C             NUM - NUMBER OF OBSERVATION WELLS OR PIEZOMETERS TIMES    00005170
C             NUMBER OF VALUES OF KZ/KR.                                00005180
C             SMALL - SMALLEST VALUE OF 1/U FOR WHICH COMPUTATION IS    00005190
C             DESIRED                                                   00005200
C             LARGE - LARGEST  VALUE OF 1/U FOR WHICH COMPUTATION IS    00005210
C             DESIRED                                                   00005220
C        2 CARDS - FORMAT(8E10.5)                                       00005230
C             BDAT - 12 VALUES OF R/BR, NON ZERO VALUES SHOULD BE       00005240
C             FIRST, WILL TERMINATE AT FIRST ZERO (OR BLANK) VALUE.     00005250
C        NUM CARDS (ONE FOR EACH OBS. WELL OR PIEZOMETER AND FOR EACH   00005260
C             VALUE OF R*SQRT(KZ/KR) - FORMAT (3F5.1)                   00005270
C             R - RADIAL DISTANCE FROM PUMPED WELL TIMES SQRT(KZ/KR).   00005280
C             LPRIME - DEPTH, BELOW TOP OF AQUIFER, TO BOTTOM OF OBS.   00005290
C             WELL SCREEN (ZERO FOR PIEZOMETER)                         00005300
C             DPRIME - DEPTH, BELOW TOP OF AQUIFER, TO TOP OF OBS. WELL 00005310
C             SCREEN (TOTAL DEPTH FOR PIEZOMETER)                       00005320
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00005330
C        DQL12,SERIES,BESK,FCT,L,FL                                     00005340
C                                                                       00005350
C     ******************************************************************00005360
      REAL*8 U,V                                                        00005370
      REAL*4 L,LB,LPB,LPRIME,LARGE                                      00005380
      DIMENSION ARRAY(55,12), ARG(6), BDAT(12), Y(55)                   00005390
      DATA ARG/1.,1.5,2.,3.,5.,7./                                      00005400
      DATA ARRAY/660*0./,Y/55*0./                                       00005410
      IRD=5                                                             00005420
      IPT=6                                                             00005430
      READ (IRD,9) B,E,D,NUM,SMALL,LARGE                                00005440
      READ (IRD,14) BDAT                                                00005450
      DO 1 I=1,12                                                       00005460
      IF (BDAT(I).EQ.0.) GO TO 2                                        00005470
    1 CONTINUE                                                          00005480
      NB=12                                                             00005490
      GO TO 3                                                           00005500
    2 NB=I-1                                                            00005510
    3 LB=E/B                                                            00005520
      DB=D/B                                                            00005530
      IBEGIN=ALOG10(SMALL)                                              00005540
      IEND=ALOG10(LARGE)+.1                                             00005550
      JLIMIT=IEND-IBEGIN                                                00005560
      IF (JLIMIT.GT.9) JLIMIT=9                                         00005570
      ILIMIT=6*JLIMIT+1                                                 00005580
      DO 8 K=1,NUM                                                      00005590
      READ (IRD,9) R,LPRIME,DPRIME                                      00005600
      RB=R/B                                                            00005610
      LPB=LPRIME/B                                                      00005620
      DPB=DPRIME/B                                                      00005630
      DO 4 I=1,ILIMIT                                                   00005640
      INDEX=(I-1)/6                                                     00005650
      IEXP=IBEGIN+INDEX                                                 00005660
      II=I-INDEX*6                                                      00005670
      Y(I)=ARG(II)*10.**IEXP                                            00005680
      U=1./Y(I)                                                         00005690
      DO 4 J=1,NB                                                       00005700
      BETA=BDAT(J)                                                      00005710
      V=BETA/2.                                                         00005720
    4 ARRAY(I,J)=L(U,V)+FL(U,RB,BETA,LB,DB,LPB,DPB)                     00005730
      IF (LPB-0.) 5,5,6                                                 00005740
    5 WRITE (IPT,10) DPB,RB,LB,DB                                       00005750
      GO TO 7                                                           00005760
    6 WRITE (IPT,11) LPB,DPB,RB,LB,DB                                   00005770
    7 WRITE (IPT,12) (BDAT(I),I=1,NB)                                   00005780
      DO 8 I=1,ILIMIT                                                   00005790
      WRITE (IPT,13) Y(I),(ARRAY(I,J),J=1,NB)                           00005800
    8 CONTINUE                                                          00005810
      STOP                                                              00005820
C                                                                       00005830
C                                                                       00005840
    9 FORMAT (3F5.1,I5,2E10.4)                                          00005850
   10 FORMAT ('1','W(U,R/BR)+F(U,R/B,R/BR,L/B,D/B,Z/B), Z/B=',F5.2,', SQ00005860
     1RT(KZ/KR)*R/B=',F5.2,', L/B=',F5.2,', D/B=',F5.2)                 00005870
   11 FORMAT ('1','W(U,R/BR)+F(U,R/B,R/BR,L/B,D/B,L''/B,D''/B), L''/B=',00005880
     1F5.2,', D''/B=',F5.2,', SQRT(KZ/KR)*R/B=',F5.2,', L/B=',F5.2,', D/00005890
     2B=',F5.2)                                                         00005900
   12 FORMAT ('0',9X,'| R/BR'/' ',5X,'1/U |',12E10.2)                   00005910
   13 FORMAT (' ',E10.3,12F10.4)                                        00005920
   14 FORMAT (8E10.5)                                                   00005930
      END                                                               00005940
      REAL FUNCTION FL*4(U,RB,BETA,LB,DB,LPB,DPB)                       00005950
C     ******************************************************************00005960
C                                                                       00005970
C     FUNCTION FL                                                       00005980
C                                                                       00005990
C     PURPOSE                                                           00006000
C        TO COMPUTE DEPARTURES FROM HANTUSH-JACOB LEAKY AQUIFER CURVE   00006010
C        CAUSED BY PARTIAL PENETRATION OF PUMPED WELL.                  00006020
C     USAGE                                                             00006030
C        FL(U,RB,BETA,LB,DB,LPB,DPB)                                    00006040
C     DESCRIPTION OF PARAMETERS                                         00006050
C        ALL REAL, U DOUBLE PRECISION                                   00006060
C        U - R**2*S/4*T*TIME (RADIAL DISTANCE SQUARED * STORAGE         00006070
C            COEFFICIENT / 4*TRANSMISSIVITY * TIME                      00006080
C        RB - R/B ( RADIAL DISTANCE / AQUIFER THICKNESS )               00006090
C        BETA - R*SQRT(K'/B'T) - (RADIAL DISTANCE * SQUARE ROOT         00006100
C               (HYD. COND. OF CONFINING BED/THICKNESS OF CONFINING     00006110
C               BED * TRANSMISSIVITY OF AQUIFER))                       00006120
C        LB - L/B ( FRACTION OF AQUIFER PENETRATED BY PUMPED WELL)      00006130
C        DB - D/B ( FRACTION OF AQUIFER ABOVE PUMPED WELL SCREEN)       00006140
C        LPB - L'/B (FRACTION OF AQUIFER PENETRATED BY OBS. WELL, ZERO  00006150
C              FOR PIEZOMETER)                                          00006160
C        DPB - D'/B (FRACTION OF AQUIFER ABOVE OBS. WELL SCREEN, TOTAL  00006170
C              DEPTH FOR PIEZOMETER)                                    00006180
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00006190
C        DQL12,SERIES,BESK,FCT,L                                        00006200
C     METHOD                                                            00006210
C        SUMS THE SERIES THROUGH N*PI*R/B EQ 20                         00006220
C                                                                       00006230
C     ******************************************************************00006240
      REAL*8 U,V,DSQRT                                                  00006250
      REAL*4 L,N,LB,LPB                                                 00006260
      SUM=0.                                                            00006270
      N=0.                                                              00006280
      BETSQ=BETA*BETA                                                   00006290
      PIRBSQ=9.869604*RB*RB                                             00006300
      PILB=3.141593*LB                                                  00006310
      PIDB=3.141593*DB                                                  00006320
      IF (LPB-0.) 1,1,4                                                 00006330
C     CHECKS FOR WELL OR PIEZOMETER                                     00006340
    1 PIZB=3.141593*DPB                                                 00006350
    2 N=N+1.                                                            00006360
      V=SQRT(BETSQ+N*N*PIRBSQ)/2.                                       00006370
      IF (V.GT.10.) GO TO 3                                             00006380
C     TRUNCATES SERIES WHEN V>10                                        00006390
      X=L(U,V)/N                                                        00006400
      SUM=SUM+(SIN(N*PILB)-SIN(N*PIDB))*COS(N*PIZB)*X                   00006410
      GO TO 2                                                           00006420
    3 FL=.6366198*SUM/(LB-DB)                                           00006430
      GO TO 7                                                           00006440
    4 PILPB=3.141593*LPB                                                00006450
      PIDPB=3.141593*DPB                                                00006460
    5 N=N+1                                                             00006470
      V=SQRT(BETSQ+N*N*PIRBSQ)/2.                                       00006480
      IF (V.GT.10.) GO TO 6                                             00006490
C     TRUNCATES SERIES WHEN V>10                                        00006500
      X=L(U,V)/N                                                        00006510
      SUM=SUM+(SIN(N*PILB)-SIN(N*PIDB))*(SIN(N*PILPB)-SIN(N*PIDPB))*X/N 00006520
      GO TO 5                                                           00006530
    6 FL=.2026424*SUM/((LB-DB)*(LPB-DPB))                               00006540
    7 RETURN                                                            00006550
      END                                                               00006560
C     ******************************************************************00006570
C                                                                       00006580
C     PURPOSE                                                           00006590
C        TO COMPUTE A TABLE OF FUNCTION VALUES FOR DRAWDOWN IN A        00006600
C        LEAKY ARTESIAN AQUIFER IN RESPONSE TO A STEP CHANGE IN         00006610
C        WATER LEVEL IN THE CONTROL WELL. FUNCTION VALUES ARE           00006620
C        EXPRESSED AS A FRACTION OF DRAWDOWN IN CONTROL WELL (S/SW).    00006630
C        REFERENCE - HANTUSH,M.S., 1959, NONSTEADY FLOW TO FLOWING      00006640
C        WELLS IN LEAKY AQUIFERS: JOUR. GEOPHYS. RESEARCH, V. 64,       00006650
C        NO. 8, P. 1043-1052.                                           00006660
C     INPUT DATA                                                        00006670
C        1 CARD - FORMAT(2E10.5)                                        00006680
C             TSMALL - SMALLEST VALUE OF ALPHA FOR WHICH COMPUTATION    00006690
C             IS DESIRED.                                               00006700
C             TLARGE -  LARGEST VALUE OF ALPHA FOR WHICH COMPUTATION    00006710
C             IS DESIRED.                                               00006720
C        1 CARD - FORMAT(13F5.0)                                        00006730
C             BDAT - 13 VALUES OF RW/B. NON ZERO VALUES SHOULD BE GE 1  00006740
C             AND LT 10. FIRST ZERO (OR BLANK) WILL TERMINATE THE       00006750
C             LIST. AT LEAST ONE NON ZERO VALUE MUST BE CODED. INPUT    00006760
C             VALUES ARE MULTIPLIED BY POWER OF TEN DETERMINED BY       00006770
C             PROGRAM FROM ALPHA.                                       00006780
C        1 CARD - FORMAT(10F8.2)                                        00006790
C             RW - RADIUS OF CONTROL WELL.                              00006800
C             RDAT - 9 VALUES OF RADIAL DISTANCE OF OBSERVATION POINTS  00006810
C             FROM CONTROL WELL. SHOULD BE CODED WITH SMALLEST NUMBER   00006820
C             FIRST, THEN BY INCREASING DISTANCE. THE FIRST ZERO        00006830
C             (OR BLANK) VALUE WILL TERMINATE COMPUTATION.              00006840
C     METHOD                                                            00006850
C        EVALUATES EQ. 13 OF HANTUSH. EVALUATION OF BESSEL FUNCTIONS    00006860
C        BY SUBROUTINES BESK AND BESY AND FUNCTION J0. EVALUATES        00006870
C        INTEGRAL BY SUM, I=1 TO 8000, F((DELTA U)*(I-.5))*(DELTA U).   00006880
C        CHOOSES INITIAL DELTA U = .001/SQRT(SMALLEST ALPHA) AND USES   00006890
C        THIS VALUE FOR ALL RW/B GE 10*(DELTA U). FOR SMALLER RW/B,     00006900
C        DIVIDES DELTA U BY 10 AND MULTIPLIES SMALLEST ALPHA BY 100.    00006910
C     REMARKS                                                           00006920
C        SMALLEST RW/B GE .01/SQRT(SMALLEST ALPHA)                      00006930
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00006940
C        BESK,BESY,J0                                                   00006950
C                                                                       00006960
C     ******************************************************************00006970
      REAL*8 SUM1,SUM2                                                  00006980
      REAL*4 KOBP,KOB,J0,JOPU,JOU,Y(8000),J(8000),F(8000),FT(8000),     00006990
     1   FB(8000),RDAT(9),TDAT(6),BDAT(13),ARRAY(25,9,13),B(13),T(25)   00007000
      DATA FT/8000*0./,FB/8000*0./                                      00007010
      DATA RDAT/9*1./                                                   00007020
      DATA ARRAY/2925*0./,TDAT/1.,1.5,2.,3.,5.,7./                      00007030
      IRD=5                                                             00007040
      IPT=6                                                             00007050
      READ (IRD,24) TSMALL,TLARGE                                       00007060
      READ (IRD,23) BDAT                                                00007070
      READ (IRD,22) RW,RDAT                                             00007080
      IBEGIN=ALOG10(TSMALL)                                             00007090
      IEND=ALOG10(TLARGE)+.99999                                        00007100
      IF ((IBEGIN/2*2).LT.IBEGIN) IBEGIN=IBEGIN-1                       00007110
      ISPAN=IEND-IBEGIN                                                 00007120
      MLIMIT=(ISPAN+1)/2                                                00007130
C     COMPUTES INITIAL DELTA U (DU) = .001/SQRT(SMALLEST ALPHA)         00007140
      DU=.001/SQRT(TDAT(1)*10.**IBEGIN)                                 00007150
C     EXPONENT (JBEGIN) OF SMALLEST RW/B IS COMPUTED FROM EXPONENT      00007160
C     (IBEGIN) OF SMALLEST ALPHA.                                       00007170
      JBEGIN=-IBEGIN/2-2                                                00007180
      DO 1 I=1,13                                                       00007190
      IF (BDAT(I).EQ.0.) GO TO 2                                        00007200
    1 CONTINUE                                                          00007210
      NB=13                                                             00007220
      GO TO 3                                                           00007230
    2 NB=I-1                                                            00007240
    3 CONTINUE                                                          00007250
      DO 4 I=1,9                                                        00007260
      IF (RDAT(I).EQ.0.) GO TO 5                                        00007270
    4 RDAT(I)=RDAT(I)/RW                                                00007280
      NR=9                                                              00007290
      GO TO 6                                                           00007300
    5 NR=I-1                                                            00007310
    6 DO 21 M=1,MLIMIT                                                  00007320
      NUM=8000                                                          00007330
      START=-DU/2.                                                      00007340
      U=START                                                           00007350
      DO 7 I=1,NUM                                                      00007360
      U=U+DU                                                            00007370
      CALL BESY(U,0,Y(I),IDUMY)                                         00007380
    7 J(I)=J0(U)                                                        00007390
      DO 19 IR=1,NR                                                     00007400
      RHO=RDAT(IR)                                                      00007410
      U=START                                                           00007420
      DO 8 I=1,NUM                                                      00007430
      U=U+DU                                                            00007440
      CALL BESY(RHO*U,0,YOPU,IDUMY)                                     00007450
      JOPU=J0(RHO*U)                                                    00007460
      JOU=J(I)                                                          00007470
      YOU=Y(I)                                                          00007480
    8 F(I)=(JOPU*YOU-YOPU*JOU)/(JOU*JOU+YOU*YOU)                        00007490
      DO 19 IT=1,25                                                     00007500
      INDEX=(IT-1)/6                                                    00007510
      IEXP=IBEGIN+INDEX                                                 00007520
      II=IT-INDEX*6                                                     00007530
      TAU=TDAT(II)*10.**IEXP                                            00007540
      T(IT)=TAU                                                         00007550
      U=START                                                           00007560
      NUMT=NUM                                                          00007570
      DO 9 I=1,NUMT                                                     00007580
      U=U+DU                                                            00007590
      FTEST=F(I)                                                        00007600
      IF (ABS(FTEST).LT.1.E-30) GO TO 10                                00007610
      XTEST=-TAU*U*U                                                    00007620
      IF (XTEST+69.) 10,10,9                                            00007630
    9 FT(I)=FTEST*EXP(XTEST)                                            00007640
      GO TO 11                                                          00007650
   10 NUMT=I-1                                                          00007660
      FT(I)=0.                                                          00007670
   11 DO 19 IB=1,13                                                     00007680
      JNDEX=(IB-1)/NB                                                   00007690
      JEXP=JBEGIN+JNDEX                                                 00007700
      JJ=IB-JNDEX*NB                                                    00007710
      BETA=BDAT(JJ)*10.**JEXP                                           00007720
      B(IB)=BETA                                                        00007730
      U=START                                                           00007740
      BSQ=BETA*BETA                                                     00007750
      NUMB=NUMT                                                         00007760
      DO 12 I=1,NUMB                                                    00007770
      U=U+DU                                                            00007780
      FTEST=FT(I)                                                       00007790
      IF (ABS(FTEST).LT.1.E-30) GO TO 13                                00007800
   12 FB(I)=FTEST/(U+BSQ/U)                                             00007810
      GO TO 14                                                          00007820
   13 NUMB=I-1                                                          00007830
      FB(I)=0.                                                          00007840
   14 SUM1=0.                                                           00007850
      SUM2=0.                                                           00007860
      DO 15 I=1,NUMB,2                                                  00007870
      SUM1=SUM1+FB(I)                                                   00007880
   15 SUM2=SUM2+FB(I+1)                                                 00007890
      XINT=(SUM1+SUM2)*DU                                               00007900
      CALL BESK(RHO*BETA,0,KOBP,IDUMY)                                  00007910
      CALL BESK(BETA,0,KOB,IDUMY)                                       00007920
      RATIO=0.                                                          00007930
      IF (KOBP.GT.0.) RATIO=KOBP/KOB                                    00007940
      XTEST=-TAU*BSQ                                                    00007950
      IF (XTEST+30.) 16,17,17                                           00007960
   16 XPT=0.                                                            00007970
      GO TO 18                                                          00007980
   17 XPT=EXP(XTEST)                                                    00007990
   18 Z=RATIO+.6366198*XPT*XINT                                         00008000
      IF ((Z.LT.0.).AND.(Z.GT.-5.E-5)) Z=0.E0                           00008010
   19 ARRAY(IT,IR,IB)=Z                                                 00008020
      DO 20 K=1,NR                                                      00008030
      WRITE (IPT,25) RDAT(K),B                                          00008040
      WRITE (IPT,26) (T(I),(ARRAY(I,K,L),L=1,13),I=1,25)                00008050
   20 CONTINUE                                                          00008060
C     EXPONENT OF SMALLEST RW/B DECREASED BY ONE EACH TIME THROUGH LOOP 00008070
      JBEGIN=JBEGIN-1                                                   00008080
C     EXPONENT OF SMALLEST ALPHA INCREASED BY TWO EACH TIME THROUGH LOOP00008090
      IBEGIN=IBEGIN+2                                                   00008100
C     DELTA U (DU) IS DIVIDED BY 10 EACH TIME THROUGH THE LOOP          00008110
   21 DU=.1*DU                                                          00008120
      STOP                                                              00008130
C                                                                       00008140
   22 FORMAT (10F8.2)                                                   00008150
   23 FORMAT (13F5.0)                                                   00008160
   24 FORMAT (2E10.5)                                                   00008170
   25 FORMAT ('1','Z(ALPHA,R/RW,RW/B), R/RW=',F6.0/'0',9X,'| RW/B'/(' ',00008180
     13X,'ALPHA |',13E9.2))                                             00008190
   26 FORMAT (' ',E10.3,13F9.3)                                         00008200
      END                                                               00008210
      REAL FUNCTION J0*4(X)                                             00008220
C     ******************************************************************00008230
C                                                                       00008240
C     FUNCTION J0                                                       00008250
C                                                                       00008260
C     PURPOSE                                                           00008270
C        TO COMPUTE THE ZERO ORDER J BESSEL FUNCTION FOR A GIVEN        00008280
C        ARGUMENT.                                                      00008290
C     USAGE                                                             00008300
C        J0(X)                                                          00008310
C     DESCRIPTION OF PARAMETER                                          00008320
C        X - REAL*4, ARGUMENT OF J0 BESSEL FUNCTION DESIRED.            00008330
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00008340
C        NONE.                                                          00008350
C     METHOD                                                            00008360
C        POLYNOMIAL APPROXIMATION FOR X<4 AND ASYMPTOTIC SERIES FOR     00008370
C        X GE 4. THE POLYNOMIAL APPROXIMATION IS THE FIRST 10 TERMS OF  00008380
C        THE POWER SERIES FOR J0(X) (MILLER, K.S., 1957,                00008390
C        ENGINEERING MATHEMATICS: RINEHART AND CO., INC., NEW YORK,     00008400
C        P. 120). THE ASYMPTOTIC EXPANSION OF J0(X) IS GIVEN ON P. 82   00008410
C        OF BOWMAN, FRANK, 1958, INTRODUCTION  TO BESSEL FUNCTIONS:     00008420
C        DOVER PUBLICATIONS INC., NEW YORK. THE TERMS P ('A*P0') AND    00008430
C        Q ('-B*Q0') OF THE ASYMPTOTIC EXPANSION ARE COMPUTED BY AN     00008440
C        ALGORITHM FROM IBM SUBROUTINE BESY.                            00008450
C                                                                       00008460
C     ******************************************************************00008470
      IF (X-4.) 1,3,3                                                   00008480
C     COMPUTE J0 BY FIRST 10 TERMS OF POWER SERIES                      00008490
    1 A=-X*X/4.                                                         00008500
      B=1.                                                              00008510
      DO 2 I=1,10                                                       00008520
      C=11.-I                                                           00008530
    2 B=1.+B*(A/(C*C))                                                  00008540
      J0=B                                                              00008550
      GO TO 4                                                           00008560
C     COMPUTE J0 BY ASYMPTOTIC SERIES                                   00008570
    3 T1=4./X                                                           00008580
      T2=T1*T1                                                          00008590
      P0=((((-.0000037043*T2+.0000173565)*T2-.0000487613)*T2+.00017343)*00008600
     1T2-.001753062)*T2+.3989423                                        00008610
      Q0=((((.0000032312*T2-.0000142078)*T2+.0000342468)*T2-.0000869791)00008620
     1*T2+.0004564324)*T2-.01246694                                     00008630
      A=2.0/SQRT(X)                                                     00008640
      B=A*T1                                                            00008650
      C=X-.7853982                                                      00008660
      J0=A*P0*COS(C)-B*Q0*SIN(C)                                        00008670
    4 RETURN                                                            00008680
      END                                                               00008690
      SUBROUTINE BESY(X,N,BY,IER)                                       00008700
C                                                                       00008710
C     ..................................................................00008720
C                                                                       00008730
C        SUBROUTINE BESY                                                00008740
C                                                                       00008750
C        PURPOSE                                                        00008760
C           COMPUTE THE Y BESSEL FUNCTION FOR A GIVEN ARGUMENT AND ORDER00008770
C                                                                       00008780
C        USAGE                                                          00008790
C           CALL BESY(X,N,BY,IER)                                       00008800
C                                                                       00008810
C        DESCRIPTION OF PARAMETERS                                      00008820
C           X  -THE ARGUMENT OF THE Y BESSEL FUNCTION DESIRED           00008830
C           N  -THE ORDER OF THE Y BESSEL FUNCTION DESIRED              00008840
C           BY -THE RESULTANT Y BESSEL FUNCTION                         00008850
C           IER-RESULTANT ERROR CODE WHERE                              00008860
C              IER=0  NO ERROR                                          00008870
C              IER=1  N IS NEGATIVE                                     00008880
C              IER=2  X IS NEGATIVE OR ZERO                             00008890
C              IER=3  BY HAS EXCEEDED MAGNITUDE OF 10**70               00008900
C                                                                       00008910
C        REMARKS                                                        00008920
C           VERY SMALL VALUES OF X MAY CAUSE THE RANGE OF THE LIBRARY   00008930
C           FUNCTION ALOG TO BE EXCEEDED                                00008940
C           X MUST BE GREATER THAN ZERO                                 00008950
C           N MUST BE GREATER THAN OR EQUAL TO ZERO                     00008960
C                                                                       00008970
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  00008980
C           NONE                                                        00008990
C                                                                       00009000
C        METHOD                                                         00009010
C           RECURRENCE RELATION AND POLYNOMIAL APPROXIMATION TECHNIQUE  00009020
C           AS DESCRIBED BY A.J.M.HITCHCOCK,'POLYNOMIAL APPROXIMATIONS  00009030
C           TO BESSEL FUNCTIONS OF ORDER ZERO AND ONE AND TO RELATED    00009040
C           FUNCTIONS', M.T.A.C., V.11,1957,PP.86-88, AND G.N. WATSON,  00009050
C           'A TREATISE ON THE THEORY OF BESSEL FUNCTIONS', CAMBRIDGE   00009060
C           UNIVERSITY PRESS, 1958, P. 62                               00009070
C                                                                       00009080
C     ..................................................................00009090
C                                                                       00009100
C                                                                       00009110
C     CHECK FOR ERRORS IN N AND X                                       00009120
C                                                                       00009130
      IF(N)180,10,10                                                    00009140
   10 IER=0                                                             00009150
      IF(X)190,190,20                                                   00009160
C                                                                       00009170
C     BRANCH IF X LESS THAN OR EQUAL 4                                  00009180
C                                                                       00009190
   20 IF(X-4.0)40,40,30                                                 00009200
C                                                                       00009210
C       COMPUTE Y0 AND Y1 FOR X GREATER THAN 4                          00009220
C                                                                       00009230
   30 T1=4.0/X                                                          00009240
      T2=T1*T1                                                          00009250
      P0=((((-.0000037043*T2+.0000173565)*T2-.0000487613)*T2            00009260
     1  +.00017343)*T2-.001753062)*T2+.3989423                          00009270
      Q0=((((.0000032312*T2-.0000142078)*T2+.0000342468)*T2             00009280
     1  -.0000869791)*T2+.0004564324)*T2-.01246694                      00009290
      P1=((((.0000042414*T2-.0000200920)*T2+.0000580759)*T2             00009300
     1  -.000223203)*T2+.002921826)*T2+.3989423                         00009310
      Q1=((((-.0000036594*T2+.00001622)*T2-.0000398708)*T2              00009320
     1  +.0001064741)*T2-.0006390400)*T2+.03740084                      00009330
      A=2.0/SQRT(X)                                                     00009340
      B=A*T1                                                            00009350
      C=X-.7853982                                                      00009360
      Y0=A*P0*SIN(C)+B*Q0*COS(C)                                        00009370
      Y1=-A*P1*COS(C)+B*Q1*SIN(C)                                       00009380
      GO TO 90                                                          00009390
C                                                                       00009400
C       COMPUTE Y0 AND Y1 FOR X LESS THAN OR EQUAL TO 4                 00009410
C                                                                       00009420
   40 XX=X/2.                                                           00009430
      X2=XX*XX                                                          00009440
      T=ALOG(XX)+.5772157                                               00009450
      SUM=0.                                                            00009460
      TERM=T                                                            00009470
      Y0=T                                                              00009480
      DO 70 L=1,15                                                      00009490
      IF(L-1)50,60,50                                                   00009500
   50 SUM=SUM+1./FLOAT(L-1)                                             00009510
   60 FL=L                                                              00009520
      TS=T-SUM                                                          00009530
      IF(ABS(TERM).LE.1.E-40) TERM=0.                                   00009540
      TERM=(TERM*(-X2)/FL**2)*(1.-1./(FL*TS))                           00009550
   70 Y0=Y0+TERM                                                        00009560
      TERM = XX*(T-.5)                                                  00009570
      SUM=0.                                                            00009580
      Y1=TERM                                                           00009590
      DO 80 L=2,16                                                      00009600
      SUM=SUM+1./FLOAT(L-1)                                             00009610
      FL=L                                                              00009620
      FL1=FL-1.                                                         00009630
      TS=T-SUM                                                          00009640
      IF(ABS(TERM).LE.1.E-40) TERM=0.                                   00009650
      TERM=(TERM*(-X2)/(FL1*FL))*((TS-.5/FL)/(TS+.5/FL1))               00009660
   80 Y1=Y1+TERM                                                        00009670
      PI2=.6366198                                                      00009680
      Y0=PI2*Y0                                                         00009690
      Y1=-PI2/X+PI2*Y1                                                  00009700
C                                                                       00009710
C     CHECK IF ONLY Y0 OR Y1 IS DESIRED                                 00009720
C                                                                       00009730
   90 IF(N-1)100,100,130                                                00009740
C                                                                       00009750
C     RETURN EITHER Y0 OR Y1 AS REQUIRED                                00009760
C                                                                       00009770
  100 IF(N)110,120,110                                                  00009780
  110 BY=Y1                                                             00009790
      GO TO 170                                                         00009800
  120 BY=Y0                                                             00009810
      GO TO 170                                                         00009820
C                                                                       00009830
C    PERFORM RECURRENCE OPERATIONS TO FIND YN(X)                        00009840
C                                                                       00009850
  130 YA=Y0                                                             00009860
      YB=Y1                                                             00009870
      K=1                                                               00009880
  140 T=FLOAT(2*K)/X                                                    00009890
      YC=T*YB-YA                                                        00009900
      IF(ABS(YC)-1.0E70)145,145,141                                     00009910
  141 IER=3                                                             00009920
      RETURN                                                            00009930
  145 K=K+1                                                             00009940
      IF(K-N)150,160,150                                                00009950
  150 YA=YB                                                             00009960
      YB=YC                                                             00009970
      GO TO 140                                                         00009980
  160 BY=YC                                                             00009990
  170 RETURN                                                            00010000
  180 IER=1                                                             00010010
      RETURN                                                            00010020
  190 IER=2                                                             00010030
      RETURN                                                            00010040
      END                                                               00010050
C***********************************************************************00010060
C                                                                       00010070
C     PURPOSE                                                           00010080
C        COMPUTES FUNCTION VALUES OF F(U,ALPHA,RHO) FOR RHO > 1 -       00010090
C        PAPADOPULOS,I.S. AND COOPER,H.H.,JR., 1967, DRAWDOWN IN        00010100
C        A WELL OF LARGE DIAMETER: WATER RESOURCES RESEARCH, V. 3,      00010110
C        NO. 1, P. 241-244.                                             00010120
C        PROGRAM BY S.S.PAPADOPULOS.                                    00010130
C     INPUT DATA - ONE OR MORE GROUPS, EACH GROUP CODED AS FOLLOWS      00010140
C        1 CARD - FORMAT(2E10.5)                                        00010150
C             ALPHA  -      RW**2*S/RC**2 - RADIUS OF WELL (SCREEN      00010160
C                  OR OPEN BORE IN AQUIFER) SQUARED * STORAGE           00010170
C                  COEFFICIENT / RADIUS OF CASING (OVER INTERVAL OF     00010180
C                  WATER LEVEL CHANGE) SQUARED.                         00010190
C             RHO - R/RW - DISTANCE FROM PUMPED WELL / RADIUS OF        00010200
C                  WELL (SCREEN OR OPEN BORE IN AQUIFER), MUST BE       00010210
C                  GREATER THAN ONE.                                    00010220
C        1 CARD - FORMAT(16E5.0)                                        00010230
C             U- 16 VALUES OF U - R**2*S/(4*T*TIME) - DISTANCE FROM     00010240
C                  PUMPED WELL SQUARED * STORAGE COEFFICIENT /          00010250
C                  4 * TRANSMISSIVITY * TIME. IF LESS THAN 16 DESIRED,  00010260
C                  BLANK OR ZERO VALUES MAY BE CODED FOR THE REST.      00010270
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00010280
C        PEAK,SIMP,APEKE,EXBSL1,JY0,JY1,ROOTS - MUST BE IN DECK.        00010290
C                                                                       00010300
C***********************************************************************00010310
      DIMENSION V(40,40),U(16)                                          00010320
      COMMON XPK,YPK                                                    00010330
      COMMON/PBLK/A,B,RHO                                               00010340
      EXTERNAL EXBSL1                                                   00010350
    1 READ (5,16,END=15) ALPHA,RHO                                      00010360
      IF (ALPHA) 15,15,2                                                00010370
    2 WRITE (6,17) ALPHA,RHO                                            00010380
      WRITE (6,18)                                                      00010390
    3 READ (5,19) U                                                     00010400
      DO 14 II=1,16                                                     00010410
      IF (U(II)) 1,1,4                                                  00010420
    4 A=ALPHA+ALPHA                                                     00010430
      B=0.25/U(II)                                                      00010440
      CALL APEKE(EXBSL1)                                                00010450
      CALL PEAK(EXBSL1)                                                 00010460
      IF (XPK-1.0E-8) 5,6,6                                             00010470
    5 WRITE (6,20) XPK,U                                                00010480
      GO TO 3                                                           00010490
    6 IF (XPK-3.0) 8,7,7                                                00010500
    7 WRITE (6,21) XPK,U                                                00010510
      GO TO 3                                                           00010520
    8 EPS=0.000001                                                      00010530
      HBAR=0.007*XPK                                                    00010540
      CALL SIMPS(0.0,XPK,EPS,HBAR,SUM,DEL,EXBSL1)                       00010550
      XM1=((3.14159265*7.0)/(8.0*(RHO-1.))+1.E-6)*RHO/2.                00010560
      DX1=XM1-(1.0E-6)*RHO                                              00010570
      DXN=(2.0*3.14159265*RHO)/(5.*(RHO-1.))                            00010580
      DL=3.14159265*RHO/(RHO-1.)                                        00010590
      CALL ROOTS(XM1,DX1,RT1,EXBSL1)                                    00010600
      HBAR=0.007*(RT1-XPK)                                              00010610
      CALL SIMPS(XPK,RT1,EPS,HBAR,TRM1,ERR1,EXBSL1)                     00010620
      SUM=SUM+TRM1                                                      00010630
      DEL=DEL+ERR1                                                      00010640
      X1=RT1                                                            00010650
      I=1                                                               00010660
    9 XM=X1+DL                                                          00010670
      CALL ROOTS(XM,DXN,X2,EXBSL1)                                      00010680
      HBAR=0.007*(X2-X1)                                                00010690
      CALL SIMPS(X1,X2,EPS,HBAR,TRM,ERR,EXBSL1)                         00010700
      V(1,I)=ABS(TRM)                                                   00010710
      DEL=DEL+ERR                                                       00010720
      I=I+1                                                             00010730
      IF (I-40) 10,10,11                                                00010740
   10 X1=X2                                                             00010750
      GO TO 9                                                           00010760
   11 EST=0.0                                                           00010770
      DO 12 K=2,40                                                      00010780
      M=41-K                                                            00010790
      DO 12 J=1,M                                                       00010800
   12 V(K,J)=V(K-1,J+1)-V(K-1,J)                                        00010810
      DO 13 N=1,40                                                      00010820
      L=N-1                                                             00010830
      DELV=(-0.5)**L*V(N,1)                                             00010840
   13 EST=EST+(0.5)*DELV                                                00010850
      SUM=SUM-EST                                                       00010860
      PUAR=4.0*A*RHO*SUM/3.14159265                                     00010870
      WRITE (6,22) U(II),SUM,DEL,PUAR                                   00010880
   14 CONTINUE                                                          00010890
      GO TO 1                                                           00010900
   15 STOP                                                              00010910
C                                                                       00010920
   16 FORMAT (2E10.5)                                                   00010930
   17 FORMAT ('1','F(U,ALPHA,RHO) FOR ALPHA=',1PE13.5,', RHO=',1E13.5)  00010940
   18 FORMAT (1H0,12X,1HU,16X,8HINTEGRAL,9X,14HINTEGRAL ERROR,6X,14HF(U,00010950
     1ALPHA,RHO)/1H )                                                   00010960
   19 FORMAT (16E5.0)                                                   00010970
   20 FORMAT (5H XPK=,E15.8,3X,16HTOO SMALL FOR U=,E10.3)               00010980
   21 FORMAT (5H XPK=,E15.8,3X,16HTOO LARGE FOR U=,E10.3)               00010990
   22 FORMAT (1H ,1P4E20.8)                                             00011000
      END                                                               00011010
      FUNCTION EXBSL1(X)                                                00011020
C***********************************************************************00011030
C                                                                       00011040
C     PURPOSE                                                           00011050
C        COMPUTES VALUES OF THE INTEGRAND FOR F(U,ALPHA,RHO)            00011060
C     DESCRIPTION OF PARAMETER                                          00011070
C        X- REAL - ARGUMENT OF INTEGRAND                                00011080
C                                                                       00011090
C***********************************************************************00011100
      COMMON/PBLK/A,B,R                                                 00011110
      IF (X) 1,1,2                                                      00011120
    1 EXBSL1=0.                                                         00011130
      GO TO 8                                                           00011140
    2 W=X/R                                                             00011150
      IF (W-1.0E7) 4,4,3                                                00011160
    3 FNU=A*COS(W*(R-1.0))-W*SIN(W*(R-1.0))                             00011170
      DE=(W*W*SQRT(R))*(W*W+A*A)                                        00011180
      EXBSL1=FNU/DE                                                     00011190
      GO TO 8                                                           00011200
    4 Y=B*X*X                                                           00011210
      IF (Y-0.01) 5,5,6                                                 00011220
    5 EXPO=Y*(1.0-Y*(0.5-Y*((1.0/6.0)-Y*(1.0/24.0))))                   00011230
      GO TO 7                                                           00011240
    6 EXPO=1.0-EXP(-Y)                                                  00011250
    7 CALL JY0(W,WJ0,WY0)                                               00011260
      CALL JY1(W,WJ1,WY1)                                               00011270
      AW=W*WY0-A*WY1                                                    00011280
      BW=W*WJ0-A*WJ1                                                    00011290
      CALL JY0(X,BJ0,BY0)                                               00011300
      FNUM=EXPO*(AW*BJ0-BW*BY0)                                         00011310
      DEN=X*X*(AW*AW+BW*BW)                                             00011320
      EXBSL1=FNUM/DEN                                                   00011330
    8 RETURN                                                            00011340
      END                                                               00011350
      SUBROUTINE ROOTS(XM,DX,ROOT,F)                                    00011360
C***********************************************************************00011370
C                                                                       00011380
C     PURPOSE                                                           00011390
C        SEARCHES FOR ROOT OF F IN THE INTERVAL XM-DX TO XM+DX.         00011400
C     DESCRIPTION OF PARAMETERS - ALL REAL                              00011410
C        XM - CENTER OF INTERVAL SEARCHED.                              00011420
C        DX - HALF WIDTH OF INTERVAL SEARCHED.                          00011430
C        ROOT - RETURNED ROOT LOCATION.                                 00011440
C        F - FUNCTION REFERENCE.                                        00011450
C                                                                       00011460
C***********************************************************************00011470
      XL=XM-DX                                                          00011480
      XR=XM+DX                                                          00011490
      YL=F(XL)                                                          00011500
      YR=F(XR)                                                          00011510
      EP=0.000001*ABS(YL)                                               00011520
      DO 9 I=1,200                                                      00011530
      YM=F(XM)                                                          00011540
      UP=ABS(YM)                                                        00011550
      IF (UP.LT.EP.AND.UP.LT.1.0D-7) GO TO 1                            00011560
      IF (YM) 2,1,2                                                     00011570
    1 ROOT=XM                                                           00011580
      GO TO 10                                                          00011590
    2 IF (YM*YL) 7,3,4                                                  00011600
    3 ROOT=XL                                                           00011610
      GO TO 10                                                          00011620
    4 IF (YM*YR) 8,5,6                                                  00011630
    5 ROOT=XR                                                           00011640
      GO TO 10                                                          00011650
    6 WRITE (6,11) XL,XR                                                00011660
      STOP                                                              00011670
    7 XR=XM                                                             00011680
      YR=YM                                                             00011690
      GO TO 9                                                           00011700
    8 XL=XM                                                             00011710
      YL=YM                                                             00011720
    9 XM=(XL+XR)/2.0                                                    00011730
      ROOT=XM                                                           00011740
   10 RETURN                                                            00011750
C                                                                       00011760
   11 FORMAT (1H ,10X,27HNO ROOT IN INTERVAL XM-DX =,1PE20.8,5X,11HAND X00011770
     1M+DX =,1PE20.8/)                                                  00011780
      END                                                               00011790
      SUBROUTINE APEKE(EXBSL)                                           00011800
C***********************************************************************00011810
C                                                                       00011820
C     PURPOSE                                                           00011830
C        GETS FIRST APPROXIMATION TO PEAK POSITION                      00011840
C                                                                       00011850
C***********************************************************************00011860
      COMMON XPK,YPK                                                    00011870
      XPK=0.0                                                           00011880
      YPK=0.0                                                           00011890
      DO 2 I=1,17                                                       00011900
      X=10.0**(I-9)                                                     00011910
      Y=EXBSL(X)                                                        00011920
      IF (Y-YPK) 3,3,1                                                  00011930
    1 XPK=X                                                             00011940
      YPK=Y                                                             00011950
    2 CONTINUE                                                          00011960
    3 RETURN                                                            00011970
      END                                                               00011980
      SUBROUTINE PEAK(EXBSL)                                            00011990
C***********************************************************************00012000
C                                                                       00012010
C     PURPOSE                                                           00012020
C        ATTEMPTS TO FIND POSITION OF MAXIMUM FOR INTEGRAND             00012030
C                                                                       00012040
C***********************************************************************00012050
      COMMON XPK,YPK                                                    00012060
      YPK=EXBSL(XPK)                                                    00012070
      DO 13 L=1,200                                                     00012080
      DX=0.01*XPK                                                       00012090
      XL=XPK-DX                                                         00012100
      YL=EXBSL(XL)                                                      00012110
      XR=XPK+DX                                                         00012120
      YR=EXBSL(XR)                                                      00012130
      DEN=YR+YL-YPK-YPK                                                 00012140
      IF (DEN) 1,9,1                                                    00012150
    1 X=XPK-0.5*(YR-YL)*DX/DEN                                          00012160
    2 IF (X) 3,4,4                                                      00012170
    3 X=0.0                                                             00012180
    4 Y=EXBSL(X)                                                        00012190
      IF (YR-Y) 6,6,5                                                   00012200
    5 Y=YR                                                              00012210
      X=XR                                                              00012220
    6 IF (YL-Y) 8,8,7                                                   00012230
    7 Y=YL                                                              00012240
      X=XL                                                              00012250
    8 IF (Y-YPK) 14,14,12                                               00012260
    9 IF (YR-YPK) 11,10,10                                              00012270
   10 X=XPK+DX+DX                                                       00012280
      GO TO 2                                                           00012290
   11 X=XPK-DX-DX                                                       00012300
      GO TO 2                                                           00012310
   12 YPK=Y                                                             00012320
      XPK=X                                                             00012330
   13 CONTINUE                                                          00012340
   14 RETURN                                                            00012350
      END                                                               00012360
      SUBROUTINE SIMPS(Q,R,EPS,HBAR,AREA,DEL,F)                         00012370
C***********************************************************************00012380
C                                                                       00012390
C     PURPOSE                                                           00012400
C        TO DETERMINE THE INTEGRAL OF A FUNCTION, F, FROM Q TO R,       00012410
C        USING SIMPSON'S RULE.                                          00012420
C     DESCRIPTION OF PARAMETERS                                         00012430
C        ALL REAL                                                       00012440
C        Q - LOWER LIMIT OF INTEGRAL                                    00012450
C        R - UPPER LIMIT OF INTEGRAL                                    00012460
C        EPS - DESIRED ACCURACY                                         00012470
C        HBAR - MINIMUM DIVISION OF THE INTERVAL                        00012480
C        AREA - COMPUTED VALUE OF INTEGRAL BETWEEN Q AND R              00012490
C        DEL - COMPUTED ESTIMATE OF ERROR                               00012500
C        F- THE INTEGRAND (FUNCTION REFERENCE)                          00012510
C     METHOD                                                            00012520
C        USES SIMPSON'S RULE TO COMPUTE A SUM APPROXIMATING THE INTEGRAL00012530
C        USES INITIAL H=(R-Q)/2. COMPUTES A SEQUENCE OF SUMS BY HALVING 00012540
C        H EACH TIME. COMPUTES ESTIMATE OF ERROR (DEL) AS (PREVIOUS     00012550
C        SUM - CURRENT SUM)/15. COMPUTATION STOPS WHEN 1) H<HBAR,       00012560
C        2) ABS(DEL)<ABS(EPS*CURRENT SUM). IF HBAR IS LE 0,             00012570
C        THEN HBAR=.007*(R-Q).                                          00012580
C                                                                       00012590
C***********************************************************************00012600
      H=R-Q                                                             00012610
      IF (H) 1,1,2                                                      00012620
    1 AREA=0.0                                                          00012630
      DEL=0.0                                                           00012640
      GO TO 10                                                          00012650
C R MUST BE GREATER THAN Q                                              00012660
    2 SP=1.0E35                                                         00012670
      S3=0.0                                                            00012680
      S1=F(Q)+F(R)                                                      00012690
      IF (HBAR) 3,3,4                                                   00012700
    3 HBAR=0.007*H                                                      00012710
    4 S2=0.0                                                            00012720
      X=Q+0.5*H                                                         00012730
    5 S2=S2+4.0*F(X)                                                    00012740
      X=X+H                                                             00012750
      IF (X-R) 5,5,6                                                    00012760
    6 SC=(S1+S2+S3)*H*0.16666667                                        00012770
      DEL=0.066666667*(SP-SC)                                           00012780
      IF (ABS(DEL)-ABS(EPS*SC)) 7,8,8                                   00012790
    7 AREA=SC-DEL                                                       00012800
      GO TO 10                                                          00012810
    8 S3=S3+0.5*S2                                                      00012820
      H=0.5*H                                                           00012830
      IF (H-HBAR) 7,9,9                                                 00012840
    9 SP=SC                                                             00012850
      GO TO 4                                                           00012860
   10 RETURN                                                            00012870
      END                                                               00012880
      SUBROUTINE JY0(X,J0,Y0)                                           00012890
C***********************************************************************00012900
C                                                                       00012910
C     PURPOSE                                                           00012920
C        COMPUTES BESSEL FUNCTIONS OF THE FIRST AND SECOND KIND,        00012930
C        ZERO ORDER, FOR POSITIVE ARGUMENTS.                            00012940
C        SEE NBS AMS 55, P. 369-370.                                    00012950
C     DESCRIPTION OF PARAMETERS - ALL REAL                              00012960
C        X- ARGUMENT, MUST BE >0                                        00012970
C        J0 - RETURNED FUNCTION VALUE, J0(X)                            00012980
C        Y0 - RETURNED FUNCTION VALUE, Y0(X)                            00012990
C                                                                       00013000
C***********************************************************************00013010
      REAL J0                                                           00013020
      IF (X-3.0) 1,2,3                                                  00013030
    1 IF (X) 4,4,2                                                      00013040
    2 Z=(0.33333333*X)**2                                               00013050
      J0=1.0-Z*(2.2499997-Z*(1.2656208-Z*(0.3163866-Z*(0.0444479-Z*(0.0000013060
     139444-0.00021*Z)))))                                              00013070
      Y0=0.63661977*ALOG(0.5*X)*J0+0.36746691+Z*(0.60559366-Z*(0.743503800013080
     14-Z*(0.25300117-Z*(0.04261214-Z*(0.00427916-0.00024846*Z)))))     00013090
      RETURN                                                            00013100
    3 Z=3.0/X                                                           00013110
      F=0.79788456-Z*(0.77E-6+Z*(0.0055274+Z*(0.00009512-Z*(0.00137237-Z00013120
     1*(0.00072805-0.00014476*Z)))))                                    00013130
      P=0.78539816+Z*(0.04166397+Z*(0.00003954-Z*(0.00262573-Z*(0.00054100013140
     125+Z*(0.00029333-0.00013558*Z)))))                                00013150
      Q=SQRT(1.0/X)                                                     00013160
      J0=Q*F*COS(X-P)                                                   00013170
      Y0=Q*F*SIN(X-P)                                                   00013180
    4 RETURN                                                            00013190
      END                                                               00013200
      SUBROUTINE JY1(X,J1,Y1)                                           00013210
C***********************************************************************00013220
C                                                                       00013230
C     PURPOSE                                                           00013240
C        COMPUTES BESSEL FUNCTIONS OF THE FIRST AND SECOND KIND,        00013250
C        FIRST ORDER, FOR POSITIVE ARGUMENTS.                           00013260
C        SEE NBS AMS 55, P. 370.                                        00013270
C     DESCRIPTION OF PARAMETERS - ALL REAL                              00013280
C        X- ARGUMENT, MUST BE >0                                        00013290
C        J1 - RETURNED FUNCTION VALUE, J1(X)                            00013300
C        Y1 - RETURNED FUNCTION VALUE, Y1(X)                            00013310
C                                                                       00013320
C***********************************************************************00013330
      REAL J1                                                           00013340
      IF (X-3.0) 1,2,3                                                  00013350
    1 IF (X) 4,4,2                                                      00013360
    2 Z=(0.33333333*X)**2                                               00013370
      J1=X*(0.5-Z*(0.56249985-Z*(0.21093573-Z*(0.03954289-Z*(0.00443319-00013380
     1Z*(0.00031761-0.00001109*Z))))))                                  00013390
      Y1=0.63661977*ALOG(0.5*X)*J1+(-0.6366198+Z*(0.2212091+Z*(2.168270900013400
     1-Z*(1.3164827-Z*(0.3123951-Z*(0.0400976-0.0027873*Z))))))/X       00013410
      RETURN                                                            00013420
    3 Z=3.0/X                                                           00013430
      F=0.79788456+Z*(0.156E-5+Z*(0.01659667+Z*(0.00017105-Z*(0.0024951100013440
     1-Z*(0.00113653-0.00020033*Z)))))                                  00013450
      P=0.78539816-Z*(0.12499612+Z*(0.0000565-Z*(0.00637879-Z*(0.000743400013460
     18+Z*(0.00079824-0.00029166*Z)))))                                 00013470
      Q=SQRT(1.0/X)                                                     00013480
      J1=Q*F*SIN(X-P)                                                   00013490
      Y1=-Q*F*COS(X-P)                                                  00013500
    4 RETURN                                                            00013510
      END                                                               00013520
C***********************************************************************00013530
C                                                                       00013540
C     PURPOSE                                                           00013550
C        COMPUTES FUNCTION VALUES OF F(UW,ALPHA)  -                     00013560
C        PAPADOPULOS,I.S. AND COOPER,H.H.,JR., 1967, DRAWDOWN IN        00013570
C        A WELL OF LARGE DIAMETER: WATER RESOURCES RESEARCH, V. 3,      00013580
C        NO. 1, P. 241-244.                                             00013590
C        PROGRAM BY S.S.PAPADOPULOS.                                    00013600
C     INPUT DATA - ONE OR MORE GROUPS, EACH GROUP CODED AS FOLLOWS      00013610
C        1 CARD - FORMAT (E10.5)                                        00013620
C             S - (ALPHA) - RW**2*S/RC**2 - RADIUS OF WELL (SCREEN      00013630
C                  OR OPEN BORE IN AQUIFER) SQUARED * STORAGE           00013640
C                  COEFFICIENT / RADIUS OF CASING (OVER INTERVAL OF     00013650
C                  WATER LEVEL CHANGE) SQUARED.                         00013660
C        1 CARD - FORMAT(16E5.0)                                        00013670
C             U- 16 VALUES OF UW - RW**2*S/(4*T*TIME) - RADIUS OF       00013680
C                  PUMPED WELL SQUARED * STORAGE COEFFICIENT /          00013690
C                  4 * TRANSMISSIVITY * TIME. IF LESS THAN 16 DESIRED,  00013700
C                  BLANK OR ZERO VALUES MAY BE CODED FOR THE REST.      00013710
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00013720
C        PEAK,SIMP,APEKE,EXBSL2,JY0,JY1 - MUST BE INCLUDED IN DECK.     00013730
C                                                                       00013740
C***********************************************************************00013750
      COMMON XPK,YPK                                                    00013760
      COMMON/PBLK/A,B                                                   00013770
      EXTERNAL EXBSL2                                                   00013780
      DIMENSION U(16)                                                   00013790
      EPS=0.0001                                                        00013800
    1 READ (5,13,END=12) S                                              00013810
      IF (S) 1,1,2                                                      00013820
    2 READ (5,14) U                                                     00013830
      WRITE (6,15) S                                                    00013840
      DO 11 I=1,16                                                      00013850
      UW=U(I)                                                           00013860
      IF (UW) 1,1,3                                                     00013870
    3 B=0.25/UW                                                         00013880
      A=S+S                                                             00013890
      CALL APEKE(EXBSL2)                                                00013900
      CALL PEAK(EXBSL2)                                                 00013910
      IF (XPK-1.0E-8) 4,5,5                                             00013920
    4 WRITE (6,16) UW,S,XPK,YPK                                         00013930
      GO TO 11                                                          00013940
    5 IF (XPK-1.0E8) 7,7,6                                              00013950
    6 WRITE (6,17) UW,S,XPK,YPK                                         00013960
      GO TO 11                                                          00013970
    7 HBAR=0.007*XPK                                                    00013980
      CALL SIMPS(0.0,XPK,EPS,HBAR,SUM,DEL,EXBSL2)                       00013990
      X2=XPK                                                            00014000
      DX=XPK                                                            00014010
    8 DX=10.0*DX                                                        00014020
      X1=X2                                                             00014030
      X2=X1+DX                                                          00014040
      Y=EXBSL2(X2)                                                      00014050
      HBAR=0.007*DX                                                     00014060
      CALL SIMPS(X1,X2,EPS,HBAR,TRM,ERR,EXBSL2)                         00014070
      SUM=SUM+TRM                                                       00014080
      DEL=DEL+ERR                                                       00014090
      IF (X2-1.0E9) 9,10,10                                             00014100
    9 YT=1.5707963/X2**4                                                00014110
      IF (ABS(Y-YT)/YT-0.5E-6) 10,8,8                                   00014120
   10 EST=0.52359878/X2**3                                              00014130
      SUM=SUM+EST                                                       00014140
      FUWS=3.2422779*S*S*SUM                                            00014150
      WRITE (6,18) UW,SUM,DEL,FUWS,XPK,YPK                              00014160
   11 CONTINUE                                                          00014170
      GO TO 1                                                           00014180
   12 STOP                                                              00014190
C                                                                       00014200
   13 FORMAT (E10.5)                                                    00014210
   14 FORMAT (16E5.0)                                                   00014220
   15 FORMAT ('1','F(UW,ALPHA) FOR ALPHA=',1PE14.5/'0',7X,'UW',12X,'INTE00014230
     1GRAL',5X,'INTEGRAL ERROR',5X,'F(UW,ALPHA)',8X,'X(PEAK)',10X,'Y(PEA00014240
     2K)'/' ')                                                          00014250
   16 FORMAT (1H ,1PE14.7,9X,34HVALUES OF DUMMY VARIABLE TOO SMALL,1PE2500014260
     1.7,1PE17.7)                                                       00014270
   17 FORMAT (1H ,1PE14.7,9X,34HVALUES OF DUMMY VARIABLE TOO LARGE,1PE2500014280
     1.7,1PE17.7)                                                       00014290
   18 FORMAT (1H ,1PE14.5,1P5E17.5)                                     00014300
      END                                                               00014310
      FUNCTION EXBSL2(X)                                                00014320
C***********************************************************************00014330
C                                                                       00014340
C     PURPOSE                                                           00014350
C        COMPUTES VALUES OF THE INTEGRAND FOR F(UW,ALPHA)               00014360
C     DESCRIPTION OF PARAMETER                                          00014370
C        X- REAL - ARGUMENT OF INTEGRAND                                00014380
C                                                                       00014390
C***********************************************************************00014400
      COMMON/PBLK/A,B                                                   00014410
      IF (X) 1,1,2                                                      00014420
    1 EXBSL2=0.                                                         00014430
      GO TO 8                                                           00014440
    2 IF (X-1.E+7) 4,4,3                                                00014450
    3 EXBSL2=1.5707963/X**4                                             00014460
      GO TO 8                                                           00014470
    4 Y=B*X*X                                                           00014480
      IF (Y-.01) 5,5,6                                                  00014490
    5 FNUM=Y*(1.-Y*(.5-Y*((1./6.)-Y*(1./24.))))                         00014500
      GO TO 7                                                           00014510
    6 FNUM=1.-EXP(-Y)                                                   00014520
    7 CALL JY0(X,BJ0,BY0)                                               00014530
      CALL JY1(X,BJ1,BY1)                                               00014540
      DEN=((X*BJ0-A*BJ1)**2+(X*BY0-A*BY1)**2)*X**3                      00014550
      EXBSL2=FNUM/DEN                                                   00014560
    8 RETURN                                                            00014570
      END                                                               00014580
C********************************************************************** 00014590
C                                                                       00014600
C     PURPOSE                                                           00014610
C        COMPUTES CHANGES IN WATER LEVEL, H(R,T), IN RESPONSE TO        00014620
C        VARYING DISCHARGE USING THE CONVOLUTION INTEGRAL FOR           00014630
C        LEAKY AQUIFERS - EQ. 3 OF MOENCH,ALLEN,1971, GROUND-WATER      00014640
C        FLUCTUATIONS IN RESPONSE TO ARBITRARY PUMPAGE: GROUND WATER,   00014650
C        V.9, NO.2,P.4-8.                                               00014660
C     INPUT DATA - ONE OR MORE GROUPS, EACH GROUP CODED AS FOLLOWS      00014670
C        1 CARD - FORMAT(2E10.5,4X,I1,5X,E10.5)                         00014680
C             TBEGIN - SMALLEST VALUE OF TIME FOR OUTPUT.               00014690
C             TEND - LARGEST VALUE OF TIME FOR OUTPUT.                  00014700
C             IQ - INDICATES FORM OF DISCHARGE FUNCTION, Q(T).          00014710
C                  IQ=1,2,3 REFER TO DISCHARGE FUNCTIONS IN             00014720
C                  HANTUSH,M.S., 1964, HYDRAULICS OF WELLS IN CHOW,     00014730
C                  VEN TE, ED., ADVANCES IN HYDROSCIENCE, VOL. 1:       00014740
C                  ACADEMIC PRESS INC., NEW YORK, P. 281-442.           00014750
C                  IQ=1, Q(T) IS AN EXPONENTIAL FUNCTION, CASE A,       00014760
C                       P. 343 OF HANTUSH.                              00014770
C                  IQ=2, Q(T) IS A HYPERBOLIC FUNCTION, CASE B,         00014780
C                       P. 344 OF HANTUSH.                              00014790
C                  IQ=3, Q(T) IS AN INVERSE SQUARE ROOT FUNCTION,       00014800
C                       CASE C, P. 344 OF HANTUSH.                      00014810
C                  IQ=4, Q(T) IS A FIFTH-DEGREE POLYNOMIAL.             00014820
C                  IQ=5, Q(T) IS A PIECEWISE LINEAR FUNCTION OF         00014830
C                       TIME (EIGHT SEGMENTS).                          00014840
C             QR - REFERENCE DISCHARGE, ZERO OR BLANK FOR PROJECTION,   00014850
C        1 OR 4 CARDS, DEPENDING ON IQ.                                 00014860
C             IF IQ=1,2,3 - 1 CARD - FORMAT(3E10.3)                     00014870
C                  QST - EVENTUAL CONSTANT DISCHARGE.                   00014880
C                  DELTA - RATE PARAMETER.                              00014890
C                  TSTAR - TIME PARAMETER.                              00014900
C             IF IQ=4 - 1 CARD - FORMAT(6E10.3)                         00014910
C                  AQ(6) - 6 VALUES - THE POLYNOMIAL COEFFICIENTS       00014920
C                       WITH A0 FIRST AND A5 LAST.                      00014930
C             IF IQ=5 - 4 CARDS - FORMAT(6E10.3)                        00014940
C                  TI(I),AI(I),BI(I),TI(I+1),AI(I+1),BI(I+1),I=1,3,5,7  00014950
C                  PARAMETERS OF THE PIECEWISE LINEAR FUNCTION          00014960
C                  (8 SEGMENTS), CODED 2 SEGMENTS PER CARD. FIRST       00014970
C                  AND SECOND SEGMENTS ON FIRST CARD, THEN SEQUENTIALLY 00014980
C                  ON SUCCEEDING CARDS. EACH SEGMENT HAS THREE          00014990
C                  PARAMETERS WHICH ARE IN CODING ORDER                 00015000
C                       TI - ENDING TIME OF THE SEGMENT.                00015010
C                       AI - DISCHARGE AT BEGINNING OF SEGMENT.         00015020
C                       BI - RATE OF CHANGE IN DISCHARGE DURING SEG.    00015030
C                  THE DISCHARGE FUNCTION IN EACH SEGMENT HAS THE       00015040
C                  FORM Q(T) = AI(I)+BI(I)*(T-TI(I-1)). IF LESS THAN 8  00015050
C                  SEGMENTS ARE NEEDED, BLANKS CAN BE CODED FOR         00015060
C                  SUCCEEDING SEGMENTS.                                 00015070
C        2 OR MORE CARDS - FORMAT(4E10.3)                               00015080
C             R - RADIAL DISTANCE FROM PUMPED WELL, BLANK OR ZERO       00015090
C                  SIGNALS PROGRAM AS END TO GROUP OF DATA.             00015100
C             S - STORAGE COEFFICIENT                                   00015110
C             T - TRANSMISSIVITY                                        00015120
C             PM - (P'/M') - HYD. COND. OF CONFINING BED DIVIDED        00015130
C                  BY THICKNESS OF CONFINING BED.                       00015140
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00015150
C        CONVOL,Q - MUST BE INCLUDED IN DECK.                           00015160
C                                                                       00015170
C********************************************************************** 00015180
      DIMENSION       D(12),IEX(12),X(6),H(12,6),QS(12,6),CP(12),CT(12) 00015190
      DIMENSION H1(12),H2(12),Q1(12),Q2(12)                             00015200
      DIMENSION H3(12),H4(12),Q3(12),Q4(12)                             00015210
      COMMON AQ(6),TI(9),AI(9),BI(9),QST,DELTA,TSTAR                    00015220
      DATA CP/12*'  T*'/,CT/12*'1/U*'/,D/12*'10**'/                     00015230
      DATA H1/12*'  S('/,H2/12*'R,T)'/,Q1/12*'    '/,Q2/12*'Q(T)'/      00015240
      DATA H3/12*'   S'/,H4/12*'O(T)'/,Q3/12*' Q(T'/,Q4/12*')/QR'/      00015250
      DATA X/1.,1.5,2.,3.,5.,7./                                        00015260
      TI(1)=0.                                                          00015270
      N=500                                                             00015280
    1 READ (5,18,END=17) TBEGIN,TEND,IQ,QR                              00015290
      IF (IQ.LT.4) READ (5,19) QST,DELTA,TSTAR                          00015300
      IF (IQ.EQ.4) READ (5,19) AQ                                       00015310
      IF (IQ.EQ.5) READ (5,19) (TI(I),AI(I),BI(I),I=2,9)                00015320
      WRITE (6,24)                                                      00015330
    2 READ (5,19) R,S,T,PM                                              00015340
      IF (R.EQ.0.) GO TO 1                                              00015350
      A=R*R*S/(4.*T)                                                    00015360
      B=PM/S                                                            00015370
      Y=ALOG10(TBEGIN)                                                  00015380
      IF (Y) 3,5,4                                                      00015390
    3 Y=Y-.001                                                          00015400
      GO TO 5                                                           00015410
    4 Y=Y+.001                                                          00015420
    5 IBEGIN=Y                                                          00015430
      Y=ALOG10(TEND)                                                    00015440
      IF (Y) 6,8,7                                                      00015450
    6 Y=Y-.001                                                          00015460
      GO TO 8                                                           00015470
    7 Y=Y+.001                                                          00015480
    8 IEND=Y                                                            00015490
      M=IEND-IBEGIN+1                                                   00015500
      IF (M.GT.12) M=12                                                 00015510
      DO 10 I=1,M                                                       00015520
      IEX(I)=IBEGIN+I-1                                                 00015530
      Y=10.**(IBEGIN+I-1)                                               00015540
      DO 10 J=1,6                                                       00015550
      TIME=X(J)*Y                                                       00015560
      IF (QR.GT.0.) TIME=A*TIME                                         00015570
      CALL CONVOL(TIME,A,B,N,IQ,SUM)                                    00015580
      IF (QR.GT.0.) GO TO 9                                             00015590
      H(I,J)=SUM/(12.5664*T)                                            00015600
      QS(I,J)=Q(TIME,IQ)                                                00015610
      GO TO 10                                                          00015620
    9 H(I,J)=SUM/QR                                                     00015630
      QS(I,J)=Q(TIME,IQ)/QR                                             00015640
   10 CONTINUE                                                          00015650
      K=M                                                               00015660
      IF (M.GT.6) K=6                                                   00015670
      IF (QR.GT.0.) GO TO 11                                            00015680
      WRITE (6,20) A,B,(CP(I),D(I),IEX(I),I=1,K)                        00015690
      WRITE (6,21) (H1(I),H2(I),Q1(I),Q2(I),I=1,K)                      00015700
      GO TO 12                                                          00015710
   11 WRITE (6,25) A,B,QR,(CT(I),D(I),IEX(I),I=1,K)                     00015720
      WRITE (6,21) (H3(I),H4(I),Q3(I),Q4(I),I=1,K)                      00015730
   12 DO 13 J=1,6                                                       00015740
      WRITE (6,22) X(J),(H(I,J),QS(I,J),I=1,K)                          00015750
   13 CONTINUE                                                          00015760
      IF (M.LE.6) GO TO 2                                               00015770
      K1=K+1                                                            00015780
      IF (QR.GT.0.) GO TO 14                                            00015790
      WRITE (6,23) (CP(I),D(I),IEX(I),I=K1,M)                           00015800
      WRITE (6,21) (H1(I),H2(I),Q1(I),Q2(I),I=K1,M)                     00015810
      GO TO 15                                                          00015820
   14 WRITE (6,26) (CT(I),D(I),IEX(I),I=K1,M)                           00015830
      WRITE (6,21) (H3(I),H4(I),Q3(I),Q4(I),I=K1,M)                     00015840
   15 DO 16 J=1,6                                                       00015850
      WRITE (6,22) X(J),(H(I,J),QS(I,J),I=K1,M)                         00015860
   16 CONTINUE                                                          00015870
      GO TO 2                                                           00015880
   17 STOP                                                              00015890
C                                                                       00015900
   18 FORMAT (2E10.5,4X,I1,5X,E10.5)                                    00015910
   19 FORMAT (6E10.3)                                                   00015920
   20 FORMAT ('0','R**2*S/(4*TRANS)=',1PE10.3,', K''/(S*B'')=',E10.3/'0'00015930
     1,2X,'T',5X,6(2A4,I2,9X))                                          00015940
   21 FORMAT (' ',4X,6(2A4,2X,2A4,1X))                                  00015950
   22 FORMAT (' ',F4.1,6(0PF8.3,1PE11.3))                               00015960
   23 FORMAT ('0',2X,'T',5X,6(2A4,I2,9X))                               00015970
   24 FORMAT (1H1)                                                      00015980
   25 FORMAT ('0','R**2*S/(4*TRANS)=',1PE10.3,', K''/(S*B'')=',E10.3,', 00015990
     1QR=',E10.3/'0',1X,'1/U',4X,6(2A4,I2,9X))                          00016000
   26 FORMAT ('0',1X,'1/U',4X,6(2A4,I2,9X))                             00016010
      END                                                               00016020
C********************************************************************** 00016030
C                                                                       00016040
C     PURPOSE                                                           00016050
C        COMPUTES FUNCTION VALUES OF F(BETA,ALPHA) - THE SLUG TEST      00016060
C        FUNCTION - COOPER,H.H.,JR., BREDEHOEFT,J.D., AND PAPADOPULOS,  00016070
C        I.S., 1967, RESPONSE OF A FINITE-DIAMETER WELL TO AN           00016080
C        INSTANTANEOUS CHARGE OF WATER: WATER RESOURCES RESEARCH,       00016090
C        V. 3, NO. 1, P. 263-269.                                       00016100
C        PROGRAM BY S.S.PAPADOPULOS.                                    00016110
C     INPUT DATA                                                        00016120
C        1 OR MORE CARDS - FORMAT(F16.5)                                00016130
C             A - (ALPHA) - RW**2*S/RC**2 - RADIUS OF WELL (SCREEN OR   00016140
C                  OPEN BORE IN AQUIFER) SQUARED * STORAGE COEFFICIENT  00016150
C                  / RADIUS OF CASING (OVER INTERVAL OF WATER LEVEL     00016160
C                  CHANGE) SQUARED.                                     00016170
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00016180
C        PRX,DJY0,DJY1,DSIMPS - MUST BE INCLUDED IN DECK                00016190
C     METHOD                                                            00016200
C        THIS PROGRAM CALCULATES THE SLUG TEST FUNCTION, F(BETA,ALPHA), 00016210
C        FOR VALUES OF BETA RANGING FROM 0.001 TO 1000.0 BY INCREMENTING00016220
C        BETA ACCORDING TO DATA ARRAY BB(I). AVERAGE COMPUTATION TIME   00016230
C        IS ABOUT 30 SECONDS PER VALUE OF ALPHA ON IBM 360/155.         00016240
C                                                                       00016250
C********************************************************************** 00016260
      DOUBLE PRECISION A,B,PI,ZZ,EPS,Y,X1,X2,TERM,FAB,DATAN,DEL,HBAR    00016270
      DIMENSION ZZ(40), BB(39)                                          00016280
      COMMON A,B,PI                                                     00016290
      EXTERNAL PRX                                                      00016300
      DATA ZZ/0.D+0,1.D-10,1.D-9,1.D-8,1.D-7,1.D-6,1.D-5,1.D-4,         00016310
     1   1.D-3,1.D-2,1.D-1,2.D-1,3.D-1,4.D-1,5.D-1,6.D-1,7.D-1,8.D-1,   00016320
     2   9.D-1,1.D+0,2.D+0,3.D+0,4.D+0,5.D+0,6.D+0,7.D+0,8.D+0,         00016330
     3   9.D+0,1.D+1,2.D+1,3.D+1,4.D+1,5.D+1,6.D+1,7.D+1,8.D+1,         00016340
     4   9.D+1,1.D+2,1.25D+2,1.5D+2/                                    00016350
      DATA BB/.001,.002,.004,.006,.008,.01,.02,.04,.06,.08,.1,.2,.4,.6,.00016360
     18,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100016370
     200.,200.,400.,600.,800.,1000./                                    00016380
      PI=4.*DATAN(1.0D+00)                                              00016390
      EPS=0.00001                                                       00016400
    1 READ (5,6) A                                                      00016410
      IF (A.LE.0.0) GO TO 5                                             00016420
      WRITE (6,7) A                                                     00016430
      WRITE (6,8)                                                       00016440
      DO 4 I=1,39                                                       00016450
      B=BB(I)                                                           00016460
      Y=0.0                                                             00016470
      DO 2 L=1,39                                                       00016480
      X1=ZZ(L)                                                          00016490
      X2=ZZ(L+1)                                                        00016500
      HBAR=0.                                                           00016510
      CALL DSIMPS(X1,X2,EPS,HBAR,TERM,DEL,PRX)                          00016520
      Y=Y+TERM                                                          00016530
      IF (L.GT.20.AND.TERM.LT.EPS) GO TO 3                              00016540
    2 CONTINUE                                                          00016550
    3 FAB=4.*A*Y/(PI*PI)                                                00016560
    4 WRITE (6,9) B,FAB                                                 00016570
      GO TO 1                                                           00016580
    5 STOP                                                              00016590
C                                                                       00016600
C                                                                       00016610
    6 FORMAT (F16.5)                                                    00016620
    7 FORMAT ('1',41X,'F(BETA,ALPHA) FOR ALPHA=',1PD9.2/)               00016630
    8 FORMAT ('0',53X,'BETA',13X,'H/H0'/)                               00016640
    9 FORMAT (' ',51X,1PD8.2,10X,0PF6.4)                                00016650
      END                                                               00016660
      DOUBLE PRECISION FUNCTION PRX(X)                                  00016670
C********************************************************************** 00016680
C                                                                       00016690
C     PURPOSE                                                           00016700
C        COMPUTE VALUES OF THE INTEGRAND FOR F(BETA,ALPHA)              00016710
C     DESCRIPTION OF PARAMETER                                          00016720
C        X - DOUBLE PRECISION - ARGUMENT OF INTEGRAND                   00016730
C                                                                       00016740
C********************************************************************** 00016750
      DOUBLE PRECISION A,B,PI,XX,X,C,F1,F2,J0,Y0,J1,Y1                  00016760
      DOUBLE PRECISION DLOG,DSQRT,DEXP                                  00016770
      COMMON A,B,PI                                                     00016780
      XX=DSQRT(A*X/B)                                                   00016790
      IF (X) 6,1,2                                                      00016800
    1 PRX=(PI*PI)/(16.*A*B)                                             00016810
      GO TO 6                                                           00016820
    2 IF (X.LT.150.) GO TO 3                                            00016830
      PRX=0.0                                                           00016840
      GO TO 6                                                           00016850
    3 IF (XX.GT.0.0001) GO TO 4                                         00016860
      C=DEXP(5.772156649D-01)/2.                                        00016870
      F1=PI*X*(1.-A)                                                    00016880
      F2=X*DLOG(C*C*A*X/B)+4.*B                                         00016890
      PRX=(B*PI*PI*DEXP(-X))/(A*(F1*F1+F2*F2))                          00016900
      GO TO 6                                                           00016910
    4 IF (XX.LT.50.) GO TO 5                                            00016920
      PRX=(PI*DEXP(-X))/(2.*XX*(X+4.*A*B))                              00016930
      GO TO 6                                                           00016940
    5 CALL DJY0(XX,J0,Y0)                                               00016950
      CALL DJY1(XX,J1,Y1)                                               00016960
      F1=(XX*J0-2.*A*J1)                                                00016970
      F2=(XX*Y0-2.*A*Y1)                                                00016980
      PRX=DEXP(-X)/(X*(F1*F1+F2*F2))                                    00016990
    6 RETURN                                                            00017000
      END                                                               00017010
      SUBROUTINE DJY0(X,J0,Y0)                                          00017020
C***********************************************************************00017030
C                                                                       00017040
C     PURPOSE                                                           00017050
C        COMPUTES BESSEL FUNCTIONS OF THE FIRST AND SECOND KIND,        00017060
C        ZERO ORDER, FOR POSITIVE ARGUMENTS.                            00017070
C     DESCRIPTION OF PARAMETERS - ALL DOUBLE PRECISION                  00017080
C        X- ARGUMENT, MUST BE >0                                        00017090
C        J0 - RETURNED FUNCTION VALUE, J0(X)                            00017100
C        Y0 - RETURNED FUNCTION VALUE, Y0(X)                            00017110
C                                                                       00017120
C***********************************************************************00017130
      DOUBLE PRECISION Z,J0,Y0,F,P,Q,U,W,X,DLOG,DCOS,DSIN,DSQRT         00017140
      IF (X-3.0) 1,2,3                                                  00017150
    1 IF (X) 4,4,2                                                      00017160
    2 Z=(X/3.0)**2                                                      00017170
      J0=1.0-Z*(2.2499997-Z*(1.2656208-Z*(0.3163866-Z*(0.0444479-Z*(0.0000017180
     139444-0.00021*Z)))))                                              00017190
      W=(0.5D0)*X                                                       00017200
      Y0=0.63661977*DLOG(W)*J0+0.36746691+Z*(0.60559366-Z*(0.74350384-Z*00017210
     1(0.25300117-Z*(0.04261214-Z*(0.00427916-0.00024846*Z)))))         00017220
      RETURN                                                            00017230
    3 Z=3.0/X                                                           00017240
      F=0.79788456-Z*(0.77D-6+Z*(0.0055274+Z*(0.00009512-Z*(0.00137237-Z00017250
     1*(0.00072805-0.00014476*Z)))))                                    00017260
      P=0.78539816+Z*(0.04166397+Z*(0.00003954-Z*(0.00262573-Z*(0.00054100017270
     125+Z*(0.00029333-0.00013558*Z)))))                                00017280
      U=(1.0D0)/X                                                       00017290
      Q=DSQRT(U)                                                        00017300
      J0=Q*F*DCOS(X-P)                                                  00017310
      Y0=Q*F*DSIN(X-P)                                                  00017320
    4 RETURN                                                            00017330
      END                                                               00017340
      SUBROUTINE DJY1(X,J1,Y1)                                          00017350
C***********************************************************************00017360
C                                                                       00017370
C     PURPOSE                                                           00017380
C        COMPUTES BESSEL FUNCTIONS OF THE FIRST AND SECOND KIND,        00017390
C        FIRST ORDER, FOR POSITIVE ARGUMENTS.                           00017400
C     DESCRIPTION OF PARAMETERS - ALL DOUBLE PRECISION                  00017410
C        X- ARGUMENT, MUST BE >0                                        00017420
C        J1 - RETURNED FUNCTION VALUE, J1(X)                            00017430
C        Y1 - RETURNED FUNCTION VALUE, Y1(X)                            00017440
C                                                                       00017450
C***********************************************************************00017460
      DOUBLE PRECISION X,J1,Y1,Z,W,DLOG,F,P,U,Q,DSQRT,DSIN,DCOS         00017470
      IF (X-3.0) 1,2,3                                                  00017480
    1 IF (X) 4,4,2                                                      00017490
    2 Z=(X/3.0)**2                                                      00017500
      J1=X*(0.5-Z*(0.56249985-Z*(0.21093573-Z*(0.03954289-Z*(0.00443319-00017510
     1Z*(0.00031761-0.00001109*Z))))))                                  00017520
      W=(0.5D0)*X                                                       00017530
      Y1=0.63661977*DLOG(W)*J1+(-0.6366198+Z*(0.2212091+Z*(2.1682709-Z*(00017540
     11.3164827-Z*(0.3123951-Z*(0.0400976-0.0027873*Z))))))/X           00017550
      RETURN                                                            00017560
    3 Z=3.0/X                                                           00017570
      F=0.79788456+Z*(0.156D-5+Z*(0.01659667+Z*(0.00017105-Z*(0.0024951100017580
     1-Z*(0.00113653-0.00020033*Z)))))                                  00017590
      P=0.78539816-Z*(0.12499612+Z*(0.0000565-Z*(0.00637879-Z*(0.000743400017600
     18+Z*(0.00079824-0.00029166*Z)))))                                 00017610
      U=(1.0D0)/X                                                       00017620
      Q=DSQRT(U)                                                        00017630
      J1=Q*F*DSIN(X-P)                                                  00017640
      Y1=-Q*F*DCOS(X-P)                                                 00017650
    4 RETURN                                                            00017660
      END                                                               00017670
      SUBROUTINE DSIMPS(A,B,EPS,HBAR,AREA,DEL,F)                        00017680
C***********************************************************************00017690
C                                                                       00017700
C     PURPOSE                                                           00017710
C        TO DETERMINE THE INTEGRAL OF A FUNCTION, F, FROM A TO B,       00017720
C        USING SIMPSON'S RULE.                                          00017730
C     DESCRIPTION OF PARAMETERS                                         00017740
C        ALL DOUBLE PRECISION                                           00017750
C        A - LOWER LIMIT OF INTEGRAL                                    00017760
C        B - UPPER LIMIT OF INTEGRAL                                    00017770
C        EPS - DESIRED ACCURACY                                         00017780
C        HBAR - MINIMUM DIVISION OF THE INTERVAL                        00017790
C        AREA - COMPUTED VALUE OF INTEGRAL BETWEEN Q AND R              00017800
C        DEL - COMPUTED ESTIMATE OF ERROR                               00017810
C        F- THE INTEGRAND (FUNCTION REFERENCE)                          00017820
C     METHOD                                                            00017830
C        USES SIMPSON'S RULE TO COMPUTE A SUM APPROXIMATING THE INTEGRAL00017840
C        USES INITIAL H=(B-A)/2. COMPUTES A SEQUENCE OF SUMS BY HALVING 00017850
C        H EACH TIME. COMPUTES ESTIMATE OF ERROR (DEL) AS (PREVIOUS     00017860
C        SUM - CURRENT SUM)/15. COMPUTATION STOPS WHEN 1) H<HBAR,       00017870
C        2) ABS(DEL)<ABS(EPS*CURRENT SUM). IF HBAR IS LE 0,             00017880
C        THEN HBAR=.007*(B-A).                                          00017890
C                                                                       00017900
C***********************************************************************00017910
      DOUBLE PRECISION H,HBAR,AREA,DEL,S1,S2,S3,SC,SP,X,A,B,EPS,F,DABS  00017920
C  AREA OF F FROM A TO B,EPS IS DESIRED ACCURACY, HBAR THE MINIMUM      00017930
C  ALLOWABLE INTERVAL, DEL THE ESTIMATE OF THE ERROR                    00017940
      H=B-A                                                             00017950
      IF (H) 1,1,2                                                      00017960
    1 AREA=0.0                                                          00017970
      DEL=0.0                                                           00017980
      GO TO 10                                                          00017990
    2 SP=1.0D35                                                         00018000
      S3=0.0                                                            00018010
      S1=F(A)+F(B)                                                      00018020
      IF (HBAR) 3,3,4                                                   00018030
    3 HBAR=0.007*H                                                      00018040
    4 S2=0.0                                                            00018050
      X=A+0.5*H                                                         00018060
    5 S2=S2+4.0*F(X)                                                    00018070
      X=X+H                                                             00018080
      IF (X-B) 5,5,6                                                    00018090
    6 SC=(S1+S2+S3)*H*0.16666666667                                     00018100
      DEL=0.06666666667*(SP-SC)                                         00018110
      IF (DABS(DEL)-DABS(EPS*SC)) 7,8,8                                 00018120
    7 AREA=SC-DEL                                                       00018130
      GO TO 10                                                          00018140
    8 S3=S3+0.5*S2                                                      00018150
      H=0.5*H                                                           00018160
      IF (H-HBAR) 7,9,9                                                 00018170
    9 SP=SC                                                             00018180
      GO TO 4                                                           00018190
   10 RETURN                                                            00018200
      END                                                               00018210
      REAL FUNCTION L*4(U,V)                                            00018220
C     ******************************************************************00018230
C                                                                       00018240
C     FUNCTION L                                                        00018250
C                                                                       00018260
C     PURPOSE                                                           00018270
C        TO COMPUTE THE INTEGRAL( EXP(-Y-V**2/Y)/Y) SUMMED OVER Y FROM  00018280
C        U TO INFINITY(WELL FUNCTION FOR LEAKY AQUIFERS).               00018290
C     DESCRIPTION OF PARAMETERS                                         00018300
C        BOTH DOUBLE PRECISION                                          00018310
C        U - R**2*S/4*T*TIME (RADIAL DISTANCE SQUARED * STORAGE         00018320
C            COEFFICIENT / 4*TRANSMISSIVITY * TIME                      00018330
C        V - R/2*SQRT(K'/(T*B'))--ONE-HALF RADIAL DISTANCE*SQUARE ROOT  00018340
C            (HYD. COND. OF CONFINING BED/TRANSMISSIVITY*THICKNESS      00018350
C            OF CONFINING BED)                                          00018360
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00018370
C        DQL12,SERIES,BESK,FCT                                          00018380
C     METHOD                                                            00018390
C        IN THE FOLLOWING F=EXP(-Y-V**2/Y)/Y                            00018400
C      (1) U>=1, USES A GAUSSIAN-LAGUERRE QUADRATURE FORMULA TO         00018410
C        EVALUATE INTEGRAL(F) FROM U TO INF.                            00018420
C      (2) V**2<U<1, USES THE G-L QUADRATURE TO EVALUATE INTEGRAL(F)    00018430
C        FROM ONE TO INF AND A SERIES EXPANSION TO EVALUATE INTEGRAL(F) 00018440
C        FROM U TO ONE.                                                 00018450
C      (3) U<1, U<=V**2, USES THE REPRESENTATION INTEGRAL(F) FROM U     00018460
C        TO INF. = 2*K0(2*V)-INTEGRAL(F) FROM V**2/U TO INF.            00018470
C        EVALUATES THE ZERO ORDER MODIFIED BESSEL FUNCTION OF SECOND    00018480
C        KIND WITH IBM SUBROUTINE, EVALUATES INTEGRAL BY G-L QUAD.      00018490
C                                                                       00018500
C     ******************************************************************00018510
      EXTERNAL FCT                                                      00018520
      REAL*8 U,V,Z,F,VV,SERIES                                          00018530
      COMMON /C1/ VV,Z                                                  00018540
      VV=V                                                              00018550
      IF (U-1.) 1,2,2                                                   00018560
C     CHECKS IF U<1                                                     00018570
    1 Z=V*V/U                                                           00018580
      IF (Z-1.) 3,4,4                                                   00018590
C     CHECKS IF V**2/U < 1                                              00018600
    2 Z=U                                                               00018610
      CALL DQL12(FCT,F)                                                 00018620
      L=F                                                               00018630
C     INTEGRAL U TO INF. EVALUATED BY GAUSS-LAGUERRE QUADRATURE         00018640
      GO TO 5                                                           00018650
    3 Z=1.                                                              00018660
      CALL DQL12(FCT,F)                                                 00018670
      L=F+SERIES(U,V)                                                   00018680
C     INTEGRAL 1 TO INF. BY G-L QUAD., INTEGRAL U TO 1 BY SERIES EXP.   00018690
      GO TO 5                                                           00018700
    4 TWOV=2.*V                                                         00018710
      CALL BESK(TWOV,0,BK,IER)                                          00018720
      CALL DQL12(FCT,F)                                                 00018730
      L=2.*BK-F                                                         00018740
C     2K0(2V)-INTEGRAL V**2/U TO INF.                                   00018750
    5 RETURN                                                            00018760
      END                                                               00018770
      REAL FUNCTION SERIES*8(U,V)                                       00018780
C     ******************************************************************00018790
C                                                                       00018800
C     FUNCTION SERIES                                                   00018810
C                                                                       00018820
C     PURPOSE                                                           00018830
C        TO EVALUATE S(1)-S(U), WHERE S IS A SERIES EXPANSION OF        00018840
C        INTEGRAL(EXP(-Y-V**2/Y)DY/Y) GIVEN BY: S= SUM, M=0 TO INFINITY,00018850
C        (F(M)*SUM, N=0 TO INF.,(V**(2*N)/((N!)*(M+N)!)) WHERE F(M)=    00018860
C        LOG(U) IF M=0 AND = ((-1)**M/M)*(U**M-(V**2/U)**M) IF M>0.     00018870
C     DESCRIPTION OF PARAMETERS                                         00018880
C        BOTH DOUBLE PRECISION                                          00018890
C        U - R**2*S/4*T*TIME (RADIAL DISTANCE SQUARED * STORAGE         00018900
C            COEFFICIENT / 4*TRANSMISSIVITY * TIME                      00018910
C        V - R/2*SQRT(K'/(T*B'))--ONE-HALF RADIAL DISTANCE*SQUARE ROOT  00018920
C            (HYD. COND. OF CONFINING BED/TRANSMISSIVITY*THICKNESS      00018930
C            OF CONFINING BED)                                          00018940
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00018950
C        NONE                                                           00018960
C     METHOD                                                            00018970
C        SUMMATION IS TERMINATED FOR THE INNER SERIES WHEN A TERM       00018980
C        BECOMES LESS THAN 5.E-7/N AND FOR OUTER SERIES WHEN A TERM     00018990
C        BECOMES LESS THAN 5.E-7                                        00019000
C                                                                       00019010
C     ******************************************************************00019020
      REAL*8 DLOG,DABS,S(2),VUM,UU                                      00019030
      REAL*8 TEST,U,UM,EM,EN,SUM1,SUM,SIGN,V,VSQ,VSQU,RMUL,TERM,TERM1   00019040
      TEST=5.D-07                                                       00019050
      VSQ=V*V                                                           00019060
      UU=U                                                              00019070
      DO 6 I=1,2                                                        00019080
C     EVALUATES SERIES FOR LOWER LIMIT = U AND UPPER LIMIT = 1          00019090
      IF (I.EQ.2) U=1.                                                  00019100
      UM=1.                                                             00019110
      EM=-1.                                                            00019120
      SUM1=0.                                                           00019130
      SIGN=-1.                                                          00019140
      VUM=1.                                                            00019150
      VSQU=VSQ/U                                                        00019160
    1 EM=EM+1.                                                          00019170
      IF (EM-.1) 2,3,3                                                  00019180
C     CHECKS FOR M=0                                                    00019190
    2 RMUL=DLOG(U)                                                      00019200
      TERM1=1.                                                          00019210
      GO TO 4                                                           00019220
    3 UM=UM*U                                                           00019230
      IF (VUM.LT.1.D-30) VUM=0.                                         00019240
      VUM=VUM*VSQU                                                      00019250
      RMUL=(UM-VUM)/EM                                                  00019260
      TERM1=TERM1/EM                                                    00019270
    4 SIGN=-SIGN                                                        00019280
      SUM=TERM1                                                         00019290
      TERM=TERM1                                                        00019300
      EN=0.                                                             00019310
    5 EN=EN+1.                                                          00019320
      TERM=TERM*VSQ/(EN*(EN+EM))                                        00019330
      SUM=SUM+TERM                                                      00019340
      IF (TEST.LE.DABS(RMUL*EN*TERM)) GO TO 5                           00019350
C     TRUNCATES INNER SERIES IF OUTER TERM*N*INNER TERM < 5.E-7         00019360
      SUM1=SUM1+SIGN*RMUL*SUM                                           00019370
      IF (EM.LT..1) GO TO 1                                             00019380
      IF (TEST.LE.DABS(RMUL*SUM)) GO TO 1                               00019390
C     TRUNCATES OUTER SERIES IF OUTER TERM*INNER SUM < 5.E-7            00019400
    6 S(I)=SUM1                                                         00019410
      U=UU                                                              00019420
      SERIES=S(2)-S(1)                                                  00019430
      RETURN                                                            00019440
      END                                                               00019450
      REAL FUNCTION FCT*8(X)                                            00019460
C     ******************************************************************00019470
C                                                                       00019480
C     FUNCTION FCT                                                      00019490
C                                                                       00019500
C     PURPOSE                                                           00019510
C        TO COMPUTE FCT(X)=EXP(-Z-V**2/(X+Z))/(X+Z)                     00019520
C     DESCRIPTION OF PARAMETERS                                         00019530
C        X - THE DOUBLE PRECISION VALUE OF X FOR WHICH FCT IS COMPUTED  00019540
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00019550
C        NONE                                                           00019560
C     METHOD                                                            00019570
C        FORTRAN EVALUATION OF FUNCTION                                 00019580
C                                                                       00019590
C     ******************************************************************00019600
      REAL*8 X,V,Z,P,DEXP                                               00019610
      COMMON /C1/ V,Z                                                   00019620
      IF (X) 1,2,2                                                      00019630
    1 FCT=0.                                                            00019640
      GO TO 4                                                           00019650
    2 P=Z+V**2/(X+Z)                                                    00019660
      IF (P-5.D1) 3,3,1                                                 00019670
    3 FCT=DEXP(-P)/(X+Z)                                                00019680
    4 RETURN                                                            00019690
      END                                                               00019700
      SUBROUTINE DQL12(FCT,Y)                                           00019710
C                                                                       00019720
C     ..................................................................00019730
C                                                                       00019740
C        SUBROUTINE DQL12                                               00019750
C                                                                       00019760
C        PURPOSE                                                        00019770
C           TO COMPUTE INTEGRAL(EXP(-X)*FCT(X), SUMMED OVER X           00019780
C                               FROM 0 TO INFINITY).                    00019790
C                                                                       00019800
C        USAGE                                                          00019810
C           CALL DQL12 (FCT,Y)                                          00019820
C           PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT                00019830
C                                                                       00019840
C        DESCRIPTION OF PARAMETERS                                      00019850
C           FCT    - THE NAME OF AN EXTERNAL DOUBLE PRECISION FUNCTION  00019860
C                    SUBPROGRAM USED.                                   00019870
C           Y      - THE RESULTING DOUBLE PRECISION INTEGRAL VALUE.     00019880
C                                                                       00019890
C        REMARKS                                                        00019900
C           NONE                                                        00019910
C                                                                       00019920
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  00019930
C           THE EXTERNAL DOUBLE PRECISION FUNCTION SUBPROGRAM FCT(X)    00019940
C           MUST BE FURNISHED BY THE USER.                              00019950
C                                                                       00019960
C        METHOD                                                         00019970
C           EVALUATION IS DONE BY MEANS OF 12-POINT GAUSSIAN-LAGUERRE   00019980
C           QUADRATURE FORMULA, WHICH INTEGRATES EXACTLY,               00019990
C           WHENEVER FCT(X) IS A POLYNOMIAL UP TO DEGREE 23.            00020000
C           FOR REFERENCE, SEE                                          00020010
C           SHAO/CHEN/FRANK, TABLES OF ZEROS AND GAUSSIAN WEIGHTS OF    00020020
C           CERTAIN ASSOCIATED LAGUERRE POLYNOMIALS AND THE RELATED     00020030
C           GENERALIZED HERMITE POLYNOMIALS, IBM TECHNICAL REPORT       00020040
C           TR00.1100 (MARCH 1964), PP.24-25.                           00020050
C                                                                       00020060
C     ..................................................................00020070
C                                                                       00020080
C                                                                       00020090
C                                                                       00020100
      DOUBLE PRECISION X,Y,FCT                                          00020110
C                                                                       00020120
      X=.3709912104446692 D2                                            00020130
      Y=.814807746742624 D-15*FCT(X)                                    00020140
      X=.2848796725098400 D2                                            00020150
      Y=Y+.3061601635035021 D-11*FCT(X)                                 00020160
      X=.2215109037939701 D2                                            00020170
      Y=Y+.1342391030515004 D-8*FCT(X)                                  00020180
      X=.1711685518746226 D2                                            00020190
      Y=Y+.1668493876540910 D-6*FCT(X)                                  00020200
      X=.1300605499330635 D2                                            00020210
      Y=Y+.836505585681980 D-5*FCT(X)                                   00020220
      X=.962131684245687 D1                                             00020230
      Y=Y+.2032315926629994 D-3*FCT(X)                                  00020240
      X=.6844525453115177 D1                                            00020250
      Y=Y+.2663973541865316 D-2*FCT(X)                                  00020260
      X=.4599227639418348 D1                                            00020270
      Y=Y+.2010238115463410 D-1*FCT(X)                                  00020280
      X=.2833751337743507 D1                                            00020290
      Y=Y+.904492222116809 D-1*FCT(X)                                   00020300
      X=.1512610269776419 D1                                            00020310
      Y=Y+.2440820113198776 D0*FCT(X)                                   00020320
      X=.6117574845151307 D0                                            00020330
      Y=Y+.3777592758731380 D0*FCT(X)                                   00020340
      X=.1157221173580207 D0                                            00020350
      Y=Y+.2647313710554432 D0*FCT(X)                                   00020360
      RETURN                                                            00020370
      END                                                               00020380
      SUBROUTINE BESK(X,N,BK,IER)                                       00020390
C                                                                       00020400
C     ..................................................................00020410
C                                                                       00020420
C        SUBROUTINE BESK                                                00020430
C                                                                       00020440
C           COMPUTE THE K BESSEL FUNCTION FOR A GIVEN ARGUMENT AND ORDER00020450
C                                                                       00020460
C        USAGE                                                          00020470
C           CALL BESK(X,N,BK,IER)                                       00020480
C                                                                       00020490
C        DESCRIPTION OF PARAMETERS                                      00020500
C           X  -THE ARGUMENT OF THE K BESSEL FUNCTION DESIRED           00020510
C           N  -THE ORDER OF THE K BESSEL FUNCTION DESIRED              00020520
C           BK -THE RESULTANT K BESSEL FUNCTION                         00020530
C           IER-RESULTANT ERROR CODE WHERE                              00020540
C              IER=0  NO ERROR                                          00020550
C              IER=1  N IS NEGATIVE                                     00020560
C              IER=2  X IS ZERO OR NEGATIVE                             00020570
C              IER=3  X .GT. 170, MACHINE RANGE EXCEEDED                00020580
C              IER=4  BK .GT. 10**70                                    00020590
C                                                                       00020600
C        REMARKS                                                        00020610
C           N MUST BE GREATER THAN OR EQUAL TO ZERO                     00020620
C                                                                       00020630
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  00020640
C           NONE                                                        00020650
C                                                                       00020660
C        METHOD                                                         00020670
C           COMPUTES ZERO ORDER AND FIRST ORDER BESSEL FUNCTIONS USING  00020680
C           SERIES APPROXIMATIONS AND THEN COMPUTES N TH ORDER FUNCTION 00020690
C           USING RECURRENCE RELATION.                                  00020700
C           RECURRENCE RELATION AND POLYNOMIAL APPROXIMATION TECHNIQUE  00020710
C           AS DESCRIBED BY A.J.M.HITCHCOCK,'POLYNOMIAL APPROXIMATIONS  00020720
C           TO BESSEL FUNCTIONS OF ORDER ZERO AND ONE AND TO RELATED    00020730
C           FUNCTIONS', M.T.A.C., V.11,1957,PP.86-88, AND G.N. WATSON,  00020740
C           'A TREATISE ON THE THEORY OF BESSEL FUNCTIONS', CAMBRIDGE   00020750
C           UNIVERSITY PRESS, 1958, P. 62                               00020760
C                                                                       00020770
C     ..................................................................00020780
C                                                                       00020790
      DIMENSION T(12)                                                   00020800
      BK=.0                                                             00020810
      IF(N)10,11,11                                                     00020820
   10 IER=1                                                             00020830
      RETURN                                                            00020840
   11 IF(X)12,12,20                                                     00020850
   12 IER=2                                                             00020860
      RETURN                                                            00020870
   20 IF(X-170.0)22,22,21                                               00020880
   21 IER=3                                                             00020890
      RETURN                                                            00020900
   22 IER=0                                                             00020910
      IF(X-1.)36,36,25                                                  00020920
   25 A=EXP(-X)                                                         00020930
      B=1./X                                                            00020940
      C=SQRT(B)                                                         00020950
      T(1)=B                                                            00020960
      DO 26 L=2,12                                                      00020970
   26 T(L)=T(L-1)*B                                                     00020980
      IF(N-1)27,29,27                                                   00020990
C                                                                       00021000
C     COMPUTE KO USING POLYNOMIAL APPROXIMATION                         00021010
C                                                                       00021020
   27 G0=A*(1.2533141-.1566642*T(1)+.08811128*T(2)-.09139095*T(3)       00021030
     2+.1344596*T(4)-.2299850*T(5)+.3792410*T(6)-.5247277*T(7)          00021040
     3+.5575368*T(8)-.4262633*T(9)+.2184518*T(10)-.06680977*T(11)       00021050
     4+.009189383*T(12))*C                                              00021060
      IF(N)20,28,29                                                     00021070
   28 BK=G0                                                             00021080
      RETURN                                                            00021090
C                                                                       00021100
C     COMPUTE K1 USING POLYNOMIAL APPROXIMATION                         00021110
C                                                                       00021120
   29 G1=A*(1.2533141+.4699927*T(1)-.1468583*T(2)+.1280427*T(3)         00021130
     2-.1736432*T(4)+.2847618*T(5)-.4594342*T(6)+.6283381*T(7)          00021140
     3-.6632295*T(8)+.5050239*T(9)-.2581304*T(10)+.07880001*T(11)       00021150
     4-.01082418*T(12))*C                                               00021160
      IF(N-1)20,30,31                                                   00021170
   30 BK=G1                                                             00021180
      RETURN                                                            00021190
C                                                                       00021200
C     FROM KO,K1 COMPUTE KN USING RECURRENCE RELATION                   00021210
C                                                                       00021220
   31 DO 35 J=2,N                                                       00021230
      GJ=2.*(FLOAT(J)-1.)*G1/X+G0                                       00021240
      IF(GJ-1.0E70)33,33,32                                             00021250
   32 IER=4                                                             00021260
      GO TO 34                                                          00021270
   33 G0=G1                                                             00021280
   35 G1=GJ                                                             00021290
   34 BK=GJ                                                             00021300
      RETURN                                                            00021310
   36 B=X/2.                                                            00021320
      A=.5772157+ALOG(B)                                                00021330
      C=B*B                                                             00021340
      IF(N-1)37,43,37                                                   00021350
C                                                                       00021360
C     COMPUTE KO USING SERIES EXPANSION                                 00021370
C                                                                       00021380
   37 G0=-A                                                             00021390
      X2J=1.                                                            00021400
      FACT=1.                                                           00021410
      HJ=.0                                                             00021420
      DO 40 J=1,6                                                       00021430
      RJ=1./FLOAT(J)                                                    00021440
      IF(X2J.LT.1.E-40) X2J=0.                                          00021450
C     PREVIOUS STATEMENT ADDED TO IBM SUBROUTINE TO CORRECT UNDERFLOW   00021460
C     PROBLEM ON WATFOR COMPILER                                        00021470
      X2J=X2J*C                                                         00021480
      FACT=FACT*RJ*RJ                                                   00021490
      HJ=HJ+RJ                                                          00021500
   40 G0=G0+X2J*FACT*(HJ-A)                                             00021510
      IF(N)43,42,43                                                     00021520
   42 BK=G0                                                             00021530
      RETURN                                                            00021540
C                                                                       00021550
C     COMPUTE K1 USING SERIES EXPANSION                          ~      00021560
C                                                                       00021570
   43 X2J=B                                                             00021580
      FACT=1.                                                           00021590
      HJ=1.                                                             00021600
      G1=1./X+X2J*(.5+A-HJ)                                             00021610
      DO 50 J=2,8                                                       00021620
      X2J=X2J*C                                                         00021630
      RJ=1./FLOAT(J)                                                    00021640
      FACT=FACT*RJ*RJ                                                   00021650
      HJ=HJ+RJ                                                          00021660
   50 G1=G1+X2J*FACT*(.5+(A-HJ)*FLOAT(J))                               00021670
      IF(N-1)31,52,31                                                   00021680
   52 BK=G1                                                             00021690
      RETURN                                                            00021700
      END                                                               00021710
      SUBROUTINE CONVOL(TIME,A,B,N,IQ,SUM)                              00021720
C********************************************************************** 00021730
C                                                                       00021740
C     PURPOSE                                                           00021750
C        COMPUTES VALUES OF THE CONVOLUTION INTEGRAL FOR LEAKY          00021760
C        AQUIFERS. THE INTEGRAL IS, FROM 0 TO T, OF                     00021770
C        Q(T-T')/T'*EXP(-A/T'-B*T')*DT'.                                00021780
C     DESCRIPTION OF PARAMETERS                                         00021790
C        A,B,SUM ARE REAL; N,IQ ARE INTEGER.                            00021800
C        A - R**2*S/(4*T) - RADIAL DISTANCE SQUARED * STORAGE           00021810
C             COEFFICIENT / 4 * TRANSMISSIVITY.                         00021820
C        B - P'/(S*M') - HYD. COND. OF CONFINING BED DIVIDED BY         00021830
C             AQUIFER STORAGE COEFFICIENT * THICKNESS OF CONF. BED.     00021840
C        N - NUMBER OF INCREMENTS FOR EACH INTERVAL OF THE SUM.         00021850
C        IQ - INDICATES FORM OF DISCHARGE FUNCTION.                     00021860
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     00021870
C        Q                                                              00021880
C     METHOD                                                            00021890
C        APPROXIMATES INTEGRAL BY SUMMING THE TRAPEZOIDAL RULE APPLIED  00021900
C        TO A SEQUENCE OF SEGMENTS. LOWER LIMIT OF FIRST SEGMENT IS     00021910
C        PICKED AT POINT WHERE EXPONENT > -100 .                        00021920
C        IF SUCH A POINT DOES NOT EXIST (A*B > 2500) A FUNCTION VALUE   00021930
C        OF 0 IS RETURNED. UPPER LIMIT = 10 * LOWER LIMIT FOR EACH      00021940
C        SEGMENT. USES INCREMENT OF DELTA T' = (U-L)/N WHERE N IS THE   00021950
C        NUMBER OF INCREMENTS IN THE CALL. CEASES SUMMATION WHEN        00021960
C        EXPONENT < -101 .                                              00021970
C                                                                       00021980
C********************************************************************** 00021990
      REAL*8 DSUM                                                       00022000
      REAL*4 NEWT,NEWTP,NEWX,NEWF                                       00022010
      DSUM=0.D+0                                                        00022020
      IS=0                                                              00022030
C     INITIAL T' COMPUTED FROM A,B                                      00022040
      AB=A*B                                                            00022050
      IF (AB.GE.2500.) GO TO 7                                          00022060
      IF (B.GT.0.) GO TO 2                                              00022070
    1 OLDT=.01*A                                                        00022080
      GO TO 3                                                           00022090
    2 OLDT=(1.-SQRT(1.-AB/2500.))*50./B                                 00022100
      IF (OLDT.EQ.0.) GO TO 1                                           00022110
C     INITIAL T-T'                                                      00022120
    3 OLDTP=TIME-OLDT                                                   00022130
      OLDX=-A/OLDT-B*OLDT                                               00022140
      OLDF=Q(OLDTP,IQ)*EXP(OLDX)/OLDT                                   00022150
C     END OF SUMMATION SEGMENT IS 10 TIMES THE BEGINNING                00022160
    4 ENDT=10.*OLDT                                                     00022170
      IF (ENDT.LT.TIME) GO TO 5                                         00022180
      IF (OLDT.GE.TIME) GO TO 7                                         00022190
      IS=1                                                              00022200
      ENDT=TIME                                                         00022210
C     DELTA T' IS COMPUTED FROM LENGTH AND NUMBER OF INCREMENTS         00022220
    5 DELT=(ENDT-OLDT)/N                                                00022230
      DO 6 I=1,N                                                        00022240
C     T' IS INCREMENTED BY DELTA T'                                     00022250
      NEWT=OLDT+DELT                                                    00022260
      NEWX=-A/NEWT-B*NEWT                                               00022270
C     TERMINATES SUMMATION WHEN EXP(-A/T'-B*T') < 1.37E-44              00022280
      IF (NEWX.LT.-101.) GO TO 7                                        00022290
      NEWTP=TIME-NEWT                                                   00022300
      NEWF=Q(NEWTP,IQ)*EXP(NEWX)/NEWT                                   00022310
      DSUM=DSUM+(NEWF+OLDF)*DELT                                        00022320
      OLDT=NEWT                                                         00022330
      OLDF=NEWF                                                         00022340
    6 CONTINUE                                                          00022350
      IF (IS.GT.0) GO TO 7                                              00022360
C     IF T' < T, BEGINS A NEW SEGMENT                                   00022370
      GO TO 4                                                           00022380
    7 SUM=DSUM/2.D+0                                                    00022390
      RETURN                                                            00022400
      END                                                               00022410
      FUNCTION Q(TIME,IQ)                                               00022420
C********************************************************************** 00022430
C                                                                       00022440
C     PURPOSE                                                           00022450
C        COMPUTES THE DISCHARGE FUNCTION, Q(T)                          00022460
C     DESCRIPTION OF PARAMETERS                                         00022470
C        TIME - REAL - ELAPSED TIME SINCE BEGINNING OF DISCHARGE.       00022480
C        IQ - INTEGER - INDICATES FORM OF DISCHARGE FUNCTION.           00022490
C             IQ=1,2,3. CASES A,B,C, RESPECTIVELY, OF HANTUSH,M.S.,     00022500
C                  1964, HYDRAULICS OF WELLS IN CHOW, VEN TE, ED.,      00022510
C                  ADVANCES IN HYDROSCIENCE, VOL. 1: ACADEMIC PRESS,    00022520
C                  NEW YORK, P. 343,344.                                00022530
C             IQ=4. DISCHARGE IS A FIFTH DEGREE POLYNOMIAL OF TIME.     00022540
C             IQ=5. DISCHARGE IS A PIECEWISE LINEAR FUNCTION OF UP TO   00022550
C                  8 SEGMENTS.                                          00022560
C     METHOD                                                            00022570
C        FORTRAN EVALUATION OF FUNCTIONS.                               00022580
C                                                                       00022590
C********************************************************************** 00022600
      COMMON AQ(6),TI(9),AI(9),BI(9),QST,DELTA,TSTAR                    00022610
      GO TO (1,2,3,4,5), IQ                                             00022620
    1 Q=QST*(1.+DELTA*EXP(-TIME/TSTAR))                                 00022630
      RETURN                                                            00022640
    2 Q=QST*(1.+DELTA/(1.+TIME/TSTAR))                                  00022650
      RETURN                                                            00022660
    3 Q=QST*(1.+DELTA/SQRT(1.+TIME/TSTAR))                              00022670
      RETURN                                                            00022680
    4 Q=AQ(1)+TIME*(AQ(2)+TIME*(AQ(3)+TIME*(AQ(4)+TIME*(AQ(5)+TIME*AQ(6)00022690
     1))))                                                              00022700
      RETURN                                                            00022710
    5 DO 6 I=2,9                                                        00022720
      IF (TIME.LE.TI(I)) GO TO 7                                        00022730
    6 CONTINUE                                                          00022740
      I=9                                                               00022750
    7 Q=AI(I)+BI(I)*(TIME-TI(I-1))                                      00022760
      RETURN                                                            00022770
      END                                                               00022780
