From 6dbc8c7a2bc155fb44ba60d95e757d4179241981 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 5 Feb 2022 17:38:21 -0500 Subject: [PATCH 1/5] bug fix: units of RUNSRF; cleanup: units of THRU[x] (catchment.F90, lsm_routines.F90) --- .../GEOScatch_GridComp/catchment.F90 | 4 +++- .../GEOSland_GridComp/Shared/lsm_routines.F90 | 19 ++++++++++++++++--- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 index 6c333c939..aff9f204b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 @@ -1772,7 +1772,9 @@ SUBROUTINE RZDRAIN ( & ENDIF IF(CATDEF(N) .LT. 0.) THEN - RUNSRF(N)=RUNSRF(N)-CATDEF(N) + ! bug fix: RUNSRF in flux units (kg m-2 s-1) for consistency with partition() + ! reichle, 5 Feb 2022 + RUNSRF(N)=RUNSRF(N)-CATDEF(N)/DTSTEP CATDEF(N)=0. ENDIF diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 index f3cf0af6f..4e934fa94 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 @@ -199,7 +199,12 @@ SUBROUTINE INTERC ( & THRUL(CHNO)=AMAX1(0., THRUL(CHNO)) THRUC(CHNO)=AMAX1(0., THRUC(CHNO)) - 100 CONTINUE + 100 CONTINUE ! ends do loop??? +!**** + ! return throughfall as a flux (kg m-2 s-1) for consistency with units of precip inputs + ! - reichle, 5 Feb 2022 + THRUL = THRUL/DTSTEP + THRUC = THRUC/DTSTEP !**** RETURN END SUBROUTINE INTERC @@ -213,7 +218,7 @@ END SUBROUTINE INTERC !**** =================================================== SUBROUTINE SRUNOFF ( & - NCH,DTSTEP,UFW4RO, FWETC, FWETL, AR1,ar2,ar4, THRUL,THRUC, & + NCH,DTSTEP,UFW4RO, FWETC, FWETL, AR1,ar2,ar4, THRUL_IN,THRUC_IN, & frice,tp1,srfmx, BUG, & SRFEXC,RUNSRF, & QINFIL & @@ -226,7 +231,7 @@ SUBROUTINE SRUNOFF ( & REAL, INTENT(IN) :: DTSTEP, FWETC, FWETL LOGICAL, INTENT (IN):: UFW4RO REAL, INTENT(IN), DIMENSION(NCH) :: AR1, ar2, ar4, frice, tp1, & - srfmx, THRUL, THRUC + srfmx, THRUL_IN, THRUC_IN LOGICAL, INTENT(IN) :: BUG REAL, INTENT(INOUT), DIMENSION(NCH) :: SRFEXC ,RUNSRF @@ -235,9 +240,16 @@ SUBROUTINE SRUNOFF ( & INTEGER N REAL deficit,srun0,frun,qin, qinfil_l, qinfil_c, qcapac, excess_infil, srunc, srunl, ptotal + REAL, DIMENSION(NCH) :: THRUL, THRUC !**** - - - - - - - - - - - - - - - - - - - - - - - - - + ! calculations throughout srunoff() are in "volumes" (kg m-2) + ! convert input fluxes to volumes, reichle, 5 Feb 2022 + THRUL = THRUL_IN * DTSTEP + THRUC = THRUC_IN * DTSTEP + RUNSRF = RUNSRF * DTSTEP + DO N=1,NCH if(.not.UFW4RO) then @@ -299,6 +311,7 @@ SUBROUTINE SRUNOFF ( & endif SRFEXC(N)=SRFEXC(N)+QIN + ! convert outputs back to flux units (kg m-2 s-1) RUNSRF(N)=RUNSRF(N)/DTSTEP QINFIL(N)=QIN/DTSTEP From ebaa89884b5ff58bd6bdceb64caabec3b493a8a1 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 5 Feb 2022 18:53:34 -0500 Subject: [PATCH 2/5] undoing unit conversions for THRU[x] (hoping to avoid non-zero-diff result from previous commit) --- .../GEOSland_GridComp/Shared/lsm_routines.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 index 4e934fa94..c72af4979 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 @@ -203,8 +203,8 @@ SUBROUTINE INTERC ( & !**** ! return throughfall as a flux (kg m-2 s-1) for consistency with units of precip inputs ! - reichle, 5 Feb 2022 - THRUL = THRUL/DTSTEP - THRUC = THRUC/DTSTEP + !!THRUL = THRUL/DTSTEP + !!THRUC = THRUC/DTSTEP !**** RETURN END SUBROUTINE INTERC @@ -246,8 +246,8 @@ SUBROUTINE SRUNOFF ( & ! calculations throughout srunoff() are in "volumes" (kg m-2) ! convert input fluxes to volumes, reichle, 5 Feb 2022 - THRUL = THRUL_IN * DTSTEP - THRUC = THRUC_IN * DTSTEP + !!THRUL = THRUL_IN * DTSTEP + !!THRUC = THRUC_IN * DTSTEP RUNSRF = RUNSRF * DTSTEP DO N=1,NCH From 30a6556c1ab489241ce8c15236383cc0cf72d2b6 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sun, 6 Feb 2022 09:32:41 -0500 Subject: [PATCH 3/5] fixed bug in previous commit --- .../GEOSland_GridComp/Shared/lsm_routines.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 index c72af4979..1a43bf00c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 @@ -246,8 +246,8 @@ SUBROUTINE SRUNOFF ( & ! calculations throughout srunoff() are in "volumes" (kg m-2) ! convert input fluxes to volumes, reichle, 5 Feb 2022 - !!THRUL = THRUL_IN * DTSTEP - !!THRUC = THRUC_IN * DTSTEP + THRUL = THRUL_IN !!* DTSTEP + THRUC = THRUC_IN !!* DTSTEP RUNSRF = RUNSRF * DTSTEP DO N=1,NCH From 0ec04eafe8bf164b3757dfbf4d99184f518c97ca Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sun, 6 Feb 2022 13:07:40 -0500 Subject: [PATCH 4/5] additional comments and cleanup: - clarified units of surface runoff, infiltration, and throughfall in subroutines interc(), rzdrain(), srunof() - moved subroutine rzdrain() to lsm_routines.F90 (identical for Catchment and CatchmentCN) --- .../Shared/catchmentCN.F90 | 135 +---------- .../GEOScatch_GridComp/catchment.F90 | 152 +----------- .../GEOSland_GridComp/Shared/lsm_routines.F90 | 226 +++++++++++++++--- 3 files changed, 200 insertions(+), 313 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchmentCN.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchmentCN.F90 index 3078b2814..5d27e5847 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchmentCN.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchmentCN.F90 @@ -101,12 +101,12 @@ MODULE CATCHMENT_CN_MODEL USE SURFPARAMS, ONLY: CSOIL_2, RSWILT, & LAND_FIX, FLWALPHA - USE lsm_routines, only : & - INTERC, BASE, PARTITION, RZEQUIL, gndtp0, & - catch_calc_soil_moist, gndtmp, & + INTERC, RZDRAIN, BASE, PARTITION, RZEQUIL,& + gndtp0, gndtmp, & + catch_calc_soil_moist, & catch_calc_wtotl, dampen_tc_oscillations, & - SRUNOFF + SRUNOFF USE SIBALB_COEFF, ONLY: coeffsib @@ -1642,132 +1642,7 @@ END SUBROUTINE CATCHCN !**** =================================================== !**** /////////////////////////////////////////////////// !**** =================================================== - - SUBROUTINE RZDRAIN ( & - NCH,DTSTEP,VGWMAX,SATCAP,RZEQ,AR1,WPWET, & - tsa1,tsa2,tsb1,tsb2,atau,btau,CDCR2,poros,BUG, & - CAPAC,RZEXC,SRFEXC,CATDEF,RUNSRF & - ) - -!----------------------------------------------------------------- -! defines drainage timescales: -! - tsc0, between srfex and rzex -! - tsc2, between rzex and catdef -! then defines correponding drainages -! and updates the water contents -!----------------------------------------------------------------- - - IMPLICIT NONE - INTEGER, INTENT(IN) :: NCH - REAL, INTENT(IN) :: DTSTEP - REAL, INTENT(IN), DIMENSION(NCH) :: VGWMAX, SATCAP, RZEQ, AR1, wpwet, & - tsa1, tsa2, tsb1, tsb2, atau, btau, CDCR2, poros - LOGICAL, INTENT(IN) :: BUG - - REAL, INTENT(INOUT), DIMENSION(NCH) :: RZEXC, SRFEXC, CATDEF, CAPAC, & - RUNSRF - - - INTEGER N - REAL srflw,rzflw,FLOW,EXCESS,TSC0,tsc2,rzave,rz0,wanom,rztot, & - rzx,btaux,ax,bx,rzdif - - -!**** - - - - - - - - - - - - - - - - - - - - - - - - - - - DO 100 N=1,NCH - -!**** Compute equivalent of root zone excess in non-saturated area: - rztot=rzeq(n)+rzexc(n) - if(ar1(n).ne.1.) then - !!! rzave=(rztot-ar1(n)*vgwmax(n))/(1.-ar1(n)) - !!! rzave=rzave*poros(n)/vgwmax(n) - rzave=rztot*poros(n)/vgwmax(n) - else - rzave=poros(n) - endif - -! updated warning statement, reichle+koster, 12 Aug 2014 -! -! Impose minimum of 1.e-4, rather than leaving positive values <1.e-4 unchanged. -! -reichle, 15 Jan 2016 - if (rzave .le. 1.e-4) then - rzave=1.e-4 - print*,'problem: rzave <= 1.e-4 in catchment',n - end if - - btaux=btau(n) - if (srfexc(n) .lt. 0.) btaux=btau(n)*(poros(n)/rzave) - rz0=amax1(0.001,rzave-srfexc(n)/(1000.*(-btaux))) - tsc0=atau(n)/(rz0**3.) - - tsc0=tsc0*3600. - if(tsc0.lt.dtstep) tsc0=dtstep - -! --------------------------------------------------------------------- - - SRFLW=SRFEXC(N)*DTSTEP/TSC0 - IF(SRFLW < 0. ) SRFLW = FLWALPHA * SRFLW! C05 change -!rr following inserted by koster Sep 22, 2003 - rzdif=rzave/poros(n)-wpwet(n) -!**** No moisture transport up if rz at wilting; employ ramping. - if(rzdif.le.0. .and. srflw.lt.0.) srflw=0. - if(rzdif.gt.0. .and. rzdif.lt.0.01 & - .and. srflw.lt.0.) srflw=srflw*(rzdif/0.01) - RZEXC(N)=RZEXC(N)+SRFLW - SRFEXC(N)=SRFEXC(N)-SRFLW - -!**** Topography-dependent tsc2, between rzex and catdef - - rzx=rzexc(n)/vgwmax(n) - - if(rzx .gt. .01) then - ax=tsa1(n) - bx=tsb1(n) - elseif(rzx .lt. -.01) then - ax=tsa2(n) - bx=tsb2(n) - else - ax=tsa2(n)+(rzx+.01)*(tsa1(n)-tsa2(n))/.02 - bx=tsb2(n)+(rzx+.01)*(tsb1(n)-tsb2(n))/.02 - endif - - tsc2=exp(ax+bx*catdef(n)) - rzflw=rzexc(n)*tsc2*dtstep/3600. - - IF (CATDEF(N)-RZFLW .GT. CDCR2(N)) then - RZFLW=CATDEF(N)-CDCR2(N) - end if - - CATDEF(N)=CATDEF(N)-RZFLW - RZEXC(N)=RZEXC(N)-RZFLW - -!**** REMOVE ANY EXCESS FROM MOISTURE RESERVOIRS: - - IF(CAPAC(N) .GT. SATCAP(N)) THEN - RZEXC(N)=RZEXC(N)+CAPAC(N)-SATCAP(N) - CAPAC(N)=SATCAP(N) - ENDIF - - IF(RZEQ(N) + RZEXC(N) .GT. VGWMAX(N)) THEN - EXCESS=RZEQ(N)+RZEXC(N)-VGWMAX(N) - RZEXC(N)=VGWMAX(N)-RZEQ(N) - CATDEF(N)=CATDEF(N)-EXCESS - ENDIF - - IF(CATDEF(N) .LT. 0.) THEN - RUNSRF(N)=RUNSRF(N)-CATDEF(N) - CATDEF(N)=0. - ENDIF - - 100 ENDDO - - RETURN - END SUBROUTINE RZDRAIN - -!**** ----------------------------------------------------------------- -!**** ///////////////////////////////////////////////////////////////// -!**** ----------------------------------------------------------------- +!**** !**** [ BEGIN FLUXES ] !**** SUBROUTINE FLUXES ( & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 index aff9f204b..ba710c9d5 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90 @@ -98,12 +98,13 @@ MODULE CATCHMENT_MODEL LAND_FIX, ASTRFR, STEXP, RSWILT, & FLWALPHA, CSOIL_2 - USE lsm_routines, only : & - INTERC, BASE, PARTITION, RZEQUIL, gndtp0, & - catch_calc_soil_moist, gndtmp, & - catch_calc_wtotl, dampen_tc_oscillations, & - SRUNOFF - + USE lsm_routines, only : & + INTERC, RZDRAIN, BASE, PARTITION, RZEQUIL,& + gndtp0, gndtmp, & + catch_calc_soil_moist, & + catch_calc_wtotl, dampen_tc_oscillations, & + SRUNOFF + USE SIBALB_COEFF, ONLY: coeffsib USE STIEGLITZSNOW, ONLY: & @@ -1649,145 +1650,6 @@ END SUBROUTINE CATCHMENT !**** =================================================== !**** /////////////////////////////////////////////////// !**** =================================================== - - SUBROUTINE RZDRAIN ( & - NCH,DTSTEP,VGWMAX,SATCAP,RZEQ,AR1,WPWET, & - tsa1,tsa2,tsb1,tsb2,atau,btau,CDCR2,poros,BUG, & - CAPAC,RZEXC,SRFEXC,CATDEF,RUNSRF & - ) - -!----------------------------------------------------------------- -! defines drainage timescales: -! - tsc0, between srfex and rzex -! - tsc2, between rzex and catdef -! then defines correponding drainages -! and updates the water contents -!----------------------------------------------------------------- - - IMPLICIT NONE - INTEGER, INTENT(IN) :: NCH - REAL, INTENT(IN) :: DTSTEP - REAL, INTENT(IN), DIMENSION(NCH) :: VGWMAX, SATCAP, RZEQ, AR1, wpwet, & - tsa1, tsa2, tsb1, tsb2, atau, btau, CDCR2, poros - LOGICAL, INTENT(IN) :: BUG - - REAL, INTENT(INOUT), DIMENSION(NCH) :: RZEXC, SRFEXC, CATDEF, CAPAC, & - RUNSRF - - - INTEGER N - REAL srflw,rzflw,FLOW,EXCESS,TSC0,tsc2,rzave,rz0,wanom,rztot, & - rzx,btaux,ax,bx,rzdif, rzavemin - - -!**** - - - - - - - - - - - - - - - - - - - - - - - - - - - DO 100 N=1,NCH - -!**** Compute equivalent of root zone excess in non-saturated area: - rztot=rzeq(n)+rzexc(n) - if(ar1(n).ne.1.) then - !!! rzave=(rztot-ar1(n)*vgwmax(n))/(1.-ar1(n)) - !!! rzave=rzave*poros(n)/vgwmax(n) - rzave=rztot*poros(n)/vgwmax(n) - else - rzave=poros(n) - endif - -! updated warning statement, reichle+koster, 12 Aug 2014 -! -! Impose minimum of 1.e-4, rather than leaving positive values <1.e-4 unchanged. -! -reichle, 15 Jan 2016 - - if (LAND_FIX) then - rzavemin = 1.e-4 - else - rzavemin = 0. - end if - - if (rzave .le. rzavemin) then ! JP: could put rzavemin in catch_constants - rzave=1.e-4 - print*,'problem: rzave <= 1.e-4 in catchment',n - end if - - btaux=btau(n) - if (srfexc(n) .lt. 0.) btaux=btau(n)*(poros(n)/rzave) - rz0=amax1(0.001,rzave-srfexc(n)/(1000.*(-btaux))) - tsc0=atau(n)/(rz0**3.) - - tsc0=tsc0*3600. - if(tsc0.lt.dtstep) tsc0=dtstep - -! --------------------------------------------------------------------- - - SRFLW=SRFEXC(N)*DTSTEP/TSC0 - - IF(SRFLW < 0. ) SRFLW = FLWALPHA * SRFLW ! C05 change - -!rr following inserted by koster Sep 22, 2003 - rzdif=rzave/poros(n)-wpwet(n) -!**** No moisture transport up if rz at wilting; employ ramping. - if(rzdif.le.0. .and. srflw.lt.0.) srflw=0. - if(rzdif.gt.0. .and. rzdif.lt.0.01 & - .and. srflw.lt.0.) srflw=srflw*(rzdif/0.01) - RZEXC(N)=RZEXC(N)+SRFLW - SRFEXC(N)=SRFEXC(N)-SRFLW - -!**** Topography-dependent tsc2, between rzex and catdef - - rzx=rzexc(n)/vgwmax(n) - - if(rzx .gt. .01) then - ax=tsa1(n) - bx=tsb1(n) - elseif(rzx .lt. -.01) then - ax=tsa2(n) - bx=tsb2(n) - else - ax=tsa2(n)+(rzx+.01)*(tsa1(n)-tsa2(n))/.02 - bx=tsb2(n)+(rzx+.01)*(tsb1(n)-tsb2(n))/.02 - endif - - tsc2=exp(ax+bx*catdef(n)) - rzflw=rzexc(n)*tsc2*dtstep/3600. - - IF (CATDEF(N)-RZFLW .GT. CDCR2(N)) then - RZFLW=CATDEF(N)-CDCR2(N) - end if - - CATDEF(N)=CATDEF(N)-RZFLW - RZEXC(N)=RZEXC(N)-RZFLW - -!**** REMOVE ANY EXCESS FROM MOISTURE RESERVOIRS: - - IF(CAPAC(N) .GT. SATCAP(N)) THEN - RZEXC(N)=RZEXC(N)+CAPAC(N)-SATCAP(N) - CAPAC(N)=SATCAP(N) - ENDIF - - IF(RZEQ(N) + RZEXC(N) .GT. VGWMAX(N)) THEN - EXCESS=RZEQ(N)+RZEXC(N)-VGWMAX(N) - RZEXC(N)=VGWMAX(N)-RZEQ(N) - CATDEF(N)=CATDEF(N)-EXCESS - ENDIF - - IF(CATDEF(N) .LT. 0.) THEN - ! bug fix: RUNSRF in flux units (kg m-2 s-1) for consistency with partition() - ! reichle, 5 Feb 2022 - RUNSRF(N)=RUNSRF(N)-CATDEF(N)/DTSTEP - CATDEF(N)=0. - ENDIF - - 100 ENDDO - - RETURN - END SUBROUTINE RZDRAIN - - -!**** -!**** ----------------------------------------------------------------- -!**** ///////////////////////////////////////////////////////////////// -!**** ----------------------------------------------------------------- !**** !**** [ BEGIN RCUNST ] !**** diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 index 1a43bf00c..6abc35d7c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 @@ -52,7 +52,7 @@ MODULE lsm_routines PRIVATE - PUBLIC :: INTERC, SRUNOFF, BASE, PARTITION, RZEQUIL, gndtp0 + PUBLIC :: INTERC, SRUNOFF, RZDRAIN, BASE, PARTITION, RZEQUIL, gndtp0 PUBLIC :: SIBALB, catch_calc_soil_moist, catch_calc_subtile2tile PUBLIC :: gndtmp, catch_calc_tp, catch_calc_wtotl, catch_calc_ght, catch_calc_FT PUBLIC :: dampen_tc_oscillations, irrigation_rate @@ -89,27 +89,34 @@ MODULE lsm_routines !**** [ BEGIN INTERC ] !**** - SUBROUTINE INTERC ( & - NCH, DTSTEP, FWETC, FWETL, TRAINL, TRAINC,SMELT, & - SATCAP,BUG, & - CAPAC, & - THRUL, THRUC & + SUBROUTINE INTERC ( & + NCH, DTSTEP, FWETC, FWETL, & + TRAINL, TRAINC, SMELT, & ! [kg m-2 s-1] + SATCAP,BUG, & + CAPAC, & + THRUL, THRUC & ! [kg m-2] !!!! ) !**** !**** THIS ROUTINE USES THE PRECIPITATION FORCING TO DETERMINE !**** CHANGES IN INTERCEPTION AND SOIL MOISTURE STORAGE. !**** +!**** +!**** NOTE: Input precip & snow melt are in flux units [kg m-2 s-1] +!**** Output throughfall is in volume units [kg m-2] +!**** (added comments for clarity but did not change units to preserve 0-diff, reichle, 6 Feb 2022) + IMPLICIT NONE !**** - INTEGER, INTENT(IN) :: NCH - REAL, INTENT(IN) :: DTSTEP, FWETC, FWETL - REAL, INTENT(IN), DIMENSION(NCH) :: TRAINL, TRAINC, SMELT, SATCAP - LOGICAL, INTENT(IN) :: BUG + INTEGER, INTENT(IN) :: NCH + REAL, INTENT(IN) :: DTSTEP, FWETC, FWETL + REAL, INTENT(IN), DIMENSION(NCH) :: TRAINL, TRAINC, SMELT ! [kg m-2 s-1] (flux units) + REAL, INTENT(IN), DIMENSION(NCH) :: SATCAP + LOGICAL, INTENT(IN) :: BUG - REAL, INTENT(INOUT), DIMENSION(NCH) :: CAPAC + REAL, INTENT(INOUT), DIMENSION(NCH) :: CAPAC - REAL, INTENT(OUT), DIMENSION(NCH) :: THRUC, THRUL + REAL, INTENT(OUT), DIMENSION(NCH) :: THRUC, THRUL ! [kg m-2] ("volume" units) !!! INTEGER CHNO REAL WETINT, WATADD, CAVAIL, THRU1, THRU2, XTCORR,SMPERS @@ -199,12 +206,7 @@ SUBROUTINE INTERC ( & THRUL(CHNO)=AMAX1(0., THRUL(CHNO)) THRUC(CHNO)=AMAX1(0., THRUC(CHNO)) - 100 CONTINUE ! ends do loop??? -!**** - ! return throughfall as a flux (kg m-2 s-1) for consistency with units of precip inputs - ! - reichle, 5 Feb 2022 - !!THRUL = THRUL/DTSTEP - !!THRUC = THRUC/DTSTEP + 100 CONTINUE !**** RETURN END SUBROUTINE INTERC @@ -217,26 +219,33 @@ END SUBROUTINE INTERC !**** /////////////////////////////////////////////////// !**** =================================================== - SUBROUTINE SRUNOFF ( & - NCH,DTSTEP,UFW4RO, FWETC, FWETL, AR1,ar2,ar4, THRUL_IN,THRUC_IN, & - frice,tp1,srfmx, BUG, & - SRFEXC,RUNSRF, & - QINFIL & + SUBROUTINE SRUNOFF ( & + NCH,DTSTEP,UFW4RO, FWETC, FWETL, AR1,ar2,ar4, & + THRUL_VOL, THRUC_VOL, & ! [kg m-2] ("volume" units) !!! + frice,tp1,srfmx, BUG, & + SRFEXC, & + RUNSRF, & ! [kg m-2 s-1] (flux units) + QINFIL & ! [kg m-2 s-1] (flux units) ) +!**** NOTE: Input throughfall is in volume units, as are calcs throughout this subroutine [kg m-2] +!**** Input-output surface runoff and output infiltration are in flux units [kg m-2 s-1] +!**** (added comments and clarified variable names, left throughfall units +!**** unchanged to preserve 0-diff, reichle, 6 Feb 2022) + IMPLICIT NONE + INTEGER, INTENT(IN) :: NCH + REAL, INTENT(IN) :: DTSTEP, FWETC, FWETL + LOGICAL, INTENT(IN) :: UFW4RO + REAL, INTENT(IN), DIMENSION(NCH) :: AR1, ar2, ar4, frice, tp1, srfmx + REAL, INTENT(IN), DIMENSION(NCH) :: THRUL_VOL, THRUC_VOL ! [kg m-2] + LOGICAL, INTENT(IN) :: BUG - INTEGER, INTENT(IN) :: NCH - REAL, INTENT(IN) :: DTSTEP, FWETC, FWETL - LOGICAL, INTENT (IN):: UFW4RO - REAL, INTENT(IN), DIMENSION(NCH) :: AR1, ar2, ar4, frice, tp1, & - srfmx, THRUL_IN, THRUC_IN - LOGICAL, INTENT(IN) :: BUG - - REAL, INTENT(INOUT), DIMENSION(NCH) :: SRFEXC ,RUNSRF + REAL, INTENT(INOUT), DIMENSION(NCH) :: SRFEXC + REAL, INTENT(INOUT), DIMENSION(NCH) :: RUNSRF ! [kg m-2 s-1] - REAL, INTENT(OUT), DIMENSION(NCH) :: QINFIL + REAL, INTENT(OUT), DIMENSION(NCH) :: QINFIL ! [kg m-2 s-1] INTEGER N REAL deficit,srun0,frun,qin, qinfil_l, qinfil_c, qcapac, excess_infil, srunc, srunl, ptotal @@ -244,11 +253,12 @@ SUBROUTINE SRUNOFF ( & !**** - - - - - - - - - - - - - - - - - - - - - - - - - - ! calculations throughout srunoff() are in "volumes" (kg m-2) - ! convert input fluxes to volumes, reichle, 5 Feb 2022 - THRUL = THRUL_IN !!* DTSTEP - THRUC = THRUC_IN !!* DTSTEP - RUNSRF = RUNSRF * DTSTEP + ! calculations throughout srunoff() are in "volume" units [kg m-2] + + THRUL = THRUL_VOL ! introduced *_VOL variables for clarity + THRUC = THRUC_VOL + + RUNSRF = RUNSRF * DTSTEP ! convert input surface runoff from flux to "volume" units DO N=1,NCH @@ -311,7 +321,8 @@ SUBROUTINE SRUNOFF ( & endif SRFEXC(N)=SRFEXC(N)+QIN - ! convert outputs back to flux units (kg m-2 s-1) + + ! convert outputs to flux units [kg m-2 s-1] RUNSRF(N)=RUNSRF(N)/DTSTEP QINFIL(N)=QIN/DTSTEP @@ -325,6 +336,145 @@ END SUBROUTINE SRUNOFF !**** ///////////////////////////////////////////////////////////////// !**** ----------------------------------------------------------------- !**** + SUBROUTINE RZDRAIN ( & + NCH,DTSTEP,VGWMAX,SATCAP,RZEQ,AR1,WPWET, & + tsa1,tsa2,tsb1,tsb2,atau,btau,CDCR2,poros,BUG, & + CAPAC,RZEXC,SRFEXC,CATDEF, & + RUNSRF & ! [kg m-2 s-1] + ) + +!----------------------------------------------------------------- +! defines drainage timescales: +! - tsc0, between srfex and rzex +! - tsc2, between rzex and catdef +! then defines correponding drainages +! and updates the water contents +!----------------------------------------------------------------- + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: NCH + REAL, INTENT(IN) :: DTSTEP + REAL, INTENT(IN), DIMENSION(NCH) :: VGWMAX, SATCAP, RZEQ, AR1, wpwet + REAL, INTENT(IN), DIMENSION(NCH) :: tsa1, tsa2, tsb1, tsb2, atau, btau, CDCR2, poros + LOGICAL, INTENT(IN) :: BUG + + REAL, INTENT(INOUT), DIMENSION(NCH) :: RZEXC, SRFEXC, CATDEF, CAPAC + REAL, INTENT(INOUT), DIMENSION(NCH) :: RUNSRF ! [kg m-2 s-1] + + INTEGER N + REAL srflw,rzflw,FLOW,EXCESS,TSC0,tsc2,rzave,rz0,wanom,rztot, & + rzx,btaux,ax,bx,rzdif, rzavemin + + +!**** - - - - - - - - - - - - - - - - - - - - - - - - - + + DO 100 N=1,NCH + +!**** Compute equivalent of root zone excess in non-saturated area: + rztot=rzeq(n)+rzexc(n) + if(ar1(n).ne.1.) then + !!! rzave=(rztot-ar1(n)*vgwmax(n))/(1.-ar1(n)) + !!! rzave=rzave*poros(n)/vgwmax(n) + rzave=rztot*poros(n)/vgwmax(n) + else + rzave=poros(n) + endif + +! updated warning statement, reichle+koster, 12 Aug 2014 +! +! Impose minimum of 1.e-4, rather than leaving positive values <1.e-4 unchanged. +! -reichle, 15 Jan 2016 + + if (LAND_FIX) then + rzavemin = 1.e-4 + else + rzavemin = 0. + end if + + if (rzave .le. rzavemin) then ! JP: could put rzavemin in catch_constants + rzave=1.e-4 + print*,'problem: rzave <= 1.e-4 in catchment',n + end if + + btaux=btau(n) + if (srfexc(n) .lt. 0.) btaux=btau(n)*(poros(n)/rzave) + rz0=amax1(0.001,rzave-srfexc(n)/(1000.*(-btaux))) + tsc0=atau(n)/(rz0**3.) + + tsc0=tsc0*3600. + if(tsc0.lt.dtstep) tsc0=dtstep + +! --------------------------------------------------------------------- + + SRFLW=SRFEXC(N)*DTSTEP/TSC0 + + IF(SRFLW < 0. ) SRFLW = FLWALPHA * SRFLW ! C05 change + +!rr following inserted by koster Sep 22, 2003 + rzdif=rzave/poros(n)-wpwet(n) +!**** No moisture transport up if rz at wilting; employ ramping. + if(rzdif.le.0. .and. srflw.lt.0.) srflw=0. + if(rzdif.gt.0. .and. rzdif.lt.0.01 & + .and. srflw.lt.0.) srflw=srflw*(rzdif/0.01) + RZEXC(N)=RZEXC(N)+SRFLW + SRFEXC(N)=SRFEXC(N)-SRFLW + +!**** Topography-dependent tsc2, between rzex and catdef + + rzx=rzexc(n)/vgwmax(n) + + if(rzx .gt. .01) then + ax=tsa1(n) + bx=tsb1(n) + elseif(rzx .lt. -.01) then + ax=tsa2(n) + bx=tsb2(n) + else + ax=tsa2(n)+(rzx+.01)*(tsa1(n)-tsa2(n))/.02 + bx=tsb2(n)+(rzx+.01)*(tsb1(n)-tsb2(n))/.02 + endif + + tsc2=exp(ax+bx*catdef(n)) + rzflw=rzexc(n)*tsc2*dtstep/3600. + + IF (CATDEF(N)-RZFLW .GT. CDCR2(N)) then + RZFLW=CATDEF(N)-CDCR2(N) + end if + + CATDEF(N)=CATDEF(N)-RZFLW + RZEXC(N)=RZEXC(N)-RZFLW + +!**** REMOVE ANY EXCESS FROM MOISTURE RESERVOIRS: + + IF(CAPAC(N) .GT. SATCAP(N)) THEN + RZEXC(N)=RZEXC(N)+CAPAC(N)-SATCAP(N) + CAPAC(N)=SATCAP(N) + ENDIF + + IF(RZEQ(N) + RZEXC(N) .GT. VGWMAX(N)) THEN + EXCESS=RZEQ(N)+RZEXC(N)-VGWMAX(N) + RZEXC(N)=VGWMAX(N)-RZEQ(N) + CATDEF(N)=CATDEF(N)-EXCESS + ENDIF + + IF(CATDEF(N) .LT. 0.) THEN + ! bug fix: RUNSRF in flux units [kg m-2 s-1] for consistency with partition() + ! reichle, 6 Feb 2022 + RUNSRF(N)=RUNSRF(N)-CATDEF(N)/DTSTEP + CATDEF(N)=0. + ENDIF + + 100 ENDDO + + RETURN + END SUBROUTINE RZDRAIN + +!**** ----------------------------------------------------------------- +!**** ///////////////////////////////////////////////////////////////// +!**** ----------------------------------------------------------------- +!**** + SUBROUTINE BASE ( & NCH,DTSTEP,BF1,BF2,BF3,CDCR1,FRICE,COND,GNU, & CATDEF, & From ff719ead1d9168cf2bcb14f6e71e0af8e4a2658b Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sun, 6 Feb 2022 13:17:28 -0500 Subject: [PATCH 5/5] clarification of RUNSRF units in subroutine partition() --- .../GEOSland_GridComp/Shared/lsm_routines.F90 | 38 ++++++++++--------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 index 6abc35d7c..b2766d3a6 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 @@ -529,33 +529,35 @@ END SUBROUTINE BASE !**** ----------------------------------------------------------------- !**** - SUBROUTINE PARTITION ( & - NCH,DTSTEP,DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2, & - PSIS,BEE,poros,WPWET, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4,BUG, & - srfexc,catdef,runsrf, & - AR1, AR2, AR4, srfmx, srfmn, & - SWSRF1,SWSRF2,SWSRF4,RZI & + SUBROUTINE PARTITION ( & + NCH,DTSTEP,DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2, & + PSIS,BEE,poros,WPWET, & + ars1,ars2,ars3,ara1,ara2,ara3,ara4, & + arw1,arw2,arw3,arw4,BUG, & + srfexc,catdef, & + runsrf, & ! [kg m-2 s-1] + AR1, AR2, AR4, srfmx, srfmn, & + SWSRF1,SWSRF2,SWSRF4,RZI & ) IMPLICIT NONE ! ------------------------------------------------------------------- - INTEGER, INTENT(IN) :: NCH + INTEGER, INTENT(IN) :: NCH - REAL, INTENT(IN) :: DTSTEP - REAL, INTENT(IN), DIMENSION(NCH) :: DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2, & - PSIS,BEE,poros,WPWET, & - ars1,ars2,ars3,ara1,ara2,ara3,ara4, & - arw1,arw2,arw3,arw4 + REAL, INTENT(IN) :: DTSTEP + REAL, INTENT(IN), DIMENSION(NCH) :: DZSF,RZEXC,RZEQ,VGWMAX,CDCR1,CDCR2, & + PSIS,BEE,poros,WPWET, & + ars1,ars2,ars3,ara1,ara2,ara3,ara4, & + arw1,arw2,arw3,arw4 - LOGICAL, INTENT(IN) :: BUG + LOGICAL, INTENT(IN) :: BUG ! ------------------------------------------------------------------- - REAL, INTENT(INOUT), DIMENSION(NCH) :: srfexc,catdef,runsrf + REAL, INTENT(INOUT), DIMENSION(NCH) :: srfexc,catdef + REAL, INTENT(INOUT), DIMENSION(NCH) :: runsrf ! [kg m-2 s-1] ! ------------------------------------------------------------------- - REAL, INTENT(OUT), DIMENSION(NCH) :: AR1, AR2, AR4, srfmx, srfmn, & - SWSRF1, SWSRF2, SWSRF4, RZI + REAL, INTENT(OUT), DIMENSION(NCH) :: AR1, AR2, AR4, srfmx, srfmn, & + SWSRF1, SWSRF2, SWSRF4, RZI ! ------------------------------------------------------------------- INTEGER :: N