Download Program Thies - Biosystems Engineering for Fuel and Chemicals | BSysE 595 and more Study notes Engineering in PDF only on Docsity! Program Theis c This code is modified from Dr. J.P. Fairley's (UI) original code by S. Dun. c One main change is in the I/O format. In addition, the modified code c provides a more explict way of handling varing pumping rates. Also, c one can choose to use latitude/longitude coordinate system. c Note: the units in the calculation and output are consistent c with units of the input. c c C---+----1----+----2----+----3----+----4----+----5----+----6 C The input file has the following format: C C NW,NB,NT,DT,SY,ERR c TRANS,STORE c XB(1),XB(1) c ........... c XB(I),YB(I) c XB(NB),YB(NB) c C XL(I),YL(I),RAD(I),NQT(I) (I loop NW times) c ST(I,J),QL(I,J) (J loop NQT(I) times) C C---+----1----+----2----+----3----+----4----+----5----+----6 c c Modified by S. Dun, June 15, 2003. c c C Program file "theis.f". This file reads data from C an input file called "theis.dat", and uses the data C and the Theis solution to caluculate drawdowns in C a hypothetical well field described by the input C data. The program was written by JPF on 11/11/02, C and last edited by JPF on 11/13/02. C---+----1----+----2----+----3----+----4----+----5----+----6 C The input file has the following format: C C NX NY NT NW (4I10) C DX DY DT RE (4E10) C TRANS STORE STATIC RAD2 (4E10) C QL TL XL YL (4E10) C QL TL XL YL C ....Etc., until QL(NW) C QL TL XL YL C C The input parameters are defined in the data dictionary below. C As currently set up, the program can handle up to 58 timesteps. C The number of timesteps can be increased up to 1000 by changing C the FORMAT statement 1100. For more timesteps, increase C both the FORMAT statement and the dimensioned size of array C "DRAWD", which will currently hold a maximum of 1000 values. C C---+----1----+----2----+----3----+----4----+----5----+----6 C C **********DATA DICTIONARY******************* C C Reals C C DX spacing between points in the X direction C DY spacing between points in the Y direction C DT timestep C ERR convergence criterion (absolute) C TRANS aquifer transmissivity C STORE aquifer storativity C STATIC aquifer static (initial) water level C QL well #? pumping rate (withdrawl positive) C TL time well #? begins pumping C XL X location of well #? C YL Y location of well #? C RAD2 effective well radius C XLOC, YLOC, TLOC x,y,t at current calculational point C U Theis' u, equals U=U1/U2 C RDIST square of radial distance from the well C WF Theis' well function (exponential integral) C REM remainder (for error checking) C WFNEW intermediate infinite series W(u) values C POWER term number for infinite series C FACT factorial of POWER C NEXT dummy number for calculating FACTorial C ASY dummy variable for calculating asymptotic series C DRAWD amount of drawdown; also, pumping surface elevation C C Integers C C NX number of points in the X (easting) direction C NY number of points in the Y (northing) direction C NT number of timesteps C NW number of wells C I,J,K,L, counting integers (indices) C C Characters C C ********************************************** c include 'location.inc' C C---+----1----+----2----+----3----+----4----+----5----+----6 C C DECLARE VARIABLES AND INITIALIZE ARRAYS c IMPLICIT NONE REAL DT,ERR,TRANS,STORE,TLOC,U1,U2,U REAL RDIST,WF,REM,WFNEW,POWER,FACT,NEXT REAL XL(20),YL(20),RAD(20),QL(20,200),ST(20,200) REAL XB(20),YB(20) REAL DRAWD(20),ASY,SY INTEGER NT,NW,NB,I,J,K,L,NQT(20),iflag C C---+----1----+----2----+----3----+----4----+----5----+----6 c c Choose the coordinate system you will use c (1) Cartisian X,Y c (2) X,Y in longitude/latitude c write (6,1350) c read (5,1500,err=10) iflag c 10 if (iflag.lt.1.or.iflag.gt.2) then write (6,1650) iflag = 1 end if 1350 format (/,' Confirm the coordinate system used in your input:',/,