From b9d73072888a780487b56a637502240d7342785a Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Thu, 11 Jun 2026 23:27:31 +0200 Subject: [PATCH 1/5] reduce log-file output avoid spamming in log file, thereby improve performance --- src/clm5/biogeophys/BalanceCheckMod.F90 | 4 ++++ src/clm5/main/lnd2atmMod.F90 | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/clm5/biogeophys/BalanceCheckMod.F90 b/src/clm5/biogeophys/BalanceCheckMod.F90 index 7ca412955a..a664391d99 100644 --- a/src/clm5/biogeophys/BalanceCheckMod.F90 +++ b/src/clm5/biogeophys/BalanceCheckMod.F90 @@ -356,6 +356,7 @@ subroutine BalanceCheck( bounds, & end do found = .false. +#ifndef COUP_OAS_PFL do c = bounds%begc, bounds%endc if (abs(errh2o(c)) > 1.e-9_r8) then found = .true. @@ -446,6 +447,7 @@ subroutine BalanceCheck( bounds, & #endif end if end if +#endif ! Snow balance check @@ -691,6 +693,7 @@ subroutine BalanceCheck( bounds, & ! Soil energy balance check +#ifndef COUP_OAS_PFL found = .false. do c = bounds%begc,bounds%endc if (col%active(c)) then @@ -714,6 +717,7 @@ subroutine BalanceCheck( bounds, & end if #endif end if +#endif end associate diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index 7d1d68602d..341235029c 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -472,7 +472,7 @@ subroutine lnd2atm(bounds, & do j = 1, nlevsoi if (lnd2atm_inst%qflx_parflow_grc(g,j) == spval) then lnd2atm_inst%qflx_parflow_grc(g,j) = 0._r8 - write(iulog,*)'WARNING: qflx_parflow_grc is nan at grid point ',g,' level',j,' replaced with 0.' +! write(iulog,*)'WARNING: qflx_parflow_grc is nan at grid point ',g,' level',j,' replaced with 0.' end if end do enddo From 1206290d670a07961e45d00a7de8a63c688f66a4 Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Tue, 11 Aug 2026 18:07:12 +0200 Subject: [PATCH 2/5] eCLM-ParFlow: account for soil ice when receiving ParFlow water content - ParFlow has no ice phase, thus pfl_h2osoi_liq is the total soil water. Cap the diagnosed ice content and take the remaining as liquid. --- src/clm5/biogeophys/SoilWaterMovementMod.F90 | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index f86212882a..2ab4fc48c7 100644 --- a/src/clm5/biogeophys/SoilWaterMovementMod.F90 +++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90 @@ -1487,21 +1487,24 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & integer :: c,j,fc ! indices associate(& - smp_l => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix potential [mm] - h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) + smp_l => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix potential [mm] + h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) + 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_psi => waterstate_inst%pfl_psi_col & ! Input: [real(r8) (:,:) ] ParFlow pressure head (mm) ) ! end associate statement - ! COUP_OAS_PFL - ! TODO + ! Exchange of soil water and pressure head between ParFlow and eCLM + ! ParFlow has no ice phase, so that the received water it the total water content. + ! Keep the ice what PhaseChanged diagnosed earlier in the time step. do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) do j = 1, nlevgrnd - h2osoi_liq(c,j) = pfl_h2osoi_liq(c,j) + 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 From 772b24f73f7d86ad2e8006d7d5570c36ef7598dc Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Tue, 18 Aug 2026 16:49:13 +0200 Subject: [PATCH 3/5] eCLM-ParFlow: Calculate hydraulic conductivity based on ParFlow's water content - bugfix: before hk_l was never set for soil water movement method 4, value was the initial value =0 all time - recompute hk_l from ParFlow water content similar to standalone eCLM - see details in #128 --- src/clm5/biogeophys/SoilWaterMovementMod.F90 | 35 +++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index 2ab4fc48c7..13a7fae132 100644 --- a/src/clm5/biogeophys/SoilWaterMovementMod.F90 +++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90 @@ -1454,7 +1454,7 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & use abortutils , only : endrun use decompMod , only : bounds_type use clm_varctl , only : iulog, use_hydrstress - use clm_varcon , only : denh2o, denice + use clm_varcon , only : denh2o, denice, e_ice use clm_varpar , only : nlevgrnd use clm_time_manager , only : get_step_size, get_nstep use SoilStateType , only : soilstate_type @@ -1485,8 +1485,18 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & class(soil_water_retention_curve_type), intent(in) :: soil_water_retention_curve integer :: c,j,fc ! indices + integer :: nlayers ! lower boundary index + real(r8) :: s1 ! "s" at interface of layer + real(r8) :: s2(1:nlevgrnd) ! "s" at layer node + real(r8) :: imped ! ice impedance + real(r8) :: hk ! hydraulic conductivity (mm/s) associate(& + dz => col%dz , & ! Input: [real(r8) (:,:) ] layer thickness (m) + nbedrock => col%nbedrock , & ! Input: [integer (:) ] index of shallowest bedrock layer + icefrac => soilhydrology_inst%icefrac_col , & ! Input: [real(r8) (:,:) ] fraction of ice + watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) + hk_l => soilstate_inst%hk_l_col , & ! Output: [real(r8) (:,:) ] hydraulic conductivity (mm/s) smp_l => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix potential [mm] h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) @@ -1509,6 +1519,29 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & smp_l(c,j) = max(smpmin(c), pfl_psi(c,j)) end if end do + + ! 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) + do j = 1, nlayers + s2(j) = max(h2osoi_liq(c,j), 1.0e-6_r8) & + / (dz(c,j) * denh2o * watsat(c,j)) + s2(j) = min(max(s2(j), 0.01_r8), 1._r8) + end do + do j = 1, nlayers + ! hk is evaluated at the layer interface, as in the standalone model + if (j == nlayers) then + s1 = s2(j) + call IceImpedance(icefrac(c,j), e_ice, imped) + else + s1 = 0.5_r8 * (s2(j) + s2(j+1)) + call IceImpedance(0.5_r8*(icefrac(c,j) + icefrac(c,j+1)), e_ice, imped) + end if + s1 = min(max(s1, 0.01_r8), 1._r8) + call soil_water_retention_curve%soil_hk(c, j, s1, imped, & + soilstate_inst, hk) + hk_l(c,j) = hk + end do end do end associate From 910d5d0b97e3733a3e808e64d328c5d5fe41b823 Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Tue, 18 Aug 2026 17:11:33 +0200 Subject: [PATCH 4/5] eCLM-ParFlow: Update ice fraction after capping soil ice - icefrac was last set in infiltration, so the ice impedance in hk_l used a stale value - refresh icefrac from the capped h2osoi_ice, as in Infiltration --- src/clm5/biogeophys/SoilWaterMovementMod.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index 13a7fae132..f119b28fb3 100644 --- a/src/clm5/biogeophys/SoilWaterMovementMod.F90 +++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90 @@ -1488,13 +1488,14 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & integer :: nlayers ! lower boundary index real(r8) :: s1 ! "s" at interface of layer real(r8) :: s2(1:nlevgrnd) ! "s" at layer node + real(r8) :: vol_ice ! partial volume of ice real(r8) :: imped ! ice impedance real(r8) :: hk ! hydraulic conductivity (mm/s) associate(& dz => col%dz , & ! Input: [real(r8) (:,:) ] layer thickness (m) nbedrock => col%nbedrock , & ! Input: [integer (:) ] index of shallowest bedrock layer - icefrac => soilhydrology_inst%icefrac_col , & ! Input: [real(r8) (:,:) ] fraction of ice + icefrac => soilhydrology_inst%icefrac_col , & ! Output: [real(r8) (:,:) ] fraction of ice watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) hk_l => soilstate_inst%hk_l_col , & ! Output: [real(r8) (:,:) ] hydraulic conductivity (mm/s) smp_l => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix potential [mm] @@ -1524,6 +1525,10 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & ! ParFlow water content using the same formulation as compute_hydraulic_properties. nlayers = nbedrock(c) do j = 1, nlayers + ! icefrac was last set in Infiltration, before the ice limiter above + ! refresh it so the impedance matches the ParFlow state + vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) + icefrac(c,j) = min(1._r8, vol_ice/watsat(c,j)) s2(j) = max(h2osoi_liq(c,j), 1.0e-6_r8) & / (dz(c,j) * denh2o * watsat(c,j)) s2(j) = min(max(s2(j), 0.01_r8), 1._r8) From 858513b54d32d3fa24ea56309688c63a74b23b88 Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Tue, 18 Aug 2026 17:14:28 +0200 Subject: [PATCH 5/5] Revert "reduce log-file output" This reverts commit b9d73072888a780487b56a637502240d7342785a. --- src/clm5/biogeophys/BalanceCheckMod.F90 | 4 ---- src/clm5/main/lnd2atmMod.F90 | 2 +- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/clm5/biogeophys/BalanceCheckMod.F90 b/src/clm5/biogeophys/BalanceCheckMod.F90 index a664391d99..7ca412955a 100644 --- a/src/clm5/biogeophys/BalanceCheckMod.F90 +++ b/src/clm5/biogeophys/BalanceCheckMod.F90 @@ -356,7 +356,6 @@ subroutine BalanceCheck( bounds, & end do found = .false. -#ifndef COUP_OAS_PFL do c = bounds%begc, bounds%endc if (abs(errh2o(c)) > 1.e-9_r8) then found = .true. @@ -447,7 +446,6 @@ subroutine BalanceCheck( bounds, & #endif end if end if -#endif ! Snow balance check @@ -693,7 +691,6 @@ subroutine BalanceCheck( bounds, & ! Soil energy balance check -#ifndef COUP_OAS_PFL found = .false. do c = bounds%begc,bounds%endc if (col%active(c)) then @@ -717,7 +714,6 @@ subroutine BalanceCheck( bounds, & end if #endif end if -#endif end associate diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index 341235029c..7d1d68602d 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -472,7 +472,7 @@ subroutine lnd2atm(bounds, & do j = 1, nlevsoi if (lnd2atm_inst%qflx_parflow_grc(g,j) == spval) then lnd2atm_inst%qflx_parflow_grc(g,j) = 0._r8 -! write(iulog,*)'WARNING: qflx_parflow_grc is nan at grid point ',g,' level',j,' replaced with 0.' + write(iulog,*)'WARNING: qflx_parflow_grc is nan at grid point ',g,' level',j,' replaced with 0.' end if end do enddo