diff --git a/docs/_toc.yml b/docs/_toc.yml index 1ff64785..2dc70d6a 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -22,6 +22,7 @@ parts: - file: users_guide/running_tsmp_pdaf/input_enkfpf - file: users_guide/running_tsmp_pdaf/tsmp_pdaf_omi - file: users_guide/running_tsmp_pdaf/grace_tws_da + - file: users_guide/running_tsmp_pdaf/lai_da - file: users_guide/misc/README title: Misc diff --git a/docs/users_guide/running_tsmp_pdaf/README.md b/docs/users_guide/running_tsmp_pdaf/README.md index 18dfd8f0..fd8753e1 100644 --- a/docs/users_guide/running_tsmp_pdaf/README.md +++ b/docs/users_guide/running_tsmp_pdaf/README.md @@ -14,6 +14,8 @@ COSMO](cos)). Additionally, a control file for the data assimilation ([Observation files](obs)) need to be present in the run directory. Furthermore, some command line options ([Command line options](cmd)) need to be specified when TSMP-PDAF is executed. +For assimilation of LAI observations into eCLM, see +[LAI Data Assimilation](laida). For assimilation of GRACE/GRACE-FO terrestrial water storage observations into eCLM, see [GRACE/TWS Data Assimilation](gracetws). diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md index 78c83b3e..aacf8857 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md +++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md @@ -54,6 +54,9 @@ nprocs = update_swc = update_texture = update_T = +update_lai = +update_lai_params = +update_lai_incr_w = print_swc = print_et = print_inc = @@ -499,6 +502,70 @@ 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. + +(enkfpf:clm:update_lai_incr_w)= +### CLM:update_lai_incr_w ### + +`CLM:update_lai_incr_w`: (real) Weight of the **additive** LAI increment when mapping the gridcell analysis back to patches. Only used if `CLM:update_lai` is non-zero. Default: `1.0`. + +After the EnKF update, the gridcell-mean LAI increment is applied to each patch in two ways and then blended: + +- **Additive** (weight `w`): each patch receives the same increment + $\Delta = \text{LAI}_\text{a} - \text{LAI}_\text{f}$: + $\text{LAI}'(p) = \text{LAI}_\text{f}(p) + \Delta$. + +- **Multiplicative** (weight $1-w$): each patch is scaled by the gridcell analysis-to-forecast ratio: + $\text{LAI}'(p) = \text{LAI}_\text{f}(p)\,\text{LAI}_\text{a} / \text{LAI}_\text{f}$. + +The final patch LAI is +$\text{LAI}'(p) = w\,\text{LAI}'_\text{add}(p) + (1-w)\,\text{LAI}'_\text{mult}(p)$. + +- `1.0`: purely additive. +- `0.0`: purely multiplicative. + ### CLM:print_swc ### `CLM:print_swc`: (integer) If set to `1`, the updated soil moisture @@ -998,6 +1065,9 @@ Default: 0, output turned off. | | `problemname` | \- | | | `nprocs` | 0 | | | `update_swc` | 1 | + | | `update_lai` | 0 | + | | `update_lai_params` | 0 | + | | `update_lai_incr_w` | 1.0 | | | `print_swc` | 0 | | | `print_et` | 0 | | | `print_inc` | 0 | diff --git a/docs/users_guide/running_tsmp_pdaf/lai_da.md b/docs/users_guide/running_tsmp_pdaf/lai_da.md new file mode 100644 index 00000000..d6b624d8 --- /dev/null +++ b/docs/users_guide/running_tsmp_pdaf/lai_da.md @@ -0,0 +1,173 @@ +(laida)= +# LAI Data Assimilation # + +Leaf Area Index (LAI) data assimilation allows observed LAI fields to +be used to update the vegetation state and, optionally, the specific +leaf area parameters of eCLM. + +**eCLM only.** Requires BGC (Biogeochemical Mode) to be active. +Not available for CLM3.5 or CLM5.0 without BGC. + +## Configuration ## + +LAI DA is controlled by three parameters in `enkfpf.par`: + +- [`CLM:update_lai`](enkfpf:clm:update_lai) — selects the state vector + layout and the location of the forward observation operator. +- [`CLM:update_lai_params`](enkfpf:clm:update_lai_params) — enables + joint estimation of specific leaf area parameters alongside LAI. +- [`CLM:update_lai_incr_w`](enkfpf:clm:update_lai_incr_w) — blend between + additive and multiplicative patch LAI increments. + +## Background ## + +eCLM with BGC does not carry LAI as an independent prognostic +variable. Instead, LAI is a diagnostic quantity derived from leaf +carbon (`leafc`) and plant-functional-type (PFT) constants for +specific leaf area: + +$$ +\text{LAI}(p) = +\begin{cases} + \dfrac{s_{\text{top}}(p)\,\bigl(\exp\!\bigl(C_\ell(p)\,s'(p)\bigr) - 1\bigr)}{s'(p)} + & s'(p) > 0 \\[6pt] + s_{\text{top}}(p)\, C_\ell(p) & s'(p) = 0 +\end{cases} +$$ + +where +- $C_\ell$ is `leafc` (leaf carbon, gC m$^{-2}$) +- $s_{\text{top}}$ is `slatop` (specific leaf area at the canopy top, +m$^2$ gC$^{-1}$) +- $s'$ is `dsladlai` (slope of specific leaf +area with respect to LAI, m$^2$ gC$^{-1}$) + +Both $s_{\text{top}}$ and $s'$ are PFT constants read from the CLM +surface dataset. The formula follows Eq. 3 of Thornton and Zimmermann +(2007, J. Clim., 20, 3902–3923). + +Because LAI is not a direct model state, assimilation updates `leafc` +(and, consistently, `leafn = leafc / leafcn`) rather than LAI +itself. + +## State Vector ## + +The choice of [`CLM:update_lai`](enkfpf:clm:update_lai) determines the +state vector layout: + +| `update_lai` | `update_lai_params` | State vector contents | Size | +|:---:|:---:|---|---| +| 1 | 0 | Gridcell-averaged LAI | $N_g$ | +| 1 | 1 | Gridcell-averaged LAI + `slatop` per patch | $N_g + N_p$ | +| 2 | 0 | `leafc`, `slatop`, `dsladlai` per patch | $3\,N_p$ | +| 2 | 2 | `leafc`, `slatop`, `dsladlai` per patch (parameters updated after assimilation) | $3\,N_p$ | + +$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. Per-patch forecast LAI $\text{LAI}_\text{f}(p)$ is stored for the + update phase; + +**PDAF update**: PDAF operates on the gridcell-level LAI values. + +**Update phase** (`update_clm`): + +1. The gridcell analysis LAI $\text{LAI}'_g$ is mapped to patches using + a blend of additive and multiplicative increments controlled by + [`CLM:update_lai_incr_w`](enkfpf:clm:update_lai_incr_w) + (default `1.0` = purely additive). +2. The updated per-patch LAI is inverted to obtain the new `leafc`: + +$$ + C'_\ell(p) = + \begin{cases} + \dfrac{\ln\!\bigl((\text{LAI}'(p)\,s'(p)\,/\,s_{\text{top}}(p)) + 1\bigr)}{s'(p)} + & s'(p) > 0 \\[6pt] + \text{LAI}'(p)\,/\,s_{\text{top}}(p) & s'(p) = 0 + \end{cases} +$$ + +3. `leafc` is clipped to zero (no negative leaf carbon). +4. `leafn` is set consistently: $N'_\ell = C'_\ell\,/\,\text{leafcn}$. + +**Joint parameter estimation** (`update_lai_params=1`): `slatop` for +each patch is appended to the state vector and updated by PDAF. The +updated `slatop` values are written back to the PFT constants before +the LAI inversion in step 2. + +## Option 2 — Patch-Level State Vector with Observation Operator ## + +With `CLM:update_lai=2`, the state vector holds `leafc`, `slatop`, +and `dsladlai` for every patch in three consecutive blocks: + +$$ +\mathbf{x} = \bigl[\underbrace{C_\ell(1),\ldots,C_\ell(N_p)}_{\text{block 1}}, + \underbrace{s_{\text{top}}(1),\ldots,s_{\text{top}}(N_p)}_{\text{block 2}}, + \underbrace{s'(1),\ldots,s'(N_p)}_{\text{block 3}}\bigr] +$$ + +The observation operator (`obs_op_pdaf`, active when +`CLM:update_lai=2`) maps this state to observed gridcell LAI: + +$$ +\mathcal{H}(\mathbf{x})_i = \sum_{p:\,g(p)=g_i} w_p\,\text{LAI}(p) +$$ + +where $g_i$ is the grid cell corresponding to observation $i$. + +After the PDAF update: + +- Patch LAI is updated with the same additive/multiplicative blend as in + option 1 (using weighted gridcell-mean LAI for the ratio and increment). +- `leafc` is obtained by inverting the blended patch LAI and clipped to + zero; `leafn` is updated consistently. +- If `CLM:update_lai_params=2`, `slatop` and `dsladlai` are also + read back from blocks 2 and 3 and written to the PFT constants. + The diagnostic LAI update is then left to eCLM's own phenology + routines. + +## Configuration Examples ## + +**Example 1** — gridcell LAI assimilation, no parameter estimation: + +```ini +[CLM] +update_lai = 1 +update_lai_params = 0 +update_lai_incr_w = 1.0 +``` + +**Example 2** — gridcell LAI assimilation with joint `slatop` estimation: + +```ini +[CLM] +update_lai = 1 +update_lai_params = 1 +``` + +**Example 3** — patch-level state vector with custom observation operator: + +```ini +[CLM] +update_lai = 2 +update_lai_params = 0 +``` + +**Example 4** — patch-level state vector with joint `slatop`/`dsladlai` estimation: + +```ini +[CLM] +update_lai = 2 +update_lai_params = 2 +``` diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index bfd08a43..a73e3cb1 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 @@ -71,6 +71,10 @@ 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 + USE pftconMod, only : pftcon #ifdef CLMFIVE USE clm_instMod, & ONLY : soilstate_inst @@ -100,8 +104,9 @@ 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 +character(len=48) :: fn_hx ! end of hcp #ifndef PARFLOW_STAND_ALONE @@ -128,6 +133,74 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) ! If no special observation operator is compiled, use point observations lpointobs = .true. +! LAI assimilation assuming statevec contains leafc, slatop, dsladlai for each patch +! two component op to get to LAI obs +! 1 calculate tlai from statevec per patch +! state_p contains leafc of varsize then slatop of varsize and then dsadlai of varsize +! 2 average patches according to weight per gridcell +if (clmupdate_lai==2) then + lpointobs = .false. ! disable general obs_op + ! Note: This fails if there ever is an observation pointing to a gridcell with 0 patches + ! this should not happen, but would cause div by zero. + + do i = 1, dim_obs_p ! loop over all observations + avelai = 0.0 ! use avelai as grid cell average collecter + do j = 1, clm_endp-clm_begp+1 ! loop over all patches + write(*,*) 'DEBUG LAI : obs_index_p(i), clm_patch2gc(j)', obs_index_p(i), clm_patch2gc(j) + if (obs_index_p(i)==clm_patch2gc(j)) then + write(*,*) 'DEBUG LAI : state dsladlai', state_p(j+2*clm_varsize) + if (state_p(j+2*clm_varsize)>0.0) then + avelai = avelai + clm_patchwt(j) * ((state_p(j+1*clm_varsize)& + *(exp(state_p(j)*state_p(j+2*clm_varsize)) - 1.0))/state_p(j+2*clm_varsize)) ! formula for tlai from leafc + else + avelai = avelai + clm_patchwt(j) *(state_p(j+1*clm_varsize)*state_p(j)) ! 2nd formula + endif ! dsladlai decider + endif ! this patch is in gridcell of observation + write(*,*) 'DEBUG LAI : average ', avelai + enddo ! all patches for obs checked + m_state_p(i) = avelai ! == lai of gridcell where the observation is + enddo ! all observations +endif ! clmupdate_lai == 2 : m_state_p now contains lai for each gridcell with an observation. + +! LAI assimilation assuming statevec contains leafc, livestemc, deadstemc +! for each vegetated patch (bare-ground PFT excluded). +! H uses 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 if (clmupdate_T==1) then diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h index a734c335..0b01df7f 100755 --- a/interface/model/common/enkf.h +++ b/interface/model/common/enkf.h @@ -91,6 +91,9 @@ GLOBAL int clmupdate_swc; GLOBAL int clmupdate_T; GLOBAL int clmupdate_texture; GLOBAL int clmupdate_tws; +GLOBAL int clmupdate_lai; +GLOBAL int clmupdate_lai_params; +GLOBAL double clmupdate_lai_incr_w; GLOBAL int clmprint_swc; GLOBAL int clmprint_et; GLOBAL int clmprint_inc; diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c index 27b26c46..4b006c50 100755 --- a/interface/model/common/read_enkfpar.c +++ b/interface/model/common/read_enkfpar.c @@ -79,6 +79,9 @@ void read_enkfpar(char *parname) clmupdate_swc = iniparser_getint(pardict,"CLM:update_swc",1); clmupdate_T = iniparser_getint(pardict,"CLM:update_T",0); clmupdate_texture = iniparser_getint(pardict,"CLM:update_texture",0); + clmupdate_lai = iniparser_getint(pardict,"CLM:update_lai",0); + clmupdate_lai_params = iniparser_getint(pardict,"CLM:update_lai_params",0); + clmupdate_lai_incr_w = iniparser_getdouble(pardict,"CLM:update_lai_incr_w",1); clmprint_swc = iniparser_getint(pardict,"CLM:print_swc",0); clmprint_et = iniparser_getint(pardict,"CLM:print_et",0); clmprint_inc = iniparser_getint(pardict,"CLM:print_inc",0); diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 68bd8c54..89a5dd00 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -52,10 +52,17 @@ module enkf_clm_mod ! LST in LST assimilation (clmupdate_T) real(r8),allocatable :: clm_paramarr(:) !hcp CLM parameter vector (f.e. LAI) integer, allocatable :: state_clm2pdaf_p(:,:) !Index of column in hydraulic active state vector (nlevsoi,endc-begc+1) + real(r8),allocatable :: clm_patch2gc(:) ! index array to simplify patch to gridcell averaging + real(r8),allocatable :: clm_patchwt(:) ! weight array to simplify patch to gridcell averaging + integer,allocatable :: clm_lai_patch_idx(:) ! mode-3: CLM patch index per state slot + integer,allocatable :: clm_lai_patch_itype(:) ! mode-3: PFT type per state slot integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc + integer(c_int),bind(C,name="clmupdate_lai") :: clmupdate_lai + integer(c_int),bind(C,name="clmupdate_lai_params") :: clmupdate_lai_params + real(c_double),bind(C,name="clmupdate_lai_incr_w") :: clmupdate_lai_incr_w ! Yorck integer(c_int),bind(C,name="clmupdate_tws") :: clmupdate_tws @@ -125,12 +132,18 @@ module enkf_clm_mod ! (currently not used for eclm) logical :: newgridcell !only eclm + ! Arrays for LAI Assimilation + real(r8),allocatable :: tlai(:) + real(r8),allocatable :: tlai_orig(:) + real(r8),allocatable :: lai_wtsum2_gc(:) + contains #if defined CLMSA subroutine define_clm_statevec(mype) use decompMod , only : get_proc_bounds use clm_varpar , only : nlevsoi + use PatchType , only : patch use clm_varcon, only: set_averaging_to_zero use PDAF_interfaces_module, only: PDAF_reset_dim_p @@ -142,6 +155,7 @@ subroutine define_clm_statevec(mype) integer :: begc, endc ! per-proc beginning and ending column indices integer :: begl, endl ! per-proc beginning and ending landunit indices integer :: begg, endg ! per-proc gridcell ending gridcell indices + integer :: i, cc call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp) @@ -191,6 +205,80 @@ subroutine define_clm_statevec(mype) clm_statevecsize = clm_statevecsize + 3*((endg-begg+1)*nlevsoi) end if + ! LAI assimilation + if(clmupdate_lai==1) then + ! Allocate array to store total leaf area index per patch + if (allocated(tlai)) deallocate(tlai) + allocate(tlai(clm_begp:clm_endp)) + if (allocated(tlai_orig)) deallocate(tlai_orig) + allocate(tlai_orig(clm_begp:clm_endp)) + if (allocated(lai_wtsum2_gc)) deallocate(lai_wtsum2_gc) + allocate(lai_wtsum2_gc(clm_begg:clm_endg)) + + ! #patches values per grid-cell ! clm_varsize not accurate due to variable #patches/grid cell + ! set sizes + clm_varsize = 1 ! each grid cell gets 1 value of LAI + ! Currently no combination of SWC and LAI DA + + clm_statevecsize = (clm_endg - clm_begg + 1) ! statevector is the size of the num of gridcells + !(clm_endp-clm_begp+1) ! So like this if lai is set it takes priority + + if (clmupdate_lai_params==1) then ! add parameters to the state vector. + ! clm_varsize = (clm_endg - clm_begg+1) ! Currently no combination of SWC and LAI DA + ! clm_statevecsize = (clm_endg - clm_begg+1)*1 ! 1 parameters + clm_varsize = (clm_endp-clm_begp+1)+1 ! Currently no combination of SWC and LAI DA + clm_statevecsize = (clm_endp-clm_begp+1)+1 ! 1 parameters + ! So like this if lai is set it takes priority + end if + + ! if (clmupdate_lai_params==3) then + ! clm_varsize = (clm_endp-clm_begp+1) + ! clm_statevecsize = 3*clm_varsize + ! end if + + end if + + ! for obs_op version allocate additional helper index array + if(clmupdate_lai==2) then + clm_varsize = (clm_endp-clm_begp+1) ! equals number of patches + clm_statevecsize = 3*clm_varsize ! 3 var/params per patch + + if (allocated(clm_patch2gc)) deallocate(clm_patch2gc) + allocate(clm_patch2gc(clm_varsize)) + if (allocated(clm_patchwt)) deallocate(clm_patchwt) + allocate(clm_patchwt(clm_varsize)) + + endif + + ! LAI assimilation mode 3: leafc, livestemc, deadstemc per vegetated patch + if(clmupdate_lai==3) then + clm_varsize = 0 + do i=clm_begp,clm_endp + if (.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" @@ -240,17 +328,16 @@ subroutine define_clm_statevec(mype) #endif IF (allocated(clm_statevec)) deallocate(clm_statevec) - if ((clmupdate_swc/=0) .or. (clmupdate_T/=0) .or. (clmupdate_texture/=0) .or. (clmupdate_tws/=0)) then + if ((clmupdate_swc/=0) .or. (clmupdate_T/=0) .or. (clmupdate_texture/=0) .or. (clmupdate_tws/=0)) .or. (clmupdate_lai/=0)) then !hcp added condition allocate(clm_statevec(clm_statevecsize)) end if - ! Allocate statevector-duplicate for saving original column mean - ! values used in computing increments during updating the state - ! vector in column-mean-mode. + ! Allocate statevector-duplicate for saving the forecast state 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) .or. clmupdate_tws/=0 ) then + if ( (clmupdate_swc/=0 .and. clmstatevec_colmean/=0) .or. clmupdate_tws/=0 .or. (clmupdate_lai/=0)) then allocate(clm_statevec_orig(clm_statevecsize)) end if @@ -701,16 +788,28 @@ subroutine cleanup_clm_statevec() IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig) + IF (allocated(clm_patch2gc)) deallocate(clm_patch2gc) + IF (allocated(clm_patchwt)) deallocate(clm_patchwt) + IF (allocated(clm_lai_patch_idx)) deallocate(clm_lai_patch_idx) + IF (allocated(clm_lai_patch_itype)) deallocate(clm_lai_patch_itype) + IF (allocated(lai_wtsum2_gc)) deallocate(lai_wtsum2_gc) + IF (allocated(tlai)) deallocate(tlai) + IF (allocated(tlai_orig)) deallocate(tlai_orig) end subroutine cleanup_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : bgc_vegetation_inst, canopystate_inst use clm_varpar , only : nlevsoi ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col + use PatchType , only : patch + use pftconMod , only : pftcon use shr_kind_mod, only: r8 => shr_kind_r8 + use PhotosynthesisMod, only : params_inst + implicit none integer,intent(in) :: tstartcycle integer,intent(in) :: mype @@ -718,10 +817,14 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) - integer :: i,j,jj,g,c,cc,offset + real(r8), pointer :: leafc(:) + real(r8), pointer :: livestemc(:) + real(r8), pointer :: deadstemc(:) + integer :: i,j,jj,g,c,cc,offset,igc integer :: n_c character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output + character (len = 48) :: fn3 !TSMP-PDAF: function name for lai_leafc_* integrate output cc = 0 offset = 0 @@ -730,6 +833,8 @@ subroutine set_clm_statevec(tstartcycle, mype) psand => soilstate_inst%cellsand_col pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN @@ -754,6 +859,126 @@ subroutine set_clm_statevec(tstartcycle, mype) error stop "Not implemented clmupdate_swc.eq.2" end if + !> @author Lukas Strebel, Haojin Zhao + !> @date 04.2026 + ! LAI assimilation + ! Case 1: LAI all patches + ! Case 1: Transformations in Set and Update functions + if(clmupdate_lai==1) then + clm_statevec(:) = 0._r8 ! reset statevector to 0 because we add to it for the average + lai_wtsum2_gc(:) = 0._r8 + do i=clm_begp,clm_endp + lai_wtsum2_gc(patch%gridcell(i)) = lai_wtsum2_gc(patch%gridcell(i)) & + + patch%wtgcell(i)**2 + end do + do i=clm_begp,clm_endp ! iterate through patches + ! update the leaf area index based on leafC and SLA + ! Eq 3 from Thornton and Zimmerman, 2007, J Clim, 20, 3902-3923. + ! (slatop(ivt(p))*(exp(leafc(p)*dsladlai(ivt(p))) - 1._r8))/dsladlai(ivt(p)) + if (pftcon%dsladlai(patch%itype(i)) > 0._r8) then + tlai(i) = (pftcon%slatop(patch%itype(i))*(exp(leafc(i)*pftcon%dsladlai(patch%itype(i))) - 1._r8)) & + /pftcon%dsladlai(patch%itype(i)) + else + tlai(i) = pftcon%slatop(patch%itype(i)) * leafc(i) + end if + tlai(i) = max(0._r8, tlai(i)) ! don't allow negative LAI + tlai_orig(i) = tlai(i) + ! Map CLM gridcell index (begg:endg) to PE-local state index (1:clm_statevecsize) + igc = patch%gridcell(i) - clm_begg + 1 + if (igc >= 1 .and. igc <= clm_statevecsize) then + clm_statevec(igc) = clm_statevec(igc) + patch%wtgcell(i)*tlai(i) + write(*,'(A,I6,A,I6,A,F10.4,A,F10.4,A,F10.4,A,F10.4,A,F10.4,A,F10.4,A,F10.4,A,F10.4)') & + "LAI1 DEBUG setvec ", i, ";", igc, ";", patch%wtgcell(i), ";", & + tlai(i), ";", patch%wtgcell(i)*tlai(i), ";", clm_statevec(igc), ";", & + canopystate_inst%tlai_patch(i), ";", leafc(i), ";", & + pftcon%dsladlai(patch%itype(i)), ";", pftcon%slatop(patch%itype(i)) + else + write(*,*) 'DEBUG LAI : igc out of bounds in LAI mapping: ', igc, & + ' valid range: 1-', clm_statevecsize, ' i=', i + end if + end do + if (allocated(clm_statevec_orig)) then + do igc = 1, min(clm_endg - clm_begg + 1, clm_statevecsize) + clm_statevec_orig(igc) = clm_statevec(igc) + end do + end if + +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN + IF(clmupdate_lai/=0) THEN + ! Same pattern as swcstate_*.integrate.*: one file, consecutive es22.15 blocks. + ! Order: leafc(patch) tlai_calc(patch) tlai_clm(patch) L_gc(gridcell) lai_wtsum2_gc(gridcell) + WRITE(fn3, "(a,i5.5,a,i5.5,a)") "lai_leafc_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn3, action="write", status="replace") + WRITE(71,"(a)") "# blocks: leafc(patch) tlai_calc(patch) tlai_clm(patch) L_gc(gridcell) lai_wtsum2_sq(gridcell)" + WRITE(71,"(es22.15)") leafc(clm_begp:clm_endp) + WRITE(71,"(es22.15)") tlai(clm_begp:clm_endp) + WRITE(71,"(es22.15)") canopystate_inst%tlai_patch(clm_begp:clm_endp) + WRITE(71,"(es22.15)") clm_statevec(1:clm_statevecsize) + WRITE(71,"(es22.15)") lai_wtsum2_gc(clm_begg:clm_endg) + CLOSE(71) + END IF + END IF +#endif + + endif + + ! Case 2: statevec with leafc, slatop, dsladlai for transform in obs_op + if (clmupdate_lai==2) then + cc = 1 + do i=clm_begp,clm_endp + clm_statevec(cc) = leafc(i) + clm_statevec(cc + 1*clm_varsize) = pftcon%slatop(patch%itype(i)) + clm_statevec(cc + 2*clm_varsize) = pftcon%dsladlai(patch%itype(i)) + clm_patch2gc(cc) = patch%gridcell(i) + clm_patchwt(cc) = patch%wtgcell(i) + cc = cc + 1 + enddo + write(*,*) 'DEBUG LAI : statevec ', clm_statevec(:) + write(*,*) 'DEBUG LAI : patches ', clm_patch2gc(:), clm_patchwt(:) + endif + + ! Case 3: leafc, livestemc, deadstemc per vegetated patch; H(x) in obs_op + if (clmupdate_lai==3) then + livestemc => bgc_vegetation_inst%cnveg_carbonstate_inst%livestemc_patch + deadstemc => bgc_vegetation_inst%cnveg_carbonstate_inst%deadstemc_patch + cc = 1 + do j = 1, clm_varsize + i = clm_lai_patch_idx(j) + clm_statevec(cc) = leafc(i) + clm_statevec(cc + clm_varsize) = livestemc(i) + clm_statevec(cc + 2*clm_varsize) = deadstemc(i) + cc = cc + 1 + end do + if (allocated(clm_statevec_orig)) then + clm_statevec_orig(:) = clm_statevec(:) + end if + +#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 + clm_statevec(cc+1*clm_varsize+offset) = pftcon%slatop(patch%itype(i)) + cc = cc + 1 + end do + endif + + ! if (clmupdate_lai_params==3) then + ! cc = 1 + ! do i=clm_begp,clm_endp + ! clm_statevec(cc+1*clm_varsize+offset) = pftcon%slatop(patch%itype(i)) + ! clm_statevec(cc+2*clm_varsize+offset) = params_inst%kmax(patch%itype(i),1) + ! clm_statevec(cc+3*clm_varsize+offset) = params_inst%kmax(patch%itype(i),2) + ! cc = cc + 1 + ! end do + ! endif + !hcp LAI if(clmupdate_T==1) then error stop "Not implemented: clmupdate_T.eq.1" @@ -1144,10 +1369,123 @@ end subroutine set_clm_statevec_tws + subroutine lai_patch_tlai_blend(tlai_fc, lai_gc_fc, lai_gc_an, tlai_out) + use shr_kind_mod , only : r8 => shr_kind_r8 + + implicit none + + real(r8), intent(in) :: tlai_fc + real(r8), intent(in) :: lai_gc_fc + real(r8), intent(in) :: lai_gc_an + real(r8), intent(out) :: tlai_out + + real(r8) :: incr_lai_gc + real(r8) :: tlai_add + real(r8) :: tlai_mult + real(r8) :: w_add + + incr_lai_gc = lai_gc_an - lai_gc_fc + tlai_add = tlai_fc + incr_lai_gc + if (abs(lai_gc_fc) > 1.0e-12_r8) then + tlai_mult = tlai_fc * (lai_gc_an / lai_gc_fc) + else + tlai_mult = tlai_add + end if + w_add = real(clmupdate_lai_incr_w, r8) + tlai_out = w_add * tlai_add + (1._r8 - w_add) * tlai_mult + tlai_out = max(0._r8, tlai_out) + + end subroutine lai_patch_tlai_blend + + + subroutine write_lai_da_mode3_debug(phase, cycle, mype) + use clm_instMod, only : bgc_vegetation_inst, canopystate_inst + + implicit none + + character(len=*), intent(in) :: phase + integer, intent(in) :: cycle + integer, intent(in) :: mype + + real(r8), pointer :: leafc(:) + real(r8), pointer :: livestemc(:) + real(r8), pointer :: deadstemc(:) + real(r8), pointer :: leafn(:) + real(r8), pointer :: livestemn(:) + real(r8), pointer :: deadstemn(:) + + character(len=48) :: fn + integer :: j, i + + IF(clmt_printensemble /= cycle .AND. clmt_printensemble /= -1) RETURN + + leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch + livestemc => bgc_vegetation_inst%cnveg_carbonstate_inst%livestemc_patch + deadstemc => bgc_vegetation_inst%cnveg_carbonstate_inst%deadstemc_patch + leafn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%leafn_patch + livestemn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%livestemn_patch + deadstemn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%deadstemn_patch + + WRITE(fn, "(a,i5.5,a,a,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 use clm_instMod, only : waterstate_inst + use clm_instMod, only : bgc_vegetation_inst, canopystate_inst + use pftconMod , only : pftcon + use PatchType , only : patch + use PhotosynthesisMod, only : params_inst implicit none @@ -1156,16 +1494,35 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) + real(r8), pointer :: leafc(:) ! leaf carbon of patch + real(r8), pointer :: leafn(:) ! leaf nitrogen of patch + real(r8), pointer :: livestemc(:) + real(r8), pointer :: deadstemc(:) + real(r8), pointer :: livestemn(:) + real(r8), pointer :: deadstemn(:) + + real(r8) :: lai_gc_an + real(r8) :: lai_gc_fc integer :: i + integer :: j + integer :: itype + integer :: cc + integer :: offset + integer :: igc character (len = 31) :: fn !TSMP-PDAF: function name for state vector output - character (len = 32) :: fn5 !TSMP-PDAF: function name for state vector outpu - character (len = 32) :: fn6 !TSMP-PDAF: function name for state vector outpu + character (len = 31) :: fn2 !TSMP-PDAF: function name for state vector output + character (len = 32) :: fn5 !TSMP-PDAF: function name for state vector output + character (len = 32) :: fn6 !TSMP-PDAF: function name for state vector output + character (len = 32) :: fn7 !TSMP-PDAF: function name for state vector output logical :: swc_zero_before_update swc_zero_before_update = .false. + cc = 1 + offset = 0 + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN ! TSMP-PDAF: For debug runs, output the state vector in files @@ -1181,6 +1538,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col + leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch + leafn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%leafn_patch + livestemc => bgc_vegetation_inst%cnveg_carbonstate_inst%livestemc_patch + deadstemc => bgc_vegetation_inst%cnveg_carbonstate_inst%deadstemc_patch + livestemn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%livestemn_patch + deadstemn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%deadstemn_patch + + + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN @@ -1198,6 +1564,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") CLOSE(71) END IF + IF(clmupdate_lai/=0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn5, "(a,i5.5,a,i5.5,a)") "leafc", mype, ".bef_up.", tstartcycle, ".txt" + OPEN(unit=71, file=fn5, action="write") + WRITE (71,"(es22.15)") leafc(:) + CLOSE(71) + END IF + END IF #endif @@ -1233,6 +1607,159 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") + ! LAI assimilation: + ! Case 1: gridcell-mean LAI state; map analysis back to patches via + ! additive/multiplicative blend (CLM:update_lai_incr_w). + if(clmupdate_lai==1) then + do i=clm_begp,clm_endp + igc = patch%gridcell(i) - clm_begg + 1 + if (igc >= 1 .and. igc <= clm_statevecsize .and. allocated(clm_statevec_orig)) then + lai_gc_fc = clm_statevec_orig(igc) + lai_gc_an = clm_statevec(igc) + call lai_patch_tlai_blend(tlai_orig(i), lai_gc_fc, lai_gc_an, tlai(i)) + else + tlai(i) = 0._r8 + endif + ! update leafc based on Eq 3 from Thornton and Zimmerman, 2007, J Clim, 20, 3902-3923 + ! reformulate the equation to solve for leafc + ! => leafc(p) = log(((tlai(p) * dsladlai(ivt(p)))/slatop(ivt(p))) + 1._r8) / dsladlai(ivt(p)) + if (tlai(i) > 0._r8) then ! invalid log protection + if (pftcon%dsladlai(patch%itype(i)) > 0._r8) then + if (pftcon%slatop(patch%itype(i)) > 0._r8) then + leafc(i) = log(((tlai(i) * pftcon%dsladlai(patch%itype(i))) & + / pftcon%slatop(patch%itype(i))) + 1.0_r8) & + / pftcon%dsladlai(patch%itype(i)) + else + leafc(i) = 0._r8 + endif + else + ! Forward: tlai = slatop * leafc (dsladlai == 0) => leafc = tlai / slatop + if (pftcon%slatop(patch%itype(i)) > 0._r8) then + leafc(i) = tlai(i) / pftcon%slatop(patch%itype(i)) + else + leafc(i) = 0._r8 + end if + endif + else + leafc(i) = 0._r8 + endif + leafc(i) = max(0._r8, leafc(i)) ! Don't allow negative leaf carbon + if (pftcon%leafcn(patch%itype(i)) > 0._r8) then + leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i)) + else + leafn(i) = 0._r8 + endif + ! DEBUG ONLY + ! print *, "LAI1 DEBUG leafc,tlai(wt) AFTER: ", i, ";", leafc(i), ";", tlai(i) + end do + +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN + ! TSMP-PDAF: lai & leafc (legacy column vector) + WRITE(fn7, "(a,i5.5,a,i5.5,a)") "lai_leafc_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn7, action="write", status="replace") + WRITE(71,"(a)") "# patch leafc_aft tlai_patch_aft" + DO i = clm_begp, clm_endp + WRITE(71,"(i8,2es22.15)") i, leafc(i), tlai(i) + END DO + CLOSE(71) + END IF +#endif + endif + + ! Case 2: statevec contains leafc, slatop, dsladlai + ! update of diagnostic tlai left to be done by clm5 itself. + if(clmupdate_lai==2) then + cc = 1 + do i=clm_begp,clm_endp + leafc(i) = clm_statevec(cc) + leafc(i) = max(0._r8, leafc(i)) + if (pftcon%leafcn(patch%itype(i)) > 0._r8) then + leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i)) + else + leafn(i) = 0._r8 + end if + if(clmupdate_lai_params==2) then + pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize) + pftcon%dsladlai(patch%itype(i)) = clm_statevec(cc+2*clm_varsize) + endif + + cc = cc + 1 + enddo + +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn7, "(a,i5.5,a,i5.5,a)") "leafc_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn7, action="write") + WRITE (71,"(es22.15)") leafc(:) + CLOSE(71) + END IF +#endif + + endif + + ! Case 3: leafc, livestemc, deadstemc per vegetated patch + if(clmupdate_lai==3) then + cc = 1 + do j = 1, clm_varsize + i = clm_lai_patch_idx(j) + itype = clm_lai_patch_itype(j) + leafc(i) = max(0._r8, clm_statevec(cc)) + livestemc(i) = max(0._r8, clm_statevec(cc + clm_varsize)) + deadstemc(i) = max(0._r8, clm_statevec(cc + 2*clm_varsize)) + if (pftcon%leafcn(itype) > 0._r8) then + leafn(i) = leafc(i) / pftcon%leafcn(itype) + else + leafn(i) = 0._r8 + end if + if (pftcon%livewdcn(itype) > 0._r8) then + livestemn(i) = livestemc(i) / pftcon%livewdcn(itype) + else + livestemn(i) = 0._r8 + end if + if (pftcon%deadwdcn(itype) > 0._r8) then + deadstemn(i) = deadstemc(i) / pftcon%deadwdcn(itype) + else + deadstemn(i) = 0._r8 + end if + cc = cc + 1 + end do + +#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 + pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize+offset) + ! DEBUG ONLY + print *, "LST DEBUG PARAM,i: ", pftcon%slatop(patch%itype(i)), ",", i + cc = cc + 1 + end do + endif + + ! if (clmupdate_lai_params==3) then + ! cc = 1 + ! do i=clm_begp,clm_endp + ! pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize+offset) + ! params_inst%kmax(patch%itype(i),1) = clm_statevec(cc+2*clm_varsize+offset) + ! params_inst%kmax(patch%itype(i),2) = clm_statevec(cc+3*clm_varsize+offset) + + ! ! DEBUG ONLY + ! print *, "LST DEBUG PARAM slatop,i: ", pftcon%slatop(patch%itype(i)), ",", i + ! print *, "LST DEBUG PARAM kmax sun,i: ", params_inst%kmax(patch%itype(i),1), ",", i + ! print *, "LST DEBUG PARAM kmax shade,i: ", params_inst%kmax(patch%itype(i),2), ",", i + + ! cc = cc + 1 + ! end do + + ! end if + end subroutine update_clm @@ -2102,6 +2629,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 @@ -2160,6 +2689,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)) @@ -2167,7 +2701,16 @@ subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, & allocate(latixy_obs(dim_obs)) do i = 1, dim_obs - if(((lon_clmobs(i) + 180) - minlon) /= 0 .and. & + if (dlon_span <= tol_degen .and. dlat_span <= tol_degen) then + longxy_obs(i) = 1 + latixy_obs(i) = 1 + else if (dlon_span <= tol_degen) then + longxy_obs(i) = 1 + latixy_obs(i) = ceiling(((lat_clmobs(i) + 90) - minlat) * nj / (maxlat - minlat)) + else if (dlat_span <= tol_degen) then + longxy_obs(i) = ceiling(((lon_clmobs(i) + 180) - minlon) * ni / (maxlon - minlon)) + latixy_obs(i) = 1 + else if(((lon_clmobs(i) + 180) - minlon) /= 0 .and. & ((lat_clmobs(i) + 90) - minlat) /= 0) then longxy_obs(i) = ceiling(((lon_clmobs(i) + 180) - minlon) * ni / (maxlon - minlon)) !+ 1 latixy_obs(i) = ceiling(((lat_clmobs(i) + 90) - minlat) * nj / (maxlat - minlat)) !+ 1 @@ -2184,6 +2727,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) diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c index 87c2d093..16892c79 100644 --- a/interface/model/wrapper_tsmp.c +++ b/interface/model/wrapper_tsmp.c @@ -198,7 +198,7 @@ void integrate_tsmp() { void update_tsmp(){ #if defined CLMSA - if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_tws != 0))){ + if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_tws != 0) || (clmupdate_lai != 0))){ update_clm(&tstartcycle, &mype_world); #ifndef CLMFIVE