MODULE presProfileOperators_mod (prefix=’ppo’ category=’8. Low-level utilities and constants’)
Purpose
Vertical interpolation, integration, and layer averaging subroutines.
Includes the special routines designed to interpolate to the
(widely spaced) RTTOV pressure levels.
To perform the vertical interpolation in log of pressure and
constant-value extrapolation of one-dimensional vectors. Input
pressure levels can vary for each profile.
Determination of interpolation weights for interpolation to points
in a profile. Applies interpolation in log(Pressure).
Input
pressInput
pressure on reference column levels assumed to be in ascending order
pressTarget
target pressure levels assumed to be in ascending order
numInputLevs
number of input/reference levels
numTargetLevs
number of target levels
Options
method_opt [character ,in,] ::
Specified interpolation method
:skipType_opt: Skipping processing of specific target levels depending on case:
’default’ - extrapolation allowed and skipping application via input success_opt only
‘noExtrap’ - no extrapolation as well as additional skipping via input success_opt
‘doAll&noExtrap’ - no extrapolation only (all other levels processed)
skiptype_opt [character ,in,] :: Skipping processing of specific target levels depending on case
Inout
Arguments
outbound_opt (numtargetlevs) [integer ,inout,] :: Flag indicating if obs outside model vertical range (0 for no)
success_opt (numtargetlevs) [logical ,inout,] :: success of interpolation
Determine profile interpolation weights by considering integrations
over of series of segments [PX1(KI-1),PX1(KI+1)] using piecewise
linear weighting with weights of zero at KI-1 and KI+1 and
max weight at KI. KI ranges from 1 to KN1.
Can also provide gradient contributions from both
linear and non-linear components of interpolator. The non-linear
components (case PZS* is present) stem from vertical coordinate
independent variables (e.g. dependency on Ps). The gradients
contributions from the linear components are the interpolation weights.
array of gradients from linearized component.
(dF/dv) or (dF/dv)*(v-vo) provided when pzs* is present
Method:
Piecewise weighted interpolation in ln(P).
Journal reference:
Rochon, Y.J., L. Garand, D.S. Turner, and S. Polavarapu.
Jacobian mapping between coordinate systems in data assimilation,
Q. J. of the Royal Met. Soc., vol 133, 1547-1558, 2007.
(www.interscience.wiley.com) DOI:10.1002/qj.117
The input profile is now extrapolated at constant value.
The impact is of most practical significance for instruments where
the weighting function peaks at or near the surface such as SSMI.
This approach increases the weights of contribution from
the lowest and highest input domain levels for output
layers intersecting these input domain boundaries. It consists
of applying contant value extrapolation by introducing
‘fake’ or ‘virtual’ layers. For the lowest level of the input domain,
as example, this implies creating a virtual surface layer which
extends to the lower boundary of the output domain layer which
contains the input domain surface. This increases the
contributing weight of the surface which would be otherwise
understimated in the original code due to the interpolator
actually doing piecewise weighted averaging.
COmment out use of ‘zb’ for consistency with RTTOV-9 when
Options
rttov_opt [logical ,in,] ::
Commented out use of ‘zb’ when .true. for consistency with RTTOV-9
See the four lines ending with !C1 and version 7 comment above.
3) A major reduction in computational time results from only
assigning values to the non-zero ranges of the 2D output
arrays. These ranges of the 2D areas are identified by ‘kstart’
and ‘kend’. Initialization to zero for values within these ranges
is done using 1D work arrays ‘zpz’ and ‘zpzd’, with the resulting
values then being assigned to the related elements of ‘PZ’ and
‘PZDPS’.
Therefore, elements of ‘PZ’ and ‘PZDPS’ outside these ranges
could be undefined (i.e. NaN if not 0.0) and should not be used.
Other notable reductions in computational time stem (a) from inlining
of ‘sublayer’ code (applied to a reduced degree in ‘layeravg’ as
compared to ‘rttov_layeravg_*) and (b) from updating the start
position ‘istart’ of the loops over J. The latter was faciliated
by moving the loop over KI inside ‘layeravg’.
These improvements were originally devised and implemented by
Deborah Salmond and Mats Hamrud (ECMWF) in the RTTOV-9 routines
‘rttov_layeravg*’.
4) Contributors to improvements and changes to the original version
of ‘layeravg’ and to ‘rttov_layeravg*’: members of the RTTOV9
development team, namely Niels Bormann, Alan Geer, Deborah Salmond,
and Mats Hamrud of ECMWF, Peter Rayer and Roger Saunders of the
Met Office and Pascal Brunel of Meteo-France, and Y.J. Rochon (EC).
Arguments
px1 (kn1) [real ,in] :: Levels of output domain (e.g. lnP; in increasing values)
px2 (kn2) [real ,in] :: Levels of input domain (e.g. lnP; in increasing values)
kn1 [integer ,in,] :: Dimension of PX1
kn2 [integer ,in,] :: Dimension of other arrays
pz (kn1,kn2) [real ,out] :: Resultant accumulated weighting factors for interp from input to output domain
kstart (kn1) [integer ,out] :: Start index for relevant PZ row
kend (kn1) [integer ,out] :: End index of relevant PZ row
pzs1_opt (kn1) [real ,in,] :: dPX1/dv or perturbation dPX1/dv * delta(v) where PX1(v), e.g. v=Ps
pzs2_opt (kn2) [real ,in,] :: dPX2/dv or perturbation dPX2/dv * delta(v) where PX2(v), e.g. v=Ps
pzdps_opt (kn1,kn2) [real ,out,] :: dF/dv or perturbations (dF/dv)*(v-vo): Resultant accumulated factors for the gradients w.r.t. v (or perturbations) associated to coordinates
validlevel_opt (kn1) [logical ,in,] :: Logical indicating validity of each output level
Determine weight coefficient contributions to w1 and w2 to assign
to input domain (e.g. NWP model) variables at x1 and x2. Weights
are determined from integration over the intersecting segment (y1,y2)
of the ranges (x1,x2) for the input domain and (z1,z2) for the
output domain. Integrals are approximated via the trapezoidal rule:
This is synonomous to having an integrand linear in x.
In the above (and below) equation(s), w1 and w2 are contributions to
the input values.
The weights for linearized contributions of non-linear interpolator
components, i.e. gradient w.r.t. the vertical coordinate
independent variable (e.g. v*=Ps), are calculated
when LGRADP* is .true.:
This routine provides terms on the right-hand-side.
Note: pxs* and pzs* can be provided either as gradients or
perturbations.
Method:
- Piecewise weighted interpolation in ln(P).
- Journal reference:
Rochon, Y.J., L. Garand, D.S. Turner, and S. Polavarapu.
Jacobian mapping between coordinate systems in data assimilation,
Q. J. of the Royal Met. Soc., vol 133, 1547-1558, 2007.
(www.interscience.wiley.com) DOI:10.1002/qj.117
z2 [real ,in] :: Lower level of output domain half-layer
dzd [real ,in] :: 1.0/(z2-z1)=1.0/dz
wgt2 [real ,in] :: Weight at z2 (1.0 or 0.0 or …)
x1 [real ,in] :: Upper boundary of input layer (x1<x2)
x2 [real ,in] :: Lower boundary of input layer
dxd [real ,in] :: 1.0/(x2-x1)=1.0/dx
w1 [real ,inout] :: Starting (in) and updated (out) weight assigned to variable at upper level x1
w2 [real ,inout] :: Starting (in) and updated (out) weight assigned to variable at upper level x2
pzs1 [real ,in] :: dz1/dvz or perturbation dz1/dvz * delta(vz)
pzs2 [real ,in] :: dz2/dvz or perturbation dz2/dvz * delta(vz)
pxs1 [real ,in] :: dx1/dvx or perturbation dx1/dvx * delta(vx) (required when LGRADPx=.true.)
pxs2 [real ,in] :: dx2/dvx or perturbation dx2/dvx * delta(vx)
wps1 [real ,inout] :: Starting (in) and updated (out) value of (pxs1*dw1/dx1 + pxs2*dw1/dx2 +pzs1*dw1/dz1 + pzs2*dw1/dz2)
wps2 [real ,inout] :: Starting (in) and updated (out) value of pxs1*dw2/dx1 + pxs2*dw2/dx2 +pzs1*dw2/dz1 + pzs2*dw2/dz2) (required when LGRADP*=.true.)
lgradpx [logical ,in] :: Input domain logical indicating if gradient w.r.t. vertical coordinate independent variable if required i.e. d/dv where P(v) (e.g. v=Ps).
lgradpz [logical ,in] :: Output domain logical indicating if gradient w.r.t. vertical coordinate independent variable if required i.e. d/dv where P(v) (e.g. v=Ps).
Preliminary calculations for producing components required for
vertical integration (or averaging) w.r.t. pressure for the full
vertical rangeor a set of target layers. To be called before
routine ppo_vertIntegWgts or ppo_vertAvgWgts.
Integration calculations are performed appling quadratic interpolation
in P between level.
Input
operatorType
‘integ’ for intergration; ‘avg’ for averaging
pressInput
Reference input levels
numInputLevs
# of model vertical levels
Output
boundaries(numInputLevs+1)
Input layer boundaries assuming
provided input levels can be taken as
mid-layer values.
weights
Second order Lagrangian
interp integration weights or
unity for averaging weights
Comments
This subroutine does the following:
Setting of layer boundaries
If integration, determining integration weights associated
to second order Lagrangian interpolation. Otherwise, initialize
weights to unity.
Layer boundaries are taken as mid-point between provided levels in
lnP coordinate. Layer values are set to be the values
interpolated to the mid-point in P within the various layers.
Interpolation in P is done quadratically.
Arguments
operatortype [character ,in] :: ‘integ’ for integration; ‘avg’ for averaging
To calculate integration weights “wgts” required for vertical integration w.r.t.
pressure for the full vertical range or a set of target layers.
Given the calculated weights and a user intergrand array vector X, the integral
for a given layer i would be given by sum(wgts(i,:)*X(:))
Integration calculations are performed applying quadratic interpolation
in P between level.
Routine ppo_vertLayersSetup to be called beforehand to generate Lagrangian weights
and related layer boundaries (arrays ‘weights’ and ‘boundaries’)
Input
targetLayersTop
top of target layers
targetLayersBot
bottom of target layers
numInputLevs
# of original/input vertical levels
numTargetLevs
# of target vertical levels
Arguments
kstart (numtargetlevs) [integer ,inout] :: Index of first relevant original/input level for each target level
kend (numtargetlevs) [integer ,inout] ::
Index of last relevant original/input level for each target level
If kstart and kend are non-zero on input,
the input are initial estimates of the values.
weights
See routine ppo_vertLayersSetup
boundaries
Boundaries of input layers
skipType_opt
Skipping processing of specific target layers depending on case:
‘default’ - skipping application via input success_opt only
‘doAll&noExtrap’ - application of both success_opt and outbound_opt
outbound_opt (numtargetlevs) [integer ,inout,] :: Flag indicating if obs outside input vertical range (0 for no)
success_opt (numtargetlevs) [logical ,inout,] :: success of interpolation
To calculate averaging weights “wgts” required for vertical averaging
w.r.t. ln(pressure) for the full vertical range or a set of target layers.
Given the calculated weights and a user input array vector X, the average
for a given layer i would be given by sum(wgts(i,:)*X(:))
Routine ppo_vertLayersSetup to be called beforehand to initial weigths
and related layer boundaries (arrays ‘weights’ and ‘boundaries’)
Input
targetLayersTop
top of target layers
targetLayersBot
bottom of target layers
numInputLevs
# of original/input vertical levels
numTargetLevs
# of target vertical levels
Arguments
kstart (numtargetlevs) [integer ,inout] :: Index of first relevant original/input level for each target level
kend (numtargetlevs) [integer ,inout] ::
Index of last relevant original/input level for each target level
If kstart and kend are non-zero on input,
the input are initial estimates of the values.
weights
See routine ppo_vertLayersSetup
boundaries
Boundaries of input layers
skipType_opt
Skipping processing of specific target layers depending on case:
‘default’ - skipping application via input success_opt only
‘doAll&noExtrap’ - application of both success_opt and outbound_opt
outbound_opt (numtargetlevs) [integer ,inout,] :: Flag indicating if obs outside input vertical range (0 for no)
success_opt (numtargetlevs) [logical ,inout,] :: success of interpolation
