Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

bug fix: units of RUNSRF in Catchment #535

Merged
merged 6 commits into from
Feb 10, 2022
Merged

bug fix: units of RUNSRF in Catchment #535

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
6dbc8c7
bug fix: units of RUNSRF; cleanup: units of THRU[x] (catchment.F90,…
gmao-rreichle Feb 5, 2022
ebaa898
undoing unit conversions for THRU[x] (hoping to avoid non-zero-diff r…
gmao-rreichle Feb 5, 2022
30a6556
fixed bug in previous commit
gmao-rreichle Feb 6, 2022
0ec04ea
additional comments and cleanup:
gmao-rreichle Feb 6, 2022
ff719ea
clarification of RUNSRF units in subroutine partition()
gmao-rreichle Feb 6, 2022
0b57244
Merge branch 'develop' into bugfix/rreichle/Catchment_RUNSRF_units
sdrabenh Feb 10, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
135 changes: 5 additions & 130 deletions ...idComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/Shared/catchmentCN.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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 ( &
Expand Down
150 changes: 7 additions & 143 deletions ...Sphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/catchment.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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: &
Expand Down Expand Up @@ -1649,143 +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
RUNSRF(N)=RUNSRF(N)-CATDEF(N)
CATDEF(N)=0.
ENDIF

100 ENDDO

RETURN
END SUBROUTINE RZDRAIN


!****
!**** -----------------------------------------------------------------
!**** /////////////////////////////////////////////////////////////////
!**** -----------------------------------------------------------------
!****
!**** [ BEGIN RCUNST ]
!****
Expand Down
Loading