[Dart-dev] [6885] DART/trunk/models/mpas_atm: add in code to support interpolation of locations outside the

nancy at ucar.edu nancy at ucar.edu
Fri Apr 11 15:41:08 MDT 2014


Revision: 6885
Author:   nancy
Date:     2014-04-11 15:41:07 -0600 (Fri, 11 Apr 2014)
Log Message:
-----------
add in code to support interpolation of locations outside the
mpas grid but very close to the bottom or top levels.  allow
user to say how close is close enough with a namelist item.
leave a placeholder for real extrapolation, but for now just
evaluate the locations as if they are on level 1 or N.

this satisfies a functionality soyoung wanted when she needed
to evaluate cell edges which were higher or lower than their
corresponding cell centers. 

also support an rbf width of 0, meaning only use the vertices 
in the current cell to compute the value.  

both of these functions may be most useful for debugging, but
still seem of value.

Modified Paths:
--------------
    DART/trunk/models/mpas_atm/model_mod.f90
    DART/trunk/models/mpas_atm/model_mod.html
    DART/trunk/models/mpas_atm/model_mod.nml

-------------- next part --------------
Modified: DART/trunk/models/mpas_atm/model_mod.f90
===================================================================
--- DART/trunk/models/mpas_atm/model_mod.f90	2014-04-11 21:20:50 UTC (rev 6884)
+++ DART/trunk/models/mpas_atm/model_mod.f90	2014-04-11 21:41:07 UTC (rev 6885)
@@ -191,6 +191,13 @@
 ! assimilation, and use only the increments to update the edge winds
 logical :: use_increments_for_u_update = .true.
 
+! if > 0, amount of distance in fractional model levels that a vertical
+! point can be above or below the top or bottom of the grid and still be
+! evaluated without error.  if extrapolate is true, extrapolate from the
+! first or last model level.  if extrapolate is false, use level 1 or N.
+real(r8) :: outside_grid_level_tolerance = -1.0_r8
+logical  :: extrapolate = .false.
+
 namelist /model_nml/             &
    model_analysis_filename,      &
    grid_definition_filename,     &
@@ -208,6 +215,8 @@
    update_u_from_reconstruct,    &
    use_increments_for_u_update,  &
    highest_obs_pressure_mb,      &
+   outside_grid_level_tolerance, &
+   extrapolate,                  &
    sfc_elev_max_diff
 
 ! DART state vector contents are specified in the input.nml:&mpas_vars_nml namelist.
@@ -741,6 +750,12 @@
 endif
 
 
+if (extrapolate) then
+   call error_handler(E_MSG,'static_init_model',&
+                      'extrapolate not supported yet; will use level 1 or N values', &
+                      source, revision, revdate)
+endif
+
 allocate( ens_mean(model_size) )
 
 end subroutine static_init_model
@@ -4605,6 +4620,7 @@
 ! models get to be huge.
 
 integer :: i
+real(r8) :: hlevel
 
 ! Initialization
 ier = 0
@@ -4612,9 +4628,51 @@
 lower = -1
 upper = -1
 
-! Bounds check
-if(height < bounds(1)) ier = 998
-if(height > bounds(nbounds)) ier = 9998
+! FIXME:
+! difficult to allow extrapolation because this routine only returns
+! an upper and lower level, plus a fraction.  it would have to return
+! a fraction < 0 or > 1 and then the calling code would have to test for
+! and support extrapolation at some later point in the code.  so ignore
+! extrapolation for now.
+!
+! do try to convert heights to grid levels and if within the limits
+! return the top or bottom level as the value.  convert lower locations
+! based on bounds(1)/(2) and upper on bounds(nbounds-1)/bounds(nbounds)
+
+! Bounds checks:
+
+! too low?
+if(height < bounds(1)) then
+   if(outside_grid_level_tolerance <= 0.0_r8) then
+      ier = 998
+   else
+      ! let points outside the grid but "close enough" use the
+      ! bottom level as the height.  hlevel should come back < 1
+      hlevel = extrapolate_level(height, bounds(1), bounds(2))
+      if (hlevel + outside_grid_level_tolerance >= 1.0_r8) then
+         fract = 0.0_r8
+         lower = 1
+         upper = 2
+      endif
+   endif
+endif
+
+! too high?
+if(height > bounds(nbounds)) then
+   if(outside_grid_level_tolerance <= 0.0_r8) then
+      ier = 9998
+   else
+      ! let points outside the grid but "close enough" use the
+      ! top level as the height.  hlevel should come back > nbounds
+      hlevel = extrapolate_level(height, bounds(nbounds-1), bounds(nbounds))
+      if (hlevel - outside_grid_level_tolerance <= real(nbounds, r8)) then
+         fract = 1.0_r8
+         lower = nbounds-1
+         upper = nbounds
+      endif
+   endif
+endif
+
 if (ier /= 0) return
 
 ! Search two adjacent layers that enclose the given point
@@ -5441,6 +5499,12 @@
 edge_list = 0
 
 select case (use_rbf_option)
+  case (0)
+      ! use edges on the closest cell only.
+      nedges = nEdgesOnCell(cellid)
+      do i=1, nedges
+         edge_list(i) = edgesOnCell(i, cellid)
+      enddo
   case (1)
       ! fill in the number of unique edges and fills the
       ! edge list with the edge ids.  the code that detects
@@ -6095,6 +6159,42 @@
 
 !------------------------------------------------------------
 
