PROGRAM THEIS
parameter(mmax=300,nmax=300,isize=mmax*nmax)
parameter(mx=mmax+1,ny=nmax+1)
dimension v(0:mx,0:ny)
integer iu(mmax,nmax)
character*1 cu(mmax,nmax)
character fn*10
C
C THIS PROGRAM SOLVES THE THEIS SOLUTION
C
COMMON /BLOCK1/ XBEG,XEND,NDIV,XINC,S,TRANS,TBEG,TEND,TINC,Q,HEAD
C
CALL INPUT
C
C
icnt=1
DO 150 TIM = TBEG, TEND, TINC
k=icnt
DO 100 IXDIST = 0,mx
DO 100 IYDIST = 0,ny
v(ixdist,iydist)=0.0
100 CONTINUE
DO 200 XDIST = XBEG, XEND, XINC
DO 200 YDIST = XBEG, XEND, XINC
R=(XDIST**2 + YDIST**2)**.5
U=R**2*S/4./TRANS/TIM
DRWDN= Q/4/3.1416/TRANS*WU(U)
IXDIST=XDIST
IYDIST=YDIST
C Conceptually, assign the matrix coordinates
C with the radius and the value is the drawdown
C v(IXDIST,IYDIST)=DRWDN
IR=R
IS=DRWDN*1.5
if(IR .gt. mx) IR=mx
if(IS .gt. mx) IS=mx
v(IR,IS)=HEAD-DRWDN
WRITE(6,900)IR,IS,HEAD-DRWDN,TIM
900 format ('IR=',I10,'IS=',I10,'H-D=',f10.2,'time=',f10.2)
200 CONTINUE
C
if(k .lt. 10) write(fn,102) k
if(k .ge. 10 .and. k .lt. 100) write(fn,103) k
if(k .ge. 100 ) write(fn,104) k
c
102 format('img00',i1,'.raw')
103 format('img0',i2,'.raw')
104 format('img',i3,'.raw')
do 33 m = 1,mmax
do 33 n= 1,nmax
iu(m,n)=255*(v(m,n)+1.)*.5
if(iu(m,n).lt.0) iu(m,n)=0
if(iu(m,n).gt.255) iu(m,n) = 255
cu(m,n) = char(iu(m,n))
33 continue
open(unit=2,file=fn,status='unknown',access='direct',
1 form='unformatted',recl=isize)
write(2,rec=1) ((cu(m,n),n=1,nmax),m=1,mmax)
close(2)
icnt=icnt+1
150 CONTINUE
c
77 continue
C
STOP
END
SUBROUTINE INPUT
C
C THIS SUBROUTINE USED TO ENTER THE DATA
C
COMMON /BLOCK1/ XBEG,XEND,NDIV,XINC,S,TRANS,TBEG,TEND,TINC,Q,HEAD
CHARACTER ANS*1, FNAME*12
C
150 WRITE(*,200)
200 FORMAT(//' ENTER BEGINNING DISTANCE (FEET):' $)
READ(*,*)XBEG
C
WRITE(*,250)
250 FORMAT(//' ENTER ENDING DISTANCE (FEET): '$)
READ(*,*)XEND
C
WRITE(*,300)
300 FORMAT(//' ENTER NUMBER OF DIVISIONS (XINC=[XEND-XBEG]/NDIV): '$)
READ(*,*)NDIV
XINC=(XEND-XBEG)/NDIV
C
WRITE(*,350)
350 FORMAT(//' ENTER STORAGE COEFFICIENT (DIMENSIONLESS): '$)
READ(*,*)S
C
WRITE(*,400)
400 FORMAT(//' ENTER AQUIFER TRANSMISSIVITY (FEET**2/DAY): '$)
READ(*,*)TRANS
C
WRITE(*,450)
450 FORMAT(//' ENTER WELL DISCHARGE (FEET**3/DAY): '$)
READ(*,*)Q
C
WRITE(*,470)
470 FORMAT(//' ENTER INITIAL HEAD (FEET): '$)
READ(*,*)HEAD
C
WRITE(*,500)
500 FORMAT(//' ENTER BEGINING TIME (DAYS): '$)
READ(*,*)TBEG
C
WRITE(*,550)
550 FORMAT(//' ENTER ENDING TIME (DAYS): '$)
READ(*,*)TEND
C
WRITE(*,600)
600 FORMAT(//' ENTER TIME INCREMENT (DAYS): '$)
READ(*,*)TINC
C
WRITE(*,650)XBEG,XEND,NDIV,S,TRANS,Q,HEAD,TBEG,TEND,TINC
650 FORMAT(//////
1' BEGINNING DISTANCE FROM PUMPED WELL (FEET) =',F8.1,/
2' ENDING DISTANCE FROM PUMPED WELL (FEET) =',F8.1,/
3' NUMBER OF DIVISIONS (XINC=[XEND-XBEG]/NDIV) =',I6,/
4' STORAGE COEFFICIENT (DIMENSIONLESS) =',E9.2,/
5' AQUIFER TRANSMISSIVITY (FEET**2/DAY) =',F8.1,/
6' WELL DISCHARGE (FEET**3/DAY) =',F8.1,/
6' INITIAL HEAD (FEET) =',F8.1,/
7' BEGINING TIME (DAYS) =',F8.1,/
8' ENDING TIME (DAYS) =',F8.1,/
9' TIME INCREMENT (DAYS) =',F8.1,/)
675 WRITE(*,700)
700 FORMAT(' IS THIS DATA CORRECT? (Y/N): '$)
C
READ(*,750)ANS
750 FORMAT(A)
IF(ANS .EQ. 'N' .OR. ANS .EQ. 'n')THEN
WRITE(*,800)
800 FORMAT(//' PLEASE RE-ENTER THE DATA.')
GOTO 150
ELSE
IF(TBEG .EQ. 0.)THEN
TBEG=.00001
TEND=TEND+.00001
ENDIF
IF(XBEG .EQ. 0.)THEN
XBEG=.00001
XEND=XEND+.00001
XINC=(XEND-XBEG)/NDIV
ENDIF
RETURN
ENDIF
c
RETURN
END
FUNCTION WU(U)
C
C THIS FUNCTION SUBPROGRAM SOLVES FOR THE WELL FUNCTION
C
IF(U .GE. 100.)THEN
WU=0.
GOTO 999
ENDIF
C
IF(U .LE. 1.0)THEN
A0 = -.57721566
A1 = .99999193
A2 = -.24991055
A3 = .05519968
A4 = -.00976004
A5 = .00107857
C
WU = A0 + A1*U + A2*U**2 + A3*U**3 + A4*U**4 + A5*U**5 -ALOG(U)
C
ELSE
A1= 8.5733287401
A2=18.0590169730
A3= 8.6347608925
A4= 0.2677737343
B1= 9.5733223454
B2=25.6329561486
B3=21.0996530827
B4= 3.9584969228
C
WU1 = U**4 + A1*U**3 + A2*U**2 + A3*U + A4
WU2 = U**4 + B1*U**3 + B2*U**2 + B3*U + B4
WU = WU1/WU2/U/EXP(U)
ENDIF
C
999 RETURN
END