[Dart-dev] DART/branches Revision: 13206

[dart at ucar.edu]

nancy at ucar.edu
2019-07-25 09:52:03 -0600 (Thu, 25 Jul 2019)
887
add an option to allow a 1D latitude array to 
start at 90 and to go -90 instead of starting 
at -90 and going to 90.  this entails iterating
the 1d array in the opposite order, but leaves
the rest of the code alone because the return values
are the corner indices of the right quad, and the
fraction across the quad.  as long as we're returning
the right corner indicies, the interpolation code
doesn't care which order the latitudes are in.

also cleaned up some comments and start to explore
exactly what would happen with a grid that is specified
as -180 to 180 in longitudes.  i believe adding 360 to
all the negative values makes this valid, since we
already have code that knows how to interpolate at
the "seam" between a line of grid cells that ends
around 360 on the left and starts around 0 on the right.
(add 360 to the "0" side and compute the fraction across
like normal).




Modified: DART/branches/basic_earth/models/utilities/quad_utils_mod.f90
===================================================================
--- DART/branches/basic_earth/models/utilities/quad_utils_mod.f90	2019-07-23 22:48:07 UTC (rev 13205)
+++ DART/branches/basic_earth/models/utilities/quad_utils_mod.f90	2019-07-25 15:52:03 UTC (rev 13206)
@@ -77,7 +77,7 @@
 
 logical, save :: module_initialized = .false.
 
-integer  :: debug = 0   ! turn up for more and more debug messages
+integer  :: debug = 0               ! turn up for more and more debug messages
 integer  :: interpolation_type = 1  ! add cases for different strategies
 logical  :: do_rotate = .false.     ! rotate edge from pts 1,2 to horizontal before interp
 
@@ -105,8 +105,7 @@
 integer, parameter :: QUAD_LOCATED_LAT_EDGES     =   3
 integer, parameter :: QUAD_LOCATED_CELL_CORNERS  =   4
 
-! data struct question - does this go directly into handle?
-! right now it's in each option.
+! options that control how the interpolation handles specific cases
 type quad_grid_options
    private
 
@@ -119,12 +118,35 @@
    ! if not, either or both could still be true but we can't assume.
    logical :: global_grid = .false.
 
-   ! does the grid cross the 360 -> 0 boundary?
-   logical :: spans_lon_zero = .false.   ! true if any lon(:,1) > lon(:,nlons)
+   ! separate out the single logical 'spans' flag into two cases?
+   ! case 1:  the longitude grid is cyclic; either a global grid or a band between 
+   !  two latitude lines that circles the globe.  all longitude values are valid.
+   ! case 2: a regional grid that crosses the prime meridian.
+   !  this will contain a discontinuity around 360 -> 0 that should be a valid region.
+   ! can a single flag handle both of these cases?  
 
+   ! are there valid values between array(N) and array(1) that should be interpolated?
+   ! always true for global grids. 
+   logical :: lon_cyclic = .false.  
+
+   ! are there valid values between array(X) and array(X+1) in the interior of
+   ! the longitude array where a(X) > a(X+1) (e.g. around where a(X) is near 360 
+   ! and a(X+1) is near 0) that should be interpolated?
+
+   ! always true for global grids. 
+   ! for partially regular grids (1D lon and 1D lat arrays) this can be
+   ! true if a regional grid crosses the prime meridian.
+   ! for irregular grids (2D lons and 2D lats) this is true if 
+   ! any lon(:,1) > lon(:,nlons)
+   logical :: spans_lon_zero = .false.   
+
    ! do we handle wrap over the poles?
    logical :: pole_wrap = .false.
 
+   ! are the latitudes specified from smallest (south) to largest (north)
+   ! or largest (north) to smallest (south)?
+   logical :: north_to_south = .false.
+
    ! i don't want to know this, but apparently we might
    ! have to know if the points given are cell-centered or
    ! on the edges.  (god forbid we need to know U stagger vs V stagger!)
@@ -348,6 +370,14 @@
                interp_handle%ir%lons_1D(num_lons))
       interp_handle%ir%lats_1D(num_lats) = MISSING_R8
       interp_handle%ir%lons_1D(num_lons) = MISSING_R8
+      if (debug > 10) then
+         write(string1, *) 'nlats, nlons: ', num_lats, num_lons
+         call log_it(string1)
+         write(string1, *) 'first lat: ', interp_handle%ir%lats_1D(1)
+         call log_it(string1)
+         write(string1, *) 'first lon: ', interp_handle%ir%lons_1D(1)
+         call log_it(string1)
+      endif
 
    case(GRID_QUAD_FULLY_IRREGULAR)
       allocate(interp_handle%ii%lats_2D(num_lons, num_lats), &
@@ -560,9 +590,47 @@
                       source, revision, revdate)
 endif
 
+! lons and lats are declared intent(in).  any code that just needs to
+! test values can use them.  any code that is going to modify values
+! has to use the xxx_1D arrays in the structures.
+
 interp_handle%ir%lons_1D(:) = lons
 interp_handle%ir%lats_1D(:) = lats
 
+! inverted order for latitudes?
+if (lats(1) > lats(interp_handle%nlat)) then
+   interp_handle%opt%north_to_south = .true.
+endif
+
+! -180 to 180 instead of 0 to 360?  add 360, which makes all
+! the values valid (we require longitudes to be between 0 and 360)
+! but it also makes a partially regular grid start at 180, go up
+! to 360, then back to 0, then up to 180.  set the 'spans 0' flag.
+if (any(lons < 0.0_r8)) then
+  where(interp_handle%ir%lons_1D < 0.0_r8) &
+     interp_handle%ir%lons_1D = interp_handle%ir%lons_1D + 360.0_r8
+     interp_handle%opt%spans_lon_zero = .true. 
+endif