Weekly TSMP-PDAF update (#188)

* doc: update entries for `ix_interp_d`, `iy_interp_d`

- added more explanatory text
- added simple example for how to determine the input values

* doc: updates concerning the merge of the FOR2131-branch

* doc: setting of environment variable FOR2131

* doc: explanation of computation of observation interpolation input

doc/content/build_tsmp/build_environment_variables.md

@@ -53,3 +53,22 @@ This will remove unnecessary code during observation reading, when
 ParFlow-CLM-PDAF is built, but no ParFlow-type observations are
 included.
 
+## FOR2131 (PDAF / ParFlow build, experimental) ##
+
+The environment variable `FOR2131` can affect
+
+1. TSMP-PDAF (`enkf_parflow.c`), then it has to be set in the
+   PDAF-build script
+   `bldsva/intf_DA/pdaf1_1/arch/<machine>/build_interface_pdaf1_1_<machine>.ksh`
+   as mentioned above
+2. Patched ParFlow files, then it has to be set in the ParFlow build
+   script
+   `bldsva/intf_oas3/parflow/arch/<machine>/build_interface_parflow_<machine>.ksh`
+
+In `enkf_parflow.c`, there are two main effects:
+- behavior of `PF:updateflag == 2` is changed to include pressure in
+  the state vector
+- a check for saturations greater than 1 is included
+
+In `ParFlow` more outputs routines are defined and there is an
+additional saturation update.

doc/content/setup_tsmp/input_enkfpf.md

@@ -148,11 +148,26 @@ pressure data (updateflag=1) in ParFlow.
 
 -   2: Update of Mannings coefficient
 
+Experimental:
+
+- 3: Update porosity
+
+- 4: Update van Genuchten parameters
+
+- 5: Update hydraulic conductivity and porosity
+
+- 6: Update hydraulic conductivity and van Genuchten parameters
+
+- 7: Update porosity and van Genuchten parameters
+
+- 8 Update hydraulic conductivity, porosity and van Genuchten parameters
+
+
 ### PF:paramupdate_frequency ###
 `PF:paramupdate_frequency`: (integer) Frequency of parameter
 updates. Default: `1`
 
-For each assimilation cycle it is checked whether `tcylce mod
+For each assimilation cycle it is checked whether `tcycle mod
 paramupdate_frequency == 0`. When this happens, a parameter update is
 applied. For the default of `1`, each assimilation cycle contains a
 parameter update.
@@ -201,15 +216,19 @@ case of assimilated soil moisture data `PF:updateflag = 2` or
 
 `PF:paramprintensemble`: (integer) Only used in case of parameter
 update. If set to `1`, the updated parameters are printed to `pfb`
-files (similar to `[PF]printensemble`). Output files include the
-specifier `update.param`.
+files (similar to `[PF]printensemble`). Output files include a
+specifier depending on `PF:paramupdate`, f.e. `update.param.ksat`,
+`update.param.ksat`.
 
 ### PF:paramprintstat ###
 
 `PF:paramprintstat`: (integer) Only used in case of parameter
 update. If set to `1` statistics on the updated parameters are printed
 to `pfb` files (similar to `[PF]printstat`).
 
+The output files include a specifier depending on `PF:paramupdate`,
+f.e.  `param.ksat`, `param.mannings` or `param.poro`.
+
 ### PF:olfmasking ###
 
 `PF:olfmasking`: (integer) Only used in case you do not want to update

doc/content/setup_tsmp/input_obs.md

@@ -84,33 +84,118 @@ observation in x-, y- and z-direction.
 
 #### ix_interp_d ####
 
-`ix_interp_d`: (real) Distance from observation location to the grid
-position assigned to the location `ix`. Only used when
-[`DA:obs_interp_switch`](./input_enkfpf.md#daobs_interp_switch) is set to
-`1`.
+`ix_interp_d`: (real) Offset of the correct observation locations
+compared to the position of the observation in the ParFlow grid
+assigned with `ix`. Only used when
+[`DA:obs_interp_switch`](./input_enkfpf.md#daobs_interp_switch) is set
+to `1`.
 
-The observation location should be at
+The correct observation location should be located at
 
 \begin{gather*} x_{obs} =
 x_{grid}(ix) + \mathtt{ix\_interp\_d} \cdot \Delta x
 \end{gather*}
 
-and `0 <= ix_interp_d < 1`.
+where 
+- Important: `0 <= ix_interp_d < 1`.
+- $x_{obs}$: x-coordinate of the observation location
+- $x_{grid}(ix)$: x-coordinate of grid cell with index `ix` in the
+ParFlow grid. Note that $x_{grid}(ix) = \max_{ix}\{x(ix) | x(ix) <
+x_{obs}\}$ (closest grid point that is smaller than $x_{obs}$)
+- $\Delta x$: Difference of x-coordinate between grid cells
+  neighboring observation location: $\Delta x = x_{grid}(ix+1) -
+  x_{grid}(ix)$
+
+##### Example #####
+
+This example shows how to set `ix_interp_d` from a given (1) grid and
+(2) observation location.
+
+Let the observation location be at longitude $x_{obs} = 6.404$.
+
+Let the grid be simple going from `6.3` to `6.5` in steps of `0.02`.
+Example grid locations:
+-  $x_{grid}(1) = 6.3$
+-  $x_{grid}(6) = 6.4$
+-  $x_{grid}(11) = 6.5$
+
+We find `ix` and `ix_interp_d` as follows:
+
+1. Find the two grid locations with longitudes closest to the
+   observation location: $x_{grid}(6) = 6.4$, $x_{grid}(7) = 6.42$
+2. Choose the index of the grid location smaller than $x_{obs}$ for
+   observation file input `ix`: `ix=6`
+3. Compute `ix_interp_d` from the two closest grid locations as follows
+
+\begin{gather*}
+\mathtt{ix\_interp\_d} = \frac{ x_{obs} -  x_{grid}(ix)}{ \Delta x}
+\end{gather*}
+
+Here this leads to `ix_interp_d=0.2`.
+
+We can check this by computing $x_{obs}$ from `ix`, `ix_interp_d`:
+
+\begin{align*}
+x_{obs} &= x_{grid}(6) + \mathtt{ix\_interp\_d} \cdot \Delta x \\
+   &= 6.4 + 0.2 \cdot 0.02 = 6.4 + 0.004 = 6.404
+\end{align*}
 
 #### iy_interp_d ####
 
-`iy_interp_d`: (real) Distance from observation location to the grid
-position assigned to the location `iy`. Only used when
-[`DA:obs_interp_switch`](./input_enkfpf.md#daobs_interp_switch) is set to
-`1`.
+`iy_interp_d`: (real) Offset of the correct observation locations
+compared to the position of the observation in the ParFlow grid
+assigned with `iy`. Only used when
+[`DA:obs_interp_switch`](./input_enkfpf.md#daobs_interp_switch) is set
+to `1`.
 
-The observation location should be at
+The correct observation location should be located at
+
+\begin{gather*} y_{obs} =
+y_{grid}(iy) + \mathtt{iy\_interp\_d} \cdot \Delta y
+\end{gather*}
+
+where 
+- Important: `0 <= iy_interp_d < 1`.
+- $y_{obs}$: y-coordinate of the observation location
+- $y_{grid}(iy)$: y-coordinate of grid cell with index `iy` in the
+ParFlow grid. Note that $y_{grid}(iy) = \max_{iy}\{y(iy) | y(iy) <
+y_{obs}\}$ (closest grid point that is smaller than $y_{obs}$)
+- $\Delta y$: Difference of y-coordinate grid cells neighboring
+  observation location: $\Delta y = y_{grid}(iy+1) - y_{grid}(iy)$
+
+##### Example #####
+
+This example shows how to set `iy_interp_d` from a given (1) grid and
+(2) observation location.
+
+Let the observation location be at longitude $y_{obs} = 50.910$.
+
+Let the grid be simple going from `50.8` to `51.0` in steps of `0.02`.
+Example grid locations:
+-  $y_{grid}(1) = 50.8$
+-  $y_{grid}(6) = 50.9$
+-  $y_{grid}(11) = 51.0$
+
+We find `iy` and `iy_interp_d` as follows:
+
+1. Find the two grid locations with longitudes closest to the
+   observation location: $y_{grid}(6) = 50.9$, $y_{grid}(7) = 50.92$
+2. Choose the index of the grid location smaller than $y_{obs}$ for
+   observation file input `iy`: `iy=6`
+3. Compute `iy_interp_d` from the two closest grid locations as follows
 
 \begin{gather*}
-y_{obs} = y_{grid}(iy) + \mathtt{iy\_interp\_d} \cdot \Delta y
+\mathtt{iy\_interp\_d} = \frac{ y_{obs} -  y_{grid}(iy)}{ \Delta y }
 \end{gather*}
 
-and `0 <= iy_interp_d < 1`.
+Here this leads to `iy_interp_d=0.5`.
+
+We can check this by computing $y_{obs}$ from `iy`, `iy_interp_d`:
+
+\begin{align*}
+y_{obs} &= y_{grid}(6) + \mathtt{iy\_interp\_d} \cdot \Delta y \\
+   &= 50.9 + 0.5 \cdot 0.02 = 50.9 + 0.01 = 50.910
+\end{align*}
 
 ### CLM observation file variables ###