From 66db07d53dcbd6a4891c1763e47f15c01ad76d7c Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Thu, 3 Sep 2026 17:18:49 +0200 Subject: [PATCH 1/4] eCLM-ParFlow: Route impervious urban runoff into ParFlow - impervious urban runoff left the domain as river runoff without entering ParFlow - route it into ParFlow layer 1 - zeroed qflx_surf/qflx_qrgwl and write it to qflx_drain as for every other coupled column --- src/clm5/biogeophys/SoilHydrologyMod.F90 | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/clm5/biogeophys/SoilHydrologyMod.F90 b/src/clm5/biogeophys/SoilHydrologyMod.F90 index 3b20ff626..d428fc8bd 100644 --- a/src/clm5/biogeophys/SoilHydrologyMod.F90 +++ b/src/clm5/biogeophys/SoilHydrologyMod.F90 @@ -2373,6 +2373,7 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & 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) @@ -2412,6 +2413,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 From f46d798ef62398c2d62f009d5e9cfeea72ad6b28 Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Thu, 24 Sep 2026 10:39:46 +0200 Subject: [PATCH 2/4] eCLM-ParFlow: Couple runoff and overland flow - h2osfc from eCLM use ponding computed by ParFlow, adapt frach2osfc accordingly - Infiltration passes all surface input minus soil and surface-water evaporation to ParFlow without the eCLM qinmax capacity - account layer-1 water change (RenewCondensation, CombineSnowLayers), store pfl_top_sync after the sync and add the layer-1 mass change to qflx_parflow - remove qflx_h2osfc_to_ice (ponded water frozen into snow) from the ParFlow source term - route snow-capping liquid qflx_snwcp_liq to ParFlow layer 1 instead of qflx_qrgwl - SnowWater: liquid taken from soil layer 1 is taken from qflx_top_soil - smp_l set to 0 for saturated/ponded cells (pfl_psi > 0) instead of keeping the previous value --- src/clm5/biogeophys/HydrologyDrainageMod.F90 | 2 +- src/clm5/biogeophys/SnowHydrologyMod.F90 | 14 +++++++++ src/clm5/biogeophys/SoilHydrologyMod.F90 | 33 ++++++++++++++++++-- src/clm5/biogeophys/SoilWaterMovementMod.F90 | 10 ++++-- src/clm5/biogeophys/WaterStateType.F90 | 2 ++ src/clm5/main/clm_driver.F90 | 9 ++++++ 6 files changed, 64 insertions(+), 6 deletions(-) diff --git a/src/clm5/biogeophys/HydrologyDrainageMod.F90 b/src/clm5/biogeophys/HydrologyDrainageMod.F90 index e332aede5..b6c8a5f3d 100644 --- a/src/clm5/biogeophys/HydrologyDrainageMod.F90 +++ b/src/clm5/biogeophys/HydrologyDrainageMod.F90 @@ -135,7 +135,7 @@ subroutine HydrologyDrainage(bounds, & ! clm3.5/bld/usr.src/SoilHydrologyMod.F90 ! ignore drainage calculations in eCLM and instead pass these fluxes to ParFlow call ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & - num_urbanc, filter_urbanc, waterflux_inst) + num_urbanc, filter_urbanc, waterstate_inst, waterflux_inst) else if (use_aquifer_layer()) then call Drainage(bounds, num_hydrologyc, filter_hydrologyc, & diff --git a/src/clm5/biogeophys/SnowHydrologyMod.F90 b/src/clm5/biogeophys/SnowHydrologyMod.F90 index 15f3f5767..594f4b93f 100644 --- a/src/clm5/biogeophys/SnowHydrologyMod.F90 +++ b/src/clm5/biogeophys/SnowHydrologyMod.F90 @@ -287,6 +287,9 @@ 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) :: liq_from_soil(bounds%begc:bounds%endc) ! liquid taken from soil layer 1 to fill snow liquid deficit [mm] +#endif !----------------------------------------------------------------------- associate( & @@ -333,6 +336,9 @@ subroutine SnowWater(bounds, & do fc = 1,num_snowc c = filter_snowc(fc) l=col%landunit(c) +#ifdef COUP_OAS_PFL + liq_from_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 +358,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) liq_from_soil(c) = wgdif +#endif enddo end if end do @@ -548,6 +559,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) + liq_from_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 diff --git a/src/clm5/biogeophys/SoilHydrologyMod.F90 b/src/clm5/biogeophys/SoilHydrologyMod.F90 index d428fc8bd..c138230ba 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 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 @@ -2343,12 +2352,13 @@ end subroutine RenewCondensation #ifdef COUP_OAS_PFL !----------------------------------------------------------------------- subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & - num_urbanc, filter_urbanc, waterflux_inst) + num_urbanc, filter_urbanc, waterstate_inst, waterflux_inst) ! ! !DESCRIPTION: ! Calculate subsurface water fluxes which will be sent to ParFlow ! ! !USES: + use clm_time_manager , only : get_step_size use column_varcon , only : icol_road_perv use clm_varpar , only : nlevsoi ! @@ -2358,10 +2368,12 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points + type(waterstate_type) , intent(in) :: waterstate_inst type(waterflux_type) , intent(inout) :: waterflux_inst ! !LOCAL VARIABLES: character(len=32) :: subname = 'ParFlowDrainage' ! subroutine name integer :: c,j,fc ! indices + real(r8) :: dtime ! land model time step (sec) real(r8), parameter :: m_per_mm = 1.e-3_r8 ! 0.001 meters per mm real(r8), parameter :: sec_per_hr = 3600._r8 ! 3600 s in 1 hour @@ -2369,6 +2381,10 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & associate( & dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) + h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) + h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) + pfl_top_sync => waterstate_inst%pfl_top_sync_col , & ! Input: [real(r8) (:) ] liquid+ice of soil layer 1 after ParFlow exchange (kg/m2) + 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) @@ -2384,6 +2400,8 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & ! contributing fraction instead of the gridcell mean by c2g. qflx_parflow(bounds%begc:bounds%endc, 1:nlevsoi) = 0._r8 + dtime = get_step_size() + ! Calculate here the source/sink term for ParFlow do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) @@ -2396,11 +2414,22 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & qflx_parflow(c,j) = -qflx_rootsoi(c,j) !mm/s end if end do + + ! Take into account changes of 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) & + + (h2osoi_liq(c,1) + h2osoi_ice(c,1) - pfl_top_sync(c)) / dtime + ! 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)) diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index 81f86e436..f582230b7 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) + pfl_top_sync => waterstate_inst%pfl_top_sync_col , & ! Output: [real(r8) (:) ] liquid+ice of soil layer 1 after exchange (kg/m2) pfl_psi => waterstate_inst%pfl_psi_col & ! Input: [real(r8) (:,:) ] ParFlow pressure head (mm) ) ! end associate statement @@ -1517,11 +1518,14 @@ 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 + ! Remember the top layer mass after the exchange. Any later change by eCLM + ! (condensation, sublimation, snow layer combination) is passed to Parflow. + pfl_top_sync(c) = h2osoi_liq(c,1) + h2osoi_ice(c,1) + ! 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/WaterStateType.F90 b/src/clm5/biogeophys/WaterStateType.F90 index fdec7cf1b..5cf696574 100644 --- a/src/clm5/biogeophys/WaterStateType.F90 +++ b/src/clm5/biogeophys/WaterStateType.F90 @@ -108,6 +108,7 @@ module WaterstateType #ifdef COUP_OAS_PFL real(r8), pointer :: pfl_psi_col (:,:) ! ParFlow pressure head COUP_OAS_PFL real(r8), pointer :: pfl_h2osoi_liq_col (:,:) ! ParFlow soil liquid COUP_OAS_PFL + real(r8), pointer :: pfl_top_sync_col (:) ! liquid+ice of soil layer 1 right after ParFlow state was applied [kg/m2] #endif real(r8), pointer :: total_plant_stored_h2o_col(:)! col water that is bound in plants, including roots, sapwood, leaves, etc ! in most cases, the vegetation scheme does not have a dynamic @@ -247,6 +248,7 @@ subroutine InitAllocate(this, bounds) #ifdef COUP_OAS_PFL allocate(this%pfl_psi_col (begc:endc,1:nlevgrnd)) ; this%pfl_psi_col (:,:) = nan allocate(this%pfl_h2osoi_liq_col (begc:endc,1:nlevgrnd)) ; this%pfl_h2osoi_liq_col (:,:) = nan + allocate(this%pfl_top_sync_col (begc:endc)) ; this%pfl_top_sync_col (:) = nan #endif allocate(this%h2osoi_ice_tot_col (begc:endc)) ; this%h2osoi_ice_tot_col (:) = nan allocate(this%h2osoi_liq_tot_col (begc:endc)) ; this%h2osoi_liq_tot_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 From 711c0067d3c4bfcddb70d9c9581f9218ea6e5b58 Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Tue, 29 Sep 2026 20:54:36 +0200 Subject: [PATCH 3/4] CLM-ParFlow: replace mass difference by accumulated flux qflx_pfl_top - replace pfl_top_sync by the flux accumulator qflx_pfl_top_col in waterflux_inst - RenewCondensation: dew is passed to ParFlow via qflx_pfl_top instead of added to h2osoi_liq; the soil-ice change from dew/sublimation is added to qflx_pfl_top - CombineSnowLayers: when snow water goes to soil layer, it is passed to ParFlow via qflx_pfl_top instead of added to h2osoi_liq - ParFlowDrainage adds qflx_pfl_top to qflx_parflow --- src/clm5/biogeophys/HydrologyDrainageMod.F90 | 2 +- src/clm5/biogeophys/SnowHydrologyMod.F90 | 21 ++++++++++++ src/clm5/biogeophys/SoilHydrologyMod.F90 | 35 ++++++++++++-------- src/clm5/biogeophys/SoilWaterMovementMod.F90 | 7 ++-- src/clm5/biogeophys/WaterStateType.F90 | 2 -- src/clm5/biogeophys/WaterfluxType.F90 | 2 ++ 6 files changed, 48 insertions(+), 21 deletions(-) diff --git a/src/clm5/biogeophys/HydrologyDrainageMod.F90 b/src/clm5/biogeophys/HydrologyDrainageMod.F90 index b6c8a5f3d..e332aede5 100644 --- a/src/clm5/biogeophys/HydrologyDrainageMod.F90 +++ b/src/clm5/biogeophys/HydrologyDrainageMod.F90 @@ -135,7 +135,7 @@ subroutine HydrologyDrainage(bounds, & ! clm3.5/bld/usr.src/SoilHydrologyMod.F90 ! ignore drainage calculations in eCLM and instead pass these fluxes to ParFlow call ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & - num_urbanc, filter_urbanc, waterstate_inst, waterflux_inst) + num_urbanc, filter_urbanc, waterflux_inst) else if (use_aquifer_layer()) then call Drainage(bounds, num_hydrologyc, filter_hydrologyc, & diff --git a/src/clm5/biogeophys/SnowHydrologyMod.F90 b/src/clm5/biogeophys/SnowHydrologyMod.F90 index 594f4b93f..493487ea1 100644 --- a/src/clm5/biogeophys/SnowHydrologyMod.F90 +++ b/src/clm5/biogeophys/SnowHydrologyMod.F90 @@ -845,6 +845,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) @@ -888,7 +891,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 +#endif h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j) +#ifdef COUP_OAS_PFL + end if +#endif h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j) if (j == 0) then @@ -1010,7 +1022,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 +#endif h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c) +#ifdef COUP_OAS_PFL + end if +#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 c138230ba..3b34c6e81 100644 --- a/src/clm5/biogeophys/SoilHydrologyMod.F90 +++ b/src/clm5/biogeophys/SoilHydrologyMod.F90 @@ -2288,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( & @@ -2297,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) [+] ) @@ -2312,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 @@ -2320,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 @@ -2352,13 +2368,12 @@ end subroutine RenewCondensation #ifdef COUP_OAS_PFL !----------------------------------------------------------------------- subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & - num_urbanc, filter_urbanc, waterstate_inst, waterflux_inst) + num_urbanc, filter_urbanc, waterflux_inst) ! ! !DESCRIPTION: ! Calculate subsurface water fluxes which will be sent to ParFlow ! ! !USES: - use clm_time_manager , only : get_step_size use column_varcon , only : icol_road_perv use clm_varpar , only : nlevsoi ! @@ -2368,12 +2383,10 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points - type(waterstate_type) , intent(in) :: waterstate_inst type(waterflux_type) , intent(inout) :: waterflux_inst ! !LOCAL VARIABLES: character(len=32) :: subname = 'ParFlowDrainage' ! subroutine name integer :: c,j,fc ! indices - real(r8) :: dtime ! land model time step (sec) real(r8), parameter :: m_per_mm = 1.e-3_r8 ! 0.001 meters per mm real(r8), parameter :: sec_per_hr = 3600._r8 ! 3600 s in 1 hour @@ -2381,9 +2394,7 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & associate( & dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) - h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) - h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) - pfl_top_sync => waterstate_inst%pfl_top_sync_col , & ! Input: [real(r8) (:) ] liquid+ice of soil layer 1 after ParFlow exchange (kg/m2) + 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) @@ -2400,8 +2411,6 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & ! contributing fraction instead of the gridcell mean by c2g. qflx_parflow(bounds%begc:bounds%endc, 1:nlevsoi) = 0._r8 - dtime = get_step_size() - ! Calculate here the source/sink term for ParFlow do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) @@ -2415,11 +2424,9 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, & end if end do - ! Take into account changes of 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) & - + (h2osoi_liq(c,1) + h2osoi_ice(c,1) - pfl_top_sync(c)) / dtime + ! 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 diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index f582230b7..5feed75bd 100644 --- a/src/clm5/biogeophys/SoilWaterMovementMod.F90 +++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90 @@ -1503,7 +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) - pfl_top_sync => waterstate_inst%pfl_top_sync_col , & ! Output: [real(r8) (:) ] liquid+ice of soil layer 1 after exchange (kg/m2) + 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 @@ -1522,9 +1522,8 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & smp_l(c,j) = max(smpmin(c), min(0._r8, pfl_psi(c,j))) end do - ! Remember the top layer mass after the exchange. Any later change by eCLM - ! (condensation, sublimation, snow layer combination) is passed to Parflow. - pfl_top_sync(c) = h2osoi_liq(c,1) + h2osoi_ice(c,1) + ! 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. diff --git a/src/clm5/biogeophys/WaterStateType.F90 b/src/clm5/biogeophys/WaterStateType.F90 index 5cf696574..fdec7cf1b 100644 --- a/src/clm5/biogeophys/WaterStateType.F90 +++ b/src/clm5/biogeophys/WaterStateType.F90 @@ -108,7 +108,6 @@ module WaterstateType #ifdef COUP_OAS_PFL real(r8), pointer :: pfl_psi_col (:,:) ! ParFlow pressure head COUP_OAS_PFL real(r8), pointer :: pfl_h2osoi_liq_col (:,:) ! ParFlow soil liquid COUP_OAS_PFL - real(r8), pointer :: pfl_top_sync_col (:) ! liquid+ice of soil layer 1 right after ParFlow state was applied [kg/m2] #endif real(r8), pointer :: total_plant_stored_h2o_col(:)! col water that is bound in plants, including roots, sapwood, leaves, etc ! in most cases, the vegetation scheme does not have a dynamic @@ -248,7 +247,6 @@ subroutine InitAllocate(this, bounds) #ifdef COUP_OAS_PFL allocate(this%pfl_psi_col (begc:endc,1:nlevgrnd)) ; this%pfl_psi_col (:,:) = nan allocate(this%pfl_h2osoi_liq_col (begc:endc,1:nlevgrnd)) ; this%pfl_h2osoi_liq_col (:,:) = nan - allocate(this%pfl_top_sync_col (begc:endc)) ; this%pfl_top_sync_col (:) = nan #endif allocate(this%h2osoi_ice_tot_col (begc:endc)) ; this%h2osoi_ice_tot_col (:) = nan allocate(this%h2osoi_liq_tot_col (begc:endc)) ; this%h2osoi_liq_tot_col (:) = nan 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 From 6a8d52bc189488c787a0aff68f795faf6354a180 Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Thu, 1 Oct 2026 12:39:37 +0200 Subject: [PATCH 4/4] eCLM-ParFlow: Increase readability -change preprocessor statements - rename liq_from_soil to snowliq_to_soil --- src/clm5/biogeophys/SnowHydrologyMod.F90 | 21 +++++++++++---------- src/clm5/biogeophys/SoilHydrologyMod.F90 | 4 ++-- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/src/clm5/biogeophys/SnowHydrologyMod.F90 b/src/clm5/biogeophys/SnowHydrologyMod.F90 index 493487ea1..2aa2df78a 100644 --- a/src/clm5/biogeophys/SnowHydrologyMod.F90 +++ b/src/clm5/biogeophys/SnowHydrologyMod.F90 @@ -288,7 +288,8 @@ subroutine SnowWater(bounds, & 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) :: liq_from_soil(bounds%begc:bounds%endc) ! liquid taken from soil layer 1 to fill snow liquid deficit [mm] + 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 !----------------------------------------------------------------------- @@ -337,7 +338,7 @@ subroutine SnowWater(bounds, & c = filter_snowc(fc) l=col%landunit(c) #ifdef COUP_OAS_PFL - liq_from_soil(c) = 0._r8 + snowliq_to_soil(c) = 0._r8 #endif wgdif = h2osoi_ice(c,snl(c)+1) & @@ -361,7 +362,7 @@ subroutine SnowWater(bounds, & #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) liq_from_soil(c) = wgdif + if (j == 0) snowliq_to_soil(c) = wgdif #endif enddo end if @@ -560,7 +561,7 @@ 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) + liq_from_soil(c) / dtime + 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 @@ -896,10 +897,10 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, & ! 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 -#endif - h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j) -#ifdef COUP_OAS_PFL + 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) @@ -1027,10 +1028,10 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, & ! Pass the snow water to ParFlow qflx_pfl_top(c) = qflx_pfl_top(c) + zwliq(c)/dtime else -#endif - h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c) -#ifdef COUP_OAS_PFL + 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 diff --git a/src/clm5/biogeophys/SoilHydrologyMod.F90 b/src/clm5/biogeophys/SoilHydrologyMod.F90 index 3b34c6e81..385788360 100644 --- a/src/clm5/biogeophys/SoilHydrologyMod.F90 +++ b/src/clm5/biogeophys/SoilHydrologyMod.F90 @@ -451,8 +451,8 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f endif #ifdef COUP_OAS_PFL - ! All surface input minus evaporation from soil and surface water is passed - ! to ParFlow, without the eCLM infiltration capacity. + ! 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)