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) H0 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=1, USES A GAUSSIAN-LAGUERRE QUADRATURE FORMULA TO 00018410 C EVALUATE INTEGRAL(F) FROM U TO INF. 00018420 C (2) V**20. 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