Program Main Integer JSTOP character*80 filen filen = 'vs2dh.fil' JSTOP = 0 C SET UP TO READ FROM DATA FILE GENERATED BY VS2DHI call SETUP(0, 0, 0, filen) 10 CALL STEP(JSTOP) IF (JSTOP.EQ.0) GOTO 10 stop end include 'd_modules.inc' SUBROUTINE SETUP(iold, ihydr, isorp, filen) C *** THIS SUBROUTINE IS THE FIRST HALF OF THE ORIGINAL MAIN CODE C *** THAT READS THE SIMULATION DATA C *** 12/1/98 C C C C****** CVSEXEC C****** C ----------------------------------------------------------------- C ************ PROGRAM VS2DH ********************* C C PROGRAM TO SOLVE FOR: C TWO DIMENSIONAL VERTICAL SECTION OR CYLINDRICAL THREE C DIMENSIONAL FLUID FLOW AND HEAT TRANSPORT UNDER C VARIABLY SATURATED CONDITIONS C C FLUID FLOW IS SOLVED FOR BY AN IMPLICIT FINITE DIFFERENCE C FORMULATION OF THE COMBINED RICHARDS AND COOPER-JACOB C EQUATIONS FOR FLUID CONTINUITY. C C ------- VERSION 3.2 August 2004 ---------- C ................................................................. C C DEFINITION OF FUNCTIONAL RELATIONSHIPS REQUIRED C VSHKU = RELATIVE HYDRAULIC CONDUCTIVITY AS A FUNCTION OF C PRESSURE HEAD C VSTHU = VOLUMETRIC MOISTURE CONTENT AS A FUNCTION OF PRESSURE HEAD C VSDTHU = FIRST DERIVATIVE OF MOISTURE CONTENT WITH RESPECT C TO PRESSURE HEAD C VSTHNV = PRESSURE HEAD AS A FUNCTION OF VOLUMETRIC MOISTURE C CONTENT C VSRDF = ROOT ACTIVITY AS A FUNCTION OF TIME AND DEPTH C VTRET = BULK DENSITY TIMES SLOPE OF ADSORPTION ISOTHERM. C C----------------------------------------------------------- C C SPECIFICATIONS FOR ARRAYS AND SCALARS C include 'd_rspac.inc' include 'd_kcon.inc' include 'd_mprop.inc' include 'd_press.inc' include 'd_disch.inc' include 'd_hcon.inc' include 'd_equat.inc' include 'd_jtxx.inc' include 'd_dumm.inc' include 'd_ptet.inc' include 'd_trxx.inc' include 'd_trxy1.inc' include 'd_pit.inc' include 'd_sip.inc' include 'd_plott.inc' include 'd_spfc.inc' include 'd_rprop.inc' include 'd_scon.inc' IMPLICIT DOUBLE PRECISION (A-H,P-Z) c include 'c_rspac.inc' c include 'c_kcon.inc' c include 'c_mprop.inc' c include 'c_press.inc' c include 'c_disch.inc' c include 'c_hcon.inc' c include 'c_equat.inc' c include 'c_jtxx.inc' c include 'c_dumm.inc' c include 'c_ptet.inc' c include 'c_trxx.inc' c include 'c_trxy1.inc' c include 'c_pit.inc' c include 'c_sip.inc' c include 'c_plott.inc' c include 'c_spfc.inc' c include 'c_rprop.inc' c include 'c_scon.inc' COMMON/ISPAC/NLY,NLYY,NXR,NXRR,NNODES COMMON/PND/POND COMMON/WGT/WUS,WDS COMMON/SCN1/TMPX,TMLT,DLTMX,DLTMIN,TRED COMMON/TCON/STIM,DSMAX,KTIM,NIT,NIT1,KP,NIT3 COMMON/JCON/JSTOP,JFLAG,jflag1 LOGICAL TRANS,TRANS1,SORP,SSTATE COMMON/TRXY/EPS1,EPS2,TRANS,TRANS1,SORP,SSTATE,MB9(72),NMB9 LOGICAL RAD,BCIT,ETSIM,SEEP,ITSTOP,CIS,CIT LOGICAL F7P,F11P,F8P,F9P,F6P,PRNT,o9p,o11p,o13p LOGICAL THPT,SPNT,PPNT,HPNT,VPNT COMMON/LOG1/RAD,BCIT,ETSIM,SEEP,ITSTOP,CIS,CIT COMMON/LOG2/F7P,F11P,F8P,F9P,F6P,PRNT,o9p,o11p,o13p COMMON/LOG4/THPT,SPNT,PPNT,HPNT,VPNT CHARACTER*80 TITL,filen,f5,f6,f7,f8,f9,f10,f11,f12,blank CHARACTER*4 ZUNIT,TUNIT,CUNX COMMON/SCHAR/TITL,ZUNIT,TUNIT,CUNX COMMON/MASSB/ BL(72),bcmft,bcmtt,bl29I,bl29IT,bl29O,bl29OT, 1bl68I,bl68IT,bl68o,bl68OT,bl66T integer hydraulicFunctionType, adsorptionType, iuHead, &iuConcentration common/functiontype/ hydraulicFunctionType, adsorptionType COMMON/VERSION/VER1P0 LOGICAL VER1P0 common/iuNumber/ iuHead, iuConcentration COMMON/TTT/ IFET,IFET1,NITT,NITT1,NITT2 real*8 ang include 'd_kdum.inc' c include 'c_kdum.inc' c *** Set version, hydraulic function type, and sorption type if (iold.eq.1) then VER1P0 = .true. hydraulicFunctionType = ihydr adsorptionType = isorp else VER1P0 = .false. endif c c open files c iuHead = 5 iuConcentration = 5 blank = ' ' open(2,file=filen) read(2,9090)f5 read(2,9090)f6 read(2,9090)f7 read(2,9090)f8 read(2,9090)f9 read(2,9090)f11 read(2,9090,end=9088)f10 iInitial = index(f10,'#') if (iInitial.ne.0) go to 9088 if(f5.ne.f10.and.f10.ne.blank) then open(10,file=f10) iuHead = 10 end if read(2,9090,end=9089)f12 iInitial = index(f12,'#') if (iInitial.ne.0) then f12 = f10 iuConcentration = iuHead go to 9088 end if 9089 continue if(f12.eq.blank) f12 = f10 if(f5.ne.f12.and.f12.ne.blank) then if(f10.ne.f12) then open(12,file=f12) iuConcentration = 12 else iuConcentration = 10 end if end if 9088 open(5,file=f5) open(6,file=f6) open(7,file=f7) open(8,file=f8) open(9,file=f9) open(11,file=f11) 9090 format(a80) C C------------------------------------------------------------------- C C ---- READ AND WRITE PROBLEM TITLE AND SPACE AND TIME CONSTANTS C JSTOP=0 READ (05,4000) TITL READ (5,*) TMAX,STIM,ANG READ (05,4010) ZUNIT,TUNIT,CUNX READ (05,*) NXR,NLY READ (05,*) NRECH,NUMT if(NUMT.lt.0) then numt = -numt o13p = .TRUE. else o13p = .FALSE. end if WRITE (06,4060) WRITE (06,4070) TITL,TMAX,TUNIT,STIM,NRECH,NUMT,NLY,NXR WRITE(06,4080) ANG IF(ANG.GT.90.0D0.OR.ANG.LT.-90.0D0)THEN WRITE(06,4090) jstop=2 return END IF READ (05,*) RAD,ITSTOP,TRANS IF(TRANS) READ(05,*)CIS,CIT SORP=.FALSE. READ (05,*) F11P,F7P,F8P,F9P,F6P READ (05,*) THPT,SPNT,PPNT,HPNT,VPNT WRITE (06,4100) F8P,ITSTOP,F7P,F11P,F9P,F6P WRITE (06,4110) THPT,SPNT,PPNT,HPNT,VPNT NLYY=NLY-1 NXRR=NXR-1 NNODES=NLY*NXR include 'd_arrays.inc' c include 'c_arrays.inc' C *** SET ALL ARRAYS TO ZERO DO 710 I=1, NLY DELZ(I) = 0.0D0 DZZ(I) = 0.0D0 710 CONTINUE DO 715 I=1, NXR DXR(I) = 0.0D0 RX(I) = 0.0D0 715 CONTINUE DO 711 I=1, NNODES HX(I) = 0.0D0 NTYP(I) = 0 THETA(I) = 0.0D0 THLST(I) = 0.0D0 P(I) = 0.0D0 PXXX(I) = 0.0D0 Q(I) = 0.0D0 QQ(I) = 0.0D0 HCND(I) = 0.0D0 HKLL(I) = 0.0D0 HKTT(I) = 0.0D0 A(I) = 0.0D0 B(I) = 0.0D0 C(I) = 0.0D0 D(I) = 0.0D0 E(I) = 0.0D0 RHS(I) = 0.0D0 XI(I) = 0.0D0 JTEX(I) = 0 DUM(I) = 0.0D0 DPTH(I) = 0.0D0 RT(I) = 0.0D0 DX1(I) = 0.0D0 DX2(I) = 0.0D0 DZ1(I) = 0.0D0 DZ2(I) = 0.0D0 VX(I) = 0.0D0 VZ(I) = 0.0D0 CC(I) = 0.0D0 COLD(I) = 0.0D0 CS(I) = 0.0D0 QT(I) = 0.0D0 NCTYP(I) = 0 RET(I) = 0.0D0 AO(I) = 0.0D0 BO(I) = 0.0D0 CO(I) = 0.0D0 DO(I) = 0.0D0 EO(I) = 0.0D0 PITT(I) = 0.0D0 RHO(I) = 0.0D0 RHOOLD(I) = 0.0D0 711 CONTINUE DO 712 I=1,72 BL(I)=0.0D0 712 CONTINUE bcmft = 0.0D0 bcmtt = 0.0D0 bl29I = 0.0D0 bl29IT = 0.0D0 bl29O = 0.0D0 bl29OT = 0.0D0 bl68I = 0.0D0 bl68IT = 0.0D0 bl68O = 0.0D0 bl68OT = 0.0D0 bl66T = 0.0D0 C C ESTABLISH HORIZONTAL OR RADIAL SPACING C READ (05,*) IFAC,FACX IF(IFAC.GT.0) GO TO 20 C C READ IN SPACING FOR EACH COLUMN C READ (05,*) (DXR(K),K=1,NXR) DO 10 K=1,NXR 10 DXR(K)=DXR(K)*FACX GO TO 60 20 IF(IFAC.EQ.2) GO TO 40 DO 30 K=1,NXR 30 DXR(K)=FACX GO TO 60 C C IF IFAC=2, HORIZONTAL NODE SPACING IS INCREMENTED BY A CONSTANT C MULTIPLIER UNTIL A USER-SPECIFIED MAXIMUM IS REACHED, WHERE- C UPON THE SPACING BECOMES CONSTANT C 40 READ (05,*) XMULT,XMAX DXR(1)=FACX DXR(2)=FACX DO 50 K=3,NXRR DXR(K)=DXR(K-1)*XMULT IF(DXR(K) .GT. XMAX)DXR(K)=XMAX 50 CONTINUE DXR(NXR)=DXR(NXRR) C C ESTABLISH VERTICAL SPACING C 60 READ (05,*) JFAC,FACZ IF(JFAC.GT.0) GO TO 80 C C READ IN VERTICAL SPACINGS INDIVIDUALLY C READ (05,*) (DELZ(K),K=1,NLY) DO 70 K=1,NLY 70 DELZ(K)=DELZ(K)*FACZ GO TO 120 80 IF(JFAC.EQ.2) GO TO 100 DO 90 K=1,NLY 90 DELZ(K)=FACZ GO TO 120 C C ESTABLISH VERTICAL SPACING BY PROGRESSION, AS ABOVE FOR HORIZ. C 100 READ (05,*) ZMULT,ZMAX DELZ(1)=FACZ DELZ(2)=FACZ DO 110 K=3,NLYY DELZ(K)=DELZ(K-1)*ZMULT IF(DELZ(K) .GT. ZMAX)DELZ(K)=ZMAX 110 CONTINUE DELZ(NLY)=DELZ(NLYY) 120 CONTINUE C C DETERMINE HORIZONTAL AND VERTICAL COORDINATES C RX(1)=-0.5D0 *DXR(1) DO 130 N=2,NXR RX(N)=RX(N-1)+0.5D0 *(DXR(N-1)+DXR(N)) 130 CONTINUE DZZ(1)=-0.5D0 *DELZ(1) DO 140 J=2,NLY 140 DZZ(J)=DZZ(J-1)+0.5D0 *(DELZ(J-1)+DELZ(J)) WRITE (06,4120) ZUNIT,(DELZ(K),K=1,NLY) WRITE (06,4130) ZUNIT,(DXR(K),K=1,NXR) PI=3.141592654D0 PI2=PI+PI ANG=ANG/360.0D0 IF(ANG.EQ.0.0D0) THEN CS1=1.0D0 CS2=0.0D0 ELSE IF(ANG.EQ.0.25D0.OR.ANG.EQ.-0.25D0) THEN CS1=0.0D0 ELSE CS1=DCOS(ANG*PI2) END IF CS2=-DSIN(ANG*PI2) END IF C C READ DATA FOR MONITORING TIMES AND POINTS C NPLT=0 o9p = .false. o11p = o9p IF(F8P) THEN READ (05,*) NPLT include 'd_obst.inc' c include 'c_obst.inc' IF(NPLT.EQ.0)NPLT=1 READ (05,*) (PLTIM(K),K=1,NPLT) WRITE (06,4140) tunit,(PLTIM(K),K=1,NPLT) else include 'd_obst.inc' END IF IF(F11P) THEN READ (05,*) NOBS if (nobs.lt.0) then nobs = -nobs o11p = .true. end if include 'd_obsp.inc' READ (05,*) ((KDUM(K,J),J=1,2),K=1,NOBS) WRITE (06,4150) ((KDUM(K,J),J=1,2),K=1,NOBS) c include 'c_obsp.inc' DO 150 K=1,NOBS N=NLY*(KDUM(K,2)-1)+KDUM(K,1) IJOBS(K)=N 150 continue END IF IF (F9P) THEN READ(05,*)NMB9 if (nmb9.lt.0) then nmb9 = -nmb9 o9p = .true. end if if (nmb9.gt.72) nmb9 = 72 READ(05,*) (MB9(K),K=1,NMB9) WRITE(06,4160) (MB9(K),K=1,NMB9) END IF PLTIM(NPLT+1)=TMAX+TMAX IF(RAD) THEN WRITE(06,4050) ELSE WRITE (06,4040) END IF IF(TRANS) THEN WRITE(06,4240) IF(CIS) THEN WRITE(6,4200) ELSE WRITE(6,4210) END IF IF(CIT) THEN WRITE(6,4220) ELSE WRITE(6,4230) END IF IF(SORP) WRITE(06,4250) END IF IF(F11P) then if (o13p) then if(trans) then WRITE (11,4033) TITL,TUNIT,ZUNIT,ZUNIT,ZUNIT,ZUNIT else WRITE (11,4031) TITL,TUNIT,ZUNIT,ZUNIT,ZUNIT,ZUNIT end if else if(trans) then WRITE (11,4032) TITL,TUNIT,ZUNIT,ZUNIT,ZUNIT,ZUNIT else WRITE (11,4030) TITL,TUNIT,ZUNIT,ZUNIT,ZUNIT,ZUNIT end if end if end if C C INITIALIZE CONSTANTS C ITEST=0 KTIM=0 NITT=0 NITT1=0 NITT2=0 JFLAG=1 KP=0 jplt = 0 WRITE (06,4170) C C C READ AND WRITE INITIAL VALUES OF PRESSURE HEAD, TOTAL HEAD, C THETA, AND SATURATION C ------------------------------------------------------------- C CALL VSREAD CALL VSSIP IFET=0 IFET2=0 CALL VSOUTP include 'd_kdumDealloc.inc' RETURN 4000 FORMAT(A80) 4010 FORMAT(4A4) 4020 FORMAT(5X,20(1H*),1X,31HDIMENSIONS TOO LARGE FOR ARRAYS, &1X,20(1H*)/5X,6HNLY = ,I5,2X,6H,NXR = ,I5) 4030 FORMAT(A80/21HMONITORING POINT FILE/2X,6HTIME, ,A4,2X, & 6H XR, ,A4,2X,6H Z, ,A4,2X,6H H, ,A4,2X,6H P, ,A4, & 2X,6H THETA,4X,8H SAT) 4031 FORMAT(A80/21HMONITORING POINT FILE/7X,6HTIME, ,A4,7X, & 6H XR, ,A4,10X,6H Z, ,A4,10X,6H H, ,A4,10X,6H P, ,A4, & 12X,6H THETA,12X,8H SAT) 4032 FORMAT(A80/21HMONITORING POINT FILE/2X,6HTIME, ,A4,2X, & 6H XR, ,A4,2X,6H Z, ,A4,2X,6H H, ,A4,2X,6H P, ,A4, & 2X,6H THETA,4X,8H SAT,9x,'TEMP') 4033 FORMAT(A80/21HMONITORING POINT FILE/7X,6HTIME, ,A4,7X, & 6H XR, ,A4,10X,6H Z, ,A4,10X,6H H, ,A4,10X,6H P, ,A4, & 12X,6H THETA,12X,8H SAT,17x,'TEMP') 4040 FORMAT(5X,32HCOORDINATE SYSTEM IS RECTANGULAR) 4050 FORMAT(5X,27HCOORDINATE SYSTEM IS RADIAL) 4060 FORMAT(35X,60('+')/35X,'+',26X,6H VS2DH,26X,'+'/35X, &'+',23x,' VERSION 3.2',23x,'+'/35x, &'+',11X,'SIMULATION OF 2-DIMENSIONAL VARIABLY',11X,'+'/ &35X,'+',12X,'SATURATED FLOW AND ENERGY TRANSPORT',11X,'+' &/35X,'+',11X,' THROUGH POROUS MEDIA. ',11X,'+' & /35X,60('+')//) 4070 FORMAT(//,1X,100(1H*)/5X,A80/1X,100(1H*)//10X, &24HSPACE AND TIME CONSTANTS/10X,23(1H-)/ & 5X,26HMAXIMUM SIMULATION TIME = ,E14.6,1X,A4/ &5X,'STARTING TIME = ',F10.4,/ &5X,28HNUMBER OF RECHARGE PERIODS =,I10/ &4X,32H MAXIMUM NUMBER OF TIME STEPS = ,I10/ &5X,17HNUMBER OF ROWS = ,I5/5X,20HNUMBER OF COLUMNS = ,I5) 4080 FORMAT(5X,'AXES TILTED BY ANGLE = ',F8.2) 4090 FORMAT(1X,'ANGLE OF AXES TILTING MUST BE BETWEEN -90 AND 90 ', &'DEGREES',/,1X,'SIMULATION TERMINATED') 4100 FORMAT(10X,16HSOLUTION OPTIONS/10X,16(1H-)/ &5X,'WRITE ALL PRESSURE HEADS TO FILE 8', &23H AT OBSERVATION TIMES? ,L1,/ &5X,'STOP SOLUTION IF MAXIMUM NO.', &' OF ITERATIONS EXCEEDED IN ANY TIME STEP? ',L1/5X, &'WRITE MAXIMUM CHANGE IN HEAD FOR EACH ITERATION TO FILE 7? ', &L1/5X,'WRITE RESULTS AT SELECTED OBSERVATION POINTS TO ', &'FILE 11? ', L1/,5X,'WRITE MASS BALANCE RATES TO FILE 9? ',L1/ &5X,'WRITE MASS BALANCE RATES TO FILE 6? ',L1) 4110 FORMAT(1H ,4X,35HWRITE MOISTURE CONTENTS TO FILE 6? ,L1/ & 5X,29HWRITE SATURATIONS TO FILE 6? ,L1/ & 5X,32HWRITE PRESSURE HEADS TO FILE 6? ,L1/ & 5X,29HWRITE TOTAL HEADS TO FILE 6? ,L1/ &5X,'WRITE VELOCITIES TO FILE 6? ',L1) 4120 FORMAT(50X,39HGRID SPACING IN VERTICAL DIRECTION, IN ,A4/ & (10(F10.3))) 4130 FORMAT(50X,47HGRID SPACING IN HORIZONTAL OR RADIAL DIRECTION, &,3H IN,1X,A4/(10F10.3)) 4140 FORMAT(5X,'TIMES (IN ',a4,') AT WHICH H WILL BE WRITTEN TO ', &'FILE 08'/10(5X,10E12.5,/)) 4150 FORMAT(5X,37HROW AND COLUMN OF OBSERVATION POINTS:/ & 10(2X,2I4)) 4160 FORMAT(5X,'MASS BALANCE COMPONENTS WRITTEN TO FILE 9', &/,5X,24I4) 4170 FORMAT(5X,36HMATRIX EQUATIONS TO BE SOLVED BY SIP) 4180 FORMAT(5X,100(1H*)/5X,17HEND OF SIMULATION/ & 5X,100(1H*)) 4190 FORMAT(5X,'TOTAL NUMBER OF ITERATIONS FOR FLOW EQUATION = ',I6 &/5X,'TOTAL NUMBER OF ITERATIONS FOR TRANSPORT EQUATION = ',I6, &/5x,'TOTAL NUMBER OF ITERATIONS FOR OUTER LOOP = ',I6) 4200 FORMAT(5X,'CENTRAL DIFFERENCING IN SPACE USED FOR TRANSPORT', &' EQUATION') 4210 FORMAT(4X,' BACKWARD DIFFERENCING IN SPACE USED FOR TRANSPORT', &' EQUATION') 4220 FORMAT(4X,' CENTRAL DIFFERENCING IN TIME USED FOR TRANSPORT', &' EQUATION') 4230 FORMAT(4X,' BACKWARD DIFFERENCING IN TIME USED FOR TRANSPORT', &' EQUATION') 4240 FORMAT(4X,' TRANSPORT TO BE SIMULATED') 4250 FORMAT(4X,' NONLINEAR SORPTION TO BE SIMULATED') END SUBROUTINE STEP(JSTP) C *** THIS IS THE SECOND HALF OF THE ORIGINAL MAIN CODE. THIS C *** PERFORMS A SINGLE TIME STEP. THE LOOPING IS DONE IN THE CALLING C *** PROGRAM. THE SUBROUTINE RETURNS THE SIMULATION TIME, TIME STEP, C *** AND THE FLAG JSTP. WHEN JSTP EQUALS 1, THIS SIGNIFIES C *** END OF SIMULATION. include 'd_rspac.inc' include 'd_kcon.inc' include 'd_mprop.inc' include 'd_press.inc' include 'd_disch.inc' include 'd_hcon.inc' include 'd_equat.inc' include 'd_jtxx.inc' include 'd_dumm.inc' include 'd_ptet.inc' include 'd_trxx.inc' include 'd_trxy1.inc' include 'd_plott.inc' include 'd_spfc.inc' include 'd_rprop.inc' include 'd_scon.inc' IMPLICIT DOUBLE PRECISION (A-H,P-Z) c include 'c_rspac.inc' c include 'c_kcon.inc' c include 'c_mprop.inc' c include 'c_press.inc' c include 'c_disch.inc' c include 'c_hcon.inc' c include 'c_equat.inc' c include 'c_jtxx.inc' c include 'c_dumm.inc' c include 'c_ptet.inc' c include 'c_trxx.inc' c include 'c_trxy1.inc' c include 'c_plott.inc' c include 'c_spfc.inc' c include 'c_rprop.inc' c include 'c_scon.inc' COMMON/ISPAC/NLY,NLYY,NXR,NXRR,NNODES COMMON/PND/POND COMMON/WGT/WUS,WDS COMMON/SCN1/TMPX,TMLT,DLTMX,DLTMIN,TRED COMMON/TCON/STIM,DSMAX,KTIM,NIT,NIT1,KP,NIT3 COMMON/JCON/JSTOP,JFLAG,jflag1 LOGICAL TRANS,TRANS1,SORP,SSTATE COMMON/TRXY/EPS1,EPS2,TRANS,TRANS1,SORP,SSTATE,MB9(72),NMB9 LOGICAL RAD,BCIT,ETSIM,SEEP,ITSTOP,CIS,CIT LOGICAL F7P,F11P,F8P,F9P,F6P,PRNT,o9p,o11p,o13p LOGICAL THPT,SPNT,PPNT,HPNT,VPNT COMMON/LOG1/RAD,BCIT,ETSIM,SEEP,ITSTOP,CIS,CIT COMMON/LOG2/F7P,F11P,F8P,F9P,F6P,PRNT,o9p,o11p,o13p COMMON/LOG4/THPT,SPNT,PPNT,HPNT,VPNT CHARACTER*80 TITL,filen,f5,f6,f7,f8,f9,f10,f11 CHARACTER*4 ZUNIT,TUNIT,CUNX COMMON/SCHAR/TITL,ZUNIT,TUNIT,CUNX integer hydraulicFunctionType, adsorptionType common/functiontype/ hydraulicFunctionType, adsorptionType COMMON/TTT/ IFET,IFET1,NITT,NITT1,NITT2 C C ------------------------------------------------------------- C START OF TIME LOOP C ------------------------------------------------------------- C 160 IF(JFLAG.EQ.1)IFET1=1 IF (JSTOP.GT.1) THEN JSTP=JSTOP RETURN ENDIF CALL VSTMER C C *** IF NO MORE PERIODS TO SIMULATE, JSTOP IS SET TO 1 IN VSTMER, C *** WE JUST RETURN AND LET THE CALLING PROGRAM TERMINATE THE C *** THE EXECUTION IF (JSTOP.GT.1) THEN JSTP=JSTOP RETURN ENDIF C TRANS1=.FALSE. c IF(.NOT.SSTATE) THEN C C SET UP AND SOLVE MATRIX EQUATIONS FOR FLOW C NIT3=0 170 TRANS1=.FALSE. C C SET UP AND SOLVE MATRIX EQUATIONS FOR FLOW C if(nit3.gt.100) then jstop = 11 jstp = jstop return end if nit=0 CALL VSMGEN IF (JSTOP.GT.1) THEN JSTP=JSTOP RETURN ENDIF C C CHECK FOR PONDING DURING THIS TIME STEP C CALL VSPOND(IFET,IFET1,IFET2) C C IF PONDING HAS OCCURRED, EQUATIONS NEED TO BE SOLVED AGAIN C IF(IFET.NE.0) THEN IF(NIT.LT.ITMAX) THEN DO 50 II=1,NNODES IF(NTYP(II).NE.1.AND.HX(II).GT.0.0D0) P(II)=PXXX(II) 50 CONTINUE nit = 1 GO TO 170 ELSE WRITE(6,4260) END IF END IF C C REEVALUATE NONLINEAR COEFFICIENTS C CALL VSCOEF NITT=NITT+NIT c END IF c IF((TRANS.OR.VPNT).AND..NOT.SSTATE) CALL VTVELO IF((TRANS.OR.VPNT)) CALL VTVELO IF(TRANS) THEN c IF(.NOT.SSTATE) THEN C C DETERMINE VELOCITIES AND DISPERSION TENSOR C CALL VTDCOEF c END IF TRANS1=.TRUE. C C SET UP AND SOLVE MATRIX EQUATION FOR TRANSPORT C CALL VTSETUP NITT1=NITT1+NIT1 DMAX1=0.0D0 NIT3=NIT3+1 IF(NIT3.LT.2.OR.RHOMAX.GT.EPS2) GO TO 170 DO 190 N=1,NNODES RHOOLD(N)=RHO(N) 190 CONTINUE END IF NITT2=NITT2+NIT3 C C PRINT RESULTS AND COMPUTE MASS BALANCE COMPONENTS C CALL VSOUTP CALL VSFLUX c IF(JSTOP.NE.1) GO TO 160 C C------------------------------------------------------------------- C END OF TIME LOOP C------------------------------------------------------------------- C C *** ASSIGN VALUES TO SUBROUTINE ARGUMENTS TO RETURN JSTP = JSTOP RETURN 4260 FORMAT(5X,'-- WARNING -- INFILTRATION/PONDING BOUNDARY WAS NOT', &' SOLVED ACCURATELY FOR THIS TIME STEP') END SUBROUTINE VSREAD C****** CVSREAD C****** C C PURPOSE: TO READ INITIAL HEAD AND SATURATION DATA C C ---------------------------------------------------------------- C C SPECIFICATIONS FOR ARRAYS AND SCALARS C include 'd_rspac.inc' include 'd_kcon.inc' include 'd_mprop.inc' include 'd_press.inc' include 'd_jtxx.inc' include 'd_dumm.inc' include 'd_ptet.inc' include 'd_trxx.inc' include 'd_idumm.inc' include 'd_rprop.inc' include 'd_scon.inc' include 'd_equat.inc' IMPLICIT DOUBLE PRECISION (A-H,P-Z) c include 'c_rspac.inc' c include 'c_kcon.inc' c include 'c_mprop.inc' c include 'c_press.inc' c include 'c_jtxx.inc' c include 'c_dumm.inc' c include 'c_ptet.inc' c include 'c_trxx.inc' c include 'c_idumm.inc' c include 'c_rprop.inc' c include 'c_scon.inc' c include 'c_equat.inc' COMMON/ISPAC/NLY,NLYY,NXR,NXRR,NNODES COMMON/WGT/WUS,WDS COMMON/JCON/JSTOP,JFLAG,jflag1 LOGICAL TRANS,TRANS1,SORP,SSTATE COMMON/TRXY/EPS1,EPS2,TRANS,TRANS1,SORP,SSTATE,MB9(72),NMB9 LOGICAL PHRD LOGICAL RAD,BCIT,ETSIM,SEEP,ITSTOP,CIS,CIT COMMON/LOG1/RAD,BCIT,ETSIM,SEEP,ITSTOP,CIS,CIT COMMON/hinterp/ DELP,nprop,I1,I2,I3,I4,I5,I6 CHARACTER*80 TITL CHARACTER*80 IFMT,FREE CHARACTER*4 ZUNIT,TUNIT,CUNX COMMON/SCHAR/TITL,ZUNIT,TUNIT,CUNX common/functiontype/ hydraulicFunctionType, adsorptionType COMMON/VERSION/VER1P0 LOGICAL VER1P0 integer hydraulicFunctionType, adsorptionType, iuHead, &iuConcentration common/iuNumber/ iuHead,iuConcentration include 'd_idummAlloc.inc' C----------------------------------------------------------------------- C C READ AND WRITE INITIAL DATA FOR SIMULATION C FREE = 'free' IF (TRANS) THEN READ (5,*) EPS,HMAX,WUS,EPS1,EPS2 ELSE READ(5,*) EPS,HMAX,WUS EPS1=0.0D0 EPS2=0.0D0 END IF READ (5,*) MINIT,ITMAX READ (05,*) PHRD IF(TRANS) THEN READ (05,*) NTEX,NPROP,NPROP1 IF (.NOT.VER1P0) READ (05,*) hydraulicFunctionType,adsorptionType ELSE READ (05,*) NTEX,NPROP IF (.NOT.VER1P0) READ (05,*) hydraulicFunctionType NPROP1=0 END IF include 'd_rpropAlloc.inc' C C CHECK THAT SUM OF WEIGHTING FACTORS IS EQUAL TO ONE C WRITE (6,4000) EPS,ZUNIT,EPS1,HMAX,EPS2 IF(WUS.EQ.1.0D0) THEN WDS=0.0D0 WRITE(06,4020) ELSE IF(WUS.EQ.0.5D0) THEN WDS=0.5D0 WRITE(06,4070) ELSE WUS=0.0D0 WRITE(06,4010) END IF END IF WRITE (6,4080) NTEX,NPROP,NPROP1,MINIT,ITMAX include 'd_sconAlloc.inc' c include 'c_sconItmax.inc' WRITE (06,4100) IF (TRANS) WRITE(06,4110) C C READ AND WRITE MATERIAL PROPERTIES FOR EACH TEXTURAL CLASS C DO 20 J22=1,NTEX DO 10 J23=1,NPROP 10 HK(J22,J23)=0.0D0 DO 20 J23=1,NPROP1 20 HT(J22,J23)=0.0D0 DO 30 J22=1,NTEX READ (5,*) J READ (5,*) ANIZ(J),(HK(J,I),I=1,NPROP) WRITE (6,4120) J,ANIZ(J),(HK(J,I),I=1,NPROP) IF(TRANS) THEN c c modification to read in only nontrivial vs2dh parameters c read(5,*) ht(j,1),ht(j,2),ht(j,5),ht(j,9),ht(j,10),ht(j,11) write(6,4130) ht(j,1),ht(j,2),ht(j,5),ht(j,9),ht(j,10),ht(j,11) ht(j,5)=ht(j,5)*(1.-hk(j,3)) END IF 30 CONTINUE c convert alpha parameter if reading version 2.5 data set IF (VER1P0.AND.(hydraulicFunctionType.EQ.1)) THEN DO 31 J=1,NTEX IF (HK(J,4).NE.0) THEN HK(J,4) = -1.D0/HK(J,4) END IF 31 CONTINUE ENDIF WRITE (06,4140) C C READ TEXTURAL CLASS INDEX MAP C READ (05,*) IROW IF(IROW.EQ.0) THEN WRITE(06,4090) DO 50 J=1,NLY READ (05,*) (IDUM(N),N=1,NXR) WRITE (06,4150) J,(IDUM(N),N=1,NXR) DO 40 N=1,NXR IN=NLY*(N-1)+J J22=IDUM(N) HX(IN)=HK(J22,1) 40 JTEX(IN)=J22 50 CONTINUE ELSE C C READ TEXTURE CLASSES BY BLOCK--EITHER CONTINUOUS LAYERS OR C LAYERS BOUNDED BY VERTICAL DISCONTINUITIES. C WRITE (06,4040) JTP=1 60 READ (05,*) IL,IR,JBT,JRD DO 70 N=IL,IR IDUM(N)=JRD 70 CONTINUE IF(IR.LT.NXR) GO TO 60 DO 80 J=JTP,JBT if (ntex.gt.9) then write (06,4151) j,(idum(n),n=1,nxr) else WRITE (06,4150) J,(IDUM(N),N=1,NXR) end if 80 continue DO 90 J=JTP,JBT DO 90 N=1,NXR IN=NLY*(N-1)+J J22=IDUM(N) HX(IN)=HK(J22,1) JTEX(IN)=J22 90 CONTINUE IF(JBT.EQ.NLY) GO TO 100 JTP=JBT+1 GO TO 60 END IF 100 CONTINUE C C BORDERS OF DOMAIN ARE ALL SET TO NO FLOW BOUNDARIES C DO 110 I=1,NLY I1=NNODES-I+1 HX(I)=0.0D0 110 HX(I1)=0.0D0 DO 120 I=2,NXR I1=(I-1)*NLY HX(I1)=0.0D0 120 HX(I1+1)=0.0D0 C C READ INITIAL HEADS OR MOISTURE CONTENTS C READ (05,*) IREAD,FACTOR IF(IREAD.EQ.2) THEN READ (05,*) DWTX,HMIN WRITE (06,4190) DWTX,ZUNIT,HMIN,ZUNIT,DWTX,ZUNIT C C CALCULATE EQUILIBRIUM INITIAL HEAD PROFILE C DO 130 J=2,NLYY DO 130 N=2,NXRR IN=NLY*(N-1)+J IF(HX(IN).EQ.0.0D0) GO TO 130 IF(CS1.EQ.1.0D0) THEN Z1=DZZ(J) ELSE Z1=DZZ(J)*CS1+RX(N)*CS2 END IF P1=Z1-DWTX IF(P1.LT.HMIN)P1=HMIN P(IN)=P1-Z1 PXXX(IN)=P(IN) 130 CONTINUE ELSE IF(IREAD.EQ.0) THEN WRITE (6,4170) FACTOR ELSE READ (05,*) IU,IFMT WRITE (06,4180) IU,FACTOR END IF DO 160 J=1,NLY IF(IREAD.EQ.1) then IF(IFMT.eq.FREE) THEN C C READ INITIAL CONDITIONS FROM FILE IU C read (iuHead,*) (dum(n),n=1,nxr) ELSE READ (iuHead,FMT=IFMT) (DUM(N),N=1,NXR) end if else DO 140 N=1,NXR 140 DUM(N)=FACTOR end if DO 150 N=1,NXR IN=NLY*(N-1)+J IF(IREAD.ne.0)DUM(N)=DUM(N)*FACTOR IF(CS1.EQ.1.0D0) THEN Z1=DZZ(J) ELSE Z1=DZZ(J)*CS1+RX(N)*CS2 END IF IF(.NOT.PHRD) THEN IF(DUM(N).LE.0.0D0) THEN WRITE(6,4230) J,N jstop=4 return END IF C C CONVERT INITIAL MOISTURE CONTENTS TO HEADS C if(hydraulicFunctionType.eq.1) then P(IN)=VSTHNVVG(DUM(N),JTEX(IN))-Z1 else if(hydraulicFunctionType.eq.2) then P(IN)=VSTHNVHK(DUM(N),JTEX(IN))-Z1 else if(hydraulicFunctionType.eq.0) then P(IN)=VSTHNVBC(DUM(N),JTEX(IN))-Z1 else if(hydraulicFunctionType.eq.3) then P(IN)=VSTHNVTAB(DUM(N),JTEX(IN))-Z1 else P(IN)=VSTHNVOT(DUM(N),JTEX(IN))-Z1 end if end if end if end if THETA(IN)=DUM(N) ELSE P(IN)=DUM(N)-Z1 END IF PXXX(IN)=P(IN) 150 CONTINUE 160 CONTINUE if(jstop.gt.1) return C C COMPUTE INITIAL NONLINEAR COEFFICIENT VALUES C END IF CALL VSCOEF C C IF ET IS TO BE SIMULATED, ALL VARIABLES MUST BE ENTERED HERE. C READ(05,*) BCIT,ETSIM IF(BCIT .OR. ETSIM) THEN C C COMPUTE DEPTHS FOR ET CALCULATIONS C DPTH(1)=-0.5D0 *DELZ(1) DO 170 J=2,NLYY DO 170 N=2,NXRR IN=NLY*(N-1)+J JM1=IN-1 IF(HX(IN).NE.0.0D0) THEN IF(HX(JM1).EQ.0.0D0) THEN DPTH(IN)=0.0D0 ELSE DPTH(IN)=DPTH(JM1)+DELZ(J-1) END IF END IF 170 CONTINUE WRITE (6,4200) CALL VSOUT(2,DPTH) C C READ EVAPORATION VARIABLES C READ(05,*)NPV,ETCYC include 'd_ptetAlloc.inc' DO 713 I=1,NPV PEVAL(I) = 0.0D0 PTVAL(I) = 0.0D0 DO 714 J=1,6 RDC(J,I) = 0.0D0 714 CONTINUE 713 CONTINUE WRITE(6,4030) NPV,ETCYC,TUNIT IF(BCIT) THEN READ (05,*)(PEVAL(I),I=1,NPV) READ(05,*) (RDC(1,I),I=1,NPV) READ(05,*) (RDC(2,I),I=1,NPV) WRITE (06,4050)ZUNIT,TUNIT,ZUNIT,ZUNIT,(I,PEVAL(I),RDC(1,I),RDC(2, *I),I=1,NPV) END IF IF (ETSIM )THEN C C READ TRANSPIRATION VARIABLES C READ(05,*)(PTVAL(I),I=1,NPV) READ(05,*) (RDC(3,I),I=1,NPV) READ(05,*) (RDC(4,I),I=1,NPV) READ(05,*) (RDC(5,I),I=1,NPV) READ(05,*) (RDC(6,I),I=1,NPV) WRITE(06,4060)ZUNIT,TUNIT,ZUNIT,ZUNIT,ZUNIT,ZUNIT,(I,PTVAL(I), *(RDC(J,I),J=3,6),I=1,NPV) END IF END IF DO 180 IN=1,NNODES NTYP(IN)=0 NCTYP(IN)=0 IF(HX(IN).GT.0.0D0) THLST(IN)=THETA(IN) 180 CONTINUE C C READ INITIAL TEMPERATURES IF TRANSPORT EQUATION IS TO C BE SOLVED C IF (TRANS) THEN READ(05,*) IREAD,FACTOR IF(IREAD.EQ.0) THEN WRITE(6,4210) FACTOR DO 190 N=1,NNODES IF(HX(N).GT.0.0D0) THEN CC(N)=FACTOR COLD(N)=FACTOR ELSE CC(N)=0.0D0 COLD(N)=0.0D0 END IF RHO(N)=VTRHO(CC(N),JTEX(N)) RHOOLD(N)=RHO(N) 190 CONTINUE ELSE READ(05,*)IU,IFMT WRITE(06,4220) IU,FACTOR DO 200 J=1,NLY if (IFMT.EQ.FREE) then read(iuConcentration,*) (dum(n),n=1,nxr) else READ(iuConcentration,FMT=IFMT) (DUM(N),N=1,NXR) end if DO 200 N=1,NXR IN=NLY*(N-1)+J IF(HX(IN).GT.0.0D0) THEN CC(IN)=DUM(N)*FACTOR COLD(IN)=CC(IN) ELSE CC(IN)=20.0D0 COLD(IN)=20.0D0 END IF RHO(IN)=VTRHO(CC(IN),JTEX(IN)) RHOOLD(IN)=RHO(IN) 200 CONTINUE end if else do 210 n = 1,nnodes RHO(N)=VTRHO(CC(N),JTEX(N)) RHOOLD(N)=RHO(N) 210 CONTINUE C C COMPUTE INTERCELL CONDUCTANCES C END IF CALL VSHCMP RETURN 4000 FORMAT(10X,27HINITIAL MOISTURE PARAMETERS/10X,27(1H_)// &5X,'CONVERGENCE CRITERION FOR SIP FOR FLOW (EPS) =',1PE12.3,1X,A4/ &5X,'CONVERGENCE CRITERION FOR SIP FOR TRANSPORT (EPS1) =',1PE12.3, &1X,/5X,23HDAMPING FACTOR, HMAX = ,1PE12.3,1X,/5X, &'CONVERGENCE CRITERION FOR OUTER ITERATION (EPS2) = ',1PE12.3) 4010 FORMAT(5X,46HGEOMETRIC MEAN USED FOR INTERCELL CONDUCTIVITY) 4020 FORMAT(5X,45HUPSTREAM WEIGHTING USED FOR INTERCELL CONDUCT &,5HIVITY) 4030 FORMAT(//15X,'NUMBER OF EVAPORATION AND/OR EVAPOTRASPIRATION PER' &,'IODS = ',I4,/,15X,'LENGTH OF EACH PERIOD = ',F10.4,2X,A4) 4040 FORMAT(5X,'TEXTURAL CLASSES READ IN BY BLOCK') 4050 FORMAT(//5X,'EVAPORATION POTENTIAL SURFACE ATMOSHERIC', &/' PERIOD RATE RESISTANCE PRESSURE', &/19X,A4,'/',A4,3X,A4,'**(-1)',5X,A4,/,1X,90('-'), &25(/,5X,I6,4X,3E14.5)) 4060 FORMAT(//,3X,'TRANSPIRATION POTENTIAL ROOT ACTIVIT &Y ACTIVITY ROOT', &/' PERIOD RATE DEPTH AT BOTTOM A &T TOP PRESSURE',/,19X,A4,'/',A4,9X,A4,5X,A4,'**(-2)',4X,A4, &'**(-2)',8X,A4,/,1X,90('-'),25(/,5X,I6,4X,5E14.5)) 4070 FORMAT(5X,47HARITHEMTIC MEAN USED FOR INTERCELL CONDUCTIVITY) 4080 FORMAT(5X,34HNUMBER OF SOIL TEXTURAL CLASSES = ,I10/ &5X,43HNUMBER OF SOIL PARAMETERS FOR EACH CLASS = ,I10/ &5X,'NUMBER OF TRANSPORT PARAMETERS FOR EACH CLASS = ',I10/ &5X,47HMINIMUM PERMITTED NO. OF ITERATIONS/TIME STEP =,I10/ &5X,47HMAXIMUM PERMITTED NO. OF ITERATIONS/TIME STEP =,I10) 4090 FORMAT(5X,41HTEXTURAL CLASS TO BE READ IN FOR EACH ROW) 4100 FORMAT(41X,35HCONSTANTS FOR SOIL TEXTURAL CLASSES// &10X,10HANISOTROPY,7X,4HKSAT,5X,8HSPECIFIC,4X,8HPOROSITY, &4x,'HK(4)',8X,'HK(5)',7X,'HK(6)',/,36X,7HSTORAGE) 4110 FORMAT(12X,'ALPHAL',8X,'ALPHAT',6X,'Cs',9X,'KT MIN', &4X,' KT MAX ',4X,' Cw ') 4120 FORMAT(1X,7HCLASS #,I2,/9X,3(1PD12.3),14(7(1PD12.3),/)) 4130 FORMAT(9X,10(1PD12.3)) 4140 FORMAT(6X,24HTEXTURAL CLASS INDEX MAP// ) 4150 FORMAT(1H ,5X,I4,2X,100I1) 4151 FORMAT(1H ,5X,I4,2X,100I2) 4160 FORMAT(5X,24H ****** VALUE OF ITMAX =,I5,8HEXCEEDS , &44HDIMENSION OF DHMX, PROGRAM TERMINATED ******) 4170 FORMAT(5X,48HINITIAL PRESSURE HEAD OR MOISTURE CONTENT WAS SE, & 24HT TO A CONSTANT VALUE OF,1PE12.3) 4180 FORMAT(5X,48HINITIAL PRESSURE HEAD OR MOISTURE CONTENT WAS RE, & 12HAD FROM UNIT,I5, & 20H A SCALING FACTOR OF,1PE12.3,9H WAS USED) 4190 FORMAT(5X,'EQUILLIBRIUM PROFILE USED TO INITIALIZE PRESSURE', & 27H HEADS ABOVE WATER TABLE AT,F10.2,1X,A4,1X, & 12HBELOW ORIGIN/5X, & 57HEQUILLIBRIUM PROFILE ONLY USED UNTIL PRESSURE HEADS EQUAL, & F10.2,1X,A4/5X, & 20HPRESSURE HEADS BELOW,F10.2,1X,A4,16H ARE HYDROSTATIC) 4200 FORMAT(1H ,50X,18HDEPTH FROM SURFACE) 4210 FORMAT(' INITIAL TEMPERATURE SET TO A CONSTANT VALUE OF ', &1PE12.3,' DEGREES C') 4220 FORMAT(' INITIAL TEMPERATURE WAS READ FROM UNIT',I5, &' A SCALING FACTOR OF, ',1PE12.3,' WAS USED') 4230 FORMAT(' INITIAL MOISTURE CONTENT AT ROW ',I3,' COLUMN ', &I3,' IS LESS THAN OR EQUAL TO 0.',/' PROGRAM TERMINATED') END SUBROUTINE VSTMER C****** CVSTMER C****** C C PURPOSE: TO CONTROL THE TIME SEQUENCE OF SIMULATION C AND TO READ NEW BOUNDARY CONDITION DATA C C ---------------------------------------------------------------- C C SPECIFICATIONS FOR ARRAYS AND SCALARS C include 'd_rspac.inc' include 'd_kcon.inc' include 'd_mprop.inc' include 'd_press.inc' include 'd_disch.inc' include 'd_trxx.inc' include 'd_plott.inc' include 'd_idumm.inc' include 'd_spfc.inc' include 'd_scon.inc' IMPLICIT DOUBLE PRECISION (A-H,P-Z) c include 'c_rspac.inc' c include 'c_kcon.inc' c include 'c_mprop.inc' c include 'c_press.inc' c include 'c_disch.inc' c include 'c_trxx.inc' c include 'c_plott.inc' c include 'c_idumm.inc' c include 'c_spfc.inc' c include 'c_scon.inc' COMMON/ISPAC/NLY,NLYY,NXR,NXRR,NNODES COMMON/PND/POND COMMON/SCN1/TMPX,TMLT,DLTMX,DLTMIN,TRED COMMON/TCON/STIM,DSMAX,KTIM,NIT,NIT1,KP,NIT3 COMMON/JCON/JSTOP,JFLAG,jflag1 LOGICAL TRANS,TRANS1,SORP,SSTATE COMMON/TRXY/EPS1,EPS2,TRANS,TRANS1,SORP,SSTATE,MB9(72),NMB9 LOGICAL RAD,BCIT,ETSIM,SEEP,ITSTOP,CIS,CIT LOGICAL F7P,F11P,F8P,F9P,F6P,PRNT,o9p,o11p,o13p COMMON/LOG1/RAD,BCIT,ETSIM,SEEP,ITSTOP,CIS,CIT COMMON/LOG2/F7P,F11P,F8P,F9P,F6P,PRNT,o9p,o11p,o13p CHARACTER*80 TITL CHARACTER*4 ZUNIT,TUNIT,CUNX COMMON/SCHAR/TITL,ZUNIT,TUNIT,CUNX COMMON/VERSION/VER1P0 LOGICAL VER1P0 SAVE STERR C------------------------------------------------------------------- C C ADVANCE TO NEXT TIME STEP C KTIM=KTIM+1 IF (KTIM.NE.1.AND.JSTOP.NE.0) RETURN JSTOP=0 JPLT=0 NIT=0 NIT1=0 IF(KTIM.EQ.1) KPLT=1 IF(JFLAG.EQ.1) THEN C C ................................................................ C C READ DATA FOR NEW RECHARGE PERIOD C ................................................................ C READ (05,*) TPER,DELT jflag1 = 1 C C CHECK FOR END OF SIMULATION C IF(TPER.LT.0.0D0) THEN WRITE (06,4070) TMAX,STIM jstop=5 return END IF READ (05,*) TMLT,DLTMX,DLTMIN,TRED KP=KP+1 SSTATE=.FALSE. if(delt.lt.dltmin) then delt = dltmin else if(delt.gt.dltmx) delt = dltmx end if WRITE (06,4000) KP,TPER,TUNIT,DELT,TUNIT,TMLT,DLTMX,TUNIT,DLTMIN, *TUNIT,TRED READ (05,*) DSMAX,STERR READ (05,*) POND WRITE (06,4020) DSMAX,ZUNIT,STERR,ZUNIT,POND,ZUNIT READ (05,*) PRNT READ (05,*) BCIT,ETSIM,SEEP WRITE (06,4010) PRNT,BCIT,ETSIM,SEEP DSMAX=DABS(DSMAX) ETOUT=0.0D0 ETOUT1=0.0D0 C C READ SEEPAGE FACE DATA C IF(SEEP) THEN READ (05,*) NFCS include 'd_spfcAlloc.inc' DO 20 K=1,NFCS READ (05,*) JJ,JLAST(K) NFC(K)=JJ READ (05,*) ((JSPX(L,J,K),L=2,3),J=1,JJ) DO 10 J=1,JJ J1=JSPX(2,J,K) N1=JSPX(3,J,K) N2=NLY*(N1-1)+J1 JSPX(1,J,K)=N2 Q(N2)=0.0D0 QQ(N2)=0.0D0 NCTYP(N2)=0 CS(N2)=0.0D0 IF(J.GT.JLAST(K)) THEN NTYP(N2)=3 ELSE NTYP(N2)=1 IF(CS1.EQ.1.0D0) THEN Z1=DZZ(J1) ELSE Z1=DZZ(J1)*CS1+RX(N1)*CS2 END IF P(N2)=-Z1 END IF 10 CONTINUE 20 CONTINUE END IF C C READ IN NEW BOUNDARY CONDITIONS FOR RECHARGE PERIOD C IF IBC=0, POINT BOUNDARY CONDITIONS ARE READ IN. C IF IBC=1, LINE BOUNDARY CONDITIONS ARE READ IN, AND IT IS NECESSARY C TO SPECIFY FOUR POINTS ON THE LINE--THIS ALLOWS VERTICAL OR HORI- C ZONTAL LINES TO BE READ IN INDISCRIMINATELY. THE SEQUENCE IS: C TOP ROW, BOTTOM ROW, LEFT COLUMN, RIGHT COLUMN, CODE, AND FLUX OR C PRESSURE HEAD FOR BOUNDARY CONDITION. C READ (05,*) IBC IF(IBC.GT.0) GO TO 40 30 IF (TRANS) THEN READ(05,*) JJ,NN,NTX,PFDUM,NTC,CF ELSE READ (05,*) JJ,NN,NTX,PFDUM CF=0.0D0 NTC=0 END IF IF(JJ.LT.0) GO TO 90 JJT=JJ JJB=JJ NNL=NN NNR=NN GO TO 50 40 IF (TRANS) THEN READ(05,*) JJT,JJB,NNL,NNR,NTX,PFDUM,NTC,CF ELSE READ (05,*) JJT,JJB,NNL,NNR,NTX,PFDUM CF=0.0D0 NTC=0 END IF IF(JJT.LT.0) GO TO 90 50 CONTINUE DO 80 JJ=JJT,JJB DO 80 NN=NNL,NNR IN=NLY*(NN-1)+JJ CS(IN)=CF c c change made 99-12-16 so that ntc = 2 conductive heat flux is c in terms of energy per time per length of boundary (i.e., c a line source) , as c opposed to energy per time (which is for point source) c if (ntc.eq.2) then area = dxr(nn) if(rad) area = pi2*rx(nn)*dxr(nn) cs(in) = cf*area end if c IF(NTC.EQ.1) CC(IN)=CF NCTYP(IN)=NTC c if(ntx.eq.5.and.hx(in-1).gt.0.0d0) go to 80 IF(NTX.NE.6) GO TO 60 NTYP(IN)=2 QQ(IN)=PFDUM GO TO 80 60 NTYP(IN)=NTX IF(NTX .EQ. 4)NTYP(IN)=1 IF(NTX.EQ.0) WRITE (06,4030) JJ,NN IF(CS1.EQ.1.0D0) THEN Z1=DZZ(JJ) ELSE Z1=DZZ(JJ)*CS1+RX(NN)*CS2 END IF IF(NTX.EQ.1) P(IN)=PFDUM-Z1 IF(NTX.EQ.4) P(IN)=PFDUM IF(NTX.EQ.2) GO TO 70 QQ(IN)=0.0D0 GO TO 80 70 CONTINUE C C SET QQ TO RAINFALL RATE C AREA=DXR(NN) IF(RAD)AREA=PI2*RX(NN)*DXR(NN) QQ(IN)=PFDUM*AREA 80 CONTINUE IF(IBC.EQ.0) GO TO 30 GO TO 40 90 CONTINUE C C WRITE INITIAL BOUNDARY CONDITIONS FOR THIS PERIOD C WRITE (06,4040) KP DO 110 J=1,NLY DO 100 N=1,NXR IN=NLY*(N-1)+J Q(IN)=0.0D0 100 IDUM(N)=NTYP(IN) 110 WRITE (06,4050) J,(IDUM(I),I=1,NXR) TMPX=STIM+TPER IF(TMPX+0.5D0*DLTMIN.GT.TMAX) TMPX=TMAX C C CALCULATE NEW COEFFICIENTS C IF(KTIM.NE.1)CALL VSCOEF END IF C C INITIALIZE REQUIRED ARRAYS FOR NEW BOUNDARY CONDITION, UPDATE C PXXX,THLST. COMPUTE MAXIMUM HEAD CHANGE DURING LAST TIME STEP C PDIF=0.0D0 IF(.NOT.SSTATE) THEN DO 120 J=2,NLYY DO 120 N=2,NXRR IN=NLY*(N-1)+J IF(HX(IN).EQ.0.0D0) GO TO 120 P12=P(IN)-PXXX(IN) PTMP=DABS(P12) IF(PTMP.GT.PDIF)PDIF=PTMP PXXX(IN)=P(IN) THLST(IN)=THETA(IN) IF(TRANS) COLD(IN)=CC(IN) 120 CONTINUE C C CHECK FOR STEADY STATE C IF(jflag.ne.1.and.PDIF.LE.STERR) THEN c SSTATE=.TRUE. c WRITE(6,4060) STIM,KTIM END IF END IF JFLAG=0 C C INITIALIZE DHMX C DO 130 K=1,itmax 130 DHMX(K)=0.0D0 C C ADVANCE DELT AND RESET TO PROPER LENGTH IF NECESSARY C if (delt.lt.dltmin) delt = dltmin DLTOLD=DELT DELT= TMLT*DELT C C MAXIMUM PERMISSABLE HEAD CHANGE CHECK C IF(KTIM.GE.2) THEN IF((PDIF*DELT/DLTOLD).GT.DSMAX)DELT=DLTOLD*DSMAX*0.98D0/PDIF END IF c IF(DABS(TMPX-PLTIM(KPLT)).LT.DLTMIN) PLTIM(KPLT)=TMPX T1=DMIN1(TMPX,PLTIM(KPLT)) T2=T1-STIM c c changes made 6-12-02 for calculating delt. This new c version will allow delt