MAPL_LatLonGridFactory.F90 Source File


This file depends on

sourcefile~~mapl_latlongridfactory.f90~~EfferentGraph sourcefile~mapl_latlongridfactory.f90 MAPL_LatLonGridFactory.F90 sourcefile~base_base.f90 Base_Base.F90 sourcefile~mapl_latlongridfactory.f90->sourcefile~base_base.f90 sourcefile~constants.f90 Constants.F90 sourcefile~mapl_latlongridfactory.f90->sourcefile~constants.f90 sourcefile~mapl_abstractgridfactory.f90 MAPL_AbstractGridFactory.F90 sourcefile~mapl_latlongridfactory.f90->sourcefile~mapl_abstractgridfactory.f90 sourcefile~mapl_comms.f90 MAPL_Comms.F90 sourcefile~mapl_latlongridfactory.f90->sourcefile~mapl_comms.f90 sourcefile~mapl_config.f90 MAPL_Config.F90 sourcefile~mapl_latlongridfactory.f90->sourcefile~mapl_config.f90 sourcefile~mapl_errorhandling.f90 MAPL_ErrorHandling.F90 sourcefile~mapl_latlongridfactory.f90->sourcefile~mapl_errorhandling.f90 sourcefile~mapl_keywordenforcer.f90 MAPL_KeywordEnforcer.F90 sourcefile~mapl_latlongridfactory.f90->sourcefile~mapl_keywordenforcer.f90 sourcefile~mapl_minmax.f90 MAPL_MinMax.F90 sourcefile~mapl_latlongridfactory.f90->sourcefile~mapl_minmax.f90 sourcefile~pfio.f90 pFIO.F90 sourcefile~mapl_latlongridfactory.f90->sourcefile~pfio.f90

Files dependent on this one

sourcefile~~mapl_latlongridfactory.f90~~AfferentGraph sourcefile~mapl_latlongridfactory.f90 MAPL_LatLonGridFactory.F90 sourcefile~base.f90 Base.F90 sourcefile~base.f90->sourcefile~mapl_latlongridfactory.f90 sourcefile~mapl_cfio.f90 MAPL_CFIO.F90 sourcefile~mapl_cfio.f90->sourcefile~mapl_latlongridfactory.f90 sourcefile~mapl_gridmanager.f90 MAPL_GridManager.F90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_latlongridfactory.f90 sourcefile~regrid_util.f90 Regrid_Util.F90 sourcefile~regrid_util.f90->sourcefile~mapl_latlongridfactory.f90 sourcefile~test_latlon_corners.pf Test_LatLon_Corners.pf sourcefile~test_latlon_corners.pf->sourcefile~mapl_latlongridfactory.f90 sourcefile~test_latlon_gridfactory.pf Test_LatLon_GridFactory.pf sourcefile~test_latlon_gridfactory.pf->sourcefile~mapl_latlongridfactory.f90 sourcefile~test_sphericaltocartesian.pf Test_SphericalToCartesian.pf sourcefile~test_sphericaltocartesian.pf->sourcefile~mapl_latlongridfactory.f90

Source Code

#include "MAPL_ErrLog.h"

! overload set interfaces in legacy
! Document PE, PC, DC, DE, GC

! This module generates ESMF_Grids corresponding to _regular_ lat-lon coordinate grids.
! I.e., spacing between lats (lons) is constant.

module MAPL_LatLonGridFactoryMod
   use MAPL_AbstractGridFactoryMod
   use MAPL_MinMaxMod
   use MAPL_KeywordEnforcerMod
   use mapl_ErrorHandlingMod
   use MAPL_Constants
   use ESMF
   use pFIO
   use MAPL_CommsMod
   use, intrinsic :: iso_fortran_env, only: REAL32
   use, intrinsic :: iso_fortran_env, only: REAL64
   implicit none
   private

   public :: LatLonGridFactory

   integer, parameter :: NUM_DIM = 2

   type, extends(AbstractGridFactory) :: LatLonGridFactory
      private
      logical :: is_regular = .false.
      character(len=:), allocatable :: grid_name
      ! Grid dimensions
      integer :: im_world = MAPL_UNDEFINED_INTEGER
      integer :: jm_world = MAPL_UNDEFINED_INTEGER
      integer :: lm = MAPL_UNDEFINED_INTEGER
      real(kind=REAL64), allocatable :: lon_centers(:)
      real(kind=REAL64), allocatable :: lat_centers(:)
      real(kind=REAL64), allocatable :: lon_centers_degrees(:)
      real(kind=REAL64), allocatable :: lat_centers_degrees(:)
      real(kind=REAL64), allocatable :: lon_corners(:)
      real(kind=REAL64), allocatable :: lat_corners(:)
      logical :: force_decomposition = .false.

      ! Domain decomposition:
      integer :: nx = MAPL_UNDEFINED_INTEGER
      integer :: ny = MAPL_UNDEFINED_INTEGER
      integer, allocatable :: ims(:)
      integer, allocatable :: jms(:)
      ! Grid conventions:
      character(len=:), allocatable :: pole
      character(len=:), allocatable :: dateline
      ! Regional vs global:
      type (RealMinMax) :: lon_range = RealMinMax(MAPL_UNDEFINED_REAL,MAPL_UNDEFINED_REAL)
      type (RealMinMax) :: lat_range = RealMinMax(MAPL_UNDEFINED_REAL,MAPL_UNDEFINED_REAL)
      ! Used for halo
      type (ESMF_DELayout) :: layout
      integer :: px, py
      logical :: is_halo_initialized = .false.
      logical :: periodic = .true.
   contains
      procedure :: make_new_grid
      procedure :: create_basic_grid
      procedure :: get_longitudes_degrees
      procedure :: get_longitudes
      procedure :: get_latitudes_degrees
      procedure :: get_latitudes
      procedure :: compute_lon_centers
      procedure :: compute_lat_centers
      procedure :: get_lon_corners
      procedure :: get_lat_corners
      procedure :: compute_lon_corners
      procedure :: compute_lat_corners
      procedure :: add_horz_coordinates
      procedure :: init_halo
      procedure :: halo


      procedure :: initialize_from_file_metadata
      procedure :: initialize_from_config_with_prefix
      procedure :: initialize_from_esmf_distGrid

      procedure :: equals

      procedure :: check_and_fill_consistency
      procedure :: generate_grid_name
      procedure :: to_string

      procedure :: append_metadata
      procedure :: get_grid_vars
      procedure :: get_file_format_vars
      procedure :: append_variable_metadata
      procedure :: check_decomposition
      procedure :: generate_newnxy
      procedure :: generate_file_bounds
      procedure :: generate_file_corner_bounds
      procedure :: generate_file_reference2D
      procedure :: generate_file_reference3D
      procedure :: decomps_are_equal
      procedure :: physical_params_are_equal
   end type LatLonGridFactory

   character(len=*), parameter :: MOD_NAME = 'MAPL_LatLonGridFactory::'

   interface LatLonGridFactory
      module procedure LatLonGridFactory_from_parameters
   end interface LatLonGridFactory

   interface set_with_default
      module procedure set_with_default_integer
      module procedure set_with_default_real
      module procedure set_with_default_character
      module procedure set_with_default_range
      module procedure set_with_default_logical
   end interface set_with_default


