diag_grid.F90

Go to the documentation of this file.

!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************

MODULE diag_grid_mod
#include <fms_platform.h>
  ! <CONTACT EMAIL="seth.underwood@noaa.gov">
  !   Seth Underwood
  ! </CONTACT>
  ! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/" />
  ! <OVERVIEW>
  !   <TT>diag_grid_mod</TT> is a set of procedures to work with the
  !   model's global grid to allow regional output.
  ! </OVERVIEW>
  ! <DESCRIPTION>
  !   <TT>diag_grid_mod</TT> contains useful utilities for dealing
  !   with, mostly, regional output for grids other than the standard
  !   lat/lon grid.  This module contains three public procedures <TT>
  !   diag_grid_init</TT>, which is shared globably in the <TT>
  !   diag_manager_mod</TT>, <TT>diag_grid_end</TT> which will free
  !   up memory used during the register field calls, and
  !   <TT>get_local_indexes</TT>.  The <TT>send_global_grid</TT>
  !   procedure is called by the model that creates the global grid.
  !   <TT>send_global_grid</TT> needs to be called before any fields
  !   are registered that will output only regions.  <TT>get_local_indexes</TT>
  !   is to be called by the <TT>diag_manager_mod</TT> to discover the
  !   global indexes defining a subregion on the tile.
  !
  !   <B>Change Log</B>
  !   <DL>
  !     <DT>September 2009</DT>
  !     <DD>
  !       <UL>
  !         <LI>Single point region in Cubed Sphere</LI>
  !         <LI>Single tile regions in the cubed sphere</LI>
  !       </UL>
  !     </DD>
  !   </DL>
  ! </DESCRIPTION>

  ! <INFO>
  !   <FUTURE>
  !     Multi-tile regional output in the cubed sphere.
  !   </FUTURE>
  !   <FUTURE>
  !     Single grid in the tri-polar grid.
  !   </FUTURE>
  !   <FUTURE>
  !     Multi-tile regional output in the tri-polar grid.
  !   </FUTURE>
  !   <FUTURE>
  !     Regional output using array masking.  This should allow
  !     regional output to work on any current or future grid.
  !   </FUTURE>
  ! </INFO>
  USE constants_mod, ONLY: deg_to_rad, rad_to_deg, radius
  USE fms_mod, ONLY: write_version_number, error_mesg, warning, fatal,&
       & mpp_pe
  USE mpp_mod, ONLY: mpp_root_pe, mpp_npes, mpp_max, mpp_min
  USE mpp_domains_mod, ONLY: domain2d, mpp_get_tile_id,&
       & mpp_get_ntile_count, mpp_get_compute_domains

  IMPLICIT NONE

  ! Parameters
  ! Include variable "version" to be written to log file.
#include<file_version.h>

  ! Derived data types
  ! <PRIVATE>
  ! <TYPE NAME="diag_global_grid_type">
  !   <DESCRIPTION>
  !     Contains the model's global grid data, and other grid information.
  !   </DESCRIPTION>
  !   <DATA NAME="glo_lat" TYPE="REAL, _ALLOCATABLE, DIMENSION(:,:)">
  !     The latitude values on the global grid.
  !   </DATA>
  !   <DATA NAME="glo_lon" TYPE="REAL, _ALLOCATABLE, DIMENSION(:,:)">
  !     The longitude values on the global grid.
  !   </DATA>
  !   <DATA NAME="aglo_lat" TYPE="REAL, _ALLOCATABLE, DIMENSION(:,:)">
  !     The latitude values on the global a-grid.  Here we expect isc-1:iec+1 and
  !     jsc=1:jec+1 to be passed in.
  !   </DATA>
  !   <DATA NAME="aglo_lon" TYPE="REAL, _ALLOCATABLE, DIMENSION(:,:)">
  !     The longitude values on the global a-grid.  Here we expec isc-1:iec+j and
  !     jsc-1:jec+1 to be passed in.
  !   </DATA>
  !   <DATA NAME="myXbegin" TYPE="INTEGER">
  !     The starting index of the compute domain on the current PE.
  !   </DATA>
  !   <DATA NAME="myYbegin" TYPE="INTEGER">
  !     The starting index of the compute domain on the cureent PE.
  !   </DATA>
  !   <DATA NAME="dimI" TYPE="INTEGER">
  !     The dimension of the global grid in the 'i' / longitudal direction.
  !   </DATA>
  !   <DATA NAME="dimJ" TYPE="INTEGER">
  !     The dimension of the global grid in the 'j' / latitudal direction.
  !   </DATA>
  !   <DATA NAME="adimI" TYPE="INTEGER">
  !     The dimension of the global a-grid in the 'i' / longitudal direction.  Again,
  !     the expected dimension for diag_grid_mod is isc-1:iec+1.
  !   </DATA>
  !   <DATA NAME="adimJ" TYPE="INTEGER">
  !     The dimension of the global a-grid in the 'j' / latitudal direction.  Again,
  !     the expected dimension for diag_grid_mod is jsc-1:jec+1.
  !   </DATA>
  !   <DATA NAME="tile_number" TYPE="INTEGER">
  !     The tile the <TT>glo_lat</TT> and <TT>glo_lon</TT> define.
  !   </DATA>
  !   <DATA NAME="ntimes" TYPE="INTEGER">
  !     The number of tiles.
  !   </DATA>
  !   <DATA NAME="peStart" TYPE="INTEGER">
  !     The starting PE number for the current tile.
  !   </DATA>
  !   <DATA NAME="peEnd" TYPE="INTEGER">
  !     The ending PE number for the current tile.
  !   </DATA>
  !   <DATA NAME="grid_type" TYPE="CHARACTER(len=128)">
  !     The global grid type.
  !   </DATA>
  TYPE :: diag_global_grid_type
     REAL, _allocatable, DIMENSION(:,:) :: glo_lat, glo_lon
     REAL, _allocatable, DIMENSION(:,:) :: aglo_lat, aglo_lon
     INTEGER :: myxbegin, myybegin
     INTEGER :: dimi, dimj
     INTEGER :: adimi, adimj
     INTEGER :: tile_number
     INTEGER :: ntiles
     INTEGER :: pestart, peend
     CHARACTER(len=128) :: grid_type
  END TYPE diag_global_grid_type
  ! </TYPE>
  ! </PRIVATE>

  ! <PRIVATE>
  ! <TYPE NAME="point">
  !   <DESCRIPTION>
  !      Private point type to hold the (x,y,z) location for a (lat,lon)
  !      location.
  !   </DESCRIPTION>
  !   <DATA NAME="x" TYPE="REAL">
  !     The x value of the (x,y,z) coordinates.
  !   </DATA>
  !   <DATA NAME="y" TYPE="REAL">
  !     The y value of the (x,y,z) coordinates.
  !   </DATA>
  !   <DATA NAME="z" TYPE="REAL">
  !     The z value of the (x,y,z) coordinates.
  !   </DATA>
  TYPE :: point
     REAL :: x,y,z
  END TYPE point
  ! </TYPE>
  ! </PRIVATE>

  ! <PRIVATE>
  ! <DATA NAME="diag_global_grid" TYPE="TYPE(diag_global_grid_type)">
  !   Variable to hold the global grid data
  ! </DATA>
  ! </PRIVATE>
  TYPE(diag_global_grid_type) :: diag_global_grid

  ! <PRIVATE>
  ! <DATA NAME="diag_grid_initialized" TYPE="LOGICAL" DEFAULT=".FALSE.">
  !   Indicates if the diag_grid_mod has been initialized.
  ! </DATA>
  ! </PRIVATE>
  LOGICAL :: diag_grid_initialized = .false.

  PRIVATE
  PUBLIC :: diag_grid_init, diag_grid_end, get_local_indexes,  &
            get_local_indexes2

