diff --git a/docs/_toc.yml b/docs/_toc.yml
index 1ff64785..2dc70d6a 100644
--- a/docs/_toc.yml
+++ b/docs/_toc.yml
@@ -22,6 +22,7 @@ parts:
       - file: users_guide/running_tsmp_pdaf/input_enkfpf
       - file: users_guide/running_tsmp_pdaf/tsmp_pdaf_omi
       - file: users_guide/running_tsmp_pdaf/grace_tws_da
+      - file: users_guide/running_tsmp_pdaf/lai_da
 
     - file: users_guide/misc/README
       title: Misc
diff --git a/docs/users_guide/running_tsmp_pdaf/README.md b/docs/users_guide/running_tsmp_pdaf/README.md
index 18dfd8f0..fd8753e1 100644
--- a/docs/users_guide/running_tsmp_pdaf/README.md
+++ b/docs/users_guide/running_tsmp_pdaf/README.md
@@ -14,6 +14,8 @@ COSMO](cos)). Additionally, a control file for the data assimilation
 ([Observation files](obs)) need to be present in the run directory.
 Furthermore, some command line options ([Command line options](cmd))
 need to be specified when TSMP-PDAF is executed.
+For assimilation of LAI observations into eCLM, see
+[LAI Data Assimilation](laida).
 
 For assimilation of GRACE/GRACE-FO terrestrial water storage
 observations into eCLM, see [GRACE/TWS Data Assimilation](gracetws).
diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
index 78c83b3e..b4aba368 100644
--- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
+++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
@@ -54,6 +54,9 @@ nprocs      =
 update_swc  =
 update_texture  =
 update_T  =
+update_lai  =
+update_lai_params =
+update_lai_incr_w =
 print_swc   =
 print_et   =
 print_inc   =
@@ -499,6 +502,89 @@ Currently only CLM3.5
 
 -  1: Update of ground and vegetation temperature
 
+(enkfpf:clm:update_lai)=
+### CLM:update_lai ###
+
+`CLM:update_lai`: (integer) Flag for assimilation of Leaf Area Index
+(LAI) in eCLM. Default: `0`.
+
+Only available for eCLM with Biogeochemical Cycling (BGC) enabled.
+See [LAI Data Assimilation](laida) for algorithmic details.
+
+-  0: No LAI assimilation.
+
+-  1: Gridcell-averaged LAI state vector. One LAI value per grid cell
+   is stored in the state vector, computed as a patch-weight-averaged
+   LAI derived from per-patch leaf carbon (`leafc`) and the specific
+   leaf area parameters (`slatop`, `dsladlai`) using the formula of
+   Thornton and Zimmermann (2007, J. Clim., 20, 3902–3923). After the
+   PDAF update, the updated gridcell LAI is back-transformed to
+   per-patch `leafc` and `leafn`.
+
+-  2: (WARNING: This option is not tested. Rather use option 3 for
+   patch-level state vector) Patch-level state vector. `leafc`,
+   `slatop`, and `dsladlai` are stored for every patch (state vector
+   size: 3 × number of patches).  The LAI-to-observation mapping is
+   computed inside the observation operator `obs_op_pdaf`. The
+   `slatop` and `dsladlai` portions of the state vector can optionally
+   be estimated jointly; set `CLM:update_lai_params=2` to write them
+   back to eCLM after the update.
+
+-  3: Carbon-pool state vector. `leafc`, `livestemc`, and `deadstemc`
+   are stored for every vegetated patch (bare-ground PFT excluded;
+   state vector size: 3 × number of vegetated patches). The
+   observation operator `obs_op_pdaf` maps the `leafc` block to
+   observed gridcell LAI using the Thornton and Zimmermann (2007)
+   formula. After the PDAF update, the three carbon pools are written
+   back to eCLM and the corresponding nitrogen pools (`leafn`,
+   `livestemn`, `deadstemn`) are updated consistently using PFT C:N
+   ratios. Set `CLM:update_lai_params=1` to also estimate `slatop` and
+   `medlynslope` jointly.
+
+(enkfpf:clm:update_lai_params)=
+### CLM:update_lai_params ###
+
+`CLM:update_lai_params`: (integer) Flag for joint
+state–parameter estimation of specific leaf area parameters alongside
+LAI. Only takes effect if `CLM:update_lai` is non-zero. Default: `0`.
+
+-  0: No parameter estimation.
+
+-  1: For `CLM:update_lai=1`: appends `slatop` (specific leaf area at
+   the canopy top) for each patch to the state vector as an additional
+   parameter block.
+
+   For `CLM:update_lai=3`: appends `slatop` and `medlynslope`
+   (Medlyn stomatal slope) per vegetated patch as two additional
+   parameter blocks (blocks 4 and 5). Both parameters are read from
+   the PFT constants before the PDAF update and written back after.
+   The observation operator uses the in-state `slatop` values for the
+   LAI forward mapping.
+
+-  2: Only for `CLM:update_lai=2`. After the PDAF update, writes the
+   updated `slatop` and `dsladlai` values from the state vector back to
+   the eCLM plant functional type constants.
+
+(enkfpf:clm:update_lai_incr_w)=
+### CLM:update_lai_incr_w ###
+
+`CLM:update_lai_incr_w`: (real) Weight of the **additive** LAI increment when mapping the gridcell analysis back to patches. Only used if `CLM:update_lai` is non-zero. Default: `1.0`. 
+
+After the EnKF update, the gridcell-mean LAI increment is applied to each patch in two ways and then blended:
+
+- **Additive** (weight `w`): each patch receives the same increment
+  $\Delta = \text{LAI}_\text{a} - \text{LAI}_\text{f}$:
+  $\text{LAI}'(p) = \text{LAI}_\text{f}(p) + \Delta$.
+
+- **Multiplicative** (weight $1-w$): each patch is scaled by the gridcell analysis-to-forecast ratio:
+  $\text{LAI}'(p) = \text{LAI}_\text{f}(p)\,\text{LAI}_\text{a} / \text{LAI}_\text{f}$.
+
+The final patch LAI is
+$\text{LAI}'(p) = w\,\text{LAI}'_\text{add}(p) + (1-w)\,\text{LAI}'_\text{mult}(p)$.
+
+-  `1.0`: purely additive.
+-  `0.0`: purely multiplicative.
+
 ### CLM:print_swc ###
 
 `CLM:print_swc`: (integer) If set to `1`, the updated soil moisture
@@ -998,6 +1084,9 @@ Default: 0, output turned off.
  |           | `problemname`           | \-            |
  |           | `nprocs`                | 0             |
  |           | `update_swc`            | 1             |
+ |           | `update_lai`            | 0             |
+ |           | `update_lai_params`     | 0             |
+ |           | `update_lai_incr_w`     | 1.0           |
  |           | `print_swc`             | 0             |
  |           | `print_et`              | 0             |
  |           | `print_inc`             | 0             |
