From 7071fe9e314a29370aeee7ede4ddc695074f4f73 Mon Sep 17 00:00:00 2001 From: Lukas Strebel Date: Fri, 21 Nov 2025 10:05:51 +0100 Subject: [PATCH 1/8] 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 Date: Fri, 21 Nov 2025 13:54:46 +0100 Subject: [PATCH 2/8] 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 Date: Fri, 21 Nov 2025 15:26:25 +0100 Subject: [PATCH 3/8] 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 Date: Fri, 21 Nov 2025 15:41:00 +0100 Subject: [PATCH 4/8] 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 Date: Fri, 21 Nov 2025 15:42:14 +0100 Subject: [PATCH 5/8] 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 Date: Fri, 21 Nov 2025 15:47:41 +0100 Subject: [PATCH 6/8] 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 Date: Fri, 21 Nov 2025 15:48:20 +0100 Subject: [PATCH 7/8] 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 Date: Fri, 21 Nov 2025 16:40:37 +0100 Subject: [PATCH 8/8] 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, &