+function extrapolate_level(h, lb, ub)
+
+! Given a test height and two level heights, extrapolate the
+! actual level number which would correspond to that height.
+
+real(r8), intent(in) :: h
+real(r8), intent(in) :: lb
+real(r8), intent(in) :: ub
+real(r8)             :: extrapolate_level
+
+real(r8) :: thickness, fract
+
+thickness = ub - lb
+
+if (h < lb) then
+   ! extrapolate down - not in data values but in level
+   ! relative to the height difference in levels 1-2
+   extrapolate_level = lb - (((lb - h) + lb) / thickness)
+   
+else if (h > ub) then
+   ! extrapolate up - not in data values but in level
+   ! relative to the height difference in levels nb-1 to nb
+   extrapolate_level = ub + ((ub - (h - ub)) / thickness)
+
+else
+   write(string1, *) h, ' not outside range of ', lb, ' and ', ub
+   call error_handler(E_ERR, 'extrapolate_level', &
+                      'extrapolate code called for height not out of range', &
+                      source, revision, revdate, text2=string1)
+endif
+
+if (debug > 1) print *, 'extrapolate: ', h, lb, ub, extrapolate_level
+end function extrapolate_level
+
+!------------------------------------------------------------
+
 subroutine latlon_to_xyz(lat, lon, x, y, z)
 
 ! Given a lat, lon in degrees, return the cartesian x,y,z coordinate

Modified: DART/trunk/models/mpas_atm/model_mod.html
===================================================================
--- DART/trunk/models/mpas_atm/model_mod.html	2014-04-11 21:20:50 UTC (rev 6884)
+++ DART/trunk/models/mpas_atm/model_mod.html	2014-04-11 21:41:07 UTC (rev 6885)
@@ -131,6 +131,8 @@
    use_increments_for_u_update  = .true.,
    highest_obs_pressure_mb      = 100.0,
    sfc_elev_max_diff            = -1.0,
+   outside_grid_level_tolerance = -1.0,
+   extrapolate                  = .false.,
    debug                        = 0,
 /
 </pre>
@@ -245,8 +247,10 @@
 <TR><TD>use_rbf_option</TD>
     <TD>integer <em class=units>[default:&nbsp;2]</em></TD>
     <TD>If <em class=code>use_u_for_wind</em> is .true., this option controls
-        how many points will be used in the RBF interpolation. Options are available as 1, 2, and 3.
-        All the edges available in N (= 1,2, or 3) cells go into the RBF reconstruction. </TD>
+        how many points will be used in the RBF interpolation. 
+        Options are available as 0, 1, 2, and 3.
+        All the edges available in N (= 0,1,2, or 3) neighboring cells 
+        go into the RBF reconstruction. </TD>
 </TR>
 
 <TR><TD>update_u_from_reconstruct</TD>
@@ -288,6 +292,36 @@
     <TD> Character string specifying the calendar being used by MPAS.</TD>
 </TR>
 
+<TR><TD>outside_grid_level_tolerance</TD>
+    <TD>real(r8) <em class=units>[default:&nbsp;-1.0]</em></TD>
+    <TD>If greater than 0.0, amount of distance in fractional model levels 
+        that a vertical location can be above or below the top or bottom of the 
+        grid and still be evaluated without error.  
+        Since <em class=code>extrapolate</em> is not implemented yet,
+        the value of <em class=code>.false.</em> will be assumed.
+        In this case, vertical locations equivalent to level 1
+        or level N will be used.  Eventually, if
+        <em class=code>extrapolate</em> is <em class=code>.true.</em>, 
+        extrapolate from the first or last model level.  
+        If <em class=code>extrapolate</em> is <em class=code>.false.</em>, 
+        simply use the value at level 1 for low vertical locations, 
+        or at level N for high vertical locations. </TD>
+</TR> 
+
+<TR><TD>extrapolate</TD>
+    <TD>logical <em class=units>[default:&nbsp;.false.]</em></TD>
+    <TD><em>NOT IMPLEMENTED YET</em>.  Vertical locations equivalant to level 1
+        or level N will be used.  When this is implemented, it will do: <br />
+        If <em class=code>outside_grid_level_tolerance</em> is greater than 0.0,
+        then control how values are assigned to locations where the vertical is
+        exterior to the grid.  If this is <em class=code>.true.</em>, then
+        extrapolate low locations from levels 1 and 2, and high locations from
+        levels N-1 and N.  If this is <em class=code>.false.</em>, then simply 
+        use the corresponding values at level 1 or N.  This item is ignored if 
+        <em class=code>outside_grid_level_tolerance</em> is less than or
+        equal to 0.0. </TD>
+</TR> 
+
 <TR><TD>debug</TD>
     <TD>integer <em class=units>[default:&nbsp;0]</em></TD>
     <TD>The switch to specify the run-time verbosity.

Modified: DART/trunk/models/mpas_atm/model_mod.nml
===================================================================
--- DART/trunk/models/mpas_atm/model_mod.nml	2014-04-11 21:20:50 UTC (rev 6884)
+++ DART/trunk/models/mpas_atm/model_mod.nml	2014-04-11 21:41:07 UTC (rev 6885)
@@ -14,6 +14,8 @@
    use_increments_for_u_update  = .true.,
    highest_obs_pressure_mb      = 100.0,
    sfc_elev_max_diff            = -1.0,
+   outside_grid_level_tolerance = -1.0,
+   extrapolate                  = .false.,
    debug                        = 0,
 /
 


More information about the Dart-dev mailing list