Straighten out p_full_level interpolation

[Johannes K. Nielsen]

1) Understand which interpolation method is correct given that ECMWF changed to finite element discretization in 2002.

2) Implement and test possible missing methods in fm, tl and ad. 3) Implement optional interpolation choices

[Johannes K. Nielsen]

We need adjoint code of

https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_fm/model_ecmwf/ropp_fm_state2state_ecmwf_tl.f90

in order to get the S&B correction into the official ROPP.

Cut and paste from e-mail 1st of May 2018:

Apologize for the not so elegant code. We shall repair this before we code the adjoint. Still I expect that someone more trained is going to do the hard work.

joe

On 2018-05-01 18:00, Sean Healy wrote:

Yes, your TL is a bit tricky to follow, and doesn't really follow the standard TL rules.

As you comment,

!Should probably use p_hlv_tl and del_p_tl

I think is right!

and I'd expect to see the p_flv_tl computed in the same loop as p_flv. In the code I'd expect

p_flv(n_flv)=0.5_wp * (p_hlv(n_hlv-1)-p_hlv(n_hlv))

p_flv_tl(n_flv)= SHOULD BE HERE == 0.5_wp * (p_hlv_tl(n_hlv-1)-p_hlv_tl(n_hlv))

DO lvl=n_flv-1,1,-1

p_flv(lvl)=exp(1.0_wp/(p_hlv(lvl)-p_hlv(lvl+1)) &

*(p_hlv(lvl)*log(p_hlv(lvl))- p_hlv(lvl+1)*log(p_hlv(lvl+1))) - 1.0_wp)

!!!! TL bits SHOULD BE HERE here

p_flv_tl(lvl) = ...

ENDDO

!

On 01/05/18 15:52, Johannes K. Nielsen wrote:

Hello, Ian and Sean. Thanks for the input. I once followed the NWP course on ECMWF and I worked through the Lorenz exercise concluding that it is good to know that this can be done and I am glad that someone else is taking care of it. So a little help is probably needed when we dig into it.

yes we have the tangent linear

https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_fm/model_ecmwf/ropp_fm_state2state_ecmwf_tl.f90

First task is to clean it up.

As far as I recall there was a problem with the ROPP verification tool, so this was verified with a few perturbation experiments.

Second task is to write the adjoint in ROPP style. I am glad that we apparently have someone else here at DMI to do that.

joe

On 2018-05-01 16:33, Culverwell, Ian wrote:

Hi Kent,

Thanks for this. I'm happy to help Joe, but I'm no expert on adjoint coding. All I know is that they’re bloody difficult to get right - and you do have to get them exactly right. The adjoint complements the tangent linear, so you have to match whatever equations are in the TL code. As far as I can see, these are the full S&B 1981 eqns.

Ian.

---

From: Kent Bækgaard Lauritsen kbl@… Sent: 01 May 2018 15:22 To: Culverwell, Ian; sean.healy@… Cc: Johannes K. Nielsen; Hans Gleisner; Stig Syndergaard Subject: SV: ROPP 9.1 changes

Hi Ian

I understand this point. We have discussed this at DMI and we will involve our IT department and try to code this missing AD routine. Joe believes he knows what needs to be done but if there are questions it would be nice if he can discuss them with you.

We (you and I) briefly discussed this when you were here at DMI in January. And - if I recall correctly - I think we discussed that rather than using the full expression from Simmons-Burridge one could perhaps Taylor expand to second order and code the adjoint for the latter. Is my recollection of this discussion correct and would you recommend that approach or is it better to use the full equation?

Thanks Kent

---

Fra: Culverwell, Ian ian.culverwell@… Sendt: 1. maj 2018 12:26 Til: Kent Bækgaard Lauritsen; sean.healy@… Cc: Johannes K. Nielsen; Hans Gleisner; Stig Syndergaard; Culverwell, Ian Emne: RE: ROPP 9.1 changes

Hi Kent,

I'm sorry but the ML change cannot go into ROPP until the adjoint code has been written. As currently written, the changes affect the default behaviour, and will therefore cause the model to fail the adjoint tests (and, more importantly, for 1DVAR to give the wrong solution) if a general user were to run it with minROPP.

[Johannes K. Nielsen]

The format messed up so here is seans suggestion again:

As you comment,

!Should probably use p_hlv_tl and del_p_tl

I think is right!

and I'd expect to see the p_flv_tl computed in the same loop as p_flv. In the code I'd expect

p_flv(n_flv)=0.5_wp * (p_hlv(n_hlv-1)-p_hlv(n_hlv))

p_flv_tl(n_flv)= SHOULD BE HERE == 0.5_wp * (p_hlv_tl(n_hlv-1)-p_hlv_tl(n_hlv))

DO lvl=n_flv-1,1,-1

p_flv(lvl)=exp(1.0_wp/(p_hlv(lvl)-p_hlv(lvl+1)) &

*(p_hlv(lvl)*log(p_hlv(lvl))- p_hlv(lvl+1)*log(p_hlv(lvl+1))) - 1.0_wp)

