pennelise / GOOPy

2 stars 1 forks source link

General Observation Operator for Python (GOOPy)

Mass-conserving vertical interpolation and satellite observation operator for inclusion in GCPy.

Discussed in this issue for GCPy: https://github.com/geoschem/gcpy/issues/242

Installing GOOPy

Using GOOPy

Configuration file

This file describes the structure of the satellite or model files used as inputs for the satellite operator.

Local Settings:

Model files:

These should not change much between GEOS-Chem files, except for CONC_AT_PRESSURE_CENTERS, which is different for CO2 and CH4.

The description of the model files requires the following fields to be defined:

Description of method for mass-conserving vertical interpolation

We follow the method described in Keppens et al. (2019), which describes a mass-conserving interpolation function in eq. 14:

$$ W' = M{out}^* W M{in} $$

Where $M{in}$ converts from concentrations at pressure centers on the model profile to partial columns on the model profile, $M{out}$ converts from partial columns on the satellite profile to concetrations at pressure edges or pressure centers on the satellite profile, and W is the interpolation matrix defined by eq. 13:

$$ W(i,j) = \frac{1}{\Delta p{in,j}}(min(p{out,i}^U, p{in,j}^U) - max(p{out,i}^L, p_{in,j}^L)) $$

Eq. 14 can be used to interpolate to both pressure edges or pressure centers depending on the choice of $M^*_{out}$.

We define $M{out}$, which transforms from pressure centers (or edges) to partial columns, and then take the inverse $M^*{out}$ to transform from partial columns to pressure centers (or edges).

To transform from N-1 pressure centers to partial columns, $M_{out}$ is just the pressure difference for each layer, and has dimension (N-1)x(N-1).

What if the averaging kernel is defined on pressure edges?

In the case of some satellite products (such as GOSAT) where the averaging kernel is defined on pressure edges, we now must interpolate this set of N-1 partial columns to N pressure edges. To avoid an underconstrained solution, we add an additional partial column by transforming to partial columns defined between layers. Instead of using the original satellite grid $h$, we instead use one that we refer to as $h'$. The transformation is described by eq. 11 (Keppens et al., 2019):

$$ h' = \begin{bmatrix} 1 & 0 & 0 & 0 \ .5 & .5 & 0 & 0 \ 0 & \ldots & \ldots & 0 \ 0 & 0 & .5 & .5 \ 0 & 0 & 0 & 1
\end{bmatrix} h \ $$

We also illustrate the transformation from $h$ to $h'$ here:

This now allows us to use $M_{out}$ with dimension NxN, which is invertible.

Assumptions: