From 7071fe9e314a29370aeee7ede4ddc695074f4f73 Mon Sep 17 00:00:00 2001
From: Lukas Strebel <l.strebel@fz-juelich.de>
Date: Fri, 21 Nov 2025 10:05:51 +0100
Subject: [PATCH 01/25] LAI-DA: First implementation

---
 interface/model/common/enkf.h           |   2 +
 interface/model/common/read_enkfpar.c   |   2 +
 interface/model/eclm/enkf_clm_mod_5.F90 | 110 +++++++++++++++++++++++-
 3 files changed, 113 insertions(+), 1 deletion(-)

diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h
index f943ad90..1dcd2498 100755
--- a/interface/model/common/enkf.h
+++ b/interface/model/common/enkf.h
@@ -88,6 +88,8 @@ GLOBAL int nx_local,ny_local,nz_local;
 GLOBAL int clmupdate_swc;
 GLOBAL int clmupdate_T;
 GLOBAL int clmupdate_texture;
+GLOBAL int clmupdate_lai;
+GLOBAL int clmupdate_lai_params;
 GLOBAL int clmprint_swc;
 GLOBAL int clmprint_et;
 GLOBAL int clmstatevec_allcol;
diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c
index b2bfec28..2bf563d2 100755
--- a/interface/model/common/read_enkfpar.c
+++ b/interface/model/common/read_enkfpar.c
@@ -78,6 +78,8 @@ 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);
   clmprint_swc          = iniparser_getint(pardict,"CLM:print_swc",0);
   clmprint_et           = iniparser_getint(pardict,"CLM:print_et",0);
   clmstatevec_allcol    = iniparser_getint(pardict,"CLM:statevec_allcol",0);
diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 82fa382a..ee985acf 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -85,6 +85,10 @@ module enkf_clm_mod
                                 ! (currently not used for eclm)
   logical :: newgridcell        !only eclm
 
+  ! Arrays for LAI Assimilation
+  real(r8),allocatable :: tlai(:)
+  real(r8),allocatable :: old_frac(:)
+
   contains
 
 #if defined CLMSA
@@ -137,6 +141,23 @@ 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_endp-clm_begp+1))
+        ! Allocate array to store fraction of LAI the patch represents.
+        if (allocated(old_frac)) deallocate(old_frac)
+        allocate(old_frac(clm_endp-clm_begp+1))
+        ! set sizes
+        clm_varsize = 1 ! each grid cell gets 1 value of LAI
+        clm_statevecsize = (clm_endg - clm_begg + 1) ! statevector is the size of the num of gridcells
+        if (clmupdate_lai_params==1) then ! add parameters to the state vector.
+            clm_varsize = (clm_endg - clm_begg+1)
+            clm_statevecsize = (clm_endg - clm_begg+1)*1 ! 1 parameters
+        endif
+    end if
+
     !hcp LST DA
     if(clmupdate_T==1) then
       error stop "Not implemented: clmupdate_T.eq.1"
@@ -149,7 +170,7 @@ 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)) then
+    if ((clmupdate_swc/=0) .or. (clmupdate_T/=0) .or. (clmupdate_texture/=0) .or. (clmupdate_lai/=0)) then
       !hcp added condition
       allocate(clm_statevec(clm_statevecsize))
     end if
@@ -381,10 +402,13 @@ 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
     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
@@ -393,6 +417,7 @@ subroutine set_clm_statevec(tstartcycle, mype)
     real(r8), pointer :: psand(:,:)
     real(r8), pointer :: pclay(:,:)
     real(r8), pointer :: porgm(:,:)
+    real(r8), pointer :: leafc(:)
     integer :: i,j,jj,g,c,cc,offset
     integer :: n_c
     character (len = 34) :: fn    !TSMP-PDAF: function name for state vector output
@@ -405,6 +430,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 < 0) THEN
@@ -429,6 +456,41 @@ subroutine set_clm_statevec(tstartcycle, mype)
       error stop "Not implemented clmupdate_swc.eq.2"
     endif
 
+
+    ! LAI assimilation
+    ! 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
+       do i=clm_begp,clm_endp ! iterate through patches
+           ! set the leaf area index based on leafC and SLA
+           ! Eq 3 from Thornton and Zimmerman, 2007, J Clim, 20, 3902-3923.
+           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%iytpe(i)) * leafc(i)
+           endif
+           tlai(i) = max(0._r8, tlai(i)) ! don't allow negative LAI
+           ! Add to grid cell average of patch with weighted tlai
+           clm_statevec(patch%gridell(i)) = clm_statevec(patch%gridcell(i)) + patch%wtgcell(i)*tlai(i)
+       end do
+       ! Require second loop to calculate fraction from collected grid cell average
+       do i=clm_begp,clm_endp
+           if (clm_statevec(patch%gridcell(i)) > 0._r8) then ! divide by zero protection
+                old_frac(i) = patch%wtgcell(i)*tlai(i) / clm_statevec(patch%gridcell(i))
+           else
+                old_frac(i) = 0._r8
+           end if
+       end do
+    endif
+
+    if (clmupdate_lai_params==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
+
     !hcp  LAI
     if(clmupdate_T==1) then
       error stop "Not implemented: clmupdate_T.eq.1"
@@ -540,6 +602,9 @@ 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
+    use pftconMod       , only : pftcon
+    use PatchType       , only : patch
 
     implicit none
 
@@ -548,6 +613,8 @@ 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
 
     integer :: i
     character (len = 31) :: fn    !TSMP-PDAF: function name for state vector output
@@ -558,6 +625,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
 
     swc_zero_before_update = .false.
 
+    leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch
+    leafn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%leafn_patch
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN
       ! TSMP-PDAF: For debug runs, output the state vector in files
@@ -619,6 +689,44 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
       call update_clm_texture(tstartcycle, mype)
     endif
 
+    ! LAI assimilation:
+    ! Case 1:
+    if(clmupdate_lai==1) then
+      do i=clm_begp,clm_endp
+      ! At this time tlai(patches) contains the pre-DA patch level LAI
+      ! clm_statevec(gridcells) contains the post-DA grid level LAI
+      ! old_frac(patches) contains the fraction of the pre-DA patch lai / grid cell average lai
+      ! To update tlai(patches) = new gridcell average * old_frac(patches) / patchweight.
+      ! Then use that tlai to determine leafc(patches)
+        if (patch%wtgcell(i) > 0._r8) then ! divide by zero protection
+            tlai(i) = clm_statevec(patch%gridcell(i)) * old_frac(i)/patch%wtgcell(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
+        if (tlai(i) > 0._r8) then ! invalid log protection
+          if (pftcon%dsladlai(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) = tlai(i) * pftcon%slatop(patch%itype(i))
+          endif
+        else
+          leafc(i) = 0._r8
+        endif
+        leafc(i) = max(0._r8, leafc(i)) ! Don't allow negative leaf carbon
+        leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i))
+      end do
+
+      if (clmupdate_lai_params==1) then
+        cc = 1
+        do i=clm_begp,clm_endp
+            pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize+offset)
+            cc = cc + 1
+        end do
+      endif
+    end if
+
   end subroutine update_clm
 
 

From b2604ae3346628cae07c6d9fc28c32d5c4b19f46 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 21 Nov 2025 13:54:46 +0100
Subject: [PATCH 02/25] fortitude: line break overlong lines

---
 interface/model/eclm/enkf_clm_mod_5.F90 | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index ee985acf..a4a7dc57 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -465,7 +465,8 @@ subroutine set_clm_statevec(tstartcycle, mype)
            ! set the leaf area index based on leafC and SLA
            ! Eq 3 from Thornton and Zimmerman, 2007, J Clim, 20, 3902-3923.
            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))
+               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%iytpe(i)) * leafc(i)
            endif
@@ -707,7 +708,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
         ! reformulate the equation to solve for leafc
         if (tlai(i) > 0._r8) then ! invalid log protection
           if (pftcon%dsladlai(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))
+            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) = tlai(i) * pftcon%slatop(patch%itype(i))
           endif

From af61ea31dc2d2f90ab89444cff65147df74f8e04 Mon Sep 17 00:00:00 2001
From: Lukas Strebel <l.strebel@fz-juelich.de>
Date: Fri, 21 Nov 2025 15:26:25 +0100
Subject: [PATCH 03/25] LAI-DA: Updates

---
 interface/model/eclm/enkf_clm_mod_5.F90 | 187 ++++++++++++++++++++----
 1 file changed, 162 insertions(+), 25 deletions(-)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index a4a7dc57..8ff85aeb 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -51,10 +51,14 @@ 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(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
 #endif
   integer(c_int),bind(C,name="clmprint_et")       :: clmprint_et
   integer(c_int),bind(C,name="clmstatevec_allcol")       :: clmstatevec_allcol
@@ -149,15 +153,42 @@ subroutine define_clm_statevec(mype)
         ! Allocate array to store fraction of LAI the patch represents.
         if (allocated(old_frac)) deallocate(old_frac)
         allocate(old_frac(clm_endp-clm_begp+1))
+
+        ! #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
+        !(clm_endp-clm_begp+1) ! 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)
-            clm_statevecsize = (clm_endg - clm_begg+1)*1 ! 1 parameters
-        endif
+            ! 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
+
     !hcp LST DA
     if(clmupdate_T==1) then
       error stop "Not implemented: clmupdate_T.eq.1"
@@ -410,6 +441,8 @@ subroutine set_clm_statevec(tstartcycle, mype)
     use PatchType  , only : patch
     use pftconMod  , only : pftcon
     use shr_kind_mod, only: r8 => shr_kind_r8
+    use PhotosynthesisMod, only : params_inst
+
     implicit none
     integer,intent(in) :: tstartcycle
     integer,intent(in) :: mype
@@ -458,23 +491,27 @@ subroutine set_clm_statevec(tstartcycle, mype)
 
 
     ! 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
+      clm_statevec(:) = 0._r8 ! reset statevector to 0 because we add to it for the average
        do i=clm_begp,clm_endp ! iterate through patches
-           ! set the leaf area index based on leafC and SLA
+           ! 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))
+             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%iytpe(i)) * leafc(i)
-           endif
+             tlai(i) = pftcon%slatop(patch%itype(i)) * leafc(i)
+           end if
            tlai(i) = max(0._r8, tlai(i)) ! don't allow negative LAI
            ! Add to grid cell average of patch with weighted tlai
-           clm_statevec(patch%gridell(i)) = clm_statevec(patch%gridcell(i)) + patch%wtgcell(i)*tlai(i)
+           clm_statevec(patch%gricdell(i)) = clm_statevec(patch%gridcell(i)) + patch%wtgcell(i)*tlai(i)
        end do
+
        ! Require second loop to calculate fraction from collected grid cell average
+       ! second loop once grid cell average is known calculate fraction
        do i=clm_begp,clm_endp
            if (clm_statevec(patch%gridcell(i)) > 0._r8) then ! divide by zero protection
                 old_frac(i) = patch%wtgcell(i)*tlai(i) / clm_statevec(patch%gridcell(i))
@@ -484,6 +521,21 @@ subroutine set_clm_statevec(tstartcycle, mype)
        end do
     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
+
     if (clmupdate_lai_params==1) then
       cc = 1
       do i=clm_begp,clm_endp
@@ -492,6 +544,16 @@ subroutine set_clm_statevec(tstartcycle, mype)
       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"
@@ -606,6 +668,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     use clm_instMod, only : bgc_vegetation_inst
     use pftconMod       , only : pftcon
     use PatchType       , only : patch
+    use PhotosynthesisMod, only : params_inst
 
     implicit none
 
@@ -617,17 +680,24 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     real(r8), pointer :: leafc(:)         ! leaf carbon of patch
     real(r8), pointer :: leafn(:)         ! leaf nitrogen of patch
 
+!    real(r8) :: tlai(clm_begp:clm_endp)
+    real(r8) :: incr_lai
+
     integer :: i
+    integer :: cc
+    integer :: offset
     character (len = 31) :: fn    !TSMP-PDAF: function name for state vector output
+    character (len = 31) :: fn2    !TSMP-PDAF: function name for state vector outpu
     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 = 32) :: fn7    !TSMP-PDAF: function name for state vector outpu
 
     logical :: swc_zero_before_update
 
     swc_zero_before_update = .false.
 
-    leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch
-    leafn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%leafn_patch
+    cc = 1
+    offset = 0
 
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN
@@ -644,6 +714,11 @@ 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
+
+
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN
 
@@ -693,19 +768,32 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     ! LAI assimilation:
     ! Case 1:
     if(clmupdate_lai==1) then
+      ! all patches
       do i=clm_begp,clm_endp
-      ! At this time tlai(patches) contains the pre-DA patch level LAI
-      ! clm_statevec(gridcells) contains the post-DA grid level LAI
-      ! old_frac(patches) contains the fraction of the pre-DA patch lai / grid cell average lai
-      ! To update tlai(patches) = new gridcell average * old_frac(patches) / patchweight.
-      ! Then use that tlai to determine leafc(patches)
+        ! DEBUG ONLY
+        print *, "LST DEBUG leafc,tlai(p) BEFORE:",i,";",leafc(i),";",tlai(i)
+
+        ! At this time tlai(patches) contains the pre-DA patch level LAI
+        ! clm_statevec(gridcells) contains the post-DA grid level LAI
+        ! old_frac(patches) contains the fraction of the pre-DA patch lai / grid cell average lai
+        ! To update tlai(patches) = new gridcell average * old_frac(patches) / patchweight.
+        ! Then use that tlai to determine leafc(patches)
+
+        !DEBUG
+        print *, "LST DEBUG old_frac BEFORE", i, ";", old_frac(i)
+
         if (patch%wtgcell(i) > 0._r8) then ! divide by zero protection
-            tlai(i) = clm_statevec(patch%gridcell(i)) * old_frac(i)/patch%wtgcell(i)
+          tlai(i) = clm_statevec(patch%gridcell(i)) * old_frac(i)/patch%wtgcell(i)
         else
             tlai(i) = 0._r8
         endif
+        ! DEBUG
+        print *, "LST DEBUG calc tlai BEFORE", i, ";", &
+          (pftcon%slatop(patch%itype(i))*exp(leafc(i)*pftcon%dsladlai(patch%itype(i)) - 1._r8)) &
+          /pftcon%dsladlai(patch%itype(i))
         ! 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
             leafc(i) = log(((tlai(i) * pftcon%dsladlai(patch%itype(i))) / pftcon%slatop(patch%itype(i))) + 1.0_r8) &
@@ -718,15 +806,64 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
         endif
         leafc(i) = max(0._r8, leafc(i)) ! Don't allow negative leaf carbon
         leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i))
+
+        ! DEBUG ONLY
+        print *, "LST DEBUG leafc,tlai(wt) AFTER: ",i,";",leafc(i),";",tlai(i)
+      end do
+
+#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=fn2, action="write")
+          WRITE (71,"(es22.15)") leafc(:)
+          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))
+        leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i))
+        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
+    endif
+
+    if (clmupdate_lai_params==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
 
-      if (clmupdate_lai_params==1) then
-        cc = 1
-        do i=clm_begp,clm_endp
-            pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize+offset)
-            cc = cc + 1
-        end do
-      endif
     end if
 
   end subroutine update_clm

From 59f1fb517b825d98846068b4f52e5f3473fe4039 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 21 Nov 2025 15:41:00 +0100
Subject: [PATCH 04/25] LAI-DA: observation operator

