[Wrf-users] Bug in REAL interpolation code

[wrfhelp]

An alternative workaround is to set zap_close_levels = 0. We will also  
be looking at a fix so that
zap_close_levels will not affect upper levels.

wrfhelp

On Jul 30, 2009, at 1:26 PM, Kemp, Eric M (IS) wrote:

> 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)
>
>


wrfhelp
http://www.mmm.ucar.edu/wrf/users/supports/wrfhelp.html