From 8b1de2966001e6b00bc05a739c1262ef80ad020a Mon Sep 17 00:00:00 2001 From: Benjamin Sulman Date: Thu, 10 Jun 2021 13:34:56 -0400 Subject: [PATCH] Bsulman/water balance fix (#4) * 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 --- .../elm/src/biogeophys/SoilHydrologyMod.F90 | 56 +++++++++++-------- .../src/biogeophys/SoilWaterMovementMod.F90 | 4 +- components/elm/src/main/pftvarcon.F90 | 3 + 3 files changed, 39 insertions(+), 24 deletions(-) diff --git a/components/elm/src/biogeophys/SoilHydrologyMod.F90 b/components/elm/src/biogeophys/SoilHydrologyMod.F90 index 49722c589de..aea2bae70e2 100644 --- a/components/elm/src/biogeophys/SoilHydrologyMod.F90 +++ b/components/elm/src/biogeophys/SoilHydrologyMod.F90 @@ -183,7 +183,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 @@ -198,7 +198,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 @@ -211,7 +211,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 @@ -321,7 +321,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f #endif use elm_time_manager , only : get_curr_date, get_curr_time #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 elm_varcon , only : secspday ! @@ -384,7 +384,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( & @@ -615,8 +615,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 @@ -655,17 +656,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 @@ -759,9 +750,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 @@ -770,6 +762,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 @@ -1138,12 +1151,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 @@ -1678,7 +1690,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 diff --git a/components/elm/src/biogeophys/SoilWaterMovementMod.F90 b/components/elm/src/biogeophys/SoilWaterMovementMod.F90 index 143f61ee390..5c0553ef851 100644 --- a/components/elm/src/biogeophys/SoilWaterMovementMod.F90 +++ b/components/elm/src/biogeophys/SoilWaterMovementMod.F90 @@ -702,7 +702,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 @@ -760,7 +760,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 diff --git a/components/elm/src/main/pftvarcon.F90 b/components/elm/src/main/pftvarcon.F90 index 8ffc9e22d1a..b5eacdae502 100644 --- a/components/elm/src/main/pftvarcon.F90 +++ b/components/elm/src/main/pftvarcon.F90 @@ -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 @@ -1094,6 +1095,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.)