<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
<HTML>
<HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso-8859-1">
<META NAME="Generator" CONTENT="MS Exchange Server version 6.5.7654.12">
<TITLE>Bug in REAL interpolation code</TITLE>
</HEAD>
<BODY>
<!-- Converted from text/plain format -->
<BR>
<BR>
<P><FONT SIZE=2>Dear WRF Community:<BR>
<BR>
I have found a subtle bug in program REAL (released with WRF 3.1) that can be<BR>
triggered when processing closely packed pressure levels from WPS (in my case,<BR>
high altitude pressure levels from 5 to 0.2 mb).<BR>
<BR>
In subroutine vert_interp (dyn_em/module_initialize_real.F), pressure data<BR>
is copied from 3D array porig to 1D array ordered_porig for each column.<BR>
When there are no underground pressure levels (lines 2667-2716),<BR>
the program skips pressure levels that are "too close" to each other, where<BR>
"too close" is defined by namelist variable zap_close_levels (default is<BR>
500 Pa):<BR>
<BR>
DO ko = knext , generic<BR>
IF ( ordered_porig(count-1) - porig(i,ko,j) .LT. zap_close_levels ) THEN<BR>
zap = zap + 1<BR>
zap_above = zap_above + 1<BR>
CYCLE<BR>
END IF<BR>
ordered_porig(count) = porig(i,ko,j)<BR>
ordered_forig(count) = forig(i,ko,j)<BR>
count = count + 1<BR>
END DO<BR>
<BR>
Note that this can cause the level at the top of the atmosphere to be<BR>
"zapped". The logic for the case with underground pressure levels (lines<BR>
2582-2667) is different, in that "zapping" only occurs underground and<BR>
immediately above the surface. (Code omitted).<BR>
<BR>
This causes problems when the program defines the subset of<BR>
ordered_porig that contains valid data (lines 2762-2795). If both<BR>
underground pressure levels and the surface layer are used (lines 2762-2744),<BR>
the program will loop through ordered_porig and stop when (1) it finds the<BR>
pressure corresponding to the top of the atmosphere, or (2) when it reaches<BR>
the end of the array (variable generic). Condition (1) is supposed to<BR>
account for "zapped" levels, but if the top level pressure isn't copied to<BR>
ordered_porig, then condition (2) will result:<BR>
<BR>
count = 0<BR>
find_how_many_1 : DO ko = 1 , generic<BR>
IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN<BR>
count = count + 1<BR>
EXIT find_how_many_1<BR>
ELSE<BR>
count = count + 1<BR>
END IF<BR>
END DO find_how_many_1<BR>
kinterp_start = 1<BR>
kinterp_end = kinterp_start + count - 1<BR>
<BR>
Similar code is used when underground pressures are used but the surface<BR>
layer is NOT (lines 2744-2777):<BR>
<BR>
count = 0<BR>
find_how_many_2 : DO ko = 1 , generic<BR>
IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN<BR>
count = count + 1<BR>
EXIT find_how_many_2<BR>
ELSE<BR>
count = count + 1<BR>
END IF<BR>
END DO find_how_many_2<BR>
kinterp_start = 1<BR>
kinterp_end = kinterp_start + count - 1<BR>
<BR>
If condition (2) holds when data was "zapped", then invalid elements of<BR>
ordered_porig will be passed to the interpolation code (lines 2799-2811),<BR>
and cause bogus results. In some of my runs, ordered_porig contained zeros<BR>
in those invalid elements, and caused NaNs/floating point exceptions when<BR>
interpolating in log(pressure).<BR>
<BR>
IF ( interp_type .EQ. 1 ) THEN<BR>
CALL lagrange_setup ( var_type , &<BR>
ordered_porig(kinterp_start:kinterp_end) , &<BR>
ordered_forig(kinterp_start:kinterp_end) , &<BR>
count , lagrange_order , extrap_type , &<BR>
ordered_pnew(kstart:kend) , ordered_fnew , kend-kstart+1 ,i,j) <BR>
ELSE<BR>
CALL lagrange_setup ( var_type , &<BR>
LOG(ordered_porig(kinterp_start:kinterp_end)) , &<BR>
ordered_forig(kinterp_start:kinterp_end) , &<BR>
count , lagrange_order , extrap_type , &<BR>
LOG(ordered_pnew(kstart:kend)) , ordered_fnew , kend-kstart+1 ,i,j) <BR>
END IF<BR>
<BR>
I see two ways to fix this bug: either restrict "zapping" to near the<BR>
surface, or change the logic in the "subset" code to use the zap variable<BR>
that keeps track of "zapped" levels. I think the former is preferable,<BR>
since "zapping" pressure levels aloft can make it difficult to raise<BR>
the model top above 10 mb.<BR>
<BR>
-Eric<BR>
<BR>
Eric M. Kemp<BR>
Meteorologist<BR>
Northrop Grumman Information Systems<BR>
Advisory Services Division (TASC)<BR>
4801 Stonecroft Boulevard<BR>
Chantilly, VA 20151<BR>
301-286-9768 (NASA Goddard Office)<BR>
703-633-8300 x7078 (Chantilly SLEET Lab)<BR>
703-633-8300 x8278 (Chantilly Office)<BR>
703-449-3400 (Chantilly Fax)<BR>
<BR>
</FONT>
</P>
</BODY>
</HTML>