diff --git a/docs/users_guide/running_tsmp_pdaf/lai_da.md b/docs/users_guide/running_tsmp_pdaf/lai_da.md
new file mode 100644
index 00000000..06a71615
--- /dev/null
+++ b/docs/users_guide/running_tsmp_pdaf/lai_da.md
@@ -0,0 +1,251 @@
+(laida)=
+# LAI Data Assimilation #
+
+Leaf Area Index (LAI) data assimilation allows observed LAI fields to
+be used to update the vegetation state and, optionally, the specific
+leaf area parameters of eCLM.
+
+**eCLM only.** Requires BGC (Biogeochemical Mode) to be active.
+Not available for CLM3.5 or CLM5.0 without BGC.
+
+## Configuration ##
+
+LAI DA is controlled by three parameters in `enkfpf.par`:
+
+- [`CLM:update_lai`](enkfpf:clm:update_lai) — selects the state vector
+  layout and the location of the forward observation operator.
+- [`CLM:update_lai_params`](enkfpf:clm:update_lai_params) — enables
+  joint estimation of specific leaf area parameters alongside LAI.
+- [`CLM:update_lai_incr_w`](enkfpf:clm:update_lai_incr_w) — blend between
+  additive and multiplicative patch LAI increments (options 1 and 2 only).
+
+(laida:background)=
+## Background ##
+
+eCLM with BGC does not carry LAI as an independent prognostic
+variable.  Instead, LAI is a diagnostic quantity derived from leaf
+carbon (`leafc`) and plant-functional-type (PFT) constants for
+specific leaf area:
+
+$$
+\text{LAI}(p) =
+\begin{cases}
+  \dfrac{s_{\text{top}}(p)\,\bigl(\exp\!\bigl(C_\ell(p)\,s'(p)\bigr) - 1\bigr)}{s'(p)}
+    & s'(p) > 0 \\[6pt]
+  s_{\text{top}}(p)\, C_\ell(p) & s'(p) = 0
+\end{cases}
+$$
+
+where
+- $C_\ell$ is `leafc` (leaf carbon, gC m$^{-2}$)
+- $s_{\text{top}}$ is `slatop` (specific leaf area at the canopy top,
+m$^2$ gC$^{-1}$)
+- $s'$ is `dsladlai` (slope of specific leaf
+area with respect to LAI, m$^2$ gC$^{-1}$)
+
+Both $s_{\text{top}}$ and $s'$ are PFT constants read from the CLM
+surface dataset. The formula follows Eq. 3 of Thornton and Zimmermann
+(2007, J. Clim., 20, 3902–3923).
+
+Because LAI is not a direct model state, assimilation updates `leafc`
+(and, consistently, `leafn = leafc / leafcn`) rather than LAI
+itself.
+
+## State Vector ##
+
+The choice of [`CLM:update_lai`](enkfpf:clm:update_lai) determines the
+state vector layout:
+
+| `update_lai` | `update_lai_params` | State vector contents                                                           | Size        |
+|:------------:|:-------------------:|---------------------------------------------------------------------------------|-------------|
+|      1       |          0          | Gridcell-averaged LAI                                                           | $N_g$       |
+|      1       |          1          | Gridcell-averaged LAI + `slatop` per patch                                      | $N_g + N_p$ |
+|      2       |          0          | `leafc`, `slatop`, `dsladlai` per patch                                         | $3\,N_p$    |
+|      2       |          2          | `leafc`, `slatop`, `dsladlai` per patch (parameters updated after assimilation) | $3\,N_p$    |
+|      3       |          0          | `leafc`, `livestemc`, `deadstemc` per vegetated patch                           | $3\,N_v$    |
+|      3       |          1          | `leafc`, `livestemc`, `deadstemc` + `slatop`, `medlynslope` per vegetated patch | $5\,N_v$    |
+
+$N_g$ = number of local grid cells, $N_p$ = number of local patches,
+$N_v$ = number of local vegetated patches (non-bare-ground PFTs).
+
+## Option 1 — Gridcell LAI State Vector ##
+
+With `CLM:update_lai=1`, the state vector holds one LAI value per grid
+cell (a weighted average over all patches in that cell):
+
+**Set phase** (`set_clm_statevec`):
+
+1. Per-patch LAI is computed from `leafc` using the formula above and
+   clipped to zero.
+2. The gridcell LAI is the patch-weight-averaged sum:
+   $\text{LAI}_g = \sum_p w_p\,\text{LAI}(p)$, where $w_p$ is
+   `patch%wtgcell`.
+3. Per-patch forecast LAI $\text{LAI}_\text{f}(p)$ is stored for the
+   update phase;
+
+**PDAF update**: PDAF operates on the gridcell-level LAI values.
+
+**Update phase** (`update_clm`):
+
+1. The gridcell analysis LAI $\text{LAI}'_g$ is mapped to patches using
+   a blend of additive and multiplicative increments controlled by
+   [`CLM:update_lai_incr_w`](enkfpf:clm:update_lai_incr_w)
+   (default `1.0` = purely additive).
+2. The updated per-patch LAI is inverted to obtain the new `leafc`:
+
+$$
+     C'_\ell(p) =
+     \begin{cases}
+       \dfrac{\ln\!\bigl((\text{LAI}'(p)\,s'(p)\,/\,s_{\text{top}}(p)) + 1\bigr)}{s'(p)}
+         & s'(p) > 0 \\[6pt]
+       \text{LAI}'(p)\,/\,s_{\text{top}}(p) & s'(p) = 0
+     \end{cases}
+$$
+
+3. `leafc` is clipped to zero (no negative leaf carbon).
+4. `leafn` is set consistently: $N'_\ell = C'_\ell\,/\,\text{leafcn}$.
+
+**Joint parameter estimation** (`update_lai_params=1`): `slatop` for
+each patch is appended to the state vector and updated by PDAF. The
+updated `slatop` values are written back to the PFT constants before
+the LAI inversion in step 2.
+
+## Option 2 — Patch-Level State Vector with Observation Operator ##
+
+:::{warning}
+This option is not tested. Use [Option 3](laida:option3) for patch-level state vectors instead.
+:::
+
+With `CLM:update_lai=2`, the state vector holds `leafc`, `slatop`,
+and `dsladlai` for every patch in three consecutive blocks:
+
+$$
+\mathbf{x} = \bigl[\underbrace{C_\ell(1),\ldots,C_\ell(N_p)}_{\text{block 1}},
+              \underbrace{s_{\text{top}}(1),\ldots,s_{\text{top}}(N_p)}_{\text{block 2}},
+              \underbrace{s'(1),\ldots,s'(N_p)}_{\text{block 3}}\bigr]
+$$
+
+The observation operator (`obs_op_pdaf`, active when
+`CLM:update_lai=2`) maps this state to observed gridcell LAI:
+
+$$
+\mathcal{H}(\mathbf{x})_i = \sum_{p:\,g(p)=g_i} w_p\,\text{LAI}(p)
+$$
+
+where $g_i$ is the grid cell corresponding to observation $i$.
+
+After the PDAF update:
+
+- Patch LAI is updated with the same additive/multiplicative blend as in
+  option 1 (using weighted gridcell-mean LAI for the ratio and increment).
+- `leafc` is obtained by inverting the blended patch LAI and clipped to
+  zero; `leafn` is updated consistently.
+- If `CLM:update_lai_params=2`, `slatop` and `dsladlai` are also
+  read back from blocks 2 and 3 and written to the PFT constants.
+  The diagnostic LAI update is then left to eCLM's own phenology
+  routines.
+
+(laida:option3)=
+## Option 3 — Carbon-Pool State Vector ##
+
+With `CLM:update_lai=3`, the state vector holds three carbon pools for
+every vegetated patch (bare-ground PFT excluded):
+
+$$
+\mathbf{x} = \bigl[\underbrace{C_\ell(1),\ldots,C_\ell(N_v)}_{\text{block 1: leafc}},
+              \underbrace{C_s(1),\ldots,C_s(N_v)}_{\text{block 2: livestemc}},
+              \underbrace{C_d(1),\ldots,C_d(N_v)}_{\text{block 3: deadstemc}}\bigr]
+$$
+
+where $C_\ell$ = `leafc`, $C_s$ = `livestemc`, $C_d$ = `deadstemc`.
+
+The observation operator (`obs_op_pdaf`) maps this state to observed
+gridcell LAI using fixed `dsladlai` from the PFT constants and `slatop`
+either from the state vector (when `update_lai_params=1`) or from the PFT
+constants:
+
+$$
+\mathcal{H}(\mathbf{x})_i = \sum_{p:\,g(p)=g_i} w_p\,\text{LAI}(p)
+$$
+
+with $\text{LAI}(p)$ computed from $C_\ell(p)$ using the Thornton and
+Zimmermann (2007) formula described in the [Background](laida:background)
+section.
+
+**Set phase** (`set_clm_statevec`):
+
+1. `leafc`, `livestemc`, and `deadstemc` for each vegetated patch are
+   copied into the three consecutive state vector blocks.
+2. If `update_lai_params=1`, two additional blocks are appended:
+   `slatop` (block 4) and `medlynslope` (block 5) per patch, read from
+   the PFT constants for each patch's PFT type.
+
+**PDAF update**: PDAF operates on the full state vector.
+
+**Update phase** (`update_clm`):
+
+1. `leafc`, `livestemc`, and `deadstemc` are read from their respective
+   state vector blocks and clipped to zero (no negative carbon pools).
+2. Nitrogen pools are updated consistently using PFT C:N ratios from
+   `pftcon`:
+   - $N'_\ell = C'_\ell / \text{leafcn}$
+   - $N'_s = C'_s / \text{livewdcn}$
+   - $N'_d = C'_d / \text{deadwdcn}$
+3. If `update_lai_params=1`, `slatop` and `medlynslope` are read from
+   blocks 4 and 5 and written back to the PFT constants
+   (`pftcon%slatop`, `pftcon%medlynslope`).
+
+**Note**: Unlike options 1 and 2, the `CLM:update_lai_incr_w` blending
+parameter has no effect in option 3. LAI is not part of the state vector;
+only the underlying carbon pools are updated directly.
+
+## Configuration Examples ##
+
+**Example 1** — gridcell LAI assimilation, no parameter estimation:
+
+```ini
+[CLM]
+update_lai        = 1
+update_lai_params = 0
+update_lai_incr_w = 1.0
+```
+
+**Example 2** — gridcell LAI assimilation with joint `slatop` estimation:
+
+```ini
+[CLM]
+update_lai        = 1
+update_lai_params = 1
+```
+
+**Example 3** — patch-level state vector with custom observation operator:
+
+```ini
+[CLM]
+update_lai        = 2
+update_lai_params = 0
+```
+
+**Example 4** — patch-level state vector with joint `slatop`/`dsladlai` estimation:
+
+```ini
+[CLM]
+update_lai        = 2
+update_lai_params = 2
+```
+
+**Example 5** — carbon-pool state vector, no parameter estimation:
+
+```ini
+[CLM]
+update_lai        = 3
+update_lai_params = 0
+```
+
+**Example 6** — carbon-pool state vector with joint `slatop`/`medlynslope` estimation:
+
+```ini
+[CLM]
+update_lai        = 3
+update_lai_params = 1
+```
diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90
index bfd08a43..b7def5b6 100644
--- a/interface/framework/obs_op_pdaf.F90
+++ b/interface/framework/obs_op_pdaf.F90
@@ -57,6 +57,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 #endif
    USE mod_assimilation, ONLY: obs_interp_indices_p
    USE mod_assimilation, ONLY: obs_interp_weights_p
+   USE mod_parallel_pdaf, ONLY: mype_world
    use mod_tsmp, &
        only: obs_interp_switch, &
        soilay, &
@@ -71,6 +72,16 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 #if defined CLMSA
    USE enkf_clm_mod, &
         ONLY : clm_varsize, clm_paramarr, clmupdate_swc, clmupdate_T, clmcrns_bd
+   USE enkf_clm_mod, ONLY : clmupdate_lai
+   USE enkf_clm_mod, ONLY : clm_begp
+   USE enkf_clm_mod, ONLY : clm_endp
+   USE enkf_clm_mod, ONLY : clm_patch2gc
+   USE enkf_clm_mod, ONLY : clm_patchwt
+   USE enkf_clm_mod, ONLY : clm_lai_patch_itype
+   USE enkf_clm_mod, ONLY : clmt_printensemble
+   USE enkf_clm_mod, ONLY : clmupdate_lai_params
+   USE enkf_clm_mod, ONLY : lai_mode3_npool
+   USE pftconMod, only : pftcon
 #ifdef CLMFIVE
    USE clm_instMod, &
      ONLY : soilstate_inst
@@ -101,7 +112,9 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 !                0.8d0,   1.3d0,   2.d0,  3.d0, 5.d0,  12.d0/) !soil depth
 
 real :: tot, avesm, avesm_temp, Dp
+real :: avelai, slatop_val
 integer :: nsc
+character(len=48) :: fn_hx
 ! end of hcp
 
 #ifndef PARFLOW_STAND_ALONE
@@ -128,6 +141,79 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 ! If no special observation operator is compiled, use point observations
 lpointobs = .true.
 
+! LAI assimilation assuming statevec contains leafc, slatop, dsladlai for each patch
+! two component op to get to LAI obs
+! 1 calculate tlai from statevec per patch
+! state_p contains leafc of varsize then slatop of varsize and then dsadlai of varsize
+! 2 average patches according to weight per gridcell
+if (clmupdate_lai==2) then
+  lpointobs = .false. ! disable general obs_op
+  ! Note: This fails if there ever is an observation pointing to a gridcell with 0 patches
+  ! this should not happen, but would cause div by zero.
+
+  do i = 1, dim_obs_p ! loop over all observations
+    avelai = 0.0 ! use avelai as grid cell average collecter
+    do j = 1, clm_endp-clm_begp+1 ! loop over all patches
+      write(*,*) 'DEBUG LAI : obs_index_p(i), clm_patch2gc(j)', obs_index_p(i), clm_patch2gc(j)
+      if (obs_index_p(i)==clm_patch2gc(j)) then
+        write(*,*) 'DEBUG LAI : state dsladlai', state_p(j+2*clm_varsize)
+        if (state_p(j+2*clm_varsize)>0.0) then
+          avelai = avelai + clm_patchwt(j) * ((state_p(j+1*clm_varsize)&
+            *(exp(state_p(j)*state_p(j+2*clm_varsize)) - 1.0))/state_p(j+2*clm_varsize)) ! formula for tlai from leafc
+        else
+          avelai = avelai +  clm_patchwt(j) *(state_p(j+1*clm_varsize)*state_p(j)) ! 2nd formula
+        endif ! dsladlai decider
+      endif ! this patch is in gridcell of observation
+      write(*,*) 'DEBUG LAI : average ', avelai
+    enddo ! all patches for obs checked
+    m_state_p(i) = avelai ! == lai of gridcell where the observation is
+  enddo ! all observations
+endif ! clmupdate_lai == 2 : m_state_p now contains lai for each gridcell with an observation.
+
+! LAI assimilation assuming statevec contains leafc, livestemc, deadstemc
+! for each vegetated patch (bare-ground PFT excluded).
+! H uses slatop from state (update_lai_params=1) or pftcon; dsladlai from pftcon.
+if (clmupdate_lai==3) then
+  lpointobs = .false.
+
+  do i = 1, dim_obs_p
+    avelai = 0.0
+    do j = 1, clm_varsize
+      if (obs_index_p(i)==clm_patch2gc(j)) then
+        if (clmupdate_lai_params == 1) then
+          slatop_val = state_p(j + lai_mode3_npool*clm_varsize)
+        else
+          slatop_val = pftcon%slatop(clm_lai_patch_itype(j))
+        end if
+        if (pftcon%dsladlai(clm_lai_patch_itype(j)) > 0.0) then
+          avelai = avelai + clm_patchwt(j) * ((slatop_val&
+            *(exp(state_p(j)*pftcon%dsladlai(clm_lai_patch_itype(j))) - 1.0))&
+            /pftcon%dsladlai(clm_lai_patch_itype(j)))
+        else
+          avelai = avelai + clm_patchwt(j) * (slatop_val*state_p(j))
+        endif
+      endif
+    enddo
+    m_state_p(i) = avelai
+  enddo
+
+#ifdef PDAF_DEBUG
+  IF(clmt_printensemble == step .OR. clmt_printensemble == -1) THEN
+    WRITE(fn_hx, "(a,i5.5,a,i5.5,a)") "lai_hx_", mype_world, "_", step, ".txt"
+    OPEN(unit=72, file=fn_hx, action="write", status="replace")
+    WRITE(72,"(a,i8)") "# mode3 H(x) at PDAF step ", step
+    WRITE(72,"(a,i8)") "# dim_obs_p ", dim_obs_p
+    DO i = 1, dim_obs_p
+      WRITE(72,"(a,i8,a)") "# obs_index ", obs_index_p(i)
+      WRITE(72,"(es22.15)") m_state_p(i)
+    END DO
+    CLOSE(72)
+  END IF
+#endif
+endif ! clmupdate_lai == 3
+
+
+
 #if defined CLMSA
 if (clmupdate_T==1) then
 
diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h
index a734c335..0b01df7f 100755
--- a/interface/model/common/enkf.h
+++ b/interface/model/common/enkf.h
@@ -91,6 +91,9 @@ GLOBAL int clmupdate_swc;
 GLOBAL int clmupdate_T;
 GLOBAL int clmupdate_texture;
 GLOBAL int clmupdate_tws;
+GLOBAL int clmupdate_lai;
+GLOBAL int clmupdate_lai_params;
+GLOBAL double clmupdate_lai_incr_w;
 GLOBAL int clmprint_swc;
 GLOBAL int clmprint_et;
 GLOBAL int clmprint_inc;
diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c
index 27b26c46..4b006c50 100755
--- a/interface/model/common/read_enkfpar.c
+++ b/interface/model/common/read_enkfpar.c
@@ -79,6 +79,9 @@ void read_enkfpar(char *parname)
   clmupdate_swc         = iniparser_getint(pardict,"CLM:update_swc",1);
   clmupdate_T           = iniparser_getint(pardict,"CLM:update_T",0);
   clmupdate_texture     = iniparser_getint(pardict,"CLM:update_texture",0);
+  clmupdate_lai         = iniparser_getint(pardict,"CLM:update_lai",0);
+  clmupdate_lai_params  = iniparser_getint(pardict,"CLM:update_lai_params",0);
+  clmupdate_lai_incr_w  = iniparser_getdouble(pardict,"CLM:update_lai_incr_w",1);
   clmprint_swc          = iniparser_getint(pardict,"CLM:print_swc",0);
   clmprint_et           = iniparser_getint(pardict,"CLM:print_et",0);
   clmprint_inc           = iniparser_getint(pardict,"CLM:print_inc",0);
diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 68bd8c54..3809b60f 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -52,10 +52,17 @@ module enkf_clm_mod
   ! LST in LST assimilation (clmupdate_T)
   real(r8),allocatable :: clm_paramarr(:)  !hcp CLM parameter vector (f.e. LAI)
   integer, allocatable :: state_clm2pdaf_p(:,:) !Index of column in hydraulic active state vector (nlevsoi,endc-begc+1)
+  real(r8),allocatable :: clm_patch2gc(:) ! index array to simplify patch to gridcell averaging
+  real(r8),allocatable :: clm_patchwt(:) ! weight array to simplify patch to gridcell averaging
+  integer,allocatable :: clm_lai_patch_idx(:)   ! mode-3: CLM patch index per state slot
+  integer,allocatable :: clm_lai_patch_itype(:) ! mode-3: PFT type per state slot
   integer(c_int),bind(C,name="clmupdate_swc")     :: clmupdate_swc
   integer(c_int),bind(C,name="clmupdate_T")     :: clmupdate_T  ! by hcp
   integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture
   integer(c_int),bind(C,name="clmprint_swc")      :: clmprint_swc
+  integer(c_int),bind(C,name="clmupdate_lai")    :: clmupdate_lai
+  integer(c_int),bind(C,name="clmupdate_lai_params") :: clmupdate_lai_params
+  real(c_double),bind(C,name="clmupdate_lai_incr_w") :: clmupdate_lai_incr_w
 
   ! Yorck
   integer(c_int),bind(C,name="clmupdate_tws") :: clmupdate_tws
@@ -76,6 +83,11 @@ module enkf_clm_mod
   logical :: first_cycle = .TRUE.
   logical :: use_omi_model = .FALSE.
 
+  ! Mode-3 state layout: 3 carbon pools; +2 PFT params when update_lai_params=1
+  ! (kmax is private/protected in PhotosynthesisMod and cannot be updated from PDAF)
+  integer, parameter :: lai_mode3_npool = 3
+  integer, parameter :: lai_mode3_nparam = 2   ! slatop, medlynslope
+
 
   ! OMI --> I want to update the observation type after each observation comes in.
   ! problem: the observation type is updated before the update of the assimilation
@@ -125,12 +137,19 @@ module enkf_clm_mod
                                 ! (currently not used for eclm)
   logical :: newgridcell        !only eclm
 
+  ! Arrays for LAI Assimilation
+  real(r8),allocatable :: tlai(:)
+  real(r8),allocatable :: tlai_orig(:)
+  real(r8),allocatable :: lai_wtsum2_gc(:)
+
   contains
 
 #if defined CLMSA
   subroutine define_clm_statevec(mype)
     use decompMod , only : get_proc_bounds
     use clm_varpar   , only : nlevsoi
+    use PatchType  , only : patch
+    use pftconMod  , only : noveg
     use clm_varcon, only: set_averaging_to_zero
     use PDAF_interfaces_module, only: PDAF_reset_dim_p
 
@@ -142,6 +161,7 @@ subroutine define_clm_statevec(mype)
     integer :: begc, endc   ! per-proc beginning and ending column indices
     integer :: begl, endl   ! per-proc beginning and ending landunit indices
     integer :: begg, endg   ! per-proc gridcell ending gridcell indices
+    integer :: i, cc
 
 
     call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp)