contains


   function LatLonGridFactory_from_parameters(unusable, grid_name, &
        & im_world, jm_world, lm, nx, ny, ims, jms, &
        & pole, dateline, lon_range, lat_range, force_decomposition, rc) result(factory)
      type (LatLonGridFactory) :: factory
      class (KeywordEnforcer), optional, intent(in) :: unusable
      character(len=*), optional, intent(in) :: grid_name

      ! grid details:
      integer, optional, intent(in) :: im_world
      integer, optional, intent(in) :: jm_world
      integer, optional, intent(in) :: lm
      character(len=2), optional, intent(in) :: pole
      character(len=2), optional, intent(in) :: dateline
      type (RealMinMax), optional, intent(in) :: lon_range
      type (RealMinMax), optional, intent(in) :: lat_range

      ! decomposition:
      integer, optional, intent(in) :: nx
      integer, optional, intent(in) :: ny
      integer, optional, intent(in) :: ims(:)
      integer, optional, intent(in) :: jms(:)
      logical, optional, intent(in) :: force_decomposition

      integer, optional, intent(out) :: rc

      integer :: status

      _UNUSED_DUMMY(unusable)

      factory%is_regular = .true.
      call set_with_default(factory%grid_name, grid_name, MAPL_GRID_NAME_DEFAULT)

      call set_with_default(factory%nx, nx, MAPL_UNDEFINED_INTEGER)
      call set_with_default(factory%ny, ny, MAPL_UNDEFINED_INTEGER)

      call set_with_default(factory%im_world, im_world, MAPL_UNDEFINED_INTEGER)
      call set_with_default(factory%jm_world, jm_world, MAPL_UNDEFINED_INTEGER)
      call set_with_default(factory%lm, lm, MAPL_UNDEFINED_INTEGER)

      ! default is unallocated
      if (present(ims)) factory%ims = ims
      if (present(jms)) factory%jms = jms

      call set_with_default(factory%pole, pole, MAPL_UNDEFINED_CHAR)
      call set_with_default(factory%dateline, dateline, MAPL_UNDEFINED_CHAR)

      call set_with_default(factory%lon_range, lon_range, RealMinMax(MAPL_UNDEFINED_REAL,MAPL_UNDEFINED_REAL))
      call set_with_default(factory%lat_range, lat_range, RealMinMax(MAPL_UNDEFINED_REAL,MAPL_UNDEFINED_REAL))
      call set_with_default(factory%force_decomposition, force_decomposition, .false.)

      call factory%check_and_fill_consistency(rc=status)
      _VERIFY(status)

      ! Compute the centers and corners
      factory%lon_centers = factory%compute_lon_centers(factory%dateline, rc=status)
      _VERIFY(status)
      factory%lat_centers = factory%compute_lat_centers(factory%pole, rc=status)
      _VERIFY(status)
      factory%lon_centers_degrees = factory%compute_lon_centers(factory%dateline, &
            convert_to_radians = .false.,  rc=status)
      _VERIFY(status)
      factory%lat_centers_degrees = factory%compute_lat_centers(factory%pole, &
            convert_to_radians = .false.,  rc=status)
      _VERIFY(status)
      factory%lon_corners = factory%compute_lon_corners(factory%dateline, rc=status)
      _VERIFY(status)
      factory%lat_corners = factory%compute_lat_corners(factory%pole, rc=status)
      _VERIFY(status)

      _RETURN(_SUCCESS)

   end function LatLonGridFactory_from_parameters


   function make_new_grid(this, unusable, rc) result(grid)
      type (ESMF_Grid) :: grid
      class (LatLonGridFactory), intent(in) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      integer :: status

      _UNUSED_DUMMY(unusable)
      grid = this%create_basic_grid(rc=status)
      _VERIFY(status)

      call this%add_horz_coordinates(grid, rc=status)
      _VERIFY(status)

      _RETURN(_SUCCESS)

   end function make_new_grid



   function create_basic_grid(this, unusable, rc) result(grid)
      type (ESMF_Grid) :: grid
      class (LatLonGridFactory), intent(in) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      integer :: status

      _UNUSED_DUMMY(unusable)

      if (this%periodic) then
         grid = ESMF_GridCreate1PeriDim( &
              & name = this%grid_name, &
              & countsPerDEDim1=this%ims, &
              & countsPerDEDim2=this%jms, &
              & indexFlag=ESMF_INDEX_DELOCAL, &
              & gridEdgeLWidth=[0,0], &
              & gridEdgeUWidth=[0,1], &
              & coordDep1=[1,2], &
              & coordDep2=[1,2], &
              & coordSys=ESMF_COORDSYS_SPH_RAD, &
              & rc=status)
         _VERIFY(status)
      else
         grid = ESMF_GridCreateNoPeriDim( &
              & name = this%grid_name, &
              & countsPerDEDim1=this%ims, &
              & countsPerDEDim2=this%jms, &
              & indexFlag=ESMF_INDEX_DELOCAL, &
              & gridEdgeLWidth=[0,0], &
              & gridEdgeUWidth=[1,1], &
              & coordDep1=[1,2], &
              & coordDep2=[1,2], &
              & coordSys=ESMF_COORDSYS_SPH_RAD, &
              & rc=status)
         _VERIFY(status)
      end if

      ! Allocate coords at default stagger location
      call ESMF_GridAddCoord(grid, rc=status)
      _VERIFY(status)
      call ESMF_GridAddCoord(grid, staggerloc=ESMF_STAGGERLOC_CORNER, rc=status)
      _VERIFY(status)


      if (this%lm /= MAPL_UNDEFINED_INTEGER) then
         call ESMF_AttributeSet(grid, name='GRID_LM', value=this%lm, rc=status)
         _VERIFY(status)
      end if

      call ESMF_AttributeSet(grid, 'GridType', 'LatLon', rc=status)
      _VERIFY(status)
      if (.not.this%periodic) then
         call ESMF_AttributeSet(grid, 'Global', .false., rc=status)
         _VERIFY(status)
      end if

      _RETURN(_SUCCESS)
   end function create_basic_grid

   ! in radians
   function get_longitudes(this, unusable, rc) result(longitudes)
      use MAPL_BaseMod
      class (LatLonGridFactory), intent(in) :: this
      real(kind=REAL64), allocatable :: longitudes(:)
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      _UNUSED_DUMMY(unusable)

      longitudes = this%lon_centers
      _RETURN(_SUCCESS)
   end function get_longitudes

   function get_longitudes_degrees(this, unusable, rc) result(longitudes)
      use MAPL_BaseMod
      class (LatLonGridFactory), intent(in) :: this
      real(kind=REAL64), allocatable :: longitudes(:)
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      _UNUSED_DUMMY(unusable)

      longitudes = this%lon_centers_degrees
      _RETURN(_SUCCESS)
   end function get_longitudes_degrees

   ! in radians
   function get_latitudes(this, unusable, rc) result(latitudes)
      use MAPL_BaseMod
      class (LatLonGridFactory), intent(in) :: this
      real(kind=REAL64), allocatable :: latitudes(:)
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      _UNUSED_DUMMY(unusable)

      latitudes = this%lat_centers
      _RETURN(_SUCCESS)
   end function get_latitudes

   function get_latitudes_degrees(this, unusable, rc) result(latitudes)
      use MAPL_BaseMod
      class (LatLonGridFactory), intent(in) :: this
      real(kind=REAL64), allocatable :: latitudes(:)
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      _UNUSED_DUMMY(unusable)

      latitudes = this%lat_centers_degrees
      _RETURN(_SUCCESS)
   end function get_latitudes_degrees

   ! in radians
   function compute_lon_centers(this, dateline, unusable, convert_to_radians, rc) result(lon_centers)
      use MAPL_Constants, only:MAPL_DEGREES_TO_RADIANS_R8
      use MAPL_BaseMod
      real(kind=REAL64), allocatable :: lon_centers(:)
      class (LatLonGridFactory), intent(in) :: this
      character(2), intent(in) :: dateline
      class (KeywordEnforcer), optional, intent(in) :: unusable
      logical, optional, intent(in)  :: convert_to_radians
      integer, optional, intent(out) :: rc

      real(kind=REAL64) :: delta, min_coord, max_coord
      logical :: local_convert_to_radians
      logical :: regional
      integer :: status

      _UNUSED_DUMMY(unusable)
      if (present(convert_to_radians)) then
         local_convert_to_radians = convert_to_radians
      else
         local_convert_to_radians = .true.
      end if

      allocate(lon_centers(this%im_world))

      regional  = (dateline == 'XY')
      if (regional) then
         delta = (this%lon_range%max - this%lon_range%min) / this%im_world
         min_coord = this%lon_range%min + delta/2
         max_coord = this%lon_range%max - delta/2
      else
         delta = 360.d0 / this%im_world
         select case (dateline)
         case ('DC')
            min_coord = -180.d0
            max_coord = +180.d0 - delta
         case ('DE')
            min_coord = -180.d0 + delta/2
            max_coord = +180.d0 - delta/2
         case ('GC')
            min_coord = 0.d0
            max_coord = 360.d0 - delta
         case ('GE')
            min_coord = delta/2
            max_coord = 360.d0 - delta/2
         end select
      end if

      if (local_convert_to_radians) then
         lon_centers = MAPL_Range(min_coord, max_coord, this%im_world, &
              & conversion_factor=MAPL_DEGREES_TO_RADIANS_R8, rc=status)
         _VERIFY(status)
      else
         lon_centers = MAPL_Range(min_coord, max_coord, this%im_world, rc=status)
         _VERIFY(status)
      end if

      _RETURN(_SUCCESS)
   end function compute_lon_centers

   function compute_lon_corners(this, dateline, unusable, rc) result(lon_corners)
      use MAPL_Constants, only:MAPL_DEGREES_TO_RADIANS_R8
      use MAPL_BaseMod
      real(kind=REAL64), allocatable :: lon_corners(:)
      class (LatLonGridFactory), intent(in) :: this
      character(2), intent(in) :: dateline
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      real(kind=REAL64) :: delta, min_coord, max_coord
      logical :: regional
      integer :: status

      _UNUSED_DUMMY(unusable)

      allocate(lon_corners(this%im_world+1))

      regional  = (dateline == 'XY')
      if (regional) then
         delta = (this%lon_range%max - this%lon_range%min) / this%im_world
         min_coord = this%lon_range%min
         max_coord = this%lon_range%max
      else
         delta = 360.d0 / this%im_world
         select case (dateline)
         case ('DC')
            min_coord = -180.d0 - delta/2
            max_coord = +180.d0 - delta/2
         case ('DE')
            min_coord = -180.d0
            max_coord = +180.d0
         case ('GC')
            min_coord = 0.d0-delta/2
            max_coord = 360.d0-delta/2
         case ('GE')
            min_coord = 0.d0
            max_coord = 360.d0 - delta
         end select
      end if

      lon_corners = MAPL_Range(min_coord, max_coord, this%im_world+1, &
           & conversion_factor=MAPL_DEGREES_TO_RADIANS_R8, rc=status)
      _VERIFY(status)

      _RETURN(_SUCCESS)
   end function compute_lon_corners


   ! in radians
   function get_lon_corners(this, unusable, rc) result(lon_corners)
      use MAPL_BaseMod
      class (LatLonGridFactory), intent(in) :: this
      real(kind=REAL64), allocatable :: lon_corners(:)
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      _UNUSED_DUMMY(unusable)

      lon_corners = this%lon_corners
      _RETURN(_SUCCESS)

   end function get_lon_corners


   ! in radians
   function get_lat_corners(this, unusable, rc) result(lat_corners)
      use MAPL_BaseMod
      class (LatLonGridFactory), intent(in) :: this
      real(kind=REAL64), allocatable :: lat_corners(:)
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      _UNUSED_DUMMY(unusable)

      lat_corners = this%lat_corners
      _RETURN(_SUCCESS)

   end function get_lat_corners


   function compute_lat_centers(this, pole, unusable, convert_to_radians, rc) result(lat_centers)
      use MAPL_Constants, only: MAPL_DEGREES_TO_RADIANS_R8
      use MAPL_BaseMod
      real(kind=REAL64), allocatable :: lat_centers(:)
      class (LatLonGridFactory), intent(in) :: this
      character(2), intent(in) :: pole
      class (KeywordEnforcer), optional, intent(in) :: unusable
      logical, optional, intent(in)  :: convert_to_radians
      integer, optional, intent(out) :: rc

      real(kind=REAL64) :: delta, min_coord, max_coord
      logical :: regional
      logical :: local_convert_to_radians
      integer :: status

      _UNUSED_DUMMY(unusable)
      if (present(convert_to_radians)) then
         local_convert_to_radians = convert_to_radians
      else
         local_convert_to_radians = .true.
      end if

      allocate(lat_centers(this%jm_world))

      regional  = (pole == 'XY')
      if (regional) then
         delta = (this%lat_range%max - this%lat_range%min) / this%jm_world
         min_coord = this%lat_range%min + delta/2
         max_coord = this%lat_range%max - delta/2
      else ! global grid

         select case (pole)
         case ('PE')
            delta = 180.d0 / this%jm_world
            min_coord = -90.d0 + delta/2
            max_coord = +90.d0 - delta/2
         case ('PC')
            _ASSERT(this%jm_world > 1,'degenerate grid')
            min_coord = -90.d0
            max_coord = +90.d0
         end select
      end if

      if (local_convert_to_radians) then
         lat_centers = MAPL_Range(min_coord, max_coord, this%jm_world, &
              & conversion_factor=MAPL_DEGREES_TO_RADIANS_R8, rc=status)
      else
         lat_centers = MAPL_Range(min_coord, max_coord, this%jm_world, rc=status)
      end if

      _RETURN(_SUCCESS)

   end function compute_lat_centers

   function compute_lat_corners(this, pole, unusable, rc) result(lat_corners)
      use MAPL_Constants, only: MAPL_DEGREES_TO_RADIANS_R8
      use MAPL_BaseMod
      real(kind=REAL64), allocatable :: lat_corners(:)
      class (LatLonGridFactory), intent(in) :: this
      character(2), intent(in) :: pole
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      real(kind=REAL64) :: delta, min_coord, max_coord
      logical :: regional

      integer :: status

      _UNUSED_DUMMY(unusable)

      allocate(lat_corners(this%jm_world+1))

      regional  = (pole == 'XY')
      if (regional) then
         delta = (this%lat_range%max - this%lat_range%min) / this%jm_world
         min_coord = this%lat_range%min
         max_coord = this%lat_range%max
      else ! global grid

         select case (pole)
         case ('PE')
            delta = 180.d0 / this%jm_world
            min_coord = -90.d0
            max_coord = +90.d0
         case ('PC')
            _ASSERT(this%jm_world > 1, 'degenerate grid')
            delta = 180.d0 / (this%jm_world-1)
            min_coord = -90.d0-delta/2
            max_coord = +90.d0+delta/2
         end select
      end if

      lat_corners = MAPL_Range(min_coord, max_coord, this%jm_world+1, &
           & conversion_factor=MAPL_DEGREES_TO_RADIANS_R8, rc=status)
      if (pole == 'PC') then
         lat_corners(1)=-90.d0*MAPL_DEGREES_TO_RADIANS_R8
         lat_corners(this%jm_world+1)=90.d0*MAPL_DEGREES_TO_RADIANS_R8
      end if

      _RETURN(_SUCCESS)

   end function compute_lat_corners


   subroutine add_horz_coordinates(this, grid, unusable, rc)
      use MAPL_BaseMod, only: MAPL_grid_interior
      class (LatLonGridFactory), intent(in) :: this
      type (ESMF_Grid), intent(inout) :: grid
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      integer :: i_1, i_n, j_1, j_n ! regional array bounds
      integer :: ic_1,ic_n,jc_1,jc_n ! regional corner bounds
      real(kind=ESMF_KIND_R8), pointer :: centers(:,:)
      real(kind=ESMF_KIND_R8), pointer :: corners(:,:)
      integer :: status
      integer :: i, j, ij(4)

      _UNUSED_DUMMY(unusable)

      call MAPL_grid_interior(grid, i_1, i_n, j_1, j_n)
      ij(1)=i_1
      ij(2)=i_n
      ij(3)=j_1
      ij(4)=j_n
      if (.not. any(ij == -1)) then
         if (this%periodic) then
            ic_1=i_1
            ic_n=i_n
         else
            ic_1=i_1
            if (i_n == this%im_world) then
               ic_n=i_n+1
            else
               ic_n=i_n
            end if
         end if

         jc_1=j_1
         if (j_n == this%jm_world) then
            jc_n=j_n+1
         else
            jc_n=j_n
         end if

         ! First we handle longitudes:
         call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, &
              staggerloc=ESMF_STAGGERLOC_CENTER, &
              farrayPtr=centers, rc=status)
         _VERIFY(status)
         call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, &
              staggerloc=ESMF_STAGGERLOC_CORNER, &
              farrayPtr=corners, rc=status)
         _VERIFY(status)
         do j = 1, size(centers,2)
            centers(:,j) = this%lon_centers(i_1:i_n)
         end do
         do j = 1, size(corners,2)
            corners(:,j) = this%lon_corners(ic_1:ic_n)
         end do

         ! Now latitudes
         call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, &
              staggerloc=ESMF_STAGGERLOC_CENTER, &
              farrayPtr=centers, rc=status)
         _VERIFY(status)
         call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, &
              staggerloc=ESMF_STAGGERLOC_CORNER, &
              farrayPtr=corners, rc=status)
         _VERIFY(status)

         do i = 1, size(centers,1)
            centers(i,:) = this%lat_centers(j_1:j_n)
         end do
         do i = 1, size(corners,1)
            corners(i,:) = this%lat_corners(jc_1:jc_n)
         end do
      end if

      _RETURN(_SUCCESS)

   end subroutine add_horz_coordinates


   subroutine initialize_from_file_metadata(this, file_metadata, unusable, force_file_coordinates, rc)
      use MAPL_KeywordEnforcerMod
      use MAPL_BaseMod, only: MAPL_DecomposeDim

      class (LatLonGridFactory), intent(inout)  :: this
      type (FileMetadata), target, intent(in) :: file_metadata
      class (KeywordEnforcer), optional, intent(in) :: unusable
      logical, optional, intent(in) :: force_file_coordinates
      integer, optional, intent(out) :: rc

      integer :: status

      class (CoordinateVariable), pointer :: v
      class (*), pointer :: ptr(:)

      character(:), allocatable :: lon_name
      character(:), allocatable :: lat_name
      character(:), allocatable :: lev_name
      integer :: i
      logical :: hasLon, hasLat, hasLongitude, hasLatitude, hasLev,hasLevel,regLat,regLon
      real(kind=REAL64) :: del12,delij

      integer :: i_min, i_max
      real(kind=REAL64) :: d_lat, d_lat_temp, extrap_lat
      logical :: is_valid, use_file_coords, compute_lons, compute_lats

      _UNUSED_DUMMY(unusable)

      if (present(force_file_coordinates)) then
         use_file_coords = force_file_Coordinates
      else
         use_file_coords = .false.
      end if

      ! Cannot assume that lats and lons are evenly spaced
      this%is_regular = .false.

      associate (im => this%im_world, jm => this%jm_world, lm => this%lm)
         lon_name = 'lon'
         hasLon = file_metadata%has_dimension(lon_name)
         if (hasLon) then
            im = file_metadata%get_dimension(lon_name, rc=status)
            _VERIFY(status)
         else
            lon_name = 'longitude'
            hasLongitude = file_metadata%has_dimension(lon_name)
            if (hasLongitude) then
               im = file_metadata%get_dimension(lon_name, rc=status)
               _VERIFY(status)
            else
               _FAIL('no longitude coordinate')
            end if
         end if
         lat_name = 'lat'
         hasLat = file_metadata%has_dimension(lat_name)
         if (hasLat) then
            jm = file_metadata%get_dimension(lat_name, rc=status)
            _VERIFY(status)
         else
            lat_name = 'latitude'
            hasLatitude = file_metadata%has_dimension(lat_name)
            if (hasLatitude) then
               jm = file_metadata%get_dimension(lat_name, rc=status)
               _VERIFY(status)
            else
               _FAIL('no latitude coordinate')
            end if
         end if
         hasLev=.false.
         hasLevel=.false.
         lev_name = 'lev'
         hasLev = file_metadata%has_dimension(lev_name)
         if (hasLev) then
            lm = file_metadata%get_dimension(lev_name,rc=status)
            _VERIFY(status)
         else
            lev_name = 'levels'
            hasLevel = file_metadata%has_dimension(lev_name)
            if (hasLevel) then
               lm = file_metadata%get_dimension(lev_name,rc=status)
               _VERIFY(status)
            end if
         end if

        ! TODO: if 'lat' and 'lon' are not present then
        ! assume ... pole/dateline are ?

        ! TODO: check radians vs degrees.  Assume degrees for now.


        ! TODO: modify CoordinateVariable so that get_coordinate_data() is overloaded
        ! for different types (as subroutine) to avoid casting here.
        ! TODO: add get_coordinate_variable() interface to avoid the need to cast
        v => file_metadata%get_coordinate_variable(lon_name, rc=status)
        _VERIFY(status)
        ptr => v%get_coordinate_data()
        _ASSERT(associated(ptr),'coordinate data not allocated')
        select type (ptr)
        type is (real(kind=REAL64))
           this%lon_centers = ptr
        type is (real(kind=REAL32))
           this%lon_centers = ptr
        class default
           _FAIL('unsuppoted type of data; must be REAL32 or REAL64')
        end select

        if (any((this%lon_centers(2:im)-this%lon_centers(1:im-1))<0)) then
           where(this%lon_centers > 180) this%lon_centers=this%lon_centers-360
        end if


        v => file_metadata%get_coordinate_variable(lat_name, rc=status)
        _VERIFY(status)
        ptr => v%get_coordinate_data()
        _ASSERT(associated(ptr),'coordinate data not allocated')
        select type (ptr)
        type is (real(kind=REAL64))
           this%lat_centers = ptr
        type is (real(kind=REAL32))
           this%lat_centers = ptr
        class default
           _FAIL('unsupported type of data; must be REAL32 or REAL64')
        end select


        ! Check: is this a "mis-specified" pole-centered grid?
        if (size(this%lat_centers) >= 4) then
           ! Assume lbound=1 and ubound=size for now
           i_min = 1 !lbound(this%lat_centers)
           i_max = size(this%lat_centers) !ubound(this%lat_centers)
           d_lat = (this%lat_centers(i_max-1) - this%lat_centers(i_min+1))/&
                    (size(this%lat_centers)-3)
           is_valid = .True.
           ! Check: is this a regular grid (i.e. constant spacing away from the poles)?
           do i=(i_min+1),(i_max-2)
              d_lat_temp = this%lat_centers(i+1) - this%lat_centers(i)
              is_valid = (is_valid.and.(abs((d_lat_temp/d_lat)-1.0) < 1.0e-5))
              if (.not. is_valid) then
                 exit
              end if
           end do
           if (is_valid) then
              ! Should the southernmost point actually be at the pole?
              extrap_lat = this%lat_centers(i_min+1) - d_lat
              if (extrap_lat <= ((d_lat/20.0)-90.0)) then
                 this%lat_centers(i_min) = -90.0
              end if
              ! Should the northernmost point actually be at the pole?
              extrap_lat = this%lat_centers(i_max-1) + d_lat
              if (extrap_lat >= (90.0-(d_lat/20.0))) then
                 this%lat_centers(i_max) =  90.0
              end if
           end if
        end if


         ! Corners are the midpoints of centers (and extrapolated at the
         ! poles for lats.)
         allocate(this%lon_corners(im+1), this%lat_corners(jm+1))

         this%lon_corners(1) = (this%lon_centers(im) + this%lon_centers(1))/2 - 180
         this%lon_corners(2:im) = (this%lon_centers(1:im-1) + this%lon_centers(2:im))/2
         this%lon_corners(im+1) = (this%lon_centers(im) + this%lon_centers(1))/2 + 180

         ! This section about pole/dateline is probably not needed in file data case.
         if (abs(this%lon_centers(1) + 180) < 1000*epsilon(1.0)) then
            this%dateline = 'DC'
         else if (abs(this%lon_centers(1)) < 1000*epsilon(1.0)) then
            this%dateline = 'GC'
         else if (abs(this%lon_corners(1) + 180) < 1000*epsilon(1.0)) then
            this%dateline = 'DE'
         else if (abs(this%lon_corners(1)) < 1000*epsilon(1.0)) then
            this%dateline = 'GE'
         else ! assume 'XY'
            this%dateline = 'XY'
            this%lon_range = RealMinMax(this%lon_centers(1), this%lon_centers(jm))
         end if

         this%lat_corners(1) = this%lat_centers(1) - (this%lat_centers(2)-this%lat_centers(1))/2
         this%lat_corners(2:jm) = (this%lat_centers(1:jm-1) + this%lat_centers(2:jm))/2
         this%lat_corners(jm+1) = this%lat_centers(jm) - (this%lat_centers(jm-1)-this%lat_centers(jm))/2

         if (abs(this%lat_centers(1) + 90) < 1000*epsilon(1.0)) then
            this%pole = 'PC'
         else if (abs(this%lat_corners(1) + 90) < 1000*epsilon(1.0)) then
            this%pole = 'PE'
         else ! assume XY
            this%pole = 'XY'
            this%lat_range = RealMinMax(this%lat_centers(1), this%lat_centers(jm))
         end if
         if (this%lat_corners(1) < -90) this%lat_corners(1)=-90
         if (this%lat_corners(jm+1) > 90) this%lat_corners(jm+1)=90

         ! check if evenly spaced
         regLon=.true.
         do i=2,size(this%lon_centers)
            del12=this%lon_centers(2)-this%lon_centers(1)
            delij=this%lon_centers(i)-this%lon_centers(i-1)
            if ((del12-delij)>epsilon(1.0)) regLon=.false.
         end do
         regLat=.true.
         do i=2,size(this%lat_centers)
            del12=this%lat_centers(2)-this%lat_centers(1)
            delij=this%lat_centers(i)-this%lat_centers(i-1)
            if ((del12-delij)>epsilon(1.0)) regLat=.false.
         end do
         this%is_regular = (regLat .and. regLon)

         if (use_file_coords) then
            this%is_regular = .false.
            this%lon_centers = MAPL_DEGREES_TO_RADIANS_R8 * this%lon_centers
            this%lat_centers = MAPL_DEGREES_TO_RADIANS_R8 * this%lat_centers
            this%lon_corners = MAPL_DEGREES_TO_RADIANS_R8 * this%lon_corners
            this%lat_corners = MAPL_DEGREES_TO_RADIANS_R8 * this%lat_corners
         else
            compute_lons=.false.
            compute_lats=.false.
            if (regLon .and. (this%dateline.ne.'XY')) then
               compute_lons=.true.
            end if
            if (regLat .and. (this%pole.ne.'XY')) then
               compute_lats=.true.
            end if
            if (compute_lons .and. compute_lats) then
               this%lon_centers = this%compute_lon_centers(this%dateline, rc=status)
               _VERIFY(status)
               this%lon_centers_degrees = this%compute_lon_centers(this%dateline, &
                      convert_to_radians=.false., rc=status)
               _VERIFY(status)
               this%lon_corners = this%compute_lon_corners(this%dateline, rc=status)
               _VERIFY(status)
               this%lat_centers_degrees = this%compute_lat_centers(this%pole, &
                      convert_to_radians=.false., rc=status)
               this%lat_centers = this%compute_lat_centers(this%pole, rc=status)
               _VERIFY(status)
               this%lat_corners = this%compute_lat_corners(this%pole, rc=status)
               _VERIFY(status)
            else
               this%lon_centers_degrees = this%lon_centers
               this%lat_centers_degrees = this%lat_centers
               this%lon_centers = MAPL_DEGREES_TO_RADIANS_R8 * this%lon_centers
               this%lat_centers = MAPL_DEGREES_TO_RADIANS_R8 * this%lat_centers
               this%lon_corners = MAPL_DEGREES_TO_RADIANS_R8 * this%lon_corners
               this%lat_corners = MAPL_DEGREES_TO_RADIANS_R8 * this%lat_corners
            end if
         end if

    end associate

    call this%make_arbitrary_decomposition(this%nx, this%ny, rc=status)
    _VERIFY(status)

    ! Determine IMS and JMS with constraint for ESMF that each DE has at least an extent
    ! of 2.  Required for ESMF_FieldRegrid().
    allocate(this%ims(0:this%nx-1))
    allocate(this%jms(0:this%ny-1))
    call MAPL_DecomposeDim(this%im_world, this%ims, this%nx, min_DE_extent=2)
    call MAPL_DecomposeDim(this%jm_world, this%jms, this%ny, min_DE_extent=2)

    call this%check_and_fill_consistency(rc=status)
    _VERIFY(status)

    _RETURN(_SUCCESS)

   end subroutine initialize_from_file_metadata



   subroutine initialize_from_config_with_prefix(this, config, prefix, unusable, rc)
      use esmf
      class (LatLonGridFactory), intent(inout) :: this
      type (ESMF_Config), intent(inout) :: config
      character(len=*), intent(in) :: prefix  ! effectively optional due to overload without this argument
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      integer :: status
      character(len=ESMF_MAXSTR) :: tmp
      type(ESMF_VM) :: VM

      _UNUSED_DUMMY(unusable)

      call ESMF_VmGetCurrent(VM, rc=status)
      _VERIFY(status)

      this%is_regular = .true.
      call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'GRIDNAME:', default=MAPL_GRID_NAME_DEFAULT)
      this%grid_name = trim(tmp)

      call ESMF_ConfigGetAttribute(config, this%nx, label=prefix//'NX:', default=MAPL_UNDEFINED_INTEGER)
      call ESMF_ConfigGetAttribute(config, this%ny, label=prefix//'NY:', default=MAPL_UNDEFINED_INTEGER)

      call ESMF_ConfigGetAttribute(config, this%im_world, label=prefix//'IM_WORLD:', default=MAPL_UNDEFINED_INTEGER)
      call ESMF_ConfigGetAttribute(config, this%jm_world, label=prefix//'JM_WORLD:', default=MAPL_UNDEFINED_INTEGER)

      call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'IMS_FILE:', rc=status)
      if ( status == _SUCCESS ) then
         call get_ims_from_file(this%ims, trim(tmp),this%nx, rc=status)
         _VERIFY(status)
      else
         call get_multi_integer(this%ims, 'IMS:', rc=status)
         _VERIFY(status)
      endif
      call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'JMS_FILE:', rc=status)
      if ( status == _SUCCESS ) then
         call get_ims_from_file(this%jms, trim(tmp),this%ny, rc=status)
         _VERIFY(status)
      else
         call get_multi_integer(this%jms, 'JMS:', rc=status)
         _VERIFY(status)
      endif

      call ESMF_ConfigGetAttribute(config, this%lm, label=prefix//'LM:', default=MAPL_UNDEFINED_INTEGER)

      call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'POLE:', default=MAPL_UNDEFINED_CHAR, rc=status)
      if (status == _SUCCESS) then
         this%pole = trim(tmp)
      end if
      call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'DATELINE:', default=MAPL_UNDEFINED_CHAR, rc=status)
      if (status == _SUCCESS) then
         this%dateline = trim(tmp)
      end if

      call get_range(this%lon_range, 'LON_RANGE:', rc=status); _VERIFY(status)
      call get_range(this%lat_range, 'LAT_RANGE:', rc=status); _VERIFY(status)
      call this%check_and_fill_consistency(rc=status); _VERIFY(status)

      ! Compute the centers and corners
      this%lon_centers = this%compute_lon_centers(this%dateline, rc=status)
      _VERIFY(status)
      this%lon_centers_degrees = this%compute_lon_centers(this%dateline, &
               convert_to_radians = .false., rc=status)
      _VERIFY(status)
      this%lat_centers = this%compute_lat_centers(this%pole, rc=status)
      _VERIFY(status)
      this%lat_centers_degrees = this%compute_lat_centers(this%pole, &
               convert_to_radians = .false., rc=status)
      this%lon_corners = this%compute_lon_corners(this%dateline, rc=status)
      _VERIFY(status)
      this%lat_corners = this%compute_lat_corners(this%pole, rc=status)
      _VERIFY(status)

      _RETURN(_SUCCESS)

   contains

      subroutine get_multi_integer(values, label, rc)
         integer, allocatable, intent(out) :: values(:)
         character(len=*) :: label
         integer, optional, intent(out) :: rc

         integer :: i
         integer :: n
         integer :: tmp
         integer :: status
         logical :: isPresent

         call ESMF_ConfigFindLabel(config, label=prefix//label, isPresent=isPresent, rc=status)
         _VERIFY(status)
         if (.not. isPresent) then
            _RETURN(_SUCCESS)
         end if

         ! First pass:  count values
         n = 0
         do
            call ESMF_ConfigGetAttribute(config, tmp, rc=status)
            if (status /= _SUCCESS) then
               exit
            else
               n = n + 1
            end if
         end do

         ! Second pass: allocate and fill
         allocate(values(n), stat=status) ! no point in checking status
         _VERIFY(status)
         call ESMF_ConfigFindLabel(config, label=prefix//label,rc=status)
         _VERIFY(status)
         do i = 1, n
            call ESMF_ConfigGetAttribute(config, values(i), rc=status)
            _VERIFY(status)
         end do

         _RETURN(_SUCCESS)

      end subroutine get_multi_integer

      subroutine get_ims_from_file(values, file_name, n, rc)
         integer, allocatable, intent(out) :: values(:)
         character(len=*), intent(in) :: file_name
         integer, intent(in) :: n
         integer, optional, intent(out) :: rc

         logical :: FileExists
         integer :: i, total, unit
         integer :: status

         inquire(FILE = trim(file_name), EXIST=FileExists)
         allocate(values(n), stat=status) ! no point in checking status
         _VERIFY(status)

         if ( .not. FileExists) then
             print*, file_name // "   not found"
             _RETURN(_FAILURE)

         elseif (MAPL_AM_I_Root(VM)) then

            open(newunit=UNIT, file=trim(file_name), form="formatted", iostat=status )
            _VERIFY(STATUS)
            read(UNIT,*) total
            if (total /= n) then
                print*, file_name // " n is different from ", total
                _RETURN(_FAILURE)
            endif
            do i = 1,total
                read(UNIT,*) values(i)
            enddo
            close(UNIT)
         endif

         call MAPL_CommsBcast(VM, values, n=N, ROOT=MAPL_Root, rc=status)
         _VERIFY(STATUS)
         _RETURN(_SUCCESS)

      end subroutine get_ims_from_file

      subroutine get_range(range, label, rc)
         type(RealMinMax), intent(out) :: range
         character(len=*) :: label
         integer, optional, intent(out) :: rc

         integer :: status
         logical :: isPresent

         call ESMF_ConfigFindLabel(config, label=prefix//label,isPresent=isPresent,rc=status)
         _VERIFY(status)
         if (.not. isPresent) then
            _RETURN(_SUCCESS)
         end if

         ! Must be 2 values: min and max
         call ESMF_ConfigGetAttribute(config, range%min, rc=status)
         _VERIFY(status)
         call ESMF_ConfigGetAttribute(config, range%max, rc=status)
         _VERIFY(status)

         _RETURN(_SUCCESS)

      end subroutine get_range


   end subroutine initialize_from_config_with_prefix



   function to_string(this) result(string)
      character(len=:), allocatable :: string
      class (LatLonGridFactory), intent(in) :: this

      _UNUSED_DUMMY(this)
      string = 'LatLonGridFactory'

   end function to_string



   subroutine check_and_fill_consistency(this, unusable, rc)
      use MAPL_BaseMod, only: MAPL_DecomposeDim
      class (LatLonGridFactory), intent(inout) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      integer :: status
      logical :: verify_decomp

      _UNUSED_DUMMY(unusable)

      if (.not. allocated(this%grid_name)) then
         this%grid_name = MAPL_GRID_NAME_DEFAULT
      end if

      ! Check decomposition/bounds
      ! WY notes: should not have this assert
      !_ASSERT(allocated(this%ims) .eqv. allocated(this%jms), 'inconsistent options')
      call verify(this%nx, this%im_world, this%ims, rc=status)
      call verify(this%ny, this%jm_world, this%jms, rc=status)

      ! Check regional vs global
      if (this%pole == 'XY') then ! regional
         this%periodic = .false.
         _ASSERT(this%lat_range%min /= MAPL_UNDEFINED_REAL, 'uninitialized min for lat_range')
         _ASSERT(this%lat_range%max /= MAPL_UNDEFINED_REAL, 'uninitialized min for lat_range')
      else ! global
         _ASSERT(any(this%pole == ['PE', 'PC']), 'unsupported option for pole:'//this%pole)
         _ASSERT(this%lat_range%min == MAPL_UNDEFINED_REAL, 'inconsistent min for lat_range')
         _ASSERT(this%lat_range%max == MAPL_UNDEFINED_REAL, 'inconsistent max for lat_range')
      end if
      if (this%dateline == 'XY') then
         this%periodic = .false.
         _ASSERT(this%lon_range%min /= MAPL_UNDEFINED_REAL, 'uninitialized min for lon_range')
         _ASSERT(this%lon_range%max /= MAPL_UNDEFINED_REAL, 'uninitialized max for lon_range')
      else
         _ASSERT(any(this%dateline == ['DC', 'DE', 'GC', 'GE']), 'unsupported option for dateline')
         _ASSERT(this%lon_range%min == MAPL_UNDEFINED_REAL, 'inconsistent min for lon_range')
         _ASSERT(this%lon_range%max == MAPL_UNDEFINED_REAL, 'inconsistent max for lon_range')
      end if
      if (.not.this%force_decomposition) then
         verify_decomp = this%check_decomposition(rc=status)
         _VERIFY(status)
         if ( (.not.verify_decomp) ) then
            call this%generate_newnxy(rc=status)
            _VERIFY(status)
         end if
      end if

      _RETURN(_SUCCESS)

   contains

      subroutine verify(n, m_world, ms, rc)
         integer, intent(inout) :: n
         integer, intent(inout) :: m_world
         integer, allocatable, intent(inout) :: ms(:)
         integer, optional, intent(out) :: rc

         integer :: status

         if (allocated(ms)) then
            _ASSERT(size(ms) > 0, 'degenerate topology')

            if (n == MAPL_UNDEFINED_INTEGER) then
               n = size(ms)
            else
               _ASSERT(n == size(ms), 'inconsistent topology')
            end if

            if (m_world == MAPL_UNDEFINED_INTEGER) then
               m_world = sum(ms)
            else
               _ASSERT(m_world == sum(ms), 'inconsistent decomponsition')
            end if

         else

            _ASSERT(n /= MAPL_UNDEFINED_INTEGER, 'uninitialized topology')
            _ASSERT(m_world /= MAPL_UNDEFINED_INTEGER,'uninitialized dimension')
            allocate(ms(n), stat=status)
            _VERIFY(status)
            !call MAPL_DecomposeDim(m_world, ms, n, min_DE_extent=2)
            call MAPL_DecomposeDim(m_world, ms, n)

         end if

         _RETURN(_SUCCESS)

      end subroutine verify

   end subroutine check_and_fill_consistency


   elemental subroutine set_with_default_integer(to, from, default)
      integer, intent(out) :: to
      integer, optional, intent(in) :: from
      integer, intent(in) :: default

      if (present(from)) then
         to = from
      else
         to = default
      end if

   end subroutine set_with_default_integer


   elemental subroutine set_with_default_real(to, from, default)
      real, intent(out) :: to
      real, optional, intent(in) :: from
      real, intent(in) :: default

      if (present(from)) then
         to = from
      else
         to = default
      end if

   end subroutine set_with_default_real

   subroutine set_with_default_character(to, from, default)
      character(len=:), allocatable, intent(out) :: to
      character(len=*), optional, intent(in) :: from
      character(len=*), intent(in) :: default

      if (present(from)) then
         to = from
      else
         to = default
      end if

   end subroutine set_with_default_character


   elemental subroutine set_with_default_range(to, from, default)
      type (RealMinMax), intent(out) :: to
      type (RealMinMax), optional, intent(in) :: from
      type (RealMinMax), intent(in) :: default

      if (present(from)) then
         to = from
      else
         to = default
      end if

   end subroutine set_with_default_range

   subroutine set_with_default_logical(to, from, default)
      logical, intent(out) :: to
      logical, optional, intent(in) :: from
      logical, intent(in) :: default

      if (present(from)) then
         to = from
      else
         to = default
      end if

   end subroutine set_with_default_logical

   ! MAPL uses values in lon_array and lat_array only to determine the
   ! general positioning.  Actual coordinates are then recomputed.
   ! This helps to avoid roundoff differences from slightly different
   ! input files.
   subroutine initialize_from_esmf_distGrid(this, dist_grid, lon_array, lat_array, unusable, rc)
      use MAPL_ConfigMod
      use MAPL_Constants, only: PI => MAPL_PI_R8
      class (LatLonGridFactory), intent(inout)  :: this
      type (ESMF_DistGrid), intent(in) :: dist_grid
      type (ESMF_LocalArray), intent(in) :: lon_array
      type (ESMF_LocalArray), intent(in) :: lat_array
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      integer :: dim_count, tile_count
      integer, allocatable :: max_index(:,:)
      integer :: status
      character(len=2) :: pole ,dateline

      type (ESMF_Config) :: config
      type (ESMF_VM) :: vm
      integer :: nPet
      real(kind=REAL32), pointer :: lon(:)
      real(kind=REAL32), pointer :: lat(:)
      integer :: nx_guess,nx,ny
      integer :: i

      real, parameter :: tiny = 1.e-4

      _UNUSED_DUMMY(unusable)

      this%is_regular = .true.
      call ESMF_DistGridGet(dist_grid, dimCount=dim_count, tileCount=tile_count)
      allocate(max_index(dim_count, tile_count))
      call ESMF_DistGridGet(dist_grid, maxindexPTile=max_index)

      config = MAPL_ConfigCreate(rc=status)
      _VERIFY(status)
      call MAPL_ConfigSetAttribute(config, max_index(1,1), 'IM_WORLD:', rc=status)
      _VERIFY(status)
      call MAPL_ConfigSetAttribute(config, max_index(2,1), 'JM_WORLD:', rc=status)
      _VERIFY(status)
      call MAPL_ConfigSetAttribute(config, max_index(3,1), 'LM:', rc=status)
      _VERIFY(status)

      lon => null()
      lat => null()
      call ESMF_LocalArrayGet(lon_array, farrayPtr=lon, rc=status)
      _VERIFY(status)
      call ESMF_LocalArrayGet(lat_array, farrayPtr=lat, rc=status)
      _VERIFY(status)


      if (abs(lat(1) + PI/2) < tiny) then
         pole = 'PC'
      elseif (abs(lat(1) + PI/2 - 0.5*(lat(2)-lat(1))) < tiny) then
         pole = 'PE'
      else
         pole = 'PC'
      end if

      ! the code below is kluge to return DE/DC wheither or not the file lons are -180 to 180 or 0 360
      ! it detects whether the first longitudes which are cell centers
      ! If first longitude is 0 or -180 (DC) it is dateline center in that 0 or -180 is
      ! in the center of a grid cell.
      ! or shifted by half a grid box (DE) so 0 or -180 is the edge of a cell
      ! really should have 4 options dateline edge (DE), dateline center(DC)
      ! grenwich center (GC) and grenwich edge (GE) but the last 2 are not supported
      ! if it is GC or GE we will shift the data on the usage so that it is DE or DC for now
      do i=0,1
         if (abs(lon(1) + PI*i) < tiny) then
            dateline = 'DC'
            exit
         elseif (abs(lon(1) + PI*i - 0.5*(lon(2)-lon(1))) < tiny) then
            dateline = 'DE'
            exit
         end if
      end do
      !if (abs(lon(1) + PI) < tiny) then
      !dateline = 'DC'
      !elseif (abs(lon(1) + PI - 0.5*(lon(2)-lon(1))) < tiny) then
      !dateline = 'DE'
      !elseif (abs(lon(1)) < tiny) then
      !dateline = 'GC'
      !elseif (abs(lon(1) - 0.5*(lon(2)-lon(1))) < tiny) then
      !dateline = 'GE'
      !end if

      call MAPL_ConfigSetAttribute(config, pole, 'POLE:')
      call MAPL_ConfigSetAttribute(config, dateline, 'DATELINE:')

      call ESMF_VMGetCurrent(vm, rc=status)
      _VERIFY(status)
      call ESMF_VMGet(vm, PETcount=nPet, rc=status)
      _VERIFY(status)

      nx_guess = nint(sqrt(real(nPet)))
      do nx = nx_guess,1,-1
         ny=nPet/nx
         if (nx*ny==nPet) then
            call MAPL_ConfigSetAttribute(config, nx, 'NX:')
            call MAPL_ConfigSetAttribute(config, ny, 'NY:')
            exit
         end if
      enddo

      call this%initialize(config, rc=status)
      _VERIFY(status)


   end subroutine initialize_from_esmf_distGrid

   function decomps_are_equal(this,a) result(equal)
      class (LatLonGridFactory), intent(in) :: this
      class (AbstractGridFactory), intent(in) :: a
      logical :: equal

      select type (a)
         class default
         equal = .false.
         return
      class is (LatLonGridFactory)
         equal = .true.


         equal = size(a%ims)==size(this%ims) .and. size(a%jms)==size(this%jms)
         if (.not. equal) return

         ! same decomposition
         equal = all(a%ims == this%ims) .and. all(a%jms == this%jms)
         if (.not. equal) return

      end select

   end function decomps_are_equal


   function physical_params_are_equal(this, a) result(equal)
      class (LatLonGridFactory), intent(in) :: this
      class (AbstractGridFactory), intent(in) :: a
      logical :: equal

      select type (a)
         class default
         equal = .false.
         return
      class is (LatLonGridFactory)
         equal = .true.

         equal = (a%im_world == this%im_world) .and. (a%jm_world == this%jm_world)
         if (.not. equal) return

         equal = (a%is_regular .eqv. this%is_regular)
         if (.not. equal) return

         if (a%is_regular) then
            equal = (a%pole == this%pole)
            if (.not. equal) return

            equal = (a%dateline == this%dateline)
            if (.not. equal) return

            if (a%pole == 'XY') then
               equal = (a%lat_range == this%lat_range)
               if (.not. equal) return
            end if

            if (a%dateline == 'XY') then
               equal = (a%lon_range == this%lon_range)
               if (.not. equal) return
            end if
         else
            equal = &
                 & all(a%lon_centers == this%lon_centers) .and. &
                 & all(a%lon_corners == this%lon_corners) .and. &
                 & all(a%lat_centers == this%lat_centers) .and. &
                 & all(a%lat_corners == this%lat_corners)
         end if
      end select

   end function physical_params_are_equal

   logical function equals(a, b)
      class (LatLonGridFactory), intent(in) :: a
      class (AbstractGridFactory), intent(in) :: b

      select type (b)
         class default
         equals = .false.
         return
      class is (LatLonGridFactory)
         equals = .true.

         equals = (a%lm == b%lm)
         if (.not. equals) return

         equals = a%decomps_are_equal(b)
         if (.not. equals) return

         equals = a%physical_params_are_equal(b)
         if (.not. equals) return

      end select

   end function equals


   function generate_grid_name(this) result(name)
      character(len=:), allocatable :: name
      class (LatLonGridFactory), intent(in) :: this

      character(len=4) :: im_string, jm_string

      write(im_string,'(i4.4)') this%im_world
      write(jm_string,'(i4.4)') this%jm_world

      name = this%dateline // im_string // 'x' // this%pole // jm_string

   end function generate_grid_name

   function check_decomposition(this,unusable,rc) result(can_decomp)
      class (LatLonGridFactory), target, intent(inout) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc
      logical :: can_decomp
      integer :: n
      _UNUSED_DUMMY(unusable)

      can_decomp = .true.
      if (this%im_world==1 .and. this%jm_world==1) then
         _RETURN(_SUCCESS)
      end if
      n = this%im_world/this%nx
      if (n < 2) can_decomp = .false.
      n = this%jm_world/this%ny
      if (n < 2) can_decomp = .false.
      _RETURN(_SUCCESS)
   end function check_decomposition

   subroutine generate_newnxy(this,unusable,rc)
      use MAPL_BaseMod, only: MAPL_DecomposeDim
      class (LatLonGridFactory), target, intent(inout) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc
      integer :: n

      _UNUSED_DUMMY(unusable)

      n = this%im_world/this%nx
      if (n < 2) then
         this%nx = generate_new_decomp(this%im_world,this%nx)
         deallocate(this%ims)
         allocate(this%ims(0:this%nx-1))
         call MAPL_DecomposeDim(this%im_world, this%ims, this%nx)
      end if
      n = this%jm_world/this%ny
      if (n < 2) then
         this%ny = generate_new_decomp(this%jm_world,this%ny)
         deallocate(this%jms)
         allocate(this%jms(0:this%ny-1))
         call MAPL_DecomposeDim(this%jm_world, this%jms, this%ny)
      end if

      _RETURN(_SUCCESS)

   end subroutine generate_newnxy

   function generate_new_decomp(im,nd) result(n)
      integer, intent(in) :: im, nd
      integer :: n
      logical :: canNotDecomp

      canNotDecomp = .true.
      n = nd
      do while(canNotDecomp)
         if ( (im/n) < 2) then
            n = n/2
         else
            canNotDecomp = .false.
         end if
      enddo
   end function generate_new_decomp

   subroutine init_halo(this, unusable, rc)
      class (LatLonGridFactory), target, intent(inout) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      type (ESMF_Grid), pointer :: grid
      integer :: dim_count
      integer :: pet
      integer :: ndes
      type (ESMF_DistGrid) :: dist_grid
      type (ESMF_VM) :: vm

      integer :: status

      _UNUSED_DUMMY(unusable)

      if (this%is_halo_initialized) return

      grid => this%get_grid(rc=status)
      _VERIFY(status)

      call ESMF_GridGet(grid,   distGrid=dist_grid, dimCount=dim_count, rc=status)
      _VERIFY(status)
      call ESMF_DistGridGet(dist_grid, delayout=this%layout, rc=status)
      _VERIFY(status)
      call ESMF_DELayoutGet (this%layout, vm=vm, rc=status)
      _VERIFY(status)

      call ESMF_VMGet(vm, localPet=pet, petCount=ndes, rc=status)
      _VERIFY(status)

      this%px = mod(pet, this%nx)
      this%py = pet / this%nx

      this%is_halo_initialized = .true.

      _RETURN(_SUCCESS)

   end subroutine init_halo


   subroutine halo(this, array, unusable, halo_width, rc)
      use MAPL_CommsMod
      class (LatLonGridFactory), intent(inout) :: this
      real(kind=REAL32), intent(inout) :: array(:,:)
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(in) :: halo_width
      integer, optional, intent(out) :: rc

      integer :: status

      integer :: pet_north
      integer :: pet_south
      integer :: pet_east
      integer :: pet_west
      integer :: pet_N_E, pet_N_W, pet_S_E, pet_S_W
      integer :: last_lon, last_lat

      _UNUSED_DUMMY(unusable)
      ! not yet implmented, default is 1
      _UNUSED_DUMMY(halo_width)

      if (.not. this%is_halo_initialized) then
         call this%init_halo(rc=status)
         _VERIFY(status)
      end if

      last_lon = size(array,1)
      last_lat = size(array,2)

      associate (nx => this%nx, ny => this% ny, px => this%px, py => this%py)
        ! Nearest neighbors processor' ids
        pet_north = get_pet(px, py+1, nx, ny)
        pet_south = get_pet(px, py-1, nx, ny)
        pet_east  = get_pet(px+1, py, nx, ny)
        pet_west  = get_pet(px-1, py, nx, ny)

        call fill_north(array, rc=status)
        _VERIFY(status)
        call fill_south(array, rc=status)
        _VERIFY(status)

        call fill_east(array, rc=status)
        _VERIFY(status)
        call fill_west(array, rc=status)
        _VERIFY(status)

        pet_N_E   = get_pet(px+1, py+1, nx, ny)
        pet_N_W   = get_pet(px-1, py+1, nx, ny)
        pet_S_E   = get_pet(px+1, py-1, nx, ny)
        pet_S_W   = get_pet(px-1, py-1, nx, ny)

        !fill north east
        call MAPL_CommsSendRecv(this%layout,              &
              array(2,      2         ), 1,  pet_S_W,  &
              array(last_lon,last_lat ), 1,  pet_N_E,  &
              rc=status)
        _VERIFY(status)

        !fill north west
        call MAPL_CommsSendRecv(this%layout,              &
              array(last_lon-1, 2), 1,  pet_S_E,  &
              array(1,   last_lat), 1,  pet_N_W,  &
              rc=status)
        _VERIFY(status)

        ! north pol corner
        if(this%py== this%ny-1) then
           array(last_lon,last_lat) = array(last_lon-1,last_lat-1)
           array(1,last_lat)        = array(2,last_lat-1)
        endif

        !fill south east
        call MAPL_CommsSendRecv(this%layout,              &
              array(2, last_lat-1), 1,  pet_N_W,  &
              array(last_lon,1),    1,  pet_S_E,  &
              rc=status)
        _VERIFY(status)

        !fill south west
        call MAPL_CommsSendRecv(this%layout,              &
              array(last_lon-1,last_lat-1), 1,  pet_N_E,  &
              array(1,1                  ), 1,  pet_S_W,  &
              rc=status)
        _VERIFY(status)

        ! south pole corner
        if(this%py==0) then
           array(last_lon,1   ) = array(last_lon-1,2 )
           array(1,1   )        = array(2,2 )
        endif

      end associate

      _RETURN(ESMF_SUCCESS)

   contains

      integer function get_pet(px, py, nx, ny) result(pet)
         integer, intent(in) :: px, py  ! rank in x/y directions
         integer, intent(in) :: nx, ny  ! npets in x/y directions

         pet = mod(px+nx,nx) + nx*mod(py+ny,ny)

      end function get_pet


      subroutine fill_north(array, rc)
         real(kind=REAL32), intent(inout) :: array(:,:)
         integer, optional, intent(out) :: rc

         integer :: status

         integer :: len, last

         last = size(array,2)-1
         len = size(array,1)

         call MAPL_CommsSendRecv(this%layout,        &
              array(:,2        ),  len,  pet_south,  &
              array(:,last+1   ),  len,  pet_north,  &
              rc=status)
         _VERIFY(status)
         if(this%py==this%ny-1) then
            array(:,last+1   ) = array(:,last )
         end if

         _RETURN(_SUCCESS)

      end subroutine fill_north


      subroutine fill_south(array, rc)
         real(kind=REAL32), intent(inout) :: array(:,:)
         integer, optional, intent(out) :: rc

         integer :: status

         integer :: len, last

         last = size(array,2)-1
         len = size(array,1)

         call MAPL_CommsSendRecv(this%layout,     &
              array(:,last     ),  len,  pet_north,  &
              array(:,1        ),  len,  pet_south,  &
              rc=status)
         _VERIFY(status)

         if(this%py==0) then
            array(:,1   ) = array(:,2 )
         endif

         _RETURN(_SUCCESS)

      end subroutine fill_south


      subroutine fill_east(array, rc)
         real(kind=REAL32), intent(inout) :: array(:,:)
         integer, optional, intent(out) :: rc

         integer :: status

         integer :: len, last

         last = size(array,1)-1
         len = size(array,2)

         call MAPL_CommsSendRecv(this%layout,      &
              array(2     , : ),  len,  pet_west,  &
              array(last+1, : ),  len,  pet_east,  &
              rc=status)
         _VERIFY(status)

         _RETURN(_SUCCESS)

      end subroutine fill_east


      subroutine fill_west(array, rc)
         real(kind=REAL32), intent(inout) :: array(:,:)
         integer, optional, intent(out) :: rc

         integer :: status

         integer :: len, last

         last = size(array,1)-1
         len = size(array,2)

         call MAPL_CommsSendRecv(this%layout,   &
              array(last  , : ),  len,  pet_east,  &
              array(1     , : ),  len,  pet_west,  &
              rc=status)
         _VERIFY(status)

         _RETURN(_SUCCESS)

      end subroutine fill_west


   end subroutine halo


   subroutine append_metadata(this, metadata)
      use MAPL_Constants
      class (LatLonGridFactory), intent(inout) :: this
      type (FileMetadata), intent(inout) :: metadata

      type (Variable) :: v
      real(kind=REAL64), allocatable :: temp_coords(:)

      ! Horizontal grid dimensions
      call metadata%add_dimension('lon', this%im_world)
      call metadata%add_dimension('lat', this%jm_world)

      ! Coordinate variables
      v = Variable(type=PFIO_REAL64, dimensions='lon')
      call v%add_attribute('long_name', 'longitude')
      call v%add_attribute('units', 'degrees_east')
      temp_coords = this%get_longitudes_degrees()
      call v%add_const_value(UnlimitedEntity(temp_coords))
      call metadata%add_variable('lon', v)
      deallocate(temp_coords)

      v = Variable(type=PFIO_REAL64, dimensions='lat')
      call v%add_attribute('long_name', 'latitude')
      call v%add_attribute('units', 'degrees_north')
      temp_coords=this%get_latitudes_degrees()
      call v%add_const_value(UnlimitedEntity(temp_coords))
      call metadata%add_variable('lat', v)

   end subroutine append_metadata

   function get_grid_vars(this) result(vars)
      class (LatLonGridFactory), intent(inout) :: this

      character(len=:), allocatable :: vars
      _UNUSED_DUMMY(this)

      vars = 'lon,lat'

   end function get_grid_vars

   function get_file_format_vars(this) result(vars)
      class (LatLonGridFactory), intent(inout) :: this

      character(len=:), allocatable :: vars
      _UNUSED_DUMMY(this)

      vars = 'lon,lat'

   end function get_file_format_vars

   subroutine append_variable_metadata(this,var)
      class (LatLonGridFactory), intent(inout) :: this
      type(Variable), intent(inout) :: var
      _UNUSED_DUMMY(this)
      _UNUSED_DUMMY(var)
   end subroutine append_variable_metadata

   subroutine generate_file_bounds(this,grid,local_start,global_start,global_count,metadata,rc)
      use MAPL_BaseMod
      class(LatLonGridFactory), intent(inout) :: this
      type(ESMF_Grid),      intent(inout) :: grid
      integer, allocatable, intent(out) :: local_start(:)
      integer, allocatable, intent(out) :: global_start(:)
      integer, allocatable, intent(out) :: global_count(:)
      type(FileMetaData), intent(in), optional :: metaData
      integer, optional, intent(out) :: rc

      integer :: status
      integer :: global_dim(3), i1,j1,in,jn

      call MAPL_GridGet(grid,globalCellCountPerDim=global_dim,rc=status)
      _VERIFY(status)
      call MAPL_GridGetInterior(grid,i1,in,j1,jn)
      allocate(local_start,source=[i1,j1])
      allocate(global_start,source=[1,1])
      allocate(global_count,source=[global_dim(1),global_dim(2)])

      _RETURN(_SUCCESS)

      _UNUSED_DUMMY(this)
      _UNUSED_DUMMY(metadata)
   end subroutine generate_file_bounds

   subroutine generate_file_corner_bounds(this,grid,local_start,global_start,global_count,rc)
      use esmf
      class (LatLonGridFactory), intent(inout) :: this
      type(ESMF_Grid), intent(inout)      :: grid
      integer, allocatable, intent(out) :: local_start(:)
      integer, allocatable, intent(out) :: global_start(:)
      integer, allocatable, intent(out) :: global_count(:)
      integer, optional, intent(out) :: rc

      _UNUSED_DUMMY(this)
      _UNUSED_DUMMY(grid)
      _UNUSED_DUMMY(local_start)
      _UNUSED_DUMMY(global_start)
      _UNUSED_DUMMY(global_count)

      _FAIL('unimplemented')
      _RETURN(_SUCCESS)
   end subroutine generate_file_corner_bounds

   function generate_file_reference2D(this,fpointer) result(ref)
      use pFIO
      type(ArrayReference) :: ref
      class(LatLonGridFactory), intent(inout) :: this
      real, pointer, intent(in) :: fpointer(:,:)
      ref = ArrayReference(fpointer)
      _UNUSED_DUMMY(this)
   end function generate_file_reference2D

   function generate_file_reference3D(this,fpointer,metaData) result(ref)
      use pFIO
      type(ArrayReference) :: ref
      class(LatLonGridFactory), intent(inout) :: this
      real, pointer, intent(in) :: fpointer(:,:,:)
      type(FileMetaData), intent(in), optional :: metaData
      ref = ArrayReference(fpointer)

      _UNUSED_DUMMY(this)
      _UNUSED_DUMMY(metaData)
   end function generate_file_reference3D


end module MAPL_LatLonGridFactoryMod