[Wrf-users] Bug in REAL interpolation code

[Kemp, Eric M (IS)]

Dear WRF Community:

I have found a subtle bug in program REAL (released with WRF 3.1) that can be
triggered when processing closely packed pressure levels from WPS (in my case,
high altitude pressure levels from 5 to 0.2 mb).

In subroutine vert_interp (dyn_em/module_initialize_real.F), pressure data
is copied from 3D array porig to 1D array ordered_porig for each column.
When there are no underground pressure levels (lines 2667-2716),
the program skips pressure levels that are "too close" to each other, where
"too close" is defined by namelist variable zap_close_levels (default is
500 Pa):

               DO ko = knext , generic
                  IF ( ordered_porig(count-1) - porig(i,ko,j) .LT. zap_close_levels ) THEN
                     zap = zap + 1
                     zap_above = zap_above + 1
                     CYCLE
                  END IF
                  ordered_porig(count) = porig(i,ko,j)
                  ordered_forig(count) = forig(i,ko,j)
                  count = count + 1
               END DO

Note that this can cause the level at the top of the atmosphere to be
"zapped".  The logic for the case with underground pressure levels (lines
2582-2667) is different, in that "zapping" only occurs underground and
immediately above the surface.  (Code omitted).

This causes problems when the program defines the subset of
ordered_porig that contains valid data (lines 2762-2795).  If both
underground pressure levels and the surface layer are used (lines 2762-2744),
the program will loop through ordered_porig and stop when (1) it finds the
pressure corresponding to the top of the atmosphere, or (2) when it reaches
the end of the array (variable generic).  Condition (1) is supposed to
account for "zapped" levels, but if the top level pressure isn't copied to
ordered_porig, then condition (2) will result:

               count = 0
               find_how_many_1 : DO ko = 1 , generic
                  IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN
                     count = count + 1
                     EXIT find_how_many_1
                  ELSE
                     count = count + 1
                  END IF
               END DO find_how_many_1
               kinterp_start = 1
               kinterp_end = kinterp_start + count - 1

Similar code is used when underground pressures are used but the surface
layer is NOT (lines 2744-2777):

               count = 0
               find_how_many_2 : DO ko = 1 , generic
                  IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN
                     count = count + 1
                     EXIT find_how_many_2
                  ELSE
                     count = count + 1
                  END IF
               END DO find_how_many_2
               kinterp_start = 1
               kinterp_end = kinterp_start + count - 1

If condition (2) holds when data was "zapped", then invalid elements of
ordered_porig will be passed to the interpolation code (lines 2799-2811),
and cause bogus results.  In some of my runs, ordered_porig contained zeros
in those invalid elements, and caused NaNs/floating point exceptions when
interpolating in log(pressure).

            IF ( interp_type .EQ. 1 ) THEN
               CALL lagrange_setup ( var_type , &
                    ordered_porig(kinterp_start:kinterp_end) , &
                    ordered_forig(kinterp_start:kinterp_end) , &
                    count , lagrange_order , extrap_type , &
                    ordered_pnew(kstart:kend)  , ordered_fnew  , kend-kstart+1 ,i,j)                                                  
            ELSE
               CALL lagrange_setup ( var_type , &
                LOG(ordered_porig(kinterp_start:kinterp_end)) , &
                    ordered_forig(kinterp_start:kinterp_end) , &
                    count , lagrange_order , extrap_type , &
                LOG(ordered_pnew(kstart:kend)) , ordered_fnew  , kend-kstart+1 ,i,j)                                                  
            END IF

I see two ways to fix this bug:  either restrict "zapping" to near the
surface, or change the logic in the "subset" code to use the zap variable
that keeps track of "zapped" levels.  I think the former is preferable,
since "zapping" pressure levels aloft can make it difficult to raise
the model top above 10 mb.

-Eric

Eric M. Kemp
Meteorologist
Northrop Grumman Information Systems
Advisory Services Division (TASC)
4801 Stonecroft Boulevard
Chantilly, VA 20151
301-286-9768              (NASA Goddard Office)
703-633-8300 x7078 (Chantilly SLEET Lab)
703-633-8300 x8278 (Chantilly Office)
703-449-3400              (Chantilly Fax)

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/wrf-users/attachments/20090730/7630981d/attachment.html