PROGRAM MAIN * -------------------------------------------------- * Calls the Menu shell to activate and control all * processes involved in using this program * -------------------------------------------------- CALL MENU() END * ------------------- * END OF MAIN PROGRAM * ------------------- * -------------------------------- * PRETTY PRINT SUPPORT SUBROUTINES * -------------------------------- SUBROUTINE STARS(LUN) * -------------------------------------------------- * This allows to always make same size space lines * of stars. * -------------------------------------------------- INTEGER LUN WRITE(LUN,*) '**********************************************', + '**********************' RETURN END SUBROUTINE SPACES(sp,LUN) * -------------------------------------------------- * This allows to always make same space lines * -------------------------------------------------- INTEGER i,sp,LUN DO i=1,sp WRITE(LUN,*) ENDDO RETURN END * -------------------------- * CHARACTER SUPPORT ROUTINES * -------------------------- SUBROUTINE GETLEN(LINE,length) * -------------------------------------------------- * This routine gets the length of a line read in * from file. This then allows other subroutines to * parse the line with a known length. * -------------------------------------------------- CHARACTER*(*) LINE CHARACTER space,tab INTEGER i,length LOGICAL finished space=' ' tab='\t' i=1 length=0 finished=.false. DO WHILE ((i .le. LEN(LINE)) .and. .not. finished) IF (LINE((LEN(LINE)-i):LEN(LINE)-i+1) .eq. space .or. + LINE((LEN(LINE)-i):LEN(LINE)-i+1) .eq. tab) THEN length=LEN(LINE)-i ELSE length=LEN(LINE)-i finished=.true. ENDIF i=i+1 ENDDO RETURN END SUBROUTINE PARSELINE(LUN,columns,error) * -------------------------------------------------- * This subroutine takes a line from the common block * named txtblk (which was read in from file) and * breaks it into smaller strings. The delimiters * are spaces and tabs. Multiple space and/or tabs * automagically accounted and removed. * -------------------------------------------------- INCLUDE "text_block.inc" CHARACTER*255 Line INTEGER LUN,error,LinLen,space_index,tab_index INTEGER placement,columns LOGICAL finished CALL READLINE(LUN,Line,error) IF (error .ne. 0) THEN columns=0 GOTO 999 ELSE finished = .false. columns=0 CALL GETLEN(Line,LinLen) IF (LinLen .eq. 0) THEN columns=0 GOTO 999 ELSE DO WHILE (.not. finished) 10 CALL GETLEN(Line,LinLen) tab_index=INDEX(Line,'\t') space_index=INDEX(Line,' ') IF (LinLen .eq. 0) THEN columns=0 GOTO 999 ELSE IF (space_index .eq. 0 .and. tab_index .eq. 0) THEN IF (LinLen .gt. 0) THEN columns=columns+1 txtblk(columns)=Line(1:LinLen) ENDIF ELSE IF (space_index .eq. 1 .or. + tab_index .eq. 1) THEN Line=Line(2:LinLen) GOTO 10 ELSE IF (space_index .eq. 0 .or. + tab_index .eq. 0) THEN placement=MAX(tab_index,space_index) columns=columns+1 txtblk(columns)=Line(1:placement-1) IF ((placement-1) .ge. LinLen) THEN GOTO 999 ENDIF Line=Line(placement+1:LinLen) GOTO 10 ELSE placement=MIN(tab_index,space_index) columns=columns+1 txtblk(columns)=Line(1:placement-1) Line=Line(placement+1:LinLen) GOTO 10 ENDIF ENDIF finished=.true. ENDDO ENDIF ENDIF 999 RETURN END * -------------------------- * NUMERICAL SUPPORT ROUTINES * -------------------------- SUBROUTINE FCDECI(strpiece,value) * -------------------------------------------------- * This subroutine takes a string of digits and * converts it into a REAL numerical value. It allows * for exponential writing (E) and accounts for * decimals as well. * -------------------------------------------------- INTEGER bgpt,endpt,cutoff,decimal,lenpiece,exponent,mantissa INTEGER whpt,frpt,j,tempV REAL pw,pw1,value CHARACTER*(*) strpiece CALL GETLEN(strpiece,lenpiece) endpt=lenpiece bgpt=1 cutoff=0 pw=1.0 pw1=1.0 mantissa=0 decimal=INDEX(strpiece,'.') exponent=INDEX(strpiece,'e') IF (exponent .eq. 0) THEN exponent=INDEX(strpiece,'E') ENDIF IF (decimal .eq. 0 .and. exponent .eq. 0) THEN IF (lenpiece .gt. 8) THEN bgpt=MAX(endpt-8,1) ENDIF CALL ICDECI(strpiece,bgpt,endpt,tempV) value=tempV*1.0 ELSE IF (exponent .eq. 0) THEN IF ((decimal-9) .gt. 8) THEN bgpt=MAX(decimal-9,1) ENDIF IF ((lenpiece-decimal) .gt. 8) THEN endpt=MIN(decimal+9,lenpiece) ENDIF CALL ICDECI(strpiece,bgpt,decimal-1,whpt) CALL ICDECI(strpiece,decimal+1,endpt,frpt) DO j=1,(endpt-decimal) pw=pw*10.0 ENDDO value=whpt+frpt/pw ELSE IF ((decimal-9) .gt. 8) THEN bgpt=MAX(decimal-9,1) ENDIF IF ((decimal+9) .gt. 9) THEN cutoff=MIN(decimal+9,exponent-1) ELSE cutoff=exponent-1 ENDIF CALL ICDECI(strpiece,bgpt,decimal-1,whpt) CALL ICDECI(strpiece,decimal+1,cutoff,frpt) CALL ICDECI(strpiece,exponent+1,lenpiece,mantissa) DO j=1,(exponent-decimal-1) pw1=pw1*10 ENDDO DO j=1,ABS(mantissa) pw=pw*10 ENDDO IF (mantissa .le. 0) THEN value=(whpt+frpt/pw1)/pw ELSE value=(whpt+frpt/pw1)*pw ENDIF ENDIF RETURN END SUBROUTINE ICDECI(strpiece,start,stop,value) * -------------------------------------------------- INTEGER start,stop,value,j,useneg CHARACTER*(*) strpiece CHARACTER*15 ASCII REAL number,total total=0.0 useneg=0 DO j=start,stop ASCII=strpiece(j:j) IF (ASCII .eq. '-') THEN useneg=1 ELSE call GETNUM(ASCII,value) number=value*(10**(stop-j)) total=total+number ENDIF ENDDO IF (useneg .eq. 1) THEN value=-1*ANINT(total) ELSE value=ANINT(total) ENDIF RETURN END SUBROUTINE GETNUM(str,value) CHARACTER str INTEGER value IF (str .eq. '1') THEN value=1 ELSE IF (str .eq. '2') THEN value=2 ELSE IF (str .eq. '3') THEN value=3 ELSE IF (str .eq. '4') THEN value=4 ELSE IF (str .eq. '5') THEN value=5 ELSE IF (str .eq. '6') THEN value=6 ELSE IF (str .eq. '7') THEN value=7 ELSE IF (str .eq. '8') THEN value=8 ELSE IF (str .eq. '9') THEN value=9 ELSE IF (str .eq. '0') THEN value=0 ELSE value=0 ENDIF RETURN END * ----------------------------------------- * SPECIFIC NUMBER FUNCTION SUPPORT ROUTINES * ----------------------------------------- SUBROUTINE FILLFACTORIAL(value,result) INCLUDE "global_variables.inc" INTEGER value,i REAL result,temp result=1.0 DO i=2,value result=result*i ENDDO RETURN END SUBROUTINE BESSY1(x,result) REAL x,result,temp REAL xx,z REAL p1,p2,p3,p4,p5,q1,q2,q3,q4,q5, + r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,s7,y SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5, + r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,s7 DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4, + .2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0, + -.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/ DATA r1,r2,r3,r4,r5,r6/-.4900604943d13,.1275274390d13, + -.5153438139d11,.7349264551d9,-.4237922726d7,.8511937935d4/, + s1,s2,s3,s4,s5,s6,s7/.2499580570d14,.4244419664d12, + .3733650367d10,.2245904002d8,.1020426050d6,.3549632885d3,1.d0/ IF (x.lt.8.) THEN y=x**2 CALL BESSJ1(x,temp) result=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+ + y*(s4+y*(s5+y*(s6+y*s7))))))+.636619772* + (temp*ALOG(x)-1./x) ELSE z=8./x y=z**2 xx=x-2.356194491 result=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*(p3+y*(p4+y* * p5))))+z*cos(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5))))) ENDIF RETURN END SUBROUTINE BESSJ1(x,result) REAL x,result,temp REAL ax,xx,z REAL p1,p2,p3,p4,p5,q1,q2,q3,q4,q5, + r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,y SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5, + r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6 DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0, + 242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/, + s1,s2,s3,s4,s5,s6/144725228442.d0,2300535178.d0, + 18583304.74d0,99447.43394d0,376.9991397d0,1.d0/ DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4, + .2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0, + -.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/ IF (abs(x).lt.8.) THEN y=x**2 result=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y* + (s2+y*(s3+y*(s4+y*(s5+y*s6))))) ELSE ax=abs(x) z=8./ax y=z**2 xx=ax-2.356194491 result=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* + p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))* + sign(1.,x) ENDIF RETURN END SUBROUTINE BESSY0(x,result) REAL x,result,temp CU USES BESSJ0 REAL xx,z REAL p1,p2,p3,p4,p5,q1,q2,q3,q4,q5, + r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,y SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5, + r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6 DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4, + -.2073370639d-5,.2093887211d-6/, + q1,q2,q3,q4,q5/-.1562499995d-1,.1430488765d-3, + -.6911147651d-5,.7621095161d-6,-.934945152d-7/ DATA r1,r2,r3,r4,r5,r6/-2957821389.d0,7062834065.d0, + -512359803.6d0,10879881.29d0,-86327.92757d0,228.4622733d0/, + s1,s2,s3,s4,s5,s6/40076544269.d0,745249964.8d0,7189466.438d0, + 47447.26470d0,226.1030244d0,1.d0/ IF (x.lt.8.) THEN y=x**2 CALL BESSJ0(x,temp) result=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y* + (s3+y*(s4+y*(s5+y*s6)))))+.636619772*temp*ALOG(x) ELSE z=8./x y=z**2 xx=x-.785398164 result=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*(p3+y*(p4+y* + p5))))+z*cos(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5))))) ENDIF RETURN END SUBROUTINE BESSJ0(x,result) REAL x,result REAL ax,xx,z REAL p1,p2,p3,p4,p5,q1,q2,q3,q4,q5, + r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,y SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5, + r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6 DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4, + -.2073370639d-5,.2093887211d-6/, + q1,q2,q3,q4,q5/-.1562499995d-1,.1430488765d-3,-.6911147651d-5, + .7621095161d-6,-.934945152d-7/ DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0, + 651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/, + s1,s2,s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0, + 59272.64853d0,267.8532712d0,1.d0/ IF (abs(x).lt.8.) THEN y=x**2 result=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y* + (s3+y*(s4+y*(s5+y*s6))))) ELSE ax=abs(x) z=8./ax y=z**2 xx=ax-.785398164 result=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* + p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5))))) ENDIF RETURN END SUBROUTINE WELLFUNC(x,alpha,u,result) REAL x,alpha,u,result REAL PI,temp1,temp2,temp3,temp4 PARAMETER (PI=3.141592653) CALL BESSJ0(x,temp1) CALL BESSJ1(x,temp2) CALL BESSY0(x,temp3) CALL BESSY1(x,temp4) delx = (x*temp1 - 2*alpha*temp2)**2 + + (x*temp3 - 2*alpha*temp4)**2 result = 32*alpha**2/PI**2*(1 - EXP(-x**2/4/u))/ + (x**3*delx) RETURN END * -------------------- * I/O SUPPORT ROUTINES * -------------------- SUBROUTINE INITLUN INCLUDE "lun_block.inc" INTEGER loop IF (lunint .ne. 1) THEN DO loop=1,maxlun lunblk(loop)=-1 ENDDO lunint=1 ENDIF RETURN END SUBROUTINE GETLUN(LUN) INCLUDE "lun_block.inc" INTEGER i,LUN LOGICAL found found=.false. i=1 DO WHILE ((i .le. maxlun) .and. (.not. found)) IF (lunblk(i) .eq. -1) THEN found=.true. LUN=i ELSE i=i+1 ENDIF ENDDO IF (.not. found) THEN LUN=0 ENDIF RETURN END SUBROUTINE RESLUN(LUN,Reser,error,fn) INCLUDE "lun_block.inc" INTEGER LUN,Reser,error CHARACTER*(*) fn error=0 IF (Reser .eq. 1) THEN 10 CALL GETLUN(LUN) IF (LUN .eq. 0) THEN error=-1 ELSE IF (LUN .eq. 5 .or. LUN .eq. 6) THEN lunblk(LUN) = 1 GOTO 10 ELSE lunblk(LUN)=1 luntxt(LUN)=fn nmflop=nmflop+1 ENDIF ELSE lunblk(LUN)=-1 luntxt(LUN)=' ' ENDIF RETURN END SUBROUTINE READLINE(LUN,Line,error) INTEGER LUN,error CHARACTER*255 Line IF (LUN .gt. 0) THEN READ (LUN,10,IOSTAT=error) Line 10 FORMAT(a) ELSE error=-1 ENDIF RETURN END * --------------------- * FILE SUPPORT ROUTINES * --------------------- SUBROUTINE CLOSEFILE(LUN,fn,error) INCLUDE "lun_block.inc" CHARACTER*(*) fn INTEGER LUN,error error=0 IF (LUN .gt. 0) THEN CLOSE(LUN) CALL RESLUN(LUN,-1,error,fn) nmflop=nmflop-1 ELSE error=-1 ENDIF RETURN END SUBROUTINE OPENPUMPFILE(tdpts,xdata,ydata,sigma,error) INCLUDE "text_block.inc" INCLUDE "lun_block.inc" INCLUDE "global_variables.inc" CHARACTER*50 file CHARACTER*255 lindat INTEGER tdpts,LinLen,Length INTEGER error,cols,datent REAL xdata(maxdat),ydata(maxdat),sigma(maxdat) CALL SPACES(1,6) CALL STARS(6) WRITE(6,*) 'Enter the name of the file to use as the pump data.' WRITE(6,*) 'It must reside in the same directory as this program.' CALL STARS(6) CALL SPACES(1,6) READ(5,*) file CALL GETLEN(file,LinLen) CALL INITLUN CALL OPENEXISTFILE(file,LUN,error) IF (error .ne. 0) THEN CALL SPACES(1,6) CALL GETLEN(file,Length) WRITE(6,*) 'Problem opening file: ',file(1:Length) CALL REPORTERROR(error) GOTO 999 ENDIF error=0 datent=0 DO WHILE (error .eq. 0) CALL PARSELINE(LUN,cols,error) datent=datent+1 IF (error .lt. 0) THEN cols=0 ELSE CALL PUTPUMPDATA(xdata,ydata,sigma,datent) ENDIF ENDDO tdpts=datent-1 CALL SPACES(1,6) 999 CALL CLOSEFILE(LUN,file,error) RETURN END c SUBROUTINE OPENDATACURVE(time,head,error) c INCLUDE "text_block.inc" c INCLUDE "lun_block.inc" c INCLUDE "global_variables.inc" c CHARACTER*50 file c CHARACTER*255 lindat c INTEGER LinLen,Length c INTEGER error,cols,datent c REAL time(maxdat),head(maxdat,2) c CALL SPACES(1,6) c CALL STARS(6) c WRITE(6,*) 'Enter name of the file to use for data.' c WRITE(6,*) 'It must reside in the same directory as this program.' c CALL STARS(6) c CALL SPACES(1,6) c READ(5,*) file c CALL GETLEN(file,LinLen) c CALL INITLUN c CALL OPENEXISTFILE(file,LUN,error) c IF (error .ne. 0) THEN c CALL SPACES(1,6) c CALL GETLEN(file,Length) c WRITE(6,*) 'Problem opening file: ',file(1:Length) c CALL REPORTERROR(error) c GOTO 999 c ENDIF c error=0 c datent=0 c DO WHILE (error .eq. 0) c CALL PARSELINE(LUN,cols,error) c datent=datent+1 c IF (error .lt. 0) THEN c cols=0 c ELSE c CALL PUTDATATABLE(time,head,datent) c ENDIF c ENDDO c thmax=datent-1 c CALL SPACES(1,6) c 999 CALL CLOSEFILE(LUN,file,error) c IF (error .ne. 0 .and. thname .ne. ' ') THEN c CALL OPENEXISTFILE(thname,LUN,error) c ENDIF c thname=file c RETURN c END SUBROUTINE OPENALPHAVALUES(Wu,error) INCLUDE "text_block.inc" INCLUDE "lun_block.inc" INCLUDE "global_variables.inc" CHARACTER*50 file CHARACTER*255 lindat INTEGER LinLen,Length INTEGER error,cols,datent REAL Wu(maxudat,maxalph) CALL SPACES(1,6) CALL STARS(6) WRITE(6,*) 'Enter name of the file to use for alpha values.' WRITE(6,*) 'It must reside in the same directory as this program.' CALL STARS(6) CALL SPACES(1,6) READ(5,*) file CALL GETLEN(file,LinLen) CALL INITLUN CALL OPENEXISTFILE(file,LUN,error) IF (error .ne. 0) THEN CALL SPACES(1,6) CALL GETLEN(file,Length) WRITE(6,*) 'Problem opening file: ',file(1:Length) CALL REPORTERROR(error) GOTO 999 ENDIF error=0 datent=0 DO WHILE (error .eq. 0) CALL PARSELINE(LUN,cols,error) datent=datent+1 IF (error .lt. 0) THEN cols=0 ELSE CALL PUTALPHATABLE(Wu,datent) ENDIF ENDDO CALL SPACES(1,6) 999 CALL CLOSEFILE(LUN,file,error) IF (error .ne. 0 .and. alname .ne. ' ') THEN CALL OPENEXISTFILE(alname,LUN,error) ENDIF alname=file RETURN END SUBROUTINE OPENTYPECURVE(Wu,error) INCLUDE "text_block.inc" INCLUDE "lun_block.inc" INCLUDE "global_variables.inc" CHARACTER*50 file CHARACTER*255 lindat INTEGER LinLen,Length INTEGER error,cols,datent REAL Wu(maxudat,maxalph) CALL SPACES(1,6) CALL STARS(6) WRITE(6,*) 'Enter name of the file to use for type curves.' WRITE(6,*) 'It must reside in the same directory as this program.' CALL STARS(6) CALL SPACES(1,6) READ(5,*) file CALL GETLEN(file,LinLen) CALL INITLUN CALL OPENEXISTFILE(file,LUN,error) IF (error .ne. 0) THEN CALL SPACES(1,6) CALL GETLEN(file,Length) WRITE(6,*) 'Problem opening file: ',file(1:Length) CALL REPORTERROR(error) GOTO 999 ENDIF error=0 datent=0 DO WHILE (error .eq. 0) CALL PARSELINE(LUN,cols,error) datent=datent+1 IF (error .lt. 0) THEN cols=0 ELSE CALL PUTCURVETABLE(Wu,cols,datent) ENDIF ENDDO umax=datent-1 CALL SPACES(1,6) 999 CALL CLOSEFILE(LUN,file,error) IF (error .ne. 0 .and. tcname .ne. ' ') THEN CALL OPENEXISTFILE(tcname,LUN,error) ENDIF tcname=file RETURN END SUBROUTINE OPENEXISTFILE(File_Name,LUN,error) INCLUDE "lun_block.inc" LOGICAL file_exists CHARACTER*(*) File_Name INTEGER error,LUN,Length error=0 CALL GETLEN(File_Name,Length) INQUIRE(FILE=File_Name(1:Length),EXIST=file_exists) IF (.not. file_exists) THEN error=-1 ELSE CALL INITLUN CALL RESLUN(LUN,1,error,File_Name) IF (error .eq. 0) THEN OPEN(UNIT=LUN,FILE=File_Name(1:Length), + STATUS='UNKNOWN',IOSTAT=error) ELSE CALL CLOSEFILE(LUN,File_Name,error) ENDIF ENDIF RETURN END SUBROUTINE GETFILENAME(fn) CHARACTER*(*) fn CALL SPACES(1,6) CALL STARS(6) WRITE(6,*) 'Please enter the file name you wish to save to:' CALL STARS(6) CALL SPACES(1,6) READ(5,*) fn RETURN END SUBROUTINE OPENNEWFILE(fn,LUN,error) INCLUDE "lun_block.inc" LOGICAL file_exists CHARACTER*(*) fn INTEGER error,LUN,Ferror,Length error=0 Ferror=0 CALL INITLUN CALL RESLUN(LUN,1,error,fn) IF (error .eq. 0) THEN CALL GETLEN(fn,Length) OPEN (UNIT=LUN,FILE=fn(1:Length), + STATUS='UNKNOWN',IOSTAT=error) IF (error .ne. 0) THEN Ferror=error CALL CLOSEFILE(LUN,fn,error) ENDIF ELSE CALL CLOSEFILE(LUN,fn,error) ENDIF error=MAX(error,Ferror) RETURN END * --------------------------------- * PROGRAM SPECIFIC SUPPORT ROUTINES * --------------------------------- SUBROUTINE REPORTERROR(error) INTEGER error WRITE (6,*) 'There was an error encountered with error code',error RETURN END SUBROUTINE PUTALPHATABLE(Wu,datent) INCLUDE "text_block.inc" INCLUDE "global_variables.inc" INTEGER datent,cols,i REAL Wu(maxudat,maxalph) CALL FCDECI(txtblk(1),Wu(1,1+datent)) alpmax=datent+1 RETURN END SUBROUTINE PUTCURVETABLE(Wu,cols,datent) INCLUDE "text_block.inc" INCLUDE "global_variables.inc" INTEGER datent,cols,i REAL Wu(maxudat,maxalph) IF (datent .eq. 1) THEN DO i=1,cols CALL FCDECI(txtblk(i),Wu(datent,i+1)) ENDDO alpmax=cols+1 ELSE DO i=1,cols CALL FCDECI(txtblk(i),Wu(datent,i)) ENDDO ENDIF RETURN END SUBROUTINE PUTDATATABLE(time,head,datent) INCLUDE "text_block.inc" INCLUDE "global_variables.inc" INTEGER datent,cols,i REAL time(maxdat),head(maxdat,2) IF (datent .gt. 1) THEN CALL FCDECI(txtblk(1),time(datent)) CALL FCDECI(txtblk(2),head(datent,1)) ELSE CALL FCDECI(txtblk(1),time(1)) ENDIF RETURN END SUBROUTINE PUTPUMPDATA(xdata,ydata,sigma,datent) INCLUDE "text_block.inc" INCLUDE "global_variables.inc" INTEGER datent,cols,i REAL xdata(maxdat),ydata(maxdat),sigma(maxdat) CALL FCDECI(txtblk(1),xdata(datent)) CALL FCDECI(txtblk(2),ydata(datent)) CALL FCDECI(txtblk(3),sigma(datent)) RETURN END SUBROUTINE PRINTTYPECURVES(Wu) INCLUDE "text_block.inc" INCLUDE "global_variables.inc" INTEGER i,j REAL Wu(maxudat,maxalph) DO i=2,umax WRITE(6,*) Wu(i,1) DO j=2,alpmax WRITE(6,*) Wu(i,j),i,j ENDDO ENDDO RETURN END SUBROUTINE PRINTALPHAVALUES(Wu) INCLUDE "text_block.inc" INCLUDE "global_variables.inc" INTEGER i,j REAL Wu(maxudat,maxalph) WRITE (6,*) 'Here are the alphas' DO j=2,alpmax WRITE(6,*) Wu(1,j),j ENDDO RETURN END c SUBROUTINE PRINTDATACURVE(time,head) c INCLUDE "global_variables.inc" c INTEGER i c REAL time(maxdat),head(maxdat,2) c DO i=2,thmax c WRITE(6,*) time(i),head(i,1),head(i,2) c ENDDO c RETURN c END c SUBROUTINE GETPUMPTIME(time,timcp) c INCLUDE "global_variables.inc" c REAL time(maxdat),timcp c timcp=time(1) c RETURN c END SUBROUTINE MAKETYPECURVES(Wu,radius,trnsmt,tfinal,tpump) INCLUDE "global_variables.inc" INTEGER numdata REAL Wu(maxudat,maxalph),tfinal,tpump,increment REAL value,trnsmt,radius,time(maxudat) IF (tfinal .gt. 0 .and. tfinal .le. 10) THEN DO i=0,ANINT(tfinal/1.0) time(i+1) = i*1.0 ENDDO ELSE IF (tfinal .gt. 10 .and. tfinal .le. 100) THEN DO i=0,9 time(i+1) = i*1.0 ENDDO DO i=0,ANINT((tfinal-10.0)/5.0) time(i+1+10) = i*5.0+10.0 ENDDO ELSE IF (tfinal .gt. 100 .and. tfinal .le. 500) THEN DO i=0,9 time(i+1) = i*1.0 ENDDO DO i=0,17 time(i+1+10) = i*5.0+10.0 ENDDO DO i=0,ANINT((tfinal-100.0)/10.0) time(i+1+28) = i*10.0+100.0 ENDDO ELSE IF (tfinal .gt. 500 .and. tfinal .le. 1000) THEN DO i=0,9 time(i+1) = i*1.0 ENDDO DO i=0,17 time(i+1+10) = i*5.0+10.0 ENDDO DO i=0,39 time(i+1+28) = i*10.0+100.0 ENDDO DO i=0,ANINT((tfinal-500.0)/25.0) time(i+1+68) = i*25.0+500.0 ENDDO ELSE IF (tfinal .gt. 1000 .and. tfinal .le. 5000) THEN DO i=0,9 time(i+1) = i*1.0 ENDDO DO i=0,17 time(i+1+10) = i*5.0+10.0 ENDDO DO i=0,39 time(i+1+28) = i*10.0+100.0 ENDDO DO i=0,19 time(i+1+68) = i*25.0+500.0 ENDDO DO i=0,ANINT((tfinal-1000.0)/100.0) time(i+1+88) = i*100.0+1000.0 ENDDO ELSE DO i=0,9 time(i+1) = i*1.0 ENDDO DO i=0,17 time(i+1+10) = i*5.0+10.0 ENDDO DO i=0,39 time(i+1+28) = i*10.0+100.0 ENDDO DO i=0,19 time(i+1+68) = i*25.0+500.0 ENDDO DO i=0,39 time(i+1+88) = i*100.0+1000.0 ENDDO DO i=0,ANINT((tfinal-5000.0)/250.0) time(i+1+128) = i*250.0+5000.0 ENDDO ENDIF CALL SPACES(2,6) WRITE (6,*) 'Beginning type curve calculations.', + ' This may take awhile...' WRITE (6,*) '\n Number of points per type curve is ',thmax CALL SPACES(2,6) CALL CALCU(time,tpump,radius,trnsmt,Wu) RETURN END * ------------------------------------ * NUMERICAL TECHNIQUE SUPPORT ROUTINES * ------------------------------------ SUBROUTINE TRAPZN(upper,lower,N,alpha,u,result) INTEGER N,TrapIt,j REAL x,tnm,sum,del,upper,lower,result REAL alpha,u,above,below,temp,temp2 IF (upper .eq. lower) THEN result=0.0 RETURN ELSE IF (N .eq. 1) THEN CALL WELLFUNC(upper,alpha,u,above) CALL WELLFUNC(lower,alpha,u,below) result = 0.5*(upper-lower)* + (above-below) ELSE temp2 = 2**N TrapIt = ANINT(temp2) tnm = TrapIt del = (upper-lower)/tnm x = 0.5*del sum = 0.0 DO j=1,TrapIt temp=0 CALL WELLFUNC(x,alpha,u,temp) sum = sum + temp x = x + del ENDDO result = (upper-lower)*sum/tnm RETURN ENDIF ENDIF END SUBROUTINE TRAPZ(upper,lower,alpha,u,result) INTEGER j,TrapIt REAL upper,lower,alpha,u,maxintervals REAL storedInt,result,largereal REAL temp PARAMETER (largereal=1e30) storedInt=-largereal result = 0.0 TrapIt = 30 DO j=0,TrapIt CALL TRAPZN(upper,lower,j,alpha,u,result) IF (ABS(result-storedInt) .lt. + abs(storedInt)/1e4) THEN GOTO 10 ELSE storedInt=result ENDIF ENDDO 10 result=storedInt RETURN END c SUBROUTINE NORMALIZEHEAD(head) c INCLUDE "global_variables.inc" c REAL head(maxudat,2) c head(1,1)=head(2,1) c DO i=2,thmax c head(i,2)=head(i,1)/head(1,1) c ENDDO c RETURN c END SUBROUTINE CALCU(time,timctp,radius,trsmt,Wu) INCLUDE "global_variables.inc" INTEGER i REAL Wu(maxudat,maxalph),time(maxudat) REAL ucrtz,ucrtg,ucrtp,timctz,timctg,timctp REAL Wucrtz,Wucrtg,Wucrtp REAL radius,trsmt,factor DO i=2,alpmax alpcrt=Wu(1,i) factor=radius**2.0*alpcrt/4/trsmt ucrtp=factor/timctp IF (ucrtp .gt. UPPER) THEN IF (ucrtp .gt. 1e7) THEN Wucrtp = 0.0 ELSE CALL TRAPZ(10.0,0.0,alpcrt,ucrtp,Wucrtp) ENDIF ELSE IF (ucrtp .lt. LOWER) THEN Wucrtp=-EULER_C-ALOG(ucrtp) ELSE CALL TRAPZ(10.0,0.0,alpcrt,ucrtp,Wucrtp) ENDIF DO j=1,thmax timctg=time(j) timctz=time(j)+timctp ucrtz=factor/timctz ucrtg=factor/timctg IF (ucrtz .gt. UPPER) THEN IF (ucrtz .gt. 1e7) THEN Wucrtz = 0.0 ELSE CALL TRAPZ(10.0,0.0,alpcrt,ucrtz,Wucrtz) ENDIF ELSE IF (ucrtz .lt. LOWER) THEN Wucrtz=-EULER_C-ALOG(ucrtz) ELSE CALL TRAPZ(10.0,0.0,alpcrt,ucrtz,Wucrtz) ENDIF IF (ucrtg .gt. UPPER) THEN IF (ucrtg .gt. 1e7 .or. timctg .eq. 0.0) THEN Wucrtg = 0.0 ELSE CALL TRAPZ(10.0,0.0,alpcrt,ucrtg,Wucrtg) ENDIF ELSE IF (ucrtg .lt. LOWER) THEN Wucrtg=-EULER_C-ALOG(ucrtg) ELSE CALL TRAPZ(10.0,0.0,alpcrt,ucrtg,Wucrtg) ENDIF Wu(j+1,1) = time(j) Wu(j+1,i) = (Wucrtz - Wucrtg)/Wucrtp ENDDO WRITE (6,*) 'Finished calculating for alpha = ',alpcrt ENDDO RETURN END c SUBROUTINE EXTRAPOLATE(time,head,radius,trsmt,Wu,NewWu) c INCLUDE "global_variables.inc" c INTEGER i c REAL Wu(maxudat,maxalph),time(maxdat),head(maxdat,2) c REAL ulow,uhigh,alpcrt c REAL Wulow,Wuhigh c REAL ucrtz,ucrtg,ucrtp,timctz,timctg,timctp c REAL Wucrtz,Wucrtg,Wucrtp c REAL radius,trsmt,NewWu(maxdat,maxalph,6) c CALL GETPUMPTIME(time,timctp) c DO i=2,alpmax c alpcrt = Wu(1,i) c DO j=2,thmax c timctg=time(j) c timctz=time(j)+timctp c ucrtz=radius**2.0*alpcrt/4/trsmt/timctz c ucrtg=radius**2.0*alpcrt/4/trsmt/timctg c ucrtp=radius**2.0*alpcrt/4/trsmt/timctp c IF (ucrtz .gt. UPPER) THEN c Wucrtz=alpcrt/ucrtz c ELSE IF (ucrtz .lt. LOWER) THEN c Wucrtz=-EULER_C-ALOG(ucrtz) c ELSE c CALL FINDUVALUE(Wu,alpcrt,ucrtz,uhigh,ulow,Wuhigh,Wulow) c CALL CONNECTLINE(ucrtz,ulow,uhigh,Wulow,Wuhigh,Wucrtz) c ENDIF c IF (ucrtg .gt. UPPER) THEN c Wucrtg=alpcrt/ucrtg c ELSE IF (ucrtg .lt. LOWER) THEN c Wucrtg=-EULER_C-ALOG(ucrtg) c ELSE c CALL FINDUVALUE(Wu,alpcrt,ucrtg,uhigh,ulow,Wuhigh,Wulow) c CALL CONNECTLINE(ucrtg,ulow,uhigh,Wulow,Wuhigh,Wucrtg) c ENDIF c IF (ucrtp .gt. UPPER) THEN c Wucrtp=alpcrt/ucrtp c ELSE IF (ucrtp .lt. LOWER) THEN c Wucrtp=-EULER_C-ALOG(ucrtp) c ELSE c CALL FINDUVALUE(Wu,alpcrt,ucrtp,uhigh,ulow,Wuhigh,Wulow) c CALL CONNECTLINE(ucrtp,ulow,uhigh,Wulow,Wuhigh,Wucrtp) c ENDIF c NewWu(1,i,1)=alpcrt c NewWu(j,i,1)=ucrtz c NewWu(j,i,2)=Wucrtz c NewWu(j,i,3)=ucrtg c NewWu(j,i,4)=Wucrtg c NewWu(j,i,5)=ucrtp c NewWu(j,i,6)=Wucrtp c ENDDO c ENDDO c RETURN c END c SUBROUTINE FINDUVALUE(Wu,alpcrt,ucrt,uhigh,ulow,Wuhigh,Wulow) c INCLUDE "global_variables.inc" c INTEGER i c REAL Wu(maxudat,maxalph),alpcrt,ucrt,uhigh,ulow,Wuhigh,Wulow c uhigh=0.0 c ulow=0.0 c Wuhigh=0.0 c Wulow=0.0 c DO i=2,alpmax c IF (alpcrt .eq. Wu(1,i)) THEN c DO j=2,umax c IF (Wu(j,1) .ge. ucrt .and. ucrt .ge. Wu(j+1,1)) THEN c Wuhigh=Wu(j,i) c Wulow=Wu(j+1,i) c uhigh=Wu(j,1) c ulow=Wu(j+1,1) c ENDIF c ENDDO c ENDIF c ENDDO c RETURN c END c SUBROUTINE CONNECTLINE(ucrt,ulow,uhigh,Wulow,Wuhigh,extrap) c REAL ulow,uhigh,Wulow,Wuhigh,ucrt,temp,extrap c temp=0.0 c temp=ALOG10(Wuhigh)- c + ALOG10(Wuhigh/Wulow)*ALOG10(ucrt/uhigh)/ALOG10(ulow/uhigh) c extrap=10.0**(temp) c RETURN c END SUBROUTINE GETRADIUS(radius) INCLUDE "global_variables.inc" REAL radius CALL SPACES(1,6) CALL STARS(6) WRITE(6,*) 'Enter radius:' CALL STARS(6) CALL SPACES(1,6) READ(5,*) radius RETURN END SUBROUTINE GETPUMPTIME(tpump) INCLUDE "global_variables.inc" REAL tpump CALL SPACES(1,6) CALL STARS(6) WRITE(6,*) 'Enter pumping time:' CALL STARS(6) CALL SPACES(1,6) READ(5,*) tpump RETURN END SUBROUTINE GETTIMEPARMS(increment,tfinal) INCLUDE "global_variables.inc" REAL increment,tfinal c 10 CALL SPACES(1,6) c CALL STARS(6) c WRITE(6,*) 'Enter time increment:' c CALL STARS(6) c CALL SPACES(1,6) c READ(5,*) increment 10 CALL SPACES(1,6) CALL STARS(6) WRITE(6,*) 'Enter final time:' CALL STARS(6) CALL SPACES(1,6) READ(5,*) tfinal c IF (tfinal/increment .ge. maxudat) THEN c WRITE (6,*) 'Too many time steps' c GOTO 10 c ELSE c thmax = ANINT(tfinal/increment) c IF (thmax*increment .lt. tfinal) THEN c WRITE (6,*) 'Due to rounding, new time increment is ', c + tfinal/increment c ENDIF c increment = tfinal/thmax c ENDIF IF (tfinal .ge. 0. .and. tfinal .lt. 10.) THEN thmax=ANINT(tfinal/1.0)+1 ELSE IF (tfinal .ge. 10. .and. tfinal .lt. 100.) THEN thmax=10 thmax=thmax+ANINT((tfinal-10.0)/5.0)+1 ELSE IF (tfinal .ge. 100. .and. tfinal .lt. 500.) THEN thmax=10+18 thmax=thmax+ANINT((tfinal-100.0)/10.0)+1 ELSE IF (tfinal .ge. 500. .and. tfinal .lt. 1000.) THEN thmax=10+18+40 thmax=thmax+ANINT((tfinal-500.0)/25.0)+1 ELSE IF (tfinal .ge. 1000. .and. tfinal .lt. 5000.) THEN thmax=10+18+40+20 thmax=thmax+ANINT((tfinal-1000.0)/100.0)+1 ELSE IF (tfinal .ge. 5000.) THEN thmax=10+18+40+20+40 thmax=thmax+ANINT((tfinal-5000.0)/250.0)+1 ELSE thmax=0 ENDIF increment=0.0 RETURN END SUBROUTINE INITDATACURVE(time) INCLUDE "global_variables.inc" INTEGER i REAL time(maxudat) DO i=1,thmax time(i)=0.0 ENDDO RETURN END SUBROUTINE INITTYPECURVE(Wu) INCLUDE "global_variables.inc" INTEGER i REAL Wu(maxudat,maxalph) DO i=1,umax DO j=1,alpmax Wu(i,j)=0.0 ENDDO ENDDO RETURN END SUBROUTINE CHECKRESPONSE(choice,answer) INTEGER length,answer CHARACTER*(*) choice CALL GETLEN(choice,length) IF (length .gt. 1) THEN answer=0 ELSE CALL ICDECI(choice,1,1,answer) ENDIF RETURN END SUBROUTINE MENU INCLUDE "lun_block.inc" INCLUDE "text_block.inc" INCLUDE "global_variables.inc" INTEGER FLUN INTEGER error,value,length REAL xdata(maxdat),ydata(maxdat),sigma(maxdat) REAL Wu(maxudat,maxalph),time(maxdat) REAL alpcrt,radius,trsmt,result,deltat,tfinal,tpump CHARACTER*50 fn,choice value=1 alname=' ' FLUN=0 DO WHILE (value .ne. 8) length=0 CALL STARS(6) CALL SPACES(2,6) IF (FLUN .ne. 0) THEN CALL GETLEN(luntxt(FLUN),length) WRITE (6,*) 'Output file open: ', + luntxt(FLUN)(1:length) ELSE WRITE (6,*) 'Output file open: ' ENDIF CALL GETLEN(alname,length) WRITE (6,*) 'Alpha list loaded: ', + alname(1:length) IF (radius .gt. 0.0) THEN WRITE (6,100) radius ELSE WRITE (6,*) 'Radius value: ' ENDIF IF (tpump .gt. 0.0) THEN WRITE (6,110) tpump ELSE WRITE (6,*) 'Pumping time value: ' ENDIF WRITE (6,*) 'Time increment value: variable' IF (tfinal .gt. 0.0) THEN WRITE (6,130) tfinal ELSE WRITE (6,*) 'Final time value: ' ENDIF CALL SPACES(3,6) WRITE (6,*) 'This is the MAIN MENU of the PICKINGmodel code.' CALL SPACES(2,6) WRITE (6,*) 'Please select from the options below...' CALL SPACES(3,6) WRITE (6,*) '\t1.) OPEN output file' WRITE (6,*) '\t2.) CLOSE output file' WRITE (6,*) '\t3.) DESIGN a type-curve table' WRITE (6,*) '\t4.) Supply radius' WRITE (6,*) '\t5.) Supply pump time' WRITE (6,*) '\t6.) Supply time parameters' WRITE (6,*) '\t7.) CALCULATE type-curves' WRITE (6,*) '\t8.) END PICKINGmodel session' CALL SPACES(3,6) CALL STARS(6) CALL SPACES(1,6) WRITE (6,*) 'Please enter your choice:' READ(5,*) choice CALL CHECKRESPONSE(choice,value) IF (value .eq. 1) THEN CALL GETFILENAME(fn) CALL OPENNEWFILE(fn,FLUN,error) ELSE IF (value .eq. 2) THEN CALL CLOSEFILE(FLUN,fn,error) ELSE IF (value .eq. 3) THEN CALL INITDATACURVE(time) CALL OPENALPHAVALUES(Wu,error) ELSE IF (value .eq. 4) THEN CALL GETRADIUS(radius) ELSE IF (value .eq. 5) THEN CALL GETPUMPTIME(tpump) ELSE IF (value .eq. 6) THEN CALL GETTIMEPARMS(deltat,tfinal) ELSE IF (value .eq. 7) THEN CALL SPACES(1,6) CALL STARS(6) WRITE (6,*) 'Enter transmissivity value to use: ' CALL STARS(6) CALL SPACES(1,6) READ(5,*) trsmt CALL MAKETYPECURVES(Wu,radius,trsmt,tfinal,tpump) CALL OUTPUTTOFILE(FLUN,tpump,trsmt,Wu) ELSE IF (value .eq. 8) THEN CALL SPACES(3,6) WRITE (6,*) 'Ending program...' STOP ENDIF ENDDO 100 FORMAT (' ','Radius value: ',F10.6,' in') 110 FORMAT (' ','Pumping time value: ',F10.3,' sec') 120 FORMAT (' ','Time increment value: ',F10.3,' sec') 130 FORMAT (' ','Final time value: ',F10.3,' sec') RETURN END SUBROUTINE OUTPUTTOFILE(FLUN,tpump,trsmt,Wu) INCLUDE "global_variables.inc" INTEGER FLUN REAL Wu(maxudat,maxalph),trsmt REAL ap(maxalph+1) WRITE(FLUN,80) tpump WRITE(FLUN,*) '-------------------------------------------------', + '-------------' WRITE(FLUN,90) trsmt WRITE(FLUN,*) '-------------------------------------------------', + '-------------' CALL SPACES(1,FLUN) WRITE(FLUN,*) 'Here are the alpha values used' WRITE(FLUN,*) '-----------------------------------' WRITE(FLUN,100) Wu(1,2),Wu(1,3),Wu(1,4), + Wu(1,5),Wu(1,6) CALL SPACES(1,FLUN) WRITE(FLUN,*) 'Time values and calculated', + ' (W(u_z) - W(u_g))/W(u_p) values:' WRITE(FLUN,*) '------------------------------------------------', + '---------------------------' CALL SPACES(1,FLUN) CALL SPACES(1,FLUN) DO j=1,thmax ap(1)=Wu(j+1,1) DO i=2,alpmax ap(i)=(Wu(j+1,i)) ENDDO WRITE(FLUN,110) ap(1),ap(2),ap(3),ap(4),ap(5),ap(6) ENDDO 80 FORMAT (' ','\n Here is the Pumping time: ',F7.4) 90 FORMAT (' ','\n Here is the transmissivity value used: ', + (E15.7,3X)) 100 FORMAT (' ',13X,5(E10.3,3X)) 110 FORMAT (' ',6(F10.5,3X)) RETURN END