Merge Joe's fix to 1dvar with q<0

[Stig Syndergaard]

Fixed here: https://trac.romsaf.org/ropp/changeset/6425/ropp_src/branches/dev/Share/jkn_dmi_trunk_9.0.3

But corresponding fix in adjoint is missing.

[Stig Syndergaard]

As is, the fix is only correct for the choice of logq = false. Some work is necessary so that the code still works correctly for the choice of use_logq = true.

This work should be done before merging to dmi_trunk. Therefore assigning ticket to jkn. When code is ready, the ticket can be reassigned to ssy.

[Johannes K. Nielsen]

A new option, "-let_qlt0", has been added to ropp_1dvar_refrac and ropp_fm_bg2ro, and a new control variable "let_qlt0" has been introduced in State1dFM. The philosophy is that default is not to allow q<0. "let_qlt0" can also be set in the 1DV_config_.cf provided at the command line for ropp_1dvar_refrac. Command line option "-let_qlt0" overrules configuration file setting.

When state%let_qlt0 is .TRUE. ropp_fm_refrac1d applies linear interpolation of humidity as function of pressure to observation levels, instead of the default exponential interpolation as function of geopotential height. Some notes on the nature of this interpolation will follow below.

[Johannes K. Nielsen]

Illustration of interpolation bias.

[Johannes K. Nielsen]

Regarding the -use_logq=.TRUE. situation: roprof2state checks whether q<0 and sets x%state_ok = .FALSE. in that case, thus there is not a conflict; if a bgr-file with q<0 is imported with -use_logq=.TRUE. ROPP throws an error, and no applications will attempt to run the forward model. If -use_logq=.TRUE., -let_qlt0=.TRUE. and the background file has only positive q then the ropp_1dvar_refrac will go into the minimization, and in each step the minimizer will do increments in log-space, while the forward model will map from log to lin, and refractivity will evaluated with linear q interpolation (with pressure as free variable). In that case negative q solutions can never be reached. Note that the use_logq option has nothing to do with the interpolation of q, it only concerns the representation of q in the minimizers incrementation-step. If -use_logq=.TRUE. and -let_qlt0=.FALSE., then the forward model maps the state vector-q to linear space and then back to log-space (for interpolation), and then back to linear space for refractivity evaluation. While the combination of -use_logq=.TRUE., -let_qlt0=.TRUE. is not meaningful, it will not do any harm, so we do not implement any check for that.

##########################################################

Interpolation issue!

##########################################################

