Tsmp pdaf patched lstda cap#49
Tsmp pdaf patched lstda cap#49visakhraja wants to merge 61 commits intoHPSCTerrSys:tsmp-pdaf-patchedfrom
Conversation
from TSMP repo, branch `dev-lstda`, commit 143ea800
``` eclm/enkf_clm_mod_5.F90(413): error #6404: This name does not have a type, and must have an explicit type. [T_SKIN] ```
this was wrongly placed outside the condition checking `newgridcell` leading to `obs_index_p` not being set at all, if not at the very first patch.
`clmupdate_T.eq.3`
Now the update is written, according to the input `CLM:t_printensemble` from `enkfpf.par`.
Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
if the absolute value of the increment is larger than clmT_max_increment, use the clmT_max_increment as increment with the sign of the actual increment.
Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
|
Closing this PR, duplicate of #47 |
There was a problem hiding this comment.
Pull request overview
This PR extends TSMP-PDAF’s CLM Land Surface Temperature (LST) data assimilation support (primarily for CLM5) by implementing new clmupdate_T modes and associated state-vector mappings, update logic, observation indexing, and new runtime parameters for masking/limiting temperature increments.
Changes:
- Implement
define_clm_statevec_T,set_clm_statevec_T, andupdate_clm_Tinenkf_clm_mod_5.F90forclmupdate_Tmodes 1–3. - Add new CLM DA runtime parameters (
T_mask_snow,T_mask_T,increment_type,T_max_increment) wired through C globals and Fortranbind(C). - Update observation indexing and observation operator behavior for
clmupdate_T==2/3, and document new parameters in the user guide.
Reviewed changes
Copilot reviewed 9 out of 9 changed files in this pull request and generated 15 comments.
Show a summary per file
| File | Description |
|---|---|
| interface/model/eclm/enkf_clm_mod_5.F90 | Adds CLM5 LST-DA state-vector definitions, state extraction, and model update logic for clmupdate_T 1–3, plus domain mapping changes. |
| interface/model/common/read_enkfpar.c | Reads new CLM temperature DA parameters from the ini file. |
| interface/model/common/enkf.h | Declares new global parameters backing the C/Fortran bindings. |
| interface/model/clm3_5/enkf_clm_mod.F90 | Adds debug output changes and a new allocatable (but introduces a compile issue). |
| interface/framework/obs_op_pdaf.F90 | Extends CLM observation operator for new clmupdate_T modes. |
| interface/framework/localize_covar_pdaf.F90 | Removes an unused clmupdate_T import. |
| interface/framework/init_dim_obs_pdaf.F90 | Updates CLM observation-to-state indexing for clmupdate_T (patch-based path for CLM5). |
| interface/framework/init_dim_obs_f_pdaf.F90 | Same as above for the “f” variant. |
| docs/users_guide/running_tsmp_pdaf/input_enkfpf.md | Documents new CLM temperature DA parameters and extends update_T description (partially). |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| if(clmupdate_T==1) then | ||
|
|
||
| IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) | ||
| allocate(state_clm2pdaf_p(begp:endp,1)) | ||
|
|
||
| do p=clm_begp,clm_endp | ||
| state_clm2pdaf_p(p,1) = (p - clm_begp + 1) | ||
| end do | ||
|
|
||
| clm_varsize = endp-begp+1 | ||
| clm_paramsize = endp-begp+1 !LAI | ||
| clm_statevecsize = (endp-begp+1)*2 !TG, then TV | ||
|
|
There was a problem hiding this comment.
In the clmupdate_T==1 branch, the state vector is sized/indexed per-patch (clm_varsize=endp-begp+1), but elsewhere (e.g., domain definition and obs indexing) LST-DA is treated as gridcell-based. This mismatch can lead to out-of-bounds or only updating a subset of patches. Please either (a) make clmupdate_T==1 consistently gridcell-mean (size/mapping based on gridcells) or (b) treat it consistently as patch-based across init_n_domains_clm/init_dim_obs_*/local analysis routines.
| allocate(col_updated(clm_begc:clm_endc)) | ||
| col_updated = .false. |
There was a problem hiding this comment.
col_updated is allocated unconditionally, but only deallocated in the clmupdate_T==2 branch. For clmupdate_T==1 or ==3 this leaks memory each DA cycle. Allocate col_updated only inside the clmupdate_T==2 block or ensure it is deallocated on all control-flow paths (including ==1 and ==3).
| allocate(col_updated(clm_begc:clm_endc)) | |
| col_updated = .false. | |
| ! Allocate col_updated only when needed (clmupdate_T == 2) to avoid leaks | |
| if (clmupdate_T == 2) then | |
| allocate(col_updated(clm_begc:clm_endc)) | |
| col_updated = .false. | |
| endif |
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = 0.0 | ||
| n_p = 0 | ||
| do p=clm_begp,clm_endp | ||
| if(patch%gridcell(p)==g) then | ||
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) + t_veg(p) | ||
| n_p = n_p + 1 | ||
| end if | ||
| end do | ||
| if(n_p > 0) then | ||
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) / real(n_p, r8) | ||
| else | ||
| write(*,*) "ERROR: Gridcell g=", g, " has no patches for TVEG averaging" | ||
| error stop "Gridcell without patches in set_clm_statevec_T (TVEG)" |
There was a problem hiding this comment.
For clmupdate_T==3, the TVEG gridcell mean is computed including patches where t_veg can be the CLM spval (~1e36) (you already guard against this in the clmupdate_T==2 branch). This will contaminate the mean state vector and can generate extreme increments. Filter out non-vegetated/bare patches the same way as in the clmupdate_T==2 TVEG averaging (e.g., ignore t_veg(p) above a large threshold) and define a safe fallback if no valid vegetated patches exist.
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = 0.0 | |
| n_p = 0 | |
| do p=clm_begp,clm_endp | |
| if(patch%gridcell(p)==g) then | |
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) + t_veg(p) | |
| n_p = n_p + 1 | |
| end if | |
| end do | |
| if(n_p > 0) then | |
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) / real(n_p, r8) | |
| else | |
| write(*,*) "ERROR: Gridcell g=", g, " has no patches for TVEG averaging" | |
| error stop "Gridcell without patches in set_clm_statevec_T (TVEG)" | |
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = 0.0_r8 | |
| n_p = 0 | |
| do p=clm_begp,clm_endp | |
| if(patch%gridcell(p)==g) then | |
| ! Ignore non-vegetated/bare patches where t_veg holds a large CLM spval (~1e36) | |
| if (abs(t_veg(p)) < 1.0e30_r8) then | |
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) + t_veg(p) | |
| n_p = n_p + 1 | |
| end if | |
| end if | |
| end do | |
| if(n_p > 0) then | |
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) / real(n_p, r8) | |
| else | |
| ! No valid vegetated patches found; use a safe fallback instead of aborting | |
| clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = 0.0_r8 |
| real(r8),allocatable :: clm_statevec(:) | ||
| real(r8),allocatable :: clm_statevec_orig(:) | ||
| integer,allocatable :: state_pdaf2clm_c_p(:) | ||
| integer,allocatable :: state_pdaf2clm_p_p(:) | ||
| integer,allocatable :: state_pdaf2clm_j_p(:) | ||
| integer,allocatable :: state_loc2clm_c_p(:) | ||
| integer,allocatable :: state_loc2clm_p_p(:) | ||
| ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model | ||
| ! LST in LST assimilation (clmupdate_T) | ||
| real(r8),allocatable :: clm_paramarr(:) !hcp CLM parameter vector (f.e. LAI) |
There was a problem hiding this comment.
The newly introduced allocatables state_pdaf2clm_p_p and state_loc2clm_p_p (and clm_paramarr when clmupdate_T==1) are not deallocated in cleanup_clm_statevec, which can leave allocations around across reinitialization or test runs. Consider extending cleanup_clm_statevec to deallocate these arrays as well (and any other new allocatables introduced for LST-DA).
| t_update = t_skin(p) + sign(clmT_max_increment, increment_factor) | ||
| incr_warn_count_skin = incr_warn_count_skin + 1 | ||
| if(incr_warn_switch_skin) then | ||
| print *, "WARNING: t_skin increment is larger then T_max_increment at p=", p |
There was a problem hiding this comment.
Spelling in warning message: "larger then" should be "larger than" to make logs easier to read/search.
| print *, "WARNING: t_skin increment is larger then T_max_increment at p=", p | |
| print *, "WARNING: t_skin increment is larger than T_max_increment at p=", p |
| ! Set index in state vector, LST will be computed | ||
| ! for first patch appearing here | ||
| obs_index_p(cnt) = state_clm2pdaf_p(p,1) | ||
|
|
||
| !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) | ||
| obs_p(cnt) = clm_obs(i) | ||
| if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) | ||
| cnt = cnt + 1 | ||
|
|
There was a problem hiding this comment.
In the clmupdate_T observation-indexing branch, the code only selects the first patch per gridcell (newgridcell guard). If clmupdate_T==1 is patch-based, this will under-observe/under-update the state. Please align obs_index_p generation with the actual state-vector layout for each clmupdate_T mode (gridcell-mean vs patch-level).
| below this threshold temperature, no T is not updated. | ||
|
|
||
| Only takes effect if `CLM:update_T``is switched on. | ||
|
|
||
| Default setting is `0.`: Masking updates below freezing temperaturs. |
There was a problem hiding this comment.
There are several formatting/typo issues in the new CLM:T_mask_* docs: extra backticks in "`CLM:update_T``is switched on", grammar in "no T is not updated", and spelling in "freezing temperaturs". Please clean these up to avoid confusing rendered documentation.
| below this threshold temperature, no T is not updated. | |
| Only takes effect if `CLM:update_T``is switched on. | |
| Default setting is `0.`: Masking updates below freezing temperaturs. | |
| below this threshold temperature, no T is updated. | |
| Only takes effect if `CLM:update_T` is switched on. | |
| Default setting is `0.`: Masking updates below freezing temperatures. |
| ! Set index in state vector, LST will be computed | ||
| ! for first patch appearing here | ||
| obs_index_p(cnt) = state_clm2pdaf_p(p,1) | ||
|
|
||
| !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) | ||
| obs_p(cnt) = clm_obs(i) | ||
| if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) | ||
| cnt = cnt + 1 |
There was a problem hiding this comment.
For clmupdate_T==1 the state vector appears to be patch-based, but this branch picks only the first patch per gridcell and assigns obs_index_p for that single patch. This likely leaves many patches never observed/updated. Please clarify the intended domain (gridcell vs patch) for update_T==1 and adjust observation indexing accordingly (e.g., use a gridcell-mean state layout, or generate indices for all relevant patches).
| - 2: Update of soil and vegetation temperature (prognostic variables) | ||
| based on comparison of satellite LST with skin temperature | ||
| (diagnostic variable) | ||
|
|
There was a problem hiding this comment.
Docs for CLM:update_T describe modes 0/1/2, but the code now implements clmupdate_T==3 as well. Please update this section to document mode 3 (and any differences between 2 and 3) so users don’t accidentally enable an undocumented mode.
| - 3: Alternative (experimental) implementation of the soil and | |
| vegetation temperature update. For details on the algorithm and | |
| differences to mode 2, please refer to the CLM update routine | |
| in the TSMP-PDAF source code. |
| | | `t_printensemble` | -2 | | ||
| | | `watmin_switch` | 0 | | ||
| | | `T_mask_snow` | 0 | | ||
| | | `T_mask_T` | 0.0 | |
There was a problem hiding this comment.
The defaults table at the end was updated for T_mask_snow/T_mask_T, but it doesn’t list the new increment_type and T_max_increment keys that were added above. Please extend the table so it stays complete/consistent with the documented parameters.
| | | `T_mask_T` | 0.0 | | |
| | | `T_mask_T` | 0.0 | | |
| | | `increment_type` | 0 | | |
| | | `T_max_increment` | 0.0 | |
No description provided.