!!!! TL bits SHOULD BE HERE here

p_flv_tl(lvl) = ...

ENDDO

!

[Johannes K. Nielsen]

Regarding the ad coding. I am rewriting the t_fascod_ad script to also test the state2state part, which has actually never been tested. The problem has been that the t_fascod* scripts starts with a perturbation of the forward modeled state (x_tl%tem, x_tl%geop etc.), instead of a perturbed state x_tl%state. I.e. the state2state-part has been bypassed. I am writing a new test t_fascod_ecmwf_tl and t_fascod_ecmwf_ad to test state2state_ecmwf properly. I am not doing this for the MetOffice model.

The first problem seems to be that extended t_fascod_ecmwf_tl cannot detect the numerical difference between the old p_flv calculation and the Simmons and Burridge calculation!

The t_fascod_ecmwf_ad does not work immediately, even for the old (pre-Simmons and Burridge) ropp_state2state_ecmwf_ad code.

[Johannes K. Nielsen]

Exchanged 0.61 with ((R_vap/R_dry)-1.d0) in ropp_fm_state2state_ecmwf_ad.f90 - and ropp_fm_state2state_ecmwf.f90

[Johannes K. Nielsen]

The old simple full level calculation is now consistent. Both TL and AD. I had to relax the AD tolerance parameter from tol=1.0e-9_wp to tol=1.0e-7_wp.

Also, I believe t_fascod_ad does not do the right thing, even with these modifications. What is being used as test vectors is two idendical vectors. That is a bit too specific:

norm1 = DOT_PRODUCT(y_tl, y_tl)

The adjoint code has to be able to fulfill <K^{T*dy,dx>=<dy,K*dx> for any pair (dy,dx). Thus, I have rewritten the code in t_fascod_ecmwf_ad to test a randomly chosen pair of perturbation vectors instead. Further more I have reduced the scope to check only ropp_fm_state2state_ecmwf_ad.f90. The test of ropp_fm_refrac_1d_new_ad and ropp_fm_bangle_1d_new_ad, starting from [temp,shum,pres,geop], can in principle be tested with t_fascod_ad (if this code is cleaned up and debugged).}

Next step is to re-implement Simmons and Burridge and run the t_fascod_ecmwf on that.

[Johannes K. Nielsen]

Noting t_fascod_ecmwf_tl cannot detect difference between LIN and Simmon and Burridge. With this in mind we can consider the problem solved; Following the logic that the ad code should be traceable line line by line when reversing the tl code I decided to move section 6.6 to right after section 6.8. In conclusion both _fascod_ecmwf_tl and _fascod_ecmwf_ad is now passed with the new Simmons an Burridge code. Keeping in mind that the t_fascod_ecmwf_tl cannot not see the difference between the new and the old code. (t_fascod_ecmwf_ad fails at the old code, as it should.) Checking into jkn branch

[Johannes K. Nielsen]

Documentation and test output.

[Johannes K. Nielsen]

We discovered some numerical instability in the 1D-Var output of the new tl/ad, which has been resolved by revision ​https://trac.grassaf.org/ropp/changeset/6065/ropp_src/branches/dev/Share/jkn_dmi_trunk_9.0.2/ropp_fm/model_ecmwf and ​https://trac.grassaf.org/ropp/changeset/6067/ropp_src/branches/dev/Share/jkn_dmi_trunk_9.0.2/ropp_fm/model_ecmwf6067. At the same time a the TL code has been changed such that x%bk is used instead of p_hlv_tl, because the old implementation did multiply with p_sfc_tl twice. The threshold for AD test has been lowered to 10e-8. A documentation and test document has been attached to this ticket.

[Johannes K. Nielsen]

Oh and Stig found yet another bug, corrected here:

Changeset 6111

Timestamp:

12/19/19 12:49:47

Author:

jkn

Message:

Stig found an error in the adjoint code,

p_hlv_ad(1:n_hlv-1) = p_hlv_ad(1:n_hlv-1) - alpha_ad/del_p*ln_prflv --> p_hlv_ad(2:n_hlv) = p_hlv_ad(2:n_hlv) - alpha_ad/del_p*ln_prflv

-much better test result

[Stig Syndergaard]

A series of revisions has been made to ropp_fm_state2state_ecmwf.f90 and associated AD/TL codes. They are r6112, r6113, r6114, r6115 (see also #661), r6116, r6117, r6120, r6121, r6122, r6123, r6124, and r6130. The last of these can be considered the final solution, and tests have shown that TL and AD codes are now fully consistent (see #661).

The solution in Section 7 of https://trac.romsaf.org/ropp/changeset/6130/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_fm/model_ecmwf/ropp_fm_state2state_ecmwf.f90 is the same as what was used to generate the ROM SAF CDR and ICDR, just written in a different way using already defined variables. The changes has been included in GPAC 2.4.0.1. Verification using the Vagrant environment showed that these changes only affect results in the 9th significant digit or beyond (in all relevant variables in all GPAC file types) compared to what was implemented for GPAC 2.4.0.0.