@@ -191,6 +211,83 @@ subroutine define_clm_statevec(mype)
       clm_statevecsize = clm_statevecsize + 3*((endg-begg+1)*nlevsoi)
     end if
 
+    ! LAI assimilation
+    if(clmupdate_lai==1) then
+        ! Allocate array to store total leaf area index per patch
+        if (allocated(tlai)) deallocate(tlai)
+        allocate(tlai(clm_begp:clm_endp))
+        if (allocated(tlai_orig)) deallocate(tlai_orig)
+        allocate(tlai_orig(clm_begp:clm_endp))
+        if (allocated(lai_wtsum2_gc)) deallocate(lai_wtsum2_gc)
+        allocate(lai_wtsum2_gc(clm_begg:clm_endg))
+
+        ! #patches values per grid-cell ! clm_varsize not accurate due to variable #patches/grid cell
+        ! set sizes
+        clm_varsize = 1 ! each grid cell gets 1 value of LAI
+        ! Currently no combination of SWC and LAI DA
+
+        clm_statevecsize = (clm_endg - clm_begg + 1) ! statevector is the size of the num of gridcells
+        !(clm_endp-clm_begp+1) ! So like this if lai is set it takes priority
+
+        if (clmupdate_lai_params==1) then ! add parameters to the state vector.
+            ! clm_varsize = (clm_endg - clm_begg+1) ! Currently no combination of SWC and LAI DA
+            ! clm_statevecsize = (clm_endg - clm_begg+1)*1 ! 1 parameters
+            clm_varsize      =  (clm_endp-clm_begp+1)+1 ! Currently no combination of SWC and LAI DA
+            clm_statevecsize =  (clm_endp-clm_begp+1)+1 ! 1 parameters
+            ! So like this if lai is set it takes priority
+        end if
+
+        ! if (clmupdate_lai_params==3) then
+        !     clm_varsize      =  (clm_endp-clm_begp+1)
+        !     clm_statevecsize =  3*clm_varsize
+        ! end if
+
+    end if
+
+    ! for obs_op version allocate additional helper index array
+    if(clmupdate_lai==2) then
+      clm_varsize      =  (clm_endp-clm_begp+1) ! equals number of patches
+      clm_statevecsize =  3*clm_varsize ! 3 var/params per patch
+
+      if (allocated(clm_patch2gc)) deallocate(clm_patch2gc)
+      allocate(clm_patch2gc(clm_varsize))
+      if (allocated(clm_patchwt)) deallocate(clm_patchwt)
+      allocate(clm_patchwt(clm_varsize))
+
+    endif
+
+    ! LAI assimilation mode 3: leafc, livestemc, deadstemc per vegetated patch
+    if(clmupdate_lai==3) then
+      clm_varsize = 0
+      do i=clm_begp,clm_endp
+        if (patch%itype(i) /= noveg) clm_varsize = clm_varsize + 1
+      end do
+      clm_statevecsize = lai_mode3_npool*clm_varsize
+      if (clmupdate_lai_params == 1) then
+        clm_statevecsize = (lai_mode3_npool + lai_mode3_nparam)*clm_varsize
+      end if
+
+      if (allocated(clm_patch2gc)) deallocate(clm_patch2gc)
+      allocate(clm_patch2gc(clm_varsize))
+      if (allocated(clm_patchwt)) deallocate(clm_patchwt)
+      allocate(clm_patchwt(clm_varsize))
+      if (allocated(clm_lai_patch_idx)) deallocate(clm_lai_patch_idx)
+      allocate(clm_lai_patch_idx(clm_varsize))
+      if (allocated(clm_lai_patch_itype)) deallocate(clm_lai_patch_itype)
+      allocate(clm_lai_patch_itype(clm_varsize))
+
+      cc = 1
+      do i=clm_begp,clm_endp
+        if (patch%itype(i) /= noveg) then
+          clm_lai_patch_idx(cc) = i
+          clm_lai_patch_itype(cc) = patch%itype(i)
+          clm_patch2gc(cc) = real(patch%gridcell(i), r8)
+          clm_patchwt(cc) = patch%wtgcell(i)
+          cc = cc + 1
+        end if
+      end do
+    endif
+
     !hcp LST DA
     if(clmupdate_T==1) then
       error stop "Not implemented: clmupdate_T.eq.1"
