Support multiple soil water retention curve methods

**The notes here were moved from a page on the old CESM wiki (https://wiki.ucar.edu/display/ccsm/Supporting+multiple+soil+water+retention+curve+methods+in+CLM), created by me Aug 19, 2015. I think that much of this is still relevant, though I haven't reviewed this to determine if everything here is still relevant.**

### Background

Currently (in clm4_5_1_r118, and quite a few earlier tags), soil_water_retention_curve_type is set up to (theoretically) support having multiple parameterizations via polymorphism. However, we haven't plugged any alternatives in yet. Sean Swenson is about to plug in a skeleton for the VanGenuchten method, but it won't work at this point.

The main challenge to supporting alternative formulations is that there are implicit assumptions throughout the code that we are using the Clapp-Hornberg formulation. Specifically, there are a number of uses of bsw (and sucsat, which I think is also specific to Clapp-Hornberg). An additional challenge is that the two formulations have different inputs (e.g., bsw vs. nsw & msw), but that challenge may go away if we handle things correctly, as described below.

This page lays out a path forward.

### Overall strategy

The overall strategy for supporting multiple implementations will consist of the following steps; each step is described in more detail below:

1. Move all equations that implicitly assume use of Clapp-Hornberg into new methods in SoilWaterRetentionCurveClappHornberg1978Mod

2. Write alternative subroutines for the VanGenuchten method

3. Move variables that are specific to the soil water retention curve equations into soil_water_retention_curve_type or one of its subclasses, and change the interface to the subroutines in here accordingly

4. Handle initialization of bsw, etc.

5. Optional: Make routines work on vector inputs for efficiency

### 1. Move all equations that implicitly assume use of Clapp-Hornberg into new subroutines in SoilWaterRetentionCurveClappHornberg1978Mod

- Identify all uses of bsw in the code

- Identify all uses of sucsat in the code?

- I think, but am not sure, that sucsat is also specific to Clapp-Hornberg. Most, but not quite all, uses of sucsat currently are in equations that also reference bsw.

- If these uses are already captured by existing subroutines (soil_hk, soil_suction, soil_suction_inverse), then call these subroutines instead of duplicating them inline.

- If these uses are not captured by existing subroutines, then introduce new subroutines in SoilWaterRetentionCurveClappHornberg1978Mod, and use them.

One major challenge is how / where to set bsw (and sucsat). For now, leave that where it is (in SoilStateInitTimeConstMod); we'll move this code later.

appHornberg1978Mod, and use them.
One major challenge is how / where to set bsw (and sucsat). For now, leave that where it is (in SoilStateInitTimeConstMod); we'll move this code later.

### 2. Write alternative equations for the VanGenuchten method

The key point here is to identify:

- For Clapp-Hornberg subroutines that have a VanGenuchten alternative, how does the interface (specifically, input requirements) differ between the two methods?

- Which Clapp-Hornberg equations, if any, do not have a VanGenuchten alternative?

We need to do the following:

- For Clapp-Hornberg subroutines that have a VanGenuchten alternative: Write that alternative in SoilWaterRetentionCurveVanGenuchten1980Mod

- For Clapp-Hornberg subroutines that do not have a VanGenuchten alternative, we will need to make these available regardless of which parameterization we use.

- This will involve pulling up some subroutines into the soil_water_retention_curve_type base class (in SoilWaterRetentionCurveMod)

### 3. Move variables that are specific to the soil water retention curve equations into soil_water_retention_curve_type or one of its subclasses, and change the interface to the subroutines in here accordingly

There are a number of variables that are really specific to the soil water retention curve equations. Some or all of these are currently in soilstate_type. For example, this includes bsw and possibly sucsat.

Previous steps, laid out above, have ensured that all uses of these variables (e.g., bsw) have been moved into SoilWaterRetentionCurveMod or one of the subclasses (except for the initial setting of bsw and sucsat: we'll return to those in a bit).

So now, move bsw (and sucsat?) out of soilstate_type, and into soil_water_retention_curve_type or one of its subclasses.

If possible, we want to move these variables into the subclass that uses it – for example, moving bsw into soil_water_retention_curve_clapp_hornberg_1978_type. However, this may not be possible if we need to use some Clapp-Hornberg equations in all cases – i.e., if there were some Clapp-Hornberg subroutines that did not have a VanGenuchten alternative. In that case, we would need to leave bsw in the base class (soil_water_retention_curve_type).

Note that we will need InitAllocate and InitHistory routines for these new variables. The following patch – based off of https://svn-ccsm-models.cgd.ucar.edu/clm2/branch_tags/clm50hydro_tags/clm50hydro_n04_clm4_5_1_r111 ... note that this patch is untested!) illustrates the changes needed:

```diff
Index: src/biogeophys/SoilStateType.F90
===================================================================
--- src/biogeophys/SoilStateType.F90 (revision 72431)
+++ src/biogeophys/SoilStateType.F90 (working copy)
@@ -37,7 +37,6 @@
real(r8), pointer :: hk_l_col (:,:) ! col hydraulic conductivity (mm/s)
real(r8), pointer :: smp_l_col (:,:) ! col soil matric potential (mm)
real(r8), pointer :: smpmin_col (:) ! col restriction for min of soil potential (mm)
- real(r8), pointer :: bsw_col (:,:) ! col Clapp and Hornberger "b" (nlevgrnd)
real(r8), pointer :: watsat_col (:,:) ! col volumetric soil water at saturation (porosity)
real(r8), pointer :: watdry_col (:,:) ! col btran parameter for btran = 0
real(r8), pointer :: watopt_col (:,:) ! col btran parameter for btran = 1
@@ -132,7 +131,6 @@
allocate(this%smp_l_col (begc:endc,nlevgrnd)) ; this%smp_l_col (:,:) = nan
allocate(this%smpmin_col (begc:endc)) ; this%smpmin_col (:) = nan

- allocate(this%bsw_col (begc:endc,nlevgrnd)) ; this%bsw_col (:,:) = nan
allocate(this%watsat_col (begc:endc,nlevgrnd)) ; this%watsat_col (:,:) = nan
allocate(this%watdry_col (begc:endc,nlevgrnd)) ; this%watdry_col (:,:) = spval
allocate(this%watopt_col (begc:endc,nlevgrnd)) ; this%watopt_col (:,:) = spval
@@ -205,13 +203,6 @@
ptr_col=this%smp_l_col, set_spec=spval, l2g_scale_type='veg', default=active)

if (use_cn) then
- this%bsw_col(begc:endc,:) = spval
- call hist_addfld2d (fname='bsw', units='unitless', type2d='levgrnd', &
- avgflag='A', long_name='clap and hornberger B', &
- ptr_col=this%bsw_col, default='inactive')
- end if
-
- if (use_cn) then
this%rootfr_patch(begp:endp,:) = spval
call hist_addfld2d (fname='ROOTFR', units='proportion', type2d='levgrnd', &
avgflag='A', long_name='fraction of roots in each soil layer', &
Index: src/biogeophys/SoilWaterRetentionCurveClappHornberg1978Mod.F90
===================================================================
--- src/biogeophys/SoilWaterRetentionCurveClappHornberg1978Mod.F90 (revision 72431)
+++ src/biogeophys/SoilWaterRetentionCurveClappHornberg1978Mod.F90 (working copy)
@@ -19,6 +19,7 @@
soil_water_retention_curve_clapp_hornberg_1978_type
private
contains
+ procedure :: Init
procedure :: soil_hk ! compute hydraulic conductivity
procedure :: soil_suction ! compute soil suction potential
procedure :: soil_suction_inverse ! compute relative saturation at which soil suction is equal to a target value
@@ -39,9 +40,35 @@
! For now, this is simply a place-holder.
!-----------------------------------------------------------------------

+ ! Eventually this should call the Init routine (or replace the Init routine
+ ! entirely). But I think it would be confusing to do that until we switch everything
+ ! to use a constructor rather than the init routine.
+
end function constructor

!-----------------------------------------------------------------------
+ subroutine Init(this, bounds)
+ !
+ ! !DESCRIPTION:
+ ! Initialize soil_water_retention_curve_clapp_hornberg_1978_type object
+ !
+ ! !USES:
+ !
+ ! !ARGUMENTS:
+ class(soil_water_retention_curve_clapp_hornberg_1978_type), intent(inout) :: this
+ type(bounds_type), intent(in) :: bounds
+ !
+ ! !LOCAL VARIABLES:
+
+ character(len=*), parameter :: subname = 'Init'
+ !-----------------------------------------------------------------------
+
+ call this%InitBase(bounds)
+
+ end subroutine Init
+
+
+ !-----------------------------------------------------------------------
subroutine soil_hk(this, c, j, s, imped, soilstate_inst, hk, dhkds)
!
! !DESCRIPTION:
Index: src/biogeophys/SoilWaterRetentionCurveMod.F90
===================================================================
--- src/biogeophys/SoilWaterRetentionCurveMod.F90 (revision 72431)
+++ src/biogeophys/SoilWaterRetentionCurveMod.F90 (working copy)
@@ -5,6 +5,10 @@
! Abstract base class for functions to compute soil water retention curve
!
! !USES:
+ use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
+ use clm_varpar , only : nlevgrnd
+ use clm_varcon , only : spval
+ !
implicit none
save
private
@@ -14,7 +18,34 @@

type, abstract :: soil_water_retention_curve_type
private
+
+ ! NOTE(wjs, 2015-08-19) bsw is really specific to the Clapp-Hornberg method.
+ ! However, there are a number of places in the code that use it directly - i.e.,
+ ! implicitly assuming that we're using Clapp-Hornberg. Thus, for now we need bsw to
+ ! be available regardless of which soil water retention curve implementation is being
+ ! used.
+ !
+ ! Ideally, all uses of bsw would be removed from elsewhere in the code, so that it is
+ ! just used in the Clapp-Hornberg implementation. If that were done, then bsw could
+ ! be pushed down into SoilWaterRetentionCurveClappHornberg1978Mod.
+ !
+ ! If that ideal isn't possible (e.g., because we can't come up with a VanGenuchten
+ ! version of some of the methods), then we could at least make it more clear what's
+ ! going on by adding new subroutines to compute things based on bsw which are
+ ! currently being computed elsewhere. Those would go in this base class. Then bsw
+ ! would only be used in this base class and in the ClappHornberg subclass. Achieving
+ ! the ideal would then mean pushing these methods (along with bsw) down into
+ ! ClappHornberg, and introducing parallel methods in the other subclasses.
+
+ real(r8), pointer, public :: bsw_col(:,:) ! col Clapp and Hornberger "b" (nlevgrnd)
+
contains
+ ! ------------------------------------------------------------------------
+ ! The following routines need to be implemented by all type extensions
+ ! ------------------------------------------------------------------------
+
+ procedure(Init_interface), deferred :: Init
+
! compute hydraulic conductivity
procedure(soil_hk_interface), deferred :: soil_hk

@@ -23,25 +54,27 @@

! compute relative saturation at which soil suction is equal to a target value
procedure(soil_suction_inverse_interface), deferred :: soil_suction_inverse
+
+ ! ------------------------------------------------------------------------
+ ! The following routine should only be called by extensions of this base type
+ ! ------------------------------------------------------------------------
+ procedure, public :: InitBase ! initialize the base type
+
+ procedure, private :: InitAllocate
+ procedure, private :: InitHistory
+
end type soil_water_retention_curve_type

abstract interface

- ! Note: The following interfaces are set up based on the arguments needed for the
- ! clapphornberg1978 implementations. It's likely that these interfaces are not
- ! totally general for all desired implementations. In that case, we'll need to think
- ! about how to support different interfaces. Some possible solutions are:
- !
- ! - Make the interfaces contain all possible inputs that are needed by any
- ! implementation; each implementation will then ignore the inputs it doesn't need.
- !
- ! - For inputs that are needed only by particular implementations - and particularly
- ! for inputs that are constant in time (e.g., this is the case for bsw, I think):
- ! pass these into the constructor, and save pointers to these inputs as components
- ! of the child type that needs them. Then they aren't needed as inputs to the
- ! individual routines, allowing the interfaces for these routines to remain more
- ! consistent between different implementations.
+ subroutine Init_interface(this, bounds)
+ use decompMod, only : bounds_type
+ import :: soil_water_retention_curve_type

+ class(soil_water_retention_curve_type), intent(inout) :: this
+ type(bounds_type), intent(in) :: bounds
+ end subroutine Init_interface
+
subroutine soil_hk_interface(this, c, j, s, imped, soilstate_inst, &
hk, dhkds)

@@ -108,4 +141,80 @@

end interface

+contains
+
+ !-----------------------------------------------------------------------
+ subroutine InitBase(this, bounds)
+ !
+ ! !DESCRIPTION:
+ ! Initialize base soil_water_retention_curve_type object
+ !
+ ! !USES:
+ !
+ ! !ARGUMENTS:
+ class(soil_water_retention_curve_type), intent(inout) :: this
+ type(bounds_type), intent(in) :: bounds
+ !
+ ! !LOCAL VARIABLES:
+
+ character(len=*), parameter :: subname = 'InitBase'
+ !-----------------------------------------------------------------------
+
+ call this%InitAllocate(bounds)
+ call this%InitHistory(bounds)
+
+ end subroutine InitBase
+
+
+ !-----------------------------------------------------------------------
+ subroutine InitAllocate(this, bounds)
+ !
+ ! !DESCRIPTION:
+ ! Initialize data structure
+ !
+ ! !USES:
+ !
+ ! !ARGUMENTS:
+ class(soil_water_retention_curve_type), intent(inout) :: this
+ type(bounds_type), intent(in) :: bounds
+ !
+ ! !LOCAL VARIABLES:
+ integer :: begc, endc
+
+ character(len=*), parameter :: subname = 'InitAllocate'
+ !-----------------------------------------------------------------------
+
+ begc = bounds%begc; endc= bounds%endc
+
+ allocate(this%bsw_col(begc:endc,nlevgrnd)) ; this%bsw_col(:,:) = nan
+
+ end subroutine InitAllocate
+
+ !-----------------------------------------------------------------------
+ subroutine InitHistory(this, bounds)
+ !
+ ! !DESCRIPTION:
+ !
+ !
+ ! !USES:
+ use histFileMod , only: hist_addfld1d, hist_addfld2d, no_snow_normal
+ !
+ ! !ARGUMENTS:
+ class(soil_water_retention_curve_type), intent(inout) :: this
+ type(bounds_type), intent(in) :: bounds
+ !
+ ! !LOCAL VARIABLES:
+
+ character(len=*), parameter :: subname = 'InitHistory'
+ !-----------------------------------------------------------------------
+
+ this%bsw_col(begc:endc,:) = spval
+ call hist_addfld2d (fname='bsw', units='unitless', type2d='levgrnd', &
+ avgflag='A', long_name='clap and hornberger B', &
+ ptr_col=this%bsw_col, default='inactive')
+
+ end subroutine InitHistory
+
+
+
end module SoilWaterRetentionCurveMod
```

We also need to insert a call to the Init method; this would best be done in SoilWaterRetentionCurveFactoryMod, following the example of what's done for Ozone (where OzoneFactoryMod calls the Init method after creating the object).

Once these variables have been moved, the interface for the subroutines in soil_water_retention_curve_type can be cleaned up. These will no longer need bsw (and possibly sucsat) passed as arguments, because now they will be obtained from 'this'. If there are still any inputs needed from soilstate_inst at this point, I would recommend passing them directly, rather than passing the whole soilstate_inst (e.g., the caller will pass soilstate_inst%hksat_col(c,j) rather than passing soilstate_inst). That will make the interface more consistent (all inputs are passed as scalars), more clear (it's apparent what the inputs are to each routine), and will remove the dependence on soilStateType.

If you want to compile and test at this point, you will need to pass the soil_water_retention_curve_type object into the routine in SoilStateInitTimeConstMod, so that bsw (and sucsat, etc.?) can be set there... we'll fix that in a later step.

### 4. Handle initialization of bsw, etc.

Currently, bsw (and sucsat) are initialized in SoilStateInitTimeConstMod, partly via a call to pedotransf (in FuncPedotransferMod). Once bsw is moved into SoilWaterRetentionCurveMod, its initialization should really be done there, via a routine like set_texture_related_params(sand, clay): this would take sand and clay as inputs, and set whatever variables are needed by that class (bsw, sucsat, nsw, msw, etc.). There would likely be a separate version of this for each subclass - e.g., the VanGenuchten version would set nsw and msw but (maybe) not bsw and sucsat. (Or, if bsw and sucsat are still needed by some routines, then the VanGenuchten method could set nsw and msw and then call the base class's version of set_texture_related_params.)

The main question I have about this is how this initialization should interact with pedotransf. This is partly a science question: Is pedotransf entirely specific to Clapp-Hornberg, or are parts of it (e.g., the setting of watsat) more general? If this is really specific to Clapp-Hornberg, then its functionality should be moved into SoilWaterRetentionCurveClappHornberg1978Mod, or at least used directly by that, rather than used from some higher-level place. (Although, if bsw is needed in the base class, then we would also need to invoke pedotransf from the base class, rather than pushing it down into the Clapp-Hornberg subclass.) On the other hand, if parts are general and parts are specific to Clapp-Hornberg, then these general vs. specific parts could be split up – e.g., there could be a subroutine that sets bsw (and sucsat?), and another that sets the general parts; the subroutine to set bsw would be called from SoilWaterRetentionCurveClappHornberg1978Mod, whereas the subroutine that sets the general parts would still be called from SoilStateInitTimeConstMod. (Alternatively, we could keep the single pedotransf function that sets everything that depends on soil texture, and then pass the computed bsw, etc. to set_texture_related_params - which could use or ignore whichever parameters it wanted, based on what that subclass needed. If we have time-varying parameters, that would be nice in that we would only need to make one call when sand and clay are updated. But that puts a burden on pedotransf to compute all texture-related parameters regardless of the soil water retention curve implementation: e.g., nsw and msw as well as bsw, and any other parameters that are needed by methods added in the future. That feels awkward, and it's hard to imagine that the same pedotransfer function options would exist for the VanGenuchten method as exist for the ClappHornberg method - e.g., table4, table5, etc.)

Once that is determined, the main remaining issue is the order dependence of the setting of various variables. bsw depends on sand, clay, om_frac and om_b, which are set in SoilStateInitTimeConstMod. So the setting of bsw would need to be done after SoilStateInitTimeConst. For that to work, om_frac and om_b should be turned into arrays - or else the new set_texture_related_params would need to be called once per point, rather than on a whole vector (which may be okay). Currently, there are other variables set later in SoilStateInitTimeConst that depend on bsw. These equations will likely have been moved into SoilWaterRetentionCurveMod (as per the above tasks), but there may be a need to split up SoilStateInitTimeConst into multiple pieces: the prerequisites for computing bsw (etc.) and then later computing things that depend on having a fully initialized SoilWaterRetentionCurveMod.

As a side-note: note that currently bsw is set via pedotransf, but then that value is ignored a few lines later:

```Fortran
soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b
```

Sean Swenson points out that the second parenthesized term should be replaced with `bsw_col(c,lev)` - i.e., the value already computed by pedotransf. Note that this change may change answers due to compiler optimizations, so should probably be done separately from the other refactoring.

If we (eventually) make bsw time-varying, the situation doesn't change much – the main thing that changes is that the call to set bsw (etc.) will need to be done once per time step rather than just once at initialization. Then the CLM driver would call soil_water_retention_curve_inst%set_texture_related_params(sand, clay) at some point in the driver loop, and users of the other subroutines in this class would not need to know that now bsw (etc.) is time-varying (those other users would never access bsw (etc.) directly).

**However, there may also be places in the code that assume that bsw (and sucsat) are constant in time. For example, IrrigationMod calls soil_suction_inverse in IrrigationInitCold; that call would need to be moved into CalcIrrigationNeeded. We should search through the code for all calls made to routines in soil_water_retention_curve_type; if any of these are called in initialization, they should be moved into the run loop.** It would be best to do this move before actually making bsw time-varying, so that we can confirm that moving the code into the run loop is done correctly, by checking that it is bit-for-bit.

### 5. Optional: Make routines work on vector inputs for efficiency
Currently, the routines in soil_water_retention_curve_type operate on scalars, and are called from within loops. This isn't good for efficiency, and is particularly bad for vectorization. If this seems to be a performance issue (which can be determined by profiling), then we should rework these routines so that they operate on vectors of points. In at least some cases, they could even be reworked to operate on 2-d arrays (c,j), but that might lead to a loss of generality, and may not be as important for performance.

In some cases, this will mean splitting some loops: currently, some loops look like:

```
do ...

Do some stuff...

call soil_water_retention_curve_inst%some_routine(scalar inputs...)

Do more stuff...

end do
```

These loops would need to be split into:

```
do ...

Do some stuff...

end do

call soil_water_retention_curve_inst%some_routine(vector inputs...)

do ...

Do more stuff...

end do
```

For this to work, temporary scalars in these loops would need to be turned into vectors.

This transformation will aid vectorizability.