diff --git a/docs/_toc.yml b/docs/_toc.yml index 19e34413..cfa774df 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -28,6 +28,7 @@ parts: - file: users_guide/running_tsmp_pdaf/input_cmd - file: users_guide/running_tsmp_pdaf/input_obs - file: users_guide/running_tsmp_pdaf/input_enkfpf + - file: users_guide/running_tsmp_pdaf/snow_da - file: users_guide/debugging_tsmp_pdaf/README title: Debugging TSMP-PDAF diff --git a/docs/users_guide/running_tsmp_pdaf/README.md b/docs/users_guide/running_tsmp_pdaf/README.md index 7d1809ab..69648a12 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 snow observations into CLM, see [Snow Data +Assimilation](snowda). + See the Virtual Machine download on webpage . diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md index 790b4aa0..8239fcde 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md +++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md @@ -54,6 +54,8 @@ nprocs = update_swc = update_texture = update_T = +update_snow = +update_snow_repartitioning = print_swc = print_et = statevec_allcol = @@ -601,6 +603,92 @@ Only takes effect if `CLM:update_swc``is switched on. Default setting is `0`: No masking of columns with snow cover. +(enkfpf:clm:update_snow)= +### CLM:update_snow ### + +`CLM:update_snow`: (integer) Flag for snow data assimilation (Snow-DA) +in CLM. Only implemented for eCLM; not available for CLM3.5. + +For a detailed description of the Snow-DA workflow see [Snow Data +Assimilation](snowda). + +Default: `0` + +- `0`: No snow update. + +- `1`: Assimilation of snow depth (SD). The state vector contains one + SD value per grid cell (`snow_depth_col`). After the PDAF analysis + step, the updated SD is written back to eCLM and the snow layer + variables are redistributed according to + `CLM:update_snow_repartitioning`. + +- `2`: Assimilation of snow water equivalent (SWE). The state vector + contains one SWE value per grid cell (`h2osno_col`). After the PDAF + analysis step, the updated SWE is written back to CLM and the snow + layer variables are redistributed. + +- `3`: Both SD and SWE are placed in the state vector (2 values per + grid cell). The ice content of each snow layer (`h2osoi_ice`) is + scaled by the SWE increment. Only compatible with + `CLM:update_snow_repartitioning=3`. + +- `4`: As `3`, and additionally scales liquid water content + (`h2osoi_liq`) by the SWE increment and layer thickness (`dz`) by + the SD increment. Only compatible with + `CLM:update_snow_repartitioning=3`. + +- `5`: As `3`, and additionally scales layer thickness (`dz`) by the + SD increment (liquid water content not updated). Only compatible + with `CLM:update_snow_repartitioning=3`. + +- `6`: As `3`, and additionally scales liquid water content + (`h2osoi_liq`) by the SWE increment (layer thickness not + updated). Only compatible with `CLM:update_snow_repartitioning=3`. + +- `7`: As `3`, with only `h2osoi_ice` scaled by the SWE increment + (no `dz`, no `h2osoi_liq` update). Intended to be functionally + equivalent to case `3` after refactoring. Only compatible with + `CLM:update_snow_repartitioning=3`. + + +(enkfpf:clm:update_snow_repartitioning)= +### CLM:update_snow_repartitioning ### + +`CLM:update_snow_repartitioning`: (integer) Method for distributing +the DA-updated bulk snow quantity across the individual CLM snow +layers. + +After the PDAF analysis step, the bulk snow variables (SD or SWE) in +the state vector are updated, but the individual snow layer variables +(`h2osoi_ice`, `h2osoi_liq`, `dz`) must be adjusted +consistently. This parameter controls the redistribution method. + +Default: `3` + +- `1`: Adjust bottom snow layer only (based on the DART + repartitioning approach). The full SWE change is applied to the + bottom layer; upper layers remain unchanged. Only available for + `CLM:update_snow=1` or `2`. + +- `2`: Adjust all active snow layers proportionally to their current + water content. Each layer receives a fraction of the SWE change + proportional to `(h2osoi_liq + h2osoi_ice) / SWE_prior`. Only + available for `CLM:update_snow=1` or `2`. + +- `3` (default): Increment-based update. Layer variables are scaled + by the ratio of the updated (posterior) to the prior bulk snow + value: + + - For `CLM:update_snow=1`: `h2osoi_ice` is scaled by `SD_out / + SD_in`. + - For `CLM:update_snow=2,3`: `h2osoi_ice` is scaled by `SWE_out / + SWE_in`. + - For `CLM:update_snow=4..7`: see [Snow Data Assimilation](snow_da) + for the full variable update table. + + Columns where either the prior or the posterior bulk snow value is + below `1e-6` are excluded from the increment update. + (enkfpf:cosmo)= ## [COSMO] ## @@ -861,61 +949,63 @@ Default: 0, output turned off. (enkfpf:summary)= ## Parameter Summary ## - | section | parameter | default value | - |:---------:|:-----------------------:|:-------------:| - | `[PF]` | | | - | | `problemname` | \- | - | | `nprocs` | 0 | - | | `starttime` | 0.0 | - | | `endtime` | 0 | - | | `simtime` | 0 | - | | `dt` | 0.0 | - | | `updateflag` | 1 | - | | `paramupdate` | 0 | - | | `paramupdate_frequency` | 1 | - | | `dampingfactor_param` | 1.0 | - | | `dampingfactor_state` | 1.0 | - | | `damping_switch_sm` | 1 | - | | `aniso_perm_y` | 1.0 | - | | `aniso_perm_z` | 1.0 | - | | `aniso_use_parflow` | 0 | - | | | | - | | `printensemble` | 1 | - | | `t_printensemble` | -2 | - | | `printstat` | 1 | - | | `paramprintensemble` | 1 | - | | `paramprintstat` | 1 | - | | `olfmasking` | 0 | - | | `olfmasking_param` | 0 | - | | `olfmasking_depth` | 0 | - | `[CLM]` | | | - | | `problemname` | \- | - | | `nprocs` | 0 | - | | `update_swc` | 1 | - | | `print_swc` | 0 | - | | `print_et` | 0 | - | | `statevec_allcol` | 0 | - | | `statevec_colmean` | 0 | - | | `statevec_only_active` | 0 | - | | `statevec_max_layer` | 25 | - | | `t_printensemble` | -2 | - | | `watmin_switch` | 0 | - | `[COSMO]` | | | - | | `nprocs` | 0 | - | | `dtmult` | 0 | - | `[DA]` | | | - | | `nreal` | 0 | - | | `outdir` | \- | - | | `da_interval` | 1 | - | | `da_interval_final` | 1 | - | | `flexible_da_interval` | 0 | - | | `stat_dumpoffset` | 0 | - | | `screen_wrapper` | 1 | - | | `point_obs` | 1 | - | | `obs_interp_switch` | 0 | - | | `crns_flag` | 1 | - | | `da_crns_depth_tol` | 0.01 | - | | `print_obs_index` | 0 | + | section | parameter | default value | + |:----------|:-----------------------------|:--------------| + | `[PF]` | | | + | | `problemname` | \- | + | | `nprocs` | 0 | + | | `starttime` | 0.0 | + | | `endtime` | 0 | + | | `simtime` | 0 | + | | `dt` | 0.0 | + | | `updateflag` | 1 | + | | `paramupdate` | 0 | + | | `paramupdate_frequency` | 1 | + | | `dampingfactor_param` | 1.0 | + | | `dampingfactor_state` | 1.0 | + | | `damping_switch_sm` | 1 | + | | `aniso_perm_y` | 1.0 | + | | `aniso_perm_z` | 1.0 | + | | `aniso_use_parflow` | 0 | + | | | | + | | `printensemble` | 1 | + | | `t_printensemble` | -2 | + | | `printstat` | 1 | + | | `paramprintensemble` | 1 | + | | `paramprintstat` | 1 | + | | `olfmasking` | 0 | + | | `olfmasking_param` | 0 | + | | `olfmasking_depth` | 0 | + | `[CLM]` | | | + | | `problemname` | \- | + | | `nprocs` | 0 | + | | `update_swc` | 1 | + | | `update_snow` | 0 | + | | `update_snow_repartitioning` | 3 | + | | `print_swc` | 0 | + | | `print_et` | 0 | + | | `statevec_allcol` | 0 | + | | `statevec_colmean` | 0 | + | | `statevec_only_active` | 0 | + | | `statevec_max_layer` | 25 | + | | `t_printensemble` | -2 | + | | `watmin_switch` | 0 | + | `[COSMO]` | | | + | | `nprocs` | 0 | + | | `dtmult` | 0 | + | `[DA]` | | | + | | `nreal` | 0 | + | | `outdir` | \- | + | | `da_interval` | 1 | + | | `da_interval_final` | 1 | + | | `flexible_da_interval` | 0 | + | | `stat_dumpoffset` | 0 | + | | `screen_wrapper` | 1 | + | | `point_obs` | 1 | + | | `obs_interp_switch` | 0 | + | | `crns_flag` | 1 | + | | `da_crns_depth_tol` | 0.01 | + | | `print_obs_index` | 0 | Default values for parameter file `enkfpf.par`. diff --git a/docs/users_guide/running_tsmp_pdaf/snow_da.md b/docs/users_guide/running_tsmp_pdaf/snow_da.md new file mode 100644 index 00000000..8b6814db --- /dev/null +++ b/docs/users_guide/running_tsmp_pdaf/snow_da.md @@ -0,0 +1,193 @@ +(snowda)= +# Snow Data Assimilation # + +Snow data assimilation (Snow-DA) in TSMP-PDAF enables the assimilation +of snow-related observations—snow depth (SD) or snow water equivalent +(SWE)—into the eCLM land-surface model. + +## Configuration ## + +Snow-DA is controlled by two parameters in the [`[CLM]` section of +`enkfpf.par`](enkfpf:clm): + +- [`CLM:update_snow`](enkfpf:clm:update_snow): selects which snow + variable(s) are placed in the state vector and which CLM layer + variables are updated after the analysis step. +- [`CLM:update_snow_repartitioning`](enkfpf:clm:update_snow_repartitioning): + selects the method used to redistribute the updated bulk snow + quantity across the individual CLM snow layers. + +When Snow-DA is active (`CLM:update_snow != 0`), the snow state +variables replace the soil water content (SWC) in the state +vector. Simultaneous assimilation of SWC and snow in the same PDAF +update step is not supported. + +## State Vector ## + +The state vector size and content depend on `CLM:update_snow`. Snow +variables are always stored at the grid-cell level (one value per grid +cell); snow depth and SWE are taken from the first column of each grid +cell. + +| `update_snow` | State vector size | Variables in state vector | +|:--------------|:------------------|:-------------------------------------------| +| `1` | 1 × n_gridcells | SD (`snow_depth_col`) | +| `2` | 1 × n_gridcells | SWE (`h2osno_col`) | +| `3` | 2 × n_gridcells | SD (`snow_depth_col`) + SWE (`h2osno_col`) | +| `4` | 2 × n_gridcells | SD (`snow_depth_col`) + SWE (`h2osno_col`) | +| `5` | 2 × n_gridcells | SD (`snow_depth_col`) + SWE (`h2osno_col`) | +| `6` | 2 × n_gridcells | SD (`snow_depth_col`) + SWE (`h2osno_col`) | +| `7` | 2 × n_gridcells | SD (`snow_depth_col`) + SWE (`h2osno_col`) | + +For cases 3–7, SD occupies the first block of the state vector and SWE +the second block (offset by `clm_varsize`). + +## Snow Layer Variable Update ## + +After the PDAF analysis step updates the bulk snow quantities in the +state vector, the individual CLM snow layer variables must be adjusted +to remain physically consistent. The following layer-indexed arrays +(indexed over active snow layers `snlsno(j)+1 : 0`) may be updated: + +| CLM variable | Description | +|-------------------------------|----------------------------------------------| +| `h2osoi_ice_col` | Ice content per snow layer (kg m⁻²) | +| `h2osoi_liq_col` | Liquid water content per snow layer (kg m⁻²) | +| `dz` (column layer thickness) | Thickness of each snow layer (m) | + +The redistribution method is selected via +`CLM:update_snow_repartitioning`. + +### Repartitioning methods 1 and 2 ### + +Methods 1 and 2 are only available for `CLM:update_snow=1` or `2`. +They operate by computing the posterior SWE and distributing the +SWE gain/loss across the snow layers, then deriving layer thickness +changes from the local snow density. + +**Method 1** (DART-style): The entire SWE change is applied to the +bottom snow layer (`i=0`). All other layers are left unchanged. Layer +thickness is adjusted based on the local snow density. This approach +follows the DART/CLM implementation +(). + +**Method 2**: The SWE change is distributed across all active snow +layers proportionally to their current water content. Each layer +receives a fraction `(h2osoi_liq + h2osoi_ice) / SWE_prior` of the +total SWE change. + +For both methods, the posterior SWE (`h2osno_po`) is derived from: +- For `update_snow=1`: `h2osno_po = snow_depth_out × ρ_avg × frac_sno`, + where `ρ_avg = min(800, h2osno / snowdp)` is the layer-averaged snow + density. +- For `update_snow=2`: `h2osno_po = h2osno_out` directly from the state + vector. + +### Repartitioning method 3 (default) ### + +Method 3 uses a multiplicative increment to scale the layer +variables. The scaling factor is the ratio of the posterior to the +prior bulk snow value. It is available for all values of +`CLM:update_snow`. + +The following table summarises which layer variables are scaled and by +which increment for each `update_snow` option: + +| `update_snow` | `h2osoi_ice` scaled by | `h2osoi_liq` scaled by | `dz` scaled by | +|:---:|:---:|:---:|:---:| +| `1` | `SD_out / SD_in` | — | — | +| `2` | `SWE_out / SWE_in` | — | — | +| `3` | `SWE_out / SWE_in` | — | — | +| `4` | `SWE_out / SWE_in` | `SWE_out / SWE_in` | `SD_out / SD_in` | +| `5` | `SWE_out / SWE_in` | — | `SD_out / SD_in` | +| `6` | `SWE_out / SWE_in` | `SWE_out / SWE_in` | — | +| `7` | `SWE_out / SWE_in` | — | — | + +`SD_in` / `SD_out`: prior / posterior snow depth +`SWE_in` / `SWE_out`: prior / posterior snow water equivalent + +For `update_snow=1`, the column-level SWE (`h2osno`) is additionally +updated after the layer loop to match the adjusted layer ice contents. + +The increment update is only applied when both the prior and posterior +bulk snow values exceed `1e-6` to avoid numerical instabilities with +near-zero snow amounts. + +## Observation Files ## + +Snow-DA uses CLM-type observation files (see [Observation files](obs)). + +- For `CLM:update_snow=1` and `3`–`7`, the observation variable should + be snow depth (`SNOWDP` / `SNOW_DEPTH`). +- For `CLM:update_snow=2`, the observation variable should be SWE. + +Observations are expected at the grid-cell level (one observation per +grid cell, no layer index). In the TSMP-PDAF observation operator +(`init_dim_obs_f_pdaf.F90`), when `CLM:update_snow != 0`, the +observation index `obs_index_p` is set to the grid-cell index rather +than a layer-specific state vector index. + +## Safety Checks ## + +The following protective measures are applied during the Snow-DA update: + +- If the posterior SD (case 1) or SWE (case 2) is negative or zero, + a warning is printed to stdout and no update is applied for that + column. +- For cases 4–7, the increment update is skipped unless all four + quantities (prior SD, posterior SD, prior SWE, posterior SWE) + exceed `1e-6`. +- NaN checks are applied to all updated layer variables; a warning is + printed if NaN values are detected. +- The SWC masking flag [`CLM:swc_mask_snow`](enkfpf:clm) is + independent of Snow-DA and only affects SWC updates. + +## Configuration Examples ## + +### Assimilate snow depth (SD) only ### + +```text +[CLM] +update_snow = 1 +update_snow_repartitioning = 3 +``` + +The state vector contains one SD value per grid cell. After the PDAF +update, `snow_depth_col` is written back to CLM and `h2osoi_ice` in +each snow layer is scaled by `SD_out / SD_in`. + +### Assimilate SWE only ### + +```text +[CLM] +update_snow = 2 +update_snow_repartitioning = 3 +``` + +The state vector contains one SWE value per grid cell. After the PDAF +update, `h2osno_col` is written back to CLM and `h2osoi_ice` in each +snow layer is scaled by `SWE_out / SWE_in`. + +### Assimilate SD, update ice and layer thickness ### + +```text +[CLM] +update_snow = 5 +update_snow_repartitioning = 3 +``` + +Both SD and SWE are placed in the state vector. After the PDAF update, +`h2osoi_ice` is scaled by the SWE increment and layer thickness `dz` +is scaled by the SD increment. + +### Assimilate SD, update ice, liquid water, and layer thickness ### + +```text +[CLM] +update_snow = 4 +update_snow_repartitioning = 3 +``` + +Both SD and SWE are placed in the state vector. After the PDAF update, +`h2osoi_ice` and `h2osoi_liq` are scaled by the SWE increment, and +layer thickness `dz` is scaled by the SD increment. diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 11547a4e..46ff67a7 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -145,6 +145,7 @@ 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_snow !hcp end #endif #endif @@ -987,7 +988,12 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) if(clmstatevec_only_active==1) then obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) else - obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + if(clmupdate_snow/=0) then + ! Snow-DA: no layer in state vector variables + obs_index_p(cnt) = g-begg+1 + else + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if end if #else obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 48f29c08..d47dc52b 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -141,6 +141,7 @@ 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_snow !hcp end #endif #endif @@ -980,7 +981,12 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) if(clmstatevec_only_active==1) then obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) else - obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + if(clmupdate_snow/=0) then + ! Snow-DA: no layer in state vector variables + obs_index_p(cnt) = g-begg+1 + else + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if end if #else obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index 7462a6f1..36960080 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -67,6 +67,8 @@ module enkf_clm_mod integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture + integer(c_int),bind(C,name="clmupdate_snow") :: clmupdate_snow + integer(c_int),bind(C,name="clmupdate_snow_repartitioning") :: clmupdate_snow_repartitioning integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et @@ -149,6 +151,25 @@ subroutine define_clm_statevec(mype) error stop "Not implemented: clmupdate_texture.eq.2" endif + ! Snow assimilation + ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5 + ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell + ! Case 2: Assimilation of snow water equivalent same as above + if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then + error stop "Not implemented: clmupdate_snow.eq.1 or clmupdate_snow.eq.1" + endif + ! Case 3: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice + if(clmupdate_snow.eq.3) then + error stop "Not implemented: clmupdate_snow.eq.3" + endif + ! Case 4: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq + ! and dz + if(clmupdate_snow.eq.4) then + error stop "Not implemented: clmupdate_snow.eq.4" + endif + !hcp LST DA if(clmupdate_T.eq.1) then clm_varsize = endg-begg+1 @@ -289,6 +310,14 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + ! Snow assimilation state vector + ! Case 1: Snow depth + ! Case 2: SWE + ! Case 3: Snow depth + SWE + if(clmupdate_snow.ne.0) then + error stop "Not implemented: clmupdate_snow.ne.0" + endif + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN ! TSMP-PDAF: For debug runs, output the state vector in files @@ -536,6 +565,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call clm_texture_to_parameters endif + ! Snow assimilation: + ! Case 1: Snow depth + ! Case 2: Snow water equivalent + ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector + ! Write updated snow variable back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow.ne.0) then + error stop "Not implemented: clmupdate_snow.ne.0" + endif + end subroutine update_clm subroutine clm_correct_texture() diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h index f943ad90..e8fac043 100755 --- a/interface/model/common/enkf.h +++ b/interface/model/common/enkf.h @@ -88,6 +88,8 @@ GLOBAL int nx_local,ny_local,nz_local; GLOBAL int clmupdate_swc; GLOBAL int clmupdate_T; GLOBAL int clmupdate_texture; +GLOBAL int clmupdate_snow; +GLOBAL int clmupdate_snow_repartitioning; GLOBAL int clmprint_swc; GLOBAL int clmprint_et; GLOBAL int clmstatevec_allcol; diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c index e654b419..01ac39af 100755 --- a/interface/model/common/read_enkfpar.c +++ b/interface/model/common/read_enkfpar.c @@ -79,6 +79,8 @@ void read_enkfpar(char *parname) clmupdate_swc = iniparser_getint(pardict,"CLM:update_swc",1); clmupdate_T = iniparser_getint(pardict,"CLM:update_T",0); clmupdate_texture = iniparser_getint(pardict,"CLM:update_texture",0); + clmupdate_snow = iniparser_getint(pardict,"CLM:update_snow",0); + clmupdate_snow_repartitioning = iniparser_getint(pardict,"CLM:update_snow_repartitioning",3); clmprint_swc = iniparser_getint(pardict,"CLM:print_swc",0); clmprint_et = iniparser_getint(pardict,"CLM:print_et",0); clmstatevec_allcol = iniparser_getint(pardict,"CLM:statevec_allcol",0); diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 1e43450c..b2544b27 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -55,6 +55,8 @@ module enkf_clm_mod integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture + integer(c_int),bind(C,name="clmupdate_snow") :: clmupdate_snow + integer(c_int),bind(C,name="clmupdate_snow_repartitioning") :: clmupdate_snow_repartitioning integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et @@ -138,6 +140,48 @@ subroutine define_clm_statevec(mype) clm_statevecsize = clm_statevecsize + 3*((endg-begg+1)*nlevsoi) end if + ! Snow assimilation + ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5 + ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell + ! Case 2: Assimilation of snow water equivalent same as above + if(clmupdate_snow==1 .or. clmupdate_snow==2) then + clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA + clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority + endif + ! Case 3: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice + if(clmupdate_snow==3) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif + ! Case 4: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq + ! and dz + if(clmupdate_snow==4) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif + ! Case 5: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice + ! and dz + if(clmupdate_snow==5) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif + ! Case 6: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq + if(clmupdate_snow==6) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif + ! Case 7: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice + ! Should reproduce case 3 + if(clmupdate_snow==7) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif + !hcp LST DA if(clmupdate_T==1) then error stop "Not implemented: clmupdate_T.eq.1" @@ -394,6 +438,8 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: snow_depth(:) + real(r8), pointer :: h2osno(:) integer :: i,j,jj,g,c,cc,offset integer :: n_c character (len = 34) :: fn !TSMP-PDAF: function name for state vector output @@ -407,6 +453,9 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN @@ -457,6 +506,53 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + ! Snow assimilation state vector + ! Case 1: Snow depth + ! Case 2: SWE + ! Case 3: Snow depth + SWE + if(clmupdate_snow/=0) then + cc = 1 + do j=clm_begg,clm_endg + ! Only get the snow_depth/swe from the first column of each gridcell + ! and add it to the clm_statevec at the position of the gridcell (cc) + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g == j) then + if (newgridcell) then + newgridcell = .false. + + if(clmupdate_snow==1) then + clm_statevec(cc+offset) = snow_depth(jj) + else if(clmupdate_snow==2) then + clm_statevec(cc+offset) = h2osno(jj) + else if(clmupdate_snow==3) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow==4) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow==5) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow==6) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow==7) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else + error stop "Wrong input for clmupdate_snow" + end if + + endif + endif + end do + cc = cc + 1 + end do + + endif + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN ! TSMP-PDAF: For debug runs, output the state vector in files @@ -540,6 +636,7 @@ end subroutine set_clm_statevec_swc subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 + use ColumnType , only : col use clm_instMod, only : waterstate_inst implicit none @@ -547,16 +644,36 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") integer,intent(in) :: tstartcycle integer,intent(in) :: mype + real(r8), pointer :: snow_depth(:) + real(r8), pointer :: h2osno(:) + + real(r8), pointer :: dz(:,:) ! layer thickness depth (m) + real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) + real(r8) :: incr_sno + real(r8) :: incr_sd + real(r8) :: incr_swe + real(r8) :: h2osno_in(clm_begc:clm_endc) + real(r8) :: snow_depth_in(clm_begc:clm_endc) + real(r8) :: h2osno_out(clm_begc:clm_endc) + real(r8) :: snow_depth_out(clm_begc:clm_endc) + integer :: i + integer :: j + integer :: cc + integer :: offset character (len = 31) :: fn !TSMP-PDAF: function name for state vector output + character (len = 32) :: fn4 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn5 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn6 !TSMP-PDAF: function name for state vector outpu + integer, pointer :: snlsno(:) logical :: swc_zero_before_update + cc = 0 + offset = 0 swc_zero_before_update = .false. #ifdef PDAF_DEBUG @@ -571,8 +688,12 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") END IF #endif + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) + h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col + snlsno => col%snl ! number of snow layers (negative) #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN @@ -583,7 +704,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") OPEN(unit=71, file=fn5, action="write") WRITE (71,"(es22.15)") h2osoi_liq(:,:) CLOSE(71) + END IF + IF(clmupdate_swc/=0 .OR. clmupdate_snow/=0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn6, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".bef_up.", tstartcycle, ".txt" OPEN(unit=71, file=fn6, action="write") @@ -620,6 +743,190 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call update_clm_texture(tstartcycle, mype) endif + ! Snow assimilation: + ! Case 1: Snow depth + ! Case 2: Snow water equivalent + ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector + ! Write updated snow variable back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow/=0) then + do j=clm_begc,clm_endc + ! Iterate through the columns + + ! For each column index, copy the gridcell-linked state vector + ! value at state vector index `cc` + + ! Set cc (the state vector index) from the + ! CLM5-grid-index + cc = (col%gridcell(j) - clm_begg + 1) + + if (clmupdate_snow==1) then + snow_depth_in(j) = snow_depth(j) + else if (clmupdate_snow==2) then + h2osno_in(j) = h2osno(j) + else if (clmupdate_snow==3) then + h2osno_in(j) = h2osno(j) + else if (clmupdate_snow==4) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + else if (clmupdate_snow==5) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + else if (clmupdate_snow==6) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + else if (clmupdate_snow==7) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + end if + + if (clmupdate_snow==1) then + snow_depth_out(j) = clm_statevec(cc+offset) + else if (clmupdate_snow==2) then + h2osno_out(j) = clm_statevec(cc+offset) + else if (clmupdate_snow==3) then + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow==4) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow==5) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow==6) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow==7) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + end if + + if(clmupdate_snow==1) then + ! Update state variable to CLM + ! Needed for Case 1/2 if they use repartioning function + if (snow_depth_out(j)>0.0) then + snow_depth(j) = snow_depth_out(j) + else + ! Catch negative or 0 values from DA + print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", & + cc, j, offset, snow_depth_out(j) + end if + else if(clmupdate_snow==2) then + if (h2osno_out(j)>0.0) then + h2osno(j) = h2osno_out(j) + else + ! Catch negative or 0 values from DA + print *, "WARNING: SWE is negative/zero at cc. cc, j, offset, h2osno(j): ", cc, j, offset, h2osno(j) + end if + end if + + end do + + ! Repartitioning + if ( clmupdate_snow_repartitioning==1 .or. clmupdate_snow_repartitioning==2) then + + if (clmupdate_snow==1) then + if ( SUM(ABS(snow_depth_in(:) - snow_depth_out(:)))>0.000001) then + call clm_repartition_snow(snow_depth_in(:)) + end if + end if + + if (clmupdate_snow==2) then + if ( SUM(ABS(h2osno_in(:) - h2osno_out(:)))>0.000001) then + call clm_repartition_snow(h2osno_in(:)) + end if + end if + + if (clmupdate_snow>=3) then + error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.ge.3" + end if + + end if + + if ( clmupdate_snow_repartitioning==3) then + + if (clmupdate_snow==1) then + do j=clm_begc,clm_endc + if (snow_depth_out(j)>1.0d-6) then + if (snow_depth_in(j)>1.0d-6) then + ! Update h2osoi_ice with increment + incr_sno = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno + end do + end if + end if + end do + end if + + if (clmupdate_snow==2 .or. clmupdate_snow==3) then + do j=clm_begc,clm_endc + if (h2osno_out(j)>1.0d-6) then + if (h2osno_in(j)>1.0d-6) then + ! Update h2osoi_ice with increment + incr_sno = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno + end do + end if + end if + end do + end if + + if (clmupdate_snow==4 .or. clmupdate_snow==5 .or. clmupdate_snow==6 .or. clmupdate_snow==7) then + do j=clm_begc,clm_endc + if (h2osno_out(j)>1.0d-6) then + if (h2osno_in(j)>1.0d-6) then + if (snow_depth_out(j)>1.0d-6) then + if (snow_depth_in(j)>1.0d-6) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + ! Update snow_depth with increment + incr_sd = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + + if (clmupdate_snow==4 .or. clmupdate_snow==5 .or. clmupdate_snow==6 .or. clmupdate_snow==7) then + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + end if + + if (clmupdate_snow==4 .or. clmupdate_snow==6) then + h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_liq(j,i))) then + print *, "WARNING: h2osoi_liq at j,i is nan: ", j, i + endif + end if + + if (clmupdate_snow==4 .or. clmupdate_snow==5) then + dz(j,i) = dz(j,i) * incr_sd + if (isnan(dz(j,i))) then + print *, "WARNING: dz at j,i is nan: ", j, i + endif + end if + + end do + end if + end if + end if + end if + end do + end if + + end if + +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) 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) + + END IF +#endif + endif + end subroutine update_clm @@ -841,6 +1148,195 @@ subroutine update_clm_texture(tstartcycle, mype) end subroutine update_clm_texture + subroutine clm_repartition_snow(h2osno_in) + use ColumnType, only : col + use clm_instMod, only : waterstate_inst + use clm_varpar, only : nlevsno, nlevsoi + use clm_varcon, only : bdsno, denice + use shr_kind_mod, only : r8 => shr_kind_r8 + + implicit none + + real(r8), intent(in), optional :: h2osno_in(clm_begc:clm_endc) + + real(r8), pointer :: snow_depth(:) + real(r8), pointer :: h2osno(:) + real(r8), pointer :: h2oliq(:,:) + real(r8), pointer :: h2oice(:,:) + real(r8), pointer :: frac_sno(:) + real(r8), pointer :: snowdp(:) + integer, pointer :: snlsno(:) + + real(r8) :: dzsno(clm_begc:clm_endc,-nlevsno+1:0) + real(r8) :: h2osno_po(clm_begc:clm_endc) + real(r8) :: h2osno_pr(clm_begc:clm_endc) + real(r8) :: snowden, frac_swe, frac_liq, frac_ice + real(r8) :: gain_h2osno, gain_h2oliq, gain_h2oice + real(r8) :: gain_dzsno(-nlevsno+1:0) + real(r8) :: rho_avg, z_avg + integer :: i,ii,j,jj,g,cc,offset + + cc = 1 + offset = 0 + + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + snowdp => waterstate_inst%snowdp_col ! area-averaged snow height (m) + h2osno => waterstate_inst%h2osno_col ! col snow water (mm H2O) + h2oliq => waterstate_inst%h2osoi_liq_col ! col liquid water (kg/m2) (-nlevsno+1:nlevgrnd) + h2oice => waterstate_inst%h2osoi_ice_col ! col ice lens (kg/m2) (-nlevsno+1:nlevgrnd) + + snlsno => col%snl ! number of snow layers (negative) + + frac_sno => waterstate_inst%frac_sno_eff_col ! fraction of ground covered by snow + ! dz for snow layers is defined like in the history output as col%dz for the snow layers + dzsno(clm_begc:clm_endc, -nlevsno+1:0) = col%dz(clm_begc:clm_endc,-nlevsno+1:0) + ! Iterate through all columns + do jj=clm_begc,clm_endc + if (h2osno(jj)<0.0) then ! No snow in column + print *, "WARNING: negative snow in col: ", jj, h2osno +! ! Set existing layers to near zero and let CLM do the layer aggregation + do i=0,snlsno(jj)+1,-1 + h2oliq(jj,i) = 0.0_r8 + h2oice(jj,i) = 0.00000001_r8 + dzsno(jj,i) = 0.00000001_r8 + end do + snow_depth(jj) = sum(dzsno(jj,:)) + h2osno(jj) = sum(h2oice(jj,:)) + + if (isnan(h2osno(jj))) then + print *, "WARNING: h2osno at jj is nan: ", jj + endif + if (isnan(snow_depth(jj))) then + print *, "WARNING: snow_depth at jj is nan: ", jj + endif + + else ! snow (h2osno) present + if (snlsno(jj)<0) then ! snow layers in the column + if (clmupdate_snow == 1) then + ! DART source: https://github.com/NCAR/DART/blob/main/models/clm/dart_to_clm.f90 + ! Formulas below from DART use h2osno_po / h2osno_pr for after / before DA SWE + ! clmupdate_snow == 1 has snow_depth after and h2osno before DA snow depth + ! Therefore need to have a transform to get h2osno_po + ! v1 init + ! h2osno_po(jj) = (snow_depth(jj) * bdsno) ! calculations from Init using constant SBD + ! v2 SoilTemperatureMod + if (snowdp(jj)>0.0_r8) then + rho_avg = min(800.0_r8, h2osno(jj)/snowdp(jj)) + else + rho_avg = 200.0_r8 + end if + if (frac_sno(jj)>0.0_r8 .and. snlsno(jj)<0.0_r8) then + h2osno_po(jj) = snow_depth(jj) * (rho_avg*frac_sno(jj)) + else + h2osno_po(jj) = snow_depth(jj) * denice + end if + h2osno_pr(jj) = h2osno(jj) + else if (clmupdate_snow == 2) then + ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already + h2osno_po(jj) = h2osno(jj) + h2osno_pr(jj) = h2osno_in(jj) + end if + + do ii=0,snlsno(jj)+1,-1 ! iterate through the snow layers + ! DART VERSION: ii = nlevsoi - i + 1 + ! snow density prior for each layer + if (dzsno(jj,ii)>0.0_r8) then + snowden = (h2oliq(jj,ii) + h2oice(jj,ii) / dzsno(jj,ii)) + else + snowden = 0.0_r8 + endif + ! fraction of SWE in each active layers + if(h2osno_pr(jj)>0.0_r8) then + ! repartition Option 1: Adjust bottom layer only (set weight to 1 for bottom 0 else) + if(clmupdate_snow_repartitioning==1) then + if (ii == 0) then ! DART version indexing check against nlevsno but for us 0 + frac_swe = 1.0_r8 + ! JK: Let CLM repartitioning do the job + ! afterwards. Provide all the snow in a single layer + else + frac_swe = 0.0_r8 + end if + ! repartition Option 2: Adjust all active layers + else if (clmupdate_snow_repartitioning==2) then + frac_swe = (h2oliq(jj,ii) + h2oice(jj,ii)) / h2osno_pr(jj) + end if + else + frac_swe = 0.0_r8 ! no fraction SWE if no snow is present in column + end if ! end SWE fraction if + + ! fraction of liquid and ice + if ((h2oliq(jj,ii) + h2oice(jj,ii))>0.0_r8) then + frac_liq = h2oliq(jj,ii) / (h2oliq(jj,ii) + h2oice(jj,ii)) + frac_ice = 1.0_r8 - frac_liq + else + frac_liq = 0.0_r8 + frac_ice = 0.0_r8 + end if + + ! SWE adjustment per layer + ! assumes identical layer distribution of liq and ice than before DA (frac_*) + if (abs(h2osno_po(jj) - h2osno_pr(jj))>0.0_r8) then + gain_h2osno = (h2osno_po(jj) - h2osno_pr(jj)) * frac_swe + gain_h2oliq = gain_h2osno * frac_liq + gain_h2oice = gain_h2osno * frac_ice + else + gain_h2osno = 0.0_r8 + gain_h2oliq = 0.0_r8 + gain_h2oice = 0.0_r8 + end if + ! layer level adjustments + if (snowden>0.0_r8) then + gain_dzsno(ii) = gain_h2osno / snowden + else + gain_dzsno(ii) = 0.0_r8 + end if + h2oliq(jj,ii) = h2oliq(jj,ii) + gain_h2oliq + h2oice(jj,ii) = h2oice(jj,ii) + gain_h2oice + + ! Adjust snow layer dimensions so that CLM5 can calculate compaction / aggregation + ! in the DART code dzsno is adjusted directly but in CLM5 dzsno is local and diagnostic + ! i.e. calculated / assigned from frac_sno and dz(:, snow_layer) in SnowHydrologyMod + ! therefore we adjust dz(:, snow_layer) here + if (abs(h2osno_po(jj) - h2osno_pr(jj))>0.0_r8) then + col%dz(jj, ii) = col%dz(jj, ii) + gain_dzsno(ii) + ! mid point and interface adjustments + ! i.e. zsno (col%z(:, snow_layers)) and zisno (col%zi(:, snow_layers)) + ! DART version the sum goes from ilevel:nlevsno to fit with our indexing: + col%zi(jj, ii) = sum(col%dz(jj,ii:0))*-1.0_r8 + ! In DART the check is ilevel == nlevsno but here + if (ii==0) then + col%z(jj,ii) = col%zi(jj,ii) / 2.0_r8 + else + col%z(jj,ii) = sum(col%zi(jj, ii:ii+1)) / 2.0_r8 + end if + end if + + end do ! End iteration of snow layers + endif ! End of snow layers present check + + ! Column level variables + if (clmupdate_snow == 1) then + ! Finally adjust SWE (h2osno) since the prior value is no longer needed + ! column level variables so we can adjust it outside the layer loop + h2osno(jj) = h2osno_po(jj) + else if (clmupdate_snow == 2 .or. clmupdate_snow == 3) then + ! Update the total snow depth to match udpates to layers for active snow layers + if (abs(h2osno_po(jj) - h2osno_pr(jj)) > 0.0_r8 .and. snlsno(jj) < 0.0_r8) then + snow_depth(jj) = snow_depth(jj) + sum(gain_dzsno(-nlevsno+1:0)) + end if + end if + + if (isnan(h2osno(jj))) then + print *, "WARNING2: h2osno at jj is nan: ", jj + endif + if (isnan(snow_depth(jj))) then + print *, "WARNING2: snow_depth at jj is nan: ", jj + endif + + end if ! End of snow present check + end do ! End of column iteration + + end subroutine clm_repartition_snow subroutine clm_correct_texture() diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c index 56687e03..19e2dc2c 100644 --- a/interface/model/wrapper_tsmp.c +++ b/interface/model/wrapper_tsmp.c @@ -198,7 +198,7 @@ void integrate_tsmp() { void update_tsmp(){ #if defined CLMSA - if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0))){ + if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_snow!=0))){ update_clm(&tstartcycle, &mype_world); if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2){ print_update_clm(&tcycle, &total_steps);