diff --git a/docs/_toc.yml b/docs/_toc.yml index 1ff64785..864629bb 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -21,6 +21,7 @@ parts: - file: users_guide/running_tsmp_pdaf/input_obs - file: users_guide/running_tsmp_pdaf/input_enkfpf - file: users_guide/running_tsmp_pdaf/tsmp_pdaf_omi + - file: users_guide/running_tsmp_pdaf/lst_da - file: users_guide/running_tsmp_pdaf/grace_tws_da - file: users_guide/misc/README diff --git a/docs/users_guide/running_tsmp_pdaf/README.md b/docs/users_guide/running_tsmp_pdaf/README.md index 18dfd8f0..b44f4aca 100644 --- a/docs/users_guide/running_tsmp_pdaf/README.md +++ b/docs/users_guide/running_tsmp_pdaf/README.md @@ -15,6 +15,9 @@ COSMO](cos)). Additionally, a control file for the data assimilation Furthermore, some command line options ([Command line options](cmd)) need to be specified when TSMP-PDAF is executed. +For assimilation of land surface temperature observations into CLM, see +[Land Surface Temperature Data Assimilation](lstda). + 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 a3902a84..273f405c 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md +++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md @@ -64,6 +64,11 @@ statevec_max_layer = t_printensemble = watmin_switch = swc_mask_snow = +T_mask_snow = +T_mask_snow_depth = +increment_type = +T_mask_T = +T_max_increment = update_tws = [COSMO] @@ -488,16 +493,48 @@ CLM (standalone only). manual ) +(enkfpf:clm:update_T)= ### CLM:update_T ### -`CLM:update_T`: (integer) Flag for updating of ground and vegetation -temperature. +`CLM:update_T`: (integer) Flag for updating temperature variables in +eCLM via LST data assimilation. -Currently only CLM3.5 +State vector variables updated for each option: -- 0: No update of ground and vegetation temperature +- 0: No update of temperature variables. -- 1: Update of ground and vegetation temperature +- 1: Update of ground temperature (`t_grnd`) and vegetation + temperature (`t_veg`) directly. The simulated LST is computed from + `t_grnd` and `t_veg` using a radiometric mixing formula + ([Kustas & Anderson, 2009](https://doi.org/10.1016/j.agrformet.2009.05.016), + Eq. 7) with LAI. + +- 2: Gridcell-mean update of skin temperature (`t_skin`, not + prognostic), soil/snow temperatures (`t_soisno`, + `min(nlevgrnd, CLM:statevec_max_layer)` layers), and vegetation + temperature (`t_veg`). Each patch/column is updated based on + gridcell-mean increments and according to selected + [increment type](enkfpf:clm:increment_type). The observation + operator uses the skin temperature (TSKIN) as the simulated LST + equivalent. + +- 3: Like `2`, but additionally updates ground temperature + (`t_grnd`). State vector: `t_skin`, `t_soisno` + (`min(nlevgrnd, CLM:statevec_max_layer)` layers), `t_veg`, + `t_grnd`. + +- 4: Like `2`, but additionally updates surface water temperature + (`t_h2osfc`). State vector: `t_skin`, `t_soisno` + (`min(nlevgrnd, CLM:statevec_max_layer)` layers), `t_veg`, + `t_h2osfc`. + +- 5: Like `3`, but additionally updates surface water temperature + (`t_h2osfc`). State vector: `t_skin`, `t_soisno` + (`min(nlevgrnd, CLM:statevec_max_layer)` layers), `t_veg`, + `t_grnd`, `t_h2osfc`. + +See [Land Surface Temperature Data Assimilation](lstda) for a detailed +description. ### CLM:print_swc ### @@ -567,25 +604,36 @@ If `0` (default): Use all columns and all layers. If `1`: Use only hydrologically active columns and only layers until bedrock. +(enkfpf:clm:statevec_max_layer)= ### CLM:statevec_max_layer ### **Not yet in main branch** -`CLM:statevec_max_layer`: (integer) Number of layers to add in the -state vector. +`CLM:statevec_max_layer`: (integer) Maximum number of soil layers +included in the state vector. + +Used in two contexts: -Only used, when `CLM:statevec_allcol` and `CLM:statevec_only_active` -are switched on. +- **SWC state vector**: when `CLM:statevec_allcol` and + `CLM:statevec_only_active` are both switched on, limits the number + of soil layers included per column. +- **T state vector**: when `CLM:update_T` is `2`, `3`, `4`, or `5`, + limits the number of `t_soisno` layers included. The effective number + of layers is `min(nlevgrnd, CLM:statevec_max_layer)`. -If `25` (default): All layers are in state vector. +If `25` (default): All layers are in the state vector (CLM5 has 25 +soil layers by default, so this effectively means no restriction). -If `9`: Only the first nine layers in state vector (corresponds to 1.2 -meter). +If `9`: Only the first nine layers are included (corresponds to +approximately 1.2 m depth). For a depth profile of CLM layers, see [CLM Technical Note: 2.2.2.1 Soil Layers](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#soil-layers). +See [Land Surface Temperature Data Assimilation](lstda) for context on +the T state vector use. + ### CLM:t_printensemble ### `CLM:t_printensemble`: (integer) The timestep for the state ensemble @@ -616,12 +664,76 @@ are allowed. `CLM:swc_mask_snow`: (integer) Switch for masking columns with snow cover from SWC updates. -Snow covers larger than 1mm are switched off for the update. +Columns with snow depth ≥ 1 mm are excluded from the update. Only takes effect if `CLM:update_swc``is switched on. Default setting is `0`: No masking of columns with snow cover. +(enkfpf:clm:T_mask_snow)= +### CLM:T_mask_snow ### + +`CLM:T_mask_snow`: (integer) Switch for masking columns with snow +cover from T updates. + +When set to `1`, columns with a snow depth exceeding +`CLM:T_mask_snow_depth` are excluded from the temperature update. + +Only takes effect if `CLM:update_T` is switched on. + +Default setting is `0`: No masking of columns with snow cover. + +(enkfpf:clm:T_mask_snow_depth)= +### CLM:T_mask_snow_depth ### + +`CLM:T_mask_snow_depth`: (double) Snow depth threshold (m) used by +the snow masking condition. Columns with `snow_depth >= +CLM:T_mask_snow_depth` are excluded from the temperature update when +`CLM:T_mask_snow = 1`. + +Only takes effect if `CLM:update_T` and `CLM:T_mask_snow` are both +switched on. + +Default setting is `0.001` (1 mm). + +(enkfpf:clm:increment_type)= +### CLM:increment_type ### + +`CLM:increment_type`: (integer) Switch for changing increment type in +T-update. + +- `0`: Multiplicative increment +- `1`: Additive increment + +Only takes effect if `CLM:update_T` is switched on. + +Default setting is `0`: Multiplicative increment. + +(enkfpf:clm:T_max_increment)= +### CLM:T_max_increment ### + +`CLM:T_max_increment`: (double) Maximum allowed magnitude of the +additive temperature increment (K). + +Only takes effect if `CLM:update_T` is switched on and +`CLM:increment_type` is set to `1`. + +Default setting is `5.0`: Updates larger than 5K are not applied. + +(enkfpf:clm:T_mask_T)= +### CLM:T_mask_T ### + +`CLM:T_mask_T`: (double) Offset above freezing (K) used as the +lower-temperature masking threshold. The update is suppressed whenever + +``` +t_soisno(:,1) < 273.15 K + CLM:T_mask_T. +``` + +Only takes effect if `CLM:update_T` is switched on. + +Default setting is `0.`: Masking updates below freezing temperatures. + (enkfpf:clm:update_tws)= ### CLM:update_tws ### @@ -994,7 +1106,7 @@ Default: 0, output turned off. ## Parameter Summary ## | section | parameter | default value | - |:---------:|:-----------------------:|:-------------:| + |:----------|:------------------------|:--------------| | `[PF]` | | | | | `problemname` | \- | | | `nprocs` | 0 | @@ -1024,6 +1136,8 @@ Default: 0, output turned off. | | `problemname` | \- | | | `nprocs` | 0 | | | `update_swc` | 1 | + | | `update_texture` | 0 | + | | `update_T` | 0 | | | `print_swc` | 0 | | | `print_et` | 0 | | | `print_inc` | 0 | @@ -1033,6 +1147,11 @@ Default: 0, output turned off. | | `statevec_max_layer` | 25 | | | `t_printensemble` | -2 | | | `watmin_switch` | 0 | + | | `T_mask_snow` | 0 | + | | `T_mask_snow_depth` | 0.001 | + | | `increment_type` | 0 | + | | `T_max_increment` | 5.0 | + | | `T_mask_T` | 0.0 | | | `update_tws` | 0 | | `[COSMO]` | | | | | `nprocs` | 0 | diff --git a/docs/users_guide/running_tsmp_pdaf/input_obs.md b/docs/users_guide/running_tsmp_pdaf/input_obs.md index fbad93be..282ca622 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_obs.md +++ b/docs/users_guide/running_tsmp_pdaf/input_obs.md @@ -290,7 +290,7 @@ means the actual types should be padded with whitespaces (before or after). So, currently possible entries would be - `"GRACE "` - `"SM "` - +- `"LST "` (obs:files:clm:multi-scale-da)= ### Multi-Scale Data Assimilation observation file variables ### diff --git a/docs/users_guide/running_tsmp_pdaf/lst_da.md b/docs/users_guide/running_tsmp_pdaf/lst_da.md new file mode 100644 index 00000000..f82a81fd --- /dev/null +++ b/docs/users_guide/running_tsmp_pdaf/lst_da.md @@ -0,0 +1,196 @@ +(lstda)= +# Land Surface Temperature Data Assimilation # + +Land Surface Temperature data assimilation (LST-DA) in TSMP-PDAF enables +the assimilation of remotely-sensed land surface temperature (LST) +observations into the eCLM land-surface model. The analysis step updates +one or more temperature variables in the CLM state, and auxiliary masking +and increment-clipping parameters allow the update to be restricted to +physically meaningful regimes. + +## Configuration ## + +LST-DA is controlled by the following parameters in the [`[CLM]` section of +`enkfpf.par`](enkfpf:clm): + +- [`CLM:update_T`](enkfpf:clm:update_T): selects which CLM temperature + variables are placed in the state vector and updated after the analysis + step. +- [`CLM:statevec_max_layer`](enkfpf:clm:statevec_max_layer): limits the + number of `t_soisno` layers included in the state vector for options 2 + and 3 (default 25, i.e. all layers). +- [`CLM:T_mask_snow`](enkfpf:clm:T_mask_snow): optionally masks out + columns with snow cover from the temperature update. +- [`CLM:increment_type`](enkfpf:clm:increment_type): switches between + multiplicative (default) and additive increments. +- [`CLM:T_max_increment`](enkfpf:clm:T_max_increment): clips additive + increments to a maximum magnitude (only relevant when + `CLM:increment_type=1`). +- [`CLM:T_mask_T`](enkfpf:clm:T_mask_T): masks out columns whose + surface-layer soil temperature falls below a threshold close to + freezing. + +When LST-DA is active (`CLM:update_T != 0`), the temperature state +variables replace the soil water content (SWC) in the state vector. +Simultaneous assimilation of SWC and LST in the same PDAF update step is +not supported. + +## State Vector ## + +The state vector content depends on `CLM:update_T`. For options 2–5 +the state vector is defined at the gridcell level (one value per grid +cell); for option 1 it is defined at the patch level. + +| `update_T` | State vector variables | +|:-----------|:----------------------------------------------------------------------------------------------| +| `1` | `t_grnd` (per patch) + `t_veg` (per patch) | +| `2` | `t_skin` + `t_soisno` (`min(nlevgrnd, statevec_max_layer)` layers) + `t_veg` (per gridcell) | +| `3` | Same as `2`, plus `t_grnd` (per gridcell) | +| `4` | Same as `2`, plus `t_h2osfc` (per gridcell) | +| `5` | Same as `3`, plus `t_h2osfc` (per gridcell) | + +## Observation Operator ## + +The observation operator maps the CLM state to a simulated LST: + +- **Option 1**: The simulated LST is computed from `t_grnd` and + `t_veg` using the radiometric mixing formula of [Kustas & Anderson + (2009)](https://doi.org/10.1016/j.agrformet.2009.05.016) (Eq. 7), + weighted by leaf area index (LAI). +- **Options 2–5**: The skin temperature `t_skin` is used directly as + the simulated LST. Observations are indexed at the gridcell level + (one observation per grid cell). + +## Increment Application ## + +After the PDAF analysis step, the updated state vector values are written +back to the CLM temperature arrays. Two increment modes are available, +selected via [`CLM:increment_type`](enkfpf:clm:increment_type): + +**Multiplicative increment** (`CLM:increment_type=0`, default): +Each temperature variable is scaled by the ratio of the posterior to the +prior gridcell-mean skin temperature: + +``` +t_update = t_prior × (TSKIN_out / TSKIN_in) +``` + +**Additive increment** (`CLM:increment_type=1`): +The difference between posterior and prior gridcell-mean TSKIN is added +directly to each temperature variable: + +``` +t_update = t_prior + (TSKIN_out - TSKIN_in) +``` + +When `CLM:increment_type=1`, the increment is clipped to +`±CLM:T_max_increment` (default 5 K). Increments that exceed this +threshold are replaced by the clipped value and a warning counter is +incremented; the total number of clipped increments is printed after each +analysis step. + +## Masking ## + +Two independent masking conditions can suppress the update for individual +columns or grid cells: + +**Snow masking** (`CLM:T_mask_snow`): When set to `1`, columns with a +snow depth exceeding `CLM:T_mask_snow_depth` (default 1 mm) are excluded +from the temperature update. This avoids applying a bare-soil LST +increment to snow-covered grid cells. The threshold can be adjusted via +[`CLM:T_mask_snow_depth`](enkfpf:clm:T_mask_snow_depth). + +**Freeze masking** (`CLM:T_mask_T`): The update is suppressed whenever +the soil/snow temperature of the surface layer (`t_soisno(:,1)`) falls +below `T_freeze + CLM:T_mask_T`. With the default value of `0`, this +masks all grid cells at or below the freezing point (273.15 K). + +Both masks are evaluated independently for each patch or column. + +## Safety Checks ## + +- NaN checks are applied to all updated temperature variables; a warning + is printed to stdout if a NaN value is detected. +- Clipped additive increments are counted per analysis step and the total + is printed for `t_skin`, `t_soisno`, `t_veg`, and `t_grnd` separately. +- The update is skipped entirely for a given patch/column if the + gridcell-mean TSKIN change falls below `1e-7 K` (numerical tolerance). + +## Configuration Examples ## + +### Assimilate LST into ground and vegetation temperature (option 1) ### + +```text +[CLM] +update_T = 1 +increment_type = 0 +``` + +The simulated LST is computed from `t_grnd` and `t_veg` via the +[Kustas & Anderson (2009)](https://doi.org/10.1016/j.agrformet.2009.05.016) +radiometric mixing formula. Both variables are updated with a +multiplicative increment. + +### Assimilate LST into skin and soil temperature (option 2) ### + +```text +[CLM] +update_T = 2 +increment_type = 0 +T_mask_snow = 1 +T_mask_T = 0.0 +``` + +TSKIN is used as the simulated LST. After analysis, `t_skin`, all +`t_soisno` layers, and `t_veg` are scaled by the TSKIN increment factor. +Snow-covered columns are excluded from the update. + +### Assimilate LST with additive increment and clipping (option 2) ### + +```text +[CLM] +update_T = 2 +increment_type = 1 +T_max_increment = 3.0 +T_mask_T = 2.0 +``` + +The TSKIN increment is applied additively and clipped to ±3 K. Grid cells +whose surface-layer temperature falls below 275.15 K (freezing + 2 K) are +excluded. + +### Assimilate LST including ground temperature (option 3) ### + +```text +[CLM] +update_T = 3 +increment_type = 0 +T_mask_snow = 1 +``` + +Like option 2, but `t_grnd` is additionally included in the state vector +and updated. Snow-covered columns are masked out. + +### Assimilate LST including surface water temperature (option 4) ### + +```text +[CLM] +update_T = 4 +increment_type = 0 +T_mask_snow = 1 +``` + +Like option 2, but `t_h2osfc` is additionally included in the state vector +and updated. Snow-covered columns are masked out. + +### Assimilate LST including ground and surface water temperature (option 5) ### + +```text +[CLM] +update_T = 5 +increment_type = 0 +T_mask_snow = 1 +``` + +Like option 3, but `t_h2osfc` is additionally included in the state vector +and updated. Snow-covered columns are masked out. diff --git a/interface/framework/Makefile b/interface/framework/Makefile index 35719909..17d74b96 100644 --- a/interface/framework/Makefile +++ b/interface/framework/Makefile @@ -77,7 +77,8 @@ MOD_ASSIM = mod_parallel_pdaf.o \ # Routines of observation handling (PDAF-OMI) MOD_USER_PDAFOMI = obs_GRACE_pdafomi.o \ - obs_SM_pdafomi.o + obs_SM_pdafomi.o \ + obs_LST_pdafomi.o # Model routines used with PDAF OBJ_MODEL_PDAF =pdaf_terrsysmp.o\ diff --git a/interface/framework/callback_obs_pdafomi.F90 b/interface/framework/callback_obs_pdafomi.F90 index 35400936..65e68730 100644 --- a/interface/framework/callback_obs_pdafomi.F90 +++ b/interface/framework/callback_obs_pdafomi.F90 @@ -36,12 +36,15 @@ SUBROUTINE init_dim_obs_pdafomi(step, dim_obs) use enkf_clm_mod, only: clmupdate_swc, clmupdate_tws + use enkf_clm_mod, only: clmupdate_T ! Include functions for different observations USE obs_GRACE_pdafomi, ONLY: assim_GRACE USE obs_GRACE_pdafomi, ONLY: init_dim_obs_GRACE USE obs_SM_pdafomi, ONLY: assim_SM USE obs_SM_pdafomi, ONLY: init_dim_obs_SM + USE obs_LST_pdafomi, ONLY: assim_LST + USE obs_LST_pdafomi, ONLY: init_dim_obs_LST !USE obs_ST_pdafomi, ONLY: assim_C !USE obs_ST_pdafomi, ONLY: init_dim_obs_C @@ -57,6 +60,7 @@ SUBROUTINE init_dim_obs_pdafomi(step, dim_obs) ! *** Local variables *** INTEGER :: dim_obs_GRACE ! Observation dimensions INTEGER :: dim_obs_SM ! Observation dimensions + INTEGER :: dim_obs_LST ! Observation dimensions !INTEGER :: dim_obs_C ! Observation dimensions @@ -67,6 +71,7 @@ SUBROUTINE init_dim_obs_pdafomi(step, dim_obs) ! Initialize number of observations dim_obs_GRACE = 0 dim_obs_SM = 0 + dim_obs_LST = 0 !dim_obs_C = 0 @@ -82,15 +87,17 @@ SUBROUTINE init_dim_obs_pdafomi(step, dim_obs) if (mype_world==0) then write(*,*)'PDAF-OMI-DEBUG: assim_GRACE=', assim_GRACE write(*,*)'PDAF-OMI-DEBUG: assim_SM=', assim_SM + write(*,*)'PDAF-OMI-DEBUG: assim_LST=', assim_LST ! write(*,*)'PDAF-OMI assim_C=', assim_C end if #endif IF (assim_GRACE) CALL init_dim_obs_GRACE(step, dim_obs_GRACE) IF (assim_SM) CALL init_dim_obs_SM(step, dim_obs_SM) + IF (assim_LST) CALL init_dim_obs_LST(step, dim_obs_LST) !IF (assim_C) CALL init_dim_obs_C(step, dim_obs_C) - dim_obs = dim_obs_GRACE + dim_obs_SM! + dim_obs_C + dim_obs = dim_obs_GRACE + dim_obs_SM + dim_obs_LST! + dim_obs_C END SUBROUTINE init_dim_obs_pdafomi @@ -109,6 +116,8 @@ SUBROUTINE obs_op_pdafomi(step, dim_p, dim_obs, state_p, ostate) USE obs_GRACE_pdafomi, ONLY: obs_op_GRACE USE obs_SM_pdafomi, ONLY: assim_SM USE obs_SM_pdafomi, ONLY: obs_op_SM + USE obs_LST_pdafomi, ONLY: assim_LST + USE obs_LST_pdafomi, ONLY: obs_op_LST !USE obs_C_pdafomi, ONLY: assim_C !USE obs_C_pdafomi, ONLY: obs_op_C @@ -139,6 +148,7 @@ SUBROUTINE obs_op_pdafomi(step, dim_p, dim_obs, state_p, ostate) IF (assim_GRACE) CALL obs_op_GRACE(dim_p, dim_obs, state_p, ostate) IF (assim_SM) CALL obs_op_SM(dim_p, dim_obs, state_p, ostate) + IF (assim_LST) CALL obs_op_LST(dim_p, dim_obs, state_p, ostate) !IF (assim_C) CALL obs_op_C(dim_p, dim_obs, state_p, ostate) END SUBROUTINE obs_op_pdafomi @@ -158,6 +168,8 @@ SUBROUTINE init_dim_obs_l_pdafomi(domain_p, step, dim_obs, dim_obs_l) USE obs_GRACE_pdafomi, ONLY: init_dim_obs_l_GRACE USE obs_SM_pdafomi, ONLY: assim_SM USE obs_SM_pdafomi, ONLY: init_dim_obs_l_SM + USE obs_LST_pdafomi, ONLY: assim_LST + USE obs_LST_pdafomi, ONLY: init_dim_obs_l_LST !USE obs_C_pdafomi, ONLY: assim_C !USE obs_C_pdafomi, ONLY: init_dim_obs_l_C @@ -178,6 +190,7 @@ SUBROUTINE init_dim_obs_l_pdafomi(domain_p, step, dim_obs, dim_obs_l) IF (assim_GRACE) CALL init_dim_obs_l_GRACE(domain_p, step, dim_obs, dim_obs_l) IF (assim_SM) CALL init_dim_obs_l_SM(domain_p, step, dim_obs, dim_obs_l) + IF (assim_LST) CALL init_dim_obs_l_LST(domain_p, step, dim_obs, dim_obs_l) !IF (assim_C) CALL init_dim_obs_l_C(domain_p, step, dim_obs, dim_obs_l) END SUBROUTINE init_dim_obs_l_pdafomi @@ -198,6 +211,8 @@ SUBROUTINE localize_covar_pdafomi(dim_p, dim_obs, HP_p, HPH) USE obs_GRACE_pdafomi, ONLY: localize_covar_GRACE USE obs_SM_pdafomi, ONLY: assim_SM USE obs_SM_pdafomi, ONLY: localize_covar_SM + USE obs_LST_pdafomi, ONLY: assim_LST + USE obs_LST_pdafomi, ONLY: localize_covar_LST !USE obs_C_pdafomi, ONLY: assim_C !USE obs_C_pdafomi, ONLY: localize_covar_C @@ -232,6 +247,7 @@ SUBROUTINE localize_covar_pdafomi(dim_p, dim_obs, HP_p, HPH) IF (assim_GRACE) CALL localize_covar_GRACE(dim_p, dim_obs, HP_p, HPH, coords_p) IF (assim_SM) CALL localize_covar_SM(dim_p, dim_obs, HP_p, HPH, coords_p) + IF (assim_LST) CALL localize_covar_LST(dim_p, dim_obs, HP_p, HPH, coords_p) !IF (assim_C) CALL localize_covar_C(dim_p, dim_obs, HP_p, HPH, coords_p) @@ -250,6 +266,8 @@ SUBROUTINE add_obs_err_pdafomi(step, dim_obs, C) USE obs_GRACE_pdafomi, ONLY: add_obs_err_GRACE USE obs_SM_pdafomi, ONLY: assim_SM USE obs_SM_pdafomi, ONLY: add_obs_err_SM + USE obs_LST_pdafomi, ONLY: assim_LST + USE obs_LST_pdafomi, ONLY: add_obs_err_LST !USE obs_C_pdafomi, ONLY: assim_C !USE obs_C_pdafomi, ONLY: add_obs_err_C @@ -265,6 +283,7 @@ SUBROUTINE add_obs_err_pdafomi(step, dim_obs, C) REAL :: variance_obs ! variance of observations IF (assim_GRACE) CALL add_obs_err_GRACE(step, dim_obs, C) IF (assim_SM) CALL add_obs_err_SM(step, dim_obs, C) + IF (assim_LST) CALL add_obs_err_LST(step, dim_obs, C) !IF (assim_C) CALL add_obs_err_C(step, dim_obs, C) END SUBROUTINE add_obs_err_pdafomi @@ -276,6 +295,8 @@ SUBROUTINE init_obscovar_pdafomi(step, dim_obs, dim_obs_p, covar, m_state_p, isd USE obs_GRACE_pdafomi, ONLY: init_obscovar_GRACE USE obs_SM_pdafomi, ONLY: assim_SM USE obs_SM_pdafomi, ONLY: init_obscovar_SM + USE obs_LST_pdafomi, ONLY: assim_LST + USE obs_LST_pdafomi, ONLY: init_obscovar_LST !USE obs_C_pdafomi, ONLY: assim_C !USE obs_C_pdafomi, ONLY: init_obscovar_C IMPLICIT NONE @@ -293,6 +314,7 @@ SUBROUTINE init_obscovar_pdafomi(step, dim_obs, dim_obs_p, covar, m_state_p, isd IF (assim_GRACE) CALL init_obscovar_GRACE(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) IF (assim_SM) CALL init_obscovar_SM(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) + IF (assim_LST) CALL init_obscovar_LST(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) !IF (assim_C) CALL init_obscovar_C(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) END SUBROUTINE init_obscovar_pdafomi @@ -304,6 +326,8 @@ SUBROUTINE prodRinvA_pdafomi(step, dim_obs_p, rank, obs_p, A_p, C_p) use obs_GRACE_pdafomi, ONLY: prodRinvA_GRACE use obs_SM_pdafomi, ONLY: assim_SM use obs_SM_pdafomi, ONLY: prodRinvA_SM + use obs_LST_pdafomi, ONLY: assim_LST + use obs_LST_pdafomi, ONLY: prodRinvA_LST !use obs_C_pdafomi, ONLY: assim_C !use obs_C_pdafomi, ONLY: prodRinvA_C @@ -317,6 +341,7 @@ SUBROUTINE prodRinvA_pdafomi(step, dim_obs_p, rank, obs_p, A_p, C_p) IF (assim_GRACE) CALL prodRinvA_GRACE(step, dim_obs_p, rank, obs_p, A_p, C_p) IF (assim_SM) CALL prodRinvA_SM(step, dim_obs_p, rank, obs_p, A_p, C_p) + IF (assim_LST) CALL prodRinvA_LST(step, dim_obs_p, rank, obs_p, A_p, C_p) !IF (assim_C) CALL prodRinvA_C(step, dim_obs_p, rank, obs_p, A_p, C_p) END SUBROUTINE prodRinvA_pdafomi @@ -328,6 +353,8 @@ SUBROUTINE prodRinvA_l_pdafomi(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) use obs_GRACE_pdafomi, ONLY: prodRinvA_l_GRACE use obs_SM_pdafomi, ONLY: assim_SM use obs_SM_pdafomi, ONLY: prodRinvA_l_SM + use obs_LST_pdafomi, ONLY: assim_LST + use obs_LST_pdafomi, ONLY: prodRinvA_l_LST !use obs_C_pdafomi, ONLY: assim_C !use obs_C_pdafomi, ONLY: prodRinvA_l_C @@ -343,6 +370,7 @@ SUBROUTINE prodRinvA_l_pdafomi(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) IF (assim_GRACE) CALL prodRinvA_l_GRACE(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) IF (assim_SM) CALL prodRinvA_l_SM(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) + IF (assim_LST) CALL prodRinvA_l_LST(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) !IF (assim_C) CALL prodRinvA_l_C(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) END SUBROUTINE prodRinvA_l_pdafomi @@ -353,6 +381,8 @@ SUBROUTINE deallocate_obs_pdafomi() use obs_GRACE_pdafomi, ONLY: deallocate_obs_GRACE use obs_SM_pdafomi, ONLY: assim_SM use obs_SM_pdafomi, ONLY: deallocate_obs_SM + use obs_LST_pdafomi, ONLY: assim_LST + use obs_LST_pdafomi, ONLY: deallocate_obs_LST !use obs_C_pdafomi, ONLY: assim_C !use obs_C_pdafomi, ONLY: deallocate_obs_C @@ -360,6 +390,7 @@ SUBROUTINE deallocate_obs_pdafomi() IF (assim_GRACE) CALL deallocate_obs_GRACE() IF (assim_SM) CALL deallocate_obs_SM() + IF (assim_LST) CALL deallocate_obs_LST() !IF (assim_C) CALL deallocate_obs_C() END SUBROUTINE deallocate_obs_pdafomi diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 11547a4e..6cf4b5e8 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -129,6 +129,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) #ifdef CLMFIVE use GridcellType, only: grc use ColumnType, only : col + use PatchType, only : patch ! use GetGlobalValuesMod, only: GetGlobalWrite ! use clm_varcon, only: nameg use enkf_clm_mod, only: state_clm2pdaf_p @@ -145,6 +146,8 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) USE enkf_clm_mod, only: domain_def_clm USE enkf_clm_mod, only: get_interp_idx use enkf_clm_mod, only: clmstatevec_allcol + use enkf_clm_mod, only: clmupdate_swc + use enkf_clm_mod, only: clmupdate_T !hcp end #endif #endif @@ -164,9 +167,11 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) integer :: ierror INTEGER :: max_var_id ! Multi-scale DA INTEGER :: sum_dim_obs_p + INTEGER :: p ! CLM Patch index INTEGER :: c ! CLM Column index INTEGER :: g ! CLM Gridcell index INTEGER :: cg + INTEGER :: pg INTEGER :: i,j,k ! Counters INTEGER :: cnt ! Counters INTEGER :: cnt_interp ! Counter for interpolation grid cells @@ -935,6 +940,8 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) do i = 1, dim_obs obs(i) = clm_obs(i) + if(clmupdate_swc==1) then + do g = begg,endg newgridcell = .true. @@ -1005,9 +1012,100 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end if end if - end do end do + + else if(clmupdate_T/=0) then +#ifdef CLMFIVE + ! patch loop + do g = begg,endg + newgridcell = .true. + + do p = begp,endp + + pg = patch%gridcell(p) + + if(pg == g) then + if(newgridcell) then + ! Sets first patch/column in a gridcell. TODO: Make + ! patch / column information part of the observation + ! file + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + + ! 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 + + end if + + newgridcell = .false. + + end if + end if + + end do + end do +#else + ! gridcell loop + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + end if + + !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 + + end do +#endif + else + + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING unsupported update in setting obs_index_p." + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING using default gridcell loop." + ! call abort_parallel() + + ! gridcell loop with layer information as default! + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if + + !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 + + end do + + end if + end do if(obs_interp_switch==1) then diff --git a/interface/framework/init_dim_obs_l_pdaf.F90 b/interface/framework/init_dim_obs_l_pdaf.F90 index 04c25803..969b1d11 100755 --- a/interface/framework/init_dim_obs_l_pdaf.F90 +++ b/interface/framework/init_dim_obs_l_pdaf.F90 @@ -319,6 +319,11 @@ SUBROUTINE init_dim_obs_l_pdaf(domain_p, step, dim_obs_f, dim_obs_l) #endif !------------------------------------------------------------------------ +#ifdef PDAF_DEBUG + WRITE(*, '(a,x,a,i5,x,a,i10,x,a,i10)') "TSMP-PDAF-debug", "mype(f)=", mype_filter, & + "init_dim_obs_l_pdaf: domain_p=", domain_p, "dim_obs_l=", dim_obs_l +#endif + ! kuw: allocate and determine local observation index and distance !#ifndef CLMSA ! Initialize index array for local observations in full observed vector diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 48f29c08..541be7ff 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -125,6 +125,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) #ifdef CLMFIVE use GridcellType, only: grc use ColumnType, only : col + use PatchType, only : patch ! use GetGlobalValuesMod, only: GetGlobalWrite ! use clm_varcon, only: nameg use enkf_clm_mod, only: state_clm2pdaf_p @@ -141,6 +142,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) USE enkf_clm_mod, only: domain_def_clm USE enkf_clm_mod, only: get_interp_idx use enkf_clm_mod, only: clmstatevec_allcol + use enkf_clm_mod, only: clmupdate_swc + use enkf_clm_mod, only: clmupdate_T !hcp end #endif #endif @@ -160,9 +163,11 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) integer :: ierror INTEGER :: max_var_id ! Multi-scale DA INTEGER :: sum_dim_obs_p + INTEGER :: p ! CLM Patch index INTEGER :: c ! CLM Column index INTEGER :: g ! CLM Gridcell index INTEGER :: cg + INTEGER :: pg INTEGER :: i,j,k ! Counters INTEGER :: cnt ! Counters INTEGER :: cnt_interp ! Counter for interpolation grid cells @@ -928,6 +933,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) do i = 1, dim_obs obs(i) = clm_obs(i) + if(clmupdate_swc==1) then + do g = begg,endg newgridcell = .true. @@ -998,9 +1005,100 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end if end if - end do end do + + else if(clmupdate_T/=0) then +#ifdef CLMFIVE + ! patch loop + do g = begg,endg + newgridcell = .true. + + do p = begp,endp + + pg = patch%gridcell(p) + + if(pg == g) then + if(newgridcell) then + ! Sets first patch/column in a gridcell. TODO: Make + ! patch / column information part of the observation + ! file + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + + ! 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 + + end if + + newgridcell = .false. + + end if + end if + + end do + end do +#else + ! gridcell loop + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + end if + + !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 + + end do +#endif + else + + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING unsupported update in setting obs_index_p." + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING using default gridcell loop." + ! call abort_parallel() + + ! gridcell loop with layer information as default! + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if + + !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 + + end do + + end if + end do if(obs_interp_switch==1) then diff --git a/interface/framework/init_pdaf.F90 b/interface/framework/init_pdaf.F90 index eb68b077..6d689b16 100644 --- a/interface/framework/init_pdaf.F90 +++ b/interface/framework/init_pdaf.F90 @@ -102,10 +102,12 @@ SUBROUTINE init_pdaf() #ifdef CLMFIVE USE obs_GRACE_pdafomi, ONLY: assim_GRACE USE obs_SM_pdafomi, ONLY: assim_SM + USE obs_LST_pdafomi, ONLY: assim_LST !USE obs_ST_pdafomi, ONLY: assim_C USE enkf_clm_mod, ONLY: clmupdate_tws USE enkf_clm_mod, ONLY: clmupdate_swc + USE enkf_clm_mod, ONLY: clmupdate_T ! use enkf_clm_mod, only: clmupdate_C #endif #endif @@ -247,6 +249,7 @@ SUBROUTINE init_pdaf() #ifdef CLMFIVE assim_GRACE = (clmupdate_tws /= 0) assim_SM = (clmupdate_swc /= 0) + assim_LST = (clmupdate_T /= 0) ! assim_C = (clmupdate_C /= 0) #endif #endif diff --git a/interface/framework/init_pdaf_parse.F90 b/interface/framework/init_pdaf_parse.F90 index 1a15295a..c9715825 100644 --- a/interface/framework/init_pdaf_parse.F90 +++ b/interface/framework/init_pdaf_parse.F90 @@ -59,9 +59,12 @@ SUBROUTINE init_pdaf_parse() use mod_assimilation,& only: cradius_GRACE, sradius_GRACE, & cradius_SM, sradius_SM + use mod_assimilation, only: cradius_LST + use mod_assimilation, only: sradius_LST #ifdef CLMFIVE use obs_GRACE_pdafomi, only: rms_obs_GRACE use obs_SM_pdafomi, only: rms_obs_SM + use obs_LST_pdafomi, only: rms_obs_LST #endif #if defined CLMSA #ifdef CLMFIVE @@ -106,6 +109,9 @@ SUBROUTINE init_pdaf_parse() rms_obs_SM = rms_obs ! backward compatibility handle = 'rms_obs_SM' ! RMS error for SM observations CALL parse(handle, rms_obs_SM) + rms_obs_LST = rms_obs ! backward compatibility + handle = 'rms_obs_LST' ! RMS error for LST observations + CALL parse(handle, rms_obs_LST) ! rms_obs_C = rms_obs ! backward compatibility ! handle = 'rms_obs_C' ! RMS error for C observations ! CALL parse(handle, rms_obs_C) @@ -172,6 +178,12 @@ SUBROUTINE init_pdaf_parse() sradius_SM = sradius ! For backward compatibility handle = 'sradius_SM' ! Set support radius for SM observations call parse(handle, sradius_SM) + cradius_LST = cradius ! For backward compatibility + handle = 'cradius_LST' ! Set cut-off radius for LST observations + call parse(handle, cradius_LST) + sradius_LST = sradius ! For backward compatibility + handle = 'sradius_LST' ! Set support radius for LST observations + call parse(handle, sradius_LST) ! Setting for file output handle = 'filename' ! Set name of output file diff --git a/interface/framework/localize_covar_pdaf.F90 b/interface/framework/localize_covar_pdaf.F90 index b0865fb8..eb93a453 100644 --- a/interface/framework/localize_covar_pdaf.F90 +++ b/interface/framework/localize_covar_pdaf.F90 @@ -40,7 +40,6 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH) USE shr_kind_mod , only : r8 => shr_kind_r8 USE mod_read_obs, ONLY: clmobs_lon USE mod_read_obs, ONLY: clmobs_lat - USE enkf_clm_mod, ONLY: clmupdate_T USE enkf_clm_mod, ONLY: clm_begc USE enkf_clm_mod, ONLY: clm_endc USE enkf_clm_mod, ONLY: state_pdaf2clm_c_p diff --git a/interface/framework/mod_assimilation.F90 b/interface/framework/mod_assimilation.F90 index 7e0f92ec..3abf5538 100755 --- a/interface/framework/mod_assimilation.F90 +++ b/interface/framework/mod_assimilation.F90 @@ -274,6 +274,10 @@ MODULE mod_assimilation REAL :: sradius_GRACE REAL :: cradius_SM REAL :: sradius_SM + REAL :: cradius_LST + REAL :: sradius_LST + ! REAL :: cradius_C + ! REAL :: sradius_C ! ! SEIK-subtype4/LSEIK-subtype4/ESTKF/LESTKF INTEGER :: type_sqrt ! Type of the transform matrix square-root ! (0) symmetric square root diff --git a/interface/framework/mod_read_obs.F90 b/interface/framework/mod_read_obs.F90 index 46c0d0f9..ae23a408 100755 --- a/interface/framework/mod_read_obs.F90 +++ b/interface/framework/mod_read_obs.F90 @@ -974,6 +974,13 @@ subroutine update_obs_type(obs_type_str) clmupdate_T = 0 clmupdate_texture = 0 + case ('LST') + ! TODO: General values for clmupdate_T (saving previous values) + clmupdate_tws = 0 + clmupdate_swc = 0 + clmupdate_T = 4 + clmupdate_texture = 0 + ! case ('C') ! clmupdate_tws = 0 ! clmupdate_swc = 0 diff --git a/interface/framework/obs_LST_pdafomi.F90 b/interface/framework/obs_LST_pdafomi.F90 new file mode 100644 index 00000000..a5a05fbb --- /dev/null +++ b/interface/framework/obs_LST_pdafomi.F90 @@ -0,0 +1,995 @@ +!> PDAF-OMI observation module for type LST observations +!! +!! This module handles operations for one data type (called 'module-type' below): +!! OBSTYPE = LST +!! +!! __Observation type LST:__ +!! The observation type LST are Land Surface Temperature observations. +!! +!! The subroutines in this module are for the particular handling of +!! a single observation type. +!! The routines are called by the different call-back routines of PDAF +!! usually by callback_obs_pdafomi.F90 +!! Most of the routines are generic so that in practice only 2 routines +!! need to be adapted for a particular data type. These are the routines +!! for the initialization of the observation information (init_dim_obs) +!! and for the observation operator (obs_op). +!! +!! The module and the routines are named according to the observation type. +!! This allows to distinguish the observation type and the routines in this +!! module from other observation types. +!! +!! The module uses two derived data types (obs_f and obs_l), which contain +!! all information about the full and local observations. Only variables +!! of the type obs_f need to be initialized in this module. The variables +!! in the type obs_l are initilized by the generic routines from PDAFomi. +!! +!! Three observation operator modes are supported via clmupdate_T: +!! * clmupdate_T==1: Kustas (2009) formula combining ground (TG) and +!! vegetation (TV) temperature: +!! LST = (exp(-tau/2)*TG^4 + (1-exp(-tau/2))*TV^4)^0.25 +!! where tau = clm_paramarr(obs_index) (vegetation optical depth) +!! * clmupdate_T==2,3: Direct TSKIN mapping (simulated LST = TSKIN) +!! +!! The state vector for LST is patch-based (not column-based as for SM). +!! For clmupdate_T==1 the state vector contains TG (first) and TV (second +!! half, offset by clm_varsize). +!! +!! These 2 routines need to be adapted for the particular observation type: +!! * init_dim_obs_LST \n +!! Count number of process-local and full observations; +!! initialize vector of observations and their inverse variances; +!! initialize coordinate array and index array for indices of +!! observed elements of the state vector. +!! * obs_op_LST \n +!! observation operator to get full observation vector of this type. +!! +!! __Revision history:__ +!! * 2024 - Initial code based on obs_SM_pdafomi.F90 (Yorck Ewerdwalbesloh) +!! adapted for LST-DA (clmupdate_T) +!! + +! Author: Based on obs_SM_pdafomi.F90 by Yorck Ewerdwalbesloh; +! adapted for LST-DA (clmupdate_T) + +#ifdef CLMFIVE +MODULE obs_LST_pdafomi + + USE mod_parallel_pdaf, & + ONLY: mype_filter ! Rank of filter process + USE PDAFomi, & + ONLY: obs_f, obs_l ! Declaration of observation data types + + IMPLICIT NONE + SAVE + + PUBLIC + + ! Variables which are inputs to the module (usually set in init_pdaf) + LOGICAL :: assim_LST !< Whether to assimilate this data type + REAL :: rms_obs_LST !< Observation error standard deviation (for constant errors) + + ! longitude and latitude of grid cells and observation cells + INTEGER, ALLOCATABLE :: longxy(:), latixy(:), longxy_obs(:), latixy_obs(:) + + ! Module-local obs-to-state index array (avoids conflict with global obs_index_p + ! from mod_assimilation when SM and LST are both active simultaneously) + INTEGER, ALLOCATABLE :: obs_index_p_LST(:) + + ! ********************************************************* + ! *** Data type obs_f defines the full observations by *** + ! *** internally shared variables of the module *** + ! ********************************************************* + + ! Declare instances of observation data types used here + TYPE(obs_f), TARGET, PUBLIC :: thisobs ! full observation + TYPE(obs_l), TARGET, PUBLIC :: thisobs_l ! local observation + + !$OMP THREADPRIVATE(thisobs_l) + + + !------------------------------------------------------------------------------- + + CONTAINS + + !> Initialize information on the module-type observation + !! + !! The routine is called by each filter process. + !! at the beginning of the analysis step before + !! the loop through all local analysis domains. + !! + !! It has to count the number of observations of the + !! observation type handled in this module according + !! to the current time step for all observations + !! required for the analyses in the loop over all local + !! analysis domains on the PE-local state domain. + !! +SUBROUTINE init_dim_obs_LST(step, dim_obs) + + USE mpi, ONLY: MPI_INTEGER + USE mpi, ONLY: MPI_DOUBLE_PRECISION + USE mpi, ONLY: MPI_SUM + USE mpi, ONLY: MPI_IN_PLACE + USE mod_parallel_pdaf, & + ONLY: mype_filter, comm_filter, npes_filter, abort_parallel, & + mype_world + USE mod_assimilation, & + ONLY: obs_filename, & + obs_pdaf2nc, obs_nc2pdaf, & + local_dims_obs, & + local_disp_obs, & + longxy_obs_floor, latixy_obs_floor, & + screen, cradius_LST + + USE PDAFomi, & + ONLY: PDAFomi_gather_obs, pi + + Use mod_read_obs, & + only: multierr, read_obs_nc_type + + use enkf_clm_mod, only: state_clm2pdaf_p + + use enkf_clm_mod, only: domain_def_clm + + use mod_parallel_pdaf, & + only: abort_parallel + + use shr_kind_mod, only: r8 => shr_kind_r8 + + use GridcellType, only: grc + + use clm_varcon, only: ispval + + use decompMod , only : get_proc_bounds, get_proc_global + + use PatchType , only : patch + + use enkf_clm_mod, only: get_interp_idx + + use mod_tsmp, only: da_print_obs_index + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: step !< Current time step + INTEGER, INTENT(inout) :: dim_obs !< Dimension of full observation vector + + ! *** Local variables *** + INTEGER :: i, p, g, pg ! Counters + INTEGER :: cnt ! Counter + INTEGER :: dim_obs_p ! Number of process-local observations + REAL, ALLOCATABLE :: obs_p(:) ! PE-local observation vector + REAL, ALLOCATABLE :: obs_g(:) ! Global observation vector + REAL, ALLOCATABLE :: ivar_obs_p(:) ! PE-local inverse observation error variance + REAL, ALLOCATABLE :: ocoord_p(:,:) ! PE-local observation coordinates + character (len = 110) :: current_observation_filename + + + character(len=20) :: obs_type_name ! name of observation type (e.g. GRACE, SM, ST, ...) + + ! clm-observation arrays are local (as opposed to non-OMI, where + ! they are found in module mod_read_obs) + REAL, ALLOCATABLE :: lon_obs(:) + REAL, ALLOCATABLE :: lat_obs(:) + INTEGER, ALLOCATABLE :: layer_obs(:) + REAL, ALLOCATABLE :: dr_obs(:) + REAL, ALLOCATABLE :: obserr(:) + REAL, ALLOCATABLE :: obscov(:,:) + + integer :: begp, endp ! per-proc beginning and ending pft indices + 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 :: numg ! total number of gridcells across all processors + integer :: numl ! total number of landunits across all processors + integer :: numc ! total number of columns across all processors + integer :: nump ! total number of pfts across all processors + + real(r8), pointer :: lon(:) + real(r8), pointer :: lat(:) + + integer :: ierror + + real :: deltax, deltay + + logical :: obs_snapped !Switch for checking multiple observation counts + logical :: newgridcell + + INTEGER :: sum_dim_obs_p + + character (len = 27) :: fn !TSMP-PDAF: function name for obs_index_p output + + + + ! ********************************************* + ! *** Initialize full observation dimension *** + ! ********************************************* + + IF (mype_filter==0) & + WRITE (*,*) 'PDAF-OMI: Assimilate observations - obs type LST' + + IF (assim_LST) thisobs%doassim = 1 + + ! Geographic distance with haversine formula + thisobs%disttype = 3 + thisobs%ncoord = 2 + + obs_type_name = 'LST' + + ! ********************************** + ! *** Read PE-local observations *** + ! ********************************** + + if(mype_filter==0 .and. screen > 2) then + write(*,*)'PDAF-OMI: load observations from type LST' + end if + + ! Set name of current NetCDF observation file + write(current_observation_filename, '(a, i5.5)') trim(obs_filename)//'.', step + + ! if I'm root in filter, read the nc file + if (mype_filter == 0) then + ! Read current NetCDF observation file + call read_obs_nc_type(current_observation_filename, obs_type_name, & + dim_obs, obs_g, lon_obs, lat_obs, layer_obs, & + dr_obs, obserr, obscov) + end if + + ! Broadcast first variables + ! ------------------------- + ! Dimension of observation vector + call mpi_bcast(dim_obs, 1, MPI_INTEGER, 0, comm_filter, ierror) + + ! Handle case if observation file contains no observations of type + ! LST + + ! Used in joint DA for handling the zero-size observation vector, + ! when this observation is not present. + if (dim_obs == 0) then + if (mype_filter==0 .and. screen > 2) then + write(*,*)'TSMP-PDAF mype(w) =', mype_world, & + ': No observations of type LST found in file ', & + trim(current_observation_filename) + end if + dim_obs_p = 0 + if (allocated(obs_p)) deallocate(obs_p) + if (allocated(ivar_obs_p)) deallocate(ivar_obs_p) + if (allocated(ocoord_p)) deallocate(ocoord_p) + if (allocated(thisobs%id_obs_p)) deallocate(thisobs%id_obs_p) + ALLOCATE(obs_p(1)) + ALLOCATE(ivar_obs_p(1)) + ALLOCATE(ocoord_p(2, 1)) + ALLOCATE(thisobs%id_obs_p(1, 1)) + thisobs%id_obs_p(1, 1) = 0 + thisobs%infile = 0 + CALL PDAFomi_gather_obs(thisobs, dim_obs_p, obs_p, ivar_obs_p, ocoord_p, & + thisobs%ncoord, cradius_LST, dim_obs) + if(mype_filter==0) DEALLOCATE(obs_g) + DEALLOCATE(obs_p, ocoord_p, ivar_obs_p) + return + end if + + ! Switch for vector of observation errors or observation covariances + call mpi_bcast(multierr, 1, MPI_INTEGER, 0, comm_filter, ierror) + + ! Allocate observation arrays for non-root procs + ! ---------------------------------------------- + if (mype_filter /= 0) then ! for all non-master proc + if(allocated(obs_g)) deallocate(obs_g) + allocate(obs_g(dim_obs)) + if(allocated(lon_obs)) deallocate(lon_obs) + allocate(lon_obs(dim_obs)) + if(allocated(lat_obs)) deallocate(lat_obs) + allocate(lat_obs(dim_obs)) + if(allocated(dr_obs)) deallocate(dr_obs) + allocate(dr_obs(2)) + if(allocated(layer_obs)) deallocate(layer_obs) + allocate(layer_obs(dim_obs)) + if(multierr==1) then + if(allocated(obserr)) deallocate(obserr) + allocate(obserr(dim_obs)) + end if + end if + + ! Broadcast the global observation arrays + ! --------------------------------------- + call mpi_bcast(obs_g, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + if(multierr==1) call mpi_bcast(obserr, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + call mpi_bcast(lon_obs, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + call mpi_bcast(lat_obs, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + call mpi_bcast(dr_obs, 2, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + call mpi_bcast(layer_obs, dim_obs, MPI_INTEGER, 0, comm_filter, ierror) + + thisobs%infile = 1 + + if (mype_filter==0 .and. screen > 2) then + write(*,*)'PDAF-OMI: Done: load observations from type LST' + end if + + ! CLM grid information + ! -------------------- + ! longxy/latixy index arrays are used for grid cell matching + ! when is_use_dr is .false. (index-based instead of distance-based snapping). + ! This is currently turned off for LSTDA. + + ! ! Generate CLM index arrays from lon/lat values + ! call domain_def_clm(lon_obs, lat_obs, dim_obs, longxy, latixy, longxy_obs, latixy_obs) + + ! Obtain CLM lon/lat information + lon => grc%londeg + lat => grc%latdeg + + ! Obtain CLM index information + call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp) + call get_proc_global(numg, numl, numc, nump) + + ! Number of observations in process-local domain + ! ---------------------------------------------- + dim_obs_p = 0 + + ! id_obs_p: placeholder to satisfy PDAFomi internal check; + ! actual obs-to-state mapping is in obs_index_p_LST + if(allocated(thisobs%id_obs_p)) deallocate(thisobs%id_obs_p) + allocate(thisobs%id_obs_p(1, 1)) + thisobs%id_obs_p(1, 1) = 1 + + ! *** Count PE-local observations *** + ! LST uses patch loop: one observation per gridcell, assigned to + ! the first patch found in that gridcell. + do i = 1, dim_obs + obs_snapped = .false. + do g = begg,endg + newgridcell = .true. + do p = begp,endp + pg = patch%gridcell(p) + if(pg == g) then + if(newgridcell) then + + deltax = abs(lon(g)-lon_obs(i)) + if (deltax > 180.0) then + deltax = 360.0 - deltax + end if + deltay = abs(lat(g)-lat_obs(i)) + + if((deltax<=dr_obs(1)).and.(deltay<=dr_obs(2))) then + dim_obs_p = dim_obs_p + 1 + + if(obs_snapped) then + print *, "TSMP-PDAF mype(w)=", mype_world, & + ": ERROR Observation snapped at multiple grid cells." + print *, "i=", i + call abort_parallel() + end if + obs_snapped = .true. + end if + + newgridcell = .false. + + end if + end if + end do + end do + end do + + if(screen > 2) then + print *, "TSMP-PDAF mype(w)=", mype_world, & + ": init_dim_obs_LST: dim_obs_p=", dim_obs_p + end if + + ! Dimension of full observation vector + ! ------------------------------------ + + ! Gather and check PE-local observation dimensions + call mpi_allreduce(dim_obs_p, sum_dim_obs_p, 1, MPI_INTEGER, MPI_SUM, & + comm_filter, ierror) + + ! Check sum of dimensions of PE-local observation vectors against + ! dimension of full observation vector + if(.not. sum_dim_obs_p == dim_obs) then + print *, "TSMP-PDAF mype(w)=", mype_world, & + ": ERROR Sum of PE-local observation dimensions" + print *, "sum_dim_obs_p=", sum_dim_obs_p + print *, "dim_obs=", dim_obs + call abort_parallel() + end if + + ! Gather PE-local observation dimensions and displacements in arrays + ! ---------------------------------------------------------------- + + ! Allocate array of PE-local observation dimensions + IF (ALLOCATED(local_dims_obs)) DEALLOCATE(local_dims_obs) + ALLOCATE(local_dims_obs(npes_filter)) + + ! Gather array of PE-local observation dimensions + call mpi_allgather(dim_obs_p, 1, MPI_INTEGER, local_dims_obs, 1, MPI_INTEGER, & + comm_filter, ierror) + + ! Allocate observation displacement array local_disp_obs + IF (ALLOCATED(local_disp_obs)) DEALLOCATE(local_disp_obs) + ALLOCATE(local_disp_obs(npes_filter)) + + ! Set observation displacement array local_disp_obs + local_disp_obs(1) = 0 + do i = 2, npes_filter + local_disp_obs(i) = local_disp_obs(i-1) + local_dims_obs(i-1) + end do + + if(mype_filter==0 .and. screen > 2) then + print *, "TSMP-PDAF mype(w)=", mype_world, & + ": init_dim_obs_LST: local_disp_obs=", local_disp_obs + end if + + ! Write index mapping obs_pdaf2nc / obs_nc2pdaf (used for debug output) + if(allocated(obs_pdaf2nc)) deallocate(obs_pdaf2nc) + allocate(obs_pdaf2nc(dim_obs)) + obs_pdaf2nc = 0 + if(allocated(obs_nc2pdaf)) deallocate(obs_nc2pdaf) + allocate(obs_nc2pdaf(dim_obs)) + obs_nc2pdaf = 0 + + cnt = 1 + do i = 1, dim_obs + obs_snapped = .true. + do g = begg,endg + newgridcell = .true. + do p = begp,endp + pg = patch%gridcell(p) + if(pg == g) then + if(newgridcell) then + + deltax = abs(lon(g)-lon_obs(i)) + if (deltax > 180.0) then + deltax = 360.0 - deltax + end if + deltay = abs(lat(g)-lat_obs(i)) + + if((deltax<=dr_obs(1)).and.(deltay<=dr_obs(2))) then + if(state_clm2pdaf_p(p,1)==ispval) then + obs_snapped = .false. + cycle + end if + obs_pdaf2nc(local_disp_obs(mype_filter+1)+cnt) = i + obs_nc2pdaf(i) = local_disp_obs(mype_filter+1)+cnt + cnt = cnt + 1 + obs_snapped = .true. + end if + + newgridcell = .false. + + end if + end if + end do + end do + + if(.not. obs_snapped) then + print *, "TSMP-PDAF mype(w)=", mype_world, & + ": ERROR observations exist at non-active gridcells." + print *, "Observation-index in NetCDF-file: i=", i + call abort_parallel() + end if + end do + + ! Gather obs_pdaf2nc / obs_nc2pdaf across all PEs via summation. + ! Each PE only wrote to its own slice (indices local_disp_obs(mype+1)+1 .. + ! local_disp_obs(mype+1)+dim_obs_p), leaving all other entries zero, so + ! MPI_SUM is equivalent to a gather without a dedicated gather buffer. + call mpi_allreduce(MPI_IN_PLACE,obs_pdaf2nc,dim_obs,MPI_INTEGER,MPI_SUM,comm_filter,ierror) + call mpi_allreduce(MPI_IN_PLACE,obs_nc2pdaf,dim_obs,MPI_INTEGER,MPI_SUM,comm_filter,ierror) + + if(mype_filter==0 .and. screen > 2) then + print *, "TSMP-PDAF mype(w)=", mype_world, & + ": init_dim_obs_LST: obs_pdaf2nc=", obs_pdaf2nc + end if + + ! Write process-local observation arrays + ! -------------------------------------- + + ! Use module-local obs_index_p_LST to avoid conflict with global obs_index_p + IF (ALLOCATED(obs_index_p_LST)) DEALLOCATE(obs_index_p_LST) + ALLOCATE(obs_index_p_LST(dim_obs_p)) + IF (ALLOCATED(obs_p)) DEALLOCATE(obs_p) + ALLOCATE(obs_p(dim_obs_p)) + + ! Initialize OMI arrays + IF (ALLOCATED(ivar_obs_p)) DEALLOCATE(ivar_obs_p) + ALLOCATE(ivar_obs_p(dim_obs_p)) + IF (ALLOCATED(ocoord_p)) DEALLOCATE(ocoord_p) + ALLOCATE(ocoord_p(2, dim_obs_p)) + + cnt = 1 + + do i = 1, dim_obs + + do g = begg,endg + newgridcell = .true. + + do p = begp,endp + + pg = patch%gridcell(p) + + if(pg == g) then + + if(newgridcell) then + ! Sets first patch/column in a gridcell. TODO: Make + ! patch / column information part of the observation + ! file + + deltax = abs(lon(g)-lon_obs(i)) + if (deltax > 180.0) then + deltax = 360.0 - deltax + end if + deltay = abs(lat(g)-lat_obs(i)) + + if((deltax<=dr_obs(1)).and.(deltay<=dr_obs(2))) then + + ! Convert observation coordinates to radians for haversine distance + if(thisobs%disttype==3) then + ocoord_p(1,cnt) = lon_obs(i) * pi / 180.0 + ocoord_p(2,cnt) = lat_obs(i) * pi / 180.0 + else + ocoord_p(1,cnt) = lon_obs(i) + ocoord_p(2,cnt) = lat_obs(i) + end if + + ! Set state vector index for this patch and save it in + ! the observation index array. + ! + ! LST uses patch (not column), variable 1 = TSKIN. + obs_index_p_LST(cnt) = state_clm2pdaf_p(p,1) + + obs_p(cnt) = obs_g(i) + if(multierr==1) ivar_obs_p(cnt) = 1.0/(obserr(i)*obserr(i)) + if(multierr==0) ivar_obs_p(cnt) = 1.0/(rms_obs_LST*rms_obs_LST) + cnt = cnt + 1 + + end if + + newgridcell = .false. + + end if + + end if + + end do + end do + + end do + +#ifdef PDAF_DEBUG + IF (da_print_obs_index > 0) THEN + WRITE(fn, "(a,i5.5,a,i5.5,a)") "obs_index_p_LST_", mype_world, ".", step, ".txt" + OPEN(unit=72, file=fn, action="write") + DO i = 1, dim_obs_p + WRITE (72,"(i10)") obs_index_p_LST(i) + END DO + CLOSE(72) + END IF +#endif + + ! **************************************** + ! *** Gather global observation arrays *** + ! **************************************** + + CALL PDAFomi_gather_obs(thisobs, dim_obs_p, obs_p, ivar_obs_p, ocoord_p, & + thisobs%ncoord, cradius_LST, dim_obs) + + ! ******************** + ! *** Finishing up *** + ! ******************** + + DEALLOCATE(obs_g) + DEALLOCATE(obs_p) + DEALLOCATE(ocoord_p) + DEALLOCATE(ivar_obs_p) + +END SUBROUTINE init_dim_obs_LST + + + +!------------------------------------------------------------------------------- +!> Implementation of observation operator for LST +!! +!! Applies the observation operator H for LST. Three modes: +!! clmupdate_T==1: Kustas (2009) Eq.(7) +!! LST = (exp(-tau/2)*TG^4 + (1-exp(-tau/2))*TV^4)^0.25 +!! clmupdate_T==2,3: direct TSKIN mapping, LST = state_p(obs_index) +!! +SUBROUTINE obs_op_LST(dim_p, dim_obs, state_p, ostate) + + use enkf_clm_mod, only: clm_paramarr, clm_varsize, clmupdate_T + + use PDAFomi_obs_f, only: PDAFomi_gather_obsstate + + IMPLICIT NONE + + INTEGER, INTENT(in) :: dim_p + INTEGER, INTENT(in) :: dim_obs + REAL, INTENT(in) :: state_p(dim_p) + REAL, INTENT(inout) :: ostate(dim_obs) + + real, allocatable :: ostate_p(:) + integer :: i + + IF (thisobs%dim_obs_p > 0) THEN + if(allocated(ostate_p)) deallocate(ostate_p) + ALLOCATE(ostate_p(thisobs%dim_obs_p)) + ELSE + if(allocated(ostate_p)) deallocate(ostate_p) + ALLOCATE(ostate_p(1)) + END IF + + if(clmupdate_T == 1) then + + ! Kustas (2009) Eq.(7): combined ground/vegetation brightness temperature + ! tau = clm_paramarr(idx): vegetation optical depth (LAI-based) + ! TG = state_p(obs_index_p_LST(i)): ground temperature + ! TV = state_p(clm_varsize + obs_index_p_LST(i)): vegetation temperature + DO i = 1, thisobs%dim_obs_p + ostate_p(i) = & + ( exp(-0.5*clm_paramarr(obs_index_p_LST(i))) & + * state_p(obs_index_p_LST(i))**4 & + + (1.0 - exp(-0.5*clm_paramarr(obs_index_p_LST(i)))) & + * state_p(clm_varsize + obs_index_p_LST(i))**4 )**0.25 + END DO + + else + + ! clmupdate_T==2 or 3: simulated LST equals TSKIN directly + DO i = 1, thisobs%dim_obs_p + ostate_p(i) = state_p(obs_index_p_LST(i)) + END DO + + end if + + CALL PDAFomi_gather_obsstate(thisobs, ostate_p, ostate) + + deallocate(ostate_p) + +END SUBROUTINE obs_op_LST + + + +!------------------------------------------------------------------------------- +!> Initialize local information on the module-type observation +!! +!! Called during the loop over all local analysis domains. +!! Uses patch-based localization domain coordinates. +!! +SUBROUTINE init_dim_obs_l_LST(domain_p, step, dim_obs, dim_obs_l) + + USE PDAFomi, ONLY: PDAFomi_init_dim_obs_l, pi + + USE mod_assimilation, & + ONLY: cradius_LST, locweight, sradius_LST, screen + + USE enkf_clm_mod, ONLY: state_loc2clm_p_p + + use shr_kind_mod, only: r8 => shr_kind_r8 + + use decompMod , only : get_proc_bounds + + USE GridcellType, ONLY: grc + USE PatchType, ONLY: patch + + use clm_varcon, only: spval + + use mod_parallel_pdaf, ONLY: mype_world + + IMPLICIT NONE + + INTEGER, INTENT(in) :: domain_p + INTEGER, INTENT(in) :: step + INTEGER, INTENT(in) :: dim_obs + INTEGER, INTENT(inout) :: dim_obs_l + + REAL :: coords_l(2) + + real(r8), pointer :: lon(:) + real(r8), pointer :: lat(:) + + integer :: pg_l + + lon => grc%londeg + lat => grc%latdeg + + ! ********************************************** + ! *** Initialize local observation dimension *** + ! ********************************************** + + if(thisobs%infile==1) then + + ! Get gridcell of local analysis domain via patch index + pg_l = patch%gridcell(state_loc2clm_p_p(domain_p)) + + if(lon(pg_l) > 180) then + coords_l(1) = lon(pg_l) - 360.0 + else + coords_l(1) = lon(pg_l) + end if + coords_l(2) = lat(pg_l) + + if(thisobs%disttype==3) then + coords_l(1) = coords_l(1) * pi / 180.0 + coords_l(2) = coords_l(2) * pi / 180.0 + end if + + else + + coords_l(1) = spval + coords_l(2) = spval + + end if + + ! For disttype=3, cradius and sradius are in km; multiply by 1000 for meters + if(thisobs%disttype==3) then + CALL PDAFomi_init_dim_obs_l(thisobs_l, thisobs, coords_l, & + locweight, cradius_LST*1000.0, sradius_LST*1000.0, dim_obs_l) + else + CALL PDAFomi_init_dim_obs_l(thisobs_l, thisobs, coords_l, & + locweight, cradius_LST, sradius_LST, dim_obs_l) + end if + +END SUBROUTINE init_dim_obs_l_LST + + + +!------------------------------------------------------------------------------- +!> Perform covariance localization for local EnKF on the module-type observation +!! +!! Uses patch-based state vector coordinates. +!! +SUBROUTINE localize_covar_LST(dim_p, dim_obs, HP_p, HPH, coords_p) + + USE PDAFomi, ONLY: PDAFomi_localize_covar + + USE mod_assimilation, & + ONLY: cradius_LST, locweight, sradius_LST + + use enkf_clm_mod, only: state_pdaf2clm_p_p + + use shr_kind_mod, only: r8 => shr_kind_r8 + + USE GridcellType, ONLY: grc + USE PatchType, ONLY: patch + + IMPLICIT NONE + + INTEGER, INTENT(in) :: dim_p + INTEGER, INTENT(in) :: dim_obs + REAL, INTENT(inout) :: HP_p(dim_obs, dim_p) + REAL, INTENT(inout) :: HPH(dim_obs, dim_obs) + REAL, INTENT(inout) :: coords_p(:,:) + + integer :: i, pg_i + + real(r8), pointer :: lon(:) + real(r8), pointer :: lat(:) + + lon => grc%londeg + lat => grc%latdeg + + ! Assign geographic coordinates from patch gridcells to state vector elements + do i = 1,dim_p + pg_i = patch%gridcell(state_pdaf2clm_p_p(i)) + if(lon(pg_i) > 180) then + coords_p(1,i) = lon(pg_i) - 360.0 + else + coords_p(1,i) = lon(pg_i) + end if + coords_p(2,i) = lat(pg_i) + end do + + CALL PDAFomi_localize_covar(thisobs, dim_p, locweight, cradius_LST, sradius_LST, & + coords_p, HP_p, HPH) + +END SUBROUTINE localize_covar_LST + + +subroutine add_obs_err_LST(step, dim_obs, C) + + USE mod_parallel_pdaf, ONLY: npes_filter + + use PDAFomi, only: obsdims + + implicit none + INTEGER, INTENT(in) :: step + INTEGER, INTENT(in) :: dim_obs + REAL, INTENT(inout) :: C(dim_obs,dim_obs) + + integer :: i, pe, cnt + INTEGER, ALLOCATABLE :: id_start(:) + INTEGER, ALLOCATABLE :: id_end(:) + + ALLOCATE(id_start(npes_filter), id_end(npes_filter)) + + pe = 1 + id_start(1) = 1 + IF (thisobs%obsid>1) id_start(1) = id_start(1) + sum(obsdims(1, 1:thisobs%obsid-1)) + id_end(1) = id_start(1) + obsdims(1,thisobs%obsid) - 1 + DO pe = 2, npes_filter + id_start(pe) = id_start(pe-1) + SUM(obsdims(pe-1,thisobs%obsid:)) + IF (thisobs%obsid>1) id_start(pe) = id_start(pe) + sum(obsdims(pe,1:thisobs%obsid-1)) + id_end(pe) = id_start(pe) + obsdims(pe,thisobs%obsid) - 1 + END DO + + cnt = 1 + DO pe = 1, npes_filter + DO i = id_start(pe), id_end(pe) + C(i,i) = C(i,i) + 1.0/thisobs%ivar_obs_f(cnt) + cnt = cnt + 1 + end do + end do + + DEALLOCATE(id_start, id_end) + +end subroutine add_obs_err_LST + + +subroutine init_obscovar_LST(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) + + USE mod_parallel_pdaf, ONLY: npes_filter + + use PDAFomi, only: obsdims, map_obs_id + + implicit none + INTEGER, INTENT(in) :: step + INTEGER, INTENT(in) :: dim_obs + INTEGER, INTENT(in) :: dim_obs_p + REAL, INTENT(inout) :: covar(dim_obs, dim_obs) + REAL, INTENT(in) :: m_state_p(dim_obs_p) + LOGICAL, INTENT(inout) :: isdiag + + integer :: i, pe, cnt + INTEGER, ALLOCATABLE :: id_start(:) + INTEGER, ALLOCATABLE :: id_end(:) + + ALLOCATE(id_start(npes_filter), id_end(npes_filter)) + + pe = 1 + id_start(1) = 1 + IF (thisobs%obsid>1) id_start(1) = id_start(1) + sum(obsdims(1, 1:thisobs%obsid-1)) + id_end(1) = id_start(1) + obsdims(1,thisobs%obsid) - 1 + DO pe = 2, npes_filter + id_start(pe) = id_start(pe-1) + SUM(obsdims(pe-1,thisobs%obsid:)) + IF (thisobs%obsid>1) id_start(pe) = id_start(pe) + sum(obsdims(pe,1:thisobs%obsid-1)) + id_end(pe) = id_start(pe) + obsdims(pe,thisobs%obsid) - 1 + END DO + + cnt = 1 + IF (thisobs%obsid-1 > 0) cnt = cnt + SUM(obsdims(:,1:thisobs%obsid-1)) + DO pe = 1, npes_filter + DO i = id_start(pe), id_end(pe) + map_obs_id(i) = cnt + cnt = cnt + 1 + END DO + END DO + + cnt = 1 + DO pe = 1, npes_filter + DO i = id_start(pe), id_end(pe) + covar(i, i) = covar(i, i) + 1.0/thisobs%ivar_obs_f(cnt) + cnt = cnt + 1 + ENDDO + ENDDO + + isdiag = .TRUE. + + DEALLOCATE(id_start, id_end) + +end subroutine init_obscovar_LST + + +subroutine prodRinvA_LST(step, dim_obs_p, rank, obs_p, A_p, C_p) + + INTEGER, INTENT(in) :: step + INTEGER, INTENT(in) :: dim_obs_p + INTEGER, INTENT(in) :: rank + REAL, INTENT(in) :: obs_p(dim_obs_p) + REAL, INTENT(in) :: A_p(dim_obs_p,rank) + REAL, INTENT(inout) :: C_p(dim_obs_p,rank) + + INTEGER :: i, j, off + + off = thisobs%off_obs_f + + do j = 1, rank + do i = 1, thisobs%dim_obs_f + C_p(i+off, j) = thisobs%ivar_obs_f(i) * A_p(i+off, j) + END DO + end do + +end subroutine prodRinvA_LST + + +subroutine prodRinvA_l_LST(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) + + use shr_kind_mod, only: r8 => shr_kind_r8 + USE mod_assimilation, ONLY: cradius_LST, locweight, sradius_LST + use pdafomi, only: PDAFomi_observation_localization_weights + + implicit none + + INTEGER, INTENT(in) :: domain_p + INTEGER, INTENT(in) :: step + INTEGER, INTENT(in) :: dim_obs_l + INTEGER, INTENT(in) :: rank + REAL, INTENT(in) :: obs_l(dim_obs_l) + REAL, INTENT(inout) :: A_l(dim_obs_l, rank) + REAL, INTENT(out) :: C_l(dim_obs_l, rank) + + INTEGER :: verbose + INTEGER, SAVE :: domain_save = -1 + REAL, ALLOCATABLE :: weight(:) + INTEGER :: i, j, off, idummy + + off = thisobs_l%off_obs_l + idummy = dim_obs_l + + IF ((domain_p <= domain_save .OR. domain_save < 0) .AND. mype_filter==0) THEN + verbose = 1 + ELSE + verbose = 0 + END IF + domain_save = domain_p + + IF (verbose == 1) THEN + WRITE (*, '(8x, a, f12.3)') & + '--- Use global rms for LST observations of ', rms_obs_LST + WRITE (*, '(8x, a, 1x)') & + '--- Domain localization' + WRITE (*, '(12x, a, 1x, f12.2)') & + '--- Local influence radius', cradius_LST + IF (locweight > 0) THEN + WRITE (*, '(12x, a)') & + '--- Use distance-dependent weight for observation errors' + END IF + ENDIF + + ALLOCATE(weight(thisobs_l%dim_obs_l)) + call PDAFomi_observation_localization_weights(thisobs_l, thisobs, rank, A_l, & + weight, verbose) + + do j = 1, rank + do i = 1, thisobs_l%dim_obs_l + C_l(i+off,j) = thisobs_l%ivar_obs_l(i) * weight(i) * A_l(i+off, j) + end do + end do + + deallocate(weight) + +end subroutine prodRinvA_l_LST + + +subroutine deallocate_obs_LST() + + USE PDAFomi, ONLY: PDAFomi_deallocate_obs + USE PDAFomi_obs_l, ONLY: obs_l_all, firstobs + + implicit none + + if(mype_filter==0) then + WRITE (*,*) 'Deallocating observations type LST' + end if + call PDAFomi_deallocate_obs(thisobs) + + if(allocated(thisobs_l%id_obs_l)) deallocate(thisobs_l%id_obs_l) + if(allocated(thisobs_l%ivar_obs_l)) deallocate(thisobs_l%ivar_obs_l) + if(allocated(thisobs_l%distance_l)) deallocate(thisobs_l%distance_l) + if(allocated(thisobs_l%cradius_l)) deallocate(thisobs_l%cradius_l) + if(allocated(thisobs_l%sradius_l)) deallocate(thisobs_l%sradius_l) + if(allocated(thisobs_l%dist_l_v)) deallocate(thisobs_l%dist_l_v) + if(allocated(thisobs_l%cradius)) deallocate(thisobs_l%cradius) + if(allocated(thisobs_l%sradius)) deallocate(thisobs_l%sradius) + + if(allocated(obs_l_all)) deallocate(obs_l_all) + + firstobs = 0 + + if(allocated(obs_index_p_LST)) deallocate(obs_index_p_LST) + +end subroutine deallocate_obs_LST + + +END MODULE obs_LST_pdafomi +#endif diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 1ce5db81..4f2f9577 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -133,13 +133,14 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) lpointobs = .false. + ! Equation for LST computation from Kustas2009, Eq(7) + ! http://dx.doi.org/10.1016/j.agrformet.2009.05.016 + ! + ! Comment: Fractional vegetation cover (Eq(8) from Kustas2009) + ! currently implemented with simplified settings: Vegetation + ! clumping parameter `Omega=1`; radiometer view angle `phi=0` + DO i = 1, dim_obs_p - ! Equation for LST computation from Kustas2009, Eq(7) - ! http://dx.doi.org/10.1016/j.agrformet.2009.05.016 - ! - ! Comment: Fractional vegetation cover (Eq(8) from Kustas2009) - ! currently implemented with simplified settings: Vegetation - ! clumping parameter `Omega=1`; radiometer view angle `phi=0` m_state_p(i) & = (exp(-0.5*clm_paramarr(obs_index_p(i))) & *state_p(obs_index_p(i))**4 & @@ -151,6 +152,51 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) ! write(*,*) 'model LST', m_state_p(:) ! write(*,*) 'TG', state_p(obs_index_p(:)) ! write(*,*) 'TV', state_p(clm_varsize+obs_index_p(:)) + +endif + +if (clmupdate_T==2) then + + lpointobs = .false. + + DO i = 1, dim_obs_p + ! first implementation: simulated LST equals TSKIN + m_state_p(i) = state_p(obs_index_p(i)) + END DO + +endif + +if (clmupdate_T==3) then + + lpointobs = .false. + + DO i = 1, dim_obs_p + ! first implementation: simulated LST equals TSKIN + m_state_p(i) = state_p(obs_index_p(i)) + END DO + +endif + +if (clmupdate_T==4) then + + lpointobs = .false. + + DO i = 1, dim_obs_p + ! first implementation: simulated LST equals TSKIN + m_state_p(i) = state_p(obs_index_p(i)) + END DO + +endif + +if (clmupdate_T==5) then + + lpointobs = .false. + + DO i = 1, dim_obs_p + ! first implementation: simulated LST equals TSKIN + m_state_p(i) = state_p(obs_index_p(i)) + END DO + endif #endif diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index 7462a6f1..8e3a3f3d 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -61,6 +61,7 @@ module enkf_clm_mod integer,allocatable :: state_pdaf2clm_c_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) @@ -216,6 +217,14 @@ subroutine set_clm_statevec(tstartcycle, mype) CLOSE(71) END IF + IF(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: Debug output of CLM t_soisno, first layer + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF + END IF #endif @@ -475,32 +484,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end do -#ifdef PDAF_DEBUG - IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN - - IF(clmupdate_swc.NE.0) THEN - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn3, action="write") - WRITE (71,"(es22.15)") h2osoi_liq(:,:) - CLOSE(71) - - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn4, action="write") - WRITE (71,"(es22.15)") h2osoi_ice(:,:) - CLOSE(71) - - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn2, action="write") - WRITE (71,"(es22.15)") swc(:,:) - CLOSE(71) - END IF - - END IF -#endif - endif !hcp: TG, TV @@ -536,6 +519,39 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call clm_texture_to_parameters endif +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN + + IF(clmupdate_swc.NE.0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn3, action="write") + WRITE (71,"(es22.15)") h2osoi_liq(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn4, action="write") + WRITE (71,"(es22.15)") h2osoi_ice(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") swc(:,:) + CLOSE(71) + END IF + + IF(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + END IF + + END IF +#endif + end subroutine update_clm subroutine clm_correct_texture() diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h index a734c335..f3a532f3 100755 --- a/interface/model/common/enkf.h +++ b/interface/model/common/enkf.h @@ -101,6 +101,8 @@ GLOBAL int clmstatevec_max_layer; GLOBAL int clmt_printensemble; GLOBAL int clmwatmin_switch; GLOBAL int clmswc_mask_snow; +GLOBAL int clmT_mask_snow; +GLOBAL int clmincrement_type; GLOBAL int dtmult_cosmo; GLOBAL int pf_olfmasking; GLOBAL int pf_olfmasking_param; @@ -149,4 +151,7 @@ GLOBAL double dampfac_state_time_dependent; GLOBAL double dampfac_param_time_dependent; GLOBAL double da_crns_depth_tol; GLOBAL double clmcrns_bd; +GLOBAL double clmT_mask_T; GLOBAL double max_inc; +GLOBAL double clmT_mask_snow_depth; +GLOBAL double clmT_max_increment; diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c index 27b26c46..ed75839f 100755 --- a/interface/model/common/read_enkfpar.c +++ b/interface/model/common/read_enkfpar.c @@ -89,7 +89,12 @@ void read_enkfpar(char *parname) clmt_printensemble = iniparser_getint(pardict,"CLM:t_printensemble",-2); clmwatmin_switch = iniparser_getint(pardict,"CLM:watmin_switch",0); clmswc_mask_snow = iniparser_getint(pardict,"CLM:swc_mask_snow",0); + clmT_mask_snow = iniparser_getint(pardict,"CLM:T_mask_snow",0); + clmincrement_type = iniparser_getint(pardict,"CLM:increment_type",0); + clmT_mask_T = iniparser_getdouble(pardict,"CLM:T_mask_T",0.0); clmupdate_tws = iniparser_getint(pardict,"CLM:update_tws",0); + clmT_mask_snow_depth = iniparser_getdouble(pardict,"CLM:T_mask_snow_depth",0.001); + clmT_max_increment = iniparser_getdouble(pardict,"CLM:T_max_increment",5.0); /* get settings for COSMO */ nproccosmo = iniparser_getint(pardict,"COSMO:nprocs",0); diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index bc2e4c59..48e7adcf 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -46,8 +46,10 @@ module enkf_clm_mod 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) @@ -104,7 +106,12 @@ module enkf_clm_mod integer(c_int),bind(C,name="clmt_printensemble") :: clmt_printensemble integer(c_int),bind(C,name="clmwatmin_switch") :: clmwatmin_switch integer(c_int),bind(C,name="clmswc_mask_snow") :: clmswc_mask_snow - real(c_double),bind(C,name="clmcrns_bd") :: clmcrns_bd + integer(c_int),bind(C,name="clmT_mask_snow") :: clmT_mask_snow + integer(c_int),bind(C,name="clmincrement_type") :: clmincrement_type + real(c_double),bind(C,name="clmT_mask_T") :: clmT_mask_T + real(c_double),bind(C,name="clmT_mask_snow_depth") :: clmT_mask_snow_depth + real(c_double),bind(C,name="clmcrns_bd") :: clmcrns_bd + real(c_double),bind(C,name="clmT_max_increment") :: clmT_max_increment integer :: nstep ! time step index real(r8) :: dtime ! time step increment (sec) @@ -159,6 +166,11 @@ subroutine define_clm_statevec(mype) clm_begp = begp clm_endp = endp +#ifdef PDAF_DEBUG + WRITE(*, '(a,x,a,i5,x,a,i10)') "TSMP-PDAF-debug", "mype(w)=", mype, & + "define_clm_statevec entry: clm_statevecsize=", clm_statevecsize +#endif + clm_statevecsize = 0 clm_varsize = 0 @@ -192,8 +204,8 @@ subroutine define_clm_statevec(mype) end if !hcp LST DA - if(clmupdate_T==1) then - error stop "Not implemented: clmupdate_T.eq.1" + if(clmupdate_T/=0) then + call define_clm_statevec_T(mype) end if !end hcp @@ -248,16 +260,18 @@ subroutine define_clm_statevec(mype) ! Allocate statevector-duplicate for saving original column mean ! values used in computing increments during updating the state - ! vector in column-mean-mode. + ! vector in column-mean-mode (SWC) or gridcell-mean-mode (T). 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_T==2 & + .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5 & + .or. clmupdate_tws/=0 ) then allocate(clm_statevec_orig(clm_statevecsize)) end if !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp - if ((clmupdate_T/=0)) then !hcp - error stop "Not implemented clmupdate_T.NE.0" + if ((clmupdate_T==1)) then !hcp + allocate(clm_paramarr(clm_paramsize)) end if if (allocated(gridcell_state)) deallocate(gridcell_state) @@ -461,6 +475,461 @@ subroutine define_clm_statevec_swc(mype) end subroutine define_clm_statevec_swc + subroutine define_clm_statevec_T(mype) + use decompMod , only : get_proc_bounds + use clm_varpar , only : nlevgrnd + use clm_varcon , only : ispval + use PatchType , only : patch + + implicit none + + integer,intent(in) :: mype + + integer :: j + integer :: jj + integer :: lev + integer :: p + integer :: cc + integer :: n_lev_T ! effective number of soil layers in T state vector + + integer :: begp, endp ! per-proc beginning and ending pft indices + 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 + + + call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp) + + n_lev_T = min(nlevgrnd, clmstatevec_max_layer) + +#ifdef PDAF_DEBUG + WRITE(*,"(a,i5,a,i10,a,i10,a,i10,a,i10,a,i10,a,i10,a,i10,a,i10,a)") & + "TSMP-PDAF mype(w)=", mype, " define_clm_statevec, CLM5-bounds (g,l,c,p)----",& + begg,",",endg,",",begl,",",endl,",",begc,",",endc,",",begp,",",endp," -------" +#endif + + clm_begg = begg + clm_endg = endg + clm_begc = begc + clm_endc = endc + clm_begp = begp + clm_endp = endp + + 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 + + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + ! Default: inactive + do cc=1,clm_statevecsize + state_pdaf2clm_p_p(cc) = ispval + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TG + state_pdaf2clm_c_p(cc) = patch%column(p) !TG + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + end do + + endif + !end hcp + + if(clmupdate_T==2) then + + ! 1) PATCH/GRC: CLM->PDAF + ! Allocate with full dimension for all variables/layers + ! + ! Here, the second index of STATE_CLM2PDAF_P is NOT simply the + ! layer index of CLM, but rather a variable index of the state + ! vector, in order tomake the index mapping 1:1. + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1:(2+n_lev_T))) + ! ^ + ! dimension layout: + ! 1: TSKIN + ! 2:(1+n_lev_T): TSOIL layers + ! (2+n_lev_T): TVEG + do lev=1,(2+n_lev_T) + do p=begp,endp + ! Default: inactive + state_clm2pdaf_p(p,lev) = ispval + end do + end do + + do p=clm_begp,clm_endp + ! All patches in a gridcell are assigned the index of + ! the gridcell-average in the state vector + cc = (patch%gridcell(p) - clm_begg + 1) + + ! TSKIN (variable index 1) + state_clm2pdaf_p(p,1) = cc + + ! TSOIL layers (variable indices 2 to 1+n_lev_T) + do lev=1,n_lev_T + state_clm2pdaf_p(p, 1+lev) = cc + lev*(clm_endg - clm_begg + 1) + end do + + ! TVEG (variable index 2+n_lev_T) + state_clm2pdaf_p(p, 2+n_lev_T) = cc + (1+n_lev_T)*(clm_endg - clm_begg + 1) + end do + + ! 2) PATCH/GRC: STATEVECSIZE + clm_varsize = endg-begg+1 + clm_statevecsize = (2 + n_lev_T)* (endg-begg+1) !TSKIN, then n_lev_T times TSOIL and TV + + ! 3) PATCH/GRC: PDAF->CLM + ! + ! Inverse mapping from state vector index to CLM indices. + ! STATE_PDAF2CLM_*_P is the inverse of STATE_CLM2PDAF_P: + ! state_pdaf2clm_p_p(cc) = p (patch index) + ! state_pdaf2clm_c_p(cc) = c (column index) + ! state_pdaf2clm_j_p(cc) = j (variable index, NOT CLM level) + ! + ! Variable index j layout (same as STATE_CLM2PDAF_P second dimension): + ! j = 1: TSKIN + ! j = 2:(1+n_lev_T): TSOIL layers (CLM level = j - 1) + ! j = (2+n_lev_T): TVEG + ! + ! NOTE: To get the CLM level from j for TSOIL, use: clm_level = j - 1 + ! + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + ! Default: inactive + do cc=1,clm_statevecsize + state_pdaf2clm_p_p(cc) = ispval + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do + + cc = 0 + + do cc=1,clm_statevecsize + + do p=clm_begp,clm_endp + if(state_clm2pdaf_p(p, 1) == cc) then + ! Set indices for the three temperature variables and then + ! exit loops + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 ! variable index 1: TSKIN + do lev=1,n_lev_T + ! variable index 2 to 1+n_lev_T: TSOIL layers + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p !TSOIL + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc + lev*clm_varsize) = 1 + lev ! variable index TSOIL + end do + state_pdaf2clm_p_p(cc+(1+n_lev_T)*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+(1+n_lev_T)*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+(1+n_lev_T)*clm_varsize) = 2 + n_lev_T ! variable index: TVEG + exit + end if + end do + + end do + + endif + + if(clmupdate_T==3) then + + ! 1) PATCH/GRC: CLM->PDAF + ! Allocate with full dimension for all variables/layers + ! + ! Here, the second index of STATE_CLM2PDAF_P is NOT simply the + ! layer index of CLM, but rather a variable index of the state + ! vector, in order tomake the index mapping 1:1. + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1:(3+n_lev_T))) + ! ^ + ! dimension layout: + ! 1: TSKIN + ! 2:(1+n_lev_T): TSOIL layers + ! (2+n_lev_T): TVEG + ! (3+n_lev_T): TGRND + do lev=1,(3+n_lev_T) + do p=begp,endp + ! Default: inactive + state_clm2pdaf_p(p,lev) = ispval + end do + end do + + do p=clm_begp,clm_endp + ! All patches in a gridcell are assigned the index of + ! the gridcell-average in the state vector + cc = (patch%gridcell(p) - clm_begg + 1) + + ! TSKIN (variable index 1) + state_clm2pdaf_p(p,1) = cc + + ! TSOIL layers (variable indices 2 to 1+n_lev_T) + do lev=1,n_lev_T + state_clm2pdaf_p(p, 1+lev) = cc + lev*(clm_endg - clm_begg + 1) + end do + + ! TVEG (variable index 2+n_lev_T) + state_clm2pdaf_p(p, 2+n_lev_T) = cc + (1+n_lev_T)*(clm_endg - clm_begg + 1) + + ! TGRND (variable index 3+n_lev_T) + state_clm2pdaf_p(p, 3+n_lev_T) = cc + (2+n_lev_T)*(clm_endg - clm_begg + 1) + end do + + ! 2) PATCH/GRC: STATEVECSIZE + clm_varsize = endg-begg+1 + clm_statevecsize = (3 + n_lev_T)* (endg-begg+1) !TSKIN, then n_lev_T times TSOIL, TV and TGRND + + ! 3) PATCH/GRC: PDAF->CLM + ! + ! Inverse mapping from state vector index to CLM indices. + ! STATE_PDAF2CLM_*_P is the inverse of STATE_CLM2PDAF_P: + ! state_pdaf2clm_p_p(cc) = p (patch index) + ! state_pdaf2clm_c_p(cc) = c (column index) + ! state_pdaf2clm_j_p(cc) = j (variable index, NOT CLM level) + ! + ! Variable index j layout (same as STATE_CLM2PDAF_P second dimension): + ! j = 1: TSKIN + ! j = 2:(1+n_lev_T): TSOIL layers (CLM level = j - 1) + ! j = (2+n_lev_T): TVEG + ! j = (3+n_lev_T): TGRND + ! + ! NOTE: To get the CLM level from j for TSOIL, use: clm_level = j - 1 + ! + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + ! Default: inactive + do cc=1,clm_statevecsize + state_pdaf2clm_p_p(cc) = ispval + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do + + cc = 0 + + do cc=1,clm_statevecsize + + do p=clm_begp,clm_endp + if(state_clm2pdaf_p(p, 1) == cc) then + ! Set indices for the four temperature variables and then + ! exit loops + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 ! variable index 1: TSKIN + do lev=1,n_lev_T + ! variable index 2 to 1+n_lev_T: TSOIL layers + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p !TSOIL + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc + lev*clm_varsize) = 1 + lev ! variable index TSOIL + end do + state_pdaf2clm_p_p(cc+(1+n_lev_T)*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+(1+n_lev_T)*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+(1+n_lev_T)*clm_varsize) = 2 + n_lev_T ! variable index: TVEG + state_pdaf2clm_p_p(cc+(2+n_lev_T)*clm_varsize) = p !TGRND + state_pdaf2clm_c_p(cc+(2+n_lev_T)*clm_varsize) = patch%column(p) !TGRND + state_pdaf2clm_j_p(cc+(2+n_lev_T)*clm_varsize) = 3 + n_lev_T ! variable index: TGRND + exit + end if + end do + + end do + + endif + + if(clmupdate_T==4) then + + ! clmupdate_T==4: like ==2 but with T_H2OSFC added to state vector. + ! State vector layout: + ! 1: TSKIN + ! 2:(1+n_lev_T): TSOIL layers + ! (2+n_lev_T): TVEG + ! (3+n_lev_T): T_H2OSFC + + ! 1) PATCH/GRC: CLM->PDAF + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1:(3+n_lev_T))) + do lev=1,(3+n_lev_T) + do p=begp,endp + ! Default: inactive + state_clm2pdaf_p(p,lev) = ispval + end do + end do + + do p=clm_begp,clm_endp + cc = (patch%gridcell(p) - clm_begg + 1) + ! TSKIN (variable index 1) + state_clm2pdaf_p(p,1) = cc + ! TSOIL layers (variable indices 2 to 1+n_lev_T) + do lev=1,n_lev_T + state_clm2pdaf_p(p, 1+lev) = cc + lev*(clm_endg - clm_begg + 1) + end do + ! TVEG (variable index 2+n_lev_T) + state_clm2pdaf_p(p, 2+n_lev_T) = cc + (1+n_lev_T)*(clm_endg - clm_begg + 1) + ! T_H2OSFC (variable index 3+n_lev_T) + state_clm2pdaf_p(p, 3+n_lev_T) = cc + (2+n_lev_T)*(clm_endg - clm_begg + 1) + end do + + ! 2) PATCH/GRC: STATEVECSIZE + clm_varsize = endg-begg+1 + clm_statevecsize = (3 + n_lev_T)* (endg-begg+1) !TSKIN, n_lev_T*TSOIL, TVEG, T_H2OSFC + + ! 3) PATCH/GRC: PDAF->CLM + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + ! Default: inactive + do cc=1,clm_statevecsize + state_pdaf2clm_p_p(cc) = ispval + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do + + do cc=1,clm_statevecsize + do p=clm_begp,clm_endp + if(state_clm2pdaf_p(p, 1) == cc) then + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 ! variable index 1: TSKIN + do lev=1,n_lev_T + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p !TSOIL + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc + lev*clm_varsize) = 1 + lev ! variable index TSOIL + end do + state_pdaf2clm_p_p(cc+(1+n_lev_T)*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+(1+n_lev_T)*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+(1+n_lev_T)*clm_varsize) = 2 + n_lev_T ! variable index: TVEG + state_pdaf2clm_p_p(cc+(2+n_lev_T)*clm_varsize) = p !T_H2OSFC + state_pdaf2clm_c_p(cc+(2+n_lev_T)*clm_varsize) = patch%column(p) !T_H2OSFC + state_pdaf2clm_j_p(cc+(2+n_lev_T)*clm_varsize) = 3 + n_lev_T ! variable index: T_H2OSFC + exit + end if + end do + end do + + endif + + if(clmupdate_T==5) then + + ! clmupdate_T==5: like ==3 but with T_H2OSFC added to state vector. + ! State vector layout: + ! 1: TSKIN + ! 2:(1+n_lev_T): TSOIL layers + ! (2+n_lev_T): TVEG + ! (3+n_lev_T): TGRND + ! (4+n_lev_T): T_H2OSFC + + ! 1) PATCH/GRC: CLM->PDAF + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1:(4+n_lev_T))) + do lev=1,(4+n_lev_T) + do p=begp,endp + ! Default: inactive + state_clm2pdaf_p(p,lev) = ispval + end do + end do + + do p=clm_begp,clm_endp + cc = (patch%gridcell(p) - clm_begg + 1) + ! TSKIN (variable index 1) + state_clm2pdaf_p(p,1) = cc + ! TSOIL layers (variable indices 2 to 1+n_lev_T) + do lev=1,n_lev_T + state_clm2pdaf_p(p, 1+lev) = cc + lev*(clm_endg - clm_begg + 1) + end do + ! TVEG (variable index 2+n_lev_T) + state_clm2pdaf_p(p, 2+n_lev_T) = cc + (1+n_lev_T)*(clm_endg - clm_begg + 1) + ! TGRND (variable index 3+n_lev_T) + state_clm2pdaf_p(p, 3+n_lev_T) = cc + (2+n_lev_T)*(clm_endg - clm_begg + 1) + ! T_H2OSFC (variable index 4+n_lev_T) + state_clm2pdaf_p(p, 4+n_lev_T) = cc + (3+n_lev_T)*(clm_endg - clm_begg + 1) + end do + + ! 2) PATCH/GRC: STATEVECSIZE + clm_varsize = endg-begg+1 + clm_statevecsize = (4 + n_lev_T)* (endg-begg+1) !TSKIN, n_lev_T*TSOIL, TVEG, TGRND, T_H2OSFC + + ! 3) PATCH/GRC: PDAF->CLM + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + ! Default: inactive + do cc=1,clm_statevecsize + state_pdaf2clm_p_p(cc) = ispval + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do + + do cc=1,clm_statevecsize + do p=clm_begp,clm_endp + if(state_clm2pdaf_p(p, 1) == cc) then + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 ! variable index 1: TSKIN + do lev=1,n_lev_T + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p !TSOIL + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc + lev*clm_varsize) = 1 + lev ! variable index TSOIL + end do + state_pdaf2clm_p_p(cc+(1+n_lev_T)*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+(1+n_lev_T)*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+(1+n_lev_T)*clm_varsize) = 2 + n_lev_T ! variable index: TVEG + state_pdaf2clm_p_p(cc+(2+n_lev_T)*clm_varsize) = p !TGRND + state_pdaf2clm_c_p(cc+(2+n_lev_T)*clm_varsize) = patch%column(p) !TGRND + state_pdaf2clm_j_p(cc+(2+n_lev_T)*clm_varsize) = 3 + n_lev_T ! variable index: TGRND + state_pdaf2clm_p_p(cc+(3+n_lev_T)*clm_varsize) = p !T_H2OSFC + state_pdaf2clm_c_p(cc+(3+n_lev_T)*clm_varsize) = patch%column(p) !T_H2OSFC + state_pdaf2clm_j_p(cc+(3+n_lev_T)*clm_varsize) = 4 + n_lev_T ! variable index: T_H2OSFC + exit + end if + end do + end do + + endif + + if(clmupdate_T < 1 .or. clmupdate_T > 5) then + error stop "LST-DA only implemented for clmupdate_T==1 to clmupdate_T==5." + end if + + end subroutine define_clm_statevec_T + + !> @author Yorck Ewerdwalbesloh, Anne Springer !> @date 29.10.2025 !> @brief define the state vector for TWS assimilation @@ -755,8 +1224,8 @@ subroutine set_clm_statevec(tstartcycle, mype) endif !hcp LAI - if(clmupdate_T==1) then - error stop "Not implemented: clmupdate_T.eq.1" + if(clmupdate_T/=0) then + call set_clm_statevec_T(tstartcycle,mype) endif !end hcp LAI @@ -864,65 +1333,499 @@ subroutine set_clm_statevec_swc() end subroutine set_clm_statevec_swc - !> @author Yorck Ewerdwalbesloh, Anne Springer - !> @date 29.10.2025 - !> @brief set the state vector for TWS assimilation - subroutine set_clm_statevec_tws() - use clm_instMod, only : soilstate_inst, waterstate_inst - use clm_varpar , only : nlevsoi + subroutine set_clm_statevec_T(tstartcycle,mype) + use clm_instMod, only : temperature_inst + use clm_instMod, only : canopystate_inst + use clm_varpar , only : nlevgrnd + use PatchType , only : patch use ColumnType , only : col use shr_kind_mod, only: r8 => shr_kind_r8 - use PatchType, only: patch - use GridcellType, only: grc + implicit none - integer :: j,g,cc,count,c,count_c - integer :: n_c - real(r8) :: avg_sum + integer,intent(in) :: tstartcycle + integer,intent(in) :: mype - real(r8), pointer :: liq(:,:) ! liquid water (kg/m2) - real(r8), pointer :: ice(:,:) ! ice lens (kg/m2) - real(r8), pointer :: snow(:) ! snow water (mm) - real(r8), pointer :: sfc(:) ! surface water - real(r8), pointer :: can(:) ! canopy water - real(r8), pointer :: TWS(:) ! TWS + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_h2osfc(:) + real(r8), pointer :: t_soisno(:,:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) + real(r8), pointer :: tlai(:) + integer :: j,g,cc,c,p + integer :: n_c,n_p + integer :: lev + integer :: n_lev_T ! effective number of soil layers in T state vector - real(r8), pointer :: tws_state(:) - real(r8), pointer :: liq_state(:,:) - real(r8), pointer :: ice_state(:,:) - real(r8), pointer :: snow_state(:) + character (len = 34) :: fn !TSMP-PDAF: function name for swc output - cc = 0 + n_lev_T = min(nlevgrnd, clmstatevec_max_layer) - tws_state => waterstate_inst%tws_state_before - liq_state => waterstate_inst%h2osoi_liq_state_before - ice_state => waterstate_inst%h2osoi_ice_state_before - snow_state => waterstate_inst%h2osno_state_before + ! LST variables + t_grnd => temperature_inst%t_grnd_col + t_h2osfc => temperature_inst%t_h2osfc_col + t_veg => temperature_inst%t_veg_patch + t_skin => temperature_inst%t_skin_patch + t_soisno => temperature_inst%t_soisno_col + tlai => canopystate_inst%tlai_patch - select case (TWS_smoother) - case(0) ! instantaneous values - liq => waterstate_inst%h2osoi_liq_col - ice => waterstate_inst%h2osoi_ice_col - snow => waterstate_inst%h2osno_col - sfc => waterstate_inst%h2osfc_col - can => waterstate_inst%h2ocan_patch - TWS => waterstate_inst%tws_hactive - case(1) ! monthly means - liq => waterstate_inst%h2osoi_liq_col_mean - ice => waterstate_inst%h2osoi_ice_col_mean - snow => waterstate_inst%h2osno_col_mean - sfc => waterstate_inst%h2osfc_col_mean - can => waterstate_inst%h2ocan_patch_mean - TWS => waterstate_inst%tws_hactive_mean - case default +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN - error stop "Unsupported TWS_smoother" + IF(clmupdate_T/=0) THEN + ! TSMP-PDAF: Debug output of CLM t_soisno, first layer + WRITE(fn, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF - end select + END IF +#endif - select case (state_setup) - case(0) ! all compartments in state vector + !hcp LAI + if(clmupdate_T==1) then + do cc = 1, clm_varsize + ! t_grnd iterated over cols + ! t_veg iterated over patches + clm_statevec(cc) = t_grnd(state_pdaf2clm_c_p(cc)) + clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_p_p(cc+clm_varsize)) + end do + + do cc = 1, clm_paramsize + ! Works only if clm_paramsize corresponds to clm_varsize (also + ! the order) + clm_paramarr(cc) = tlai(state_pdaf2clm_p_p(cc)) + end do + endif + !end hcp LAI + + ! Skin temperature updating state vector with skin, soil and vegetation temperature. + ! Uses gridcell mean: averages over all patches (for TSKIN, TVEG) or + ! columns (for TSOIL) within each gridcell. + if(clmupdate_T==2) then + + do cc = 1, clm_varsize + + ! Get gridcell from the reference patch + g = patch%gridcell(state_pdaf2clm_p_p(cc)) + + ! --- TSKIN: average over patches in gridcell --- + clm_statevec(cc) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + if(patch%gridcell(p)==g) then + clm_statevec(cc) = clm_statevec(cc) + t_skin(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc) = clm_statevec(cc) / real(n_p, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no patches for TSKIN averaging" + error stop "Gridcell without patches in set_clm_statevec_T (TSKIN)" + end if + + ! --- TSOIL: average over columns in gridcell (per layer) --- + do lev=1,n_lev_T + clm_statevec(cc+lev*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) + t_soisno(c,lev) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " layer=", lev, " has no columns for TSOIL averaging" + error stop "Gridcell without columns in set_clm_statevec_T (TSOIL)" + end if + end do + + ! --- TVEG: average over patches in gridcell --- + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + ! Skip bare-ground/lake patches where t_veg retains spval (1e36). + if(patch%gridcell(p)==g .and. t_veg(p) < 1.0e20_r8) then + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc+(1+n_lev_T)*clm_varsize) + t_veg(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc+(1+n_lev_T)*clm_varsize) / real(n_p, r8) + else + ! No vegetated patches in this gridcell (bare ground, glacier, lake): + ! t_veg retains spval (~1e36) for all patches, so no meaningful TVEG + ! average can be formed. TSKIN is used as a fallback to keep the value + ! physically plausible within the DA step. + ! + ! Note: this slot does not affect the model state. The distribution + ! step guards against writing back to patches where t_veg >= 1e20 + ! (see update_clm_statevec_T), so the fallback value is never applied. + ! + ! The fallback does, however, participate in ensemble covariance + ! estimation during DA, introducing an artificial TSKIN-TVEG + ! correlation for these cells. A conceptually cleaner solution would + ! be to exclude non-vegetated gridcells from the TVEG part of the + ! state vector entirely (detectable at initialisation: all patches in + ! the gridcell have t_veg >= 1e20). This would require a mask array + ! and non-uniform state vector layout, breaking the current assumption + ! that every gridcell contributes the same variable block. The mapping + ! arrays state_pdaf2clm_p_p / state_clm2pdaf_p already handle + ! non-trivial patch mappings and could serve as a template for this. + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc) + end if + + end do + + ! Save prior gridcell mean state vector for computing + ! increment in updating the state vector + do cc = 1, clm_statevecsize + clm_statevec_orig(cc) = clm_statevec(cc) + end do + + endif + + ! Skin temperature updating state vector with skin, soil and vegetation + ! temperature and ground temperature. + ! Uses gridcell mean: averages over all patches (for TSKIN, TVEG) or + ! columns (for TSOIL, TGRND) within each gridcell. + if(clmupdate_T==3) then + + do cc = 1, clm_varsize + + ! Get gridcell from the reference patch + g = patch%gridcell(state_pdaf2clm_p_p(cc)) + + ! --- TSKIN: average over patches in gridcell --- + clm_statevec(cc) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + if(patch%gridcell(p)==g) then + clm_statevec(cc) = clm_statevec(cc) + t_skin(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc) = clm_statevec(cc) / real(n_p, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no patches for TSKIN averaging" + error stop "Gridcell without patches in set_clm_statevec_T (TSKIN)" + end if + + ! --- TSOIL: average over columns in gridcell (per layer) --- + do lev=1,n_lev_T + clm_statevec(cc+lev*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) + t_soisno(c,lev) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " layer=", lev, " has no columns for TSOIL averaging" + error stop "Gridcell without columns in set_clm_statevec_T (TSOIL)" + end if + end do + + ! --- TVEG: average over patches in gridcell --- + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + ! Skip bare-ground/lake patches where t_veg retains spval (1e36). + if(patch%gridcell(p)==g .and. t_veg(p) < 1.0e20_r8) then + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc+(1+n_lev_T)*clm_varsize) + t_veg(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc+(1+n_lev_T)*clm_varsize) / real(n_p, r8) + else + ! No vegetated patches: reuse TSKIN mean as safe fallback. + ! See the analogous block in clmupdate_T==2 for a detailed discussion. + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc) + end if + + ! --- TGRND: average over columns in gridcell --- + clm_statevec(cc+(2+n_lev_T)*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+(2+n_lev_T)*clm_varsize) = clm_statevec(cc+(2+n_lev_T)*clm_varsize) + t_grnd(c) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+(2+n_lev_T)*clm_varsize) = clm_statevec(cc+(2+n_lev_T)*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no columns for TGRND averaging" + error stop "Gridcell without columns in set_clm_statevec_T (TGRND)" + end if + + end do + + ! Save prior gridcell mean state vector for computing + ! increment in updating the state vector + do cc = 1, clm_statevecsize + clm_statevec_orig(cc) = clm_statevec(cc) + end do + + endif + + ! clmupdate_T==4: like ==2 but with T_H2OSFC added. + if(clmupdate_T==4) then + + do cc = 1, clm_varsize + + g = patch%gridcell(state_pdaf2clm_p_p(cc)) + + ! --- TSKIN: average over patches in gridcell --- + clm_statevec(cc) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + if(patch%gridcell(p)==g) then + clm_statevec(cc) = clm_statevec(cc) + t_skin(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc) = clm_statevec(cc) / real(n_p, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no patches for TSKIN averaging" + error stop "Gridcell without patches in set_clm_statevec_T (TSKIN)" + end if + + ! --- TSOIL: average over columns in gridcell (per layer) --- + do lev=1,n_lev_T + clm_statevec(cc+lev*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) + t_soisno(c,lev) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " layer=", lev, " has no columns for TSOIL averaging" + error stop "Gridcell without columns in set_clm_statevec_T (TSOIL)" + end if + end do + + ! --- TVEG: average over patches in gridcell --- + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + ! Skip bare-ground/lake patches where t_veg retains spval (1e36). + if(patch%gridcell(p)==g .and. t_veg(p) < 1.0e20_r8) then + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc+(1+n_lev_T)*clm_varsize) + t_veg(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc+(1+n_lev_T)*clm_varsize) / real(n_p, r8) + else + ! No vegetated patches: reuse TSKIN mean as safe fallback. + ! See the analogous block in clmupdate_T==2 for a detailed discussion. + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc) + end if + + ! --- T_H2OSFC: average over columns in gridcell --- + clm_statevec(cc+(2+n_lev_T)*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+(2+n_lev_T)*clm_varsize) = clm_statevec(cc+(2+n_lev_T)*clm_varsize) + t_h2osfc(c) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+(2+n_lev_T)*clm_varsize) = clm_statevec(cc+(2+n_lev_T)*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no columns for T_H2OSFC averaging" + error stop "Gridcell without columns in set_clm_statevec_T (T_H2OSFC)" + end if + + end do + + ! Save prior gridcell mean state vector for computing + ! increment in updating the state vector + do cc = 1, clm_statevecsize + clm_statevec_orig(cc) = clm_statevec(cc) + end do + + endif + + ! clmupdate_T==5: like ==3 but with T_H2OSFC added. + if(clmupdate_T==5) then + + do cc = 1, clm_varsize + + g = patch%gridcell(state_pdaf2clm_p_p(cc)) + + ! --- TSKIN: average over patches in gridcell --- + clm_statevec(cc) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + if(patch%gridcell(p)==g) then + clm_statevec(cc) = clm_statevec(cc) + t_skin(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc) = clm_statevec(cc) / real(n_p, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no patches for TSKIN averaging" + error stop "Gridcell without patches in set_clm_statevec_T (TSKIN)" + end if + + ! --- TSOIL: average over columns in gridcell (per layer) --- + do lev=1,n_lev_T + clm_statevec(cc+lev*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) + t_soisno(c,lev) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " layer=", lev, " has no columns for TSOIL averaging" + error stop "Gridcell without columns in set_clm_statevec_T (TSOIL)" + end if + end do + + ! --- TVEG: average over patches in gridcell --- + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + ! Skip bare-ground/lake patches where t_veg retains spval (1e36). + if(patch%gridcell(p)==g .and. t_veg(p) < 1.0e20_r8) then + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc+(1+n_lev_T)*clm_varsize) + t_veg(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc+(1+n_lev_T)*clm_varsize) / real(n_p, r8) + else + ! No vegetated patches: reuse TSKIN mean as safe fallback. + ! See the analogous block in clmupdate_T==2 for a detailed discussion. + clm_statevec(cc+(1+n_lev_T)*clm_varsize) = clm_statevec(cc) + end if + + ! --- TGRND: average over columns in gridcell --- + clm_statevec(cc+(2+n_lev_T)*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+(2+n_lev_T)*clm_varsize) = clm_statevec(cc+(2+n_lev_T)*clm_varsize) + t_grnd(c) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+(2+n_lev_T)*clm_varsize) = clm_statevec(cc+(2+n_lev_T)*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no columns for TGRND averaging" + error stop "Gridcell without columns in set_clm_statevec_T (TGRND)" + end if + + ! --- T_H2OSFC: average over columns in gridcell --- + clm_statevec(cc+(3+n_lev_T)*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+(3+n_lev_T)*clm_varsize) = clm_statevec(cc+(3+n_lev_T)*clm_varsize) + t_h2osfc(c) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+(3+n_lev_T)*clm_varsize) = clm_statevec(cc+(3+n_lev_T)*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no columns for T_H2OSFC averaging" + error stop "Gridcell without columns in set_clm_statevec_T (T_H2OSFC)" + end if + + end do + + ! Save prior gridcell mean state vector for computing + ! increment in updating the state vector + do cc = 1, clm_statevecsize + clm_statevec_orig(cc) = clm_statevec(cc) + end do + + endif + + end subroutine set_clm_statevec_T + + !> @author Yorck Ewerdwalbesloh, Anne Springer + !> @date 29.10.2025 + !> @brief set the state vector for TWS assimilation + subroutine set_clm_statevec_tws() + use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_varpar , only : nlevsoi + use ColumnType , only : col + use shr_kind_mod, only: r8 => shr_kind_r8 + use PatchType, only: patch + use GridcellType, only: grc + implicit none + + integer :: j,g,cc,count,c,count_c + integer :: n_c + real(r8) :: avg_sum + + real(r8), pointer :: liq(:,:) ! liquid water (kg/m2) + real(r8), pointer :: ice(:,:) ! ice lens (kg/m2) + real(r8), pointer :: snow(:) ! snow water (mm) + real(r8), pointer :: sfc(:) ! surface water + real(r8), pointer :: can(:) ! canopy water + real(r8), pointer :: TWS(:) ! TWS + + real(r8), pointer :: tws_state(:) + real(r8), pointer :: liq_state(:,:) + real(r8), pointer :: ice_state(:,:) + real(r8), pointer :: snow_state(:) + + cc = 0 + + tws_state => waterstate_inst%tws_state_before + liq_state => waterstate_inst%h2osoi_liq_state_before + ice_state => waterstate_inst%h2osoi_ice_state_before + snow_state => waterstate_inst%h2osno_state_before + + select case (TWS_smoother) + + case(0) ! instantaneous values + liq => waterstate_inst%h2osoi_liq_col + ice => waterstate_inst%h2osoi_ice_col + snow => waterstate_inst%h2osno_col + sfc => waterstate_inst%h2osfc_col + can => waterstate_inst%h2ocan_patch + TWS => waterstate_inst%tws_hactive + case(1) ! monthly means + liq => waterstate_inst%h2osoi_liq_col_mean + ice => waterstate_inst%h2osoi_ice_col_mean + snow => waterstate_inst%h2osno_col_mean + sfc => waterstate_inst%h2osfc_col_mean + can => waterstate_inst%h2ocan_patch_mean + TWS => waterstate_inst%tws_hactive_mean + case default + + error stop "Unsupported TWS_smoother" + + end select + + select case (state_setup) + case(0) ! all compartments in state vector cc = 1 do j = 1,nlevsoi do count = 1, num_layer(j) @@ -1217,8 +2120,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif !hcp: TG, TV - if(obs_type_update_T==1) then - error stop "Not implemented: clmupdate_T.eq.1" + if(obs_type_update_T/=0) then + call update_clm_T(tstartcycle, mype) endif ! end hcp TG, TV @@ -1312,138 +2215,781 @@ subroutine update_clm_swc(tstartcycle, mype) snow_inc(j) = waterstate_inst%h2osno_col(j) end if - end do - end do + end do + end do + + ! cc = 0 + do i=1,nlevsoi + ! CLM3.5: iterate over grid cells + ! CLM5.0: iterate over columns + ! do j=clm_begg,clm_endg + do j=clm_begc,clm_endc + + ! If snow is masked, update only, when snow depth is less than 1mm + if( (clmswc_mask_snow == 0) .or. snow_depth(j) < 0.001 ) then + ! Update only those SWCs that are not excluded by ispval + if(state_clm2pdaf_p(j,i) /= ispval) then + + if(swc(j,i)==0.0) then + swc_zero_before_update = .true. + + ! Zero-SWC leads to zero denominator in computation of + ! rliq/rice, therefore setting rliq/rice to special + ! value + rliq = spval + rice = spval + else + swc_zero_before_update = .false. + + rliq = h2osoi_liq(j,i)/(dz(j,i)*denh2o*swc(j,i)) + rice = h2osoi_ice(j,i)/(dz(j,i)*denice*swc(j,i)) + !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) + end if + + if (clmstatevec_colmean==1) then + ! If there is no significant increment, do not + ! implement any update / check. + ! + ! Note: Computing the absolute difference here, + ! because the whole state vector should be soil + ! moistures. For variables with very small values in + ! the state vector, this would have to be adapted + ! (e.g. to relative difference). + if( abs(clm_statevec(state_clm2pdaf_p(j,i)) - clm_statevec_orig(state_clm2pdaf_p(j,i))) <= 1.0e-7) then + cycle + end if + + ! Update SWC column value with the increment-factor + ! of the state vector update (state vector updates + ! are means of cols in grc) + swc_update = swc(j,i) * clm_statevec(state_clm2pdaf_p(j,i)) / clm_statevec_orig(state_clm2pdaf_p(j,i)) + else + ! Update SWC with updated state vector + swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + end if + + if(swc_update<=watmin_check) then + swc(j,i) = watmin_set + else if(swc_update>=watsat(j,i)) then + swc(j,i) = watsat(j,i) + else + swc(j,i) = swc_update + endif + + if (ieee_is_nan(swc(j,i))) then + swc(j,i) = watmin_set + print *, "WARNING: swc at j,i is nan: ", j, i + endif + + if(swc_zero_before_update) then + ! This case should not appear for hydrologically + ! active columns/layers, where always: swc > watmin + ! + ! If you want to make sure that no zero SWCs appear in + ! the code, comment out the error stop + +#ifdef PDAF_DEBUG + ! error stop "ERROR: Update of zero-swc" + print *, "WARNING: Update of zero-swc" + print *, "WARNING: Any new H2O added to h2osoi_liq(j,i) with j,i = ", j, i +#endif + h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o + h2osoi_ice(j,i) = 0.0 + else + ! update liquid water content + h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o*rliq + ! update ice content + h2osoi_ice(j,i) = swc(j,i) * dz(j,i)*denice*rice + end if + + end if + end if + ! cc = cc + 1 + end do + end do + +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn3, action="write") + WRITE (71,"(es22.15)") h2osoi_liq(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn4, action="write") + WRITE (71,"(es22.15)") h2osoi_ice(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") swc(:,:) + CLOSE(71) + + END IF +#endif + + do i = 1,nlevsoi + do j=clm_begc,clm_endc + + liq_inc(j,i) = h2osoi_liq(j,i)-liq_inc(j,i) + ice_inc(j,i) = h2osoi_ice(j,i)-ice_inc(j,i) + + if (i==1) then + snow_inc(j) = waterstate_inst%h2osno_col(j)-snow_inc(j) + end if + + end do + end do + + end subroutine update_clm_swc + + subroutine update_clm_T(tstartcycle,mype) + + use PatchType , only : patch + use clm_varpar , only : nlevgrnd + use clm_instMod, only : temperature_inst + use clm_instMod, only : waterstate_inst + use IEEE_ARITHMETIC, only: ieee_is_nan + use shr_const_mod, only: SHR_CONST_TKFRZ + + implicit none + + integer,intent(in) :: tstartcycle + integer,intent(in) :: mype + + integer :: i + integer :: j + integer :: c + integer :: p + integer :: cc + integer :: lev + integer :: n_lev_T ! effective number of soil layers in T state vector + + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_h2osfc(:) + real(r8), pointer :: t_soisno(:,:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) + + real(r8), pointer :: snow_depth(:) + + real(r8) :: increment_factor + real(r8) :: t_update + + character (len = 31) :: fn !TSMP-PDAF: function name for swc output + + integer :: incr_warn_count_skin, incr_warn_count_soisno, & + incr_warn_count_veg, incr_warn_count_grnd, incr_warn_count_h2osfc + + ! Guard against applying t_soisno increment multiple times to the same + ! column when several patches share a column (T==2 loop is over patches). + logical, allocatable :: col_updated(:) + + ! LST + t_grnd => temperature_inst%t_grnd_col + t_h2osfc => temperature_inst%t_h2osfc_col + t_soisno => temperature_inst%t_soisno_col + t_veg => temperature_inst%t_veg_patch + t_skin => temperature_inst%t_skin_patch + ! tlai => canopystate_inst%tlai_patch + + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + + incr_warn_count_skin = 0 + incr_warn_count_soisno = 0 + incr_warn_count_veg = 0 + incr_warn_count_grnd = 0 + incr_warn_count_h2osfc = 0 + + allocate(col_updated(clm_begc:clm_endc)) + col_updated = .false. + + n_lev_T = min(nlevgrnd, clmstatevec_max_layer) + + !hcp: TG, TV + if(clmupdate_T==1) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1)) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + end do + endif + ! end hcp TG, TV + + ! Skin temperature updating skin, soil and vegetation temperature. + ! Uses gridcell mean increment factor: applies the ratio of + ! (new gridcell mean / old gridcell mean) to each patch/column value. + if(clmupdate_T==2) then + + do p = clm_begp, clm_endp + c = patch%column(p) + + ! If snow is masked, update only, when snow depth is less than 1mm + mask_snow_1: if( (clmT_mask_snow == 0) .or. snow_depth(c) < clmT_mask_snow_depth ) then + ! No update for (near-to) freezing soil temperatures + mask_freeze_1: if( t_soisno(c,1) > SHR_CONST_TKFRZ + clmT_mask_T ) then + + ! --- TSKIN: update with increment factor --- + cc = state_clm2pdaf_p(p,1) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_skin(p) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_skin increment_factor is NaN at p=", p, " - leaving t_skin unchanged" + t_update = t_skin(p) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_skin(p) + increment_factor + else + t_update = t_skin(p) + sign(clmT_max_increment, increment_factor) + incr_warn_count_skin = incr_warn_count_skin + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_skin update is NaN at p=", p + else + t_skin(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if + + ! Guard: only update each column once even if multiple patches share it. + if (.not. col_updated(c)) then + ! --- TSOIL: update with increment factor for each layer --- + do lev=1,n_lev_T + cc = state_clm2pdaf_p(p, 1+lev) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_soisno(c,lev) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_soisno increment_factor is NaN at c=", c, " lev=", lev, " - leaving t_soisno unchanged" + t_update = t_soisno(c,lev) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_soisno(c,lev) + increment_factor + else + t_update = t_soisno(c,lev) + sign(clmT_max_increment, increment_factor) + incr_warn_count_soisno = incr_warn_count_soisno + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_soisno update is NaN at c=", c, " lev=", lev + else + t_soisno(c,lev) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if + end do + + col_updated(c) = .true. + end if ! col_updated + + ! --- TVEG: update with increment factor --- + cc = state_clm2pdaf_p(p, 2+n_lev_T) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7 .and. t_veg(p) < 1.0e20_r8) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_veg(p) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_veg increment_factor is NaN at p=", p, " - leaving t_veg unchanged" + t_update = t_veg(p) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_veg(p) + increment_factor + else + t_update = t_veg(p) + sign(clmT_max_increment, increment_factor) + incr_warn_count_veg = incr_warn_count_veg + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_veg update is NaN at p=", p + else + t_veg(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if + + end if mask_freeze_1 + end if mask_snow_1 + + end do + if (incr_warn_count_skin > 0) print *, "WARNING: t_skin total increments exceeding T_max_increment:", & + incr_warn_count_skin + if (incr_warn_count_soisno > 0) print *, "WARNING: t_soisno total increments exceeding T_max_increment:", & + incr_warn_count_soisno + if (incr_warn_count_veg > 0) print *, "WARNING: t_veg total increments exceeding T_max_increment:", & + incr_warn_count_veg + deallocate(col_updated) + endif + + ! Skin temperature updating skin, soil, vegetation and ground temperature. + ! Uses gridcell mean increment factor: applies the ratio of + ! (new gridcell mean / old gridcell mean) to each patch/column value. + if(clmupdate_T==3) then + + do p = clm_begp, clm_endp + c = patch%column(p) + + ! If snow is masked, update only, when snow depth is less than 1mm + mask_snow_2: if( (clmT_mask_snow == 0) .or. snow_depth(c) < clmT_mask_snow_depth ) then + ! No update for (near-to) freezing soil temperatures + mask_freeze_2: if( t_soisno(c,1) > SHR_CONST_TKFRZ + clmT_mask_T ) then + + ! --- TSKIN: update with increment factor --- + cc = state_clm2pdaf_p(p,1) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_skin(p) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_skin increment_factor is NaN at p=", p, " - leaving t_skin unchanged" + t_update = t_skin(p) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_skin(p) + increment_factor + else + t_update = t_skin(p) + sign(clmT_max_increment, increment_factor) + incr_warn_count_skin = incr_warn_count_skin + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_skin update is NaN at p=", p + else + t_skin(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if + + ! Guard: only update each column once even if multiple patches share it. + if (.not. col_updated(c)) then + ! --- TSOIL: update with increment factor for each layer --- + do lev=1,n_lev_T + cc = state_clm2pdaf_p(p, 1+lev) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_soisno(c,lev) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_soisno increment_factor is NaN at c=", c, " lev=", lev, " - leaving t_soisno unchanged" + t_update = t_soisno(c,lev) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_soisno(c,lev) + increment_factor + else + t_update = t_soisno(c,lev) + sign(clmT_max_increment, increment_factor) + incr_warn_count_soisno = incr_warn_count_soisno + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_soisno update is NaN at c=", c, " lev=", lev + else + t_soisno(c,lev) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if + end do + end if ! col_updated + + ! --- TVEG: update with increment factor --- + cc = state_clm2pdaf_p(p, 2+n_lev_T) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7 .and. t_veg(p) < 1.0e20_r8) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_veg(p) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_veg increment_factor is NaN at p=", p, " - leaving t_veg unchanged" + t_update = t_veg(p) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_veg(p) + increment_factor + else + t_update = t_veg(p) + sign(clmT_max_increment, increment_factor) + incr_warn_count_veg = incr_warn_count_veg + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_veg update is NaN at p=", p + else + t_veg(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if + + ! Guard: only update each column once even if multiple patches share it. + if (.not. col_updated(c)) then + ! --- TGRND: update with increment factor --- + cc = state_clm2pdaf_p(p, 3+n_lev_T) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_grnd(c) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_grnd increment_factor is NaN at c=", c, " - leaving t_grnd unchanged" + t_update = t_grnd(c) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_grnd(c) + increment_factor + else + t_update = t_grnd(c) + sign(clmT_max_increment, increment_factor) + incr_warn_count_grnd = incr_warn_count_grnd + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_grnd update is NaN at c=", c + else + t_grnd(c) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if + + col_updated(c) = .true. + end if ! col_updated - ! cc = 0 - do i=1,nlevsoi - ! CLM3.5: iterate over grid cells - ! CLM5.0: iterate over columns - ! do j=clm_begg,clm_endg - do j=clm_begc,clm_endc + end if mask_freeze_2 + end if mask_snow_2 - ! If snow is masked, update only, when snow depth is less than 1mm - if( (clmswc_mask_snow == 0) .or. snow_depth(j) < 0.001 ) then - ! Update only those SWCs that are not excluded by ispval - if(state_clm2pdaf_p(j,i) /= ispval) then + end do + if (incr_warn_count_skin > 0) print *, "WARNING: t_skin total increments exceeding T_max_increment:", & + incr_warn_count_skin + if (incr_warn_count_soisno > 0) print *, "WARNING: t_soisno total increments exceeding T_max_increment:", & + incr_warn_count_soisno + if (incr_warn_count_veg > 0) print *, "WARNING: t_veg total increments exceeding T_max_increment:", & + incr_warn_count_veg + if (incr_warn_count_grnd > 0) print *, "WARNING: t_grnd total increments exceeding T_max_increment:", & + incr_warn_count_grnd + deallocate(col_updated) + endif - if(swc(j,i)==0.0) then - swc_zero_before_update = .true. + ! clmupdate_T==4: like ==2 but also updating T_H2OSFC. + if(clmupdate_T==4) then - ! Zero-SWC leads to zero denominator in computation of - ! rliq/rice, therefore setting rliq/rice to special - ! value - rliq = spval - rice = spval - else - swc_zero_before_update = .false. + do p = clm_begp, clm_endp + c = patch%column(p) - rliq = h2osoi_liq(j,i)/(dz(j,i)*denh2o*swc(j,i)) - rice = h2osoi_ice(j,i)/(dz(j,i)*denice*swc(j,i)) - !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) - end if + ! If snow is masked, update only, when snow depth is less than 1mm + mask_snow_3: if( (clmT_mask_snow == 0) .or. snow_depth(c) < clmT_mask_snow_depth ) then + ! No update for (near-to) freezing soil temperatures + mask_freeze_3: if( t_soisno(c,1) > SHR_CONST_TKFRZ + clmT_mask_T ) then - if (clmstatevec_colmean==1) then - ! If there is no significant increment, do not - ! implement any update / check. - ! - ! Note: Computing the absolute difference here, - ! because the whole state vector should be soil - ! moistures. For variables with very small values in - ! the state vector, this would have to be adapted - ! (e.g. to relative difference). - if( abs(clm_statevec(state_clm2pdaf_p(j,i)) - clm_statevec_orig(state_clm2pdaf_p(j,i))) <= 1.0e-7) then - cycle - end if + ! --- TSKIN: update with increment factor --- + cc = state_clm2pdaf_p(p,1) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_skin(p) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_skin increment_factor is NaN at p=", p, " - leaving t_skin unchanged" + t_update = t_skin(p) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_skin(p) + increment_factor + else + t_update = t_skin(p) + sign(clmT_max_increment, increment_factor) + incr_warn_count_skin = incr_warn_count_skin + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_skin update is NaN at p=", p + else + t_skin(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if - ! Update SWC column value with the increment-factor - ! of the state vector update (state vector updates - ! are means of cols in grc) - swc_update = swc(j,i) * clm_statevec(state_clm2pdaf_p(j,i)) / clm_statevec_orig(state_clm2pdaf_p(j,i)) - else - ! Update SWC with updated state vector - swc_update = clm_statevec(state_clm2pdaf_p(j,i)) - end if + ! Guard: only update each column once even if multiple patches share it. + if (.not. col_updated(c)) then + ! --- TSOIL: update with increment factor for each layer --- + do lev=1,n_lev_T + cc = state_clm2pdaf_p(p, 1+lev) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_soisno(c,lev) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_soisno increment_factor is NaN at c=", c, " lev=", lev, " - leaving t_soisno unchanged" + t_update = t_soisno(c,lev) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_soisno(c,lev) + increment_factor + else + t_update = t_soisno(c,lev) + sign(clmT_max_increment, increment_factor) + incr_warn_count_soisno = incr_warn_count_soisno + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_soisno update is NaN at c=", c, " lev=", lev + else + t_soisno(c,lev) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if + end do + end if ! col_updated + + ! --- TVEG: update with increment factor --- + cc = state_clm2pdaf_p(p, 2+n_lev_T) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7 .and. t_veg(p) < 1.0e20_r8) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_veg(p) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_veg increment_factor is NaN at p=", p, " - leaving t_veg unchanged" + t_update = t_veg(p) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_veg(p) + increment_factor + else + t_update = t_veg(p) + sign(clmT_max_increment, increment_factor) + incr_warn_count_veg = incr_warn_count_veg + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_veg update is NaN at p=", p + else + t_veg(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if - if(swc_update<=watmin_check) then - swc(j,i) = watmin_set - else if(swc_update>=watsat(j,i)) then - swc(j,i) = watsat(j,i) - else - swc(j,i) = swc_update - endif + ! Guard: only update each column once even if multiple patches share it. + if (.not. col_updated(c)) then + ! --- T_H2OSFC: update with increment factor --- + cc = state_clm2pdaf_p(p, 3+n_lev_T) + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_h2osfc(c) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_h2osfc increment_factor is NaN at c=", c, " - leaving t_h2osfc unchanged" + t_update = t_h2osfc(c) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_h2osfc(c) + increment_factor + else + t_update = t_h2osfc(c) + sign(clmT_max_increment, increment_factor) + incr_warn_count_h2osfc = incr_warn_count_h2osfc + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_h2osfc update is NaN at c=", c + else + t_h2osfc(c) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if - if (ieee_is_nan(swc(j,i))) then - swc(j,i) = watmin_set - print *, "WARNING: swc at j,i is nan: ", j, i - endif + col_updated(c) = .true. + end if ! col_updated - if(swc_zero_before_update) then - ! This case should not appear for hydrologically - ! active columns/layers, where always: swc > watmin - ! - ! If you want to make sure that no zero SWCs appear in - ! the code, comment out the error stop + end if mask_freeze_3 + end if mask_snow_3 -#ifdef PDAF_DEBUG - ! error stop "ERROR: Update of zero-swc" - print *, "WARNING: Update of zero-swc" - print *, "WARNING: Any new H2O added to h2osoi_liq(j,i) with j,i = ", j, i -#endif - h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o - h2osoi_ice(j,i) = 0.0 - else - ! update liquid water content - h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o*rliq - ! update ice content - h2osoi_ice(j,i) = swc(j,i) * dz(j,i)*denice*rice - end if + end do + if (incr_warn_count_skin > 0) print *, "WARNING: t_skin total increments exceeding T_max_increment:", & + incr_warn_count_skin + if (incr_warn_count_soisno > 0) print *, "WARNING: t_soisno total increments exceeding T_max_increment:", & + incr_warn_count_soisno + if (incr_warn_count_veg > 0) print *, "WARNING: t_veg total increments exceeding T_max_increment:", & + incr_warn_count_veg + if (incr_warn_count_h2osfc > 0) print *, "WARNING: t_h2osfc total increments exceeding T_max_increment:", & + incr_warn_count_h2osfc + deallocate(col_updated) + endif - end if - end if - ! cc = cc + 1 - end do - end do + ! clmupdate_T==5: like ==3 but also updating T_H2OSFC. + if(clmupdate_T==5) then -#ifdef PDAF_DEBUG - IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN + do p = clm_begp, clm_endp + c = patch%column(p) - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn3, action="write") - WRITE (71,"(es22.15)") h2osoi_liq(:,:) - CLOSE(71) + ! If snow is masked, update only, when snow depth is less than 1mm + mask_snow_4: if( (clmT_mask_snow == 0) .or. snow_depth(c) < clmT_mask_snow_depth ) then + ! No update for (near-to) freezing soil temperatures + mask_freeze_4: if( t_soisno(c,1) > SHR_CONST_TKFRZ + clmT_mask_T ) then - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn4, action="write") - WRITE (71,"(es22.15)") h2osoi_ice(:,:) - CLOSE(71) + ! --- TSKIN: update with increment factor --- + cc = state_clm2pdaf_p(p,1) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_skin(p) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_skin increment_factor is NaN at p=", p, " - leaving t_skin unchanged" + t_update = t_skin(p) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_skin(p) + increment_factor + else + t_update = t_skin(p) + sign(clmT_max_increment, increment_factor) + incr_warn_count_skin = incr_warn_count_skin + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_skin update is NaN at p=", p + else + t_skin(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn2, action="write") - WRITE (71,"(es22.15)") swc(:,:) - CLOSE(71) + ! Guard: only update each column once even if multiple patches share it. + if (.not. col_updated(c)) then + ! --- TSOIL: update with increment factor for each layer --- + do lev=1,n_lev_T + cc = state_clm2pdaf_p(p, 1+lev) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_soisno(c,lev) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_soisno increment_factor is NaN at c=", c, " lev=", lev, " - leaving t_soisno unchanged" + t_update = t_soisno(c,lev) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_soisno(c,lev) + increment_factor + else + t_update = t_soisno(c,lev) + sign(clmT_max_increment, increment_factor) + incr_warn_count_soisno = incr_warn_count_soisno + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_soisno update is NaN at c=", c, " lev=", lev + else + t_soisno(c,lev) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if + end do + end if ! col_updated + + ! --- TVEG: update with increment factor --- + cc = state_clm2pdaf_p(p, 2+n_lev_T) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7 .and. t_veg(p) < 1.0e20_r8) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_veg(p) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_veg increment_factor is NaN at p=", p, " - leaving t_veg unchanged" + t_update = t_veg(p) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_veg(p) + increment_factor + else + t_update = t_veg(p) + sign(clmT_max_increment, increment_factor) + incr_warn_count_veg = incr_warn_count_veg + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_veg update is NaN at p=", p + else + t_veg(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if - END IF -#endif + ! Guard: only update each column once even if multiple patches share it. + if (.not. col_updated(c)) then + ! --- TGRND: update with increment factor --- + cc = state_clm2pdaf_p(p, 3+n_lev_T) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_grnd(c) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_grnd increment_factor is NaN at c=", c, " - leaving t_grnd unchanged" + t_update = t_grnd(c) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_grnd(c) + increment_factor + else + t_update = t_grnd(c) + sign(clmT_max_increment, increment_factor) + incr_warn_count_grnd = incr_warn_count_grnd + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_grnd update is NaN at c=", c + else + t_grnd(c) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if - do i = 1,nlevsoi - do j=clm_begc,clm_endc + ! --- T_H2OSFC: update with increment factor --- + cc = state_clm2pdaf_p(p, 4+n_lev_T) + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + if( (clmincrement_type == 0)) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_h2osfc(c) * increment_factor + else + increment_factor = clm_statevec(cc) - clm_statevec_orig(cc) + if (ieee_is_nan(increment_factor)) then + print *, "WARNING: t_h2osfc increment_factor is NaN at c=", c, " - leaving t_h2osfc unchanged" + t_update = t_h2osfc(c) + else if (abs(increment_factor) < clmT_max_increment) then + t_update = t_h2osfc(c) + increment_factor + else + t_update = t_h2osfc(c) + sign(clmT_max_increment, increment_factor) + incr_warn_count_h2osfc = incr_warn_count_h2osfc + 1 + end if + end if + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_h2osfc update is NaN at c=", c + else + t_h2osfc(c) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update)) + end if + end if - liq_inc(j,i) = h2osoi_liq(j,i)-liq_inc(j,i) - ice_inc(j,i) = h2osoi_ice(j,i)-ice_inc(j,i) + col_updated(c) = .true. + end if ! col_updated - if (i==1) then - snow_inc(j) = waterstate_inst%h2osno_col(j)-snow_inc(j) - end if + end if mask_freeze_4 + end if mask_snow_4 end do - end do + if (incr_warn_count_skin > 0) print *, "WARNING: t_skin total increments exceeding T_max_increment:", & + incr_warn_count_skin + if (incr_warn_count_soisno > 0) print *, "WARNING: t_soisno total increments exceeding T_max_increment:", & + incr_warn_count_soisno + if (incr_warn_count_veg > 0) print *, "WARNING: t_veg total increments exceeding T_max_increment:", & + incr_warn_count_veg + if (incr_warn_count_grnd > 0) print *, "WARNING: t_grnd total increments exceeding T_max_increment:", & + incr_warn_count_grnd + if (incr_warn_count_h2osfc > 0) print *, "WARNING: t_h2osfc total increments exceeding T_max_increment:", & + incr_warn_count_h2osfc + deallocate(col_updated) + endif - end subroutine update_clm_swc +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF +#endif + end subroutine update_clm_T subroutine update_clm_texture(tstartcycle, mype) use clm_varpar , only : nlevsoi @@ -2285,6 +3831,7 @@ subroutine init_n_domains_clm(n_domains_p) use decompMod, only : get_proc_bounds use clm_varcon , only : ispval use ColumnType , only : col + use PatchType , only : patch implicit none @@ -2292,6 +3839,7 @@ subroutine init_n_domains_clm(n_domains_p) integer :: domain_p integer :: begg, endg ! per-proc gridcell ending gridcell indices integer :: begc, endc ! per-proc beginning and ending column indices + integer :: begp, endp ! per-proc beginning and ending pft indices integer :: g integer :: c @@ -2299,7 +3847,7 @@ subroutine init_n_domains_clm(n_domains_p) ! TODO: remove unnecessary calls of get_proc_bounds (use clm_begg, ! clm_endg, etc) - call get_proc_bounds(begg=begg, endg=endg, begc=begc, endc=endc) + call get_proc_bounds(begg=begg, endg=endg, begc=begc, endc=endc, begp=begp, endp=endp) if(clmupdate_swc==1) then if(clmstatevec_allcol==1) then @@ -2311,6 +3859,12 @@ subroutine init_n_domains_clm(n_domains_p) ! -> DIM_L: number of layers in gridcell n_domains_p = endg - begg + 1 end if + else if (clmupdate_T/=0) then + ! For LSTDA: gridcells are domains + ! Each gridcell is a local domain + ! -> DIM_L: number of temperature variables (each with gridcell + ! -> averages) + n_domains_p = endg - begg + 1 else ! Process-local number of gridcells Default, possibly not tested ! for other updates except SWC @@ -2327,86 +3881,125 @@ subroutine init_n_domains_clm(n_domains_p) ! local domain domain_p (from the column, the gridcell can be ! derived) - ! Allocate state_loc2clm_c_p with preliminary n_domains_p - IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) - allocate(state_loc2clm_c_p(n_domains_p)) - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = ispval - end do - - if(clmstatevec_only_active == 1) then - - ! Reset n_domains_p - n_domains_p = 0 - domain_p = 0 - - if(clmstatevec_allcol == 1) then - ! COLUMNS - - ! Each hydrologically active layer is a local domain - ! -> DIM_L: number of layers in hydrologically active column - do c=clm_begc,clm_endc - ! Skip state vector loop and directly check if column is - ! hydrologically active - if(col%hydrologically_active(c)) then - domain_p = domain_p + 1 - n_domains_p = n_domains_p + 1 - state_loc2clm_c_p(domain_p) = c - end if - end do + if(clmupdate_swc==1) then + ! Allocate state_loc2clm_c_p with preliminary n_domains_p + IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) + allocate(state_loc2clm_c_p(n_domains_p)) + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = ispval + end do - else - ! GRIDCELLS + if(clmstatevec_only_active == 1) then - ! For gridcells - do g = clm_begg,clm_endg + ! Reset n_domains_p + n_domains_p = 0 + domain_p = 0 - ! Search the state vector for col in grc - do cc = 1,clm_statevecsize + if(clmstatevec_allcol == 1) then + ! COLUMNS - if (col%gridcell(state_pdaf2clm_c_p(cc)) == g) then - ! Set local domain index + ! Each hydrologically active column is a local domain + ! -> DIM_L: number of layers in hydrologically active column + do c=clm_begc,clm_endc + ! Skip state vector loop and directly check if column is + ! hydrologically active + if(col%hydrologically_active(c)) then domain_p = domain_p + 1 - ! Set new number of local domains n_domains_p = n_domains_p + 1 - ! Set CLM-column-index corresponding to local domain - state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(cc) - ! Exit state vector loop, when fitting column is found - exit + state_loc2clm_c_p(domain_p) = c end if end do - end do + else + ! GRIDCELLS + + ! For gridcells + do g = clm_begg,clm_endg + + ! Search the state vector for col in grc + do cc = 1,clm_statevecsize + + if (col%gridcell(state_pdaf2clm_c_p(cc)) == g) then + ! Set local domain index + domain_p = domain_p + 1 + ! Set new number of local domains + n_domains_p = n_domains_p + 1 + ! Set CLM-column-index corresponding to local domain + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(cc) + ! Exit state vector loop, when fitting column is found + exit + end if + end do - end if + end do - else + end if - ! Set state_loc2clm_c_p for non-excluding hydrologically - ! inactive cols/grcs - if(clmstatevec_allcol == 1) then - ! COLUMNS - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = clm_begc + domain_p - 1 - end do else - ! GRIDCELLS - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = clm_begg + domain_p - 1 - end do - end if + ! Set state_loc2clm_c_p for non-excluding hydrologically + ! inactive cols/grcs + if(clmstatevec_allcol == 1) then + ! COLUMNS + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = clm_begc + domain_p - 1 + end do + else + ! GRIDCELLS + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(domain_p) + end do + end if + + end if end if ! Possibly: Warning when final n_domains_p actually excludes ! hydrologically inactive gridcells + if(clmupdate_T/=0) then + ! For LSTDA: Gridcells are domains + + ! Allocate state_loc2clm_c_p with preliminary n_domains_p + IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) + allocate(state_loc2clm_c_p(n_domains_p)) + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = ispval + end do + IF (allocated(state_loc2clm_p_p)) deallocate(state_loc2clm_p_p) + allocate(state_loc2clm_p_p(n_domains_p)) + do domain_p=1,n_domains_p + state_loc2clm_p_p(domain_p) = ispval + end do + + ! Columns from gridcells + ! Patches from gridcells + ! + ! domain_p is equal to cc for the first variable TSKIN + ! Therefore, state_pdaf2clm_c_p / state_pdaf2clm_p_p can be used + ! to derive the corresponding column / patch + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(domain_p) + state_loc2clm_p_p(domain_p) = state_pdaf2clm_p_p(domain_p) + end do + + end if + else NOGRACE n_domains_p = num_hactiveg end if NOGRACE +#ifdef PDAF_DEBUG + WRITE(*, '(a,x,a,i10,x,a,i10)') "TSMP-PDAF-debug", "begc=", begc, & + "init_n_domains_clm: n_domains_p=", n_domains_p + if (allocated(state_loc2clm_c_p)) then + WRITE(*, '(a,x,a,i10,x,a,*(i10))') "TSMP-PDAF-debug", "begc=", begc, & + "init_n_domains_clm: state_loc2clm_c_p=", state_loc2clm_c_p + end if +#endif + end subroutine init_n_domains_clm @@ -2417,6 +4010,7 @@ end subroutine init_n_domains_clm !> This routine sets DIM_L, the local state vector dimension. subroutine init_dim_l_clm(domain_p, dim_l) use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd use ColumnType , only : col implicit none @@ -2453,6 +4047,31 @@ subroutine init_dim_l_clm(domain_p, dim_l) dim_l = 3*nlevsoi + nshift endif + if(clmupdate_T==1) then + ! TG + TV: 2 temperatures per patch + dim_l = 2 + endif + + if(clmupdate_T==2) then + ! TSKIN + TSOIL(n_lev_T layers) + TV + dim_l = 2 + min(nlevgrnd, clmstatevec_max_layer) + endif + + if(clmupdate_T==3) then + ! TSKIN + TSOIL(n_lev_T layers) + TV + TGRND + dim_l = 3 + min(nlevgrnd, clmstatevec_max_layer) + endif + + if(clmupdate_T==4) then + ! TSKIN + TSOIL(n_lev_T layers) + TV + T_H2OSFC + dim_l = 3 + min(nlevgrnd, clmstatevec_max_layer) + endif + + if(clmupdate_T==5) then + ! TSKIN + TSOIL(n_lev_T layers) + TV + TGRND + T_H2OSFC + dim_l = 4 + min(nlevgrnd, clmstatevec_max_layer) + endif + if (clmupdate_tws==1) then dim_l = 0 g = hactiveg_levels(domain_p,1) @@ -2536,11 +4155,22 @@ subroutine g2l_state_clm(domain_p, dim_p, state_p, dim_l, state_l) ! ENDDO NOGRACE: if (clmupdate_tws/=1) then ! Column index inside gridcell index domain_p + if(clmupdate_swc==1) then DO i = 1, dim_l ! Column index from DOMAIN_P via STATE_LOC2CLM_C_P ! Layer index: i state_l(i) = state_p(state_clm2pdaf_p(state_loc2clm_c_p(domain_p),i)) END DO + end if + + if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then + DO i = 1, dim_l + ! Patch index from DOMAIN_P via STATE_LOC2CLM_P_P + ! Variable index: i + state_l(i) = state_p(state_clm2pdaf_p(state_loc2clm_p_p(domain_p),i)) + END DO + end if + else NOGRACE if (clm_varsize_tws(5)/=0) then @@ -2658,11 +4288,22 @@ subroutine l2g_state_clm(domain_p, dim_l, state_l, dim_p, state_p) ! ENDDO NOGRACE: if (clmupdate_tws==0) then ! Column index inside gridcell index domain_p + if(clmupdate_swc==1) then DO i = 1, dim_l ! Column index from DOMAIN_P via STATE_LOC2CLM_C_P ! Layer index i state_p(state_clm2pdaf_p(state_loc2clm_c_p(domain_p),i)) = state_l(i) END DO + end if + + if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then + DO i = 1, dim_l + ! Patch index from DOMAIN_P via STATE_LOC2CLM_P_P + ! Variable index: i + state_p(state_clm2pdaf_p(state_loc2clm_p_p(domain_p),i)) = state_l(i) + END DO + end if + else NOGRACE if (clm_varsize_tws(5)/=0) then