ropp_fm_bg2ro_1d -new_op can generate NaNs above model top

[Ian Culverwell]

Hans Gleisner (DMI) reports that the new interpolation scheme can generate NaNs above the model top:

          i         geop(i)                  refrac(i)
         773   77200.0000000000       4.595168840727702E-003
         774   77300.0000000000       4.513892752756471E-003
         775   77400.0000000000       4.433896771776462E-003
         776   77500.0000000000      NaN                    
         777   77600.0000000000      NaN                    
         778   77700.0000000000      NaN                    

This is not ideal.

I have traced this to the section of ropp_fm/refrac_1d/ropp_fm_refrac_1d_new.f90 which says

! 3.4 Find model layer (j to j+1) in which ob lies
! ------------------------------------------------
    DO j=1, x%n_lev-1
      IF (y%geop(i) < x%geop(1)) THEN
        below_lowest_lev=.TRUE.
        EXIT
      END IF
      IF (y%geop(i) > x%geop(j) .AND. y%geop(i) < x%geop(j+1)) THEN
        EXIT
      END IF
    END DO

The problem is that if y%geop(i) > x%geop(x%n_lev), ie, the ob is above the model top, then the last IF ... EXIT is never executed and j is updated to the 'next' index in the loop, ie j becomes x%n_lev. All later references to variables at index j+1 are therefore outside the array bounds. This causes, in Hans' case, alpha (=-dlog q/dz) to be minus infinity, and generates all the NaNs.

I see two possible solutions.

(1) Extrapolate the top layer variables, between j=n_lev-1 and n_lev above the model top. In other words, set j=x%n_lev-1 if the ob is above the model top. This is fine unless one is extrapolating wildly, as Hans was (up to 750 km in his case). You can then get negative temperature t_ob etc. This causes NaNs for a different reason.

(2) Safer would be to introduce a new logical above_highest_lev, analagous to below_lowest_lev, which would set refracs above the model top to ropp_MDFV. This fixes Hans's problem, but is it what we want to do?

