From 22f6de1986fa6a136cf535f2ec5b029c34489f2f Mon Sep 17 00:00:00 2001 From: kvrigor Date: Sun, 14 May 2023 16:59:49 +0200 Subject: [PATCH 1/9] Added coupling field ECLM_ICE_FRAC --- src/clm5/main/clm_driver.F90 | 4 ++-- src/clm5/main/lnd2atmMod.F90 | 22 ++++++++++++++++++---- src/clm5/main/lnd2atmType.F90 | 7 +++++++ src/clm5/oasis3/oas_defineMod.F90 | 5 +++-- src/clm5/oasis3/oas_sendReceiveMod.F90 | 1 + src/clm5/oasis3/oas_vardefMod.F90 | 2 +- 6 files changed, 32 insertions(+), 9 deletions(-) diff --git a/src/clm5/main/clm_driver.F90 b/src/clm5/main/clm_driver.F90 index befbba7a45..7959d610f2 100644 --- a/src/clm5/main/clm_driver.F90 +++ b/src/clm5/main/clm_driver.F90 @@ -1065,8 +1065,8 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call lnd2atm(bounds_proc, & atm2lnd_inst, surfalb_inst, temperature_inst, frictionvel_inst, & waterstate_inst, waterflux_inst, irrigation_inst, energyflux_inst, & - solarabs_inst, drydepvel_inst, & - vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, & + solarabs_inst, drydepvel_inst, vocemis_inst, fireemis_inst, & + dust_inst, ch4_inst, glc_behavior, soilhydrology_inst, & lnd2atm_inst, & #ifdef USE_PDAF soilhydrology_inst, soilstate_inst, & diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index 7bd00d126d..21caedb5dd 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -36,6 +36,7 @@ module lnd2atmMod use WaterstateType , only : waterstate_type use IrrigationMod , only : irrigation_type use glcBehaviorMod , only : glc_behavior_type + use SoilHydrologyType , only : soilhydrology_type use glc2lndMod , only : glc2lnd_type use ColumnType , only : col use LandunitType , only : lun @@ -135,10 +136,7 @@ subroutine lnd2atm(bounds, & waterstate_inst, waterflux_inst, irrigation_inst, energyflux_inst, & solarabs_inst, drydepvel_inst, & vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, & - lnd2atm_inst, & -#ifdef USE_PDAF - soilhydrology_inst, soilstate_inst, & -#endif + soilhydrology_inst, lnd2atm_inst, & net_carbon_exchange_grc) ! ! !DESCRIPTION: @@ -164,6 +162,7 @@ subroutine lnd2atm(bounds, & type(dust_type) , intent(in) :: dust_inst type(ch4_type) , intent(in) :: ch4_inst type(glc_behavior_type) , intent(in) :: glc_behavior + type(soilhydrology_type) , intent(in) :: soilhydrology_inst type(lnd2atm_type) , intent(inout) :: lnd2atm_inst #ifdef USE_PDAF ! Yorck @@ -578,6 +577,21 @@ subroutine lnd2atm(bounds, & averaging_var = 0 + call c2g( bounds, nlevsoi, & + soilhydrology_inst%icefrac_col (bounds%begc:bounds%endc, :), & + lnd2atm_inst%ice_frac_grc (bounds%begg:bounds%endg, :), & + c2l_scale_type= 'unity', l2g_scale_type='unity' ) + + do c = bounds%begc, bounds%endc + if (col%hydrologically_active(c)) then + if (col%itype(c) == istsoil .or. col%itype(c) == istcrop) then + g = col%gridcell(c) + do j = 1, nlevsoi + ! Convert eCLM fluxes (mm/s) to ParFlow fluxes (1/hr): + ! 1/hr = [mm/s] * [s/hr] * [m/mm] * [1/m] + lnd2atm_inst%qflx_parflow_grc(g,j) = lnd2atm_inst%qflx_parflow_grc(g,j) * sec_per_hr * m_per_mm * (1/col%dz(c,j)) + enddo + end if end if averaging_var = averaging_var+1 diff --git a/src/clm5/main/lnd2atmType.F90 b/src/clm5/main/lnd2atmType.F90 index 11e93930a1..9b6e3f1ff8 100644 --- a/src/clm5/main/lnd2atmType.F90 +++ b/src/clm5/main/lnd2atmType.F90 @@ -77,6 +77,7 @@ module lnd2atmType real(r8), pointer :: qflx_liq_from_ice_col (:) => null() ! liquid runoff from converted ice runoff #ifdef COUP_OAS_PFL real(r8), pointer :: qflx_parflow_grc (:,:) => null() ! source/sink flux per soil layer sent to ParFlow [1/hr] [- out from root] + real(r8), pointer :: ice_frac_grc (:,:) => null() ! soil ice fraction (values in [0,1]) #endif real(r8), pointer :: qirrig_grc (:) => null() ! irrigation flux @@ -191,6 +192,7 @@ subroutine InitAllocate(this, bounds) allocate(this%qflx_liq_from_ice_col(begc:endc)) ; this%qflx_liq_from_ice_col(:) =ival #ifdef COUP_OAS_PFL allocate(this%qflx_parflow_grc (begg:endg,1:nlevsoi)); this%qflx_parflow_grc (:,:) =ival + allocate(this%ice_frac_grc (begg:endg,1:nlevsoi)); this%ice_frac_grc (:,:) =ival #endif allocate(this%qirrig_grc (begg:endg)) ; this%qirrig_grc (:) =ival @@ -344,6 +346,11 @@ subroutine InitHistory(this, bounds) call hist_addfld2d (fname='QPARFLOW_TO_OASIS', units='1/hr', type2d='levsoi', & avgflag='A', long_name='source/sink flux per soil layer sent to ParFlow', & ptr_lnd=this%qflx_parflow_grc, default='inactive') + + this%ice_frac_grc(begg:endg, :) = 0._r8 + call hist_addfld2d (fname='FRAC_ICE', units='unitless', type2d='levsoi', & + avgflag='A', long_name='soil ice fraction sent to ParFlow', & + ptr_lnd=this%ice_frac_grc, default='inactive') #endif end subroutine InitHistory diff --git a/src/clm5/oasis3/oas_defineMod.F90 b/src/clm5/oasis3/oas_defineMod.F90 index cdae0fceeb..7f414fe1ff 100644 --- a/src/clm5/oasis3/oas_defineMod.F90 +++ b/src/clm5/oasis3/oas_defineMod.F90 @@ -102,8 +102,9 @@ subroutine oas_definitions_init(bounds) var_nodims(1) = 1 ! unused var_nodims(2) = nlevsoi ! number of fields in a bundle - call oasis_def_var(oas_et_loss_id, "ECLM_ET", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) - + call oasis_def_var(oas_et_loss_id, "ECLM_ET", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) + call oasis_def_var(oas_ice_frac_id, "ECLM_ICE_FRAC", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) + var_nodims(2) = nlevgrnd ! number of fields in a bundle call oasis_def_var(oas_sat_id, "ECLM_SOILLIQ", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) call oasis_def_var(oas_psi_id, "ECLM_PSI", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) diff --git a/src/clm5/oasis3/oas_sendReceiveMod.F90 b/src/clm5/oasis3/oas_sendReceiveMod.F90 index abd8db7db9..bc750d1604 100644 --- a/src/clm5/oasis3/oas_sendReceiveMod.F90 +++ b/src/clm5/oasis3/oas_sendReceiveMod.F90 @@ -48,6 +48,7 @@ subroutine oas_send(bounds, seconds_elapsed, lnd2atm_inst) integer :: info call oasis_put(oas_et_loss_id, seconds_elapsed, lnd2atm_inst%qflx_parflow_grc, info) + call oasis_put(oas_ice_frac_id, seconds_elapsed, lnd2atm_inst%ice_frac_grc, info) end subroutine oas_send #endif diff --git a/src/clm5/oasis3/oas_vardefMod.F90 b/src/clm5/oasis3/oas_vardefMod.F90 index 56ef8e0752..02bdee743d 100644 --- a/src/clm5/oasis3/oas_vardefMod.F90 +++ b/src/clm5/oasis3/oas_vardefMod.F90 @@ -3,7 +3,7 @@ module oas_vardefMod save #ifdef COUP_OAS_PFL - integer :: oas_psi_id, oas_et_loss_id, oas_sat_id + integer :: oas_psi_id, oas_et_loss_id, oas_sat_id, oas_ice_frac_id #endif #ifdef COUP_OAS_ICON From a08c9a7cfc81e9523f307a22fd2250d58c37e9f6 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Tue, 13 Jun 2023 17:31:00 +0200 Subject: [PATCH 2/9] Fixed incorrect dimensions of ice_grac_grc --- src/clm5/biogeophys/SoilHydrologyType.F90 | 1 + src/clm5/main/lnd2atmMod.F90 | 2 +- src/clm5/main/lnd2atmType.F90 | 4 ++-- src/clm5/oasis3/oas_defineMod.F90 | 4 ++-- 4 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/clm5/biogeophys/SoilHydrologyType.F90 b/src/clm5/biogeophys/SoilHydrologyType.F90 index c62ec465f1..033ca988a2 100644 --- a/src/clm5/biogeophys/SoilHydrologyType.F90 +++ b/src/clm5/biogeophys/SoilHydrologyType.F90 @@ -235,6 +235,7 @@ subroutine InitCold(this, bounds) ! averaging for the accum field do c = bounds%begc, bounds%endc this%num_substeps_col(c) = spval + this%icefrac_col(c,:) = spval end do end subroutine InitCold diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index 21caedb5dd..00b484563c 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -577,7 +577,7 @@ subroutine lnd2atm(bounds, & averaging_var = 0 - call c2g( bounds, nlevsoi, & + call c2g( bounds, nlevgrnd, & soilhydrology_inst%icefrac_col (bounds%begc:bounds%endc, :), & lnd2atm_inst%ice_frac_grc (bounds%begg:bounds%endg, :), & c2l_scale_type= 'unity', l2g_scale_type='unity' ) diff --git a/src/clm5/main/lnd2atmType.F90 b/src/clm5/main/lnd2atmType.F90 index 9b6e3f1ff8..de48694772 100644 --- a/src/clm5/main/lnd2atmType.F90 +++ b/src/clm5/main/lnd2atmType.F90 @@ -192,7 +192,7 @@ subroutine InitAllocate(this, bounds) allocate(this%qflx_liq_from_ice_col(begc:endc)) ; this%qflx_liq_from_ice_col(:) =ival #ifdef COUP_OAS_PFL allocate(this%qflx_parflow_grc (begg:endg,1:nlevsoi)); this%qflx_parflow_grc (:,:) =ival - allocate(this%ice_frac_grc (begg:endg,1:nlevsoi)); this%ice_frac_grc (:,:) =ival + allocate(this%ice_frac_grc (begg:endg,1:nlevgrnd)); this%ice_frac_grc (:,:) =ival #endif allocate(this%qirrig_grc (begg:endg)) ; this%qirrig_grc (:) =ival @@ -348,7 +348,7 @@ subroutine InitHistory(this, bounds) ptr_lnd=this%qflx_parflow_grc, default='inactive') this%ice_frac_grc(begg:endg, :) = 0._r8 - call hist_addfld2d (fname='FRAC_ICE', units='unitless', type2d='levsoi', & + call hist_addfld2d (fname='FRAC_ICE', units='unitless', type2d='levgrnd', & avgflag='A', long_name='soil ice fraction sent to ParFlow', & ptr_lnd=this%ice_frac_grc, default='inactive') #endif diff --git a/src/clm5/oasis3/oas_defineMod.F90 b/src/clm5/oasis3/oas_defineMod.F90 index 7f414fe1ff..f0d07bd548 100644 --- a/src/clm5/oasis3/oas_defineMod.F90 +++ b/src/clm5/oasis3/oas_defineMod.F90 @@ -103,12 +103,12 @@ subroutine oas_definitions_init(bounds) var_nodims(2) = nlevsoi ! number of fields in a bundle call oasis_def_var(oas_et_loss_id, "ECLM_ET", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) - call oasis_def_var(oas_ice_frac_id, "ECLM_ICE_FRAC", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) var_nodims(2) = nlevgrnd ! number of fields in a bundle call oasis_def_var(oas_sat_id, "ECLM_SOILLIQ", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) call oasis_def_var(oas_psi_id, "ECLM_PSI", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) -#endif + call oasis_def_var(oas_ice_frac_id, "ECLM_ICE_FRAC", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) +#endif #ifdef COUP_OAS_ICON From fa01ebbc14d7f2ec2f16f7cf2524d9d0e74e4349 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Tue, 29 Aug 2023 10:45:53 +0200 Subject: [PATCH 3/9] Sent ice impedance to ParFlow instead of ice fraction --- src/clm5/biogeophys/SoilHydrologyType.F90 | 9 +++++++ src/clm5/biogeophys/SoilWaterMovementMod.F90 | 27 +++++++++++--------- src/clm5/main/lnd2atmMod.F90 | 4 +-- src/clm5/main/lnd2atmType.F90 | 12 ++++----- src/clm5/oasis3/oas_defineMod.F90 | 2 +- src/clm5/oasis3/oas_sendReceiveMod.F90 | 2 +- src/clm5/oasis3/oas_vardefMod.F90 | 2 +- 7 files changed, 35 insertions(+), 23 deletions(-) diff --git a/src/clm5/biogeophys/SoilHydrologyType.F90 b/src/clm5/biogeophys/SoilHydrologyType.F90 index 033ca988a2..088b4b4b35 100644 --- a/src/clm5/biogeophys/SoilHydrologyType.F90 +++ b/src/clm5/biogeophys/SoilHydrologyType.F90 @@ -52,6 +52,9 @@ Module SoilHydrologyType real(r8), pointer :: max_infil_col (:) ! col VIC maximum infiltration rate calculated in VIC real(r8), pointer :: i_0_col (:) ! col VIC average saturation in top soil layers real(r8), pointer :: ice_col (:,:) ! col VIC soil ice (kg/m2) for VIC soil layers +#ifdef COUP_OAS_PFL + real(r8), pointer :: ice_impedance_col (:,:) ! col ice impdeance +#endif #ifdef USE_PDAF ! Yorck @@ -130,6 +133,9 @@ subroutine InitAllocate(this, bounds) allocate(this%qcharge_col (begc:endc)) ; this%qcharge_col (:) = nan allocate(this%fracice_col (begc:endc,nlevgrnd)) ; this%fracice_col (:,:) = nan allocate(this%icefrac_col (begc:endc,nlevgrnd)) ; this%icefrac_col (:,:) = nan +#ifdef COUP_OAS_PFL + allocate(this%ice_impedance_col (begc:endc,nlevgrnd)) ; this%ice_impedance_col (:,:) = nan +#endif allocate(this%fcov_col (begc:endc)) ; this%fcov_col (:) = nan allocate(this%fsat_col (begc:endc)) ; this%fsat_col (:) = nan allocate(this%h2osfc_thresh_col (begc:endc)) ; this%h2osfc_thresh_col (:) = nan @@ -236,6 +242,9 @@ subroutine InitCold(this, bounds) do c = bounds%begc, bounds%endc this%num_substeps_col(c) = spval this%icefrac_col(c,:) = spval +#ifdef COUP_OAS_PFL + this%ice_impedance_col(c,:) = 1.0_r8 +#endif end do end subroutine InitCold diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index 81f86e4366..1b58c7b36c 100644 --- a/src/clm5/biogeophys/SoilWaterMovementMod.F90 +++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90 @@ -1455,7 +1455,7 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & use decompMod , only : bounds_type use clm_varctl , only : iulog, use_hydrstress use clm_varcon , only : denh2o, denice, e_ice - use clm_varpar , only : nlevgrnd + use clm_varpar , only : nlevsoi, nlevgrnd use clm_time_manager , only : get_step_size, get_nstep use SoilStateType , only : soilstate_type use SoilHydrologyType , only : soilhydrology_type @@ -1493,17 +1493,12 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & 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 , & ! 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] - 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) + 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) + 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) + icefrac => soilhydrology_inst%icefrac_col , & ! Input: [real(r8) (:,:) ] fraction of ice + ice_impedance => soilhydrology_inst%ice_impedance_col & ! Input: [real(r8) (:,:) ] ice impedance ) ! end associate statement @@ -1520,6 +1515,14 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & if (pfl_psi(c,j) <= 0) then smp_l(c,j) = max(smpmin(c), pfl_psi(c,j)) end if + + if(j==nlevsoi)then + call IceImpedance(icefrac(c,j), e_ice, ice_impedance(c,j) ) + else if(j null() ! liquid runoff from converted ice runoff #ifdef COUP_OAS_PFL real(r8), pointer :: qflx_parflow_grc (:,:) => null() ! source/sink flux per soil layer sent to ParFlow [1/hr] [- out from root] - real(r8), pointer :: ice_frac_grc (:,:) => null() ! soil ice fraction (values in [0,1]) + real(r8), pointer :: ice_impedance_grc (:,:) => null() ! soil ice impedance (values in [0,1]) #endif real(r8), pointer :: qirrig_grc (:) => null() ! irrigation flux @@ -192,7 +192,7 @@ subroutine InitAllocate(this, bounds) allocate(this%qflx_liq_from_ice_col(begc:endc)) ; this%qflx_liq_from_ice_col(:) =ival #ifdef COUP_OAS_PFL allocate(this%qflx_parflow_grc (begg:endg,1:nlevsoi)); this%qflx_parflow_grc (:,:) =ival - allocate(this%ice_frac_grc (begg:endg,1:nlevgrnd)); this%ice_frac_grc (:,:) =ival + allocate(this%ice_impedance_grc (begg:endg,1:nlevgrnd)); this%ice_impedance_grc (:,:) =1.0_r8 #endif allocate(this%qirrig_grc (begg:endg)) ; this%qirrig_grc (:) =ival @@ -347,10 +347,10 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='source/sink flux per soil layer sent to ParFlow', & ptr_lnd=this%qflx_parflow_grc, default='inactive') - this%ice_frac_grc(begg:endg, :) = 0._r8 - call hist_addfld2d (fname='FRAC_ICE', units='unitless', type2d='levgrnd', & - avgflag='A', long_name='soil ice fraction sent to ParFlow', & - ptr_lnd=this%ice_frac_grc, default='inactive') + this%ice_impedance_grc(begg:endg, :) = 1._r8 + call hist_addfld2d (fname='ICE_IMPEDANCE', units='unitless', type2d='levgrnd', & + avgflag='A', long_name='soil ice impedance sent to ParFlow', & + ptr_lnd=this%ice_impedance_grc, default='inactive') #endif end subroutine InitHistory diff --git a/src/clm5/oasis3/oas_defineMod.F90 b/src/clm5/oasis3/oas_defineMod.F90 index f0d07bd548..db980ed608 100644 --- a/src/clm5/oasis3/oas_defineMod.F90 +++ b/src/clm5/oasis3/oas_defineMod.F90 @@ -107,7 +107,7 @@ subroutine oas_definitions_init(bounds) var_nodims(2) = nlevgrnd ! number of fields in a bundle call oasis_def_var(oas_sat_id, "ECLM_SOILLIQ", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) call oasis_def_var(oas_psi_id, "ECLM_PSI", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) - call oasis_def_var(oas_ice_frac_id, "ECLM_ICE_FRAC", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) + call oasis_def_var(oas_ice_impedance_id, "ECLM_ICE_IMPEDANCE", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) #endif #ifdef COUP_OAS_ICON diff --git a/src/clm5/oasis3/oas_sendReceiveMod.F90 b/src/clm5/oasis3/oas_sendReceiveMod.F90 index bc750d1604..ccff123aef 100644 --- a/src/clm5/oasis3/oas_sendReceiveMod.F90 +++ b/src/clm5/oasis3/oas_sendReceiveMod.F90 @@ -48,7 +48,7 @@ subroutine oas_send(bounds, seconds_elapsed, lnd2atm_inst) integer :: info call oasis_put(oas_et_loss_id, seconds_elapsed, lnd2atm_inst%qflx_parflow_grc, info) - call oasis_put(oas_ice_frac_id, seconds_elapsed, lnd2atm_inst%ice_frac_grc, info) + call oasis_put(oas_ice_impedance_id, seconds_elapsed, lnd2atm_inst%ice_impedance_grc, info) end subroutine oas_send #endif diff --git a/src/clm5/oasis3/oas_vardefMod.F90 b/src/clm5/oasis3/oas_vardefMod.F90 index 02bdee743d..033a6f1a65 100644 --- a/src/clm5/oasis3/oas_vardefMod.F90 +++ b/src/clm5/oasis3/oas_vardefMod.F90 @@ -3,7 +3,7 @@ module oas_vardefMod save #ifdef COUP_OAS_PFL - integer :: oas_psi_id, oas_et_loss_id, oas_sat_id, oas_ice_frac_id + integer :: oas_psi_id, oas_et_loss_id, oas_sat_id, oas_ice_impedance_id #endif #ifdef COUP_OAS_ICON From 8a56bf0c470301be2b8ce0f230acc29ee82f716a Mon Sep 17 00:00:00 2001 From: kvrigor Date: Fri, 4 Sep 2026 13:50:40 +0200 Subject: [PATCH 4/9] Fixed incorrect auto rebase conflict resolution at lnd2atmMod.F90 --- src/clm5/main/lnd2atmMod.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index 12fadeba49..1b849829df 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -137,6 +137,9 @@ subroutine lnd2atm(bounds, & solarabs_inst, drydepvel_inst, & vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, & soilhydrology_inst, lnd2atm_inst, & +#ifdef USE_PDAF + soilstate_inst, & +#endif net_carbon_exchange_grc) ! ! !DESCRIPTION: From f0775a2c8589a6fdc816eb6646bd70fc3630ded4 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Fri, 4 Sep 2026 13:57:27 +0200 Subject: [PATCH 5/9] Resolved rebase conflicts with 55db76e https://github.com/HPSCTerrSys/eCLM/commit/55db76e4bd97a4c52f2380a6786dca3b7c87c7c5 --- src/clm5/biogeophys/SoilWaterMovementMod.F90 | 33 +++++++++----------- 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index 1b58c7b36c..c8b14d6f3d 100644 --- a/src/clm5/biogeophys/SoilWaterMovementMod.F90 +++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90 @@ -1489,16 +1489,21 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & 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(& - 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) - 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) - icefrac => soilhydrology_inst%icefrac_col , & ! Input: [real(r8) (:,:) ] fraction of ice - ice_impedance => soilhydrology_inst%ice_impedance_col & ! Input: [real(r8) (:,:) ] ice impedance + dz => col%dz , & ! Input: [real(r8) (:,:) ] layer thickness (m) + nbedrock => col%nbedrock , & ! Input: [integer (:) ] index of shallowest bedrock layer + icefrac => soilhydrology_inst%icefrac_col , & ! Output: [real(r8) (:,:) ] fraction of ice + ice_impedance => soilhydrology_inst%ice_impedance_col, & ! Input: [real(r8) (:,:) ] ice impedance + 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) + 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 @@ -1515,14 +1520,6 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & if (pfl_psi(c,j) <= 0) then smp_l(c,j) = max(smpmin(c), pfl_psi(c,j)) end if - - if(j==nlevsoi)then - call IceImpedance(icefrac(c,j), e_ice, ice_impedance(c,j) ) - else if(j Date: Mon, 21 Sep 2026 11:28:34 +0200 Subject: [PATCH 6/9] eCLM-ParFlow: couple effective porosity to ParFlow - eCLM manages ParFlow's porosity - Send watsat - vol_ice to ParFlow as ECLM_EFFPOROSITY - Treat water in ParFlow water as liquid; subtract ice formed this time step, which ParFlow applies one coupling interval later - Average porosity per gridcell over hydrologically active columns only; set not coupled cells to -9999 --- src/clm5/biogeophys/SoilHydrologyType.F90 | 6 +++ src/clm5/biogeophys/SoilWaterMovementMod.F90 | 18 ++++--- src/clm5/biogeophys/WaterStateType.F90 | 3 ++ src/clm5/main/clm_driver.F90 | 7 ++- src/clm5/main/lnd2atmMod.F90 | 49 ++++++++++++++++++-- src/clm5/main/lnd2atmType.F90 | 6 +++ src/clm5/oasis3/oas_defineMod.F90 | 1 + src/clm5/oasis3/oas_sendReceiveMod.F90 | 3 +- src/clm5/oasis3/oas_vardefMod.F90 | 2 +- 9 files changed, 80 insertions(+), 15 deletions(-) diff --git a/src/clm5/biogeophys/SoilHydrologyType.F90 b/src/clm5/biogeophys/SoilHydrologyType.F90 index c62ec465f1..60b7e4ba29 100644 --- a/src/clm5/biogeophys/SoilHydrologyType.F90 +++ b/src/clm5/biogeophys/SoilHydrologyType.F90 @@ -29,6 +29,9 @@ Module SoilHydrologyType real(r8), pointer :: qcharge_col (:) ! col aquifer recharge rate (mm/s) real(r8), pointer :: fracice_col (:,:) ! col fractional impermeability (-) real(r8), pointer :: icefrac_col (:,:) ! col fraction of ice +#ifdef COUP_OAS_PFL + real(r8), pointer :: pfl_eff_porosity_col (:,:) ! col effective porosity sent to ParFlow = watsat - vol_ice (m3/m3), unfloored +#endif real(r8), pointer :: fcov_col (:) ! col fractional impermeable area real(r8), pointer :: fsat_col (:) ! col fractional area with water table at surface real(r8), pointer :: h2osfc_thresh_col (:) ! col level at which h2osfc "percolates" (time constant) @@ -130,6 +133,9 @@ subroutine InitAllocate(this, bounds) allocate(this%qcharge_col (begc:endc)) ; this%qcharge_col (:) = nan allocate(this%fracice_col (begc:endc,nlevgrnd)) ; this%fracice_col (:,:) = nan allocate(this%icefrac_col (begc:endc,nlevgrnd)) ; this%icefrac_col (:,:) = nan +#ifdef COUP_OAS_PFL + allocate(this%pfl_eff_porosity_col (begc:endc,nlevgrnd)) ; this%pfl_eff_porosity_col(:,:) = 0._r8 +#endif allocate(this%fcov_col (begc:endc)) ; this%fcov_col (:) = nan allocate(this%fsat_col (begc:endc)) ; this%fsat_col (:) = nan allocate(this%h2osfc_thresh_col (begc:endc)) ; this%h2osfc_thresh_col (:) = nan diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index 81f86e4366..17d1f25ee5 100644 --- a/src/clm5/biogeophys/SoilWaterMovementMod.F90 +++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90 @@ -1503,20 +1503,26 @@ 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_psi => waterstate_inst%pfl_psi_col & ! Input: [real(r8) (:,:) ] ParFlow pressure head (mm) + pfl_psi => waterstate_inst%pfl_psi_col , & ! Input: [real(r8) (:,:) ] ParFlow pressure head (mm) + h2osoi_ice_prev => waterstate_inst%h2osoi_ice_prev_col, & ! Input: [real(r8) (:,:) ] ice lens before PhaseChange (kg/m2) + pfl_eff_porosity => soilhydrology_inst%pfl_eff_porosity_col & ! Output: [real(r8) (:,:) ] effective porosity sent to ParFlow (m3/m3) ) ! end associate statement ! 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 and - ! attribute the remainder to liquid, so that liq and ice matches the received total. + ! ParFlow carries the liquid phase only, so the received water is the liquid + ! water content and h2osoi_ice stays as PhaseChange left it. ParFlow applies the + ! freeze one coupling interval later, so subtract it here to keep eCLM in sync. do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) 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)) + ! PhaseChange bounds the ice by the available water. Capping the state + ! keeps the frozen mass ParFlow derives from it exact. + h2osoi_ice(c,j) = min(h2osoi_ice(c,j), watsat(c,j)*dz(c,j)*denice) + h2osoi_liq(c,j) = max(0._r8, pfl_h2osoi_liq(c,j) & + - (h2osoi_ice(c,j) - h2osoi_ice_prev(c,j))) + pfl_eff_porosity(c,j) = watsat(c,j) - h2osoi_ice(c,j)/(dz(c,j)*denice) if (pfl_psi(c,j) <= 0) then smp_l(c,j) = max(smpmin(c), pfl_psi(c,j)) end if diff --git a/src/clm5/biogeophys/WaterStateType.F90 b/src/clm5/biogeophys/WaterStateType.F90 index fdec7cf1b3..c09f251170 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 :: h2osoi_ice_prev_col (:,:) ! ice lens at the start of the time step, before PhaseChange (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%h2osoi_ice_prev_col (begc:endc,1:nlevgrnd)) ; this%h2osoi_ice_prev_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 @@ -861,6 +863,7 @@ subroutine InitCold(this, bounds, & #ifdef COUP_OAS_PFL this%pfl_psi_col(bounds%begc:bounds%endc, 1:) = -1000._r8 this%pfl_h2osoi_liq_col(bounds%begc:bounds%endc, 1:) = spval + this%h2osoi_ice_prev_col(bounds%begc:bounds%endc, 1:) = 0._r8 #endif this%h2osoi_vol_prs_grc(bounds%begg:bounds%endg, 1:) = spval this%h2osoi_liq_col(bounds%begc:bounds%endc,-nlevsno+1:) = spval diff --git a/src/clm5/main/clm_driver.F90 b/src/clm5/main/clm_driver.F90 index befbba7a45..b43dd5cfda 100644 --- a/src/clm5/main/clm_driver.F90 +++ b/src/clm5/main/clm_driver.F90 @@ -1068,7 +1068,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro solarabs_inst, drydepvel_inst, & vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, & lnd2atm_inst, & -#ifdef USE_PDAF +#if defined(USE_PDAF) || defined(COUP_OAS_PFL) soilhydrology_inst, soilstate_inst, & #endif net_carbon_exchange_grc = net_carbon_exchange_grc(bounds_proc%begg:bounds_proc%endg)) @@ -1219,7 +1219,7 @@ subroutine clm_drv_init(bounds, & ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - use clm_varpar , only : nlevsno, nlevsoi + use clm_varpar , only : nlevsno, nlevsoi, nlevgrnd use CanopyStateType , only : canopystate_type use WaterStateType , only : waterstate_type use WaterFluxType , only : waterflux_type @@ -1257,6 +1257,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 + h2osoi_ice_prev => waterstate_inst%h2osoi_ice_prev_col , & ! Output: [real(r8) (:,:) ] COUP_OAS_PFL #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 @@ -1314,6 +1315,8 @@ subroutine clm_drv_init(bounds, & pfl_psi(c,:) = atm2lnd_inst%pfl_psi_grc(g,:) pfl_h2osoi_liq(c,:) = atm2lnd_inst%pfl_h2osoi_liq_grc(g,:) end if + ! Get ice content before PhaseChange, to make it consistent with ParFlow + h2osoi_ice_prev(c,1:nlevgrnd) = h2osoi_ice(c,1:nlevgrnd) end do #endif end associate diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index 7bd00d126d..82f1ea5114 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -43,9 +43,11 @@ module lnd2atmMod use landunit_varcon , only : istice_mec, istsoil, istcrop #ifdef USE_PDAF use clm_time_manager , only : get_nstep + use PatchType , only : patch +#endif +#if defined(USE_PDAF) || defined(COUP_OAS_PFL) use SoilHydrologyType , only : soilhydrology_type use SoilStateType , only : soilstate_type - use PatchType , only : patch #endif ! ! !PUBLIC TYPES: @@ -60,6 +62,10 @@ module lnd2atmMod ! !PRIVATE MEMBER FUNCTIONS: private :: handle_ice_runoff +#ifdef COUP_OAS_PFL + real(r8), parameter, private :: pfl_uncoupled = -9999._r8 +#endif + character(len=*), parameter, private :: sourcefile = & __FILE__ !------------------------------------------------------------------------ @@ -136,7 +142,7 @@ subroutine lnd2atm(bounds, & solarabs_inst, drydepvel_inst, & vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, & lnd2atm_inst, & -#ifdef USE_PDAF +#if defined(USE_PDAF) || defined(COUP_OAS_PFL) soilhydrology_inst, soilstate_inst, & #endif net_carbon_exchange_grc) @@ -165,16 +171,17 @@ subroutine lnd2atm(bounds, & type(ch4_type) , intent(in) :: ch4_inst type(glc_behavior_type) , intent(in) :: glc_behavior type(lnd2atm_type) , intent(inout) :: lnd2atm_inst -#ifdef USE_PDAF - ! Yorck +#if defined(USE_PDAF) || defined(COUP_OAS_PFL) type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(soilstate_type) , intent(inout) :: soilstate_inst - ! end Yorck #endif real(r8) , intent(in) :: net_carbon_exchange_grc( bounds%begg: ) ! net carbon exchange between land and atmosphere, positive for source (gC/m2/s) ! ! !LOCAL VARIABLES: integer :: c, g, j ! indices +#ifdef COUP_OAS_PFL + real(r8) :: pfl_eff_porosity_wt(bounds%begg:bounds%endg) ! summed weight of the hydrologically active columns +#endif #ifdef USE_PDAF integer :: p, l, index, counter ! indices #endif @@ -466,6 +473,38 @@ subroutine lnd2atm(bounds, & waterflux_inst%qflx_parflow_col (bounds%begc:bounds%endc, :), & lnd2atm_inst%qflx_parflow_grc (bounds%begg:bounds%endg, :), & c2l_scale_type= 'unity', l2g_scale_type='unity' ) + + ! Porosity does not scale with size, so it cannot use c2g. Average over the hydrologically + ! active columns only. Not coupled cells get out of range value, must be <0. + pfl_eff_porosity_wt(bounds%begg:bounds%endg) = 0._r8 + do c = bounds%begc, bounds%endc + if (col%hydrologically_active(c) .and. col%wtgcell(c) > 0._r8) then + g = col%gridcell(c) + pfl_eff_porosity_wt(g) = pfl_eff_porosity_wt(g) + col%wtgcell(c) + end if + end do + do g = bounds%begg, bounds%endg + if (pfl_eff_porosity_wt(g) > 0._r8) then + lnd2atm_inst%pfl_eff_porosity_grc(g,:) = 0._r8 + else + lnd2atm_inst%pfl_eff_porosity_grc(g,:) = pfl_uncoupled + end if + end do + do c = bounds%begc, bounds%endc + if (col%hydrologically_active(c) .and. col%wtgcell(c) > 0._r8) then + g = col%gridcell(c) + do j = 1, nlevgrnd + lnd2atm_inst%pfl_eff_porosity_grc(g,j) = lnd2atm_inst%pfl_eff_porosity_grc(g,j) & + + col%wtgcell(c) * soilhydrology_inst%pfl_eff_porosity_col(c,j) + end do + end if + end do + do g = bounds%begg, bounds%endg + if (pfl_eff_porosity_wt(g) > 0._r8) then + lnd2atm_inst%pfl_eff_porosity_grc(g,:) = lnd2atm_inst%pfl_eff_porosity_grc(g,:) & + / pfl_eff_porosity_wt(g) + end if + end do #endif #ifdef USE_PDAF diff --git a/src/clm5/main/lnd2atmType.F90 b/src/clm5/main/lnd2atmType.F90 index 11e93930a1..915c7b2065 100644 --- a/src/clm5/main/lnd2atmType.F90 +++ b/src/clm5/main/lnd2atmType.F90 @@ -77,6 +77,7 @@ module lnd2atmType real(r8), pointer :: qflx_liq_from_ice_col (:) => null() ! liquid runoff from converted ice runoff #ifdef COUP_OAS_PFL real(r8), pointer :: qflx_parflow_grc (:,:) => null() ! source/sink flux per soil layer sent to ParFlow [1/hr] [- out from root] + real(r8), pointer :: pfl_eff_porosity_grc (:,:) => null() ! effective porosity per soil layer sent to ParFlow [m3/m3] #endif real(r8), pointer :: qirrig_grc (:) => null() ! irrigation flux @@ -191,6 +192,7 @@ subroutine InitAllocate(this, bounds) allocate(this%qflx_liq_from_ice_col(begc:endc)) ; this%qflx_liq_from_ice_col(:) =ival #ifdef COUP_OAS_PFL allocate(this%qflx_parflow_grc (begg:endg,1:nlevsoi)); this%qflx_parflow_grc (:,:) =ival + allocate(this%pfl_eff_porosity_grc (begg:endg,1:nlevgrnd)); this%pfl_eff_porosity_grc(:,:) =-9999._r8 #endif allocate(this%qirrig_grc (begg:endg)) ; this%qirrig_grc (:) =ival @@ -344,6 +346,10 @@ subroutine InitHistory(this, bounds) call hist_addfld2d (fname='QPARFLOW_TO_OASIS', units='1/hr', type2d='levsoi', & avgflag='A', long_name='source/sink flux per soil layer sent to ParFlow', & ptr_lnd=this%qflx_parflow_grc, default='inactive') + + call hist_addfld2d (fname='EFF_POROSITY_TO_OASIS', units='m3/m3', type2d='levgrnd', & + avgflag='A', long_name='effective porosity per soil layer sent to ParFlow', & + ptr_lnd=this%pfl_eff_porosity_grc, default='inactive') #endif end subroutine InitHistory diff --git a/src/clm5/oasis3/oas_defineMod.F90 b/src/clm5/oasis3/oas_defineMod.F90 index cdae0fceeb..439ecaa3ad 100644 --- a/src/clm5/oasis3/oas_defineMod.F90 +++ b/src/clm5/oasis3/oas_defineMod.F90 @@ -107,6 +107,7 @@ subroutine oas_definitions_init(bounds) var_nodims(2) = nlevgrnd ! number of fields in a bundle call oasis_def_var(oas_sat_id, "ECLM_SOILLIQ", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) call oasis_def_var(oas_psi_id, "ECLM_PSI", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) + call oasis_def_var(oas_eff_porosity_id, "ECLM_EFFPOROSITY", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) #endif #ifdef COUP_OAS_ICON diff --git a/src/clm5/oasis3/oas_sendReceiveMod.F90 b/src/clm5/oasis3/oas_sendReceiveMod.F90 index abd8db7db9..008135fba6 100644 --- a/src/clm5/oasis3/oas_sendReceiveMod.F90 +++ b/src/clm5/oasis3/oas_sendReceiveMod.F90 @@ -48,7 +48,8 @@ subroutine oas_send(bounds, seconds_elapsed, lnd2atm_inst) integer :: info call oasis_put(oas_et_loss_id, seconds_elapsed, lnd2atm_inst%qflx_parflow_grc, info) - + call oasis_put(oas_eff_porosity_id, seconds_elapsed, lnd2atm_inst%pfl_eff_porosity_grc, info) + end subroutine oas_send #endif diff --git a/src/clm5/oasis3/oas_vardefMod.F90 b/src/clm5/oasis3/oas_vardefMod.F90 index 56ef8e0752..11639a4802 100644 --- a/src/clm5/oasis3/oas_vardefMod.F90 +++ b/src/clm5/oasis3/oas_vardefMod.F90 @@ -3,7 +3,7 @@ module oas_vardefMod save #ifdef COUP_OAS_PFL - integer :: oas_psi_id, oas_et_loss_id, oas_sat_id + integer :: oas_psi_id, oas_et_loss_id, oas_sat_id, oas_eff_porosity_id #endif #ifdef COUP_OAS_ICON From 31bab877be06a947e830c616967a611400c214fe Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Mon, 21 Sep 2026 11:45:47 +0200 Subject: [PATCH 7/9] eCLM-ParFlow: Rename oas_send/oas_receive to *_parflow - unify subroutine naming scheme --- src/clm5/cpl/lnd_comp_mct.F90 | 6 +++--- src/clm5/oasis3/oas_sendReceiveMod.F90 | 12 ++++++------ 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/clm5/cpl/lnd_comp_mct.F90 b/src/clm5/cpl/lnd_comp_mct.F90 index 56179beab6..b2a9b467c6 100644 --- a/src/clm5/cpl/lnd_comp_mct.F90 +++ b/src/clm5/cpl/lnd_comp_mct.F90 @@ -18,7 +18,7 @@ module lnd_comp_mct use oas_sendReceiveMod, only : oas_receive_icon, oas_send_icon #endif #ifdef COUP_OAS_PFL - use oas_sendReceiveMod, only : oas_receive, oas_send + use oas_sendReceiveMod, only : oas_receive_parflow, oas_send_parflow #endif #endif ! @@ -474,7 +474,7 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l) #endif #ifdef COUP_OAS_PFL - call oas_receive(bounds, time_elapsed, atm2lnd_inst) + call oas_receive_parflow(bounds, time_elapsed, atm2lnd_inst) #endif ! Run clm call t_barrierf('sync_clm_run1', mpicom) @@ -493,7 +493,7 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l) #endif #if defined(COUP_OAS_PFL) - call oas_send(bounds, time_elapsed, lnd2atm_inst) + call oas_send_parflow(bounds, time_elapsed, lnd2atm_inst) #endif ! Create l2x_l export state - add river runoff input to l2x_l if appropriate call t_startf ('lc_lnd_export') diff --git a/src/clm5/oasis3/oas_sendReceiveMod.F90 b/src/clm5/oasis3/oas_sendReceiveMod.F90 index 008135fba6..46843ee0d9 100644 --- a/src/clm5/oasis3/oas_sendReceiveMod.F90 +++ b/src/clm5/oasis3/oas_sendReceiveMod.F90 @@ -11,8 +11,8 @@ module oas_sendReceiveMod private #ifdef COUP_OAS_PFL - public :: oas_send - public :: oas_receive + public :: oas_send_parflow + public :: oas_receive_parflow #endif #ifdef COUP_OAS_ICON @@ -23,7 +23,7 @@ module oas_sendReceiveMod contains #ifdef COUP_OAS_PFL - subroutine oas_receive(bounds, seconds_elapsed, atm2lnd_inst) + subroutine oas_receive_parflow(bounds, seconds_elapsed, atm2lnd_inst) use atm2lndType, only: atm2lnd_type type(bounds_type), intent(in) :: bounds @@ -34,9 +34,9 @@ subroutine oas_receive(bounds, seconds_elapsed, atm2lnd_inst) call oasis_get(oas_psi_id, seconds_elapsed, atm2lnd_inst%pfl_psi_grc, info) call oasis_get(oas_sat_id, seconds_elapsed, atm2lnd_inst%pfl_h2osoi_liq_grc, info) - end subroutine oas_receive + end subroutine oas_receive_parflow - subroutine oas_send(bounds, seconds_elapsed, lnd2atm_inst) + subroutine oas_send_parflow(bounds, seconds_elapsed, lnd2atm_inst) use lnd2atmType, only : lnd2atm_type use spmdMod, only : mpicom use shr_mpi_mod, only: shr_mpi_barrier @@ -50,7 +50,7 @@ subroutine oas_send(bounds, seconds_elapsed, lnd2atm_inst) call oasis_put(oas_et_loss_id, seconds_elapsed, lnd2atm_inst%qflx_parflow_grc, info) call oasis_put(oas_eff_porosity_id, seconds_elapsed, lnd2atm_inst%pfl_eff_porosity_grc, info) - end subroutine oas_send + end subroutine oas_send_parflow #endif #ifdef COUP_OAS_ICON From bd4829b8f017bdc34ec222d26d8094dc95dc20af Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Tue, 22 Sep 2026 16:21:26 +0200 Subject: [PATCH 8/9] CLM-ParFlow: Remove ECLM_ICE_IMPEDANCE, ParFlow will derive it from porosity - Ice impedance is computed in ParFlow based on ECLM_EFFPOROSITY - Remove ECLM_ICE_IMPEDANCE OASIS field - Use local imped in soilwater_parflow as before --- src/clm5/biogeophys/SoilHydrologyType.F90 | 7 ------- src/clm5/biogeophys/SoilWaterMovementMod.F90 | 10 +++++----- src/clm5/main/lnd2atmMod.F90 | 6 ------ src/clm5/main/lnd2atmType.F90 | 6 ------ src/clm5/oasis3/oas_defineMod.F90 | 1 - src/clm5/oasis3/oas_sendReceiveMod.F90 | 1 - src/clm5/oasis3/oas_vardefMod.F90 | 2 +- 7 files changed, 6 insertions(+), 27 deletions(-) diff --git a/src/clm5/biogeophys/SoilHydrologyType.F90 b/src/clm5/biogeophys/SoilHydrologyType.F90 index 0ed13421d6..dbc5626b24 100644 --- a/src/clm5/biogeophys/SoilHydrologyType.F90 +++ b/src/clm5/biogeophys/SoilHydrologyType.F90 @@ -55,9 +55,6 @@ Module SoilHydrologyType real(r8), pointer :: max_infil_col (:) ! col VIC maximum infiltration rate calculated in VIC real(r8), pointer :: i_0_col (:) ! col VIC average saturation in top soil layers real(r8), pointer :: ice_col (:,:) ! col VIC soil ice (kg/m2) for VIC soil layers -#ifdef COUP_OAS_PFL - real(r8), pointer :: ice_impedance_col (:,:) ! col ice impdeance -#endif #ifdef USE_PDAF ! Yorck @@ -138,7 +135,6 @@ subroutine InitAllocate(this, bounds) allocate(this%icefrac_col (begc:endc,nlevgrnd)) ; this%icefrac_col (:,:) = nan #ifdef COUP_OAS_PFL allocate(this%pfl_eff_porosity_col (begc:endc,nlevgrnd)) ; this%pfl_eff_porosity_col(:,:) = 0._r8 - allocate(this%ice_impedance_col (begc:endc,nlevgrnd)) ; this%ice_impedance_col (:,:) = nan #endif allocate(this%fcov_col (begc:endc)) ; this%fcov_col (:) = nan allocate(this%fsat_col (begc:endc)) ; this%fsat_col (:) = nan @@ -246,9 +242,6 @@ subroutine InitCold(this, bounds) do c = bounds%begc, bounds%endc this%num_substeps_col(c) = spval this%icefrac_col(c,:) = spval -#ifdef COUP_OAS_PFL - this%ice_impedance_col(c,:) = 1.0_r8 -#endif end do end subroutine InitCold diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90 index 8e2f761127..17d1f25ee5 100644 --- a/src/clm5/biogeophys/SoilWaterMovementMod.F90 +++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90 @@ -1455,7 +1455,7 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & use decompMod , only : bounds_type use clm_varctl , only : iulog, use_hydrstress use clm_varcon , only : denh2o, denice, e_ice - use clm_varpar , only : nlevsoi, nlevgrnd + use clm_varpar , only : nlevgrnd use clm_time_manager , only : get_step_size, get_nstep use SoilStateType , only : soilstate_type use SoilHydrologyType , only : soilhydrology_type @@ -1489,13 +1489,13 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & 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 , & ! Output: [real(r8) (:,:) ] fraction of ice - ice_impedance => soilhydrology_inst%ice_impedance_col, & ! Output: [real(r8) (:,:) ] ice impedance 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] @@ -1544,13 +1544,13 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, & ! 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, ice_impedance(c,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, ice_impedance(c,j)) + 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, ice_impedance(c,j), & + call soil_water_retention_curve%soil_hk(c, j, s1, imped, & soilstate_inst, hk) hk_l(c,j) = hk end do diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index e50adf4e95..82f1ea5114 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -505,12 +505,6 @@ subroutine lnd2atm(bounds, & / pfl_eff_porosity_wt(g) end if end do - - ! Ice impedance sent to ParFlow - call c2g( bounds, nlevgrnd, & - soilhydrology_inst%ice_impedance_col (bounds%begc:bounds%endc, :), & - lnd2atm_inst%ice_impedance_grc (bounds%begg:bounds%endg, :), & - c2l_scale_type= 'unity', l2g_scale_type='unity' ) #endif #ifdef USE_PDAF diff --git a/src/clm5/main/lnd2atmType.F90 b/src/clm5/main/lnd2atmType.F90 index 7253dbdb6c..915c7b2065 100644 --- a/src/clm5/main/lnd2atmType.F90 +++ b/src/clm5/main/lnd2atmType.F90 @@ -78,7 +78,6 @@ module lnd2atmType #ifdef COUP_OAS_PFL real(r8), pointer :: qflx_parflow_grc (:,:) => null() ! source/sink flux per soil layer sent to ParFlow [1/hr] [- out from root] real(r8), pointer :: pfl_eff_porosity_grc (:,:) => null() ! effective porosity per soil layer sent to ParFlow [m3/m3] - real(r8), pointer :: ice_impedance_grc (:,:) => null() ! soil ice impedance (values in [0,1]) #endif real(r8), pointer :: qirrig_grc (:) => null() ! irrigation flux @@ -194,7 +193,6 @@ subroutine InitAllocate(this, bounds) #ifdef COUP_OAS_PFL allocate(this%qflx_parflow_grc (begg:endg,1:nlevsoi)); this%qflx_parflow_grc (:,:) =ival allocate(this%pfl_eff_porosity_grc (begg:endg,1:nlevgrnd)); this%pfl_eff_porosity_grc(:,:) =-9999._r8 - allocate(this%ice_impedance_grc (begg:endg,1:nlevgrnd)); this%ice_impedance_grc (:,:) =1.0_r8 #endif allocate(this%qirrig_grc (begg:endg)) ; this%qirrig_grc (:) =ival @@ -352,10 +350,6 @@ subroutine InitHistory(this, bounds) call hist_addfld2d (fname='EFF_POROSITY_TO_OASIS', units='m3/m3', type2d='levgrnd', & avgflag='A', long_name='effective porosity per soil layer sent to ParFlow', & ptr_lnd=this%pfl_eff_porosity_grc, default='inactive') - this%ice_impedance_grc(begg:endg, :) = 1._r8 - call hist_addfld2d (fname='ICE_IMPEDANCE', units='unitless', type2d='levgrnd', & - avgflag='A', long_name='soil ice impedance sent to ParFlow', & - ptr_lnd=this%ice_impedance_grc, default='inactive') #endif end subroutine InitHistory diff --git a/src/clm5/oasis3/oas_defineMod.F90 b/src/clm5/oasis3/oas_defineMod.F90 index 629343203b..0ecbac760a 100644 --- a/src/clm5/oasis3/oas_defineMod.F90 +++ b/src/clm5/oasis3/oas_defineMod.F90 @@ -108,7 +108,6 @@ subroutine oas_definitions_init(bounds) call oasis_def_var(oas_sat_id, "ECLM_SOILLIQ", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) call oasis_def_var(oas_psi_id, "ECLM_PSI", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) call oasis_def_var(oas_eff_porosity_id, "ECLM_EFFPOROSITY", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) - call oasis_def_var(oas_ice_impedance_id, "ECLM_ICE_IMPEDANCE", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) #endif #ifdef COUP_OAS_ICON diff --git a/src/clm5/oasis3/oas_sendReceiveMod.F90 b/src/clm5/oasis3/oas_sendReceiveMod.F90 index 2b5ca615bf..46843ee0d9 100644 --- a/src/clm5/oasis3/oas_sendReceiveMod.F90 +++ b/src/clm5/oasis3/oas_sendReceiveMod.F90 @@ -49,7 +49,6 @@ subroutine oas_send_parflow(bounds, seconds_elapsed, lnd2atm_inst) call oasis_put(oas_et_loss_id, seconds_elapsed, lnd2atm_inst%qflx_parflow_grc, info) call oasis_put(oas_eff_porosity_id, seconds_elapsed, lnd2atm_inst%pfl_eff_porosity_grc, info) - call oasis_put(oas_ice_impedance_id, seconds_elapsed, lnd2atm_inst%ice_impedance_grc, info) end subroutine oas_send_parflow #endif diff --git a/src/clm5/oasis3/oas_vardefMod.F90 b/src/clm5/oasis3/oas_vardefMod.F90 index 69a5c0bb63..11639a4802 100644 --- a/src/clm5/oasis3/oas_vardefMod.F90 +++ b/src/clm5/oasis3/oas_vardefMod.F90 @@ -3,7 +3,7 @@ module oas_vardefMod save #ifdef COUP_OAS_PFL - integer :: oas_psi_id, oas_et_loss_id, oas_sat_id, oas_eff_porosity_id, oas_ice_impedance_id + integer :: oas_psi_id, oas_et_loss_id, oas_sat_id, oas_eff_porosity_id #endif #ifdef COUP_OAS_ICON From 024a9b52f9aafcae53138672d6c52c6ef3cf9c82 Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Tue, 22 Sep 2026 16:24:34 +0200 Subject: [PATCH 9/9] eCLM-ParFlow: add possibility to output WATSAT --- src/clm5/biogeophys/SoilStateType.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/clm5/biogeophys/SoilStateType.F90 b/src/clm5/biogeophys/SoilStateType.F90 index 18baac4976..de718d9ed1 100644 --- a/src/clm5/biogeophys/SoilStateType.F90 +++ b/src/clm5/biogeophys/SoilStateType.F90 @@ -293,8 +293,12 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='urban factor limiting ground evap', & ptr_col=this%soilalpha_u_col, set_nourb=spval, default='inactive') +#ifdef COUP_OAS_PFL + if (.true.) then +#else if (use_cn) then - this%watsat_col(begc:endc,:) = spval +#endif + this%watsat_col(begc:endc,:) = spval call hist_addfld2d (fname='watsat', units='m^3/m^3', type2d='levgrnd', & avgflag='A', long_name='water saturated', & ptr_col=this%watsat_col, default='inactive')