Bug in vinth2p?

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
Received on Thu May 28 2009 - 09:55:18 MDT