From 46ea27ed43a1c1096c8293b9800ae75d3a28edda Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 20 Nov 2025 16:40:17 +0100 Subject: [PATCH 1/5] Soil Hydraulic Parameter Perturbations for eCLM --- .../biogeophys/SoilStateInitTimeConstMod.F90 | 44 +++++++++---------- 1 file changed, 20 insertions(+), 24 deletions(-) diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index b3054d47dd..7652b216cd 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -150,7 +150,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) real(r8) ,pointer :: gti (:) ! read in - fmax real(r8) ,pointer :: sand3d (:,:) ! read in - soil texture: percent sand (needs to be a pointer for use in ncdio) real(r8) ,pointer :: clay3d (:,:) ! read in - soil texture: percent clay (needs to be a pointer for use in ncdio) -#ifdef USE_PDAF +! SHP start real(r8) ,pointer :: psis_sat (:,:) ! read in - soil parameter: sucsat (needs to be a pointer for use in ncdio) real(r8) ,pointer :: shape_param (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) real(r8) ,pointer :: thetas (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) @@ -160,16 +160,16 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) real(r8) ,pointer :: shape_param_adj (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) real(r8) ,pointer :: thetas_adj (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) real(r8) ,pointer :: ks_adj (:,:) ! read in - soil parameter: xksat (needs to be a pointer for use in ncdio) -#endif +! SHP end real(r8) ,pointer :: organic3d (:,:) ! read in - organic matter: kg/m3 (needs to be a pointer for use in ncdio) character(len=256) :: locfn ! local filename integer :: ipedof integer :: begp, endp integer :: begc, endc integer :: begg, endg -#ifdef USE_PDAF +! SHP start logical :: parameters_in_file, parameters_in_file_adj -#endif +! SHP end !----------------------------------------------------------------------- begp = bounds%begp; endp= bounds%endp @@ -239,7 +239,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(sand3d(begg:endg,nlevsoifl)) allocate(clay3d(begg:endg,nlevsoifl)) -#ifdef USE_PDAF +! SHP start allocate(thetas(begg:endg,nlevsoifl)) allocate(shape_param(begg:endg,nlevsoifl)) allocate(psis_sat(begg:endg,nlevsoifl)) @@ -249,7 +249,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(shape_param_adj(begg:endg,nlevgrnd)) allocate(psis_sat_adj(begg:endg,nlevgrnd)) allocate(ks_adj(begg:endg,nlevgrnd)) -#endif +! SHP end ! Determine organic_max from parameter file @@ -287,7 +287,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call endrun(msg=' ERROR: PCT_CLAY NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if -#ifdef USE_PDAF +! SHP start ! include option to also read hydraulic parameters from file. Keep it variable so that the code also works for surface files that were ! generated without parameter perturbation and parameter as input variables @@ -340,7 +340,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) parameters_in_file_adj = .False. end if end if -#endif +! SHP end do p = begp,endp g = patch%gridcell(p) if ( sand3d(g,1)+clay3d(g,1) == 0.0_r8 )then @@ -533,7 +533,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call pedotransf(ipedof, sand, clay, & soilstate_inst%watsat_col(c,lev), soilstate_inst%bsw_col(c,lev), soilstate_inst%sucsat_col(c,lev), xksat) -#ifdef USE_PDAF +! SHP start ! if parameters are included in the file, watsat,... are overwritten with the values from there. If not, the pedotransfer ! function is used @@ -552,7 +552,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) xksat = ks(col%gridcell(c), nlevsoifl) end if end if -#endif +! SHP end om_watsat = max(0.93_r8 - 0.1_r8 *(zsoi(lev)/zsapric), 0.83_r8) om_b = min(2.7_r8 + 9.3_r8 *(zsoi(lev)/zsapric), 12.0_r8) om_sucsat = min(10.3_r8 - 0.2_r8 *(zsoi(lev)/zsapric), 10.1_r8) @@ -561,15 +561,13 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%bd_col(c,lev) = (1._r8 - soilstate_inst%watsat_col(c,lev))*2.7e3_r8 soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac) * soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac tkm = (1._r8-om_frac) * (8.80_r8*sand+2.92_r8*clay)/(sand+clay)+om_tkm*om_frac ! W/(m K) -#ifndef USE_PDAF - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b -#else +! SHP adapt start if (parameters_in_file) then soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b else - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b end if -#endif +! SHP adapt end soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac) * soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac soilstate_inst%hksat_min_col(c,lev) = xksat @@ -593,7 +591,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) end if soilstate_inst%hksat_col(c,lev) = uncon_frac*uncon_hksat + (perc_frac*om_frac)*om_hksat -#ifdef USE_PDAF +! SHP start if (parameters_in_file_adj) then ! Use values from the file for the soil layers soilstate_inst%watsat_col(c,lev) = thetas_adj(col%gridcell(c), lev) @@ -601,7 +599,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%sucsat_col(c,lev) = psis_sat_adj(col%gridcell(c), lev) soilstate_inst%hksat_col(c,lev) = ks_adj(col%gridcell(c), lev) ! mm/s end if -#endif +! SHP end soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev)) soilstate_inst%tksatu_col(c,lev) = soilstate_inst%tkmg_col(c,lev)*0.57_r8**soilstate_inst%watsat_col(c,lev) @@ -675,15 +673,13 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac)*soilstate_inst%watsat_col(c,lev) + om_watsat_lake * om_frac tkm = (1._r8-om_frac)*(8.80_r8*sand+2.92_r8*clay)/(sand+clay) + om_tkm * om_frac ! W/(m K) -#ifndef USE_PDAF - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake -#else +! SHP adapt start if (parameters_in_file .or. parameters_in_file_adj) then soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake else - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake end if -#endif +! SHP adapt end soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac)*soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac @@ -747,10 +743,10 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) deallocate(sand3d, clay3d, organic3d) deallocate(zisoifl, zsoifl, dzsoifl) -#ifdef USE_PDAF +! SHP start deallocate(thetas, shape_param, psis_sat, ks) deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj) -#endif +! SHP end end subroutine SoilStateInitTimeConst From a5bfe2aadee88346d1e93ae662a7263f64172c22 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 21 Nov 2025 16:24:00 +0100 Subject: [PATCH 2/5] namelist inputs: `parameters_in_file` and `parameters_in_file_adj` This way, additional memory is only allocated, when parameters are read from file. --- docs/users_guide/introduction/perturbation.md | 84 ++++++++++++++--- .../biogeophys/SoilStateInitTimeConstMod.F90 | 90 +++++++++++-------- 2 files changed, 122 insertions(+), 52 deletions(-) diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index 9fd8ec73b5..14fd44308f 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -159,20 +159,78 @@ with the `_adj` suffix: - If adjusted parameters are not present in the surface file, eCLM falls back to the original parameters or pedotransfer functions -#### Parameter Priority and Fallback Option - -The code in `SoilStateInitTimeConstMod.F90` follows this priority order: - -1. **First priority**: Read `_adj` parameters (if present) → applies - to all `nlevgrnd` layers → **overwrites** organic matter mixing - results -2. **Second priority**: Read original parameters → applies to first 10 - layers (`nlevsoifl`) → undergoes organic matter mixing -3. **Fallback**: Use pedotransfer functions if parameters are absent - from the file +#### Namelist Configuration + +The soil hydraulic parameter reading behavior is controlled by two +namelist settings in the `clm_soilstate_inparm` section of the `lnd_in` +namelist file: + +##### `parameters_in_file` + +**Type:** logical +**Default:** `.false.` +**Description:** Controls whether to read baseline hydraulic parameters +from the surface dataset file. + +When set to `.true.`: +- eCLM reads `THETAS`, `SHAPE_PARAM`, `PSIS_SAT`, and `KSAT` from the + surface file +- Parameters apply to the first 10 soil layers (`nlevsoifl=10`) +- Parameters undergo organic matter mixing +- If any required variable is missing, the model aborts with an error + message + +When set to `.false.` (default): +- Hydraulic parameters are computed via pedotransfer functions from + sand and clay fractions +- No parameters are read from the surface file + +##### `parameters_in_file_adj` + +**Type:** logical +**Default:** `.false.` +**Description:** Controls whether to read adjusted hydraulic parameters +from the surface dataset file. + +When set to `.true.`: +- eCLM reads `THETAS_adj`, `SHAPE_PARAM_adj`, `PSIS_SAT_adj`, and + `KSAT_adj` from the surface file +- Parameters apply to **all** `nlevgrnd` soil layers (typically 25 + layers) +- Adjusted parameters **overwrite** the results from organic matter + mixing +- If any required variable is missing, the model aborts with an error + message + +When set to `.false.` (default): +- No adjusted parameters are read from the surface file +- Organic matter mixing results are used as final parameter values + +##### Example Configuration + +```fortran +&clm_soilstate_inparm + organic_frac_squared = .false. + parameters_in_file = .false. + parameters_in_file_adj = .true. +/ +``` -This hierarchical approach ensures maximum flexibility for both -standard simulations and data assimilation applications. +This configuration: +- Reads adjusted parameters from the surface file which override the + organic matter mixing results +- Is typical for ensemble data assimilation applications with perturbed + soil parameters + +##### Error Handling + +Both namelist parameters enforce strict error checking: +- If set to `.true.`, **all** required parameters must be present in + the surface file +- Missing variables trigger an immediate model abort with a descriptive + error message +- This ensures users are explicitly aware when parameter files are + incomplete #### Note about the Brooks-Corey Shape Parameter #### diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index 7652b216cd..9a42e74185 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -22,6 +22,8 @@ module SoilStateInitTimeConstMod ! !PRIVATE DATA: ! Control variables (from namelist) logical, private :: organic_frac_squared ! If organic fraction should be squared (as in CLM4.5) + logical, private :: parameters_in_file ! If soil hydraulic parameters should be read from file + logical, private :: parameters_in_file_adj ! If adjusted soil hydraulic parameters should be read from file character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -55,11 +57,14 @@ subroutine ReadNL( nlfilename ) character(len=*), parameter :: nl_name = 'clm_soilstate_inparm' ! Namelist name ! MUST agree with name in namelist and read - namelist / clm_soilstate_inparm / organic_frac_squared + namelist / clm_soilstate_inparm / organic_frac_squared, parameters_in_file, & + parameters_in_file_adj ! preset values organic_frac_squared = .false. + parameters_in_file = .false. + parameters_in_file_adj = .false. if ( masterproc )then @@ -80,6 +85,14 @@ subroutine ReadNL( nlfilename ) end if call shr_mpi_bcast(organic_frac_squared, mpicom) + call shr_mpi_bcast(parameters_in_file, mpicom) + call shr_mpi_bcast(parameters_in_file_adj, mpicom) + + ! Check for incompatible namelist settings + if (parameters_in_file .and. parameters_in_file_adj) then + call endrun(msg=' ERROR: parameters_in_file and parameters_in_file_adj cannot both be .true.'//& + errmsg(sourcefile, __LINE__)) + end if end subroutine ReadNL @@ -167,9 +180,6 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) integer :: begp, endp integer :: begc, endc integer :: begg, endg -! SHP start - logical :: parameters_in_file, parameters_in_file_adj -! SHP end !----------------------------------------------------------------------- begp = bounds%begp; endp= bounds%endp @@ -240,15 +250,19 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(sand3d(begg:endg,nlevsoifl)) allocate(clay3d(begg:endg,nlevsoifl)) ! SHP start - allocate(thetas(begg:endg,nlevsoifl)) - allocate(shape_param(begg:endg,nlevsoifl)) - allocate(psis_sat(begg:endg,nlevsoifl)) - allocate(ks(begg:endg,nlevsoifl)) - - allocate(thetas_adj(begg:endg,nlevgrnd)) - allocate(shape_param_adj(begg:endg,nlevgrnd)) - allocate(psis_sat_adj(begg:endg,nlevgrnd)) - allocate(ks_adj(begg:endg,nlevgrnd)) + if(parameters_in_file) then + allocate(thetas(begg:endg,nlevsoifl)) + allocate(shape_param(begg:endg,nlevsoifl)) + allocate(psis_sat(begg:endg,nlevsoifl)) + allocate(ks(begg:endg,nlevsoifl)) + end if + + if(parameters_in_file_adj) then + allocate(thetas_adj(begg:endg,nlevgrnd)) + allocate(shape_param_adj(begg:endg,nlevgrnd)) + allocate(psis_sat_adj(begg:endg,nlevgrnd)) + allocate(ks_adj(begg:endg,nlevgrnd)) + end if ! SHP end ! Determine organic_max from parameter file @@ -279,65 +293,59 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call ncd_io(ncid=ncid, varname='PCT_SAND', flag='read', data=sand3d, dim1name=grlnd, readvar=readvar) if (.not. readvar) then - call endrun(msg=' ERROR: PCT_SAND NOT on surfdata file'//errMsg(sourcefile, __LINE__)) + call endrun(msg=' ERROR: PCT_SAND NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if call ncd_io(ncid=ncid, varname='PCT_CLAY', flag='read', data=clay3d, dim1name=grlnd, readvar=readvar) if (.not. readvar) then - call endrun(msg=' ERROR: PCT_CLAY NOT on surfdata file'//errMsg(sourcefile, __LINE__)) + call endrun(msg=' ERROR: PCT_CLAY NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if ! SHP start ! include option to also read hydraulic parameters from file. Keep it variable so that the code also works for surface files that were ! generated without parameter perturbation and parameter as input variables - call ncd_io(ncid=ncid, varname='THETAS', flag='read', data=thetas, dim1name=grlnd, readvar=readvar) - if (.not. readvar) then - parameters_in_file = .False. - else - parameters_in_file = .True. - end if + if (parameters_in_file) then + call ncd_io(ncid=ncid, varname='THETAS', flag='read', data=thetas, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + call endrun(msg=' ERROR: THETAS NOT on surfdata file'//errMsg(sourcefile, __LINE__)) + end if - if (parameters_in_file) then ! read also other paramters, if one of them is not included, use sand and clay to compute parameters - ! via pedotransfer function call ncd_io(ncid=ncid, varname='SHAPE_PARAM', flag='read', data=shape_param, dim1name=grlnd, readvar=readvar) if (.not. readvar) then - parameters_in_file = .False. + call endrun(msg=' ERROR: SHAPE_PARAM NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if call ncd_io(ncid=ncid, varname='PSIS_SAT', flag='read', data=psis_sat, dim1name=grlnd, readvar=readvar) if (.not. readvar) then - parameters_in_file = .False. + call endrun(msg=' ERROR: PSIS_SAT NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if call ncd_io(ncid=ncid, varname='KSAT', flag='read', data=ks, dim1name=grlnd, readvar=readvar) if (.not. readvar) then - parameters_in_file = .False. + call endrun(msg=' ERROR: KSAT NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if end if - call ncd_io(ncid=ncid, varname='THETAS_adj', flag='read', data=thetas_adj, dim1name=grlnd, readvar=readvar) - if (.not. readvar) then - parameters_in_file_adj = .False. - else - parameters_in_file_adj = .True. - end if + if (parameters_in_file_adj) then + call ncd_io(ncid=ncid, varname='THETAS_adj', flag='read', data=thetas_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + call endrun(msg=' ERROR: THETAS_ADJ NOT on surfdata file'//errMsg(sourcefile, __LINE__)) + end if - if (parameters_in_file_adj) then ! read also other paramters, if one of them is not included, use sand and clay to compute parameters - ! via pedotransfer function call ncd_io(ncid=ncid, varname='SHAPE_PARAM_adj', flag='read', data=shape_param_adj, dim1name=grlnd, readvar=readvar) if (.not. readvar) then - parameters_in_file_adj = .False. + call endrun(msg=' ERROR: SHAPE_PARAM_adj NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if call ncd_io(ncid=ncid, varname='PSIS_SAT_adj', flag='read', data=psis_sat_adj, dim1name=grlnd, readvar=readvar) if (.not. readvar) then - parameters_in_file_adj = .False. + call endrun(msg=' ERROR: PSIS_SAT_adj NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if call ncd_io(ncid=ncid, varname='KSAT_adj', flag='read', data=ks_adj, dim1name=grlnd, readvar=readvar) if (.not. readvar) then - parameters_in_file_adj = .False. + call endrun(msg=' ERROR: KSAT_adj NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if end if ! SHP end @@ -744,8 +752,12 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) deallocate(sand3d, clay3d, organic3d) deallocate(zisoifl, zsoifl, dzsoifl) ! SHP start - deallocate(thetas, shape_param, psis_sat, ks) - deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj) + if(parameters_in_file) then + deallocate(thetas, shape_param, psis_sat, ks) + end if + if(parameters_in_file_adj) then + deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj) + end if ! SHP end end subroutine SoilStateInitTimeConst From de5e97251b3d7ee0c690391da27de902785357b7 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 25 Nov 2025 10:11:29 +0100 Subject: [PATCH 3/5] more descriptive namelist entry `soil_hyd_inparm_from_surfdata` --- docs/users_guide/introduction/perturbation.md | 14 +++--- .../biogeophys/SoilStateInitTimeConstMod.F90 | 43 ++++++++++--------- 2 files changed, 29 insertions(+), 28 deletions(-) diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index 931688d941..f3b16a27af 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -168,7 +168,7 @@ The soil hydraulic parameter reading behavior is controlled by two namelist settings in the `clm_soilstate_inparm` section of the `lnd_in` namelist file: -##### `parameters_in_file` +##### `soil_hyd_inparm_from_surfdata` **Type:** logical **Default:** `.false.` @@ -188,12 +188,12 @@ When set to `.false.` (default): sand and clay fractions - No parameters are read from the surface file -##### `parameters_in_file_adj` +##### `soil_hyd_inparm_from_surfdata_adj` **Type:** logical **Default:** `.false.` -**Description:** Controls whether to read adjusted hydraulic parameters -from the surface dataset file. +**Description:** Controls whether to read organic-matter-adjusted +hydraulic parameters from the surface dataset file. When set to `.true.`: - eCLM reads `THETAS_adj`, `SHAPE_PARAM_adj`, `PSIS_SAT_adj`, and @@ -206,7 +206,7 @@ When set to `.true.`: message When set to `.false.` (default): -- No adjusted parameters are read from the surface file +- No organic-matter-adjusted parameters are read from the surface file - Organic matter mixing results are used as final parameter values ##### Example Configuration @@ -214,8 +214,8 @@ When set to `.false.` (default): ```fortran &clm_soilstate_inparm organic_frac_squared = .false. - parameters_in_file = .false. - parameters_in_file_adj = .true. + soil_hyd_inparm_from_surfdata = .false. + soil_hyd_inparm_from_surfdata_adj = .true. / ``` diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index 9a42e74185..ab75501ad9 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -22,8 +22,8 @@ module SoilStateInitTimeConstMod ! !PRIVATE DATA: ! Control variables (from namelist) logical, private :: organic_frac_squared ! If organic fraction should be squared (as in CLM4.5) - logical, private :: parameters_in_file ! If soil hydraulic parameters should be read from file - logical, private :: parameters_in_file_adj ! If adjusted soil hydraulic parameters should be read from file + logical, private :: soil_hyd_inparm_from_surfdata ! If soil hydraulic parameters should be read from file + logical, private :: soil_hyd_inparm_from_surfdata_adj ! If adjusted soil hydraulic parameters should be read from file character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -57,14 +57,14 @@ subroutine ReadNL( nlfilename ) character(len=*), parameter :: nl_name = 'clm_soilstate_inparm' ! Namelist name ! MUST agree with name in namelist and read - namelist / clm_soilstate_inparm / organic_frac_squared, parameters_in_file, & - parameters_in_file_adj + namelist / clm_soilstate_inparm / organic_frac_squared, soil_hyd_inparm_from_surfdata, & + soil_hyd_inparm_from_surfdata_adj ! preset values organic_frac_squared = .false. - parameters_in_file = .false. - parameters_in_file_adj = .false. + soil_hyd_inparm_from_surfdata = .false. + soil_hyd_inparm_from_surfdata_adj = .false. if ( masterproc )then @@ -85,12 +85,12 @@ subroutine ReadNL( nlfilename ) end if call shr_mpi_bcast(organic_frac_squared, mpicom) - call shr_mpi_bcast(parameters_in_file, mpicom) - call shr_mpi_bcast(parameters_in_file_adj, mpicom) + call shr_mpi_bcast(soil_hyd_inparm_from_surfdata, mpicom) + call shr_mpi_bcast(soil_hyd_inparm_from_surfdata_adj, mpicom) ! Check for incompatible namelist settings - if (parameters_in_file .and. parameters_in_file_adj) then - call endrun(msg=' ERROR: parameters_in_file and parameters_in_file_adj cannot both be .true.'//& + if (soil_hyd_inparm_from_surfdata .and. soil_hyd_inparm_from_surfdata_adj) then + call endrun(msg=' ERROR: soil_hyd_inparm_from_surfdata and soil_hyd_inparm_from_surfdata_adj cannot both be .true.'//& errmsg(sourcefile, __LINE__)) end if @@ -250,14 +250,14 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(sand3d(begg:endg,nlevsoifl)) allocate(clay3d(begg:endg,nlevsoifl)) ! SHP start - if(parameters_in_file) then + if(soil_hyd_inparm_from_surfdata) then allocate(thetas(begg:endg,nlevsoifl)) allocate(shape_param(begg:endg,nlevsoifl)) allocate(psis_sat(begg:endg,nlevsoifl)) allocate(ks(begg:endg,nlevsoifl)) end if - if(parameters_in_file_adj) then + if(soil_hyd_inparm_from_surfdata_adj) then allocate(thetas_adj(begg:endg,nlevgrnd)) allocate(shape_param_adj(begg:endg,nlevgrnd)) allocate(psis_sat_adj(begg:endg,nlevgrnd)) @@ -305,7 +305,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) ! include option to also read hydraulic parameters from file. Keep it variable so that the code also works for surface files that were ! generated without parameter perturbation and parameter as input variables - if (parameters_in_file) then + if (soil_hyd_inparm_from_surfdata) then call ncd_io(ncid=ncid, varname='THETAS', flag='read', data=thetas, dim1name=grlnd, readvar=readvar) if (.not. readvar) then call endrun(msg=' ERROR: THETAS NOT on surfdata file'//errMsg(sourcefile, __LINE__)) @@ -327,7 +327,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) end if end if - if (parameters_in_file_adj) then + if (soil_hyd_inparm_from_surfdata_adj) then call ncd_io(ncid=ncid, varname='THETAS_adj', flag='read', data=thetas_adj, dim1name=grlnd, readvar=readvar) if (.not. readvar) then call endrun(msg=' ERROR: THETAS_ADJ NOT on surfdata file'//errMsg(sourcefile, __LINE__)) @@ -545,7 +545,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) ! if parameters are included in the file, watsat,... are overwritten with the values from there. If not, the pedotransfer ! function is used - if (parameters_in_file) then + if (soil_hyd_inparm_from_surfdata) then if (lev <= nlevsoifl) then ! Use values from the file for the soil layers soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), lev) @@ -570,7 +570,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac) * soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac tkm = (1._r8-om_frac) * (8.80_r8*sand+2.92_r8*clay)/(sand+clay)+om_tkm*om_frac ! W/(m K) ! SHP adapt start - if (parameters_in_file) then + if (soil_hyd_inparm_from_surfdata) then soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b else soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b @@ -600,8 +600,9 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%hksat_col(c,lev) = uncon_frac*uncon_hksat + (perc_frac*om_frac)*om_hksat ! SHP start - if (parameters_in_file_adj) then - ! Use values from the file for the soil layers + if (soil_hyd_inparm_from_surfdata_adj) then + ! Ovewrite organic-matter-adjusted parameters from + ! the file for all ground layers soilstate_inst%watsat_col(c,lev) = thetas_adj(col%gridcell(c), lev) soilstate_inst%bsw_col(c,lev) = shape_param_adj(col%gridcell(c), lev) soilstate_inst%sucsat_col(c,lev) = psis_sat_adj(col%gridcell(c), lev) @@ -682,7 +683,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) tkm = (1._r8-om_frac)*(8.80_r8*sand+2.92_r8*clay)/(sand+clay) + om_tkm * om_frac ! W/(m K) ! SHP adapt start - if (parameters_in_file .or. parameters_in_file_adj) then + if (soil_hyd_inparm_from_surfdata .or. soil_hyd_inparm_from_surfdata_adj) then soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake else soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake @@ -752,10 +753,10 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) deallocate(sand3d, clay3d, organic3d) deallocate(zisoifl, zsoifl, dzsoifl) ! SHP start - if(parameters_in_file) then + if(soil_hyd_inparm_from_surfdata) then deallocate(thetas, shape_param, psis_sat, ks) end if - if(parameters_in_file_adj) then + if(soil_hyd_inparm_from_surfdata_adj) then deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj) end if ! SHP end From e12540b2f31c1042831aaafc7aca0fbfe7a655d4 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 25 Nov 2025 13:06:50 +0100 Subject: [PATCH 4/5] docs update - separate atmospheric forcing noise and soil hydraulic parameters - change shp-sections suggesting a PDAF-only usage --- docs/_toc.yml | 4 +- docs/users_guide/introduction/perturbation.md | 235 ++---------------- .../introduction/soil_hydraulic_parameters.md | 205 +++++++++++++++ 3 files changed, 222 insertions(+), 222 deletions(-) create mode 100644 docs/users_guide/introduction/soil_hydraulic_parameters.md diff --git a/docs/_toc.yml b/docs/_toc.yml index 659415e6bb..7dd640808d 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -10,7 +10,9 @@ parts: - file: users_guide/introduction/introduction title: Scientific Background - file: users_guide/introduction/perturbation - title: Perturbation routines + title: Atmospheric forcing noise + - file: users_guide/introduction/soil_hydraulic_parameters + title: Soil hydraulic properties from surface file - caption: User's Guide chapters: diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index f3b16a27af..206281675e 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -1,15 +1,11 @@ -# Perturbation routines # +# Atmospheric forcing noise # -This section describes **perturbation capabilities for eCLM-PDAF**. - -The implementation of the Soil Hydraulic Parameter perturbation was -first introduced by Fernand Eloundou and adapted by Yorck -Ewerdwalbesloh. +This section describes **atmospheric forcing noise for eCLM-PDAF**. The implementation of the atmospheric forcing perturbation routines was first introduced by Yorck Ewerdwalbesloh. -## Description +## General Description eCLM implements perturbations for atmospheric forcing variables (temperature, precipitation, longwave radiation, and shortwave @@ -46,202 +42,7 @@ is adjustable based on application requirements. Perturbation fields are read from preprocessed noise files and applied to atmospheric forcing variables within the eCLM model during runtime. -## Key Components ## - -### 1. Soil Parameter Perturbation (`SoilStateInitTimeConstMod.F90`) ### - -Allows reading perturbed hydraulic parameters from input files instead -of computing them via pedotransfer functions. - -**Parameters:** -- `THETAS` - saturated water content -- `SHAPE_PARAM` - Brooks-Corey shape parameter (`bsw`) -- `PSIS_SAT` - saturated matric potential (`sucsat`) -- `KSAT` - saturated hydraulic conductivity (`xksat`) - -**Implementation:** -- Hydraulic parameters read from surface file -- Reads EITHER original parameters (for `nlevsoifl=10` soil layers) - and applies organic matter mixing OR: Reads adjusted parameters with - `_adj` suffix (for all `nlevgrnd` layers) and **overwrites** - parameters from organic matter mixing. -- Falls back to pedotransfer functions if parameters aren't in the file -- Modifies organic matter mixing to preserve perturbed parameter - values - -#### Surface File - -Soil hydraulic parameters are read from the **surface dataset file**, -which is specified in the `lnd_in` namelist file under key `fsurdat`. - -##### Sufrace File: Parameter Structure and Dimensions - -All soil hydraulic parameters are stored as three-dimensional arrays -with the structure: - -``` - (nlevgrnd, lsmlat, lsmlon) -``` - -Where: -- `nlevgrnd`: Number of soil layers (vertical dimension, typically 25 - layers in eCLM simulations, see namelist file key - `soil_layerstruct`) -- `lsmlat`: Number of latitude grid points -- `lsmlon`: Number of longitude grid points - -This 3D structure enables depth-varying, spatially-distributed soil -properties across the model domain. - -##### Surface File: Input Parameter Format - -Two options -- original parameters that undergo organic mixing -- adjusted parameters that are used by input values - -###### Original Parameters - -The surface file must contain the following baseline parameters: - -``` text - THETAS(nlevgrnd, lsmlat, lsmlon) - long_name: "Porosity" - units: "vol/vol" - - SHAPE_PARAM(nlevgrnd, lsmlat, lsmlon) - long_name: "Shape (b) parameter" - units: "unitless" - - PSIS_SAT(nlevgrnd, lsmlat, lsmlon) - long_name: "Saturated soil matric potential" - units: "mmH2O" - - KSAT(nlevgrnd, lsmlat, lsmlon) - long_name: "Saturated hydraulic conductivity" - units: "mm/s" -``` - -These parameters apply to the first 10 soil layers (`nlevsoifl=10`) and -undergo organic matter mixing. - -###### Adjusted Parameters for Data Assimilation - -When running with perturbed soil parameters (e.g., from ensemble data -assimilation), the surface file should contain adjusted parameters -with the `_adj` suffix: - -``` text - THETAS_adj(nlevgrnd, lsmlat, lsmlon) - long_name: "Porosity" - units: "vol/vol" - _FillValue: 1.e+30 - - SHAPE_PARAM_adj(nlevgrnd, lsmlat, lsmlon) - long_name: "Shape (b) parameter" - units: "unitless" - _FillValue: 1.e+30 - - PSIS_SAT_adj(nlevgrnd, lsmlat, lsmlon) - long_name: "Saturated soil matric potential" - units: "mmH2O" - _FillValue: 1000.0 - - KSAT_adj(nlevgrnd, lsmlat, lsmlon) - long_name: "Saturated hydraulic conductivity" - units: "mm/s" - _FillValue: 1.e+30 -``` - -**Important Notes:** -- Adjusted parameters apply to **all** `nlevgrnd` layers -- When present, adjusted parameters **overwrite** the results of organic - matter mixing -- `PSIS_SAT_adj` uses a different fill value (1000.0) compared to other - parameters (1.e+30), reflecting special handling for undefined matric - potential values -- If adjusted parameters are not present in the surface file, eCLM falls - back to the original parameters or pedotransfer functions - -#### Namelist Configuration - -The soil hydraulic parameter reading behavior is controlled by two -namelist settings in the `clm_soilstate_inparm` section of the `lnd_in` -namelist file: - -##### `soil_hyd_inparm_from_surfdata` - -**Type:** logical -**Default:** `.false.` -**Description:** Controls whether to read baseline hydraulic parameters -from the surface dataset file. - -When set to `.true.`: -- eCLM reads `THETAS`, `SHAPE_PARAM`, `PSIS_SAT`, and `KSAT` from the - surface file -- Parameters apply to the first 10 soil layers (`nlevsoifl=10`) -- Parameters undergo organic matter mixing -- If any required variable is missing, the model aborts with an error - message - -When set to `.false.` (default): -- Hydraulic parameters are computed via pedotransfer functions from - sand and clay fractions -- No parameters are read from the surface file - -##### `soil_hyd_inparm_from_surfdata_adj` - -**Type:** logical -**Default:** `.false.` -**Description:** Controls whether to read organic-matter-adjusted -hydraulic parameters from the surface dataset file. - -When set to `.true.`: -- eCLM reads `THETAS_adj`, `SHAPE_PARAM_adj`, `PSIS_SAT_adj`, and - `KSAT_adj` from the surface file -- Parameters apply to **all** `nlevgrnd` soil layers (typically 25 - layers) -- Adjusted parameters **overwrite** the results from organic matter - mixing -- If any required variable is missing, the model aborts with an error - message - -When set to `.false.` (default): -- No organic-matter-adjusted parameters are read from the surface file -- Organic matter mixing results are used as final parameter values - -##### Example Configuration - -```fortran -&clm_soilstate_inparm - organic_frac_squared = .false. - soil_hyd_inparm_from_surfdata = .false. - soil_hyd_inparm_from_surfdata_adj = .true. -/ -``` - -This configuration: -- Reads adjusted parameters from the surface file which override the - organic matter mixing results -- Is typical for ensemble data assimilation applications with perturbed - soil parameters - -##### Error Handling - -Both namelist parameters enforce strict error checking: -- If set to `.true.`, **all** required parameters must be present in - the surface file -- Missing variables trigger an immediate model abort with a descriptive - error message -- This ensures users are explicitly aware when parameter files are - incomplete - -#### Note about the Brooks-Corey Shape Parameter #### - -When perturbed soil parameters are read from input files, the organic -matter mixing for `bsw` uses the file-read value instead of the -hard-coded Cosby et al. (1984) Table 5 formula (`2.91 + 0.159*clay`). - -### 2. Noise-Based Forcing Perturbation ### +## 1. Noise-Based Forcing Perturbation ## When having an ensemble run, the memory consumption is too large if the forcing data is perturbed one by one so that I get a file for each @@ -274,13 +75,13 @@ data assimilation. - Detects noise fields by checking if model field name contains `"noise"` - Adjusts time coordinate reading to account for ensemble-extended noise files -#### Preparing Atmospheric Forcing Noise #### +### Preparing Atmospheric Forcing Noise ### This section describes the workflow for generating and configuring spatiotemporally correlated noise fields for atmospheric forcing perturbations. -##### Noise Generation Tool ##### +#### Noise Generation Tool #### The `correlatedNoise.py` (https://github.com/HPSCTerrSys/eCLM_atm-forcing-generator/pull/8) @@ -299,7 +100,7 @@ correlated noise fields for all ensemble members. **Output**: One NetCDF file per month containing noise fields for all ensemble members and all perturbed variables. -##### Stream File Configuration ##### +#### Stream File Configuration #### Stream files configure how eCLM reads and applies noise fields to atmospheric forcing data. Three separate stream files are required @@ -314,7 +115,7 @@ file must specify the following metadata fields: - **Perturbation variables**: `tbot_noise` (temperature), `lwdn_noise` (longwave radiation), `precn_noise` (precipitation), `swdn_noise` (shortwave radiation) -##### Stream File Examples ##### +#### Stream File Examples #### The following examples demonstrate proper stream file configuration for precipitation, solar radiation, and other atmospheric @@ -322,7 +123,7 @@ variables. Each example must be adapted with correct file paths, temporal parameters, and ensemble member IDs for your specific application. -###### Precipitation Noise Stream File ###### +##### Precipitation Noise Stream File ##### **File**: `user_datm.streams.precip_noise.stream_0000.txt` @@ -381,7 +182,7 @@ application. ``` -###### Solar Radiation Noise Stream File ###### +##### Solar Radiation Noise Stream File ##### **File**: `user_datm.streams.solar_noise.stream_0000.txt` @@ -440,7 +241,7 @@ application. ``` -###### Temperature and Longwave Radiation Noise Stream File ###### +##### Temperature and Longwave Radiation Noise Stream File ##### **File**: `user_datm.streams.other_noise.stream_0000.txt` @@ -501,7 +302,7 @@ application. ``` -### 3. DATM Integration (`datm_comp_mod.F90`) ### +## 2. DATM Integration (`datm_comp_mod.F90`) ## Passes ensemble information (`inst_index`, `dt_option`, `ninst`) from the data atmosphere (DATM) component to the stream infrastructure, @@ -510,22 +311,14 @@ CESM's multi-instance framework. ## Design Rationale ## -1. **Dual perturbation approach**: - - Parameter perturbation for soil properties (offline preprocessing) - - Forcing perturbation via temporal noise patterns (runtime) - -2. **Backward compatibility**: All changes are guarded by `USE_PDAF` +1. **Backward compatibility**: All changes are guarded by `USE_PDAF` preprocessor directives -3. **Efficiency**: Uses a single noise file for all ensemble members +2. **Efficiency**: Uses a single noise file for all ensemble members by time-shifting, avoiding storage duplication -4. **Flexibility**: Soil parameters can be perturbed or computed - traditionally based on file contents - ## Modified Source Code ## -- `src/clm5/biogeophys/SoilStateInitTimeConstMod.F90` - `src/csm_share/streams/shr_stream_mod.F90` - `src/csm_share/streams/shr_strdata_mod.F90` - `src/csm_share/streams/shr_dmodel_mod.F90` diff --git a/docs/users_guide/introduction/soil_hydraulic_parameters.md b/docs/users_guide/introduction/soil_hydraulic_parameters.md new file mode 100644 index 0000000000..d8b788fd62 --- /dev/null +++ b/docs/users_guide/introduction/soil_hydraulic_parameters.md @@ -0,0 +1,205 @@ +# Soil Hydraulic Parameters # + +This section describes **setting Soil Hydraulic Parameters from the +surface file**. + +The implementation was first introduced by Fernand Eloundou and +adapted by Yorck Ewerdwalbesloh. + +## Soil Parameter Perturbation (`SoilStateInitTimeConstMod.F90`) ## + +Allows reading perturbed hydraulic parameters from input files instead +of computing them via pedotransfer functions. + +**Parameters:** +- `THETAS` - saturated water content +- `SHAPE_PARAM` - Brooks-Corey shape parameter (`bsw`) +- `PSIS_SAT` - saturated matric potential (`sucsat`) +- `KSAT` - saturated hydraulic conductivity (`xksat`) + +**Implementation:** +- Hydraulic parameters read from surface file +- Reads EITHER original parameters (for `nlevsoifl=10` soil layers) + and applies organic matter mixing OR: Reads adjusted parameters with + `_adj` suffix (for all `nlevgrnd` layers) and **overwrites** + parameters from organic matter mixing. +- Falls back to pedotransfer functions if parameters aren't in the file +- Modifies organic matter mixing to preserve perturbed parameter + values + +### Surface File + +Soil hydraulic parameters are read from the **surface dataset file**, +which is specified in the `lnd_in` namelist file under key `fsurdat`. + +#### Sufrace File: Parameter Structure and Dimensions + +All soil hydraulic parameters are stored as three-dimensional arrays +with the structure: + +``` + (nlevgrnd, lsmlat, lsmlon) +``` + +Where: +- `nlevgrnd`: Number of soil layers (vertical dimension, typically 25 + layers in eCLM simulations, see namelist file key + `soil_layerstruct`) +- `lsmlat`: Number of latitude grid points +- `lsmlon`: Number of longitude grid points + +This 3D structure enables depth-varying, spatially-distributed soil +properties across the model domain. + +#### Surface File: Input Parameter Format + +Two options +- parameters that undergo organic mixing +- adjusted parameters that are used by input values + +##### Parameters before organic mixing + +The surface file must contain the following baseline parameters: + +``` text + THETAS(nlevgrnd, lsmlat, lsmlon) + long_name: "Porosity" + units: "vol/vol" + + SHAPE_PARAM(nlevgrnd, lsmlat, lsmlon) + long_name: "Shape (b) parameter" + units: "unitless" + + PSIS_SAT(nlevgrnd, lsmlat, lsmlon) + long_name: "Saturated soil matric potential" + units: "mmH2O" + + KSAT(nlevgrnd, lsmlat, lsmlon) + long_name: "Saturated hydraulic conductivity" + units: "mm/s" +``` + +These parameters apply to the first 10 soil layers (`nlevsoifl=10`) and +undergo organic matter mixing. + +##### Adjusted Parameters after organic mixing + +For using adjusted soil parameters from file that represent the whole +soil including organic matter, the surface file should contain +adjusted parameters with the `_adj` suffix: + +``` text + THETAS_adj(nlevgrnd, lsmlat, lsmlon) + long_name: "Porosity" + units: "vol/vol" + _FillValue: 1.e+30 + + SHAPE_PARAM_adj(nlevgrnd, lsmlat, lsmlon) + long_name: "Shape (b) parameter" + units: "unitless" + _FillValue: 1.e+30 + + PSIS_SAT_adj(nlevgrnd, lsmlat, lsmlon) + long_name: "Saturated soil matric potential" + units: "mmH2O" + _FillValue: 1000.0 + + KSAT_adj(nlevgrnd, lsmlat, lsmlon) + long_name: "Saturated hydraulic conductivity" + units: "mm/s" + _FillValue: 1.e+30 +``` + +**Important Notes:** +- Adjusted parameters apply to **all** `nlevgrnd` layers +- When present, adjusted parameters **overwrite** the results of organic + matter mixing +- `PSIS_SAT_adj` uses a different fill value (1000.0) compared to other + parameters (1.e+30), reflecting special handling for undefined matric + potential values +- If adjusted parameters are not present in the surface file, eCLM falls + back to the original parameters or pedotransfer functions + +### Namelist Configuration + +The soil hydraulic parameter reading behavior is controlled by two +namelist settings in the `clm_soilstate_inparm` section of the `lnd_in` +namelist file: + +#### `soil_hyd_inparm_from_surfdata` + +**Type:** logical +**Default:** `.false.` +**Description:** Controls whether to read baseline hydraulic parameters +from the surface dataset file. + +When set to `.true.`: +- eCLM reads `THETAS`, `SHAPE_PARAM`, `PSIS_SAT`, and `KSAT` from the + surface file +- Parameters apply to the first 10 soil layers (`nlevsoifl=10`) +- Parameters undergo organic matter mixing +- If any required variable is missing, the model aborts with an error + message + +When set to `.false.` (default): +- Hydraulic parameters are computed via pedotransfer functions from + sand and clay fractions +- No parameters are read from the surface file + +#### `soil_hyd_inparm_from_surfdata_adj` + +**Type:** logical +**Default:** `.false.` +**Description:** Controls whether to read organic-matter-adjusted +hydraulic parameters from the surface dataset file. + +When set to `.true.`: +- eCLM reads `THETAS_adj`, `SHAPE_PARAM_adj`, `PSIS_SAT_adj`, and + `KSAT_adj` from the surface file +- Parameters apply to **all** `nlevgrnd` soil layers (typically 25 + layers) +- Adjusted parameters **overwrite** the results from organic matter + mixing +- If any required variable is missing, the model aborts with an error + message + +When set to `.false.` (default): +- No organic-matter-adjusted parameters are read from the surface file +- Organic matter mixing results are used as final parameter values + +#### Example Configuration + +```fortran +&clm_soilstate_inparm + organic_frac_squared = .false. + soil_hyd_inparm_from_surfdata = .false. + soil_hyd_inparm_from_surfdata_adj = .true. +/ +``` + +This configuration: +- Reads adjusted parameters from the surface file which override the + organic matter mixing results +- Can e.g. be used for ensemble data assimilation applications with + perturbed soil parameters specified in a group of surface files, one + surface file per ensemble member. + +#### Error Handling + +Both namelist parameters enforce strict error checking: +- If set to `.true.`, **all** required parameters must be present in + the surface file +- Missing variables trigger an immediate model abort with a descriptive + error message +- This ensures users are explicitly aware when parameter files are + incomplete + +### Note about the Brooks-Corey Shape Parameter ### + +When perturbed soil parameters are read from input files, the organic +matter mixing for `bsw` uses the file-read value instead of the +hard-coded Cosby et al. (1984) Table 5 formula (`2.91 + 0.159*clay`). + +## Modified Source Code ## + +- `src/clm5/biogeophys/SoilStateInitTimeConstMod.F90` From 66774164c32cda172a967467a130a5218bf74d41 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 25 Nov 2025 13:10:02 +0100 Subject: [PATCH 5/5] change namelist parameter name to `soil_hyd_inparm_from_file[_adj]` --- .../introduction/soil_hydraulic_parameters.md | 8 ++-- .../biogeophys/SoilStateInitTimeConstMod.F90 | 40 +++++++++---------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/docs/users_guide/introduction/soil_hydraulic_parameters.md b/docs/users_guide/introduction/soil_hydraulic_parameters.md index d8b788fd62..f0f272a80c 100644 --- a/docs/users_guide/introduction/soil_hydraulic_parameters.md +++ b/docs/users_guide/introduction/soil_hydraulic_parameters.md @@ -126,7 +126,7 @@ The soil hydraulic parameter reading behavior is controlled by two namelist settings in the `clm_soilstate_inparm` section of the `lnd_in` namelist file: -#### `soil_hyd_inparm_from_surfdata` +#### `soil_hyd_inparm_from_file` **Type:** logical **Default:** `.false.` @@ -146,7 +146,7 @@ When set to `.false.` (default): sand and clay fractions - No parameters are read from the surface file -#### `soil_hyd_inparm_from_surfdata_adj` +#### `soil_hyd_inparm_from_file_adj` **Type:** logical **Default:** `.false.` @@ -172,8 +172,8 @@ When set to `.false.` (default): ```fortran &clm_soilstate_inparm organic_frac_squared = .false. - soil_hyd_inparm_from_surfdata = .false. - soil_hyd_inparm_from_surfdata_adj = .true. + soil_hyd_inparm_from_file = .false. + soil_hyd_inparm_from_file_adj = .true. / ``` diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index ab75501ad9..935e26c76a 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -22,8 +22,8 @@ module SoilStateInitTimeConstMod ! !PRIVATE DATA: ! Control variables (from namelist) logical, private :: organic_frac_squared ! If organic fraction should be squared (as in CLM4.5) - logical, private :: soil_hyd_inparm_from_surfdata ! If soil hydraulic parameters should be read from file - logical, private :: soil_hyd_inparm_from_surfdata_adj ! If adjusted soil hydraulic parameters should be read from file + logical, private :: soil_hyd_inparm_from_file ! If soil hydraulic parameters should be read from file + logical, private :: soil_hyd_inparm_from_file_adj ! If adjusted soil hydraulic parameters should be read from file character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -57,14 +57,14 @@ subroutine ReadNL( nlfilename ) character(len=*), parameter :: nl_name = 'clm_soilstate_inparm' ! Namelist name ! MUST agree with name in namelist and read - namelist / clm_soilstate_inparm / organic_frac_squared, soil_hyd_inparm_from_surfdata, & - soil_hyd_inparm_from_surfdata_adj + namelist / clm_soilstate_inparm / organic_frac_squared, soil_hyd_inparm_from_file, & + soil_hyd_inparm_from_file_adj ! preset values organic_frac_squared = .false. - soil_hyd_inparm_from_surfdata = .false. - soil_hyd_inparm_from_surfdata_adj = .false. + soil_hyd_inparm_from_file = .false. + soil_hyd_inparm_from_file_adj = .false. if ( masterproc )then @@ -85,12 +85,12 @@ subroutine ReadNL( nlfilename ) end if call shr_mpi_bcast(organic_frac_squared, mpicom) - call shr_mpi_bcast(soil_hyd_inparm_from_surfdata, mpicom) - call shr_mpi_bcast(soil_hyd_inparm_from_surfdata_adj, mpicom) + call shr_mpi_bcast(soil_hyd_inparm_from_file, mpicom) + call shr_mpi_bcast(soil_hyd_inparm_from_file_adj, mpicom) ! Check for incompatible namelist settings - if (soil_hyd_inparm_from_surfdata .and. soil_hyd_inparm_from_surfdata_adj) then - call endrun(msg=' ERROR: soil_hyd_inparm_from_surfdata and soil_hyd_inparm_from_surfdata_adj cannot both be .true.'//& + if (soil_hyd_inparm_from_file .and. soil_hyd_inparm_from_file_adj) then + call endrun(msg=' ERROR: soil_hyd_inparm_from_file and soil_hyd_inparm_from_file_adj cannot both be .true.'//& errmsg(sourcefile, __LINE__)) end if @@ -250,14 +250,14 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(sand3d(begg:endg,nlevsoifl)) allocate(clay3d(begg:endg,nlevsoifl)) ! SHP start - if(soil_hyd_inparm_from_surfdata) then + if(soil_hyd_inparm_from_file) then allocate(thetas(begg:endg,nlevsoifl)) allocate(shape_param(begg:endg,nlevsoifl)) allocate(psis_sat(begg:endg,nlevsoifl)) allocate(ks(begg:endg,nlevsoifl)) end if - if(soil_hyd_inparm_from_surfdata_adj) then + if(soil_hyd_inparm_from_file_adj) then allocate(thetas_adj(begg:endg,nlevgrnd)) allocate(shape_param_adj(begg:endg,nlevgrnd)) allocate(psis_sat_adj(begg:endg,nlevgrnd)) @@ -305,7 +305,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) ! include option to also read hydraulic parameters from file. Keep it variable so that the code also works for surface files that were ! generated without parameter perturbation and parameter as input variables - if (soil_hyd_inparm_from_surfdata) then + if (soil_hyd_inparm_from_file) then call ncd_io(ncid=ncid, varname='THETAS', flag='read', data=thetas, dim1name=grlnd, readvar=readvar) if (.not. readvar) then call endrun(msg=' ERROR: THETAS NOT on surfdata file'//errMsg(sourcefile, __LINE__)) @@ -327,7 +327,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) end if end if - if (soil_hyd_inparm_from_surfdata_adj) then + if (soil_hyd_inparm_from_file_adj) then call ncd_io(ncid=ncid, varname='THETAS_adj', flag='read', data=thetas_adj, dim1name=grlnd, readvar=readvar) if (.not. readvar) then call endrun(msg=' ERROR: THETAS_ADJ NOT on surfdata file'//errMsg(sourcefile, __LINE__)) @@ -545,7 +545,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) ! if parameters are included in the file, watsat,... are overwritten with the values from there. If not, the pedotransfer ! function is used - if (soil_hyd_inparm_from_surfdata) then + if (soil_hyd_inparm_from_file) then if (lev <= nlevsoifl) then ! Use values from the file for the soil layers soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), lev) @@ -570,7 +570,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac) * soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac tkm = (1._r8-om_frac) * (8.80_r8*sand+2.92_r8*clay)/(sand+clay)+om_tkm*om_frac ! W/(m K) ! SHP adapt start - if (soil_hyd_inparm_from_surfdata) then + if (soil_hyd_inparm_from_file) then soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b else soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b @@ -600,7 +600,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%hksat_col(c,lev) = uncon_frac*uncon_hksat + (perc_frac*om_frac)*om_hksat ! SHP start - if (soil_hyd_inparm_from_surfdata_adj) then + if (soil_hyd_inparm_from_file_adj) then ! Ovewrite organic-matter-adjusted parameters from ! the file for all ground layers soilstate_inst%watsat_col(c,lev) = thetas_adj(col%gridcell(c), lev) @@ -683,7 +683,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) tkm = (1._r8-om_frac)*(8.80_r8*sand+2.92_r8*clay)/(sand+clay) + om_tkm * om_frac ! W/(m K) ! SHP adapt start - if (soil_hyd_inparm_from_surfdata .or. soil_hyd_inparm_from_surfdata_adj) then + if (soil_hyd_inparm_from_file .or. soil_hyd_inparm_from_file_adj) then soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake else soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake @@ -753,10 +753,10 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) deallocate(sand3d, clay3d, organic3d) deallocate(zisoifl, zsoifl, dzsoifl) ! SHP start - if(soil_hyd_inparm_from_surfdata) then + if(soil_hyd_inparm_from_file) then deallocate(thetas, shape_param, psis_sat, ks) end if - if(soil_hyd_inparm_from_surfdata_adj) then + if(soil_hyd_inparm_from_file_adj) then deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj) end if ! SHP end