---
 interface/framework/obs_op_pdaf.F90 | 32 +++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90
index 1ce5db81..1d4900e0 100644
--- a/interface/framework/obs_op_pdaf.F90
+++ b/interface/framework/obs_op_pdaf.F90
@@ -71,6 +71,8 @@ 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_mode, &
+        ONLY : clmupdate_lai, clm_begp, clm_endp, clm_patch2gc, clm_patchwt
 #ifdef CLMFIVE
    USE clm_instMod, &
      ONLY : soilstate_inst
@@ -128,6 +130,36 @@ 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
+    avesm = 0.0 ! use avesm as grid cell average collecter
+    do j = 1, clm_endp-clm_begp ! 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
+          avesm = avesm + 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
+          avesm = avesm +  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 ', avesm
+    enddo ! all patches for obs checked
+    m_state_p(i) = avesm ! == 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.
+
+
+
 #if defined CLMSA
 if (clmupdate_T==1) then
 

From e522ae95a1d7f7221f1209f1c58b979da870e2a8 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 21 Nov 2025 15:42:14 +0100
Subject: [PATCH 05/25] LAI-DA: wrapper_tsmp.c

---
 interface/model/wrapper_tsmp.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c
index 56687e03..c20d9a4b 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))){
+  if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_lai != 0))){
     update_clm(&tstartcycle, &mype_world);
     if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2){
       print_update_clm(&tcycle, &total_steps);

From a1d577d2bb3ea7b1d59ecec20959d23ff0191289 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 21 Nov 2025 15:47:41 +0100
Subject: [PATCH 06/25] typo fix

---
 interface/model/eclm/enkf_clm_mod_5.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 8ff85aeb..23ad16a0 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -507,7 +507,7 @@ subroutine set_clm_statevec(tstartcycle, mype)
            end if
            tlai(i) = max(0._r8, tlai(i)) ! don't allow negative LAI
            ! Add to grid cell average of patch with weighted tlai
-           clm_statevec(patch%gricdell(i)) = clm_statevec(patch%gridcell(i)) + patch%wtgcell(i)*tlai(i)
+           clm_statevec(patch%gridcell(i)) = clm_statevec(patch%gridcell(i)) + patch%wtgcell(i)*tlai(i)
        end do
 
        ! Require second loop to calculate fraction from collected grid cell average

From bced8e57f6f56b7d341577ba695bd4dfd6e53fb2 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 21 Nov 2025 15:48:20 +0100
Subject: [PATCH 07/25] comment clmupdate_lai_params == 3
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Error:

```
      |                                                                 1
Error: Component ‘kmax’ at (1) is a PRIVATE component of
‘photo_params_type’
```

https://github.com/HPSCTerrSys/TSMP2/actions/runs/19573560213/job/56052817373
---
 interface/model/eclm/enkf_clm_mod_5.F90 | 52 ++++++++++++-------------
 1 file changed, 26 insertions(+), 26 deletions(-)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 23ad16a0..a4821e32 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -170,10 +170,10 @@ subroutine define_clm_statevec(mype)
             ! 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
+        ! if (clmupdate_lai_params==3) then
+        !     clm_varsize      =  (clm_endp-clm_begp+1)
+        !     clm_statevecsize =  3*clm_varsize
+        ! end if
 
     end if
 
@@ -544,15 +544,15 @@ subroutine set_clm_statevec(tstartcycle, mype)
       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
+    ! 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
@@ -849,22 +849,22 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
       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)
+    ! 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
+    !     ! 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
+    !     cc = cc + 1
+    !   end do
 
-    end if
+    ! end if
 
   end subroutine update_clm
 

From 887bc191572d53fcd0ce36df877c4b7396001dae Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 21 Nov 2025 16:40:37 +0100
Subject: [PATCH 08/25] typo fix

---
 interface/framework/obs_op_pdaf.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90
index 1d4900e0..5966b9fa 100644
--- a/interface/framework/obs_op_pdaf.F90
+++ b/interface/framework/obs_op_pdaf.F90
@@ -71,7 +71,7 @@ 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_mode, &
+   USE enkf_clm_mod, &
         ONLY : clmupdate_lai, clm_begp, clm_endp, clm_patch2gc, clm_patchwt
 #ifdef CLMFIVE
    USE clm_instMod, &

From 203d3962c08fe5193faee4d49d88d03e80e48f00 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Thu, 26 Mar 2026 11:33:33 +0100
Subject: [PATCH 09/25] doc for LAI-DA

Generated with Claude Code
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 docs/_toc.yml                                 |   1 +
 docs/users_guide/running_tsmp_pdaf/README.md  |   2 +
 .../running_tsmp_pdaf/input_enkfpf.md         |  48 +++++
 docs/users_guide/running_tsmp_pdaf/lai_da.md  | 164 ++++++++++++++++++
 4 files changed, 215 insertions(+)
 create mode 100644 docs/users_guide/running_tsmp_pdaf/lai_da.md

diff --git a/docs/_toc.yml b/docs/_toc.yml
index 19e34413..178fcab1 100644
--- a/docs/_toc.yml
+++ b/docs/_toc.yml
@@ -28,6 +28,7 @@ parts:
       - file: users_guide/running_tsmp_pdaf/input_cmd
       - file: users_guide/running_tsmp_pdaf/input_obs
       - file: users_guide/running_tsmp_pdaf/input_enkfpf
+      - file: users_guide/running_tsmp_pdaf/lai_da
 
     - file: users_guide/debugging_tsmp_pdaf/README
       title: Debugging TSMP-PDAF
diff --git a/docs/users_guide/running_tsmp_pdaf/README.md b/docs/users_guide/running_tsmp_pdaf/README.md
index 7d1809ab..a3b03a83 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).
 
 See the Virtual Machine download on webpage
 <https://datapub.fz-juelich.de/slts/tsmp-vm/index.html>.
diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
index 790b4aa0..926c8f78 100644
--- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
+++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
@@ -54,6 +54,8 @@ nprocs      =
 update_swc  =
 update_texture  =
 update_T  =
+update_lai  =
+update_lai_params =
 print_swc   =
 print_et   =
 statevec_allcol =