(Attached is Hans' test file, which he ran with ropp_fm_bg2ro_1d singlefile.bg.nc -f -comp -new_op --no-ranchk -d -o singlefile.bg_out.nc.)

[Ian Culverwell]

singlefile.bg.nc

[Ian Culverwell]

For the record, the above_highest_levels switch produces

refrac(773)=0.00459516884073 
refrac(774)=0.00451389275276 
refrac(775)=0.00443389677178 
refrac(776)=-99999000 
refrac(777)=-99999000 
refrac(778)=-99999000 

The bending angles are unchanged, of course: they're already ropp_MDFV above the model top. (And, I suggest, should stay that way.)

[Ian Culverwell]

Annotated code with the above_highest_levels switch is attached.

[Ian Culverwell]

ropp_fm_refrac_1d_new.f90_above_highest_lev

[Ian Culverwell]

This solution would also be fairly easy to introduce into the TL and AD routines.

[Ian Culverwell]

Note that we don't have the same problem with the old refractivity interpolation scheme, because that just interpolates/extrapolates log N from the state levels to the refrac levels via a call to ropp_fm_interpol_log. This is why it can happily generate refracs of ~ 1e-53 N-units at 750 km. This may not really be desired. We could add a fix to ropp_fm/refrac_1d/ropp_fm_refrac.f90 to set those to ropp_MDFV above model top too. There's already a fix to do that for refracs below the bottom model level.

[Ian Culverwell]

above_highest_levels has been introduced at r4909.

By temporarily overwriting ropp_fm/data/FASCOD_non_ideal.nc with Hans's input file, and hacking ropp_fm/tests/t_fascod_tl.f90 and ropp_fm/tests/t_fascod_ad.f90 so that the output levels go above the model top at ~75 km, we can test the TL and AD parts. We get:

 t_fascod_tl -comp -new_op      
  
 Checking REFRACTIVITY TL ECMWF, CAT: 2
 |dy|/|dy_init|           cos_angle              100*rel diff
    1.0000000000E+00    9.8087181900E-01    4.4603720552E+01
    1.0000000000E-01    9.9630549549E-01    1.2352855809E+01
    1.0000000000E-02    9.9689573447E-01    9.3601495481E+00
    1.0000000000E-03    9.9849092061E-01    6.1777441480E+00
    1.0000000000E-04    9.9946475341E-01    3.4267647275E+00
    1.0000000000E-05    9.9987676082E-01    1.5974507168E+00
    1.0000000000E-06    9.9999444697E-01    3.3734494288E-01
    1.0000000000E-07    9.9999992628E-01    3.8842578450E-02
    1.0000000000E-08    9.9999999924E-01    3.9462068024E-03
    1.0000000000E-09    9.9999999999E-01    4.1882172800E-04
    1.0000000000E-10    9.9999999964E-01    2.7371053586E-03
    1.0000000000E-11    9.9999999059E-01    1.4057946494E-02
    1.0000000000E-12    9.9999873925E-01    1.6235754181E-01
    1.0000000000E-13    9.9994472870E-01    1.0704744077E+00
    1.0000000000E-14    9.9999571588E-01    3.0563895230E-01
    Refractivity TL check passed
  
 Checking BENDING ANGLE TL ECMWF, CAT: 2
 |dy|/|dy_init|           cos_angle              100*rel diff
    1.0000000000E+00    3.3020723155E-02    9.9999999889E+01
    1.0000000000E-01    4.4678573738E-02    9.9999999981E+01
    1.0000000000E-02    1.8192345326E-01    9.9999999983E+01
    1.0000000000E-03    7.5606922664E-01    9.9999999980E+01
    1.0000000000E-04    7.1463362890E-01    2.6271103738E+02
    1.0000000000E-05    9.4589600557E-01    1.2745866023E+02
    1.0000000000E-06    9.9734989466E-01    2.8674053090E+01
    1.0000000000E-07    9.9995658884E-01    3.6737256744E+00
    1.0000000000E-08    9.9999953385E-01    3.8062012546E-01
    1.0000000000E-09    9.9999999504E-01    3.8940516190E-02
    1.0000000000E-10    9.9999999972E-01    7.9919117385E-03
    1.0000000000E-11    9.9999963395E-01    3.3049124887E-01
    1.0000000000E-12    9.9997291704E-01    3.0544702994E+00
    1.0000000000E-13    9.9796260964E-01    2.5340990210E+01
    1.0000000000E-14    2.6737488182E-01    4.6467898497E+02
    Bending angle TL check passed
****************************
**********  PASS  **********
****************************

and

 t_fascod_ad -comp -new_op
 Checking REFRACTIVITY AD ECMWF, CAT: 2
 ad, x%non_ideal, x%new_ref_op, x%new_bangle_op =  T T T
norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 =   560427191.182  560427191.182 1.0000000    0.00000E+00
 Checking BENDING ANGLE AD ECMWF, CAT: 2
norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 =        1358.642       1358.642 1.0000000   -0.66941E-15
****************************
**********  PASS  **********
****************************

And I checked that by deliberately doing it wrong, by omitting the lines

      IF (y%geop(i) > x%geop(x%n_lev)) THEN
        above_highest_lev=.TRUE.
        EXIT
      END IF

from ropp_fm_refrac_1d_new_tl.f90 and ropp_fm_refrac_1d_new_ad.f90, these tests fail, either by generating NaNs in the TL tests or a seg fault in the AD tests.

Looking good, so passing to CB for review.

[cburrows]

The code change looks fine and the tests show that the TL and AD are consistent in cases where there are output levels above the top level in the background profile.

As noted, we should probably add a similar check for the "old" version at some point.

Chris

[Ian Culverwell]

The attached file does the analogous thing for the old refrac scheme. (Not committed until the corresponding changes to the TL and AD have been implemented.) It (ropp_fm_refrac_1d.f90) generates the following refracs:

refrac(773)=0.00458012287887 
refrac(774)=0.00450250270677 
refrac(775)=0.00442619797779 
refrac(776)=-99999000 
refrac(777)=-99999000 
refrac(778)=-99999000 

which show the right sort of differences to the results of the new operator (see above).

[Ian Culverwell]

ropp_fm_refrac_1d.f90_above_highest_lev

[Ian Culverwell]

The corresponding changes to the old refractivity operator have been made at r4914. Note that this includes some rationalisation of the adjoints of the 1d refrac and bangle operators, to nullify y_ad if the weights are zero, as they (now) are where the obs levels are outside the model domain. A change to the 1dvar routine ropp_1dvar_cost.f90 entails.

[Ian Culverwell]

This fix has caused no problems, so closing the ticket.