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
This archive was generated by hypermail 2.2.0 : Wed Jun 03 2009 - 12:52:48 MDT