@@ -490,6 +492,50 @@ 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: 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.
+
+(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: Only 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.
+
+-  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.
+
 ### CLM:print_swc ###
 
 `CLM:print_swc`: (integer) If set to `1`, the updated soil moisture
@@ -892,6 +938,8 @@ Default: 0, output turned off.
  |           | `problemname`           | \-            |
  |           | `nprocs`                | 0             |
  |           | `update_swc`            | 1             |
+ |           | `update_lai`            | 0             |
+ |           | `update_lai_params`     | 0             |
  |           | `print_swc`             | 0             |
  |           | `print_et`              | 0             |
  |           | `statevec_allcol`       | 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..31f579ac
--- /dev/null
+++ b/docs/users_guide/running_tsmp_pdaf/lai_da.md
@@ -0,0 +1,164 @@
+(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 two 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.
+
+## 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, kg C m$^{-2}$),
+$s_{\text{top}}$ is `slatop` (specific leaf area at the canopy top,
+m$^2$ kg C$^{-1}$) and $s'$ is `dsladlai` (slope of specific leaf
+area with respect to LAI, m$^2$ kg C$^{-1}$ LAI$^{-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$ |
+
+$N_g$ = number of local grid cells, $N_p$ = number of local patches.
+
+## 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. The relative contribution of each patch is stored as
+   $f_p = w_p\,\text{LAI}(p)\,/\,\text{LAI}_g$, which is used in the
+   update phase.
+
+**PDAF update**: PDAF operates on the gridcell-level LAI values.
+
+**Update phase** (`update_clm`):
+
+1. The updated gridcell LAI from the state vector is distributed back
+   to patches: $\text{LAI}'(p) = \text{LAI}'_g\,f_p\,/\,w_p$.
+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 ##
+
+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:
+
+- `leafc` is read back from block 1 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.
+
+## Configuration Examples ##
+
+**Example 1** — gridcell LAI assimilation, no parameter estimation:
+
+```ini
+[CLM]
+update_lai        = 1
+update_lai_params = 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
+```

From 25e098fd7ca961bcb443636644ea74fe52cd6121 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Thu, 26 Mar 2026 11:45:17 +0100
Subject: [PATCH 10/25] fortitude fix

---
 interface/framework/obs_op_pdaf.F90 | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90
index 5966b9fa..4194110b 100644
--- a/interface/framework/obs_op_pdaf.F90
+++ b/interface/framework/obs_op_pdaf.F90
@@ -147,7 +147,8 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
       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
-          avesm = avesm + 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
+          avesm = avesm + 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
           avesm = avesm +  clm_patchwt(j) *(state_p(j+1*clm_varsize)*state_p(j)) ! 2nd formula
         endif ! dsladlai decider

From 69e222baf421535a38e36d4df7463ccd90be43f0 Mon Sep 17 00:00:00 2001
From: zhao7777 <zhao2@jwlogin10.juwels>
Date: Wed, 1 Apr 2026 16:01:19 +0200
Subject: [PATCH 11/25] Fix divide-by-zero in domain_def_clm for degenerate 1x1
 domains

---
 interface/model/eclm/enkf_clm_mod_5.F90 | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index b70734dc..a8324e0a 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -1313,6 +1313,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
@@ -1371,6 +1373,11 @@ subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, &
     maxlon = MAXVAL(lon(:) + 180)
     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), stat=ier)
@@ -1378,7 +1385,16 @@ subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, &
     allocate(latixy_obs(dim_obs), stat=ier)
 
     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
@@ -1395,6 +1411,7 @@ subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, &
           latixy_obs(i) = 1
        endif
     end do
+    
     ! deallocate temporary arrays
     !deallocate(longxy)
     !deallocate(latixy)

From 93ec5208687f737220e16c140080c52a21a711b5 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Thu, 2 Apr 2026 13:59:10 +0200
Subject: [PATCH 12/25] LAI-DA debug output

---
 interface/model/eclm/enkf_clm_mod_5.F90 | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index a8324e0a..8c452b84 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -737,6 +737,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
 
@@ -838,6 +846,17 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
 
         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=fn2, action="write")
+          WRITE (71,"(es22.15)") leafc(:)
+          CLOSE(71)
+        END IF
+#endif
+
     endif
 
     if (clmupdate_lai_params==1) then

From c6fb0dafdfb5ffdc2551bc9292113c18dca9ab96 Mon Sep 17 00:00:00 2001
From: Haojin Zhao <46377736+zhao7777@users.noreply.github.com>
Date: Thu, 2 Apr 2026 16:50:41 +0200
Subject: [PATCH 13/25] Merge pull request #66 from
 zhao7777/fix-lai-obs-operator

Fix LAI DA obs operator loop and leafc debug I/O
---
 interface/framework/obs_op_pdaf.F90     | 14 +++++++-------
 interface/model/eclm/enkf_clm_mod_5.F90 |  6 +++---
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90
index 4194110b..4fefd0f3 100644
--- a/interface/framework/obs_op_pdaf.F90
+++ b/interface/framework/obs_op_pdaf.F90
@@ -102,7 +102,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 ! soide=(/0.d0,  0.02d0,  0.05d0,  0.1d0,  0.17d0, 0.3d0,  0.5d0, &
 !                0.8d0,   1.3d0,   2.d0,  3.d0, 5.d0,  12.d0/) !soil depth
 
-real :: tot, avesm, avesm_temp, Dp
+real :: tot, avesm, avelai, avesm_temp, Dp
 integer :: nsc
 ! end of hcp
 
@@ -141,21 +141,21 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
   ! this should not happen, but would cause div by zero.
 
   do i = 1, dim_obs_p ! loop over all observations
-    avesm = 0.0 ! use avesm as grid cell average collecter
-    do j = 1, clm_endp-clm_begp ! loop over all patches
+    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
-          avesm = avesm + clm_patchwt(j) * ((state_p(j+1*clm_varsize)&
+          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
-          avesm = avesm +  clm_patchwt(j) *(state_p(j+1*clm_varsize)*state_p(j)) ! 2nd formula
+          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 ', avesm
+      write(*,*) 'DEBUG LAI : average ', avelai
     enddo ! all patches for obs checked
-    m_state_p(i) = avesm ! == lai of gridcell where the observation is
+    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.
 
diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 8c452b84..9608f363 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -741,7 +741,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
         ! 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(:,:)
+        WRITE (71,"(es22.15)") leafc(:)
         CLOSE(71)
       END IF
 
@@ -824,7 +824,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
         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=fn2, action="write")
+          OPEN(unit=71, file=fn7, action="write")
           WRITE (71,"(es22.15)") leafc(:)
           CLOSE(71)
         END IF
@@ -851,7 +851,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
         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=fn2, action="write")
+          OPEN(unit=71, file=fn7, action="write")
           WRITE (71,"(es22.15)") leafc(:)
           CLOSE(71)
         END IF

From 70038451add2e2aec2c62ea9715c6e838b7f15fc Mon Sep 17 00:00:00 2001
From: Haojin Zhao <46377736+zhao7777@users.noreply.github.com>
Date: Fri, 10 Apr 2026 10:05:22 +0200
Subject: [PATCH 14/25] LAI DA debug dumps + fix indexing and leafc inverse
 (#73)

Co-authored-by: zhao7777 <zhao2@jwlogin10.juwels>
---
 interface/model/eclm/enkf_clm_mod_5.F90 | 150 +++++++++++++++---------
 1 file changed, 95 insertions(+), 55 deletions(-)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 9608f363..8633da46 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -92,7 +92,7 @@ module enkf_clm_mod
 
   ! Arrays for LAI Assimilation
   real(r8),allocatable :: tlai(:)
-  real(r8),allocatable :: old_frac(:)
+  real(r8),allocatable :: lai_wtsum2_gc(:)
 
   contains
 
@@ -150,15 +150,15 @@ subroutine define_clm_statevec(mype)
     if(clmupdate_lai==1) then
         ! Allocate array to store total leaf area index per patch
         if (allocated(tlai)) deallocate(tlai)
-        allocate(tlai(clm_endp-clm_begp+1))
+        allocate(tlai(clm_begp:clm_endp))
         ! Allocate array to store fraction of LAI the patch represents.
-        if (allocated(old_frac)) deallocate(old_frac)
-        allocate(old_frac(clm_endp-clm_begp+1))
+        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
-        !(clm_endp-clm_begp+1) ! Currently no combination of SWC and LAI DA
+        ! 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
@@ -429,12 +429,14 @@ 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(lai_wtsum2_gc)) deallocate(lai_wtsum2_gc)
+    IF (allocated(tlai)) deallocate(tlai)
 
   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
+    use clm_instMod, only : bgc_vegetation_inst, canopystate_inst
     use clm_varpar   , only : nlevsoi
     ! use clm_varcon, only: nameg, namec
     ! use GetGlobalValuesMod, only: GetGlobalWrite
@@ -452,10 +454,11 @@ subroutine set_clm_statevec(tstartcycle, mype)
     real(r8), pointer :: pclay(:,:)
     real(r8), pointer :: porgm(:,:)
     real(r8), pointer :: leafc(:)
-    integer :: i,j,jj,g,c,cc,offset
+    integer :: i,j,jj,g,c,cc,offset,igc
     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
@@ -490,12 +493,18 @@ subroutine set_clm_statevec(tstartcycle, mype)
       error stop "Not implemented clmupdate_swc.eq.2"
     endif
 