The ropp_refrac_1d_new algorithm was introduced to avoid an issue with the logarithmic interpolation of refractivity at high altitudes, but it introduced a bias at low altitudes, which has not received sufficient attention. The -new_opt moves the issue with interpolating a convex exp-like function from refractivty to humidity. The problem is that the specific humidity is not really exponential at low altitude, near the boundary layer. The most characteristic feature of the specific humidity is very steep vertical gradients. Mostly negative and close to linear in shape. The default loq(q) interpolation fails to capture these, thus introducing a negative bias. This becomes evident when we compare the "-new_op" to the -let_qlt0 (linear interpolation along pressure). Pressure is chosen instead of geopotential height as independent variable, because it varies exponentially, so the specific humidity varies a little more smoothly in pressure that geopotential height, but this is of little importance. The main problem is that the logarithmic interpolation (log(q/q_0)=q_0+alpha(q_0) underestimates log(q) when the gradient is large. In figure 1. the impact on refractivity of this interpolation error is viewed for few handfuls of profiles. It is quite prominent. Figure 2. shows how the undershoot of the logarithmic interpolation of specific humidity is associated with steep gradients, and not with exponential dependency.

###################################################################

TESTS

###################################################################

Needless to say: t_fascod fails already before the fix. t_fascod_tl has been changed such that it also can be evaluated with -let_qlt0.

[Johannes K. Nielsen]

At little more about figure 2:

Upper panel, occultation 1 Left: Background specific humidity (black) interpolated in log space -new_op (blue) and in lin space (orange). Middle: difference (lin minus log interpolation). Right: Same as left plot with altitude zoomed in. Lower panel, occultation 2 Left: Background specific humidity (black) interpolated in log space -new_op (blue) and in lin space (orange). Middle: difference (lin minus log interpolation). Right: Same as left plot with altitude zoomed in.

The conclusion is in both cases that whenever there is a steep, close to linear, gradient, the log interpolation underestimates q.

More about tests: The t_fascod_ecmwf_tl.f90 _fascod_ecmwf_ad.f90 + Makefile.am that I coded never made it into the dmi_trunk, so I grabbed them from ​https://svn.grassaf.org/ropp/ropp_src/branches/dev/Share/jkn_dmi_trunk_9.0/ropp_fm/tests Here is the output of t_fascod_ecmwf_tl:

jkn@b2795a8ef76e:/dvl/jkn_dmi_trunk_9.0.3/ropp_fm/tests$ t_fascod_ecmwf_tl -new_op -let_qlt0

TL: x%non_ideal, x%new_ref_op, x%new_bangle_op = F T T

Checking REFRACTIVITY TL FASCOD_TRO |dy|/|dy_init| cos_angle 100*rel diff

1.0000000000E+00 8.0497227753E-01 7.8480548500E+01 1.0000000000E-01 8.8036029621E-01 5.1717474415E+01 1.0000000000E-02 9.9672212057E-01 8.0934967295E+00 1.0000000000E-03 9.9995714676E-01 9.5942180562E-01 1.0000000000E-04 9.9999791041E-01 2.9246038597E-01 1.0000000000E-05 9.9999838576E-01 2.7249335639E-01 1.0000000000E-06 9.9999839721E-01 2.7175029569E-01 1.0000000000E-07 9.9999839800E-01 2.7168929667E-01 1.0000000000E-08 9.9999839807E-01 2.7168358298E-01 1.0000000000E-09 9.9999839807E-01 2.7168375058E-01 1.0000000000E-10 9.9999839812E-01 2.7168112881E-01 1.0000000000E-11 9.9999839990E-01 2.7164080713E-01 1.0000000000E-12 9.9999838623E-01 2.7329884245E-01 1.0000000000E-13 9.9999811070E-01 2.9086144295E-01 1.0000000000E-14 9.9998739376E-01 5.3734120779E-01 Refractivity TL check passed

Checking BENDING ANGLE TL FASCOD_TRO |dy|/|dy_init| cos_angle 100*rel diff

1.0000000000E+00 -9.1669466336E-02 1.0000000007E+02 1.0000000000E-01 -1.1552613550E-01 1.0000000001E+02 1.0000000000E-02 -9.8666521003E-02 1.0000000000E+02 1.0000000000E-03 1.0855067677E-01 9.9999999999E+01 1.0000000000E-04 9.9267741769E-01 1.2108170279E+01 1.0000000000E-05 9.9995015563E-01 1.0711914164E+00 1.0000000000E-06 9.9999949645E-01 1.0757206473E-01 1.0000000000E-07 9.9999999497E-01 1.0753242096E-02 1.0000000000E-08 9.9999999996E-01 8.8118028242E-04 1.0000000000E-09 9.9999999980E-01 2.3095759338E-03 1.0000000000E-10 9.9999998536E-01 1.8206981892E-02 1.0000000000E-11 9.9999913002E-01 1.3358902900E-01 1.0000000000E-12 9.9983362374E-01 2.1093627856E+00 1.0000000000E-13 9.7955712546E-01 2.0232031159E+01 1.0000000000E-14 4.6949483313E-01 9.6621338789E+01 Bending angle TL check passed

PASS

[Johannes K. Nielsen]

Summary: A new switch, let_qlt0, has been introduced. The code ropp_fm_refrac_1d_new and ropp_fm_refrac_1d_new_tl code can run both in the normal mode (with log(q) interpolation) and in let_qlt0 mode (with linear q interpolation). The let_qlt0 mode has not been implemented in ropp_fm_refrac_1d_new_ad. TL test allowing let_qlt0 has been implemented in ropp_fm/tests/t_fascod_ecmwf_tl.f90. Test passes.

Comments: 1) In future opp_fm_refrac_1d_new_ad (and corresponding test) should be updated for NWP use.

2) Options should be reorganized such that ROM SAF preferences are default.

3) Tests t_fascod_ecmwf_tl.f90 and t_fascod_ecmwf_tl.f90 should be merged, but it requires a rethinking of the whole test-concept, because it has to be aligned with MetOffice model.

The following files have changed (rev. 6675 and 6682)

ropp_1dvar/tools/ropp_1dvar_refrac.f90 ropp_1dvar/common/ropp_1dvar_read_config.f90

ropp_fm/common/ropp_fm_types.f90 ropp_fm/refrac_1d/ropp_fm_refrac_1d_new.f90 ropp_fm/refrac_1d/ropp_fm_refrac_1d_new_tl.f90 ropp_fm/tools/ropp_fm_bg2ro_1d.f90 ropp_fm/tests/Makefile.am ropp_fm/tests/t_fascod_ecmwf_ad.f90 ropp_fm/tests/t_fascod_ecmwf_tl.f90

joe