Tweak interpolation of pressure in grib2bgrasc

[Ian Culverwell]

Kent Bækgaard Lauritsen writes:

We have re-processed some of our bgr-files from the operational production with
the new ROPP tools. When we compared these files to the old/current files we
noted that the surface pressure sometimes differ. This seems to be due to the
following different code logic:

1) Our operational code: reads ln(Psurf); interpolates; takes exp

2) ROPP code (i.e. grib2bgrasc): reads ln(Psurf); takes exp; interpolates

The relative difference between these two methods for two points P1 and P2 are
(to lowest order):  a = (1/8)*((P1-P2)/P1)^2

We obtain values like a = 0.1-0.5 mbar / 1000 mbar  (in about 10% of the cases).
I.e., P1-P2 = 10-20 hPa. Note, a can also be e.g. 2 mbar (i.e., a 40 hPa
difference).

These values could be from regions with strong gradients (e.g. near a low
pressure region) but more likely it is probably from regions near "mountains":

height diff / 8000 m = delta P / P = 40 / 1000,

i .e., it corresponds to a height diff of 320 m.

Since the pressure more or less decreases exponentially, it seems that one
should interpolate in ln(Psurf) and only afterwards take the exp in order to
avoid making a systematic (and sometimes large) interpolation error.

Therefore, we recommend that this code change is made in ROPP-8 such that ROPP
will follow the logic in our operational code.

We will implement this change in our DMI-ROPP (dmi_trunk) already now. As far as
we can see in the code (grib2bgrasc, attached) it is straightforward and simply
consists of moving the exp from line 882:

  P0_grb(ilat,:) = EXP(values(from:to))  ! convert to P* in Pa 

to line 931-932 (i.e. line 931):

P0_prf  = (1.0d0-tt)*(1.0d0-uu)*P0_grb(  ilat,  ilon)  + &    
                  tt*(1.0d0-uu)*P0_grb(  ilat,ilonp1)  + & 
                  tt*uu*P0_grb(ilatp1,ilonp1)  + &
          (1.0d0-tt)*uu*P0_grb(ilatp1,  ilon)

Please confirm if this is correct.

Interpolating log P rather than P makes sense to me. Agreed to include it in ROPP8.0.

[Ian Culverwell]

Done - see r4123.

[Ian Culverwell]

Close ticket.