-
+    !> @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.
@@ -507,19 +516,39 @@ subroutine set_clm_statevec(tstartcycle, mype)
              tlai(i) = pftcon%slatop(patch%itype(i)) * leafc(i)
            end if
            tlai(i) = max(0._r8, tlai(i)) ! don't allow negative LAI
-           ! Add to grid cell average of patch with weighted tlai
-           clm_statevec(patch%gridcell(i)) = clm_statevec(patch%gridcell(i)) + patch%wtgcell(i)*tlai(i)
-       end do
-
-       ! Require second loop to calculate fraction from collected grid cell average
-       ! second loop once grid cell average is known calculate fraction
-       do i=clm_begp,clm_endp
-           if (clm_statevec(patch%gridcell(i)) > 0._r8) then ! divide by zero protection
-                old_frac(i) = patch%wtgcell(i)*tlai(i) / clm_statevec(patch%gridcell(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
-                old_frac(i) = 0._r8
+             write(*,*) 'DEBUG LAI : igc out of bounds in LAI mapping: ', igc, &
+               ' valid range: 1-', clm_statevecsize, ' i=', i
            end if
        end do
+
+#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
@@ -666,7 +695,7 @@ 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
+    use clm_instMod, only : bgc_vegetation_inst, canopystate_inst
     use pftconMod       , only : pftcon
     use PatchType       , only : patch
     use PhotosynthesisMod, only : params_inst
@@ -687,11 +716,12 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     integer :: i
     integer :: cc
     integer :: offset
+    integer :: igc
     character (len = 31) :: fn    !TSMP-PDAF: function name for state vector output
-    character (len = 31) :: fn2    !TSMP-PDAF: function name for state vector outpu
-    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 = 32) :: fn7    !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
 
@@ -777,55 +807,61 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     ! LAI assimilation:
     ! Case 1:
     if(clmupdate_lai==1) then
-      ! all patches
+      ! Distribute post-DA grid-cell LAI (clm_statevec(g)) to patches using patch weights.
+      ! Conservation: sum_p wtgcell(p)*tlai(p) = L_an  if
+      !   tlai(p) = L_an * wtgcell(p) / sum_k wtgcell(k)^2  (sum over patches k in grid cell g).
       do i=clm_begp,clm_endp
-        ! DEBUG ONLY
-        print *, "LST DEBUG leafc,tlai(p) BEFORE:",i,";",leafc(i),";",tlai(i)
-
-        ! At this time tlai(patches) contains the pre-DA patch level LAI
-        ! clm_statevec(gridcells) contains the post-DA grid level LAI
-        ! old_frac(patches) contains the fraction of the pre-DA patch lai / grid cell average lai
-        ! To update tlai(patches) = new gridcell average * old_frac(patches) / patchweight.
-        ! Then use that tlai to determine leafc(patches)
-
-        !DEBUG
-        print *, "LST DEBUG old_frac BEFORE", i, ";", old_frac(i)
-
-        if (patch%wtgcell(i) > 0._r8) then ! divide by zero protection
-          tlai(i) = clm_statevec(patch%gridcell(i)) * old_frac(i)/patch%wtgcell(i)
+        ! Map CLM gridcell index to PE-local clm_statevec index (same as SWC: col%gridcell - begg + 1)
+        igc = patch%gridcell(i) - clm_begg + 1
+        if (igc >= 1 .and. igc <= clm_statevecsize .and. patch%wtgcell(i) > 0._r8 &
+            .and. lai_wtsum2_gc(patch%gridcell(i)) > 0._r8) then
+          tlai(i) = clm_statevec(igc) * patch%wtgcell(i) &
+            / lai_wtsum2_gc(patch%gridcell(i))
         else
-            tlai(i) = 0._r8
+          tlai(i) = 0._r8
         endif
-        ! DEBUG
-        print *, "LST DEBUG calc tlai BEFORE", i, ";", &
-          (pftcon%slatop(patch%itype(i))*exp(leafc(i)*pftcon%dsladlai(patch%itype(i)) - 1._r8)) &
-          /pftcon%dsladlai(patch%itype(i))
         ! 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
-            leafc(i) = log(((tlai(i) * pftcon%dsladlai(patch%itype(i))) / pftcon%slatop(patch%itype(i))) + 1.0_r8) &
-              / pftcon%dsladlai(patch%itype(i))
+            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
-              leafc(i) = tlai(i) * pftcon%slatop(patch%itype(i))
+              ! 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
-        leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i))
-
+        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 *, "LST DEBUG leafc,tlai(wt) AFTER: ",i,";",leafc(i),";",tlai(i)
+        ! print *, "LAI1 DEBUG leafc,tlai(wt) AFTER: ", i, ";", leafc(i), ";", tlai(i)
       end do
 
 #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(:)
+        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
@@ -838,7 +874,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
       do i=clm_begp,clm_endp
         leafc(i) = clm_statevec(cc)
         leafc(i) = max(0._r8, leafc(i))
-        leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(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)

From cf08175870231e76bc9550d28ca279384e77d31e Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 10 Apr 2026 13:13:01 +0200
Subject: [PATCH 15/25] docs: LAIDA docs update

---
 docs/users_guide/running_tsmp_pdaf/lai_da.md | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/docs/users_guide/running_tsmp_pdaf/lai_da.md b/docs/users_guide/running_tsmp_pdaf/lai_da.md
index 31f579ac..c0eb297f 100644
--- a/docs/users_guide/running_tsmp_pdaf/lai_da.md
+++ b/docs/users_guide/running_tsmp_pdaf/lai_da.md
@@ -33,13 +33,16 @@ $$
 \end{cases}
 $$
 
-where $C_\ell$ is `leafc` (leaf carbon, kg C m$^{-2}$),
-$s_{\text{top}}$ is `slatop` (specific leaf area at the canopy top,
-m$^2$ kg C$^{-1}$) and $s'$ is `dsladlai` (slope of specific leaf
-area with respect to LAI, m$^2$ kg C$^{-1}$ LAI$^{-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).
+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

