Add CLM ET formulation improvements: PFT photosynthesis, canopy interception, Medlyn stomata (parflow#712)

## Summary

Fixes CLM3-era ET suppression by addressing four high-priority gaps
identified in the gap analysis (Jefferson et al. 2017, Lawrence et al.
2019):
- **PFT-dependent photosynthesis:** Vcmax 39–62 for forests (vs
hardcoded 33), C4 pathway for grasses/crops, with improved stomata (dark
respiration, smooth colimitation, niter=5). Activated automatically when
`vcmx25` is present in `drv_vegp.dat`.
- **Configurable canopy interception:** `InterceptionFpiMax` (default
0.25) and `FwetExponent` (default 0.6667) keys for CLM3 scheme tuning.
- **CLM5 tanh canopy interception:** `InterceptionScheme=CLM5Tanh` with
`InterceptionTanhAlpha` replaces CLM3's 0.25 ceiling with `fpi = alpha *
tanh(LAI+SAI)`, approaching 1.0 for LAI>2.
- **Medlyn stomatal conductance:** `StomataScheme=Medlyn` uses VPD-based
conductance (Medlyn et al. 2011) with PFT-dependent `g1_medlyn` values
in `drv_vegp.dat`.

All changes are backward-compatible: default key values and existing
`drv_vegp.dat` files reproduce original CLM3 behavior bit-for-bit.

  ## New PF Keys (5)
  | Key | Default | Description |
  |-----|---------|-------------|
| `Solver.CLM.InterceptionFpiMax` | 0.25 | Max interception fraction
(CLM3 scheme) |
  | `Solver.CLM.FwetExponent` | 0.6667 | Wet canopy fraction exponent |
  | `Solver.CLM.StomataScheme` | BallBerry | `BallBerry` or `Medlyn` |
  | `Solver.CLM.InterceptionScheme` | CLM3 | `CLM3` or `CLM5Tanh` |
| `Solver.CLM.InterceptionTanhAlpha` | 1.0 | CLM5 tanh scaling
coefficient |

  ## New drv_vegp.dat parameters (7)

`vcmx25`, `c3psn`, `mp`, `bp`, `qe25`, `folnmx`, `g1_medlyn` — provided
via `drv_vegp_pft.dat`

Presence of `vcmx25` in vegp auto-sets `photosyn_custom=.true.` →
activates improved stomata physics (Rd, smooth colimitation, niter=5).
Old vegp files without these params → original CLM3 defaults →
bit-for-bit.

## Backward compatibility
- `photosyn_custom` boolean flag in `drv_readvegpf.F90` gates improved
stomata (not `isnan()`)
- Guarded fpi/fwet expressions use original Fortran literals at default
key values
  - Original CLM3 stomata loop preserved verbatim in separate code path
- `InterceptionScheme` defaults to `CLM3` → existing fpi code path
unchanged
- New struct fields at END of `clm1d` to preserve existing field offsets

  ## Commits

1. `c0001aaf` — Add CLM ET formulation improvements: PFT photosynthesis,
canopy interception, Medlyn stomata
  2. `9e09d14f` — Fix Black formatting in clm_et_interception.py
  3. `58065fa3` — Fix uninitialized folnvt on custom photosynthesis path
4. `f272dfda` — Add CLM5 tanh canopy interception scheme
(InterceptionScheme key)
  5. `b2969dbb` — Fix Uncrustify formatting in parflow_proto_f.h

  ## References

- Jefferson, J.L., Maxwell, R.M., and Constantine, P.G. (2017).
Exploring the sensitivity of photosynthesis and
stomatal resistance parameters in a land surface model. *J.
Hydrometeorol.* 18, 897–915.
- Lawrence, D.M., et al. (2019). The Community Land Model version 5.
*JAMES*, 11, 4245–4287.
- Medlyn, B.E., et al. (2011). Reconciling the optimal and empirical
approaches to modelling stomatal conductance. *Global Change Biol.* 17,
2134–2144.
- Collatz, G.J., et al. (1991). Physiological and environmental
regulation of stomatal conductance, photosynthesis and transpiration.
*Agric. For. Meteorol.* 54, 107–136.

  ## Test plan

- [x] All existing CLM tests pass bit-for-bit (clm, clm_snow_defaults, 8
snow tests)
- [x] 4 new CI tests: `clm_et_pft_photosynthesis`,
`clm_et_interception`, `clm_et_medlyn`, `clm_et_tanh_interception`
  - [x] Validated against 8 AmeriFlux tower sites (latent heat flux)
  - [x] Uncrustify (C) and Black (Python) formatting verified
  - [x] CI pipeline passes on Linux

  🤖 Generated with [Claude Code](https://claude.com/claude-code)

Create at:
parflow/parflow@master...feature/clm_et

---------

Co-authored-by: Reed <rmaxwell@mines.edu>
Co-authored-by: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Co-authored-by: Steven Smith <smith84@llnl.gov>

Author: 4 people

docs/user_manual/keys.rst

@@ -6839,6 +6839,88 @@ diurnal artifacts in fractional snow cover. Typical range 48-96 hours.
       pfset Solver.CLM.FracSnoAvgWindow 72.0           ## TCL syntax
       <runname>.Solver.CLM.FracSnoAvgWindow = 72.0     ## Python syntax
 
+**ET Formulation Improvements**
+
+These keys improve CLM's evapotranspiration calculations. Default values
+reproduce original CLM3 behavior for backward compatibility. PFT-dependent
+photosynthesis parameters (vcmx25, c3psn, mp, bp, qe25, folnmx) can be
+provided in drv_vegp.dat; when detected, improved stomata physics
+(dark respiration, smooth colimitation, niter=5) activates automatically.
+
+*double* **Solver.CLM.InterceptionFpiMax** 0.25 Maximum interception
+fraction coefficient. CLM3 default 0.25 caps canopy interception at 25%
+even for dense canopies. CLM5 uses higher values.
+Formula: fpi = InterceptionFpiMax * (1 - exp(-0.5*(LAI+SAI)))
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.InterceptionFpiMax 0.25           ## TCL syntax
+      <runname>.Solver.CLM.InterceptionFpiMax = 0.25     ## Python syntax
+
+*double* **Solver.CLM.FwetExponent** 0.6667 Power-law exponent for wet
+canopy fraction. CLM default 2/3. Lower values keep the canopy wet longer,
+increasing wet canopy evaporation.
+Formula: fwet = (h2ocan/(dewmx*LAI))^FwetExponent
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.FwetExponent 0.6667               ## TCL syntax
+      <runname>.Solver.CLM.FwetExponent = 0.6667         ## Python syntax
+
+*string* **Solver.CLM.StomataScheme** BallBerry Selects the stomatal
+conductance model.
+
+**BallBerry**:
+   CLM3 default. Uses relative humidity in the conductance equation.
+
+**Medlyn**:
+   Medlyn et al. (2011) model. Uses vapor pressure deficit (VPD) instead
+   of relative humidity. Requires g1_medlyn parameter in drv_vegp.dat for
+   PFT-dependent slope values.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.StomataScheme BallBerry            ## TCL syntax
+      <runname>.Solver.CLM.StomataScheme = "BallBerry"   ## Python syntax
+
+*string* **Solver.CLM.InterceptionScheme** CLM3 Selects the canopy
+interception scheme.
+
+**CLM3**:
+   CLM3 default exponential formulation.
+   Formula: fpi = InterceptionFpiMax * (1 - exp(-0.5*(LAI+SAI)))
+   Asymptotes to InterceptionFpiMax (default 0.25) regardless of LAI.
+
+**CLM5Tanh**:
+   CLM5 tanh formulation (Lawrence et al. 2019).
+   Formula: fpi = InterceptionTanhAlpha * tanh(LAI+SAI)
+   Approaches InterceptionTanhAlpha (default 1.0) for LAI > 2, correcting
+   the CLM3 structural low bias in canopy interception for dense canopies.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.InterceptionScheme CLM3              ## TCL syntax
+      <runname>.Solver.CLM.InterceptionScheme = "CLM3"     ## Python syntax
+
+*double* **Solver.CLM.InterceptionTanhAlpha** 1.0 Scaling coefficient for
+CLM5 tanh interception scheme. Formula: fpi = alpha * tanh(LAI+SAI).
+Default 1.0 matches CLM5. Only used when InterceptionScheme is CLM5Tanh.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.InterceptionTanhAlpha 1.0             ## TCL syntax
+      <runname>.Solver.CLM.InterceptionTanhAlpha = 1.0      ## Python syntax
+
 
 .. _ParFlow NetCDF4 Parallel I/O:
 

pf-keys/definitions/solver.yaml

@@ -987,6 +987,72 @@ Solver:
           max_value: 1.0
         RequiresModule: CLM
 
+    # ET formulation improvements @RMM 2026
+    InterceptionFpiMax:
+      help: >
+        [Type: double] Maximum interception fraction coefficient. CLM3 default
+        0.25 caps canopy interception at 25% even for dense canopies. CLM5 uses
+        higher values. Formula: fpi = InterceptionFpiMax * (1 - exp(-0.5*(LAI+SAI)))
+      default: 0.25
+      domains:
+        DoubleValue:
+          min_value: 0.05
+          max_value: 1.0
+        RequiresModule: CLM
+
+    FwetExponent:
+      help: >
+        [Type: double] Power-law exponent for wet canopy fraction. CLM default
+        2/3 (0.6667). Lower values keep the canopy wet longer, increasing wet
+        canopy evaporation. Formula: fwet = (h2ocan/(dewmx*LAI))^FwetExponent
+      default: 0.6667
+      domains:
+        DoubleValue:
+          min_value: 0.3
+          max_value: 1.0
+        RequiresModule: CLM
+
+    StomataScheme:
+      help: >
+        [Type: string] Stomatal conductance model. BallBerry uses relative
+        humidity (CLM3 default). Medlyn uses vapor pressure deficit and is
+        recommended for modern applications. Reference: Medlyn et al. (2011)
+        Global Change Biology.
+      default: BallBerry
+      domains:
+        EnumDomain:
+          enum_list:
+            - BallBerry
+            - Medlyn
+        RequiresModule: CLM
+
+    InterceptionScheme:
+      help: >
+        [Type: string] Canopy interception scheme. CLM3 uses exponential
+        formulation capped at InterceptionFpiMax (default 0.25). CLM5Tanh
+        uses tanh(LAI+SAI) scaled by InterceptionTanhAlpha, which approaches
+        1.0 for LAI>2 and corrects the CLM3 structural low bias.
+        Reference: Lawrence et al. (2019) JAMES.
+      default: CLM3
+      domains:
+        EnumDomain:
+          enum_list:
+            - CLM3
+            - CLM5Tanh
+        RequiresModule: CLM
+
+    InterceptionTanhAlpha:
+      help: >
+        [Type: double] Scaling coefficient for CLM5 tanh interception scheme.
+        Formula: fpi = alpha * tanh(LAI+SAI). Default 1.0 matches CLM5.
+        Only used when InterceptionScheme is CLM5Tanh.
+      default: 1.0
+      domains:
+        DoubleValue:
+          min_value: 0.1
+          max_value: 1.0
+        RequiresModule: CLM
+
     # -----------------------------------------------------------------------------
     # CLM input variables to write drv_clmin.dat (Input) and drv_vegp.dat (VegParams) files
     # -----------------------------------------------------------------------------

pfsimulator/clm/clm.F90

@@ -20,11 +20,13 @@ subroutine clm_lsm(pressure,saturation,evap_trans,top,bottom,porosity,pf_dz_mult
 frac_sno_typepf,frac_sno_roughnesspf,frac_sno_roughness_minpf,frac_sno_roughness_maxpf,        &
 frac_sno_gamma_szapf,frac_sno_tau_szapf,                                                          &
 snowage_tau0_vispf,snowage_tau0_nirpf,snowage_grain_growth_vispf,snowage_grain_growth_nirpf,        &
-snowage_dirt_soot_vispf,snowage_dirt_soot_nirpf,snowage_reset_factorpf)
+snowage_dirt_soot_vispf,snowage_dirt_soot_nirpf,snowage_reset_factorpf,                                    &
+interception_fpi_maxpf,fwet_exponentpf,stomata_schemepf,                                     &
+interception_schemepf,interception_tanh_alphapf)
 
   !=========================================================================
   !
-  !  CLMCLMCLMCLMCLMCLMCLMCLMCL  A community developed and sponsored, freely   
+  !  CLMCLMCLMCLMCLMCLMCLMCLMCL  A community developed and sponsored, freely
   !  L                        M  available land surface process model.  
   !  M --COMMON LAND MODEL--  C  	
   !  C                        L  CLM WEB INFO: http://clm.gsfc.nasa.gov
@@ -206,6 +208,13 @@ subroutine clm_lsm(pressure,saturation,evap_trans,top,bottom,porosity,pf_dz_mult
   real(r8) :: snowage_dirt_soot_nirpf            ! NIR dirt/soot factor [-]
   real(r8) :: snowage_reset_factorpf             ! fresh snow reset factor [-]
 
+  ! ET formulation improvements @RMM 2026
+  real(r8) :: interception_fpi_maxpf             ! max interception fraction coeff [-]
+  real(r8) :: fwet_exponentpf                    ! power-law exponent for wet canopy [-]
+  integer  :: stomata_schemepf                   ! stomatal model: 0=BallBerry, 1=Medlyn
+  integer  :: interception_schemepf              ! interception: 0=CLM3, 1=CLM5Tanh
+  real(r8) :: interception_tanh_alphapf          ! CLM5 tanh scaling coeff [-]
+
   ! local indices & counters
   integer  :: i,j,k,k1,j1,l1                     ! indices for local looping
   integer  :: bj,bl                              ! indices for local looping !BH
@@ -596,6 +605,13 @@ subroutine clm_lsm(pressure,saturation,evap_trans,top,bottom,porosity,pf_dz_mult
            clm(t)%snowage_dirt_soot_nir   = snowage_dirt_soot_nirpf
            clm(t)%snowage_reset_factor    = snowage_reset_factorpf
 
+           ! for ET formulation improvements @RMM 2026
+           clm(t)%interception_fpi_max = interception_fpi_maxpf
+           clm(t)%fwet_exponent        = fwet_exponentpf
+           clm(t)%stomata_scheme       = stomata_schemepf
+           clm(t)%interception_scheme  = interception_schemepf
+           clm(t)%interception_tanh_alpha = interception_tanh_alphapf
+
            ! for irrigation
            clm(t)%irr_type           = irr_typepf
            clm(t)%irr_cycle          = irr_cyclepf

pfsimulator/clm/clm_hydro_canopy.F90

@@ -111,7 +111,17 @@ subroutine clm_hydro_canopy (clm)
 
         ! Direct throughfall
 
-        fpi = 0.25d0*(1. - exp(-0.5*(clm%elai + clm%esai)))
+        if (clm%interception_scheme == 1) then
+           ! CLM5 tanh interception (Lawrence et al. 2019)
+           fpi = clm%interception_tanh_alpha                             &
+                 * tanh(clm%elai + clm%esai)
+        else if (clm%interception_fpi_max /= 0.25d0) then
+           fpi = clm%interception_fpi_max                                &
+                 * (1.d0 - exp(-0.5d0*(clm%elai + clm%esai)))
+        else
+           ! Original CLM3 expression for bit-for-bit backward compat
+           fpi = 0.25d0*(1. - exp(-0.5*(clm%elai + clm%esai)))
+        endif
         qflx_through  = prcp*(1.-fpi)*clm%frac_veg_nosno
 
         ! Water storage of intercepted precipitation and dew
@@ -313,7 +323,11 @@ subroutine clm_hydro_canopy (clm)
   if (clm%h2ocan > 0. .and. clm%frac_veg_nosno == 1) then
      vegt     = clm%frac_veg_nosno*(clm%elai + clm%esai)
      dewmxi   = 1.0/clm%dewmx
-     clm%fwet = ((dewmxi/vegt)*clm%h2ocan)**.666666666666
+     if (clm%fwet_exponent /= 0.6667d0) then
+        clm%fwet = ((dewmxi/vegt)*clm%h2ocan)**clm%fwet_exponent
+     else
+        clm%fwet = ((dewmxi/vegt)*clm%h2ocan)**.666666666666
+     endif
      clm%fwet = min(clm%fwet,dble(1.0))     ! Check for maximum limit of fwet
      clm%fdry = (1.-clm%fwet)*clm%elai/(clm%elai+clm%esai)
   else