@@ -240,17 +337,18 @@ subroutine define_clm_statevec(mype)
 #endif
 
     IF (allocated(clm_statevec)) deallocate(clm_statevec)
-    if ((clmupdate_swc/=0) .or. (clmupdate_T/=0) .or. (clmupdate_texture/=0) .or. (clmupdate_tws/=0)) then
+    if ((clmupdate_swc/=0) .or. (clmupdate_T/=0) .or. (clmupdate_texture/=0) .or. (clmupdate_tws/=0) .or. (clmupdate_lai/=0)) then
       !hcp added condition
       allocate(clm_statevec(clm_statevecsize))
     end if
 
 
-    ! Allocate statevector-duplicate for saving original column mean
-    ! values used in computing increments during updating the state
-    ! vector in column-mean-mode.
+    ! Allocate statevector-duplicate for saving the forecast state
+    ! vector. The forecast values may then be used during the
+    ! update-phase, when mapping analysis increments back to CLM (used
+    ! in SWC DA column-mean mode or LAI DA mode).
     IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig)
-    if ( (clmupdate_swc/=0 .and. clmstatevec_colmean/=0) .or. clmupdate_tws/=0 ) then
+    if ( (clmupdate_swc/=0 .and. clmstatevec_colmean/=0) .or. clmupdate_tws/=0  .or. (clmupdate_lai/=0)) then
       allocate(clm_statevec_orig(clm_statevecsize))
     end if
 
@@ -701,16 +799,27 @@ subroutine cleanup_clm_statevec()
     IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p)
     IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p)
     IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig)
+    IF (allocated(clm_patch2gc)) deallocate(clm_patch2gc)
+    IF (allocated(clm_patchwt)) deallocate(clm_patchwt)
+    IF (allocated(clm_lai_patch_idx)) deallocate(clm_lai_patch_idx)
+    IF (allocated(clm_lai_patch_itype)) deallocate(clm_lai_patch_itype)
+    IF (allocated(lai_wtsum2_gc)) deallocate(lai_wtsum2_gc)
+    IF (allocated(tlai)) deallocate(tlai)
+    IF (allocated(tlai_orig)) deallocate(tlai_orig)
 
   end subroutine cleanup_clm_statevec
 
   subroutine set_clm_statevec(tstartcycle, mype)
     use clm_instMod, only : soilstate_inst, waterstate_inst
+    use clm_instMod, only : bgc_vegetation_inst, canopystate_inst
     use clm_varpar   , only : nlevsoi
     ! use clm_varcon, only: nameg, namec
     ! use GetGlobalValuesMod, only: GetGlobalWrite
     use ColumnType , only : col
+    use PatchType  , only : patch
+    use pftconMod  , only : pftcon
     use shr_kind_mod, only: r8 => shr_kind_r8
+
     implicit none
     integer,intent(in) :: tstartcycle
     integer,intent(in) :: mype
@@ -718,10 +827,14 @@ subroutine set_clm_statevec(tstartcycle, mype)
     real(r8), pointer :: psand(:,:)
     real(r8), pointer :: pclay(:,:)
     real(r8), pointer :: porgm(:,:)