From 8df9a4e833e39b310c474a9e2f13624c81c2a04d Mon Sep 17 00:00:00 2001
From: Haojin Zhao <46377736+zhao7777@users.noreply.github.com>
Date: Mon, 18 May 2026 20:03:10 +0200
Subject: [PATCH 16/25] Fix LAI DA update and add additive/multiplicative
 weighting incr. (#80)

Co-authored-by: zhao7777 <zhao2@jwlogin10.juwels>
---
 .../running_tsmp_pdaf/input_enkfpf.md         | 22 ++++++
 docs/users_guide/running_tsmp_pdaf/lai_da.md  | 22 +++---
 interface/model/common/enkf.h                 |  1 +
 interface/model/common/read_enkfpar.c         |  1 +
 interface/model/eclm/enkf_clm_mod_5.F90       | 67 ++++++++++++++-----
 5 files changed, 89 insertions(+), 24 deletions(-)

diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
index 926c8f78..a12e0459 100644
--- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
+++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
@@ -56,6 +56,7 @@ update_texture  =
 update_T  =
 update_lai  =
 update_lai_params =
+update_lai_incr_w =
 print_swc   =
 print_et   =
 statevec_allcol =
@@ -536,6 +537,26 @@ LAI. Only takes effect if `CLM:update_lai` is non-zero. Default: `0`.
    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
@@ -940,6 +961,7 @@ Default: 0, output turned off.
  |           | `update_swc`            | 1             |
  |           | `update_lai`            | 0             |
  |           | `update_lai_params`     | 0             |
+ |           | `update_lai_incr_w`     | 1.0           |
  |           | `print_swc`             | 0             |
  |           | `print_et`              | 0             |
  |           | `statevec_allcol`       | 0             |
diff --git a/docs/users_guide/running_tsmp_pdaf/lai_da.md b/docs/users_guide/running_tsmp_pdaf/lai_da.md
index c0eb297f..d6b624d8 100644
--- a/docs/users_guide/running_tsmp_pdaf/lai_da.md
+++ b/docs/users_guide/running_tsmp_pdaf/lai_da.md
@@ -10,12 +10,14 @@ Not available for CLM3.5 or CLM5.0 without BGC.
 
 ## Configuration ##
 
-LAI DA is controlled by two parameters in `enkfpf.par`:
+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.
 
 ## Background ##
 
@@ -74,16 +76,17 @@ cell (a weighted average over all patches in that cell):
 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. The relative contribution of each patch is stored as
-   $f_p = w_p\,\text{LAI}(p)\,/\,\text{LAI}_g$, which is used in the
-   update phase.
+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 updated gridcell LAI from the state vector is distributed back
-   to patches: $\text{LAI}'(p) = \text{LAI}'_g\,f_p\,/\,w_p$.
+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`:
 
 $$
@@ -125,8 +128,10 @@ where $g_i$ is the grid cell corresponding to observation $i$.
 
 After the PDAF update:
 
-- `leafc` is read back from block 1 and clipped to zero;
-  `leafn` is updated consistently.
+- 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
@@ -140,6 +145,7 @@ After the PDAF update:
 [CLM]
 update_lai        = 1
 update_lai_params = 0
+update_lai_incr_w = 1.0
 ```
 
 **Example 2** — gridcell LAI assimilation with joint `slatop` estimation:
diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h
index 1dcd2498..5bb86638 100755
--- a/interface/model/common/enkf.h
+++ b/interface/model/common/enkf.h
@@ -90,6 +90,7 @@ GLOBAL int clmupdate_T;
 GLOBAL int clmupdate_texture;
 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 clmstatevec_allcol;
diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c
index 9cf5befd..2181966e 100755
--- a/interface/model/common/read_enkfpar.c
+++ b/interface/model/common/read_enkfpar.c
@@ -81,6 +81,7 @@ void read_enkfpar(char *parname)
   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);
   clmstatevec_allcol    = iniparser_getint(pardict,"CLM:statevec_allcol",0);
diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 8633da46..677d092d 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -60,6 +60,7 @@ module enkf_clm_mod
   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
 #endif
   integer(c_int),bind(C,name="clmprint_et")       :: clmprint_et
   integer(c_int),bind(C,name="clmstatevec_allcol")       :: clmstatevec_allcol
@@ -92,6 +93,7 @@ module enkf_clm_mod
 
   ! Arrays for LAI Assimilation
   real(r8),allocatable :: tlai(:)
+  real(r8),allocatable :: tlai_orig(:)
   real(r8),allocatable :: lai_wtsum2_gc(:)
 
   contains
@@ -151,7 +153,8 @@ subroutine define_clm_statevec(mype)
         ! Allocate array to store total leaf area index per patch
         if (allocated(tlai)) deallocate(tlai)
         allocate(tlai(clm_begp:clm_endp))
-        ! Allocate array to store fraction of LAI the patch represents.
+        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))
 
@@ -208,11 +211,10 @@ subroutine define_clm_statevec(mype)
     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 used when
+    ! mapping analysis increments back to CLM (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) then
+    if ((clmupdate_swc/=0 .and. clmstatevec_colmean/=0) .or. (clmupdate_lai/=0)) then
       allocate(clm_statevec_orig(clm_statevecsize))
     end if
 
@@ -431,6 +433,7 @@ subroutine cleanup_clm_statevec()
     IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig)
     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
 
@@ -516,6 +519,7 @@ subroutine set_clm_statevec(tstartcycle, mype)
              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
@@ -530,6 +534,11 @@ subroutine set_clm_statevec(tstartcycle, mype)
                ' 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
@@ -691,6 +700,35 @@ subroutine set_clm_statevec_swc()
   end subroutine set_clm_statevec_swc
 
 
+  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 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
@@ -710,8 +748,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     real(r8), pointer :: leafc(:)         ! leaf carbon of patch
     real(r8), pointer :: leafn(:)         ! leaf nitrogen of patch
 
-!    real(r8) :: tlai(clm_begp:clm_endp)
-    real(r8) :: incr_lai
+    real(r8) :: lai_gc_an
+    real(r8) :: lai_gc_fc
 
     integer :: i
     integer :: cc
@@ -805,18 +843,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     endif
 
     ! LAI assimilation:
-    ! Case 1:
+    ! 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
-      ! Distribute post-DA grid-cell LAI (clm_statevec(g)) to patches using patch weights.
-      ! Conservation: sum_p wtgcell(p)*tlai(p) = L_an  if
-      !   tlai(p) = L_an * wtgcell(p) / sum_k wtgcell(k)^2  (sum over patches k in grid cell g).
       do i=clm_begp,clm_endp
-        ! Map CLM gridcell index to PE-local clm_statevec index (same as SWC: col%gridcell - begg + 1)
         igc = patch%gridcell(i) - clm_begg + 1
-        if (igc >= 1 .and. igc <= clm_statevecsize .and. patch%wtgcell(i) > 0._r8 &
-            .and. lai_wtsum2_gc(patch%gridcell(i)) > 0._r8) then
-          tlai(i) = clm_statevec(igc) * patch%wtgcell(i) &
-            / lai_wtsum2_gc(patch%gridcell(i))
+        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

