!$Author: fredfox $ !$Date: 2006-08-21 21:12:09 $ !$Revision: 1.12 $ !$Source: /weru/cvs/wepp/wepp.watbal/grna.for,v $ SUBROUTINE GRNA( NF, DEPSTO, TR, R, RR, KS, SM, & & NS, TF, RCUM, F, FF, RE, RECUM, TP, & & RPRINT, DDEPSTO, RUNOFF, DUREXR, EFFINT, EFFDRR ) ! + + + PURPOSE + + + ! THIS PROGRAM CALCULATES INFILTRATION RATES AND DEPTHS FOR ! UNSTEADY RAIN USING THE GREEN AND AMPT INFILTRATION EQUATION ! AS MODIFIED BY MEIN AND LARSON. THE EQUATION HAS THE FORM : ! F = KS + KS*SM/FF ! AND ! KS*T = FF - SM*LN(1 + FF/SM) ! WHERE: ! F = INFILTRATION RATE (L/T) ! KS = SATURATED CONDUCTIVITY (L/T) ! SM = EFFECTIVE MATRIC POTENTIAL (L) ! FF = CUMULATIVE INFILTRATION (L) ! LN = NATURAL LOG ! INPUT: ! 1. TWO COLUMN FILE ! COL 1 = TIME (MINUTES) ! COL 2 = RAINFALL RATE (MILLIMETERS/HOUR) ! ON THE VAX THE FILE IS AN EDR FILE ! ON THE PC THE FILE HAS A FORMAT 2F10.2 ! 2. SATURATED CONDUCTIVITY (MILLIMETERS/HOUR) ! 3. EFFECTIVE MATRIC POTENTIAL (S*M) (MILLIMETERS) ! WHERE: ! S = DIFFERENCE IN AVERAGE CAPILARY POTENTIAL ! BEFORE AND AFTER WETTING ! M = DIFFERENCE IN AVERAGE SOIL MOISTURE !----------------------------------------------------------- ! PARAMETERS ! 1. R(I) ..............RAINFALL RATE ! 2. RCUM(I) ...........ACCUMULATED RAINFALL DEPTH ! 3. F(I) ..............INFILTRATION RATE ! 4. FF(I) .............ACCUMULATED INFILTRATION DEPTH ! 5. RE(I) .............RAINFALL EXCESS RATE ! 6. RECUM(I) ..........ACCUMULATED RAINFALL EXCESS DEPTH ! 7. CU ................INDICATOR OF PONDING IF NO PONDING AT BEGINNING ! CU < 0 - NO PONDING CU > 0 - PONDING ! 8. CP ................INDICATOR OF PONDING WHEN PONDED AT BEGINNING ! CP < 0 - PONDING STOPS DURING INTERVAL CP > 0 - PONDING ! 9. TR(I) .............RAINFALL TIMES ! 10. TP ...............TIME OF PONDING ! 11. TS ...............PSEUDOTIME TO ADJUST REAL TIME FOR INFILTRATION ! 12. T ................REAL TIME = TR(I)-TP+TS ! 13. PT ...............ACCUMULATED RAINFALL AT TIME OF PONDING !----------------------------------------------------------- ! CALLED FROM MAIN ! AUTHOR(S): D. FLANAGAN, J. ASCOUGH ! VERSION: THIS MODULE TAKEN FROM ASCOUGH STANDALONE IRS CODE ! DATE CODED: 3-22-2005 ! CODED BY: D. FLANAGAN ! + + + PARAMETER DECLARATIONS + + + INTEGER MXTIME, MXPOND PARAMETER (MXTIME = 1000, MXPOND = 65) ! + + + ARGUMENT DECLARATIONS + + + INTEGER, intent(in) :: NF REAL, intent(in) :: DEPSTO, TR(MXTIME), R(MXTIME), RR(MXTIME), & & KS, SM INTEGER, intent(inout) :: NS REAL, intent(inout) :: TF(MXTIME), RCUM(MXTIME), & & F(MXTIME), FF(MXTIME), RE(MXTIME), RECUM(MXTIME), TP(MXPOND),& & RPRINT(MXTIME), DDEPSTO(MXTIME), & & RUNOFF, DUREXR, EFFINT, EFFDRR ! + + + ARGUMENT DEFINITIONS + + + ! DEPSTO - Depression Storage (L) ! TR - time values in rainfall breakpoint array (T) ! TF - time values in infiltration/rainfall excess/depression storage array (T) ! R - rainfall rate (L/T) ! RCUM - accumulated rainfall depth (L) ! F - infiltration rate (L/T) ! FF - accumulated infiltration depth (L) ! RE - rainfall excess rate (L/T) ! RECUM - accumulated rainfall excess depth (L) ! RR - precalculated rainfall rate accumulation array (L) ! NF - number of values in infiltration arrays ! KS - saturated conductivity (L/T) ! SM - effective matric potential (L) ! TP - time of ponding (T) ! NS - number of element for end of runoff arrays ! RPRINT - ! DDEPSTO - depth of depression storage at infiltration timestep (L) ! RUNOFF - RUNOFF DEPTH (L) ! DUREXR - DURATION OF RAINFALL EXCESS (T) ! EFFINT - EFFECTIVE RAINFALL INTENSITY ! EFFDRR - EFFECTIVE RAINFALL DURATION ! + + + LOCAL VARIABLES + + + INTEGER I, K, POND, IT, NP, NT REAL RTEMP(MXTIME), DTIME REAL DURRE, SUMINT, KSM, TT, TS, CU, CP, FSUM REAL PT(MXPOND), PRECUM(MXPOND) DOUBLE PRECISION XX ! + + + LOCAL DEFINITIONS + + + ! DTIME - time step ! DURRE - duration of effective rainfall ! SUMINT - rainfall intensity sum ! KSM - product of KS and SM ! TT - temporary "real" time variable ! TS - pseudotime to adjust real time for infiltration ! CU - indicator of ponding if no ponding at beginning ! CU < 0 - no ponding CU > 0 - ponding ! CP - indicator of ponding when ponded at beginning ! CP < 0 - ponding stops during interval CP > 0 - ponding ! XX - local temporary variable (double precision) ! FSUM - infiltration sum value ! PT - accumulated rainfall at time of ponding ! PRECUM - accumulated rainfall excess at time of ponding ! + + + END SPECIFICATIONS + + + KSM = KS * SM POND = 0 FSUM = 0.0 K = 1 FF(1) = 0.0 RECUM(1) = 0.0 RCUM(1) = 0.0 RTEMP(1) = R(1) DDEPSTO(1) = 0.0 IT = 0 NP = 0 DO 10 I = 1, MXPOND TP(I) = 0. 10 CONTINUE ! START RUN DO 20 I = 2, NF K = K + 1 ! CHECK IF THERE IS PONDING IN PREVIOUS INTERVAL IF (POND.EQ.0) THEN ! CASE ONE: NO PONDING IN PREVIOUS INTERVAL IF (R(I-1).LE.KS) THEN IF (R(I-1).EQ.0.0) THEN FF(K) = FSUM RECUM(K) = RECUM(K-1) RCUM(K) = RR(I) RTEMP(K) = R(I) TF(K) = TR(I) POND = 0 DDEPSTO(K) = 0.0 GO TO 20 ELSE FSUM = RR(I) - RECUM(K-1) FF(K) = FSUM RECUM(K) = RECUM(K-1) RCUM(K) = RR(I) RTEMP(K) = R(I) TF(K) = TR(I) POND = 0 DDEPSTO(K) = 0.0 GO TO 20 END IF END IF ! PONDING INDICATOR WHEN NO PONDING IN PREVIOUS INTERVAL CU = RR(I) - RECUM(K-1) - KSM / (R(I-1)-KS) ! CASE ONE-A: NO PONDING IF (CU.LE.0.0) THEN FSUM = RR(I) - RECUM(K-1) FF(K) = FSUM RECUM(K) = RECUM(K-1) RCUM(K) = RR(I) RTEMP(K) = R(I) TF(K) = TR(I) POND = 0 DDEPSTO(K) = 0.0 GO TO 20 END IF ! CASE ONE-B: PONDING - GET TIME TO PONDING, TP POND = 1 IT = IT + 1 TP(IT) = (KSM/(R(I-1)-KS)-RR(I-1)+RECUM(K-1))/R(I-1)+TR(I-1) IF (TP(IT).LE.TR(I-1)) THEN ! ponding starts at beginning of present interval TP(IT) = TR(I-1) PT(IT) = RR(I-1) ELSE ! ponding starts part way through the interval ! insert new time point FF(K) = FF(K-1) + R(I-1) * (TP(IT)-TR(I-1)) RECUM(K) = RECUM(K-1) PT(IT) = RR(I-1) + (TP(IT)-TR(I-1)) * R(I-1) RCUM(K) = PT(IT) RTEMP(K) = R(I) TF(K) = TP(IT) DDEPSTO(K) = 0.0 K = K + 1 END IF ! save cumulative rainfall excess at this ponding point PRECUM(IT) = RECUM(K-1) ! CUMULATIVE RAINFALL, PT, AT TIME TO PONDING, TP IF( SM .GT. 0.0 ) THEN ! infiltration with matric potential ! PSEUDOTIME - GET TIME SHIFT DUE TO INFILTRATION, TS XX = (PT(IT)-RECUM(K-1)) / SM TS = SM / KS * (XX-DLOG(1.D0+XX)) ! REAL TIME, T TT = TR(I) - TP(IT) + TS ! CUMULATIVE INFILTRATION, NEWTONS METHOD ALA LJL CALL NEWTON(TT, FF(K-1), FF(K), KS, SM) ELSE ! infiltration with gravitational potential only FF(K) = PT(IT) - PRECUM(IT) + KS * (TR(I) - TP(IT)) END IF DDEPSTO(K) = RR(I) - FF(K) IF( DDEPSTO(K) .GT. DEPSTO ) THEN RECUM(K) = DDEPSTO(K) - DEPSTO DDEPSTO(K) = DEPSTO ELSE RECUM(K) = RECUM(K-1) END IF RCUM(K) = RR(I) RTEMP(K) = R(I) TF(K) = TR(I) ELSE ! CASE TWO: PONDING IN PREVIOUS INTERVAL IF( SM .GT. 0.0 ) THEN ! infiltration with matric potential TT = TR(I) - TP(IT) + TS CALL NEWTON(TT, FF(K-1), FF(K), KS, SM) ELSE ! infiltration with gravitational potential only FF(K) = PT(IT) - PRECUM(IT) + KS * (TR(I) - TP(IT)) END IF ! CHECK IF NO PONDING BEFORE END OF INTERVAL ! with ponding, depression storage adds to available water CP = RR(I) - FF(K) - RECUM(K-1) ! CASE TWO-A: NO PONDING BEFORE END OF INTERVAL IF (CP.LT.0.0) THEN FF(K) = RR(I) - RECUM(K-1) FSUM = FF(K) DDEPSTO(K) = 0.0 RECUM(K) = RECUM(K-1) RCUM(K) = RR(I) RTEMP(K) = R(I) TF(K) = TR(I) POND = 0 ELSE ! CASE TWO-B: PONDING CONTINUES MERRILY ON DDEPSTO(K) = RR(I) - FF(K)- RECUM(K-1) IF( DDEPSTO(K) .GT. DEPSTO ) THEN RECUM(K) = DDEPSTO(K) - DEPSTO + RECUM(K-1) DDEPSTO(K) = DEPSTO ELSE RECUM(K) = RECUM(K-1) END IF RCUM(K) = RR(I) RTEMP(K) = R(I) TF(K) = TR(I) END IF END IF 20 CONTINUE ! check for remaining depression storage IF( K .GT. 2 ) THEN DTIME = TF(K-1) - TF(K-2) DO WHILE( DDEPSTO(K) .GT. 0.0 ) K = K + 1 TF(K) = TF(K-1) + DTIME ! CASE TWO: PONDING IN PREVIOUS INTERVAL IF( SM .GT. 0.0 ) THEN ! infiltration with matric potential TT = TF(K) - TP(IT) + TS CALL NEWTON(TT, FF(K-1), FF(K), KS, SM) ELSE ! infiltration with gravitational potential only FF(K) = PT(IT) - PRECUM(IT) + KS * (TF(K) - TP(IT)) END IF ! CHECK IF NO PONDING BEFORE END OF INTERVAL ! with ponding, depression storage adds to available water CP = RR(NF) - FF(K) - RECUM(K-1) ! CASE TWO-A: NO PONDING BEFORE END OF INTERVAL IF (CP.LT.0.0) THEN FF(K) = RR(NF) - RECUM(K-1) DDEPSTO(K) = 0.0 ELSE ! CASE TWO-B: PONDING CONTINUES MERRILY ON DDEPSTO(K) = RR(NF) - FF(K) - RECUM(K-1) END IF RECUM(K) = RECUM(K-1) RCUM(K) = RCUM(K-1) RTEMP(K) = 0.0 END DO END IF NS = K NT = 0 DO 30 I = 1, NS - 1 DTIME = TF(I+1) - TF(I) F(I) = (FF(I+1)-FF(I)) / DTIME IF (F(I).LT.0.0) F(I) = 0.0 RPRINT(I) = RTEMP(I) RE(I) = (RECUM(I+1)-RECUM(I)) / DTIME IF (RE(I).LT.0.0) RE(I) = 0.0 IF (RE(I).GT.0.0) THEN NT = I NP = NP + 1 END IF 30 CONTINUE F(NS) = 0.0 RE(NS) = 0.0 RUNOFF = RECUM(NS) DUREXR = MAX( 0.0, TF(NT+1) - TP(1) ) FF(NS) = RCUM(NS) - RECUM(NS) ! GET DURATION OF RAINFALL EXCESS AND EFFECTIVE RAINFALL INTENSITY DURRE = 0.0 SUMINT = 0.0 DO I = 1, NS-1 IF (RE(I).GT.0.0) THEN DURRE = DURRE + TF(I+1) - TF(I) SUMINT = SUMINT + (TF(I+1)-TF(I)) * RTEMP(I) ENDIF END DO IF (DURRE.GT.0.0) THEN EFFINT = SUMINT/DURRE EFFDRR = DURRE ELSE EFFINT = 0.0 EFFDRR = 0.0 ENDIF RETURN END