Skip to content

Commit

Permalink
Bsulman/water balance fix (#4)
Browse files Browse the repository at this point in the history 
* Initial test for water balance error corrections

* Tentatively working fix attempt from Ben

* Cleanup write statements etc

* Speed up surface water transfer and fix mistake with tide timing

* Move surface flow rate scale to settable parameter

Co-authored-by: omearata@ornl.gov <omearata@ornl.gov>
  • Loading branch information
@bsulman @omearata
bsulman and omearata committed Jun 9, 2023
1 parent 1549c6d commit 06a189d
Show file tree
Hide file tree
Showing 3 changed files with 39 additions and 24 deletions.
56 changes: 34 additions & 22 deletions components/elm/src/biogeophys/SoilHydrologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,7 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &

#if (defined MARSH)
if (c .eq. 1) fsat(c) = 1.0 * exp(-3.0_r8/humhol_ht*(zwt(c))) !at 30cm, hummock saturated at 5% changed to 0.1 TAO
if (c .eq. 2) fsat(c) = min(1.0 * exp(-3.0_r8/humhol_ht*(zwt(c)-h2osfc(c)/1000.+0.15_r8)), 1._r8) !TAO 0.3 t0 0.1, 0.15 to 0.35 !bsulman: what does 0.15 represent?
if (c .eq. 2) fsat(c) = min(1.0 * exp(-3.0_r8/humhol_ht*(zwt(c)-h2osfc(c)/1000.+humhol_ht)), 1._r8) !TAO 0.3 t0 0.1, 0.15 to 0.35 !bsulman: what does 0.15 represent?
#endif
! use perched water table to determine fsat (if present)
if ( frost_table(c) > zwt(c)) then
Expand All @@ -195,7 +195,7 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &

#if (defined MARSH)
if (c .eq. 1) fsat(c) = 1.0 * exp(-3.0_r8/humhol_ht*(zwt(c))) !at 30cm, hummock saturated at 5%
if (c .eq. 2) fsat(c) = min(1.0 * exp(-3.0_r8/humhol_ht*(zwt(c)-h2osfc(c)/1000.+0.15_r8)), 1._r8) !TAO 0.3 t 0.1, 0.15 to 0.35
if (c .eq. 2) fsat(c) = min(1.0 * exp(-3.0_r8/humhol_ht*(zwt(c)-h2osfc(c)/1000.+humhol_ht)), 1._r8) !TAO 0.3 t 0.1, 0.15 to 0.35
#endif

else
Expand All @@ -208,7 +208,7 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &

#if (defined MARSH)
if (c .eq. 1) fsat(c) = 1.0 * exp(-3.0_r8/humhol_ht*(zwt(c))) !at 30cm, hummock saturated at 5%
if (c .eq. 2) fsat(c) = min(1.0 * exp(-3.0_r8/humhol_ht*(zwt(c)-h2osfc(c)/1000.+0.15_r8)), 1._r8) !TAO 0.3 t 1.5, 0.15 to 0.35
if (c .eq. 2) fsat(c) = min(1.0 * exp(-3.0_r8/humhol_ht*(zwt(c)-h2osfc(c)/1000.+humhol_ht)), 1._r8) !TAO 0.3 t 1.5, 0.15 to 0.35
#endif
endif
if (origflag == 1) then
Expand Down Expand Up @@ -317,7 +317,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
use pftvarcon , only : humhol_ht, humhol_dist, hum_frac, qflx_h2osfc_surfrate
#endif
#if defined MARSH
use pftvarcon , only : num_tide_comps, tide_baseline,tide_coeff_period, tide_coeff_phase, tide_coeff_amp
use pftvarcon , only : num_tide_comps, tide_baseline,tide_coeff_period, tide_coeff_phase, tide_coeff_amp,sfcflow_ratescale
#endif
use clm_time_manager , only : get_step_size, get_curr_date, get_curr_time
use elm_varcon , only : secspday
Expand Down Expand Up @@ -380,7 +380,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
integer :: yr, mon, day, tod !
integer :: days, seconds !
integer :: ii
real(r8) :: h2osfc_before
real(r8) :: h2osfc_tide
!-----------------------------------------------------------------------

associate( &
Expand Down Expand Up @@ -592,8 +592,9 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
!4. soil infiltration and h2osfc "run-on"
qflx_infl(c) = qflx_in_soil(c) - qflx_infl_excess(c)
qflx_in_h2osfc(c) = qflx_in_h2osfc(c) + qflx_infl_excess(c)
#if (defined HUM_HOL || defined MARSH)
#if (defined HUM_HOL)
if (c .eq. 1) then
! qflx_surf is surface runoff. So this means there is no surface runoff from hollow (channel) to hummock (marsh?)
qflx_surf(1) = qflx_surf(1) + qflx_in_h2osfc(c)
qflx_surf_input(2) = qflx_surf_input(2) + qflx_in_h2osfc(c)
qflx_in_h2osfc(c) = 0._r8 !TAO 22/8/2018 changing to sin function gave an error
Expand All @@ -620,17 +621,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
if (h2osfc(c) .gt. 0._r8 .and. c==1) then
qflx_h2osfc_surf(c) = min(qflx_h2osfc_surfrate*h2osfc(c)**2.0_r8,h2osfc(c) / dtime)
else if (c .eq. 2) then
call get_curr_time (days, seconds)
qflx_h2osfc_surf(c) = 0._r8
! bsulman : Changed to use flexible set of parameters up to full NOAA tidal components (37 coefficients)
! Tidal cycle is the sum of all the sinusoidal components
h2osfc_before = h2osfc(c)
h2osfc(c) = 0.0_r8
do ii=1,num_tide_comps
h2osfc(c) = h2osfc(c) + tide_coeff_amp(ii) * sin(2.0_r8*SHR_CONST_PI*(1/tide_coeff_period(ii)*(days*secspday+seconds) + tide_coeff_phase(ii)))
enddo
h2osfc(c) = max(h2osfc(c) + tide_baseline, 0.0)
qflx_tide(c) = (h2osfc(c)-h2osfc_before)/dtime

#else
if(h2osfc(c) >= h2osfc_thresh(c) .and. h2osfcflag/=0) then
Expand Down Expand Up @@ -714,9 +705,10 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
zwt_ho = zwt_ho - h2osfc(2)/1000._r8 !DMR 4/29/13
end if
!DMR 12/4/2015
if (icefrac(1,min(jwt(1)+1,nlevsoi)) .ge. 0.01_r8 .or. &
icefrac(2,min(jwt(2)+1,nlevsoi)) .ge. 0.01_r8) then
!turn off lateral transport if any ice is present at or below
if (icefrac(1,min(jwt(1)+1,nlevsoi)) .ge. 0.90_r8 .or. &
icefrac(2,min(jwt(2)+1,nlevsoi)) .ge. 0.90_r8) then
!turn off lateral transport if any ice is present at or below,
!changed from 0.01 to 0.90 TAO 6/4/2021
!water table
qflx_lat_aqu(:) = 0._r8
else
Expand All @@ -725,6 +717,27 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
qflx_lat_aqu(2) = -2._r8/(1._r8/ka_hu+1._r8/ka_ho) * (zwt_hu-zwt_ho- &
humhol_ht) / humhol_dist * sqrt(hum_frac/hol_frac)
endif
! bsulman : Changed to use flexible set of parameters up to full NOAA tidal components (37 coefficients)
! Tidal cycle is the sum of all the sinusoidal components
#ifdef MARSH
call get_curr_time(days, seconds)
h2osfc_tide = 0.0_r8
do ii=1,num_tide_comps
h2osfc_tide = h2osfc_tide + tide_coeff_amp(ii) * sin(2.0_r8*SHR_CONST_PI*(1/tide_coeff_period(ii)*(days*secspday+seconds) + tide_coeff_phase(ii)))
enddo
h2osfc_tide = max(h2osfc_tide + tide_baseline, 0.0)
! qflx_tide(c) = (h2osfc(c)-h2osfc_before)/dtime
qflx_lat_aqu(2) = qflx_lat_aqu(2) + (h2osfc_tide-h2osfc(c))/dtime

! If flooded water surface of one column is higher than the other, add faster flow since aquifer transfer (ka parameters) is slow
if(h2osfc(2)>0 .and. h2osfc(2)>(h2osfc(1)+humhol_ht*1000.0)) then
qflx_lat_aqu(2) = qflx_lat_aqu(2) - min((h2osfc(2)-(h2osfc(1)+humhol_ht*1000.0))*sfcflow_ratescale,h2osfc(2)*0.5/dtime)
qflx_lat_aqu(1) = qflx_lat_aqu(1) + min((h2osfc(2)-(h2osfc(1)+humhol_ht*1000.0))*sfcflow_ratescale,h2osfc(2)*0.5/dtime)
elseif(h2osfc(1)>0 .and. h2osfc(1)>(h2osfc(2)-humhol_ht*1000.0)) then
qflx_lat_aqu(2) = qflx_lat_aqu(2) + min((h2osfc(1)-(h2osfc(2)-humhol_ht*1000.0))*sfcflow_ratescale,h2osfc(1)*0.5/dtime)
qflx_lat_aqu(1) = qflx_lat_aqu(1) - min((h2osfc(1)-(h2osfc(2)-humhol_ht*1000.0))*sfcflow_ratescale,h2osfc(1)*0.5/dtime)
endif
#endif
endif
#endif

Expand Down Expand Up @@ -1093,12 +1106,11 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, fil
end do
else ! deepening water table (negative lateral flux)
!Remove from surface water first if available
if (h2osfc(c) .gt. 0 .and. maxval(icefrac(c,1:jwt(c)+1)) .le. 0.9) then
if (h2osfc(c) .gt. 0) then ! .and. maxval(icefrac(c,1:jwt(c)+1)) .le. 0.9) then
h2osfc(c) = h2osfc(c) + qflx_lat_aqu_tot
qflx_lat_aqu_tot = 0._r8
if (h2osfc(c) .lt. 0) then
qflx_lat_aqu_tot = h2osfc(c)
call get_curr_time(days, seconds)
h2osfc(c) = 0._r8
end if
end if
Expand Down Expand Up @@ -1633,7 +1645,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
endif

! bsulman: Does MARSH need this deep_seep stuff?
#if (defined HUM_HOL || defined MARSH )
#if (defined HUM_HOL)
deep_seep = 0._r8 !100.0_r8 / 365._r8 / 86400._r8 !rate per second
!changes for hummock hollow topography
if (c .eq. 1) then !hummock
Expand Down
4 changes: 2 additions & 2 deletions components/elm/src/biogeophys/SoilWaterMovementMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -703,7 +703,7 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
qout(c,j) = 0._r8
dqodw1(c,j) = 0._r8
rmx(c,j) = qin(c,j) - qout(c,j) - qflx_rootsoi_col(c,j)
#if (defined HUM_HOL)
#if (defined HUM_HOL || defined MARSH)
if (j > jwt(c)+1) then !Water table above this layer
rmx(c,j) = qin(c,j) - qout(c,j)
else if (j == jwt(c)+1) then !water table in this layer
Expand Down Expand Up @@ -761,7 +761,7 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
end if

rmx(c,j) = qin(c,j) - qout(c,j) - qflx_rootsoi_col(c,j)
#if (defined HUM_HOL)
#if (defined HUM_HOL || defined MARSH)
if (j > jwt(c)+1) then !Water table above this layer
rmx(c,j) = qin(c,j) - qout(c,j)
else if (j == jwt(c)+1) then !water table in this layer
Expand Down
3 changes: 3 additions & 0 deletions components/elm/src/main/pftvarcon.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -328,6 +328,7 @@ module pftvarcon
real(r8),allocatable :: tide_coeff_amp(:) ! Amplitude of tide component (mm)
real(r8),allocatable :: tide_coeff_period(:) ! Period of tide component (s)
real(r8),allocatable :: tide_coeff_phase(:) ! Phase shift of tide component (s)
real(r8) :: sfcflow_ratescale ! Rate scale for surface water flow across columns (s-1)
!endif
!phenology
real(r8) :: phen_a
Expand Down Expand Up @@ -1093,6 +1094,8 @@ subroutine pftconrd
tide_coeff_phase(1) = 513.4328
tide_coeff_phase(2) = 0.0
endif
call ncd_io('sfcflow_ratescale',sfcflow_ratescale, 'read', ncid, readvar=readv, posNOTonfile=.true.)
if (.not. readv) sfcflow_ratescale = 7.0e-5_r8 ! Probably better to have default be zero for safety
#endif

call ncd_io('phen_a', phen_a, 'read', ncid, readvar=readv, posNOTonfile=.true.)
Expand Down

0 comments on commit 06a189d

Please sign in to comment.