-    integer :: i,j,jj,g,c,cc,offset
+    real(r8), pointer :: leafc(:)
+    real(r8), pointer :: livestemc(:)
+    real(r8), pointer :: deadstemc(:)
+    integer :: i, j, jj, g, c, cc, offset, igc, itype
     integer :: n_c
     character (len = 34) :: fn    !TSMP-PDAF: function name for state vector output
     character (len = 34) :: fn2    !TSMP-PDAF: function name for swc output
+    character (len = 48) :: fn3    !TSMP-PDAF: function name for lai_leafc_* integrate output
 
     cc = 0
     offset = 0
@@ -730,6 +843,8 @@ subroutine set_clm_statevec(tstartcycle, mype)
     psand => soilstate_inst%cellsand_col
     pclay => soilstate_inst%cellclay_col
     porgm => soilstate_inst%cellorg_col
+    leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch
+
 
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN
@@ -754,6 +869,136 @@ subroutine set_clm_statevec(tstartcycle, mype)
       error stop "Not implemented clmupdate_swc.eq.2"
     end if
 
+    !> @author  Lukas Strebel, Haojin Zhao
+    !> @date    04.2026
+    ! LAI assimilation
+    ! Case 1: LAI all patches
+    ! Case 1: Transformations in Set and Update functions
+    if(clmupdate_lai==1) then
+      clm_statevec(:) = 0._r8 ! reset statevector to 0 because we add to it for the average
+      lai_wtsum2_gc(:) = 0._r8
+      do i=clm_begp,clm_endp
+        lai_wtsum2_gc(patch%gridcell(i)) = lai_wtsum2_gc(patch%gridcell(i)) &
+          + patch%wtgcell(i)**2
+      end do
+       do i=clm_begp,clm_endp ! iterate through patches
+           ! update the leaf area index based on leafC and SLA
+           ! Eq 3 from Thornton and Zimmerman, 2007, J Clim, 20, 3902-3923.
+           ! (slatop(ivt(p))*(exp(leafc(p)*dsladlai(ivt(p))) - 1._r8))/dsladlai(ivt(p))
+           if (pftcon%dsladlai(patch%itype(i)) > 0._r8) then
+             tlai(i) = (pftcon%slatop(patch%itype(i))*(exp(leafc(i)*pftcon%dsladlai(patch%itype(i))) - 1._r8)) &
+               /pftcon%dsladlai(patch%itype(i))
+           else
+             tlai(i) = pftcon%slatop(patch%itype(i)) * leafc(i)
+           end if
+           tlai(i) = max(0._r8, tlai(i)) ! don't allow negative LAI
+           tlai_orig(i) = tlai(i)
+           ! Map CLM gridcell index (begg:endg) to PE-local state index (1:clm_statevecsize)
+           igc = patch%gridcell(i) - clm_begg + 1
+           if (igc >= 1 .and. igc <= clm_statevecsize) then
+             clm_statevec(igc) = clm_statevec(igc) + patch%wtgcell(i)*tlai(i)
+             write(*,'(A,I6,A,I6,A,F10.4,A,F10.4,A,F10.4,A,F10.4,A,F10.4,A,F10.4,A,F10.4,A,F10.4)') &
+                 "LAI1 DEBUG setvec ", i, ";", igc, ";", patch%wtgcell(i), ";", &
+                 tlai(i), ";", patch%wtgcell(i)*tlai(i), ";", clm_statevec(igc), ";", &
+                 canopystate_inst%tlai_patch(i), ";", leafc(i), ";", &
+                 pftcon%dsladlai(patch%itype(i)), ";", pftcon%slatop(patch%itype(i))
+           else
+             write(*,*) 'DEBUG LAI : igc out of bounds in LAI mapping: ', igc, &
+               ' valid range: 1-', clm_statevecsize, ' i=', i
+           end if
+       end do
+       if (allocated(clm_statevec_orig)) then
+         do igc = 1, min(clm_endg - clm_begg + 1, clm_statevecsize)
+           clm_statevec_orig(igc) = clm_statevec(igc)
+         end do
+       end if
+
+#ifdef PDAF_DEBUG
+      IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN
+        IF(clmupdate_lai/=0) THEN
+          ! Same pattern as swcstate_*.integrate.*: one file, consecutive es22.15 blocks.
+          ! Order: leafc(patch) tlai_calc(patch) tlai_clm(patch) L_gc(gridcell) lai_wtsum2_gc(gridcell)
+          WRITE(fn3, "(a,i5.5,a,i5.5,a)") "lai_leafc_", mype, ".integrate.", tstartcycle + 1, ".txt"
+          OPEN(unit=71, file=fn3, action="write", status="replace")
+          WRITE(71,"(a)") "# blocks: leafc(patch) tlai_calc(patch) tlai_clm(patch) L_gc(gridcell) lai_wtsum2_sq(gridcell)"
+          WRITE(71,"(es22.15)") leafc(clm_begp:clm_endp)
+          WRITE(71,"(es22.15)") tlai(clm_begp:clm_endp)
+          WRITE(71,"(es22.15)") canopystate_inst%tlai_patch(clm_begp:clm_endp)
+          WRITE(71,"(es22.15)") clm_statevec(1:clm_statevecsize)
+          WRITE(71,"(es22.15)") lai_wtsum2_gc(clm_begg:clm_endg)
+          CLOSE(71)
+        END IF
+      END IF
+#endif
+
+    endif
+
+    ! Case 2: statevec with leafc, slatop, dsladlai for transform in obs_op
+    if (clmupdate_lai==2) then
+      cc = 1
+      do i=clm_begp,clm_endp
+        clm_statevec(cc) = leafc(i)
+        clm_statevec(cc + 1*clm_varsize) = pftcon%slatop(patch%itype(i))
+        clm_statevec(cc + 2*clm_varsize) = pftcon%dsladlai(patch%itype(i))
+        clm_patch2gc(cc) = patch%gridcell(i)
+        clm_patchwt(cc) = patch%wtgcell(i)
+        cc = cc + 1
+      enddo
+      write(*,*) 'DEBUG LAI : statevec ', clm_statevec(:)
+      write(*,*) 'DEBUG LAI : patches ', clm_patch2gc(:), clm_patchwt(:)
+    endif
+
+    ! Case 3: leafc, livestemc, deadstemc per vegetated patch; H(x) in obs_op
+    if (clmupdate_lai==3) then
+      livestemc => bgc_vegetation_inst%cnveg_carbonstate_inst%livestemc_patch
+      deadstemc => bgc_vegetation_inst%cnveg_carbonstate_inst%deadstemc_patch
+      cc = 1
+      do j = 1, clm_varsize
+        i = clm_lai_patch_idx(j)
+        clm_statevec(cc) = leafc(i)
+        clm_statevec(cc + clm_varsize) = livestemc(i)
+        clm_statevec(cc + 2*clm_varsize) = deadstemc(i)
+        cc = cc + 1
+      end do
+      if (allocated(clm_statevec_orig)) then
+        clm_statevec_orig(:) = clm_statevec(:)
+      end if
+
+      if (clmupdate_lai_params == 1) then
+        cc = 1
+        do j = 1, clm_varsize
+          itype = clm_lai_patch_itype(j)
+          clm_statevec(cc + lai_mode3_npool    *clm_varsize) = pftcon%slatop(itype)
+          clm_statevec(cc + (lai_mode3_npool+1)*clm_varsize) = pftcon%medlynslope(itype)
+          cc = cc + 1
+        end do
+      end if
+
+#ifdef PDAF_DEBUG
+      IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN
+        call write_lai_da_mode3_debug('integrate', tstartcycle + 1, mype)
+      END IF
+#endif
+    endif
+
+    if (clmupdate_lai_params==1 .and. clmupdate_lai==1) then
+      cc = 1
+      do i=clm_begp,clm_endp
+          clm_statevec(cc+1*clm_varsize+offset) = pftcon%slatop(patch%itype(i))
+          cc = cc + 1
+      end do
+    endif
+
+    ! if (clmupdate_lai_params==3) then
+    !   cc = 1
+    !   do i=clm_begp,clm_endp
+    !     clm_statevec(cc+1*clm_varsize+offset) = pftcon%slatop(patch%itype(i))
+    !     clm_statevec(cc+2*clm_varsize+offset) = params_inst%kmax(patch%itype(i),1)
+    !     clm_statevec(cc+3*clm_varsize+offset) = params_inst%kmax(patch%itype(i),2)
+    !     cc = cc + 1
+    !   end do
+    ! endif
+
     !hcp  LAI
     if(clmupdate_T==1) then
       error stop "Not implemented: clmupdate_T.eq.1"
@@ -1144,10 +1389,125 @@ end subroutine set_clm_statevec_tws
 
 
 
