KGS Home Hydrogrology Home

Kansas Geological Survey, Ground Water Series 3, originally published in 1980
Prev Page--Start


Appendix I--Listing of Program THEIS

C     PROGRAM: THEIS
C     PURPOSE: CALCULATES THE DRAWDOWN FOR A WELL PUMPING AT
C         CONSTANT DISCHARGE FROM AN INFINITE AQUIFER BY THE
C         THEIS EQUATION
C
C     ARGUMENTS (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME )
C          RS CONSTANT OBSERVATION SWITCH; YES OR NO
C          TS CONSTANT OBSERVATION TIME SWITCH; YES OR NO
C          US UNIT SWITCH, USE GAL, DAYS, FT; YES OR NO
C       PRTSW PRINT SWITCH FOR U, W, AND SENSITIVITY; YES OR NO
C          SC STORAGE COEFFICIENT FOR AQUIFER (UNITLESS)
C          KB TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OR (L**2/T)
C           Q CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C           R OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C           T OBSERVATION TIME (DAYS) OR (T)
C           S DRAWDOWN (FT) OR (L)
C        DSDT SENSITIVITY W.R.T. TRANSMISSIVITY (FT**2*DAY/GAL) OR (T/L)
C       DSDSC SENSITIVITY W.R.T. STORAGE (FT) OR (L)
C -------------------------------------------------------------
C
C READ IN ARGUMENTS
      CHARACTER RS, TS, US, PRTSW 
      REAL KB 
      ICOUNT =0 
      PRINT 1 
      PRINT 15
   15 FORMAT(' DO YOU WANT U, W, AND SENSITIVITY VALUES PRINTED ? ')
      READ:PRTSW
    1 FORMAT(1H0)
      PRINT 14
   14 FORMAT(' DO YOU WANT TO USE THE GAL, DAY, FT SYSTEM ? ')
      READ:US
      PRINT 2
    2 FORMAT(' STORAGE COEFFICIENT (UNITLESS)')
      READ:SC
      PRINT 3
    3 FORMAT(' TRANSMISSIVITY (GAL/DAY/FT) OR (L**2/T)')
      READ:KB
      PRINT 4
    4 FORMAT(' CONSTANT PUMPAGE (GAL/DAY) OR (L**3/T)')
      READ:Q
      PRINT 10
   10 FORMAT (' CONSTANT OBSERVATION POINT?
      READ:RS
      PRINT 11
   11 FORMAT (' CONSTANT OBSERVATION TIME? READ:TS
   80 IF (ICOUNT .GT. 0 AND. RS EQ. 3HYES) GO TO 70
      PRINT 5
    5 FORMAT(' OBSERVATION POINT (FT) OR (L)') 
      READ:R
   70 IF (ICOUNT .GT. 0 AND. TS EQ. 3HYES) GO TO 95
      PRINT 6
    6 FORMAT(' OBSERVATION TIME (DAYS) OR (T)')
      READ:T
C CALCULATE THEIS SOLUTION
C THEIS EQUATION FROM 'INTRO TO HYDROLOGY' BY VIESSMAN PAGE 257
C RATIONAL APPROXIMATION FROM 'HANDBOOK OF MATH FUNCTIONS'
C BY ABRAMOWITZ AND STEGUN PAGE 231
   95 CONTINUE
      U =  R*R*SC/(4.0*KB*T)
      IF (US .EQ. 3HYES ) U = 7.48*U
      IF (V .LE. 0.0) GO TO 96
      IF (U .GT. 1.0) GO TO 86
      W= -ALOG(U)-.57721566+.99999193*U-.24991055*U*U+
     &     .05519968*U*U*U-.00976004*U**4+.00107857*U**5
      GO TO 87
   86 W (EXP(-U)/U)*(U*U+2.334733*U+.250621)/
     &     (U*U+3.330657*U + 1.681534)
   87 S = (Q/(4.0*3.14159*KB))*W
      DSDT = (Q/(4.0*3.14159*KB**2))*(-W + EXP(-U))
      DSDSC=-(Q/(4.0*3.14159*KB*SC))*EXP(-U)
C OUTPUT RESULTS
   90 PRINT 7, S
    7 FORMAT(1H0, 'THE DRAWDOWN IS ', F20.8)
      IF (PRTSW EQ. 3HYES ) PRINT 8,U
    8 FORMAT (' U = ', E20.8)
      IF (PRTSW EQ. 3HYES ) PRINT 9,W
    9 FORMAT (' W = ', E20.8 )
      IF (PRTSW EQ. 3HYES ) PRINT 12,DSDT
   12 FORMAT ( ' SENSITIVITY WITH RESPECT TO TRANSMISSIVITY = ', E20.8 )
      IF (PRTSW EQ. 3HYES ) PRINT 21,DSDSC
   21 FORMAT (' SENSITIVITY WITH RESPECT TO STORAGE = ', E20.8 )
      GO TO 97
   96 PRINT 13
   13 FORMAT (' DRAWDOWN UNDEFINED U TOO SMALL ')
   97 PRINT 1
      ICOUNT ICOUNT +1
      IF (RS EQ. 3HYES AND. TS EQ. 3HYES ) GO TO 99
      GO TO 80
   99 STOP
      END

Appendix II--Listing of subroutine THEIS

      SUBROUTINE THEIS(SC)KB,Q,R,T,S,DSDT,DSDSC,UNIT)
C     PURPOSE   CALCULATES THE DRAWDOWN FOR A WELL PUMPING AT
C         CONSTANT DISCHARGE FROM AN INFINITE AQUIFER BY THE
C         THEIS EQUATION
C
C     ARGUMENTS (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C        SC    STORAGE COEFFICIENT FOR AQUIFER (UNITLESS)
C        KB    TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OF (L**2/T)
C         Q    CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C         R    OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C         T    OBSERVATION TIME (DAYS) OR (T)
C         S    DRAWDOWN (FT) OR (L)
C      DSDT    SENSITIVITY W.R.T. TRANSMISSIVITY (FT**2*DAY/GAL) OR (T/L)
C     DSDSC    SENSITIVITY W.R.T. STORAGE (FT) OR (L)
C      UNIT    SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C                 CONSISTENT UNITS ARE ASSUMED
C -------------------------------------------------------------
C
C CALCULATE THEIS SOLUTION
C THEIS EQUATION FROM 'INTRO TO HYDROLOGY' BY VIESSMAN PAGE 257
C RATIONAL APPROXIMATION FROM 'HANDBOOK OF MATH FUNCTIONS'
C   BY ABRAMOWITZ AND STEGUN PAGE 231
      REAL KB
      CHARACTER UNIT
      U = R*R*SC/(4.0*KB*T)
      IF (UNIT EQ. 3HYES U = 7.48*U
      IF (U .LE. 0.0) GO TO 96
      IF (U .GT. 1.0) GO TO 86
      W= -ALOG(U)-.57721566+.99999193*U-.24991055*U*U+
     &      .05519968*U*U*U-.00976004*U**4+.00107857*U**5
      GO TO 87
   86 W = (EXP(-U)/U)*(U*U+2.334733*U+.250621)/
     &      (U*U+3.330657*U + 1.681534)
   87 S = (Q/(4.0*3.14159*KB))*W
      DSDT = (Q/(4.0*3.14159*KB**2))*(-W + EXP(-U))
      DSDSC=-(Q/(4.0*3.14159*KB*SC))*EXP(-U)
      RETURN
   96 PRINT 13
   13 FORMAT DRAWDOWN UNDEFINED U TOO SMALL
   99 STOP
     END

Appendix III--Listing of program THEISFIT (time sharing version)

C     PROGRAM THEISFIT (TIME SHARING VERSION)
C     PURPOSE: TO CALCULATE THE BEST FIT STORAGE AND TRANSMISSIVITY
C         BY FITTING THE THEIS EQUATION TO EXPERIMENTAL PUMPTEST DATA IN
C         A LEAST SQUARES SENSE.
C
C     IMPORTANT VARIABLES: (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C         UNIT    SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C                    CONSISTENT UNITS ARE ASSUMED.
C        GUESS    DO YOU WANT TO ENTER AN INITIAL GUESS FOR SC AND KB.?
C                 YES OR NO IS THE APPROPRIATE RESPONSE.
C           SC    STORAGE COEFFICIENT OF AQUIFER (UNITLESS)
C           KB    TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OR (L"*2/T)
C            Q    CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L*113/T)
C            R    OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C            N    NUMBER OF DRAWDOWN-TIME PAIRS TO BE READ
C         S(I)    EXPERIMENTAL DRAWDOWN AT TIME I (FT) OR (L)
C         T(I)    THE ITH TIME AT WHICH AN EXPERIMENTAL MEASUREMENT
C                    FOR THE DRAWDOWN IS MADE (DAYS) OR (T)
C        SIGMA    THE RMS ERROR IN DRAWDOWN AFTER THE BEST
C                    FIT HAS BEEN OBTAINED
C ---------------------------------------------------------------------------
      DIMENSION S(100),T(100),SP(100),DSDT(100),DSDSC(100)
      REAL KB
      CHARACTER GUESS, UNIT
C READ IN THE INITIAL DATA
      PRINT 801
  801 FORMAT (' DO YOU WANT TO USE GAL, DAY, FT SYSTEM ? ')
      READ:UNIT
      PRINT 105
  105 FORMAT (' GUESS FOR STORAGE AND TRANSMISSIVITY ')
      READ: GUESS
      IF (GUESS NE. 3HYES ) GO TO 114
      PRINT 112
  112 FORMAT (' ESTIMATE FOR STORAGE ? ')
      READ: SC
      PRINT 113
  113 FORMAT(' ESTIMATE FOR TRANSMISSIVITY ? GAL/DAY/FT OR L**2/T ')
      READ: KB
  114 PRINT 115
  115 FORMAT (' CONSTANT PUMPAGE RATE ? GAL/DAY OR L**3/T ')
      READ: Q
      PRINT 116
  116 FORMAT (' OBSERVATION DISTANCE FROM PUMPING WELL ? FT OR L ')
      READ: R
      PRINT 117
  117 FORMAT(' NUMBER OF DRAWDOWN-TIME PAIRS TO BE READ ? ')
      READ: N
C ECHO PRINT THE INITIAL DATA
      WRITE (6, 106 ) SC, KB, Q, R, N
  106 FORMAT ( ' ECHO THE INITIAL DATA ' / 5H SC =, E15.8, 5H KB =, E15.8
     &       /4H Q =, E15.8, 4H R =, E15.8, 4H N =, I5 )
C TYPE IN THE DRAWDOWN-TIME PAIRS IN ORDER OF INCREASING TIME
      PRINT 118
  118 FORMAT (' TYPE IN DRAWDOWN-TIME PAIRS IN ORDER OF INCREASING TIME. ')
      DO 120 1 = 1,N
  120 READ: S(I),T(I)
C ECHO PRINT THE DRAWDOWN-TIME PAIRS
      WRITE (6, 107) (S(I),T(I), I = 1, N)
  107 FORMAT (' THE PUMP TEST DATA IN DRAWDOWN-TIME PAIRS '/(2E20.8))
  101 FORMAT ( 8F10.0 )
  100 FORMAT ( 4F10.0, I10 )
C CALCULATE THE INITIAL GUESS FOR KB AND SC IF NOT GIVEN
      NUM = 4
      IF (GUESS EQ. 6HYES GO TO 1
      SUMLNT = 0.0
      SMLNT2 = 0.0
      SUMSE = 0.0
      SSELNT = 0.0
      DO 50 1 = 1,NUM
      II = N +1-I
      ALOGT = ALOG(T(II))
      SUMLNT = SUMLNT + ALOGT
      SMLNT2 = SMLNT2 + ALOGT**2
      SUMSE = SUMSE + S(II)
   50 SSELNT = SSELNT + S(II)*ALOGT
      SLOPE = (NUM*IRSSELNT-SUMSE*SUMLNT)/(NUM*SMLNT2-SUMLNT**2)
      CON =-(SLOPE*SUMLNT - SUMSE)/NUM
      KB = Q/(SLOPE*4.0*3.14159)
      ALNTO = -CON/SLOPE
      SC = 4.0*KB/( R**C*2*EXP(.5772-ALNTO))
      IF (UNIT  EQ. 3HYES ) SC = SC/7.48
C PRINT THE CALCULATED GUESS FOR KB AND SC
      WRITE (6,110) KB, SC
  110 FORMAT (' THE CALCULATED GUESS FOR KB AND SC IS ', 2E16.8 )
    1 ICOUNT = 0
    5 TEMPKB = KB
      TEMPSC = SC
C USE SENSITIVITY ANALYSIS AND LEAST SQUARES TO FIND A BETTER KB AND SC
      DO 10 1 = 1, N
   10 CALL THEIS(SC,KB,Q,R,T(I),SP(I),DSDT(I),DSDSC(I),UNIT)
      SSUS = 0.0
      SSUT = 0.0
      SUTUS = 0.0
      SUSDIF = 0.0
      SUTDIF = 0.0
      DO 15 1 = 1, N
      SSUS = DSDSC(I)11*2+SSUS
      SSUT = DSDT(I)**2+SSUT
      SUTUS = DSDSC(I)*DSDT(I) + SUTUS
      SUSDIF = DSDSC(I)*(S(I)-SP(I)) + SUSDIF
   15 SUTDIF = DSDT(I)*(S(I)-SP(I)) + SUTDIF
      DELTSC = (SSUT*SUSDIF-SUTUS*SUTDIF)/(SSUS*SSUT-SUTUS**2)
      IF (DELTSC .LT. -.95*SC ) DELTSC = -.95*SC
      IF (DELTSC .GT. 2.0*SC )  DELTSC = 2.0*SC
      SC = SC + DELTSC
      IF (SC .GT. 1.0 ) SC = 1.0
      DELTKB = (SUTDIF-DELTSC*SUTUS)/SSIJT
      IF (DELTKB .LT. -.95*KB ) DELTKB = -.95*KB
      IF (DELTKB .GT. 2.0*KB ) DELTKB = 2.0*KB
      KB = KB + DELTKB
      WRITE (6,104) KB, SC
  104 FORMAT( BEST FIT KB AND SC THIS ITERATION ARE ', 2E16.8)
      ICOUNT ICOUNT + I
C IF IT DID NOT CONVERGE PRINT MESSAGE
      IF ( ICOUNT .GT. 20) GO TO 40
C DO ANOTHER ITERATION IF IT HASNT CONVERGED AND ITERATION LIMIT NOT
C EXCEEDED.
      IF (ABS((TEMPKB-KB)/KB) .GT. .001 .OR. ABS((TEMPSC-SC)/SC) .GT. .001)
     &   GO TO 5
C PRINT BEST FIT DRAWDOWN-TIME PAIRS IF CONVERGENCE OBTAINED
  108 FORMAT (' THE BEST FIT DRAWDOWN-TIME PAIRS     ')
      WRITE (6,108)
      DO 20 1 = 1, N
      CALL THEIS(SC,KB,Q,R,T(I),SP(I),DSDT(I),DSDSC(I),UNIT)
      WRITE (6,109) SP(I), T(I)
  109 FORMAT(2E20.8)
   20 CONTINUE
C CALCULATE THE RMS ERROR IN DRAWDOWN AFTER BEST FIT
      SUM = 0.0
      DO 30 1 = 1, N
   30 SUM = SUM + (SP(I)-S(I))**2
      SIGMA = SQRT(SUM/N)
      WRITE(6,111) SIGMA
  111 FORMAT(' THE RMS ERROR FOR DRAWDOWN IS ', E20.8)
      STOP
   40 WRITE (6,125)
  125 FORMAT (' DID NOT CONVERGE IN 20 ITERATIONS)
      STOP
      END
      SUBROUTINE THEIS(SC,KB,Q,R,T,S,DSDT,DSDSC$UNIT)
C     PURPOSE   CALCULATES THE DRAWDOWN FOR A WELL PUMPING AT
C         CONSTANT DISCHARGE FROM AN INFINITE AQUIFER BY THE
C         THEIS EQUATION
C
C      ARGUMENTS (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C         SC    STORAGE COEFFICIENT FOR AQUIFER (UNITLESS)
C         KB    TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OF (L**2/T)
C          Q    CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C          R    OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C          T    OBSERVATION TIME (DAYS) OR (T)
C          S    DRAWDOWN (FT) OR (L)
C       DSDT    SENSITIVITY W.R.T. TRANSMISSIVITY (FT**2*DAY/GAL) OR (T/L)
C      DSDSC    SENSITIVITY W.R.T. STORAGE (FT) OR (L)
C       UNIT    SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C                  CONSISTENT UNITS ARE ASSUMED
C-------------------------------------------------------------
C
C CALCULATE THEIS SOLUTION
C THEIS EQUATION FROM 'INTRO TO HYDROLOGY' BY VIESSMAN PAGE 257
C RATIONAL APPROXIMATION FROM 'HANDBOOK OF MATH FUNCTIONS'
C   BY ABRAMOWITZ AND STEGUN PAGE 231
      REAL KB
      CHARACTER UNIT
      U = R*R*SC/(4.O*KB*T)
      IF (UNIT EQ. 3HYES) U = 7.48*U
      IF (U .LE. 0.0) GO TO 96
      IF (U .GT. 1.0) GO TO 86
      W= -ALOG(U)-.57721566+.99999193*U-.24991055*U*U+
     &   .05519968*U*U*U-.00976004*U**4+.00107857*U**5
      GO TO 87
   86 W = (EXP(-U)/U)*(U*U+2.334733*U+.250621)/
     &   (U*U+3.330657*U + 1.681534)
   87 S = (Q/(4.0*3.14159*KB))*W
      DSDT = (Q/(4.0*3.14159*KB**2))*(-W + EXP(-U))
      DSDSC=-(Q/(4.0*3.14159"KB*SC))*EXP(-U)
      RETURN
   96 PRINT 13
   13 FORMAT DRAWDOWN UNDEFINED U TOO SMALL
   99 STOP
      END

Appendix IV--Listing of program THEISFIT (batch version)

C     PROGRAM THEISFIT (FORTRAN IV BATCH VERSION)
C     PURPOSE: TO CALCULATE THE BEST FIT STORAGE AND TRANSMISSIVITY
C         BY FITTING THE THEIS EQUATION TO EXPERIMENTAL PUMPTEST DATA
C         IN A LEAST SQUARES SENSE.
C
C     IMPORTANT VARIABLES: (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C         UNIT    SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C                    CONSISTENT UNITS ARE ASSUMED.
C        GUESS    DO YOU WANT TO ENTER AN INITIAL GUESS FOR SC AND KB?
C                    YES OR NO IS THE APPROPRIATE RESPONSE.
C           SC    STORAGE COEFFICIENT OF AQUIFER (UNITLESS)
C           KB    TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OR (L**2/T)
C            Q    CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C            R    OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C            N    NUMBER OF DRAWDOWN-TIME PAIRS TO BE READ
C         S(I)    EXPERIMENTAL DRAWDOWN AT TIME I (FT) OR (L)
C         T(I)    THE ITH TIME AT WHICH AN EXPERIMENTAL MEASUREMENT
C                    FOR THE DRAWDOWN IS MADE (DAYS) OR (T)
C        SIGMA    THE RMS ERROR IN DRAWDOWN AFTER THE BEST
C                    FIT HAS BEEN OBTAINED
C ---------------------------------------------------------------------------
      DIMENSION S(100), T(100),SP(100),DSDT(100),DSDSC(100)
      REAL KB
      CHARACTER GUESS, UNIT
C READ IN INITIAL DATA
      READ (5,105) GUESS, UNIT
  105 FORMAT ( 2A6 )
      READ (5,100) SC, KB, Q, R, N
C ECHO PRINT THE INITIAL DATA
      WRITE (6, 106 ) SC, KB, Q, R, N
  106 FORMAT ( ' ECHO THE INITIAL DATA ' / 5H SC =, E15.8, 5H KB =, E15.8
     + /4H Q =, E15.8, 4H R =, E15.8, 4H N =, I5 )
C READ IN DRAWDOWN-TIME PAIRS IN ORDER OF INCREASING TIME
      READ (5, 101) (S(I),T(I), I = 1, N)
C ECHO PRINT THE DRAWDOWN-TIME PAIRS
      WRITE (6, 107) (S(I),T(I), I = 1, N)
  107 FORMAT THE PUMP TEST DATA IN DRAWDOWN-TIME PAIRS   '/(2E20.8))
  101 FORMAT 8F10.0 )
  100 FORMAT 4F10.0, I10 )
C CALCULATE THE INITIAL GUESS FOR KB AND SC IF NOT GIVEN
      NUM = 4
      IF (GUESS EQ. 6HYES   ) GO TO 1
      SUMLNT = 0.0
      SMLNT2 = 0.0
      SUMSE = 0.0
      SSELNT = 0.0
      DO 50 1 = 1,NUM
      II = N + 1 -1
      ALOGT = ALOG(T(II))
      SUMLNT = SUMLNT + ALOGT
      SMLNT2 = SMLNT2 + ALOGT**2
      SUMSE = SUMSE + S(II)
   50 SSELNT = SSELNT + S(II)*ALOGT
      SLOPE = (NUM*SSELNT-SUMSE*SUMLNT)/(NUM*SMLNT2-SUMLNT**2)
      CON = (SLOPE*SUMLNT - SUMSE)/NUM
      KB = Q/(SLOPE*4.0*3.14159)
      ALNTO = -CON/SLOPE
      SC = 4.0*KB/( R**2*EXP(.5772-ALNTO))
      IF (UNIT  EQ. 3HYES ) SC = SC/7.48
C PRINT THE CALCULATED GUESS FOR KB AND SC
      WRITE (6,110) KB, SC
  110 FORMAT (' THE CALCULATED GUESS FOR KB AND SC IS ', 2E16.8)
    1 ICOUNT = 0
    5 TEMPKB = KB
      TEMPSC = SC
C USE SENSITIVITY ANALYSIS AND LEAST SQUATES TO FIND A BETTER KB AND SC DO 10 1 = 1, N
   10 CALL THEIS(SC)KB,Q,R,T(I),SP(I),DSDT(I),DSDSC(I),UNIT)
      SSUS = 0.0
      SSUT = 0.0
      SUTUS = 0.0
      SUSDIF = 0.0
      SUTDIF = 0.0
      DO 15 1 = 1, N
      SSUS = DSDSC(I)**2+SSUS
      SSUT = DSDT(I)**2+SSUT
      SUTUS = DSDSC(I)*DSDT(I) + SUTUS
      SUSDIF = DSDSC(I)*(S(I)-SP(I)) + SUSDIF
   15 SUTDIF = DSDT(I)*(S(I)-SP(I)) + SUTDIF
      DELTSC = (SSUT*SUSDIF-SUTUS*SUTDIF)/(SSUS*SSUT-SUTUS**2)
      IF (DELTSC .LT. -.95*SC ) DELTSC = -.95*SC
      IF (DELTSC .GT. 2.0*SC ) DELTSC = 2.0*SC
      SC = SC + DELTSC
      IF (SC .GT. 1.0 ) SC = 1.0
      DELTKB = (SUTDIF-DELTSC*SUTUS)/SSUT
      IF (DELTKB .LT. -.95*KB) DELTKB = -.95*KB
      IF (DELTKB .GT. 2.0*KB DELTKB = 2.0*KB
      KB = KB + DELTKB
      WRITE (6,104) KB, SC
  104 FORMAT( ' BEST FIT KB AND SC THIS ITERATION ARE ', 2E16.8)
      ICOUNT ICOUNT + 1
C IF IT DID NOT CONVERGE PRINT MESSAGE
      IF ( ICOUNT .GT. 20) GO TO 40
C DO ANOTHER ITERATION IF IT HASNT CONVERGED AND ITERATION LIMIT NOT
C EXCEEDED.
      IF(ABS((TEMPKB-KB)/KB) .GT. .001 .OR. ABS((TEMPSC-SC)/SC) .GT. .001)
     & GO TO 5
C PRINT BEST FIT DRAWDOWN-TIME PAIRS IF CONVERGENCE OBTAINED
  108 FORMAT (' THE BEST FIT DRAWDOWN-TIME PAIRS     ')
      WRITE (6,108)
      DO 20 1 = 1, N
      CALL THEIS(SC,KB,Q,R,T(I),SP(I),DSDT(I),DSDSC(I),UNIT)
      WRITE (6,109) SP(I), T(I)
  109 FORMAT(2E20.8)
   20 CONTINUE
C CALCULATE THE RMS ERROR IN DRAWDOWN AFTER BEST FIT
      SUM = 0. 0
      DO 30 1 = 1, N
   30 SUM = SUM + (SP(I)-S(I))**2
      SIGMA = SQRT(SUM/N)
      WRITE(6,111) SIGMA
  111 FORMAT(' THE RMS ERROR FOR DRAWDOWN IS ', E20.8)
      STOP
   40 WRITE (6,112)
  112 FORMAT(' DID NOT CONVERGE IN 20 ITERATIONS ')
      STOP
      END

      SUBROUTINE THEIS(SC,KB,Q,R,T,S,DSDT,DSDSC,UNIT)
C     PURPOSE CALCULATES THE DRAWDOWN FOR A WELL PUMPING AT
C         CONSTANT DISCHARGE FROM AN -INFINITE AQUIFER BY THE
C         THEIS EQUATION
C
C     ARGUMENTS (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C        SC    STORAGE COEFFICIENT FOR AQUIFER (UNITLESS)
C        KB    TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OR (L**2/T)
C         Q    CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C         R    OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C         T    OBSERVATION TIME (DAYS) OR (T)
C         S    DRAWDOWN (FT) OR (L)
C      DSDT    SENSITIVITY W.R.T. TRANSMISSIVITY (FT**2*DAY/GAL) OR (T/L)
C     DSDSC    SENSITIVITY W.R.T. STORAGE (FT) OR (L)
C      UNIT    SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C                CONSISTENT UNITS ARE ASSUMED
C -------------------------------------------------------------
C
C CALCULATE THEIS SOLUTION
C THEIS EQUATION FROM 'INTRO TO HYDROLOGY' BY VIESSMAN PAGE 257
C RATIONAL APPROXIMATION FROM 'HANDBOOK OF MATH FUNCTIONS'
C   BY ABRAMOWITZ AND STEGUN PAGE 231
      REAL KB
      CHARACTER UNIT
      U = R*R*SC/(4.0*KB*T)
      IF (UNIT EQ. 6HYES   ) U = 7.48*U
      IF (U .LE. 0.0) GO TO 96
      IF (U .GT. 1.0) GO TO 86
      W= -ALOG(U)-.57721566+.99999193*U-.24991055*U*U+
     &    1.05519968*U*U*U-.00976004*U**4+.00107857*U**5
      GO TO 87
   86 W = (EXP(-U)/U)*(U*U+2.334733*U+.250621)/
     &      (U*U+3.330657*U + 1.681534)
   87 S = (Q/(4.0*3.14159*KB))*W
      DSDT = (Q/(4.0*3.14159*KB**2))*(-W + EXP(-U))
      DSDSC = -(Q/(4.0*3.14159*KB*SC))*EXP(-U)
      RETURN
   96 PRINT 13
   13 FORMAT ( ' DRAWDOWN UNDEFINED U TOO SMALL ')
   99 STOP
      END

Prev Page--Start

Kansas Geological Survey, The Theis Equation
Placed on web June 18, 2012; originally published in 1980.
Comments to webadmin@kgs.ku.edu
The URL for this page is http://www.kgs.ku.edu/Publications/Bulletins/GW3/02_append.html