From 7a9bb7538c3a6fd7932e657d365d5fb0f1ac1c4d Mon Sep 17 00:00:00 2001
From: Haojin Zhao <46377736+zhao7777@users.noreply.github.com>
Date: Thu, 16 Jul 2026 11:35:13 +0200
Subject: [PATCH 17/25] Implement LAI-DA mode 3 carbon-pool update (#83)

* Fix LAI DA update and add additive/multiplicative weighting incr.

* Implement Option 3 to update carbon pools

---------

Co-authored-by: zhao7777 <zhao2@jwlogin10.juwels>
---
 interface/framework/obs_op_pdaf.F90     |  44 +++++-
 interface/model/eclm/enkf_clm_mod_5.F90 | 186 ++++++++++++++++++++++++
 2 files changed, 228 insertions(+), 2 deletions(-)

diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90
index 4fefd0f3..9309e1db 100644
--- a/interface/framework/obs_op_pdaf.F90
+++ b/interface/framework/obs_op_pdaf.F90
@@ -49,7 +49,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 !
 ! !USES:
    USE mod_assimilation, &
-        ONLY: obs_index_p, obs_p
+        ONLY: obs_index_p, obs_p, mype_world
 #ifndef CLMSA
 #ifndef OBS_ONLY_CLM
    USE mod_assimilation, ONLY: sc_p
@@ -72,7 +72,9 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
    USE enkf_clm_mod, &
         ONLY : clm_varsize, clm_paramarr, clmupdate_swc, clmupdate_T, clmcrns_bd
    USE enkf_clm_mod, &
-        ONLY : clmupdate_lai, clm_begp, clm_endp, clm_patch2gc, clm_patchwt
+        ONLY : clmupdate_lai, clm_begp, clm_endp, clm_patch2gc, clm_patchwt, &
+               clm_lai_patch_itype, clmt_printensemble
+   USE pftconMod, only : pftcon
 #ifdef CLMFIVE
    USE clm_instMod, &
      ONLY : soilstate_inst
@@ -104,6 +106,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 
 real :: tot, avesm, avelai, avesm_temp, Dp
 integer :: nsc
+character(len=48) :: fn_hx
 ! end of hcp
 
 #ifndef PARFLOW_STAND_ALONE
@@ -159,6 +162,43 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
   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 fixed slatop/dsladlai from pftcon to compute gridcell LAI.
+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 (pftcon%dsladlai(clm_lai_patch_itype(j)) > 0.0) then
+          avelai = avelai + clm_patchwt(j) * ((pftcon%slatop(clm_lai_patch_itype(j))&
+            *(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) * (pftcon%slatop(clm_lai_patch_itype(j))*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
diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 677d092d..28ac82c8 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -54,6 +54,8 @@ module enkf_clm_mod
   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
@@ -102,6 +104,7 @@ module enkf_clm_mod
   subroutine define_clm_statevec(mype)
     use decompMod , only : get_proc_bounds
     use clm_varpar   , only : nlevsoi
+    use PatchType  , only : patch
 
     implicit none
 
@@ -111,6 +114,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)
@@ -193,6 +197,35 @@ subroutine define_clm_statevec(mype)
 
     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 (.not. patch%is_bareground(i)) clm_varsize = clm_varsize + 1
+      end do
+      clm_statevecsize = 3*clm_varsize
+
+      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 (.not. patch%is_bareground(i)) 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"
@@ -431,6 +464,10 @@ 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)
@@ -457,6 +494,8 @@ subroutine set_clm_statevec(tstartcycle, mype)
     real(r8), pointer :: pclay(:,:)
     real(r8), pointer :: porgm(:,:)
     real(r8), pointer :: leafc(:)
+    real(r8), pointer :: livestemc(:)
+    real(r8), pointer :: deadstemc(:)
     integer :: i,j,jj,g,c,cc,offset,igc
     integer :: n_c
     character (len = 34) :: fn    !TSMP-PDAF: function name for state vector output
@@ -575,6 +614,29 @@ subroutine set_clm_statevec(tstartcycle, mype)
       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
+
+#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) then
       cc = 1
       do i=clm_begp,clm_endp
@@ -729,6 +791,86 @@ subroutine lai_patch_tlai_blend(tlai_fc, lai_gc_fc, lai_gc_an, 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,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"
+    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
@@ -747,11 +889,17 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     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
@@ -785,6 +933,10 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
 
     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
 
 
 
@@ -934,6 +1086,40 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
 
     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
+
+#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) then
       cc = 1
       do i=clm_begp,clm_endp

From 1d26d2e95f5840dd66dd95f106da1708633e7a05 Mon Sep 17 00:00:00 2001
From: Haojin Zhao <46377736+zhao7777@users.noreply.github.com>
Date: Thu, 23 Jul 2026 10:39:34 +0200
Subject: [PATCH 18/25] Bug fix and add parameter update in mode 3 (#85)

Co-authored-by: zhao7777 <zhao2@jwlogin10.juwels>
---
 interface/framework/obs_op_pdaf.F90         |  19 ++--
 interface/model/Makefile                    |   1 +
 interface/model/eclm/enkf_clm_mod_5.F90     |  50 ++++++++--
 interface/model/eclm/print_update_clm_5.F90 | 104 +++++++++++++++++++-
 interface/model/wrapper_tsmp.c              |   4 +
 5 files changed, 157 insertions(+), 21 deletions(-)

diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90
index a73e3cb1..9dce2705 100644
--- a/interface/framework/obs_op_pdaf.F90
+++ b/interface/framework/obs_op_pdaf.F90
@@ -49,7 +49,8 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 !
 ! !USES:
    USE mod_assimilation, &
-        ONLY: obs_index_p, obs_p, mype_world
+        ONLY: obs_index_p, obs_p
+   USE mod_parallel_pdaf, ONLY: mype_world
 #ifndef CLMSA
 #ifndef OBS_ONLY_CLM
    USE mod_assimilation, ONLY: sc_p
@@ -73,7 +74,8 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
         ONLY : clm_varsize, clm_paramarr, clmupdate_swc, clmupdate_T, clmcrns_bd
    USE enkf_clm_mod, &
         ONLY : clmupdate_lai, clm_begp, clm_endp, clm_patch2gc, clm_patchwt, &
-               clm_lai_patch_itype, clmt_printensemble
+               clm_lai_patch_itype, clmt_printensemble, clmupdate_lai_params, &
+               lai_mode3_npool
    USE pftconMod, only : pftcon
 #ifdef CLMFIVE
    USE clm_instMod, &
@@ -104,7 +106,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 ! soide=(/0.d0,  0.02d0,  0.05d0,  0.1d0,  0.17d0, 0.3d0,  0.5d0, &
 !                0.8d0,   1.3d0,   2.d0,  3.d0, 5.d0,  12.d0/) !soil depth
 
-real :: tot, avesm, avelai, avesm_temp, Dp
+real :: tot, avesm, avelai, avesm_temp, Dp, slatop_val
 integer :: nsc
 character(len=48) :: fn_hx
 ! end of hcp
@@ -164,7 +166,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 
 ! LAI assimilation assuming statevec contains leafc, livestemc, deadstemc
 ! for each vegetated patch (bare-ground PFT excluded).
-! H uses fixed slatop/dsladlai from pftcon to compute gridcell LAI.
+! H uses slatop from state (update_lai_params=1) or pftcon; dsladlai from pftcon.
 if (clmupdate_lai==3) then
   lpointobs = .false.
 
@@ -172,12 +174,17 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_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) * ((pftcon%slatop(clm_lai_patch_itype(j))&
+          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) * (pftcon%slatop(clm_lai_patch_itype(j))*state_p(j))
+          avelai = avelai + clm_patchwt(j) * (slatop_val*state_p(j))
         endif
       endif
     enddo
diff --git a/interface/model/Makefile b/interface/model/Makefile
index 185b906a..2c7f2585 100644
--- a/interface/model/Makefile
+++ b/interface/model/Makefile
@@ -1,5 +1,6 @@
 BASEDIR = $(PDAF_DIR)
 include $(BASEDIR)/make.arch/$(PDAF_ARCH).h
+FINCS = -I$(BASEDIR)/include -I$(BASEDIR)/build
 ## includes
 
 FFLAGS = -cpp $(TSMPPDAFIMPORTFLAGS)
diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 89a5dd00..b4754628 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -83,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
@@ -144,6 +149,7 @@ 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
 
@@ -254,9 +260,12 @@ subroutine define_clm_statevec(mype)
     if(clmupdate_lai==3) then
       clm_varsize = 0
       do i=clm_begp,clm_endp
-        if (.not. patch%is_bareground(i)) clm_varsize = clm_varsize + 1
+        if (patch%itype(i) /= noveg) clm_varsize = clm_varsize + 1
       end do
-      clm_statevecsize = 3*clm_varsize
+      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))
@@ -269,7 +278,7 @@ subroutine define_clm_statevec(mype)
 
       cc = 1
       do i=clm_begp,clm_endp
-        if (.not. patch%is_bareground(i)) then
+        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)
@@ -328,7 +337,7 @@ 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)) .or. (clmupdate_lai/=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
@@ -808,7 +817,6 @@ subroutine set_clm_statevec(tstartcycle, mype)
     use PatchType  , only : patch
     use pftconMod  , only : pftcon
     use shr_kind_mod, only: r8 => shr_kind_r8
-    use PhotosynthesisMod, only : params_inst
 
     implicit none
     integer,intent(in) :: tstartcycle
@@ -820,7 +828,7 @@ subroutine set_clm_statevec(tstartcycle, mype)
     real(r8), pointer :: leafc(:)
     real(r8), pointer :: livestemc(:)
     real(r8), pointer :: deadstemc(:)
-    integer :: i,j,jj,g,c,cc,offset,igc
+    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
@@ -954,6 +962,16 @@ subroutine set_clm_statevec(tstartcycle, mype)
         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)
@@ -961,7 +979,7 @@ subroutine set_clm_statevec(tstartcycle, mype)
 #endif
     endif
 
-    if (clmupdate_lai_params==1) then
+    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))
@@ -1426,9 +1444,12 @@ subroutine write_lai_da_mode3_debug(phase, cycle, mype)
     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,i5.5,a)") "lai_da_", mype, ".", trim(phase), ".", cycle, ".txt"