+  subroutine lai_patch_tlai_blend(tlai_fc, lai_gc_fc, lai_gc_an, tlai_out)
+    use shr_kind_mod , only : r8 => shr_kind_r8
+
+    implicit none
+
+    real(r8), intent(in) :: tlai_fc
+    real(r8), intent(in) :: lai_gc_fc
+    real(r8), intent(in) :: lai_gc_an
+    real(r8), intent(out) :: tlai_out
+
+    real(r8) :: incr_lai_gc
+    real(r8) :: tlai_add
+    real(r8) :: tlai_mult
+    real(r8) :: w_add
+
+    incr_lai_gc = lai_gc_an - lai_gc_fc
+    tlai_add = tlai_fc + incr_lai_gc
+    if (abs(lai_gc_fc) > 1.0e-12_r8) then
+      tlai_mult = tlai_fc * (lai_gc_an / lai_gc_fc)
+    else
+      tlai_mult = tlai_add
+    end if
+    w_add = real(clmupdate_lai_incr_w, r8)
+    tlai_out = w_add * tlai_add + (1._r8 - w_add) * tlai_mult
+    tlai_out = max(0._r8, tlai_out)
+
+  end subroutine lai_patch_tlai_blend
+
+
+  subroutine write_lai_da_mode3_debug(phase, cycle, mype)
+    use clm_instMod, only : bgc_vegetation_inst, canopystate_inst
+
+    implicit none
+
+    character(len=*), intent(in) :: phase
+    integer, intent(in) :: cycle
+    integer, intent(in) :: mype
+
+    real(r8), pointer :: leafc(:)
+    real(r8), pointer :: livestemc(:)
+    real(r8), pointer :: deadstemc(:)
+    real(r8), pointer :: leafn(:)
+    real(r8), pointer :: livestemn(:)
+    real(r8), pointer :: deadstemn(:)
+
+    character(len=48) :: fn
+    integer :: j, i
+
+    IF(clmt_printensemble /= cycle .AND. clmt_printensemble /= -1) RETURN
+
+    leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch
+    livestemc => bgc_vegetation_inst%cnveg_carbonstate_inst%livestemc_patch
+    deadstemc => bgc_vegetation_inst%cnveg_carbonstate_inst%deadstemc_patch
+    leafn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%leafn_patch
+    livestemn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%livestemn_patch
+    deadstemn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%deadstemn_patch
+
+    WRITE(fn, "(a,i5.5,a,a,a,i5.5,a)") "lai_da_", mype, ".", trim(phase), ".", cycle, ".txt"
+    OPEN(unit=71, file=fn, action="write", status="replace")
+    WRITE(71,"(a)") "# mode3 blocks: clm_statevec clm_lai_patch_idx clm_lai_patch_itype clm_patch2gc clm_patchwt"
+    IF (clmupdate_lai_params == 1) THEN
+      WRITE(71,"(a)") "# param blocks (if update_lai_params=1): slatop medlynslope per slot"
+    END IF
+    WRITE(71,"(a)") "#            leafc livestemc deadstemc leafn livestemn deadstemn tlai_clm"
+    DO i = 1, clm_statevecsize
+      WRITE(71,"(es22.15)") clm_statevec(i)
+    END DO
+    DO j = 1, clm_varsize
+      WRITE(71,"(i12)") clm_lai_patch_idx(j)
+    END DO
+    DO j = 1, clm_varsize
+      WRITE(71,"(i12)") clm_lai_patch_itype(j)
+    END DO
+    DO j = 1, clm_varsize
+      WRITE(71,"(es22.15)") clm_patch2gc(j)
+    END DO
+    DO j = 1, clm_varsize
+      WRITE(71,"(es22.15)") clm_patchwt(j)
+    END DO
+    DO j = 1, clm_varsize
+      i = clm_lai_patch_idx(j)
+      WRITE(71,"(es22.15)") leafc(i)
+    END DO
+    DO j = 1, clm_varsize
+      i = clm_lai_patch_idx(j)
+      WRITE(71,"(es22.15)") livestemc(i)
+    END DO
+    DO j = 1, clm_varsize
+      i = clm_lai_patch_idx(j)
+      WRITE(71,"(es22.15)") deadstemc(i)
+    END DO
+    DO j = 1, clm_varsize
+      i = clm_lai_patch_idx(j)
+      WRITE(71,"(es22.15)") leafn(i)
+    END DO
+    DO j = 1, clm_varsize
+      i = clm_lai_patch_idx(j)
+      WRITE(71,"(es22.15)") livestemn(i)
+    END DO
+    DO j = 1, clm_varsize
+      i = clm_lai_patch_idx(j)
+      WRITE(71,"(es22.15)") deadstemn(i)
+    END DO
+    DO j = 1, clm_varsize
+      i = clm_lai_patch_idx(j)
+      WRITE(71,"(es22.15)") canopystate_inst%tlai_patch(i)
+    END DO
+    CLOSE(71)
+
+  end subroutine write_lai_da_mode3_debug
+
+
   subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     use clm_time_manager  , only : update_DA_nstep
     use shr_kind_mod , only : r8 => shr_kind_r8
     use clm_instMod, only : waterstate_inst
+    use clm_instMod, only : bgc_vegetation_inst, canopystate_inst
+    use pftconMod       , only : pftcon
+    use PatchType       , only : patch
 
     implicit none
 
@@ -1156,16 +1516,35 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
 
     real(r8), pointer :: h2osoi_liq(:,:)  ! liquid water (kg/m2)
     real(r8), pointer :: h2osoi_ice(:,:)
+    real(r8), pointer :: leafc(:)         ! leaf carbon of patch
+    real(r8), pointer :: leafn(:)         ! leaf nitrogen of patch
+    real(r8), pointer :: livestemc(:)
+    real(r8), pointer :: deadstemc(:)
+    real(r8), pointer :: livestemn(:)
+    real(r8), pointer :: deadstemn(:)
+
+    real(r8) :: lai_gc_an
+    real(r8) :: lai_gc_fc
 
     integer :: i
+    integer :: j
+    integer :: itype
+    integer :: cc
+    integer :: offset
+    integer :: igc
     character (len = 31) :: fn    !TSMP-PDAF: function name for state vector output
-    character (len = 32) :: fn5    !TSMP-PDAF: function name for state vector outpu
-    character (len = 32) :: fn6    !TSMP-PDAF: function name for state vector outpu
+    character (len = 31) :: fn2    !TSMP-PDAF: function name for state vector output
+    character (len = 32) :: fn5    !TSMP-PDAF: function name for state vector output
+    character (len = 32) :: fn6    !TSMP-PDAF: function name for state vector output
+    character (len = 32) :: fn7    !TSMP-PDAF: function name for state vector output
 
     logical :: swc_zero_before_update
 
     swc_zero_before_update = .false.
 
+    cc = 1
+    offset = 0
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN
       ! TSMP-PDAF: For debug runs, output the state vector in files
@@ -1181,6 +1560,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     h2osoi_liq    => waterstate_inst%h2osoi_liq_col
     h2osoi_ice    => waterstate_inst%h2osoi_ice_col
 
+    leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch
+    leafn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%leafn_patch
+    livestemc => bgc_vegetation_inst%cnveg_carbonstate_inst%livestemc_patch
+    deadstemc => bgc_vegetation_inst%cnveg_carbonstate_inst%deadstemc_patch
+    livestemn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%livestemn_patch
+    deadstemn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%deadstemn_patch
+
+
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN
 
@@ -1198,6 +1586,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
         CLOSE(71)
       END IF
 
+      IF(clmupdate_lai/=0) THEN
+        ! TSMP-PDAF: For debug runs, output the state vector in files
+        WRITE(fn5, "(a,i5.5,a,i5.5,a)") "leafc", mype, ".bef_up.", tstartcycle, ".txt"
+        OPEN(unit=71, file=fn5, action="write")
+        WRITE (71,"(es22.15)") leafc(:)
+        CLOSE(71)
+      END IF
+
     END IF
 #endif
 
@@ -1233,6 +1629,169 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
 
 
 