CONTAINS

  ! <SUBROUTINE NAME="diag_grid_init">
  !   <OVERVIEW>
  !     Send the global grid to the <TT>diag_manager_mod</TT> for
  !     regional output.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE diag_grid_init(domain, glo_lat, glo_lon, aglo_lat, aglo_lon)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     In order for the diag_manager to do regional output for grids
  !     other than the standard lat/lon grid, the <TT>
  !     diag_manager_mod</TT> needs to know the the latitude and
  !     longitude values for the entire global grid.  This procedure
  !     is the mechanism the models will use to share their grid with
  !     the diagnostic manager.
  !
  !     This procedure needs to be called after the grid is created,
  !     and before the first call to register the fields.
  !   </DESCRIPTION>
  !   <IN NAME="domain" TYPE="INTEGER">
  !     The domain to which the grid data corresponds.
  !   </IN>
  !   <IN NAME="glo_lat" TYPE="REAL, DIMENSION(:,:)">
  !     The latitude information for the grid tile.
  !   </IN>
  !   <IN NAME="glo_lon" TYPE="REAL, DIMENSION(:,:)">
  !     The longitude information for the grid tile.
  !   </IN>
  !   <IN NAME="aglo_lat" TYPE="REAL, DIMENSION(:,:)">
  !     The latitude information for the a-grid tile.
  !   </IN>
  !   <IN NAME="aglo_lon" TYPE="REAL, DIMENSION(:,:)">
  !     The longitude information for the a-grid tile.
  !   </IN>
  SUBROUTINE diag_grid_init(domain, glo_lat, glo_lon, aglo_lat, aglo_lon)
    TYPE(domain2d), INTENT(in) :: domain
    REAL, INTENT(in), DIMENSION(:,:) :: glo_lat, glo_lon
    REAL, INTENT(in), DIMENSION(:,:) :: aglo_lat, aglo_lon

    INTEGER, DIMENSION(1) :: tile
    INTEGER :: ntiles
    INTEGER :: stat
    INTEGER :: i_dim, j_dim
    INTEGER :: ai_dim, aj_dim
    INTEGER, DIMENSION(2) :: latdim, londim
    INTEGER, DIMENSION(2) :: alatdim, alondim
    INTEGER :: mype, npes, npespertile
    INTEGER, ALLOCATABLE, DIMENSION(:) :: xbegin, xend, ybegin, yend

    ! Write the file version to the logfile
    CALL write_version_number("DIAG_GRID_MOD", version)

    ! Verify all allocatable / pointers for diag_global_grid hare not
    ! allocated / associated.
    IF ( ALLOCATED(xbegin) ) DEALLOCATE(xbegin)
    IF ( ALLOCATED(ybegin) ) DEALLOCATE(ybegin)
    IF ( ALLOCATED(xend) ) DEALLOCATE(xend)
    IF ( ALLOCATED(yend) ) DEALLOCATE(yend)

    ! What is my PE
    mype = mpp_pe() -mpp_root_pe() + 1

    ! Get the domain/pe layout, and allocate the [xy]begin|end arrays/pointers
    npes = mpp_npes()
    ALLOCATE(xbegin(npes), &
         &   ybegin(npes), &
         &   xend(npes), &
         &   yend(npes), stat=stat)
    IF ( stat .NE. 0 ) THEN
       CALL error_mesg('diag_grid_mod::diag_grid_init',&
            &'Could not allocate memory for the compute grid indices&
            &.', fatal)
    END IF

    ! Get tile information
    ntiles = mpp_get_ntile_count(domain)
    tile = mpp_get_tile_id(domain)

    ! Number of PEs per tile
    npespertile = npes / ntiles
    diag_global_grid%peEnd = npespertile * tile(1)
    diag_global_grid%peStart = diag_global_grid%peEnd - npespertile + 1

    ! Get the compute domains
    CALL mpp_get_compute_domains(domain,&
         & xbegin=xbegin, xend=xend,&
         & ybegin=ybegin, yend=yend)

    ! Module initialized
    diag_grid_initialized = .true.

    ! Get the size of the grids
    latdim = shape(glo_lat)
    londim = shape(glo_lon)
    IF (  (latdim(1) == londim(1)) .AND.&
         &(latdim(2) == londim(2)) ) THEN
       i_dim = latdim(1)
       j_dim = latdim(2)
    ELSE
       CALL error_mesg('diag_grid_mod::diag_grid_init',&
            &'glo_lat and glo_lon must be the same shape.', fatal)
    END IF

    ! Same thing for the a-grid
    alatdim = shape(aglo_lat)
    alondim = shape(aglo_lon)
    IF (  (alatdim(1) == alondim(1)) .AND. &
         &(alatdim(2) == alondim(2)) ) THEN
       IF ( tile(1) == 4 .OR. tile(1) == 5 ) THEN
          ! These tiles need to be transposed.
          ai_dim = alatdim(2)
          aj_dim = alatdim(1)
       ELSE
          ai_dim = alatdim(1)
          aj_dim = alatdim(2)
       END IF
    ELSE
       CALL error_mesg('diag_grid_mod::diag_grid_init',&
            & "a-grid's glo_lat and glo_lon must be the same shape.", fatal)
    END IF

    ! Allocate the grid arrays
    IF (   _allocated(diag_global_grid%glo_lat) .OR.&
         & _allocated(diag_global_grid%glo_lon) ) THEN
       IF ( mpp_pe() == mpp_root_pe() ) &
            & CALL error_mesg('diag_grid_mod::diag_grid_init',&
            &'The global grid has already been initialized', warning)
    ELSE
       ALLOCATE(diag_global_grid%glo_lat(i_dim,j_dim),&
            &   diag_global_grid%glo_lon(i_dim,j_dim), stat=stat)
       IF ( stat .NE. 0 ) THEN
          CALL error_mesg('diag_grid_mod::diag_grid_init',&
               &'Could not allocate memory for the global grid.', fatal)
       END IF
    END IF

    ! Same thing for the a-grid
    IF (   _allocated(diag_global_grid%aglo_lat) .OR.&
         & _allocated(diag_global_grid%aglo_lon) ) THEN
       IF ( mpp_pe() == mpp_root_pe() ) &
            & CALL error_mesg('diag_grid_mod::diag_grid_init',&
            &'The global a-grid has already been initialized', warning)
    ELSE
       ALLOCATE(diag_global_grid%aglo_lat(0:ai_dim-1,0:aj_dim-1),&
            &   diag_global_grid%aglo_lon(0:ai_dim-1,0:aj_dim-1), stat=stat)
       IF ( stat .NE. 0 ) THEN
          CALL error_mesg('diag_global_mod::diag_grid_init',&
               &'Could not allocate memory for the global a-grid', fatal)
       END IF
    END IF

    ! Set the values for diag_global_grid

    ! If we are on tile 4 or 5, we need to transpose the grid to get
    ! this to work.
    IF ( tile(1) == 4 .OR. tile(1) == 5 ) THEN
       diag_global_grid%aglo_lat = transpose(aglo_lat)
       diag_global_grid%aglo_lon = transpose(aglo_lon)
    ELSE
       diag_global_grid%aglo_lat = aglo_lat
       diag_global_grid%aglo_lon = aglo_lon
    END IF
    diag_global_grid%glo_lat = glo_lat
    diag_global_grid%glo_lon = glo_lon
    diag_global_grid%dimI = i_dim
    diag_global_grid%dimJ = j_dim
    diag_global_grid%adimI = ai_dim
    diag_global_grid%adimJ = aj_dim
    !--- For the nested model, the nested region only has 1 tile ( ntiles = 1) but
    !--- the tile_id is 7 for the nested region. In the routine get_local_indexes,
    !--- local variables ijMin and ijMax have dimesnion (ntiles) and will access
    !--- ijMin(diag_global_grid%tile_number,:). For the nested region, ntiles = 1 and
    !--- diag_global_grid%tile_number = 7 will cause out of bounds. So need to
    !--- set diag_global_grid%tile_number = 1 when ntiles = 1 for the nested model.
    if(ntiles == 1) then
       diag_global_grid%tile_number = 1
    else
       diag_global_grid%tile_number = tile(1)
    endif
    diag_global_grid%ntiles = ntiles
    diag_global_grid%myXbegin = xbegin(mype)
    diag_global_grid%myYbegin = ybegin(mype)

    ! Unallocate arrays used here
    DEALLOCATE(xbegin)
    DEALLOCATE(ybegin)
    DEALLOCATE(xend)
    DEALLOCATE(yend)
  END SUBROUTINE diag_grid_init
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="diag_grid_end">
  !   <OVERVIEW>
  !     Unallocate the diag_global_grid variable.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE diag_grid_end()
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     The <TT>diag_global_grid</TT> variable is only needed during
  !     the register field calls, and then only if there are fields
  !     requestion regional output.  Once all the register fields
  !     calls are complete (before the first <TT>send_data</TT> call
  !     this procedure can be called to free up memory.
  !   </DESCRIPTION>
  SUBROUTINE diag_grid_end()

    IF ( diag_grid_initialized ) THEN
       ! De-allocate grid
       IF ( _allocated(diag_global_grid%glo_lat) ) THEN
          DEALLOCATE(diag_global_grid%glo_lat)
       ELSE IF ( mpp_pe() == mpp_root_pe() ) THEN
          CALL error_mesg('diag_grid_mod::diag_grid_end',&
               &'diag_global_grid%glo_lat was not allocated.', warning)
       END IF

       IF ( _allocated(diag_global_grid%glo_lon) ) THEN
          DEALLOCATE(diag_global_grid%glo_lon)
       ELSE IF ( mpp_pe() == mpp_root_pe() ) THEN
          CALL error_mesg('diag_grid_mod::diag_grid_end',&
               &'diag_global_grid%glo_lon was not allocated.', warning)
       END IF
       ! De-allocate a-grid
       IF ( _allocated(diag_global_grid%aglo_lat) ) THEN
          DEALLOCATE(diag_global_grid%aglo_lat)
       ELSE IF ( mpp_pe() == mpp_root_pe() ) THEN
          CALL error_mesg('diag_grid_mod::diag_grid_end',&
               &'diag_global_grid%aglo_lat was not allocated.', warning)
       END IF

       IF ( _allocated(diag_global_grid%aglo_lon) ) THEN
          DEALLOCATE(diag_global_grid%aglo_lon)
       ELSE IF ( mpp_pe() == mpp_root_pe() ) THEN
          CALL error_mesg('diag_grid_mod::diag_grid_end',&
               &'diag_global_grid%aglo_lon was not allocated.', warning)
       END IF

       diag_grid_initialized = .false.
    END IF
  END SUBROUTINE diag_grid_end
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="get_local_indexes">
  !   <OVERVIEW>
  !     Find the local start and local end indexes on the local PE
  !     for regional output.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE get_local_indexes(latStart, latEnd, lonStart,
  !     lonEnd, istart, iend, jstart, jend)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Given a defined region, find the local indexes on the local
  !     PE surrounding the region.
  !   </DESCRIPTION>
  !   <IN NAME="latStart" TYPE="REAL">
  !     The minimum latitude value defining the region.  This value
  !     must be less than latEnd, and be in the range [-90,90]
  !   </IN>
  !   <IN NAME="latEnd" TYPE="REAL">
  !     The maximum latitude value defining the region.  This value
  !     must be greater than latStart, and be in the range [-90,90]
  !   </IN>
  !   <IN NAME="lonStart" TYPE="REAL">
  !     The western most longitude value defining the region.
  !     Possible ranges are either [-180,180] or [0,360].
  !   </IN>
  !   <IN NAME="lonEnd" TYPE="REAL">
  !     The eastern most longitude value defining the region.
  !     Possible ranges are either [-180,180] or [0,360].
  !   </IN>
  !   <OUT NAME="istart" TYPE="INTEGER">
  !     The local start index on the local PE in the 'i' direction.
  !   </OUT>
  !   <OUT NAME="iend" TYPE="INTEGER">
  !     The local end index on the local PE in the 'i' direction.
  !   </OUT>
  !   <OUT NAME="jstart" TYPE="INTEGER">
  !     The local start index on the local PE in the 'j' direction.
  !   </OUT>
  !   <OUT NAME="jend" TYPE="INTEGER">
  !     The local end index on the local PE in the 'j' direction.
  !   </OUT>
  SUBROUTINE get_local_indexes(latStart, latEnd, lonStart, lonEnd,&
       & istart, iend, jstart, jend)
    REAL, INTENT(in) :: latstart, lonstart !< lat/lon start angles
    REAL, INTENT(in) :: latend, lonend !< lat/lon end angles
    INTEGER, INTENT(out) :: istart, jstart !< i/j start indexes
    INTEGER, INTENT(out) :: iend, jend !< i/j end indexes

    REAL, ALLOCATABLE, DIMENSION(:,:) :: delta_lat, delta_lon, grid_lon

    REAL, DIMENSION(4) :: dists_lon, dists_lat
    REAL :: lonendadj, my_lonstart, my_lonend

    INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmin, ijmax
    INTEGER :: mytile, ntiles, i, j, k, dimi, dimj, istat
    INTEGER :: count

    LOGICAL :: onmype

    !For cfsite potential fix.
    INTEGER :: mini
    INTEGER :: minj
    REAL :: minimum_distance
    REAL :: global_min_distance
    INTEGER :: rank_buf

    IF ( .NOT. diag_grid_initialized )&
         & CALL error_mesg('diag_grid_mod::get_local_indexes',&
         &'Module not initialized, first initialze module with a call &
         &to diag_grid_init', fatal)

    ! Make adjustment for negative longitude values
    if ( lonstart < 0. ) then
       my_lonstart = lonstart + 360.
    else
       my_lonstart = lonstart
    end if
    if ( lonend < 0. ) then
       my_lonend = lonend + 360.
    else
       my_lonend = lonend
    end if

    IF (latstart .EQ. latend .AND. my_lonstart .EQ. my_lonend) THEN

        !For a single point, use the a-grid longitude and latitude
        !values.

        mytile = diag_global_grid%tile_number
        ntiles = diag_global_grid%ntiles

        allocate(ijmin(ntiles,2))
        ijmin = 0

        !Find the i,j indices of the a-grid point nearest to the
        !my_lonStart,latStart point.
        CALL find_nearest_agrid_index(latstart, &
                                      my_lonstart, &
                                      mini, &
                                      minj, &
                                      minimum_distance)

        !Find the minimum distance across all ranks.
        global_min_distance = minimum_distance
        CALL mpp_min(global_min_distance)

        !In the case of a tie (i.e. two ranks with exactly the same
        !minimum distance), use the i,j values from the larger rank id.
        IF (global_min_distance .EQ. minimum_distance) THEN
            rank_buf = mpp_pe()
        ELSE
            rank_buf = -1
        ENDIF
        CALL mpp_max(rank_buf)

        !Sanity check.
        IF (rank_buf .EQ. -1) THEN
            CALL error_mesg("get_local_indexes", &
                            "No rank has minimum distance.", &
                            fatal)
        ENDIF

        IF (rank_buf .EQ. mpp_pe()) THEN
            ijmin(mytile,1) = mini + diag_global_grid%myXbegin - 1
            ijmin(mytile,2) = minj + diag_global_grid%myYbegin - 1
        ENDIF

        DO i = 1,ntiles
            CALL mpp_max(ijmin(i,1))
            CALL mpp_max(ijmin(i,2))
        ENDDO

        istart = ijmin(mytile,1)
        jstart = ijmin(mytile,2)
        iend = istart
        jend = jstart

        DEALLOCATE(ijmin)
    ELSE

        mytile = diag_global_grid%tile_number
        ntiles = diag_global_grid%ntiles

        ! Arrays to home min/max for each tile
        ALLOCATE(ijmin(ntiles,2), stat=istat)
        IF ( istat .NE. 0 )&
             & CALL error_mesg('diag_grid_mod::get_local_indexes',&
             &'Cannot allocate ijMin index array', fatal)
        ALLOCATE(ijmax(ntiles,2), stat=istat)
        IF ( istat .NE. 0 )&
             & CALL error_mesg('diag_grid_mod::get_local_indexes',&
             &'Cannot allocate ijMax index array', fatal)
        ijmin = 0
        ijmax = 0

        ! There will be four points to define a region, find all four.
        ! Need to call the correct function depending on if the tile is a
        ! pole tile or not.
       dimi = diag_global_grid%dimI
       dimj = diag_global_grid%dimJ

       ! Build the delta array
       ALLOCATE(delta_lat(dimi,dimj), stat=istat)
       IF ( istat .NE. 0 )&
            & CALL error_mesg('diag_grid_mod::get_local_indexes',&
            &'Cannot allocate latitude delta array', fatal)
       ALLOCATE(delta_lon(dimi,dimj), stat=istat)
       IF ( istat .NE. 0 )&
            & CALL error_mesg('diag_grid_mod::get_local_indexes',&
            &'Cannot allocate longitude delta array', fatal)
       DO j=1, dimj
          DO i=1, dimi
             count = 0
             dists_lon = 0.
             dists_lat = 0.
             IF ( i < dimi ) THEN
                dists_lon(1) = abs(diag_global_grid%glo_lon(i+1,j) - diag_global_grid%glo_lon(i,j))
                dists_lat(1) = abs(diag_global_grid%glo_lat(i+1,j) - diag_global_grid%glo_lat(i,j))
                count = count+1
             END IF
             IF ( j < dimj ) THEN
                dists_lon(2) = abs(diag_global_grid%glo_lon(i,j+1) - diag_global_grid%glo_lon(i,j))
                dists_lat(2) = abs(diag_global_grid%glo_lat(i,j+1) - diag_global_grid%glo_lat(i,j))
                count = count+1
             END IF
             IF ( i > 1 ) THEN
                dists_lon(3) = abs(diag_global_grid%glo_lon(i,j) - diag_global_grid%glo_lon(i-1,j))
                dists_lat(3) = abs(diag_global_grid%glo_lat(i,j) - diag_global_grid%glo_lat(i-1,j))
                count = count+1
             END IF
             IF ( j > 1 ) THEN
                dists_lon(4) = abs(diag_global_grid%glo_lon(i,j) - diag_global_grid%glo_lon(i,j-1))
                dists_lat(4) = abs(diag_global_grid%glo_lat(i,j) - diag_global_grid%glo_lat(i,j-1))
                count = count+1
             END IF

             ! Fix wrap around problem
             DO k=1, 4
                IF ( dists_lon(k) > 180.0 ) THEN
                   dists_lon(k) = 360.0 - dists_lon(k)
                END IF
             END DO
             delta_lon(i,j) = sum(dists_lon)/real(count)
             delta_lat(i,j) = sum(dists_lat)/real(count)
          END DO
       END DO

       ijmin = huge(1)
       ijmax = -huge(1)

       ! Adjusted longitude array
       ALLOCATE(grid_lon(dimi,dimj), stat=istat)
       IF ( istat .NE. 0 )&
            & CALL error_mesg('diag_grid_mod::get_local_indexes',&
            &'Cannot allocate temporary longitude array', fatal)
       grid_lon = diag_global_grid%glo_lon

       ! Make adjustments where required
       IF ( my_lonstart > my_lonend ) THEN
          WHERE ( grid_lon < my_lonstart )
             grid_lon = grid_lon + 360.0
          END WHERE
          lonendadj = my_lonend + 360.0
       ELSE
          lonendadj = my_lonend
       END IF

       DO j=1, dimj-1
          DO i=1, dimi-1
             onmype = .false.
             IF ( latstart-delta_lat(i,j) <= diag_global_grid%glo_lat(i,j) .AND.&
                  & diag_global_grid%glo_lat(i,j) < latend+delta_lat(i,j) ) THEN
                ! Short-cut for the poles
                IF ( (abs(latstart)-delta_lat(i,j) <= 90.0 .AND.&
                     & 90.0 <= abs(latend)+delta_lat(i,j)) .AND.&
                     & abs(diag_global_grid%glo_lat(i,j)) == 90.0 ) THEN
                   onmype = .true.
                ELSE IF ( (my_lonstart-delta_lon(i,j) <= grid_lon(i,j) .AND.&
                     & grid_lon(i,j) < lonendadj+delta_lon(i,j)) ) THEN
                   onmype = .true.
                ELSE
                   onmype = .false.
                END IF
                IF ( onmype ) THEN
                   ijmin(mytile,1) = min(ijmin(mytile,1),i + diag_global_grid%myXbegin - 1)
                   ijmax(mytile,1) = max(ijmax(mytile,1),i + diag_global_grid%myXbegin - 1)
                   ijmin(mytile,2) = min(ijmin(mytile,2),j + diag_global_grid%myYbegin - 1)
                   ijmax(mytile,2) = max(ijmax(mytile,2),j + diag_global_grid%myYbegin - 1)
                END IF
             END IF
          END DO
       END DO
       DEALLOCATE(delta_lon)
       DEALLOCATE(delta_lat)
       DEALLOCATE(grid_lon)

       ! Global min/max reduce
       DO i=1, ntiles
          CALL mpp_min(ijmin(i,1))
          CALL mpp_max(ijmax(i,1))
          CALL mpp_min(ijmin(i,2))
          CALL mpp_max(ijmax(i,2))
       END DO

       IF ( ijmin(mytile,1) == huge(1) .OR. ijmax(mytile,1) == -huge(1) ) THEN
          ijmin(mytile,1) = 0
          ijmax(mytile,1) = 0
       END IF
       IF ( ijmin(mytile,2) == huge(1) .OR. ijmax(mytile,2) == -huge(1) ) THEN
          ijmin(mytile,2) = 0
          ijmax(mytile,2) = 0
       END IF

       istart = ijmin(mytile,1)
       jstart = ijmin(mytile,2)
       iend = ijmax(mytile,1)
       jend = ijmax(mytile,2)

       DEALLOCATE(ijmin)
       DEALLOCATE(ijmax)
    END IF

  END SUBROUTINE get_local_indexes
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="get_local_indexes2">
  !   <OVERVIEW>
  !     Find the indices of the nearest grid point of the a-grid to the
  !     specified (lon,lat) location on the local PE. if desired point not
  !     within domain of local PE, return (0,0) as the indices.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE get_local_indexes2 (lat, lon, iindex, jindex)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Given a specified location, find the nearest a-grid indices on
  !     the local PE.
  !   </DESCRIPTION>
  !   <IN NAME="lat" TYPE="REAL">
  !     The requested latitude.  This value must be in the range [-90,90]
  !   </IN>
  !   <IN NAME="lon" TYPE="REAL">
  !     The requested longitude.
  !     Possible ranges are either [-180,180] or [0,360].
  !   </IN>
  !   <OUT NAME="iindex" TYPE="INTEGER">
  !     The local index on the local PE in the 'i' direction.
  !   </OUT>
  !   <OUT NAME="jindex" TYPE="INTEGER">
  !     The local index on the local PE in the 'j' direction.
  !   </OUT>
  SUBROUTINE get_local_indexes2(lat, lon, iindex, jindex)
    REAL, INTENT(in) :: lat, lon !< lat/lon location
    INTEGER, INTENT(out) :: iindex, jindex !< i/j indexes

    INTEGER  :: indexes(2)

    IF ( .NOT. diag_grid_initialized )&
         & CALL error_mesg('diag_grid_mod::get_local_indexes2',&
         &'Module not initialized, first initialze module with a call &
         &to diag_grid_init', fatal)

    indexes = 0

    IF ( mod(diag_global_grid%tile_number,3) == 0 ) THEN
       IF ( lat > 30.0 .AND. diag_global_grid%tile_number == 3 ) THEN
          indexes(:) = find_pole_index_agrid(lat,lon)
       ELSE IF ( lat < -30.0 .AND. diag_global_grid%tile_number == 6 ) THEN
          indexes(:) = find_pole_index_agrid(lat,lon)
       ENDIF
    ELSE
       indexes(:) = find_equator_index_agrid(lat,lon)
    END IF

    iindex = indexes(1)
    jindex = indexes(2)
    IF (iindex ==  diag_global_grid%adimI -1 .OR.&
        jindex ==  diag_global_grid%adimJ -1 ) THEN
      iindex = 0
      jindex = 0
    ENDIF

  END SUBROUTINE get_local_indexes2
  ! </SUBROUTINE>

  ! <PRIVATE>
  ! <FUNCTION NAME="rad2deg">
  !   <OVERVIEW>
  !     Convert and angle in radian to degrees.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     PURE ELEMENTAL REAL FUNCTION rad2deg(angle)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Given a scalar, or an array of angles in radians this
  !     function will return a scalar or array (of the same
  !     dimension) of angles in degrees.
  !   </DESCRIPTION>
  !   <IN NAME="angle" TYPE="REAL">
  !     Scalar or array of angles in radians.
  !   </IN>
  !   <OUT NAME="rad2deg" TYPE="REAL">
  !     Scalar or array (depending on the size of angle) of angles in
  !     degrees.
  !   </OUT>
  PURE ELEMENTAL REAL FUNCTION rad2deg(angle)
    REAL, INTENT(in) :: angle

    rad2deg = rad_to_deg * angle
  END FUNCTION rad2deg
  ! </FUNCTION>
  ! </PRIVATE>

  ! <PRIVATE>
  ! <FUNCTION NAME="deg2rad">
  !   <OVERVIEW>
  !     Convert an angle in degrees to radians.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     PURE ELEMENTAL REAL FUNCTION deg2rad(angle)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Given a scalar, or an array of angles in degrees this
  !     function will return a scalar or array (of the same
  !     dimension) of angles in radians.
  !   </DESCRIPTION>
  !   <IN NAME="angle" TYPE="REAL">
  !     Scalar or array of angles in degrees.
  !   </IN>
  PURE ELEMENTAL REAL FUNCTION deg2rad(angle)
    REAL, INTENT(in) :: angle

    deg2rad = deg_to_rad * angle
  END FUNCTION deg2rad
  ! </FUNCTION>
  ! </PRIVATE>

  ! <PRIVATE>
  ! <FUNCTION NAME="find_pole_index_agrid">
  !   <OVERVIEW>
  !     Return the closest index (i,j) to the given (lat,lon) point.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     PURE FUNCTION find_pole_index_agrid(lat, lon)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     This function searches a pole a-grid tile looking for the grid point
  !     closest to the give (lat, lon) location, and returns the i
  !     and j indexes of the point.
  !   </DESCRIPTION>
  !   <IN NAME="lat" TYPE="REAL">
  !     Latitude location
  !   </IN>
  !   <IN NAME="lon" TYPE="REAL">
  !     Longitude location
  !   </IN>
  !   <OUT NAME="find_pole_index" TYPE="INTEGER, DIMENSION(2)">
  !     The (i, j) location of the closest grid to the given (lat,
  !     lon) location.
  !   </OUT>
  PURE FUNCTION find_pole_index_agrid(lat, lon)
    INTEGER, DIMENSION(2) :: find_pole_index_agrid
    REAL, INTENT(in) :: lat, lon

    INTEGER :: indxi, indxj !< Indexes to be returned.
    INTEGER :: dimi, dimj !< Size of the grid dimensions
    INTEGER :: i,j !< Count indexes
    INTEGER :: nearestcorner !< index of the nearest corner.
    INTEGER, DIMENSION(4,2) :: ijarray !< indexes of the cornerPts and pntDistances arrays
    REAL :: llat, llon
    REAL :: maxctrdist !< maximum distance to center of grid
    REAL, DIMENSION(4) :: pntdistances !< distance from origPt to corner
    TYPE(point) :: origpt !< Original point
    REAL, DIMENSION(4,2) :: cornerpts !< Corner points using (lat,lon)
    TYPE(point), DIMENSION(9) :: points !< xyz of 8 nearest neighbors
    REAL, DIMENSION(9) :: distsqrd !< distance between origPt and points(:)

    ! Set the inital fail values for indxI and indxJ
    indxi = 0
    indxj = 0

    dimi = diag_global_grid%adimI
    dimj = diag_global_grid%adimJ

    ! Since the poles have an non-unique longitude value, make a small correction if looking for one of the poles.
    IF ( lat == 90.0 ) THEN
       llat = lat - .1
    ELSE IF ( lat == -90.0 ) THEN
       llat = lat + .1
    ELSE
       llat = lat
    END IF
    llon = lon

    origpt = latlon2xyz(llat,llon)

    iloop: DO i=0, dimi-2
       jloop: DO j = 0, dimj-2
          cornerpts = reshape( (/ diag_global_grid%aglo_lat(i,  j),  diag_global_grid%aglo_lon(i,  j),&
               &                  diag_global_grid%aglo_lat(i+1,j+1),diag_global_grid%aglo_lon(i+1,j+1),&
               &                  diag_global_grid%aglo_lat(i+1,j),  diag_global_grid%aglo_lon(i+1,j),&
               &                  diag_global_grid%aglo_lat(i,  j+1),diag_global_grid%aglo_lon(i,  j+1) /),&
               &               (/ 4, 2 /), order=(/2,1/) )
          ! Find the maximum half distance of the corner points
          maxctrdist = max(gcirdistance(cornerpts(1,1),cornerpts(1,2), cornerpts(2,1),cornerpts(2,2)),&
               &           gcirdistance(cornerpts(3,1),cornerpts(3,2), cornerpts(4,1),cornerpts(4,2)))

          ! Find the distance of the four corner points to the point of interest.
          pntdistances = gcirdistance(cornerpts(:,1),cornerpts(:,2), llat,llon)

          IF ( (minval(pntdistances) <= maxctrdist) .AND. (i*j.NE.0) ) THEN
             ! Set up the i,j index array
             ijarray = reshape( (/ i, j, i+1, j+1, i+1, j, i, j+1 /), (/ 4, 2 /), order=(/2,1/) )

             ! the nearest point index
             nearestcorner = minloc(pntdistances,1)

             indxi = ijarray(nearestcorner,1)
             indxj = ijarray(nearestcorner,2)

             EXIT iloop
          END IF
       END DO jloop
    END DO iloop


    ! Make sure we have indexes in the correct range
    valid: IF (  (indxi <= 0 .OR. dimi-1 <= indxi) .OR. &
         &       (indxj <= 0 .OR. dimj-1 <= indxj) ) THEN
       indxi = 0
       indxj = 0
    ELSE ! indxI and indxJ are valid.
       ! Since we are looking for the closest grid point to the
       ! (lat,lon) point, we need to check the surrounding
       ! points.  The indexes for the variable points are as follows
       !
       ! 1---4---7
       ! |   |   |
       ! 2---5---8
       ! |   |   |
       ! 3---6---9

       ! Set the 'default' values for points(:) x,y,z to some large
       ! value.
       DO i=1, 9
          points(i)%x = 1.0e20
          points(i)%y = 1.0e20
          points(i)%z = 1.0e20
       END DO

       ! All the points around the i,j indexes
       points(1) = latlon2xyz(diag_global_grid%aglo_lat(indxi-1,indxj+1),&
            &                 diag_global_grid%aglo_lon(indxi-1,indxj+1))
       points(2) = latlon2xyz(diag_global_grid%aglo_lat(indxi-1,indxj),&
            &                 diag_global_grid%aglo_lon(indxi-1,indxj))
       points(3) = latlon2xyz(diag_global_grid%aglo_lat(indxi-1,indxj-1),&
            &                 diag_global_grid%aglo_lon(indxi-1,indxj-1))
       points(4) = latlon2xyz(diag_global_grid%aglo_lat(indxi,  indxj+1),&
            &                 diag_global_grid%aglo_lon(indxi,  indxj+1))
       points(5) = latlon2xyz(diag_global_grid%aglo_lat(indxi,  indxj),&
            &                 diag_global_grid%aglo_lon(indxi,  indxj))
       points(6) = latlon2xyz(diag_global_grid%aglo_lat(indxi,  indxj-1),&
            &                 diag_global_grid%aglo_lon(indxi,  indxj-1))
       points(7) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj+1),&
            &                 diag_global_grid%aglo_lon(indxi+1,indxj+1))
       points(8) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj),&
            &                 diag_global_grid%aglo_lon(indxi+1,indxj))
       points(9) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj-1),&
            &                 diag_global_grid%aglo_lon(indxi+1,indxj-1))


       ! Calculate the distance squared between the points(:) and the origPt
       distsqrd = distancesqrd(origpt, points)

       SELECT CASE (minloc(distsqrd,1))
       CASE ( 1 )
          indxi = indxi-1
          indxj = indxj+1
       CASE ( 2 )
          indxi = indxi-1
          indxj = indxj
       CASE ( 3 )
          indxi = indxi-1
          indxj = indxj-1
       CASE ( 4 )
          indxi = indxi
          indxj = indxj+1
       CASE ( 5 )
          indxi = indxi
          indxj = indxj
       CASE ( 6 )
          indxi = indxi
          indxj = indxj-1
       CASE ( 7 )
          indxi = indxi+1
          indxj = indxj+1
       CASE ( 8 )
          indxi = indxi+1
          indxj = indxj
       CASE ( 9 )
          indxi = indxi+1
          indxj = indxj-1
       CASE DEFAULT
          indxi = 0
          indxj = 0
       END SELECT
    END IF valid

    ! Set the return value for the funtion
    find_pole_index_agrid = (/indxi, indxj/)
  END FUNCTION find_pole_index_agrid
  ! </FUNCTION>
  ! </PRIVATE>

  ! <PRIVATE>
  ! <FUNCTION NAME="find_equator_index_agrid">
  !   <OVERVIEW>
  !     Return the closest index (i,j) to the given (lat,lon) point.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     PURE FUNCTION find_equator_index_agrid(lat, lon)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     This function searches a equator grid tile looking for the grid point
  !     closest to the give (lat, lon) location, and returns the i
  !     and j indexes of the point.
  !   </DESCRIPTION>
  !   <IN NAME="lat" TYPE="REAL">
  !     Latitude location
  !   </IN>
  !   <IN NAME="lon" TYPE="REAL">
  !     Longitude location
  !   </IN>
  !   <OUT NAME="find_equator_index" TYPE="INTEGER, DIMENSION(2)">
  !     The (i, j) location of the closest grid to the given (lat,
  !     lon) location.
  !   </OUT>
  PURE FUNCTION find_equator_index_agrid(lat, lon)
    INTEGER, DIMENSION(2) :: find_equator_index_agrid
    REAL, INTENT(in) :: lat, lon

    INTEGER :: indxi, indxj !< Indexes to be returned.
    INTEGER :: indxi_tmp !< Hold the indxI value if on tile 3 or 4
    INTEGER :: dimi, dimj !< Size of the grid dimensions
    INTEGER :: i,j !< Count indexes
    INTEGER :: jstart, jend, nextj !< j counting variables
    TYPE(point) :: origpt !< Original point
    TYPE(point), DIMENSION(4) :: points !< xyz of 8 nearest neighbors
    REAL, DIMENSION(4) :: distsqrd !< distance between origPt and points(:)

    ! Set the inital fail values for indxI and indxJ
    indxi = 0
    indxj = 0

    dimi = diag_global_grid%adimI
    dimj = diag_global_grid%adimJ

    ! check to see if the 'fix' for the latitude index is needed
    IF ( diag_global_grid%aglo_lat(1,1) > &
         &diag_global_grid%aglo_lat(1,2) ) THEN
       ! reverse the j search
       jstart = dimj-1
       jend = 1
       nextj = -1
    ELSE
       jstart = 0
       jend = dimj-2
       nextj = 1
    END IF

    ! find the I index
    iloop: DO i=0, dimi-2
       IF (   diag_global_grid%aglo_lon(i,0) >&
            & diag_global_grid%aglo_lon(i+1,0) ) THEN
          ! We are at the 0 longitudal line
          IF (   (diag_global_grid%aglo_lon(i,0) <= lon .AND. lon <= 360.) .OR.&
               & (0. <= lon .AND. lon < diag_global_grid%aglo_lon(i+1, 0)) ) THEN
             indxi = i
             EXIT iloop
          END IF
       ELSEIF ( diag_global_grid%aglo_lon(i,0) <= lon .AND.&
            &   lon <= diag_global_grid%aglo_lon(i+1,0) ) THEN
          indxi = i
          EXIT iloop
       END IF
    END DO iloop

    ! Find the J index
    IF ( indxi > 0 ) THEN
       jloop: DO j=jstart, jend, nextj
          IF (   diag_global_grid%aglo_lat(indxi,j) <= lat .AND.&
               & lat <= diag_global_grid%aglo_lat(indxi,j+nextj) ) THEN
             indxj = j
             EXIT jloop
          END IF
       END DO jloop
    END IF

    ! Make sure we have indexes in the correct range
    valid: IF ( (indxi <= 0 .OR. dimi-1 < indxi) .OR. &
         &      (indxj <= 0 .OR. dimj-1 < indxj) ) THEN
       indxi = 0
       indxj = 0
    ELSE ! indxI and indxJ are valid.
       ! Since we are looking for the closest grid point to the
       ! (lat,lon) point, we need to check the surrounding
       ! points.  The indexes for the variable points are as follows
       !
       ! 1---3
       ! |   |
       ! 2---4

       ! The original point
       origpt = latlon2xyz(lat,lon)

       ! Set the 'default' values for points(:) x,y,z to some large
       ! value.
       DO i=1, 4
          points(i)%x = 1.0e20
          points(i)%y = 1.0e20
          points(i)%z = 1.0e20
       END DO

       ! The original point
       origpt = latlon2xyz(lat,lon)

       points(1) = latlon2xyz(diag_global_grid%aglo_lat(indxi,indxj),&
            &                 diag_global_grid%aglo_lon(indxi,indxj))
       points(2) = latlon2xyz(diag_global_grid%aglo_lat(indxi,indxj+nextj),&
            &                 diag_global_grid%aglo_lon(indxi,indxj+nextj))
       points(3) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj+nextj),&
            &                 diag_global_grid%aglo_lon(indxi+1,indxj+nextj))
       points(4) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj),&
            &                 diag_global_grid%aglo_lon(indxi+1,indxj))

       ! Find the distance between the original point and the four
       ! grid points
       distsqrd = distancesqrd(origpt, points)

       SELECT CASE (minloc(distsqrd,1))
       CASE ( 1 )
          indxi = indxi;
          indxj = indxj;
       CASE ( 2 )
          indxi = indxi;
          indxj = indxj+nextj;
       CASE ( 3 )
          indxi = indxi+1;
          indxj = indxj+nextj;
       CASE ( 4 )
          indxi = indxi+1;
          indxj = indxj;
       CASE DEFAULT
          indxi = 0;
          indxj = 0;
       END SELECT

       ! If we are on tile 3 or 4, then the indxI and indxJ are
       ! reversed due to the transposed grids.
       IF (   diag_global_grid%tile_number == 4 .OR.&
            & diag_global_grid%tile_number == 5 ) THEN
          indxi_tmp = indxi
          indxi = indxj
          indxj = indxi_tmp
       END IF
    END IF valid

    ! Set the return value for the function
    find_equator_index_agrid = (/indxi, indxj/)
  END FUNCTION find_equator_index_agrid
  ! </FUNCTION>
  ! </PRIVATE>

  ! <PRIVATE>
  ! <FUNCTION NAME="latlon2xyz">
  !   <OVERVIEW>
  !     Return the (x,y,z) position of a given (lat,lon) point.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     PURE ELEMENTAL TYPE(point) FUNCTION latlon2xyz(lat, lon)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Given a specific (lat, lon) point on the Earth, return the
  !     corresponding (x,y,z) location.  The return of latlon2xyz
  !     will be either a scalar or an array of the same size as lat
  !     and lon.
  !   </DESCRIPTION>
  !   <IN NAME="lat" TYPE="REAL">
  !     The latitude of the (x,y,z) location to find.  <TT>lat</TT>
  !     can be either a scalar or array.  <TT>lat</TT> must be of the
  !     same rank / size as <TT>lon</TT>.  This function assumes
  !     <TT>lat</TT> is in the range [-90,90].
  !   </IN>
  !   <IN NAME="lon" TYPE="REAL">
  !     The longitude of the (x,y,z) location to find.  <TT>lon</TT>
  !     can be either a scalar or array.  <TT>lon</TT> must be of the
  !     same rank / size as <TT>lat</TT>.  This function assumes
  !     <TT>lon</TT> is in the range [0,360].
  !   </IN>
  PURE ELEMENTAL TYPE(point) function latlon2xyz(lat, lon)
    REAL, INTENT(in) :: lat, lon

    ! lat/lon angles in radians
    REAL :: theta, phi

    ! Convert the lat lon values to radians The lat values passed in
    ! are in the range [-90,90], but we need to have a radian range
    ! [0,pi], where 0 is at the north pole.  This is the reason for
    ! the subtraction from 90
    theta = deg2rad(90.-lat)
    phi = deg2rad(lon)

    ! Calculate the x,y,z point
    latlon2xyz%x = radius * sin(theta) * cos(phi)
    latlon2xyz%y = radius * sin(theta) * sin(phi)
    latlon2xyz%z = radius * cos(theta)
  END FUNCTION latlon2xyz
  ! </FUNCTION>
  ! </PRIVATE>

  ! <PRIVATE>
  ! <FUNCTION NAME="distanceSqrd">
  !   <OVERVIEW>
  !     Find the distance between two points in the Cartesian
  !     coordinate space.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     PURE ELEMENTAL REAL FUNCTION distanceSqrd(pt1, pt2)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     <TT>distanceSqrd</TT> will find the distance squared between
  !     two points in the xyz coordinate space.  <TT>pt1</TT> and <TT>
  !     pt2</TT> can either be both scalars, both arrays of the same
  !     size, or one a scalar and one an array.  The return value
  !     will be a scalar or array of the same size as the input array.
  !   </DESCRIPTION>
  !   <IN NAME="pt1" TYPE="TYPE(POINT)" />
  !   <IN NAME="pt2" TYPE="TYPE(POINT)" />
  PURE ELEMENTAL REAL FUNCTION distancesqrd(pt1, pt2)
    TYPE(point), INTENT(in) :: pt1, pt2

    distancesqrd = (pt1%x-pt2%x)**2 +&
         &         (pt1%y-pt2%y)**2 +&
         &         (pt1%z-pt2%z)**2
  END FUNCTION distancesqrd
  ! </FUNCTION>
  ! </PRIVATE>

  ! <FUNCTION NAME="gCirDistance">
  !   <OVERVIEW>
  !     Find the distance, along the geodesic, between two points.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     PURE ELEMENTAL REAL FUNCTION gCirDistance(lat1, lon1, lat2, lon2)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     <TT>aCirDistance</TT> will find the distance, along the geodesic, between two points defined by the (lat,lon) position of
  !     each point.
  !   </DESCRIPTION>
  !   <IN NAME="lat1" TYPE="REAL">Latitude of the first point</IN>
  !   <IN NAME="lon1" TYPE="REAL">Longitude of the first point</IN>
  !   <IN NAME="lat2" TYPE="REAL">Latitude of the second point</IN>
  !   <IN NAME="lon2" TYPE="REAL">Longitude of the second point</IN>
  PURE ELEMENTAL REAL FUNCTION gcirdistance(lat1, lon1, lat2, lon2)
    REAL, INTENT(in) :: lat1, lat2, lon1, lon2

    REAL :: theta1, theta2
    REAL :: deltalambda !< Difference in longitude angles, in radians.
    REAL :: deltatheta !< Difference in latitude angels, in radians.

    theta1 = deg2rad(lat1)
    theta2 = deg2rad(lat2)
    deltalambda = deg2rad(lon2-lon1)
    deltatheta = deg2rad(lat2-lat1)

    gcirdistance = radius * 2. * asin(sqrt((sin(deltatheta/2.))**2 + cos(theta1)*cos(theta2)*(sin(deltalambda/2.))**2))
  END FUNCTION gcirdistance
  ! </FUNCTION>

  !Find the i,j indices and distance of the a-grid point nearest to
  !the inputted lat,lon point.
  SUBROUTINE find_nearest_agrid_index(lat, &
                                      lon, &
                                      minI, &
                                      minJ, &
                                      minimum_distance)

    !Inputs/outputs
    REAL,INTENT(IN) :: lat
    REAL,INTENT(IN) :: lon
    INTEGER,INTENT(OUT) :: minI
    INTEGER,INTENT(OUT) :: minJ
    REAL,INTENT(OUT) :: minimum_distance

    !Local variables
    REAL :: llat
    REAL :: llon
    INTEGER :: j
    INTEGER :: i
    REAL :: dist

    !Since the poles have an non-unique longitude value, make a small
    !correction if looking for one of the poles.
    IF (lat .EQ. 90.0) THEN
       llat = lat - .1
    ELSEIF (lat .EQ. -90.0) THEN
       llat = lat + .1
    ELSE
       llat = lat
    END IF
    llon = lon

    !Loop through non-halo points.  Calculate the distance
    !between each a-grid point and the point that we
    !are seeking.  Store the minimum distance and its
    !corresponding i,j indices.
    mini = 0
    minj = 0
    minimum_distance = 2.0*radius*3.141592653
    DO j = 1,diag_global_grid%adimJ-2
        DO i = 1,diag_global_grid%adimI-2
            dist = gcirdistance(llat, &
                                llon, &
                                diag_global_grid%aglo_lat(i,j), &
                                diag_global_grid%aglo_lon(i,j))
            IF (dist .LT. minimum_distance) THEN

                !These number shouldn't be hardcoded, but they have to
                !match the ones in diag_grid_init.
                if (diag_global_grid%tile_number .eq. 4 .or. &
                        diag_global_grid%tile_number .eq. 5) then

                    !Because of transpose in diag_grid_init.
                    mini = j
                    minj = i

                else
                    mini = i
                    minj = j
                endif
                minimum_distance = dist
            ENDIF
        ENDDO
    ENDDO

    !Check that valid i,j indices have been found.
    IF (mini .EQ. 0 .OR. minj .EQ. 0) THEN
        call error_mesg("find_nearest_agrid_index", &
                        "A minimum distance was not found.", &
                        fatal)
    ENDIF

  END SUBROUTINE find_nearest_agrid_index

END MODULE diag_grid_mod