Re: Bug in vinth2p?

All,

Mark discovered a bug that affects the top two interpolated levels
only in the vinth2p, vinth2p_ecmwf, and vintp2p_ecmwf functions. We
fixed this bug in V5.1.1.

Mark, thanks for taking the time to look into this and report it!

Dennis Shea provided some information for me to include here.

The underlying codes were taken *directly* from the original CCM
Processor:

    http://www.cgd.ucar.edu/cms/processor/subject/vertical.interp.html

These subroutines have been used possibly since the late 1980s.

An example will be used to illustrate the bug:

Consider temperatures [TM (K)] at the following model hybrid levels
[PM (hPa)]:

            PM TM <=== Raw Data
          3.64 226.77
          7.59 225.53
         14.36 224.89
         24.61 221.39
         38.27 213.97
[snip]

Let's say that interpolation to the following pressure levels (PI) is
desired:

         PI = (/ 5, 7, 10, 20, 30, 50, ...... /)

The following shows the values interpolated by the vinth2p and
vinth2p_ecmwf functions before (v5.1.0) and after (v5.1.1) the bug
fix:

          PI v5.1.0 v5.1.1 diff
           5 225.78 226.34 0.564
           7 225.59 225.72 0.129
          10 225.30 225.30 0.000 No differences
          20 222.96 222.96 0.000
          30 218.46 218.46 0.000
[snip]

The differences could be 'small' or 'large'. It depends upon the
vertical distribution of the variable.

--Mary [using input from Dennis]

On Thu, 28 May 2009, Mark Branson wrote:

>
> I was taking a close look at the vinth2p routine recently and I found
> something odd about how it handles the vertical interpolation for the
> very top output pressure levels.
>
> I wrote a test ncl program that reads in a zonal wind field from a CAM
> history file and calls vinth2p to interpolate them. Then I created an
> analogous fortran program that calls the ni/src/lib/nfpfort/
> vinth2p_dp.f subroutine and it gives the same output values as the ncl
> test code.
>
> Next I put a print statement inside vinth2p_dp.f just to check the
> input and output data values, and I found that it seems to do an
> extrapolation for any plevo values that fall between plevi(1) and
> plevi(2). For example, for one particular grid column in my sample
> data, I have this:
>
> plevi = 3.64, 7.59, 14.36, 24.6, ...
> plevo = 5., 7., 10., 20., ...
>
> You would think that for plevo=5 and 7 you would use the plevi values
> at 3.64 and 7.59 for the interpolation, but it's using 7.59 and
> 14.36. Then it correctly uses those for plevo=10, and then moves down
> to 14.36 and 24.6 for plevo=20.
>
> This is the pertinent piece of code:
>
> C
> C**** IF BRANCH FOR MODEL INTERIOR
> C**** LOOP THROUGH INPUT LEVELS TILL YOU ARE BRACKETING
> C**** OUTPUT LEVEL
> C
> ELSE
> KP = 1
> 20 CONTINUE
> KP = KP + 1
> IF (PLEVO(K).LE.PLEVI(KP+1)) GO TO 30
> IF (KP.GT.NLEVI) THEN
> WRITE (6,FMT=25) KP,NLEVI
> 25 FORMAT (' KP.GT.KLEVI IN P2HBD. KP,KLEVI= ',
> + 2I5)
> C CALL PRNT(' PLEVI',PLEVI,NLEVIP1,1,1,1)
> C CALL PRNT(' PLEVO',PLEVO,NLEVO,1,1,1)
> C CALL ABORT(' KP.GT.NLEVI IN P2HBD')
> END IF
> GO TO 20
> END IF
> 30 CONTINUE
>
> I think the KP = 1 statement should instead by KP = 0. Because when
> you first drop into it you will then set KP = 2 when you do the
> statement just below "20 continue" and then the if test is checking if
> plevo(k) <= plevi(3), which makes no sense.
>
> I'm probably missing something obvious, but if so I'm not seeing it
> right now. Hopefully someone can set me straight.
>
> Thanks,
> Mark
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Jun 18 2009 - 07:40:49 MDT