+    ! LAI assimilation:
+    ! Case 1: gridcell-mean LAI state; map analysis back to patches via
+    ! additive/multiplicative blend (CLM:update_lai_incr_w).
+    if(clmupdate_lai==1) then
+      do i=clm_begp,clm_endp
+        igc = patch%gridcell(i) - clm_begg + 1
+        if (igc >= 1 .and. igc <= clm_statevecsize .and. allocated(clm_statevec_orig)) then
+          lai_gc_fc = clm_statevec_orig(igc)
+          lai_gc_an = clm_statevec(igc)
+          call lai_patch_tlai_blend(tlai_orig(i), lai_gc_fc, lai_gc_an, tlai(i))
+        else
+          tlai(i) = 0._r8
+        endif
+        ! update leafc based on Eq 3 from Thornton and Zimmerman, 2007, J Clim, 20, 3902-3923
+        ! reformulate the equation to solve for leafc
+        ! => leafc(p) = log(((tlai(p) * dsladlai(ivt(p)))/slatop(ivt(p))) + 1._r8) / dsladlai(ivt(p))
+        if (tlai(i) > 0._r8) then ! invalid log protection
+          if (pftcon%dsladlai(patch%itype(i)) > 0._r8) then
+            if (pftcon%slatop(patch%itype(i)) > 0._r8) then
+              leafc(i) = log(((tlai(i) * pftcon%dsladlai(patch%itype(i))) &
+                / pftcon%slatop(patch%itype(i))) + 1.0_r8) &
+                / pftcon%dsladlai(patch%itype(i))
+            else
+              leafc(i) = 0._r8
+            endif
+          else
+              ! Forward: tlai = slatop * leafc  (dsladlai == 0) => leafc = tlai / slatop
+              if (pftcon%slatop(patch%itype(i)) > 0._r8) then
+                leafc(i) = tlai(i) / pftcon%slatop(patch%itype(i))
+              else
+                leafc(i) = 0._r8
+              end if
+          endif
+        else
+          leafc(i) = 0._r8
+        endif
+        leafc(i) = max(0._r8, leafc(i)) ! Don't allow negative leaf carbon
+        if (pftcon%leafcn(patch%itype(i)) > 0._r8) then
+          leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i))
+        else
+          leafn(i) = 0._r8
+        endif
+        ! DEBUG ONLY
+        ! print *, "LAI1 DEBUG leafc,tlai(wt) AFTER: ", i, ";", leafc(i), ";", tlai(i)
+      end do
+
+#ifdef PDAF_DEBUG
+        IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN
+          ! TSMP-PDAF: lai & leafc  (legacy column vector)
+          WRITE(fn7, "(a,i5.5,a,i5.5,a)") "lai_leafc_", mype, ".update.", tstartcycle, ".txt"
+          OPEN(unit=71, file=fn7, action="write", status="replace")
+          WRITE(71,"(a)") "# patch  leafc_aft  tlai_patch_aft"
+          DO i = clm_begp, clm_endp
+            WRITE(71,"(i8,2es22.15)") i, leafc(i), tlai(i)
+          END DO
+          CLOSE(71)
+        END IF
+#endif
+    endif
+
+    ! Case 2: statevec contains leafc, slatop, dsladlai
+    ! update of diagnostic tlai left to be done by clm5 itself.
+    if(clmupdate_lai==2) then
+      cc = 1
+      do i=clm_begp,clm_endp
+        leafc(i) = clm_statevec(cc)
+        leafc(i) = max(0._r8, leafc(i))
+        if (pftcon%leafcn(patch%itype(i)) > 0._r8) then
+          leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i))
+        else
+          leafn(i) = 0._r8
+        end if
+        if(clmupdate_lai_params==2) then
+          pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize)
+          pftcon%dsladlai(patch%itype(i)) = clm_statevec(cc+2*clm_varsize)
+        endif
+
+        cc = cc + 1
+      enddo
+
+#ifdef PDAF_DEBUG
+        IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN
+          ! TSMP-PDAF: For debug runs, output the state vector in files
+          WRITE(fn7, "(a,i5.5,a,i5.5,a)") "leafc_", mype, ".update.", tstartcycle, ".txt"
+          OPEN(unit=71, file=fn7, action="write")
+          WRITE (71,"(es22.15)") leafc(:)
+          CLOSE(71)
+        END IF
+#endif
+
+    endif
+
+    ! Case 3: leafc, livestemc, deadstemc per vegetated patch
+    if(clmupdate_lai==3) then
+      cc = 1
+      do j = 1, clm_varsize
+        i = clm_lai_patch_idx(j)
+        itype = clm_lai_patch_itype(j)
+        leafc(i) = max(0._r8, clm_statevec(cc))
+        livestemc(i) = max(0._r8, clm_statevec(cc + clm_varsize))
+        deadstemc(i) = max(0._r8, clm_statevec(cc + 2*clm_varsize))
+        if (pftcon%leafcn(itype) > 0._r8) then
+          leafn(i) = leafc(i) / pftcon%leafcn(itype)
+        else
+          leafn(i) = 0._r8
+        end if
+        if (pftcon%livewdcn(itype) > 0._r8) then
+          livestemn(i) = livestemc(i) / pftcon%livewdcn(itype)
+        else
+          livestemn(i) = 0._r8
+        end if
+        if (pftcon%deadwdcn(itype) > 0._r8) then
+          deadstemn(i) = deadstemc(i) / pftcon%deadwdcn(itype)
+        else
+          deadstemn(i) = 0._r8
+        end if
+        cc = cc + 1
+      end do
+
+      if (clmupdate_lai_params == 1) then
+        cc = 1
+        do j = 1, clm_varsize
+          itype = clm_lai_patch_itype(j)
+          pftcon%slatop(itype) = clm_statevec(cc + lai_mode3_npool    *clm_varsize)
+          pftcon%medlynslope(itype) = clm_statevec(cc + (lai_mode3_npool+1)*clm_varsize)
+          cc = cc + 1
+        end do
+      end if
+
+#ifdef PDAF_DEBUG
+      IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN
+        call write_lai_da_mode3_debug('update', tstartcycle, mype)
+      END IF
+#endif
+    endif
+
+    if (clmupdate_lai_params==1 .and. clmupdate_lai==1) then
+      cc = 1
+      do i=clm_begp,clm_endp
+        pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize+offset)
+        ! DEBUG ONLY
+        print *, "LST DEBUG PARAM,i: ", pftcon%slatop(patch%itype(i)), ",", i
+        cc = cc + 1
+      end do
+    endif
+
+    ! if (clmupdate_lai_params==3) then
+    !   cc = 1
+    !   do i=clm_begp,clm_endp
+    !     pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize+offset)
+    !     params_inst%kmax(patch%itype(i),1) = clm_statevec(cc+2*clm_varsize+offset)
+    !     params_inst%kmax(patch%itype(i),2) = clm_statevec(cc+3*clm_varsize+offset)
+
+    !     ! DEBUG ONLY
+    !     print *, "LST DEBUG PARAM slatop,i: ", pftcon%slatop(patch%itype(i)), ",", i
+    !     print *, "LST DEBUG PARAM kmax sun,i: ", params_inst%kmax(patch%itype(i),1), ",", i
+    !     print *, "LST DEBUG PARAM kmax shade,i: ", params_inst%kmax(patch%itype(i),2), ",", i
+
+    !     cc = cc + 1
+    !   end do
+
+    ! end if
+
   end subroutine update_clm
 
 
@@ -2102,6 +2661,8 @@ subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, &
     integer :: ncols, counter
     integer :: npatches, ncohorts
     real :: minlon, minlat, maxlon, maxlat
+    real(r8) :: dlon_span, dlat_span
+    real(r8), parameter :: tol_degen = 1.0e-6_r8 
     real(r8), pointer :: lon(:)
     real(r8), pointer :: lat(:)
     integer :: begg, endg   ! per-proc gridcell ending gridcell indices
@@ -2161,13 +2722,27 @@ subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, &
     minlat = MINVAL(lat(:) + 90)
     maxlat = MAXVAL(lat(:) + 90)
 
+    ! Degenerate horizontal domain (e.g. 1x1 grid): min==max => zero span;
+    ! Avoid dividing by (maxlon-minlon) or (maxlat-minlat).
+    dlon_span = maxval(lon(:) + 180._r8) - minval(lon(:) + 180._r8)
+    dlat_span = maxval(lat(:) + 90._r8) - minval(lat(:) + 90._r8)
+
     if(allocated(longxy_obs)) deallocate(longxy_obs)
     allocate(longxy_obs(dim_obs))
     if(allocated(latixy_obs)) deallocate(latixy_obs)
     allocate(latixy_obs(dim_obs))
 
     do i = 1, dim_obs