+    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)
@@ -1485,7 +1506,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     use clm_instMod, only : bgc_vegetation_inst, canopystate_inst
     use pftconMod       , only : pftcon
     use PatchType       , only : patch
-    use PhotosynthesisMod, only : params_inst
 
     implicit none
 
@@ -1726,6 +1746,16 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
         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)
@@ -1733,7 +1763,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
 #endif
     endif
 
-    if (clmupdate_lai_params==1) then
+    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)
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 16892c79..aaad257f 100644
--- a/interface/model/wrapper_tsmp.c
+++ b/interface/model/wrapper_tsmp.c
@@ -208,6 +208,10 @@ void update_tsmp(){
 #endif
 
 #ifdef CLMFIVE
+    if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2 ||
+       (clmupdate_lai == 3 && clmupdate_lai_params == 1)){
+      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();

From bb9c5ac63c2c6f9571c3102ba7fd95dcdd9d9e12 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 24 Jul 2026 14:07:59 +0200
Subject: [PATCH 19/25] docs: documentation for option 3 of LAI-DA

---
 .../running_tsmp_pdaf/input_enkfpf.md         | 39 ++++++--
 docs/users_guide/running_tsmp_pdaf/lai_da.md  | 94 +++++++++++++++++--
 2 files changed, 115 insertions(+), 18 deletions(-)

diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
index aacf8857..b4aba368 100644
--- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
+++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md
@@ -521,13 +521,25 @@ See [LAI Data Assimilation](laida) for algorithmic details.
    PDAF update, the updated gridcell LAI is back-transformed to
    per-patch `leafc` and `leafn`.
 
--  2: 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.
+-  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 ###
@@ -538,9 +550,16 @@ LAI. Only takes effect if `CLM:update_lai` is non-zero. Default: `0`.
 
 -  0: No parameter estimation.
 
--  1: Only 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.
+-  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
diff --git a/docs/users_guide/running_tsmp_pdaf/lai_da.md b/docs/users_guide/running_tsmp_pdaf/lai_da.md
index d6b624d8..06a71615 100644
--- a/docs/users_guide/running_tsmp_pdaf/lai_da.md
+++ b/docs/users_guide/running_tsmp_pdaf/lai_da.md
@@ -17,8 +17,9 @@ LAI DA is controlled by three parameters in `enkfpf.par`:
 - [`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.
+  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
@@ -55,14 +56,17 @@ itself.
 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$ |
+| `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_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 ##
 
@@ -108,6 +112,10 @@ 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:
 
@@ -137,6 +145,60 @@ After the PDAF update:
   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:
@@ -171,3 +233,19 @@ update_lai_params = 0
 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
+```

From bee9d6544dbdc685e2f1b0237bd7ca8a3b843466 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 24 Jul 2026 14:09:06 +0200
Subject: [PATCH 20/25] small comment update

---
 interface/model/eclm/enkf_clm_mod_5.F90 | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index b4754628..30020016 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -343,8 +343,10 @@ subroutine define_clm_statevec(mype)
     end if
 
 
-    ! Allocate statevector-duplicate for saving the forecast state used when
-    ! mapping analysis increments back to CLM (SWC DA column-mean mode or LAI DA 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  .or. (clmupdate_lai/=0)) then
       allocate(clm_statevec_orig(clm_statevecsize))

From 0333c9b54e7f153c8196865bfd0acc6ebb851fb5 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 24 Jul 2026 17:41:59 +0200
Subject: [PATCH 21/25] obs_op_pdaf: style changes in module commands and
 declarations

---
 interface/framework/obs_op_pdaf.F90 | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90
index 9dce2705..b7def5b6 100644
--- a/interface/framework/obs_op_pdaf.F90
+++ b/interface/framework/obs_op_pdaf.F90
@@ -50,7 +50,6 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 ! !USES:
    USE mod_assimilation, &
         ONLY: obs_index_p, obs_p
-   USE mod_parallel_pdaf, ONLY: mype_world
 #ifndef CLMSA
 #ifndef OBS_ONLY_CLM
    USE mod_assimilation, ONLY: sc_p
@@ -58,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, &
@@ -72,10 +72,15 @@ 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, clm_begp, clm_endp, clm_patch2gc, clm_patchwt, &
-               clm_lai_patch_itype, clmt_printensemble, clmupdate_lai_params, &
-               lai_mode3_npool
+   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, &
@@ -106,7 +111,8 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 ! soide=(/0.d0,  0.02d0,  0.05d0,  0.1d0,  0.17d0, 0.3d0,  0.5d0, &
 !                0.8d0,   1.3d0,   2.d0,  3.d0, 5.d0,  12.d0/) !soil depth
 
-real :: tot, avesm, avelai, avesm_temp, Dp, slatop_val
+real :: tot, avesm, avesm_temp, Dp
+real :: avelai, slatop_val
 integer :: nsc
 character(len=48) :: fn_hx
 ! end of hcp

From 58e707ba3b816bbbcf422a45ae7e1507a1daba0f Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 24 Jul 2026 17:42:20 +0200
Subject: [PATCH 22/25] model/Makefile: remove FINCS

should be fixed by update of TSMP2 compilation environment
---
 interface/model/Makefile | 1 -
 1 file changed, 1 deletion(-)

diff --git a/interface/model/Makefile b/interface/model/Makefile
index 2c7f2585..185b906a 100644
--- a/interface/model/Makefile
+++ b/interface/model/Makefile
@@ -1,6 +1,5 @@
 BASEDIR = $(PDAF_DIR)
 include $(BASEDIR)/make.arch/$(PDAF_ARCH).h
-FINCS = -I$(BASEDIR)/include -I$(BASEDIR)/build
 ## includes
 
 FFLAGS = -cpp $(TSMPPDAFIMPORTFLAGS)

From 415d8e2b1eb9fa04b42e21db4494e50b9d590cc6 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 24 Jul 2026 17:42:58 +0200
Subject: [PATCH 23/25] white space

---
 interface/model/eclm/enkf_clm_mod_5.F90 | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index 30020016..d5ce7d3f 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -2721,9 +2721,9 @@ subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, &
     maxlon = MAXVAL(lon(:) + 180)
     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). 
+
+    ! 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)
 
@@ -2759,7 +2759,7 @@ subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, &
           latixy_obs(i) = 1
        end if
     end do
-    
+
     ! deallocate temporary arrays
     !deallocate(longxy)
     !deallocate(latixy)

From 52d03e1100dbfed75a212686d305978bdade4106 Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 24 Jul 2026 17:43:09 +0200
Subject: [PATCH 24/25] wrapper: call print_update_clm only for LAI-DA

---
 interface/model/wrapper_tsmp.c | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c
index aaad257f..a20ff35f 100644
--- a/interface/model/wrapper_tsmp.c
+++ b/interface/model/wrapper_tsmp.c
@@ -208,8 +208,9 @@ void update_tsmp(){
 #endif
 
 #ifdef CLMFIVE
-    if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2 ||
-       (clmupdate_lai == 3 && clmupdate_lai_params == 1)){
+    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" */

From 36ed964e009fd45afa583924ac3cd259dc17472b Mon Sep 17 00:00:00 2001
From: Johannes Keller <jo.keller@fz-juelich.de>
Date: Fri, 24 Jul 2026 17:44:35 +0200
Subject: [PATCH 25/25] style

---
 interface/model/eclm/enkf_clm_mod_5.F90 | 1 -
 1 file changed, 1 deletion(-)

diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90
index d5ce7d3f..3809b60f 100755
--- a/interface/model/eclm/enkf_clm_mod_5.F90
+++ b/interface/model/eclm/enkf_clm_mod_5.F90
@@ -2759,7 +2759,6 @@ subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, &
           latixy_obs(i) = 1
        end if
     end do
-
     ! deallocate temporary arrays
     !deallocate(longxy)
     !deallocate(latixy)