[Dart-dev] [3969] DART/trunk/models/POP: Integrated all Jeff' s dipole interpolation code.
nancy at ucar.edu
nancy at ucar.edu
Tue Jul 21 15:52:57 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090721/12c4a84d/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/models/POP/model_mod.f90
===================================================================
--- DART/trunk/models/POP/model_mod.f90 2009-07-20 17:21:37 UTC (rev 3968)
+++ DART/trunk/models/POP/model_mod.f90 2009-07-21 21:52:56 UTC (rev 3969)
@@ -56,7 +56,8 @@
get_close_maxdist_init, &
get_close_obs_init, &
get_close_obs, &
- ens_mean_for_model
+ ens_mean_for_model, &
+ test_interpolation
! generally useful routines for various support purposes.
! the interfaces here can be changed as appropriate.
@@ -181,10 +182,14 @@
! Grid parameters - the values will be read from a
! standard POP namelist and filled in here.
+! nx, ny and nz are the size of the dipole (or irregular) grids.
integer :: Nx=-1, Ny=-1, Nz=-1 ! grid counts for each field
! locations of cell centers (C) and edges (G) for each axis.
real(r8), allocatable :: XC(:), XG(:), YC(:), YG(:), ZC(:), ZG(:)
+
+! These arrays store the longitude and latitude of the lower left corner of
+! each of the dipole u quadrilaterals and t quadrilaterals.
real(r8), allocatable :: ULAT(:,:), ULON(:,:), TLAT(:,:), TLON(:,:)
@@ -217,6 +222,47 @@
END INTERFACE
+!------------------------------------------------
+
+! The regular grid used for dipole interpolation divides the sphere into
+! a set of regularly spaced lon-lat boxes. The number of boxes in
+! longitude and latitude are set by num_reg_x and num_reg_y. Making the
+! number of regular boxes smaller decreases the computation required for
+! doing each interpolation but increases the static storage requirements
+! and the initialization computation (which seems to be pretty small).
+integer, parameter :: num_reg_x = 90, num_reg_y = 90
+
+! The max_reg_list_num controls the size of temporary storage used for
+! initializing the regular grid. Four arrays
+! of size num_reg_x*num_reg_y*max_reg_list_num are needed. The initialization
+! fails and returns an error if max_reg_list_num is too small. A value of
+! 30 is sufficient for the x3 POP grid with 180 regular lon and lat boxes
+! and a value of ??? is sufficient for for the x1 grid.
+integer, parameter :: max_reg_list_num = 30
+
+! The dipole interpolation keeps a list of how many and which dipole quads
+! overlap each regular lon-lat box. The number for the u and t grids are stored
+! in u_dipole_num and t_dipole_num. The allocatable arrays u_dipole_lon(lat)_list and
+! t_dipole_lon(lat)_list list the longitude and latitude indices for the overlapping
+! dipole quads. The entry in u_dipole_start and t_dipole_start for a given
+! regular lon-lat box indicates where the list of dipole quads begins
+! in the u_dipole_lon(lat)_list and t_dipole_lon(lat)_list arrays.
+integer :: u_dipole_start(num_reg_x, num_reg_y)
+integer :: u_dipole_num (num_reg_x, num_reg_y) = 0
+integer :: t_dipole_start(num_reg_x, num_reg_y)
+integer :: t_dipole_num (num_reg_x, num_reg_y) = 0
+integer, allocatable :: u_dipole_lon_list(:), t_dipole_lon_list(:)
+integer, allocatable :: u_dipole_lat_list(:), t_dipole_lat_list(:)
+
+! Need to check for pole quads: for now we are not interpolating in them
+integer :: pole_x, t_pole_y, u_pole_y
+
+
+
+! Have a global variable saying whether this is dipole or regular lon-lat grid
+! This should be initialized static_init_model. Code to do this is below.
+logical :: dipole_grid
+
contains
!==================================================================
@@ -337,11 +383,356 @@
model_size = (n3dfields * (Nx * Ny * Nz)) + (n2dfields * (Nx * Ny))
if (do_output()) write(*,*) 'model_size = ', model_size
+! Initialize the interpolation routines
+call init_interp()
+
end subroutine static_init_model
+!------------------------------------------------------------
+subroutine init_interp()
+
+! Initializes data structures needed for POP interpolation for
+! either dipole or irregular grid.
+! This should be called at static_init_model time to avoid
+! having all this temporary storage in the middle of a run.
+
+integer :: i
+
+! Determine whether this is a irregular lon-lat grid or a dipole.
+! Do this by seeing if the lons have the same values at both
+! the first and last latitude row; this is not the case for dipole.
+
+dipole_grid = .false.
+do i = 1, nx
+ if(ulon(i, 1) /= ulon(i, ny)) then
+ dipole_grid = .true.
+ write(*, *) 'dipole grid is ', dipole_grid
+ call init_dipole_interp()
+ return
+ endif
+end do
+
+end subroutine init_interp
+
+
+!------------------------------------------------------------
+
+
+subroutine init_dipole_interp()
+
+! Build the data structure for interpolation for a dipole grid.
+
+! Need a temporary data structure to build this.
+! These arrays keep a list of the x and y indices of dipole quads
+! that potentially overlap the regular boxes. Need one for the u
+! and one for the t grid.
+integer :: ureg_list_lon(num_reg_x, num_reg_y, max_reg_list_num)
+integer :: ureg_list_lat(num_reg_x, num_reg_y, max_reg_list_num)
+integer :: treg_list_lon(num_reg_x, num_reg_y, max_reg_list_num)
+integer :: treg_list_lat(num_reg_x, num_reg_y, max_reg_list_num)
+
+
+real(r8) :: u_c_lons(4), u_c_lats(4), t_c_lons(4), t_c_lats(4), pole_row_lon
+integer :: i, j, k, pindex
+integer :: reg_lon_ind(2), reg_lat_ind(2), u_total, t_total, u_index, t_index
+logical :: is_pole
+
+! Begin by finding the quad that contains the pole for the dipole t_grid. To do this
+! locate the u quad with the pole on its right boundary. This is on the row
+! that is opposite the shifted pole and exactly follows a lon circle.
+pole_x = nx / 2;
+! Search for the row at which the longitude flips over the pole
+pole_row_lon = ulon(pole_x, 1);
+do i = 1, ny
+ pindex = i
+ if(ulon(pole_x, i) /= pole_row_lon) exit
+enddo
+
+! Pole boxes for u have indices pole_x or pole_x-1 and index - 1;
+! (it's right before the flip).
+u_pole_y = pindex - 1;
+
+! Locate the T dipole quad that contains the pole.
+! We know it is in either the same lat quad as the u pole edge or one higher.
+! Figure out if the pole is more or less than halfway along
+! the u quad edge to find the right one.
+if(ulat(pole_x, u_pole_y) > ulat(pole_x, u_pole_y + 1)) then
+ t_pole_y = u_pole_y;
+else
+ t_pole_y = u_pole_y + 1;
+endif
+
+! Loop through each of the dipole grid quads
+do i = 1, nx
+ ! There's no wraparound in y, one box less than grid boundaries
+ do j = 1, ny - 1
+ ! Is this the pole quad for the T grid?
+ is_pole = (i == pole_x .and. j == t_pole_y)
+
+ ! Set up an array of the lons and lats for the corners of these u and t quads
+ call get_quad_corners(ulon, i, j, u_c_lons)
+ call get_quad_corners(ulat, i, j, u_c_lats)
+ call get_quad_corners(tlon, i, j, t_c_lons)
+ call get_quad_corners(tlat, i, j, t_c_lats)
+
+ ! Get list of regular boxes that cover this u dipole quad
+ ! The false indicates that for the u grid there's nothing special about pole
+ call reg_box_overlap(u_c_lons, u_c_lats, .false., reg_lon_ind, reg_lat_ind)
+
+ ! Update the temporary data structures for the u quad
+ call update_reg_list(u_dipole_num, ureg_list_lon, &
+ ureg_list_lat, reg_lon_ind, reg_lat_ind, i, j)
+
+ ! Repeat for t dipole quads
+ call reg_box_overlap(t_c_lons, t_c_lats, is_pole, reg_lon_ind, reg_lat_ind)
+ call update_reg_list(t_dipole_num, treg_list_lon, &
+ treg_list_lat, reg_lon_ind, reg_lat_ind, i, j)
+ enddo
+enddo
+
+! Invert the temporary data structure. The total number of entries will be the sum
+! of the number of dipole cells for each regular cell.
+u_total = sum(u_dipole_num)
+t_total = sum(t_dipole_num)
+
+! Allocate storage for the final structures in module storage
+allocate(u_dipole_lon_list(u_total), u_dipole_lat_list(u_total))
+allocate(t_dipole_lon_list(t_total), t_dipole_lat_list(t_total))
+
+! Fill up the long list by traversing the temporary structure. Need indices
+! to keep track of where to put the next entry.
+u_index = 1
+t_index = 1
+
+! Loop through each regular grid box
+do i = 1, num_reg_x
+ do j = 1, num_reg_y
+
+ ! The list for this regular box starts at the current indices.
+ u_dipole_start(i, j) = u_index
+ t_dipole_start(i, j) = t_index
+
+ ! Copy all the close dipole quads for regular u box(i, j)
+ do k = 1, u_dipole_num(i, j)
+ u_dipole_lon_list(u_index) = ureg_list_lon(i, j, k)
+ u_dipole_lat_list(u_index) = ureg_list_lat(i, j, k)
+ u_index = u_index + 1
+ enddo
+
+ ! Copy all the close dipoles for regular t box (i, j)
+ do k = 1, t_dipole_num(i, j)
+ t_dipole_lon_list(t_index) = treg_list_lon(i, j, k)
+ t_dipole_lat_list(t_index) = treg_list_lat(i, j, k)
+ t_index = t_index + 1
+ enddo
+
+ enddo
+enddo
+
+! Confirm that the indices come out okay as debug
+if(u_index /= u_total + 1) then
+ msgstring = "Storage indices did not balance for U grid: : contact DART developers"
+ call error_handler(E_ERR, "init_dipole_interp", msgstring, source, revision, revdate)
+endif
+if(t_index /= t_total + 1) then
+ msgstring = "Storage indices did not balance for T grid: : contact DART developers"
+ call error_handler(E_ERR, "init_dipole_interp", msgstring, source, revision, revdate)
+endif
+
+end subroutine init_dipole_interp
+
+!------------------------------------------------------------
+
+subroutine get_reg_box_indices(lon, lat, x_ind, y_ind)
+
+! Given a longitude and latitude in degrees returns the index of the regular
+! lon-lat box that contains the point.
+
+real(r8), intent(in) :: lon, lat
+integer, intent(out) :: x_ind, y_ind
+
+call get_reg_lon_box(lon, x_ind)
+call get_reg_lat_box(lat, y_ind)
+
+end subroutine get_reg_box_indices
+
+!------------------------------------------------------------
+
+subroutine get_reg_lon_box(lon, x_ind)
+
+! Determine which regular longitude box a longitude is in.
+
+real(r8), intent(in) :: lon
+integer, intent(out) :: x_ind
+
+x_ind = int(num_reg_x * lon / 360.0_r8) + 1
+
+! Watch out for exactly at top; assume all lats and lons in legal range
+if(lon == 360.0_r8) x_ind = num_reg_x
+
+end subroutine get_reg_lon_box
+
+!------------------------------------------------------------
+
+subroutine get_reg_lat_box(lat, y_ind)
+
+! Determine which regular latitude box a latitude is in.
+
+real(r8), intent(in) :: lat
+integer, intent(out) :: y_ind
+
+y_ind = int(num_reg_y * (lat + 90.0_r8) / 180.0_r8) + 1
+
+! Watch out for exactly at top; assume all lats and lons in legal range
+if(lat == 90.0_r8) y_ind = num_reg_y
+
+end subroutine get_reg_lat_box
+
+!------------------------------------------------------------
+
+subroutine reg_box_overlap(x_corners, y_corners, is_pole, reg_lon_ind, reg_lat_ind)
+
+real(r8), intent(in) :: x_corners(4), y_corners(4)
+logical, intent(in) :: is_pole
+integer, intent(out) :: reg_lon_ind(2), reg_lat_ind(2)
+
+! Find a set of regular lat lon boxes that covers all of the area covered
+! by a dipole grid qaud whose corners are given by the
+! dimension four x_cornersand y_corners arrays. The two dimensional arrays reg_lon_ind and
+! reg_lat_ind return the first and last indices of the regular boxes in
+! latitude and longitude respectively. These indices may wraparound for
+! reg_lon_ind. A special computation is needed for a dipole quad that has
+! the true north pole in its interior. The logical is_pole is set to true
+! if this is the case. This can only happen for the t grid.
+! If the longitude boxes overlap 0 degrees, the indices returned are adjusted
+! by adding the total number of boxes to the second index (e.g. the indices
+! might be 88 and 93 for a case with 90 longitude boxes).
+
+real(r8) :: lat_min, lat_max, lon_min, lon_max
+integer :: i
+
+! A quad containing the pole is fundamentally different
+if(is_pole) then
+ ! Need all longitude boxes
+ reg_lon_ind(1) = 1
+ reg_lon_ind(2) = num_reg_x
+ ! Need to cover from lowest latitude to top box
+ lat_min = minval(y_corners)
+ reg_lat_ind(1) = int(num_reg_y * (lat_min + 90.0_r8) / 180.0_r8) + 1
+ call get_reg_lat_box(lat_min, reg_lat_ind(1))
+ reg_lat_ind(2) = num_reg_y
+else
+ ! All other quads do not contain the pole (pole could be on edge but no problem)
+ ! This is specific to the dipole POP grids that do not go to the south pole
+ ! Finding the range of latitudes is cake
+ lat_min = minval(y_corners)
+ lat_max = maxval(y_corners)
+
+ ! Figure out the indices of the regular boxes for min and max lats
+ call get_reg_lat_box(lat_min, reg_lat_ind(1))
+ call get_reg_lat_box(lat_max, reg_lat_ind(2))
+
+ ! Lons are much trickier. Need to make sure to wraparound the
+ ! right way. There is no guarantee on direction of lons in the
+ ! high latitude dipole rows.
+ ! All longitudes for non-pole rows have to be within 180 degrees
+ ! of one another.
+ lon_min = minval(x_corners)
+ lon_max = maxval(x_corners)
+ if((lon_max - lon_min) > 180.0_r8) then
+ ! If the max longitude value is more than 180
+ ! degrees larger than the min, then there must be wraparound.
+ ! Then, find the smallest value > 180 and the largest < 180 to get range.
+ lon_min = 360.0_r8
+ lon_max = 0.0_r8
+ do i=1, 4
+ if(x_corners(i) > 180.0_r8 .and. x_corners(i) < lon_min) lon_min = x_corners(i)
+ if(x_corners(i) < 180.0_r8 .and. x_corners(i) > lon_max) lon_max = x_corners(i)
+ end do
+ end if
+
+ ! Get the indices for the extreme longitudes
+ call get_reg_lon_box(lon_min, reg_lon_ind(1))
+ call get_reg_lon_box(lon_max, reg_lon_ind(2))
+
+ ! Watch for wraparound again; make sure that second index is greater than first
+ if(reg_lon_ind(2) < reg_lon_ind(1)) reg_lon_ind(2) = reg_lon_ind(2) + num_reg_x
+endif
+
+end subroutine reg_box_overlap
+
+!------------------------------------------------------------
+
+subroutine get_quad_corners(x, i, j, corners)
+
+real(r8), intent(in) :: x(:, :)
+integer, intent(in) :: i, j
+real(r8), intent(out) :: corners(4)
+
+! Grabs the corners for a given quadrilateral from the global array of lower
+! right corners. Note that corners go counterclockwise around the quad.
+
+integer :: ip1
+
+! Have to worry about wrapping in longitude but not in latitude
+ip1 = i + 1
+if(ip1 > nx) ip1 = 1
+
+corners(1) = x(i, j )
+corners(2) = x(ip1, j )
+corners(3) = x(ip1, j+1)
+corners(4) = x(i, j+1)
+
+end subroutine get_quad_corners
+
+!------------------------------------------------------------
+
+subroutine update_reg_list(reg_list_num, reg_list_lon, reg_list_lat, &
+ reg_lon_ind, reg_lat_ind, dipole_lon_index, dipole_lat_index)
+
+integer, intent(inout) :: reg_list_num(:, :), reg_list_lon(:, :, :), reg_list_lat(:, :, :)
+integer, intent(inout) :: reg_lon_ind(2), reg_lat_ind(2)
+integer, intent(in) :: dipole_lon_index, dipole_lat_index
+
+! Updates the data structure listing dipole quads that are in a given regular box
+integer :: ind_x, index_x, ind_y
+
+! Loop through indices for each possible regular cell
+! Have to watch for wraparound in longitude
+if(reg_lon_ind(2) < reg_lon_ind(1)) reg_lon_ind(2) = reg_lon_ind(2) + num_reg_x
+
+do ind_x = reg_lon_ind(1), reg_lon_ind(2)
+ ! Inside loop, need to go back to wraparound indices to find right box
+ index_x = ind_x
+ if(index_x > num_reg_x) index_x = index_x - num_reg_x
+
+ do ind_y = reg_lat_ind(1), reg_lat_ind(2)
+ ! Make sure the list storage isn't full
+ if(reg_list_num(index_x, ind_y) > max_reg_list_num) then
+ msgstring = 'max_reg_list_num is too small'
+ call error_handler(E_ERR, 'update_reg_list', msgstring, source, revision, revdate)
+ endif
+
+ ! Increment the count
+ reg_list_num(index_x, ind_y) = reg_list_num(index_x, ind_y) + 1
+ ! Store this quad in the list for this regular box
+ reg_list_lon(index_x, ind_y, reg_list_num(index_x, ind_y)) = dipole_lon_index
+ reg_list_lat(index_x, ind_y, reg_list_num(index_x, ind_y)) = dipole_lat_index
+ enddo
+enddo
+
+end subroutine update_reg_list
+
+!------------------------------------------------------------
+
+
+
+
+
+
subroutine init_conditions(x)
!------------------------------------------------------------------
!
@@ -513,7 +904,7 @@
! For Sea Surface Height don't need the vertical coordinate
if( vert_is_surface(location) ) then
- call lat_lon_interpolate(x(base_offset:), llon, llat, obs_type, 1, interp_val, istatus)
+ call lon_lat_interpolate(x(base_offset:), llon, llat, obs_type, 1, interp_val, istatus)
return
endif
@@ -530,7 +921,7 @@
offset = base_offset + (hgt_bot - 1) * nx * ny
if (debug > 1) print *, 'relative bot height offset = ', offset - base_offset
if (debug > 1) print *, 'absolute bot height offset = ', offset
-call lat_lon_interpolate(x(offset:), llon, llat, obs_type, hgt_bot, bot_val, istatus)
+call lon_lat_interpolate(x(offset:), llon, llat, obs_type, hgt_bot, bot_val, istatus)
! Failed istatus from interpolate means give up
if(istatus /= 0) return
@@ -539,7 +930,7 @@
offset = base_offset + (hgt_top - 1) * nx * ny
if (debug > 1) print *, 'relative top height offset = ', offset - base_offset
if (debug > 1) print *, 'absolute top height offset = ', offset
-call lat_lon_interpolate(x(offset:), llon, llat, obs_type, hgt_top, top_val, istatus)
+call lon_lat_interpolate(x(offset:), llon, llat, obs_type, hgt_top, top_val, istatus)
! Failed istatus from interpolate means give up
if(istatus /= 0) return
@@ -554,302 +945,349 @@
end subroutine model_interpolate
+! Three different types of grids are used here. The POP dipole
+! grid is referred to as a dipole grid and each region is
+! referred to as a quad, short for quadrilateral.
+! The longitude latitude rectangular grid with possibly irregular
+! spacing in latitude used for some POP applications and testing
+! is referred to as the irregular grid and each region is
+! called a box.
+! Finally, a regularly spaced longitude latitude grid is used
+! as a computational tool for interpolating from the dipole
+! grid. This is referred to as the regular grid and each region
+! is called a box.
+! All grids are referenced by the index of the lower left corner
+! of the quad or box.
-subroutine height_bounds(lheight, nheights, hgt_array, bot, top, fract, istatus)
+! The dipole grid is assumed to be global for all applications.
+! The irregular grid is also assumed to be global east
+! west for all applications.
+
+
+subroutine lon_lat_interpolate(x, lon, lat, var_type, height, interp_val, istatus)
!=======================================================================
!
-real(r8), intent(in) :: lheight
-integer, intent(in) :: nheights
-real(r8), intent(in) :: hgt_array(nheights)
-integer, intent(out) :: bot, top
-real(r8), intent(out) :: fract
-integer, intent(out) :: istatus
-! Local variables
-integer :: i
+! Subroutine to interpolate to a lon lat location given the state vector
+! for that level, x. This works just on one horizontal slice.
+! NOTE: Using array sections to pass in the x array may be inefficient on some
+! compiler/platform setups. Might want to pass in the entire array with a base
+! offset value instead of the section if this is an issue.
+! This routine works for either the dipole or a regular lat-lon grid.
+! Successful interpolation returns istatus=0.
+real(r8), intent(in) :: x(:)
+real(r8), intent(in) :: lon, lat
+integer, intent(in) :: var_type, height
+real(r8), intent(out) :: interp_val
+integer, intent(out) :: istatus
+
+! Local storage
+integer :: lat_bot, lat_top, lon_bot, lon_top, num_inds, start_ind
+integer :: x_ind, y_ind, i
+real(r8) :: p(4), x_corners(4), y_corners(4), xbot, xtop
+real(r8) :: lon_fract, lat_fract
+logical :: masked
+
if ( .not. module_initialized ) call static_init_model
-! Succesful istatus is 0
-! Make any failure here return istatus in the 20s
+! Succesful return has istatus of 0
istatus = 0
-! The zc array contains the depths of the center of the vertical grid boxes
-! In this case (unlike how we handle the MIT depths), positive is really down.
-! FIXME: in the MIT model, we're given box widths and we compute the centers,
-! and we computed them with larger negative numbers being deeper. Here,
-! larger positive numbers are deeper.
+! Get the lower left corner for either grid type
+if(dipole_grid) then
+ ! Figure out which of the regular grid boxes this is in
+ call get_reg_box_indices(lon, lat, x_ind, y_ind)
-! It is assumed that the top box is shallow and any observations shallower
-! than the depth of this box's center are just given the value of the
-! top box.
-if(lheight <= hgt_array(1)) then
- top = 1
- bot = 2
- ! NOTE: the fract definition is the relative distance from bottom to top
- fract = 1.0_r8
-if (debug > 1) print *, 'above first level in height'
-if (debug > 1) print *, 'hgt_array, top, bot, fract=', hgt_array(1), top, bot, fract
-endif
+ ! Is this on the U or T grid?
+ if(is_on_ugrid(var_type)) then
+ ! On U grid
+ num_inds = u_dipole_num (x_ind, y_ind)
+ start_ind = u_dipole_start(x_ind, y_ind)
-! Search through the boxes
-do i = 2, nheights
- ! If the location is shallower than this entry, it must be in this box
- if(lheight < hgt_array(i)) then
- top = i -1
- bot = i
- fract = (hgt_array(bot) - lheight) / (hgt_array(bot) - hgt_array(top))
-if (debug > 1) print *, 'i, hgt_array, top, bot, fract=', i, hgt_array(i), top, bot, fract
- return
- endif
-end do
+ ! If there are no quads overlapping, can't do interpolation
+ if(num_inds == 0) then
+ istatus = 1
+ return
+ endif
-! Falling off the end means the location is lower than the deepest height
-istatus = 20
+ ! Search the list of quads to see if (lon, lat) is in one
+ call get_dipole_quad(lon, lat, ulon, ulat, num_inds, start_ind, &
+ u_dipole_lon_list, u_dipole_lat_list, lon_bot, lat_bot, istatus)
+ ! Fail on bad istatus return
+ if(istatus /= 0) return
-end subroutine height_bounds
+ ! Getting corners for accurate interpolation
+ call get_quad_corners(ulon, lon_bot, lat_bot, x_corners)
+ call get_quad_corners(ulat, lon_bot, lat_bot, y_corners)
+ ! Fail if point is in one of the U boxes that go through the
+ ! pole (this could be fixed up if necessary)
+ if(lat_bot == u_pole_y .and. (lon_bot == pole_x -1 .or. &
+ lon_bot == pole_x)) then
+ istatus = 4
+ return
+ endif
-subroutine lat_lon_interpolate(x, llon, llat, var_type, hgt_bot, &
- interp_val, istatus)
-!=======================================================================
-!
+ else
+ ! On T grid
+ num_inds = t_dipole_num (x_ind, y_ind)
+ start_ind = t_dipole_start(x_ind, y_ind)
+ call get_dipole_quad(lon, lat, tlon, tlat, num_inds, start_ind, &
+ t_dipole_lon_list, t_dipole_lat_list, lon_bot, lat_bot, istatus)
+ ! Fail on bad istatus return
+ if(istatus /= 0) return
-! Subroutine to interpolate to a lat lon location given that portion of state
-! vector.
-! NOTE: Using array sections to pass in the x array may be inefficient on some
-! compiler/platform setups. Might want to pass in the entire array with a base
-! offset value instead of the section if this is an issue.
+ ! Fail if point is in T box that covers pole
+ if(lon_bot == pole_x .and. lat_bot == t_pole_y) then
+ istatus = 5
+ return
+ endif
-real(r8), intent(in) :: x(:)
-real(r8), intent(in) :: llon, llat
-integer, intent(in) :: var_type
-integer, intent(in) :: hgt_bot
-integer, intent(out) :: istatus
-real(r8), intent(out) :: interp_val
+ ! Getting corners for accurate interpolation
+ call get_quad_corners(tlon, lon_bot, lat_bot, x_corners)
+ call get_quad_corners(tlat, lon_bot, lat_bot, y_corners)
+ endif
-! Local storage
-real(r8) :: lat_array(ny), lon_array(nx)
-integer :: lat_bot, lat_top, lon_bot, lon_top
-real(r8) :: lat_fract, lon_fract
-real(r8) :: pa, pb, pc, pd, xbot, xtop
-integer :: lat_status, lon_status
+else
+ ! This is an irregular grid
+ ! U and V are on velocity grid
+ if (is_on_ugrid(var_type)) then
+ ! Get the corner indices and the fraction of the distance between
+ call get_irreg_box(lon, lat, ulon, ulat, &
+ lon_bot, lat_bot, lon_fract, lat_fract, istatus)
+ else
+ ! Eta, T and S are on the T grid
+ ! Get the corner indices
+ call get_irreg_box(lon, lat, tlon, tlat, &
+ lon_bot, lat_bot, lon_fract, lat_fract, istatus)
+ endif
-if ( .not. module_initialized ) call static_init_model
+ ! Return passing through error status
+ if(istatus /= 0) return
-! Succesful return has istatus of 0
-! Make any failure here return istatus in the 30s
-istatus = 30 ! unknown error
+endif
-! Find out what latitude box and fraction
-! The latitude grid being used depends on the variable type
-if(is_on_ugrid(var_type)) then
- ! U and V are on the YG latitude grid
- lat_array = yg
- call lat_bounds(llat, ny, lat_array, lat_bot, lat_top, lat_fract, lat_status)
-else
- ! SHGT, T and S are on the YC latitude grid
- lat_array = yc
- call lat_bounds(llat, ny, lat_array, lat_bot, lat_top, lat_fract, lat_status)
+! Find the indices to get the values for interpolating
+lat_top = lat_bot + 1
+if(lat_top > ny) then
+ istatus = 2
+ return
endif
-if (debug > 1) print *, 'lat bot, top, fract = ', lat_bot, lat_top, lat_fract
-! Check for error on the latitude interpolation
-if(lat_status /= 0) then
- istatus = 31
+! Watch for wraparound in longitude
+lon_top = lon_bot + 1
+if(lon_top > nx) lon_top = 1
+
+! Get the values at the four corners of the box or quad
+! Corners go around counterclockwise from lower left
+p(1) = get_val(lon_bot, lat_bot, nx, x, var_type, height, masked)
+if(masked) then
+ istatus = 3
return
endif
-! Find out what longitude box and fraction
-if(is_on_ugrid(var_type)) then
- ! U and V velocity is on the XG grid
- lon_array = xg
- call lon_bounds(llon, nx, lon_array, lon_bot, lon_top, lon_fract, lon_status)
-else
- ! SHGT, T, and S are on the XC grid
- lon_array = xc
- call lon_bounds(llon, nx, lon_array, lon_bot, lon_top, lon_fract, lon_status)
+p(2) = get_val(lon_top, lat_bot, nx, x, var_type, height, masked)
+if(masked) then
+ istatus = 3
+ return
endif
-if (debug > 1) print *, 'lon bot, top, fract = ', lon_bot, lon_top, lon_fract
-! Check for error on the longitude interpolation
-if(lat_status /= 0) then
- istatus = 32
+p(3) = get_val(lon_top, lat_top, nx, x, var_type, height, masked)
+if(masked) then
+ istatus = 3
return
endif
-! FIXME: is this the right place to test? make sure bottom of cell is valid
-! DO WE NEED TO TEST ALL 4 lon_bot, lon_top, lat_bot, lat_tops?
-if (is_not_ocean(var_type, lon_bot, lat_bot, hgt_bot) .or. &
- is_not_ocean(var_type, lon_top, lat_bot, hgt_bot) .or. &
- is_not_ocean(var_type, lon_bot, lat_top, hgt_bot) .or. &
- is_not_ocean(var_type, lon_top, lat_top, hgt_bot)) then
- istatus = 33
-if (debug > 1) print *, 'one corner not ocean'
+p(4) = get_val(lon_bot, lat_top, nx, x, var_type, height, masked)
+if(masked) then
+ istatus = 3
return
endif
-! Vector is laid out with lat outermost loop, lon innermost loop
-! Find the bounding points for the lat lon box
-if (debug > 1) print *, 'lon_bot, lon_top = ', lon_bot, lon_top
-if (debug > 1) print *, 'lat_bot, lat_top = ', lat_bot, lat_top
+! Full bilinear interpolation for quads
+if(dipole_grid) then
+ call quad_bilinear_interp(lon, lat, x_corners, y_corners, p, interp_val)
+else
+ ! Rectangular biliear interpolation
+ xbot = p(1) + lon_fract * (p(2) - p(1))
+ xtop = p(4) + lon_fract * (p(3) - p(4))
+ ! Now interpolate in latitude
+ interp_val = xbot + lat_fract * (xtop - xbot)
+endif
-pa = get_val(lon_bot, lat_bot, nx, x)
-if (debug > 1) print *, 'pa = ', pa
-pb = get_val(lon_top, lat_bot, nx, x)
-if (debug > 1) print *, 'pb = ', pb
-pc = get_val(lon_bot, lat_top, nx, x)
-if (debug > 1) print *, 'pc = ', pc
-pd = get_val(lon_top, lat_top, nx, x)
-if (debug > 1) print *, 'pd = ', pd
+end subroutine lon_lat_interpolate
-if (debug > 1) print *, 'pa,b,c,d = ', pa, pb, pc, pd
+!------------------------------------------------------------
-! Finish bi-linear interpolation
-! First interpolate in longitude
-if (debug > 1) print *, 'bot lon_fract, delta = ', lon_fract, (pb-pa)
-xbot = pa + lon_fract * (pb - pa)
-if (debug > 1) print *, 'xbot = ', xbot
-if (debug > 1) print *, 'top lon_fract, delta = ', lon_fract, (pd-pc)
-xtop = pc + lon_fract * (pd - pc)
-if (debug > 1) print *, 'xtop = ', xtop
-! Now interpolate in latitude
-if (debug > 1) print *, 'lat_fract, delta = ', lat_fract, (xtop - xbot)
-interp_val = xbot + lat_fract * (xtop - xbot)
-if (debug > 1) print *, 'lat_lon_interp: interp_val = ', interp_val
-! good return
+function get_val(lon_index, lat_index, nlon, x, var_type, height, masked)
+!=======================================================================
+!
+! Returns the value from a single level array given the lat and lon indices
+! 'masked' returns true if this is NOT a valid grid location (e.g. land, or
+! below the ocean floor in shallower areas).
+
+integer, intent(in) :: lon_index, lat_index, nlon, var_type, height
+real(r8), intent(in) :: x(:)
+logical, intent(out) :: masked
+real(r8) :: get_val
+
+if ( .not. module_initialized ) call static_init_model
+
+! check the land/ocean bottom map and return if not valid water cell.
+if(is_not_ocean(var_type, lon_index, lat_index, height)) then
+ masked = .true.
+ return
+endif
+
+! Layout has lons varying most rapidly
+get_val = x((lat_index - 1) * nlon + lon_index)
+
+! this is a valid ocean water cell, not land or below ocean floor
+masked = .false.
+
+end function get_val
+
+
+!------------------------------------------------------------
+
+subroutine get_irreg_box(lon, lat, lon_array, lat_array, &
+ found_x, found_y, lon_fract, lat_fract, istatus)
+!
+! Given a longitude and latitude of a point (lon and lat) and the
+! longitudes and latitudes of the lower left corner of the regular grid
+! boxes, gets the indices of the grid box that contains the point and
+! the fractions along each directrion for interpolation.
+
+real(r8), intent(in) :: lon, lat
+real(r8), intent(in) :: lon_array(nx, ny), lat_array(nx, ny)
+real(r8), intent(out) :: lon_fract, lat_fract
+integer, intent(out) :: found_x, found_y, istatus
+
+! Local storage
+integer :: lon_status, lat_status, lon_top, lat_top
+
+! Succesful return has istatus of 0
istatus = 0
-end subroutine lat_lon_interpolate
+! Get latitude box boundaries
+call lat_bounds(lat, ny, lat_array, found_y, lat_top, lat_fract, lat_status)
+! Check for error on the latitude interpolation
+if(lat_status /= 0) then
+ istatus = 1
+ return
+endif
+! Find out what longitude box and fraction
+call lon_bounds(lon, nx, lon_array, found_x, lon_top, lon_fract)
+end subroutine get_irreg_box
-subroutine lat_bounds(llat, nlats, lat_array, bot, top, fract, istatus)
+!------------------------------------------------------------
+
+subroutine lon_bounds(lon, nlons, lon_array, bot, top, fract)
+
!=======================================================================
!
-! Given a latitude llat, the array of latitudes for grid boundaries, and the
-! number of latitudes in the grid, returns the indices of the latitude
-! below and above the location latitude and the fraction of the distance
-! between. istatus is returned as 0 unless the location latitude is
-! south of the southernmost grid point (1 returned) or north of the
-! northernmost (2 returned). If one really had lots of polar obs would
-! want to worry about interpolating around poles.
+! Given a longitude lon, the array of longitudes for grid boundaries, and the
+! number of longitudes in the grid, returns the indices of the longitude
+! below and above the location longitude and the fraction of the distance
+! between. It is assumed that the longitude wraps around for a global grid.
+! Since longitude grids are going to be regularly spaced, this could be made more efficient.
+! Algorithm fails for a silly grid that has only two longitudes separated by 180 degrees.
-real(r8), intent(in) :: llat
-integer, intent(in) :: nlats
-real(r8), intent(in) :: lat_array(nlats)
+real(r8), intent(in) :: lon
+integer, intent(in) :: nlons
+real(r8), intent(in) :: lon_array(:, :)
integer, intent(out) :: bot, top
real(r8), intent(out) :: fract
-integer, intent(out) :: istatus
! Local storage
-integer :: i
+integer :: i
+real(r8) :: dist_bot, dist_top
if ( .not. module_initialized ) call static_init_model
-! Success should return 0, failure a positive number.
-! Make any failure here return istatus in the 40s
-istatus = 40
-
-! Check for too far south or north
-if(llat < lat_array(1)) then
- istatus = 41
- return
-else if(llat > lat_array(nlats)) then
- istatus = 42
- return
-endif
-
-! In the middle, search through
-do i = 2, nlats
- if(llat <= lat_array(i)) then
+! This is inefficient, someone could clean it up since longitudes are regularly spaced
+! But note that they don't have to start at 0
+do i = 2, nlons
+ dist_bot = lon_dist(lon, lon_array(i - 1, 1))
+ dist_top = lon_dist(lon, lon_array(i, 1))
+ if(dist_bot <= 0 .and. dist_top > 0) then
bot = i - 1
top = i
- fract = (llat - lat_array(bot)) / (lat_array(top) - lat_array(bot))
-if (debug > 1) print *, 'bot, top vals = ', lat_array(i-1), lat_array(i)
-if (debug > 1) print *, 'bot, top index = ', bot, top
-if (debug > 1) print *, 'fract = ', fract
- istatus = 0
+ fract = abs(dist_bot) / (abs(dist_bot) + dist_top)
return
endif
end do
-! shouldn't get here.
+! Falling off the end means it's in between; wraparound
+bot = nlons
+top = 1
+dist_bot = lon_dist(lon, lon_array(bot, 1))
+dist_top = lon_dist(lon, lon_array(top, 1))
+fract = abs(dist_bot) / (abs(dist_bot) + dist_top)
-end subroutine lat_bounds
+end subroutine lon_bounds
+!-------------------------------------------------------------
-subroutine lon_bounds(llon, nlons, lon_array, bot, top, fract, istatus)
+subroutine lat_bounds(lat, nlats, lat_array, bot, top, fract, istatus)
!=======================================================================
!
-! Given a longitude llon, the array of longitudes for grid boundaries, and the
-! number of longitudes in the grid, returns the indices of the longitude
-! below and above the location longitude and the fraction of the distance
-! between. istatus is returned as 0 unless the location longitude is
-! not between any of the longitude box boundaries. A module storage flag
-! indicates whether the grid wraps around in longitude.
-! Algorithm fails for a silly grid that
-! has only two longitudes separated by 180 degrees.
+! Given a latitude lat, the array of latitudes for grid boundaries, and the
+! number of latitudes in the grid, returns the indices of the latitude
+! below and above the location latitude and the fraction of the distance
+! between. istatus is returned as 0 unless the location latitude is
+! south of the southernmost grid point (1 returned) or north of the
+! northernmost (2 returned). If one really had lots of polar obs would
+! want to worry about interpolating around poles.
-real(r8), intent(in) :: llon
-integer, intent(in) :: nlons
-real(r8), intent(in) :: lon_array(nlons)
+real(r8), intent(in) :: lat
+integer, intent(in) :: nlats
+real(r8), intent(in) :: lat_array(:, :)
integer, intent(out) :: bot, top
real(r8), intent(out) :: fract
integer, intent(out) :: istatus
! Local storage
-integer :: i
-real(r8) :: dist_bot, dist_top
+integer :: i
if ( .not. module_initialized ) call static_init_model
! Success should return 0, failure a positive number.
-! Make any failure here return istatus in the 50s
-istatus = 50
+istatus = 0
-if (debug > 1) print *, 'computing bounds for = ', llon
-! This is inefficient, someone could clean it up for regularly spaced longitudes
-do i = 2, nlons
- dist_bot = lon_dist(llon, lon_array(i - 1))
- dist_top = lon_dist(llon, lon_array(i))
-if (debug > 3) print *, 'dist top, bot = ', dist_top, dist_bot
- if(dist_bot <= 0 .and. dist_top > 0) then
+! Check for too far south or north
+if(lat < lat_array(1, 1)) then
+ istatus = 1
+ return
+else if(lat > lat_array(1, nlats)) then
+ istatus = 2
+ return
+endif
+
+! In the middle, search through
+do i = 2, nlats
+ if(lat <= lat_array(1, i)) then
bot = i - 1
top = i
-if (debug > 1) print *, 'bot, top vals = ', lon_array(i-1), lon_array(i)
-if (debug > 1) print *, 'bot, top index = ', bot, top
-if (debug > 1) print *, 'numerator = ', dist_bot
-if (debug > 1) print *, 'denominator = ', dist_bot + abs(dist_top)
- fract = abs(dist_bot) / (abs(dist_bot) + dist_top)
-if (debug > 1) print *, 'fract = ', fract
- istatus = 0
+ fract = (lat - lat_array(1, bot)) / (lat_array(1, top) - lat_array(1, bot))
return
endif
end do
-! Falling off the end means it's in between. Check for wraparound.
-if(longitude_wrap) then
- bot = nlons
- top = 1
- dist_bot = lon_dist(llon, lon_array(bot))
- dist_top = lon_dist(llon, lon_array(top))
- fract = abs(dist_bot) / (abs(dist_bot) + dist_top)
- istatus = 0
- return
-else
- ! Not in the localized longitude grid; return istatus 1
- istatus = 51
-endif
+! Shouldn't get here. Might want to fail really hard through error handler
+istatus = 40
-end subroutine lon_bounds
+end subroutine lat_bounds
@@ -878,27 +1316,489 @@
end function lon_dist
-function get_val(lon_index, lat_index, nlon, x)
+!------------------------------------------------------------
+
+
+subroutine get_dipole_quad(lon, lat, qlons, qlats, num_inds, start_ind, &
+ x_inds, y_inds, found_x, found_y, istatus)
+
+real(r8), intent(in) :: lon, lat, qlons(:, :), qlats(:, :)
+integer, intent(in) :: num_inds, start_ind, x_inds(:), y_inds(:)
+integer, intent(out) :: found_x, found_y, istatus
+
+! Given the lon and lat of a point, and a list of the
+! indices of the quads that might contain a point at (lon, lat), determines
+! which quad contains the point. istatus is returned as 0 if all went
+! well and 1 if the point was not found to be in any of the quads.
+
+integer :: i, my_index
+real(r8) :: x_corners(4), y_corners(4)
+
+istatus = 0
+
+! Loop through all the quads and see if the point is inside
+do i = 1, num_inds
+ my_index = start_ind + i - 1
+ call get_quad_corners(qlons, x_inds(my_index), y_inds(my_index), x_corners)
+ call get_quad_corners(qlats, x_inds(my_index), y_inds(my_index), y_corners)
+
+ ! Ssearch in this individual quad
+ if(in_quad(lon, lat, x_corners, y_corners)) then
+ found_x = x_inds(my_index)
+ found_y = y_inds(my_index)
+ return
+ endif
+enddo
+
+! Falling off the end means search failed, return istatus 1
+istatus = 1
+
+end subroutine get_dipole_quad
+
+
+!------------------------------------------------------------
+function in_quad(lon, lat, x_corners, y_corners)
+
+real(r8), intent(in) :: lon, lat, x_corners(4), y_corners(4)
+logical :: in_quad
+
+! Return in_quad true if the point (lon, lat) is in the quad with
+! the given corners.
+
+! Do this by ray tracing in latitude for now. For non-pole point, want a
+! downward directed ray from (lon,lat) to intersect only one edge of the
+! quadrilateral. Have to watch out for intersecting at a vertex nearly
+! exactly, other round-off issues. This should be checked carefully.
+
+real(r8) :: x(2), y(2)
+logical :: cant_be_in_box, on_side
+integer :: intercepts(4), exact_corner(4), second_corner, num_sides, num_corners, i
+
+! Default answer is point is not in quad
+in_quad = .false.
+
+! Loop through the sides and compute intercept (if any) with a downward
+! ray from the point. This ray has equation x=lon, y<=lat.
+do i = 1, 4
+ ! Load up the sides endpoints
+ if(i <= 3) then
+ x(1:2) = x_corners(i:i+1)
+ y(1:2) = y_corners(i:i+1)
+ else
+ x(1) = x_corners(4)
+ x(2) = x_corners(1)
+ y(1) = y_corners(4)
+ y(2) = y_corners(1)
+ endif
+
+ ! Check to see how a southward going ray from the point is related to this side
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list