Skip to content

Commit 9545d4e

Browse files
jjokellajjokella
committed
Merge branch 'tsmp-pdaf-patched' into tsmp-pdaf-patched-eclm-pfl-fixes
2 parents 6e54928 + fca3bd5 commit 9545d4e

34 files changed

Lines changed: 6332 additions & 59 deletions

‎CONTRIBUTORS.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,5 +6,6 @@ Wolfgang Kurtz
66
Stefan Poll
77
Mukund Pondkule
88
Prabhakar Shrestha
9+
Anne Springer
910
Lukas Strebel
1011

‎Makefile‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -138,7 +138,8 @@ SRC_PDAF_GEN = PDAF_analysis_utils.F90 \
138138
PDAFlocalomi_put_state_nondiagR.F90 \
139139
PDAFlocalomi_put_state_nondiagR_si.F90 \
140140
PDAFlocalomi_put_state_si.F90 \
141-
PDAF_correlation_function.F90
141+
PDAF_correlation_function.F90 \
142+
PDAF_reset_dim_p.F90
142143

143144
# Specific PDAF-routines for SEIK
144145
SRC_SEIK = PDAF_seik_init.F90 \

‎docs/_toc.yml‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,8 @@ parts:
2020
- file: users_guide/running_tsmp_pdaf/input_cmd
2121
- file: users_guide/running_tsmp_pdaf/input_obs
2222
- file: users_guide/running_tsmp_pdaf/input_enkfpf
23+
- file: users_guide/running_tsmp_pdaf/tsmp_pdaf_omi
24+
- file: users_guide/running_tsmp_pdaf/grace_tws_da
2325

2426
- file: users_guide/misc/README
2527
title: Misc
@@ -40,6 +42,8 @@ parts:
4042

4143
- file: developers_guide/how_to_use_this_document
4244

45+
- file: developers_guide/omi_developer_guide
46+
4347
- caption: Reference
4448
chapters:
4549
# - file: reference/history_fields

‎docs/developers_guide/omi_developer_guide.md‎

