Update ropp_1dvar_add_refrac_error to handle missing data better

[Ian Culverwell]

Joe Nielsen (DMI) found some problems with ropp_1dvar_add_refrac_error when using missing data:

My impression is that the obs%weights variable is supposed to handle 
this, but if the missing values are passed to 
ropp_1dvar_add_refrac_error it produces a 247x247 correlation file 
which is not positive definite.  

(He's referring to the reduction of levels to accommodate missing data, as is/was done in 1DV.)

By using the following patch, he managed to get consistent output between ropp_1dvar and DMI 1DV:

===================================================================
--- ropp_1dvar/errors/ropp_1dvar_add_refrac_error.f90   (revision 5460)
+++ ropp_1dvar/errors/ropp_1dvar_add_refrac_error.f90   (working copy)
@@ -82,6 +82,7 @@

  ! For TP formula
    REAL(wp), DIMENSION(:), ALLOCATABLE   :: grad_T
+  LOGICAL,  DIMENSION(:), ALLOCATABLE   :: weight
    REAL(wp)                              :: trop

+ ALLOCATE(weight(ro_data%Lev2a%Npoints))

+
+          WHERE( (ro_data%lev2a%refrac      < ropp_MDTV) .OR. &
+               (ro_data%lev2a%geop_refrac < ropp_MDTV) )
+             ro_data%Lev2a%refrac_sigma = ropp_MDFV
+             weight=.FALSE.
+          ELSEWHERE
+             weight=.TRUE.
+          END WHERE
+
+

+
         DO j=1,ro_data%Lev2a%Npoints
           DO k=1, ro_data%Lev2a%Npoints
-           Ocorr(j,k) = EXP(-ABS(RO_data%Lev2a%geop_refrac(j) -    &
-              RO_data%Lev2a%geop_refrac(k))/3000.0_wp)
+
+            IF ((weight(j).AND.weight(k)).OR.(j.EQ.k)) THEN
+               Ocorr(j,k) = &
+                    EXP(-ABS(RO_data%Lev2a%geop_refrac(j) - &
+                    RO_data%Lev2a%geop_refrac(k))/3000.0_wp)
+            ELSE
+               Ocorr(j,k)= 0.0_wp
+            ENDIF
           ENDDO
         ENDDO

===================================================================
--- ropp_fm/refrac_1d/ropp_fm_refrac_1d_grad.f90        (revision 5460)
+++ ropp_fm/refrac_1d/ropp_fm_refrac_1d_grad.f90        (working copy)

-     K(:,i) = grad
+     K(:,i) = grad*y%weights

We should consider this for inclusion in ROPP, but not urgently.

[Johannes K. Nielsen]

Just noting that I did not look at bending angle 1dvar. jkn

[Johannes K. Nielsen]

It is possible that ropp_fm/bangle_1d/ropp_fm_bangle_1d_grad.f90 could need

K(:,i) = grad

changed to

K(:,i) = grad*y%weights

jkn Replying to jkn:

Just noting that I did not look at bending angle 1dvar. jkn

[Johannes K. Nielsen]

A few additional differences between DMI 1DV and ropp_1dvar were noticed during follow up tests:

1) Not really incorrect or crucial but leading to errors in some plotting routines: In att files where bending angle and refractivity are missing values impact parameter and alt_refac still contains values. Should be set to ropp_MDFV by roop2ropp. Nothing has been done in DMI trunk.

2) More important: In cases where ropp_1dvar_add_refrac_error fails to find the tropopause in the background file (TP-option), ropp_1dvar_add_refrac_error falls back on a 10000m tropopause. 1DV_add_obs_error falls back on a tropopause climatology. The 1D-Var results can be quite different - up to 0.5 K in some cases. Unfortunately the 1DV implementation of this is not coded nicely so it would require a little effort to implement the climatology solution in ROPP. We are considering whether it is worse the effort.

3) A minor difference between 1DV and ropp_1dvar appears in some cases due to uneven implementation of new_ref_op: In 1DV new_ref_op is not activated in the initial quality check (1DV_refrac), meaning that the OmB calculation used at this point causes the obs%weights to be slightly different in some cases. This is a minor error in 1DV 3.3 (just reported here in case the question should come up again).

jkn

[Ian Culverwell]

In response to the above:

1) I don't agree that the independent variables (impact parameter and refractivity altitude) should be nullified when the dependent variables (bending angle and refractivity) are missing. IP and H are valid things, even if they don't have valid quantities attached to them. And, for example, you need a complete set of IPs if you want to interpolate different bending angle profiles to the same set of levels before taking the difference. Plotting difficulties should be handled in plotting routines. (I suspect as many users would be fazed by missing IPs as would be pleased by them.)

2) Using some sort of tropopause climatology definitely sounds more sensible than defaulting to a flat constant TPH, especially if that climatology can be expressed as a simple formula (function of latitude and/or season?). If it needs an auxiliary data set to hold it, then things could get a bit messy.

We are open to such a development in ROPP.

3) Noted.

[Ian Culverwell]

Some preliminary tidying up has been done at r5609.

[Ian Culverwell]

Joe's changes to ropp_1dvar_add_refrac_error.f90 were made at r5462. But I don't understand how this can work in general, because weight is only calculated if the TP option is activated (sec 6.5), but weight is used in to calculate Obs_covs in Sec 7, even if you don't use TP. Asked Joe if he had ever tried it without TP - he thinks not.

Postpone this fix until ROPP10.0, as it is unlikely to be fixed anytime soon.

[Johannes K. Nielsen]

I never claimed that it was supposed to work in general. However I will fix this and see that it is merged it into DMI trunk, and the eventually it will go into ROPP that way. joe

[Johannes K. Nielsen]

1D-Var whith some data cut out in te lower stratosphere

[Johannes K. Nielsen]

1D-Var with no data cut out.

[Johannes K. Nielsen]

So the weight assignment has been moved down (outside the CASE structure), so it will now cover all choices of refractivity error models.

Changed the name of logical "weight" to "quality_mask"

Moreover I have decided to set _data%Lev2a%refrac_sigma equal to unity in cases where quality_mask == zero. This is because in the inversion of the observation error covariance diagonal elements with ropp_MDTV^{2 would possible cause a "non positive definite error" from time to time. 1.0_wp is chosen because it is in the realistic range of refrac_sigma. The philosophy here is that the O matrix is composed of blocks with valid covariances and blocks of pure (non-valid) diagonal contributions. See attached "missing data" figure. These blocks do not mix during inversion. In the minimizers, the non valid ranges are completely ignored by a weight function. The problem with positive indefiniteness has not been observed in practise though.}

I will check it into ropp_src/branches/dev/Share/jkn_dmi_trunk_9.0/ropp_1dvar/errors/ropp_1dvar_add_refrac_error.f90 for now. I believe the problem is basically solved for refractivity 1D-Var but not for bending angle 1D-Var joe

[Johannes K. Nielsen]

Made sure that in ropp_fm/bangle_1d/ropp_fm_bangle_1d_grad.f90

K(:,i) = grad

changed to

K(:,i) = grad*y%weights

ropp_1dvar_add_bangle_error does not have the feature of exporting correlation matrices. It is up to the user to provide a consistent correlation matrix if desired. Therefore there is no needed fix for ropp_1dvar_add_bangle_error. One might fear that there is a risk of mishaps in the unlikely event that a user should decide to run ropp_1dvar_bangle with correlations without paying attention to missing values.

For now I consider this issue closed, but I leave the ticket open until we have migrated the changes into the dmi_trunk. joe

[Stig Syndergaard]

Moving ticket to DMI development. Changing milestone and owner.