diff --git a/src/clm5/biogeophys/SnowHydrologyMod.F90 b/src/clm5/biogeophys/SnowHydrologyMod.F90 index 15f3f5767..2aa2df78a 100644 --- a/src/clm5/biogeophys/SnowHydrologyMod.F90 +++ b/src/clm5/biogeophys/SnowHydrologyMod.F90 @@ -287,6 +287,10 @@ subroutine SnowWater(bounds, & real(r8) :: vol_ice(bounds%begc:bounds%endc,-nlevsno+1:0) ! partial volume of ice lens in layer real(r8) :: eff_porosity(bounds%begc:bounds%endc,-nlevsno+1:0) ! effective porosity = porosity - vol_ice real(r8) :: mss_liqice(bounds%begc:bounds%endc,-nlevsno+1:0) ! mass of liquid+ice in a layer +#ifdef COUP_OAS_PFL + real(r8) :: snowliq_to_soil(bounds%begc:bounds%endc) ! liquid from the bottom snow layer to soil layer [mm] + ! <0: snow liquid deficit refilled from soil +#endif !----------------------------------------------------------------------- associate( & @@ -333,6 +337,9 @@ subroutine SnowWater(bounds, & do fc = 1,num_snowc c = filter_snowc(fc) l=col%landunit(c) +#ifdef COUP_OAS_PFL + snowliq_to_soil(c) = 0._r8 +#endif wgdif = h2osoi_ice(c,snl(c)+1) & + frac_sno_eff(c) * (qflx_dew_snow(c) - qflx_sub_snow(c)) * dtime @@ -352,6 +359,11 @@ subroutine SnowWater(bounds, & if (wgdif >= 0._r8) exit h2osoi_liq(c,j) = 0._r8 h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + wgdif +#ifdef COUP_OAS_PFL + ! Deficit moved from the bottom snow layer into the soil, which is + ! overwritten by ParFlow. Take it from ParFlow instead. + if (j == 0) snowliq_to_soil(c) = wgdif +#endif enddo end if end do @@ -548,6 +560,9 @@ subroutine SnowWater(bounds, & qflx_top_soil(c) = (qout(c) / dtime) & + (1.0_r8 - frac_sno_eff(c)) * qflx_rain_grnd(c) +#ifdef COUP_OAS_PFL + qflx_top_soil(c) = qflx_top_soil(c) + snowliq_to_soil(c) / dtime +#endif int_snow(c) = int_snow(c) + frac_sno_eff(c) & * (qflx_dew_snow(c) + qflx_dew_grnd(c) + qflx_rain_grnd(c)) * dtime end do @@ -831,6 +846,9 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, & snw_rds => waterstate_inst%snw_rds_col , & ! Output: [real(r8) (:,:) ] effective snow grain radius (col,lyr) [microns, m^-6] qflx_sl_top_soil => waterflux_inst%qflx_sl_top_soil_col , & ! Output: [real(r8) (:) ] liquid water + ice from layer above soil to top soil layer or sent to qflx_qrgwl (mm H2O/s) +#ifdef COUP_OAS_PFL + qflx_pfl_top => waterflux_inst%qflx_pfl_top_col , & ! Output: [real(r8) (:) ] water added to soil layer 1 after ParFlow exchange (mm H2O/s) +#endif snl => col%snl , & ! Output: [integer (:) ] number of snow layers dz => col%dz , & ! Output: [real(r8) (:,:) ] layer depth (m) @@ -874,7 +892,16 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, & ! use 0.01 to avoid runaway ice buildup if (h2osoi_ice(c,j) <= .01_r8) then if (ltype(l) == istsoil .or. urbpoi(l) .or. ltype(l) == istcrop) then +#ifdef COUP_OAS_PFL + if (j == 0 .and. col%hydrologically_active(c)) then + ! Pass the snow water to ParFlow. The ice phase is kept in eCLM. + qflx_pfl_top(c) = qflx_pfl_top(c) + (h2osoi_liq(c,j) + h2osoi_ice(c,j))/dtime + else + h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j) + end if +#else h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j) +#endif h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j) if (j == 0) then @@ -996,7 +1023,16 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, & ! this is where water is transfered from layer 0 (snow) to layer 1 (soil) if (ltype(l) == istsoil .or. urbpoi(l) .or. ltype(l) == istcrop) then h2osoi_liq(c,0) = 0.0_r8 +#ifdef COUP_OAS_PFL + if (col%hydrologically_active(c)) then + ! Pass the snow water to ParFlow + qflx_pfl_top(c) = qflx_pfl_top(c) + zwliq(c)/dtime + else + h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c) + end if +#else h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c) +#endif end if if (ltype(l) == istwet) then h2osoi_liq(c,0) = 0.0_r8 diff --git a/src/clm5/biogeophys/SoilHydrologyMod.F90 b/src/clm5/biogeophys/SoilHydrologyMod.F90 index 3b20ff626..385788360 100644 --- a/src/clm5/biogeophys/SoilHydrologyMod.F90 +++ b/src/clm5/biogeophys/SoilHydrologyMod.F90 @@ -450,6 +450,14 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f qflx_evap(c)=qflx_ev_soil(c) endif +#ifdef COUP_OAS_PFL + ! All surface input minus evaporation (and surface runoff) from soil and + ! surface water is passed to ParFlow, without the eCLM infiltration capacity. + qflx_infl(c) = qflx_top_soil(c) - qflx_surf(c) & + - (1.0_r8 - fsno - frac_h2osfc(c))*qflx_evap(c) & + - frac_h2osfc(c)*qflx_ev_h2osfc(c) + qflx_h2osfc_surf(c) = 0._r8 +#else !1. partition surface inputs between soil and h2osfc qflx_in_soil(c) = (1._r8 - frac_h2osfc(c)) * (qflx_top_soil(c) - qflx_surf(c)) qflx_in_h2osfc(c) = frac_h2osfc(c) * (qflx_top_soil(c) - qflx_surf(c)) @@ -545,6 +553,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f !7. remove drainage from h2osfc and add to qflx_infl h2osfc(c) = h2osfc(c) - qflx_h2osfc_drain(c) * dtime qflx_infl(c) = qflx_infl(c) + qflx_h2osfc_drain(c) +#endif else ! non-vegetated landunits (i.e. urban) use original CLM4 code if (snl(c) >= 0) then @@ -2279,6 +2288,9 @@ subroutine RenewCondensation(bounds, num_hydrologyc, filter_hydrologyc, & ! !LOCAL VARIABLES: integer :: c,j,fc,i ! indices real(r8) :: dtime ! land model time step (sec) +#ifdef COUP_OAS_PFL + real(r8) :: h2osoi_ice1_old ! ice of soil layer 1 before condensation/sublimation (kg/m2) +#endif !----------------------------------------------------------------------- associate( & @@ -2288,6 +2300,9 @@ subroutine RenewCondensation(bounds, num_hydrologyc, filter_hydrologyc, & frac_h2osfc => waterstate_inst%frac_h2osfc_col , & ! Input: [real(r8) (:) ] qflx_dew_grnd => waterflux_inst%qflx_dew_grnd_col , & ! Input: [real(r8) (:) ] ground surface dew formation (mm H2O /s) [+] qflx_dew_snow => waterflux_inst%qflx_dew_snow_col , & ! Input: [real(r8) (:) ] surface dew added to snow pack (mm H2O /s) [+] +#ifdef COUP_OAS_PFL + qflx_pfl_top => waterflux_inst%qflx_pfl_top_col , & ! Output: [real(r8) (:) ] water added to soil layer 1 after ParFlow exchange (mm H2O/s) +#endif qflx_sub_snow => waterflux_inst%qflx_sub_snow_col & ! Output: [real(r8) (:) ] sublimation rate from snow pack (mm H2O /s) [+] ) @@ -2303,7 +2318,13 @@ subroutine RenewCondensation(bounds, num_hydrologyc, filter_hydrologyc, & if (snl(c)+1 >= 1) then ! make consistent with how evap_grnd removed in infiltration +#ifdef COUP_OAS_PFL + ! Pass the dew to ParFlow. + qflx_pfl_top(c) = qflx_pfl_top(c) + (1._r8 - frac_h2osfc(c))*qflx_dew_grnd(c) + h2osoi_ice1_old = h2osoi_ice(c,1) +#else h2osoi_liq(c,1) = h2osoi_liq(c,1) + (1._r8 - frac_h2osfc(c))*qflx_dew_grnd(c) * dtime +#endif h2osoi_ice(c,1) = h2osoi_ice(c,1) + (1._r8 - frac_h2osfc(c))*qflx_dew_snow(c) * dtime if (qflx_sub_snow(c)*dtime > h2osoi_ice(c,1)) then qflx_sub_snow(c) = h2osoi_ice(c,1)/dtime @@ -2311,6 +2332,10 @@ subroutine RenewCondensation(bounds, num_hydrologyc, filter_hydrologyc, & else h2osoi_ice(c,1) = h2osoi_ice(c,1) - (1._r8 - frac_h2osfc(c)) * qflx_sub_snow(c) * dtime end if +#ifdef COUP_OAS_PFL + ! Pass the ice change to ParFlow in addition. SPo: check for eff. poro. coupling + qflx_pfl_top(c) = qflx_pfl_top(c) + (h2osoi_ice(c,1) - h2osoi_ice1_old) / dtime +#endif end if end do @@ -2369,10 +2394,13 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & associate( & dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) + qflx_pfl_top => waterflux_inst%qflx_pfl_top_col , & ! Input: [real(r8) (:) ] water added to soil layer 1 after ParFlow exchange (mm H2O/s) + qflx_h2osfc_to_ice => waterflux_inst%qflx_h2osfc_to_ice_col , & ! Input: [real(r8) (:) ] surface water converted to ice (mm H2O /s) qflx_snwcp_liq => waterflux_inst%qflx_snwcp_liq_col , & ! excess rainfall due to snow capping (mm H2O /s) [+] qflx_drain => waterflux_inst%qflx_drain_col , & ! sub-surface runoff (mm H2O /s) qflx_drain_perched => waterflux_inst%qflx_drain_perched_col , & ! perched wt sub-surface runoff (mm H2O /s) qflx_qrgwl => waterflux_inst%qflx_qrgwl_col , & ! qflx_surf at glaciers, wetlands, lakes (mm H2O /s) + qflx_surf => waterflux_inst%qflx_surf_col , & ! surface runoff (mm H2O /s) qflx_rsub_sat => waterflux_inst%qflx_rsub_sat_col , & ! soil saturation excess [mm h2o/s] qflx_infl => waterflux_inst%qflx_infl_col , & ! infiltration (mm H2O /s) qflx_rootsoi => waterflux_inst%qflx_rootsoi_col , & ! vegetation/soil water exchange (mm H2O/s) (+ = to atm) @@ -2395,11 +2423,20 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & qflx_parflow(c,j) = -qflx_rootsoi(c,j) !mm/s end if end do + + ! Water added to the top soil layer by eCLM after the ParFlow state was applied + ! (condensation/sublimation in RenewCondensation, snow water from CombineSnowLayers) + qflx_parflow(c,1) = qflx_parflow(c,1) + qflx_pfl_top(c) + ! Remove what froze into ice from surface water + qflx_parflow(c,1) = qflx_parflow(c,1) - qflx_h2osfc_to_ice(c) + ! Excess rainfall due to snow capping pass to ParFlow instead of river routing + qflx_parflow(c,1) = qflx_parflow(c,1) + qflx_snwcp_liq(c) + ! Compute subsurface run-off (mm/s) qflx_drain(c) = -sum(qflx_parflow(c,:)) qflx_drain_perched(c) = 0._r8 qflx_rsub_sat(c) = 0._r8 - qflx_qrgwl(c) = qflx_snwcp_liq(c) ! Set imbalance for snow capping + qflx_qrgwl(c) = 0._r8 ! Convert eCLM fluxes (mm/s) to ParFlow fluxes (1/hr): ! 1/hr = [mm/s] * [s/hr] * [m/mm] * [1/m] qflx_parflow(c,1:nlevsoi) = qflx_parflow(c,1:nlevsoi) * sec_per_hr * m_per_mm * (1._r8/dz(c,1:nlevsoi)) @@ -2412,6 +2449,17 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & qflx_drain(c) = 0._r8 ! This must be done for roofs and impervious road (walls will be zero) qflx_qrgwl(c) = qflx_snwcp_liq(c) + + ! Instead of leaving as river runoff, impervious urban runoff enters + ! ParFlow through layer 1. It is zeroed here so it is not counted twice. + qflx_parflow(c,1:nlevsoi) = 0._r8 + qflx_parflow(c,1) = qflx_surf(c) + qflx_qrgwl(c) + qflx_surf(c) = 0._r8 + qflx_qrgwl(c) = 0._r8 + + qflx_drain(c) = -sum(qflx_parflow(c,:)) + + qflx_parflow(c,1:nlevsoi) = qflx_parflow(c,1:nlevsoi) * sec_per_hr * m_per_mm * (1._r8/dz(c,1:nlevsoi)) end if end do end associate diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index 81f86e436..5feed75bd 100644 --- a/src/clm5/biogeophys/SoilWaterMovementMod.F90 +++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90 @@ -1503,6 +1503,7 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) smpmin => soilstate_inst%smpmin_col , & ! Input: [real(r8) (:) ] restriction for min of soil potential (mm) pfl_h2osoi_liq => waterstate_inst%pfl_h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] ParFlow soil water (mm) + qflx_pfl_top => waterflux_inst%qflx_pfl_top_col , & ! Output: [real(r8) (:) ] water added to soil layer 1 after exchange (mm H2O/s) pfl_psi => waterstate_inst%pfl_psi_col & ! Input: [real(r8) (:,:) ] ParFlow pressure head (mm) ) ! end associate statement @@ -1517,11 +1518,13 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & do j = 1, nlevgrnd h2osoi_ice(c,j) = min(h2osoi_ice(c,j), pfl_h2osoi_liq(c,j)) h2osoi_liq(c,j) = max(0._r8, pfl_h2osoi_liq(c,j) - h2osoi_ice(c,j)) - if (pfl_psi(c,j) <= 0) then - smp_l(c,j) = max(smpmin(c), pfl_psi(c,j)) - end if + ! Saturated/ponded cells (pfl_psi > 0) have zero matric potential + smp_l(c,j) = max(smpmin(c), min(0._r8, pfl_psi(c,j))) end do + ! Reset flux after passing ParFlow state to eCLM + qflx_pfl_top(c) = 0._r8 + ! ParFlow replaces the eCLM soil water solver. Recompute hk_l from the ! ParFlow water content using the same formulation as compute_hydraulic_properties. nlayers = nbedrock(c) diff --git a/src/clm5/biogeophys/WaterfluxType.F90 b/src/clm5/biogeophys/WaterfluxType.F90 index 098de9849..30c65edf5 100644 --- a/src/clm5/biogeophys/WaterfluxType.F90 +++ b/src/clm5/biogeophys/WaterfluxType.F90 @@ -74,6 +74,7 @@ module WaterfluxType real(r8), pointer :: qflx_rootsoi_col (:,:) ! col root and soil water exchange [mm H2O/s] [+ into root] #ifdef COUP_OAS_PFL real(r8), pointer :: qflx_parflow_col (:,:) ! col source/sink flux per soil layer sent to ParFlow [1/hr] [- out from root] + real(r8), pointer :: qflx_pfl_top_col (:) ! col water added to soil layer 1 by eCLM within time step, sent to ParFlow (mm H2O/s) #endif real(r8), pointer :: qflx_infl_col (:) ! col infiltration (mm H2O /s) real(r8), pointer :: qflx_surf_col (:) ! col surface runoff (mm H2O /s) @@ -219,6 +220,7 @@ subroutine InitAllocate(this, bounds) allocate(this%qflx_rootsoi_col (begc:endc,1:nlevsoi)) ; this%qflx_rootsoi_col (:,:) = nan #ifdef COUP_OAS_PFL allocate(this%qflx_parflow_col (begc:endc,1:nlevsoi)) ; this%qflx_parflow_col (:,:) = nan + allocate(this%qflx_pfl_top_col (begc:endc)) ; this%qflx_pfl_top_col (:) = nan #endif allocate(this%qflx_infl_col (begc:endc)) ; this%qflx_infl_col (:) = nan allocate(this%qflx_surf_col (begc:endc)) ; this%qflx_surf_col (:) = nan diff --git a/src/clm5/main/clm_driver.F90 b/src/clm5/main/clm_driver.F90 index befbba7a4..1a25eb490 100644 --- a/src/clm5/main/clm_driver.F90 +++ b/src/clm5/main/clm_driver.F90 @@ -1224,6 +1224,9 @@ subroutine clm_drv_init(bounds, & use WaterStateType , only : waterstate_type use WaterFluxType , only : waterflux_type use EnergyFluxType , only : energyflux_type +#ifdef COUP_OAS_PFL + use landunit_varcon , only : istsoil, istcrop +#endif ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -1257,6 +1260,7 @@ subroutine clm_drv_init(bounds, & #ifdef COUP_OAS_PFL pfl_psi => waterstate_inst%pfl_psi_col , & ! Input: [real(r8) (:,:) ] COUP_OAS_PFL pfl_h2osoi_liq => waterstate_inst%pfl_h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] COUP_OAS_PFL + h2osfc => waterstate_inst%h2osfc_col , & ! Output: [real(r8) (:) ] surface water (mm) #endif elai => canopystate_inst%elai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index with burying by snow esai => canopystate_inst%esai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index with burying by snow @@ -1313,6 +1317,11 @@ subroutine clm_drv_init(bounds, & g = col%gridcell(c) pfl_psi(c,:) = atm2lnd_inst%pfl_psi_grc(g,:) pfl_h2osoi_liq(c,:) = atm2lnd_inst%pfl_h2osoi_liq_grc(g,:) + ! ParFlow's positive pressure head of the top cell is the ponding depth. + ! Surface water is diagnosed from it to sync eCLM and ParFlow state. + if (lun%itype(col%landunit(c)) == istsoil .or. lun%itype(col%landunit(c)) == istcrop) then + h2osfc(c) = max(0._r8, pfl_psi(c,1)) + end if end if end do #endif