Lines changed: 262 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,262 @@
1+
(omi:developer-guide)=
2+
# Developer Guide: Adding a New Observation Type to the OMI Interface
3+
4+
This page describes how to add a new observation type — generically
5+
called C below — to the TSMP-PDAF OMI (Observation Model Interface)
6+
framework. It covers the architecture, the role of each source file,
7+
and the changes required, using the existing GRACE-DA and SM-DA
8+
implementations as reference.
9+
10+
## Architecture Overview
11+
12+
Observation handling sits between the PDAF library and the eCLM model
13+
state in three layers:
14+
15+
```
16+
PDAF library (src/)
17+
│ calls callback routines at each analysis step
18+
▼
19+
callback_obs_pdafomi.F90 ← single routing hub
20+
│ dispatches to one routine per observation type
21+
▼
22+
obs_GRACE_pdafomi.F90 ← one module per observation type
23+
obs_SM_pdafomi.F90
24+
obs_C_pdafomi.F90 (new)
25+
│ read observations, build H, apply H to state vector
26+
▼
27+
eCLM state vector (via enkf_clm_mod, mod_assimilation, …)
28+
```
29+
30+
Each observation-type module is self-contained. The only coupling
31+
between modules is through the central routing file and the shared
32+
`mod_assimilation` module.
33+
34+
## Files Involved
35+
36+
To add a new observation type `C`, touch exactly **eight** things:
37+
38+
| File | Action |
39+
|------------------------------------------------|-----------------------------------------------------------------------|
40+
| `interface/framework/obs_C_pdafomi.F90` | **Create** (new file, modelled on the template) |
41+
| `interface/framework/Makefile` | Add `obs_C_pdafomi.o` to `MOD_USER_PDAFOMI` |
42+
| `interface/framework/callback_obs_pdafomi.F90` | Add calls for C to every callback routine |
43+
| `interface/framework/mod_assimilation.F90` | Declare `cradius_C`, `sradius_C`, and any other C-specific parameters |
44+
| `interface/framework/mod_read_obs.F90` | Add `'C'` case to `update_obs_type` |
45+
| `interface/model/eclm/enkf_clm_mod_5.F90` | Add `clmupdate_C` flag (parallel to `clmupdate_tws`, `clmupdate_swc`) |
46+
| `interface/framework/init_pdaf.F90` | Import `assim_C` and set it from `clmupdate_C` |
47+
| `interface/framework/init_pdaf_parse.F90` | Import `rms_obs_C` and add a `parse` call for it |
48+
49+
The template at `templates/omi/obs_OBSTYPE_pdafomi_TEMPLATE.F90` is a
50+
good starting skeleton. The existing SM and GRACE modules show how the
51+
template is applied in this codebase.
52+
53+
---
54+
55+
## Step 1 — Create `obs_C_pdafomi.F90`
56+
57+
### Module-level variables
58+
59+
At minimum, declare `assim_C` (a logical switch controlled by
60+
`callback_obs_pdafomi`), `rms_obs_C` (a constant fallback error
61+
standard deviation), and the PDAFomi types `thisobs` (of type `obs_f`,
62+
thread-safe) and `thisobs_l` (of type `obs_l`, declared
63+
`THREADPRIVATE`).
64+
65+
`assim_C` and `rms_obs_C` are wired up externally:
66+
- **`init_pdaf.F90`** imports `assim_C` and sets it from `clmupdate_C`
67+
(i.e. `assim_C = (clmupdate_C /= 0)`), so the observation type is
68+
enabled automatically when the corresponding model-state update flag
69+
is active. Commented-out placeholders mark the insertion point.
70+
- **`init_pdaf_parse.F90`** imports `rms_obs_C` and populates it via a
71+
`parse` call that reads `rms_obs_C` from `enkfpf.par`. Commented-out
72+
placeholders mark the insertion point here as well.
73+
74+
Additional module-level allocatable arrays — e.g. a cached
75+
climatological mean needed for an anomaly operator — follow the same
76+
pattern as in `obs_GRACE_pdafomi.F90`.
77+
78+
### The `obs_f` type
79+
80+
`obs_f` carries all observation metadata and is populated inside
81+
`init_dim_obs_C`. Mandatory fields:
82+
83+
| Field | Type | Meaning |
84+
|-------------------------|------------------------|-------------------------------------------------------|
85+
| `thisobs%doassim` | `INTEGER` | 1 = assimilate, 0 = skip |
86+
| `thisobs%disttype` | `INTEGER` | distance metric for localization (see below) |
87+
| `thisobs%ncoord` | `INTEGER` | number of spatial coordinates (usually 2) |
88+
| `thisobs%id_obs_p(:,:)` | `INTEGER, ALLOCATABLE` | state-vector index for each process-local observation |
89+
90+
Frequently useful optional fields:
91+
- `icoeff_p`: bilinear interpolation weights
92+
- `obs_err_type`: 0 = Gaussian, 1 = Laplace
93+
- `inno_omit`: innovation outlier rejection threshold
94+
95+
**`disttype` choices:**
96+
- `0` — Cartesian distance; coordinates in any consistent unit (GRACE
97+
uses integer grid-cell indices, so `cradius_GRACE` is in number of
98+
grid cells).
99+
- `3` — geographic distance; coordinates in degrees, radii in km (used
100+
by SM).
101+
102+
Use `disttype=3` for point observations distributed geographically (C,
103+
SM) and `disttype=0` for products already mapped to grid-cell indices
104+
(e.g. GRACE).
105+
106+
---
107+
108+
### Routine 1 — `init_dim_obs_C`
109+
110+
Called once per analysis step before the filter loop. Responsibilities:
111+
112+
1. Set the mandatory `thisobs` fields.
113+
2. Read the global observation vector from the NetCDF file using
114+
`read_obs_nc_type` from `mod_read_obs`.
115+
3. Handle the `dim_obs == 0` case (no observations scheduled for this
116+
step): allocate dummy arrays of size 1 and call
117+
`PDAFomi_gather_obs` before returning.
118+
4. For each global observation that falls within the PE-local domain,
119+
record the state-vector index in `thisobs%id_obs_p` and populate
120+
the local arrays `obs_p`, `ivar_obs_p`, and `ocoord_p`.
121+
5. Call `PDAFomi_gather_obs` to finalize observation bookkeeping across
122+
MPI ranks.
123+
124+
**Reading observations:** Both SM and GRACE use `read_obs_nc_type`,
125+
which checks the `type_clm` NetCDF variable and returns `dim_obs = 0`
126+
if the file does not contain observations of the requested type.
127+
128+
**Who reads the file:** SM reads only on `mype_filter==0` and
129+
broadcasts; GRACE reads on all filter processes because the spatial
130+
averaging operator needs each PE to know its contribution. For point
131+
observations like C, the SM pattern (read + broadcast) is simpler.
132+
133+
**Snapping observations to the state vector:** For point observations
134+
(SM, C), loop over all global observations and all PE-local gridcells
135+
and record the state-vector index for matching pairs. For GRACE the
136+
logic is inverted: each observation aggregates multiple gridcells, so
137+
the loop iterates over gridcells and maps each one to the observation
138+
it contributes to.
139+
140+
---
141+
142+
### Routine 2 — `obs_op_C`
143+
144+
Applies the observation operator H, mapping the PE-local model state to
145+
the observation space. For a simple point observation (SM pattern):
146+
allocate a local output array, extract `state_p(obs_index_p(i))` for
147+
each local observation, then call `PDAFomi_gather_obsstate` to assemble
148+
the global result in `ostate`.
149+
150+
For aggregating operators (GRACE pattern), sum all gridcell
151+
contributions mapping to the same observation, apply any transformation
152+
(e.g. subtract temporal mean for TWS anomaly), and then call
153+
`PDAFomi_gather_obsstate`. Follow the SM pattern for C unless your
154+
observations represent spatial averages over a large footprint.
155+
156+
---
157+
158+
### Routines 3 & 4 — Localization (required for LESTKF/LETKF)
159+
160+
**`init_dim_obs_l_C`** counts observations within the localization
161+
radius of the current local analysis domain. It is a thin wrapper
162+
around `PDAFomi_init_dim_obs_l`, passing `thisobs_l`, `thisobs`, the
163+
local domain coordinates `coords_l` (in the same units and `disttype`
164+
as `ocoord_p`), and the localization parameters `cradius_C` and
165+
`sradius_C`.
166+
167+
**`localize_covar_C`** is required only for the LEnKF. It is an
168+
equally thin wrapper around `PDAFomi_localize_covar` with the same
169+
localization parameters.
170+
171+
---
172+
173+
## Step 2 — Update `callback_obs_pdafomi.F90`
174+
175+
This file dispatches PDAF callbacks to each observation module. Add C
176+
in the same pattern as GRACE and SM. The commented-out
177+
`!USE obs_C_pdafomi` placeholders show exactly where to insert new
178+
entries.
179+
180+
For each of the eight callback subroutines, add a `USE obs_C_pdafomi`
181+
statement and one call to the corresponding C routine:
182+
183+
- **`init_dim_obs_pdafomi`**: set `assim_C = .true.`, call
184+
`init_dim_obs_C`, and add `dim_obs_C` to the running total. The
185+
`assim_*` flags exist to disable an observation type for debugging;
186+
keep `assim_C = .true.` here because the module itself returns
187+
`dim_obs = 0` when no observations are present. (This logic is
188+
planned to be complemented in the future by setting `assim_C` from
189+
the EnKF input file)
190+
- **`obs_op_pdafomi`**: call `obs_op_C`. Call order does not affect
191+
correctness because offsets in `ostate` are determined by the order
192+
of `init_dim_obs_*` calls above.
193+
- **`init_dim_obs_l_pdafomi`**: call `init_dim_obs_l_C`. The variable
194+
`dim_obs_l` is incremented inside each call; PDAF initializes it to
195+
0 before invoking the callback.
196+
- **`localize_covar_pdafomi`**: call `localize_covar_C`.
197+
- **Remaining callbacks** (`add_obs_err`, `init_obscovar`, `prodRinvA`,
198+
`prodRinvA_l`, `deallocate_obs`): each needs a corresponding C
199+
call, all of which are thin wrappers around the matching PDAFomi
200+
generic routine, identical in structure to the existing GRACE and SM
201+
entries.
202+
203+
---
204+
205+
## Step 3 — Add Parameters to `mod_assimilation.F90`
206+
207+
Declare `cradius_C` and `sradius_C` next to the existing GRACE and
208+
SM radius parameters. These are read from the `enkfpf.par` input file
209+
via the existing namelist mechanism; add corresponding entries to the
210+
`[DA]` section parser in the same module.
211+
212+
---
213+
214+
## Step 4 — Register the New Type in `mod_read_obs.F90`
215+
216+
The routine `update_obs_type` reads the `type_clm` string from the
217+
observation file and sets the `clmupdate_*` flags for the next
218+
analysis cycle. Add a `'C'` case that sets the appropriate flags. If
219+
no existing `clmupdate_*` flag covers your new variable, declare a new
220+
`clmupdate_C` in `enkf_clm_mod.F90` and wire it to the state-vector
221+
fill and scatter logic in the eCLM interface routines (`eclm/`).
222+
223+
---
224+
225+
## Key Design Differences Between GRACE-DA and SM-DA
226+
227+
### Localization distance units
228+
229+
| Type | `disttype` | Coordinate arrays | Radius unit |
230+
|-------|---------------|----------------------------|-------------|
231+
| GRACE | 0 (Cartesian) | integer grid-cell indices | grid cells |
232+
| SM | 3 (haversine) | degrees longitude/latitude | km |
233+
234+
### Observation operator complexity
235+
236+
- **SM**: direct extraction from `state_p` using a precomputed index.
237+
- **GRACE**: aggregating — all gridcell values within the GRACE
238+
footprint are summed and averaged, then the temporal mean is
239+
subtracted to produce a TWS anomaly. `thisobs%id_obs_p(1, g)` maps
240+
each gridcell `g` to its observation (inverse of the usual
241+
convention).
242+
243+
Follow the SM pattern for C point observations; follow the GRACE
244+
pattern if observations represent spatial averages over a large
245+
footprint.
246+
247+
### Coverage checks and error specification
248+
249+
GRACE discards observations where fewer than a minimum number of
250+
gridcells lie within the support radius (important near coastlines). SM
251+
and C point observations do not need this, but consider rejecting
252+
observations over water bodies or outside the CLM domain.
253+
254+
Error specification is controlled by `multierr`: constant
255+
(`multierr=0`), per-observation from `obserr_clm` (`multierr=1`), or
256+
full covariance matrix from `obscov_clm` (`multierr=2`).
257+
258+
---
259+
260+
## Assimilating Multiple Observation Types at the Same Timestep
261+
262+
The framework generates a state vector for each type individually before the assimilation, some things would need to be adapted when mutliple observation types are assimilated at the same timestep. Currently, one observation file only consists of one observation type. As SM observations are usually assimilated at noon and GRACE observations are assimilated at the end of the month at midnight, this should not provide any problems.

‎docs/users_guide/running_tsmp_pdaf/README.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,9 @@ COSMO](cos)). Additionally, a control file for the data assimilation
1515
Furthermore, some command line options ([Command line options](cmd))
1616
need to be specified when TSMP-PDAF is executed.
1717

18+
For assimilation of GRACE/GRACE-FO terrestrial water storage
19+
observations into eCLM, see [GRACE/TWS Data Assimilation](gracetws).
20+
1821
See the Virtual Machine download on webpage
1922
<https://datapub.fz-juelich.de/slts/tsmp-vm/index.html>.
2023

0 commit comments

Comments
 (0)