namelist inputs: parameters_in_file and parameters_in_file_adj

This way, additional memory is only allocated, when parameters are
read from file.

Author: jjokella

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 ####
 

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