-       if(((lon_clmobs(i) + 180) - minlon) /= 0 .and. &
+       if (dlon_span <= tol_degen .and. dlat_span <= tol_degen) then
+          longxy_obs(i) = 1
+          latixy_obs(i) = 1
+       else if (dlon_span <= tol_degen) then
+          longxy_obs(i) = 1
+          latixy_obs(i) = ceiling(((lat_clmobs(i) + 90) - minlat) * nj / (maxlat - minlat))
+       else if (dlat_span <= tol_degen) then
+          longxy_obs(i) = ceiling(((lon_clmobs(i) + 180) - minlon) * ni / (maxlon - minlon))
+          latixy_obs(i) = 1
+       else if(((lon_clmobs(i) + 180) - minlon) /= 0 .and. &
          ((lat_clmobs(i) + 90) - minlat) /= 0) then
           longxy_obs(i) = ceiling(((lon_clmobs(i) + 180) - minlon) * ni / (maxlon - minlon)) !+ 1
           latixy_obs(i) = ceiling(((lat_clmobs(i) + 90) - minlat) * nj / (maxlat - minlat)) !+ 1
diff --git a/interface/model/eclm/print_update_clm_5.F90 b/interface/model/eclm/print_update_clm_5.F90
index 07ea27cd..227f7831 100644
--- a/interface/model/eclm/print_update_clm_5.F90
+++ b/interface/model/eclm/print_update_clm_5.F90
@@ -48,6 +48,13 @@ subroutine print_update_clm(ts,ttot) bind(C,name="print_update_clm")
     use netcdf, only : nf90_put_var
     use netcdf, only : nf90_close
     use enkf_clm_mod, only : clmupdate_swc,clmupdate_texture,clmprint_swc
+    use enkf_clm_mod, only : clmupdate_lai, clmupdate_lai_params
+    use enkf_clm_mod, only : clm_begg, clm_endg, clm_begp, clm_endp
+    use PatchType  , only : patch
+    use pftconMod  , only : pftcon
+    use mpi, only : mpi_gatherv, mpi_real8
+    use spmdmod    , only : npes, mpicom, iam
+    use decompmod  , only : get_proc_total
 
     implicit none
 
@@ -63,7 +70,7 @@ subroutine print_update_clm(ts,ttot) bind(C,name="print_update_clm")
     integer :: begc,endc      ! local beg/end columns
     integer :: begp,endp      ! local beg/end pfts
 
-    integer ::   isec, info, jn, jj, ji, g1, jx    ! temporary integer
+    integer ::   isec, info, jn, jj, ji, g1, jx, g, p    ! temporary integer
     real(r8), pointer :: swc(:,:)
     real(r8), pointer :: psand(:,:)
     real(r8), pointer :: pclay(:,:)
@@ -71,13 +78,26 @@ subroutine print_update_clm(ts,ttot) bind(C,name="print_update_clm")
     ! real(r8), pointer :: clmstate_tmp_local(:)
     ! real(r8), pointer :: clmstate_tmp_global(:)
     real(r8), allocatable :: clmstate_out(:,:,:)
+    real(r8), allocatable :: clmstate_out_2d(:,:)
+    real(r8), allocatable :: slatop_gc(:)
+    real(r8), allocatable :: medlynslope_gc(:)
+    real(r8), allocatable :: param_tmp_global(:)
     integer ,dimension(4) :: dimids
+    integer ,dimension(3) :: dimids_xyt
     integer ,dimension(1) :: il_var_id
-    integer :: il_file_id, ncvarid(4), status
+    integer :: il_file_id, ncvarid(4), slatop_varid, medlynslope_varid, status
     character(len = 300) :: update_filename
     integer :: ndlon,ndlat
-
-
+    logical :: write_lai_params
+    integer :: ier
+    integer :: beg
+    integer :: numrecvv(0:npes-1)
+    integer :: displsv(0:npes-1)
+    integer :: numsend
+    integer :: pid
+    integer :: ncells, nlunits, ncols, npfts_local, ncohorts
+
+    write_lai_params = (clmupdate_lai == 3 .and. clmupdate_lai_params == 1)
     call get_proc_global(ng=numg,nl=numl,nc=numc,np=nump)
     call get_proc_bounds(begg,endg,begl,endl,begc,endc,begp,endp)
     ! allocate(clmstate_tmp_local(nlevsoi*(-begc+endc)))
@@ -88,6 +108,7 @@ subroutine print_update_clm(ts,ttot) bind(C,name="print_update_clm")
     if (masterproc) then
       ! allocate(clmstate_tmp_global(nlevsoi*numg))
       allocate(clmstate_out(ndlon,ndlat,nlevsoi))
+      if (write_lai_params) allocate(clmstate_out_2d(ndlon,ndlat))
     end if
 
     if(masterproc) then
@@ -98,6 +119,7 @@ subroutine print_update_clm(ts,ttot) bind(C,name="print_update_clm")
         status =  nf90_def_dim(il_file_id, "y", ndlat, dimids(2))
         status =  nf90_def_dim(il_file_id, "z", nlevsoi, dimids(3))
         status =  nf90_def_dim(il_file_id, "t", ttot, dimids(4))
+        dimids_xyt = [ dimids(1), dimids(2), dimids(4) ]
 
         if(clmprint_swc==1) then
           status =  nf90_def_var(il_file_id, "swc", NF90_DOUBLE, dimids, ncvarid(1))
@@ -114,13 +136,17 @@ subroutine print_update_clm(ts,ttot) bind(C,name="print_update_clm")
           status =  nf90_def_var(il_file_id, "clay", NF90_DOUBLE, dimids, ncvarid(3))
           status =  nf90_def_var(il_file_id, "orgm", NF90_DOUBLE, dimids, ncvarid(4))
         end if
+
+        if (write_lai_params) then
+          status = nf90_def_var(il_file_id, "slatop", NF90_DOUBLE, dimids_xyt, slatop_varid)
+          status = nf90_def_var(il_file_id, "medlynslope", NF90_DOUBLE, dimids_xyt, medlynslope_varid)
+        end if
         status =  nf90_enddef(il_file_id)
       else
         status = nf90_open(update_filename,NF90_WRITE,il_file_id)
       end if
     end if
 
-
     if(clmprint_swc==1) then
       swc  => waterstate_inst%h2osoi_vol_col
       ! swc
@@ -209,9 +235,77 @@ subroutine print_update_clm(ts,ttot) bind(C,name="print_update_clm")
 
     end if
 
+    if (write_lai_params) then
+      allocate(slatop_gc(begg:endg))
+      allocate(medlynslope_gc(begg:endg))
+      slatop_gc(:) = 0._r8
+      medlynslope_gc(:) = 0._r8
+
+      do p = begp, endp
+        g = patch%gridcell(p)
+        if (g >= begg .and. g <= endg) then
+          slatop_gc(g) = slatop_gc(g) + patch%wtgcell(p) * pftcon%slatop(patch%itype(p))
+          medlynslope_gc(g) = medlynslope_gc(g) + patch%wtgcell(p) * pftcon%medlynslope(patch%itype(p))
+        end if
+      end do
+
+      call get_proc_total(iam, ncells, nlunits, ncols, npfts_local, ncohorts)
+      numsend = ncells
+      beg = begg
+      do pid = 0, npes-1
+        call get_proc_total(pid, ncells, nlunits, ncols, npfts_local, ncohorts)
+        numrecvv(pid) = ncells
+      end do
+      displsv(0) = 0
+      do pid = 1, npes-1
+        displsv(pid) = displsv(pid-1) + numrecvv(pid-1)
+      end do
+
+      if (masterproc) allocate(param_tmp_global(1:numg))
+
+      if (masterproc) then
+        call mpi_gatherv(slatop_gc(beg), numsend, MPI_REAL8, &
+             param_tmp_global, numrecvv, displsv, MPI_REAL8, 0, mpicom, ier)
+      else
+        call mpi_gatherv(slatop_gc(beg), numsend, MPI_REAL8, 0._r8, numrecvv, displsv, MPI_REAL8, 0, mpicom, ier)
+      end if
+      if (masterproc) then
+        do g1 = 1, numg
+          ji = mod(ldecomp%gdc2glo(g1)-1,ldomain%ni) + 1
+          jj = (ldecomp%gdc2glo(g1) - 1)/ldomain%ni + 1
+          clmstate_out_2d(ji,jj) = param_tmp_global(g1)
+        end do
+        status = nf90_inq_varid(il_file_id, "slatop", slatop_varid)
+        status = nf90_put_var(il_file_id, slatop_varid, clmstate_out_2d(:,:), &
+             start = [1, 1, ts], count = [ndlon, ndlat, 1])
+      end if
+
+      if (masterproc) then
+        call mpi_gatherv(medlynslope_gc(beg), numsend, MPI_REAL8, &
+             param_tmp_global, numrecvv, displsv, MPI_REAL8, 0, mpicom, ier)
+      else
+        call mpi_gatherv(medlynslope_gc(beg), numsend, MPI_REAL8, 0._r8, numrecvv, displsv, MPI_REAL8, 0, mpicom, ier)
+      end if
+      if (masterproc) then
+        do g1 = 1, numg
+          ji = mod(ldecomp%gdc2glo(g1)-1,ldomain%ni) + 1
+          jj = (ldecomp%gdc2glo(g1) - 1)/ldomain%ni + 1
+          clmstate_out_2d(ji,jj) = param_tmp_global(g1)
+        end do
+        status = nf90_inq_varid(il_file_id, "medlynslope", medlynslope_varid)
+        status = nf90_put_var(il_file_id, medlynslope_varid, clmstate_out_2d(:,:), &
+             start = [1, 1, ts], count = [ndlon, ndlat, 1])
+      end if
+
+      if (masterproc) deallocate(param_tmp_global)
+      deallocate(slatop_gc)
+      deallocate(medlynslope_gc)
+    end if
+
     if(masterproc) then
       status = nf90_close(il_file_id)
       deallocate(clmstate_out)
+      if (write_lai_params) deallocate(clmstate_out_2d)
       ! deallocate(clmstate_tmp_global)
     end if
 
diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c
index 87c2d093..a20ff35f 100644
--- a/interface/model/wrapper_tsmp.c
+++ b/interface/model/wrapper_tsmp.c
@@ -198,7 +198,7 @@ void integrate_tsmp() {
 void update_tsmp(){
 
 #if defined CLMSA
-  if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_tws != 0))){
+  if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_tws != 0) || (clmupdate_lai != 0))){
     update_clm(&tstartcycle, &mype_world);
 
 #ifndef CLMFIVE
@@ -208,6 +208,11 @@ void update_tsmp(){
 #endif
 
 #ifdef CLMFIVE
+    if((clmupdate_lai == 3 && clmupdate_lai_params == 1)){
+      /* Printing update could be implemented for SM-DA as well */
+      /* clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2 || */
+      print_update_clm(&tcycle, &total_steps);
+    }
     /* SWC: enkfpf.par input switch "CLM:print_inc" */
     if ((clmupdate_tws != 0) || (clmupdate_swc != 0 && clmprint_inc !=0)){
       print_inc_clm();