MAPL_CubedSphereGridFactory.F90 Source File


This file depends on

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

Files dependent on this one

sourcefile~~mapl_cubedspheregridfactory.f90~~AfferentGraph sourcefile~mapl_cubedspheregridfactory.f90 MAPL_CubedSphereGridFactory.F90 sourcefile~base.f90 Base.F90 sourcefile~base.f90->sourcefile~mapl_cubedspheregridfactory.f90 sourcefile~mapl_gridmanager.f90 MAPL_GridManager.F90 sourcefile~base.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_cubedspheregridfactory.f90 sourcefile~regrid_util.f90 Regrid_Util.F90 sourcefile~regrid_util.f90->sourcefile~mapl_cubedspheregridfactory.f90 sourcefile~regrid_util.f90->sourcefile~mapl_gridmanager.f90 sourcefile~cfiocollection.f90 CFIOCollection.F90 sourcefile~cfiocollection.f90->sourcefile~mapl_gridmanager.f90 sourcefile~comp_testing_driver.f90 Comp_Testing_Driver.F90 sourcefile~comp_testing_driver.f90->sourcefile~mapl_gridmanager.f90 sourcefile~datacollection.f90 DataCollection.F90 sourcefile~datacollection.f90->sourcefile~mapl_gridmanager.f90 sourcefile~esmfl_mod.f90 ESMFL_Mod.F90 sourcefile~esmfl_mod.f90->sourcefile~mapl_gridmanager.f90 sourcefile~extdatagridcompmod.f90 ExtDataGridCompMod.F90 sourcefile~extdatagridcompmod.f90->sourcefile~mapl_gridmanager.f90 sourcefile~extdatagridcompng.f90 ExtDataGridCompNG.F90 sourcefile~extdatagridcompng.f90->sourcefile~mapl_gridmanager.f90 sourcefile~fieldbundleread.f90 FieldBundleRead.F90 sourcefile~fieldbundleread.f90->sourcefile~mapl_gridmanager.f90 sourcefile~filemetadatautilities.f90 FileMetadataUtilities.F90 sourcefile~filemetadatautilities.f90->sourcefile~mapl_gridmanager.f90 sourcefile~griddedio.f90 GriddedIO.F90 sourcefile~griddedio.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl.f90 MAPL.F90 sourcefile~mapl.f90->sourcefile~base.f90 sourcefile~mapl_bundleio_test.f90 mapl_bundleio_test.F90 sourcefile~mapl_bundleio_test.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_capgridcomp.f90 MAPL_CapGridComp.F90 sourcefile~mapl_capgridcomp.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_cfio.f90 MAPL_CFIO.F90 sourcefile~mapl_cfio.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_epochswathmod.f90 MAPL_EpochSwathMod.F90 sourcefile~mapl_epochswathmod.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_esmfregridder.f90 MAPL_EsmfRegridder.F90 sourcefile~mapl_esmfregridder.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_generic.f90 MAPL_Generic.F90 sourcefile~mapl_generic.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_historycollection.f90 MAPL_HistoryCollection.F90 sourcefile~mapl_historycollection.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_historygridcomp.f90 MAPL_HistoryGridComp.F90 sourcefile~mapl_historygridcomp.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_nuopcwrappermod.f90 MAPL_NUOPCWrapperMod.F90 sourcefile~mapl_nuopcwrappermod.f90->sourcefile~base.f90 sourcefile~mapl_regriddermanager.f90 MAPL_RegridderManager.F90 sourcefile~mapl_regriddermanager.f90->sourcefile~mapl_gridmanager.f90 sourcefile~mapl_tilingregridder.f90 MAPL_TilingRegridder.F90 sourcefile~mapl_tilingregridder.f90->sourcefile~mapl_gridmanager.f90 sourcefile~newregriddermanager.f90 NewRegridderManager.F90 sourcefile~newregriddermanager.f90->sourcefile~mapl_gridmanager.f90 sourcefile~regridderspec.f90 RegridderSpec.F90 sourcefile~regridderspec.f90->sourcefile~mapl_gridmanager.f90 sourcefile~test_cfio_bundle.pf Test_CFIO_Bundle.pf sourcefile~test_cfio_bundle.pf->sourcefile~base.f90 sourcefile~test_gridmanager.pf Test_GridManager.pf sourcefile~test_gridmanager.pf->sourcefile~mapl_gridmanager.f90 sourcefile~tstqsat.f90 tstqsat.F90 sourcefile~tstqsat.f90->sourcefile~base.f90 sourcefile~ut_extdata.f90 ut_ExtData.F90 sourcefile~ut_extdata.f90->sourcefile~base.f90 sourcefile~utcfio_bundle.f90 utCFIO_Bundle.F90 sourcefile~utcfio_bundle.f90->sourcefile~base.f90

Source Code

!-----------------------------------------------------
! Note that this implementation only supports
! "square" faces on the cube.   I.e. the number of
! cells along each axis (of each face) are the same.
! IM_WORLD is used for this quantity, and there is no
! equivalent for the "other" axis.
!-----------------------------------------------------

#include "MAPL_Generic.h"


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

   public :: CubedSphereGridFactory

   integer, parameter :: ndims = 2

   integer, parameter :: FV_GRID_TYPE_DEFAULT = 0

   integer, parameter :: NUM_CUBE_FACES = 6

   type, extends(AbstractGridFactory) :: CubedSphereGridFactory
      private


      character(len=:), allocatable :: grid_name
      integer :: grid_type = MAPL_UNDEFINED_INTEGER

      ! Grid dimensions - Note that we only support "square" grids
      integer :: im_world = MAPL_UNDEFINED_INTEGER
      integer :: lm = MAPL_UNDEFINED_INTEGER
      integer :: ntiles = NUM_CUBE_FACES

      ! Domain decomposition: - note that we only support "square" dec
      integer :: nx = MAPL_UNDEFINED_INTEGER
      integer :: ny = MAPL_UNDEFINED_INTEGER
      integer, allocatable :: ims(:)
      integer, allocatable :: jms(:)
      ! rectangle decomposition
      integer, allocatable :: jms_2d(:,:)
      ! stretching parameters
      real :: stretch_factor = MAPL_UNDEFINED_REAL
      real :: target_lon = MAPL_UNDEFINED_REAL
      real :: target_lat = MAPL_UNDEFINED_REAL
      real :: target_lon_degrees = MAPL_UNDEFINED_REAL
      real :: target_lat_degrees = MAPL_UNDEFINED_REAL
      logical :: stretched_cube = .false.

      ! For halo
      type(ESMF_RouteHandle) :: rh

      logical :: halo_initialized = .false.

   contains

      procedure :: make_new_grid
      procedure :: create_basic_grid

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

      procedure :: halo_init
      procedure :: halo

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

      procedure :: append_metadata
      procedure :: get_grid_vars
      procedure :: get_file_format_vars
      procedure :: append_variable_metadata
      procedure :: generate_file_bounds
      procedure :: generate_file_corner_bounds
      procedure :: generate_file_reference2D
      procedure :: generate_file_reference3D
      procedure :: get_fake_longitudes
      procedure :: get_fake_latitudes
      procedure :: decomps_are_equal
      procedure :: physical_params_are_equal
   end type CubedSphereGridFactory

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

   interface CubedSphereGridFactory
      module procedure CubedSphereGridFactory_from_parameters
   end interface CubedSphereGridFactory

   interface set_with_default
      module procedure set_with_default_integer
      module procedure set_with_default_real
      module procedure set_with_default_real64
      module procedure set_with_default_character
      module procedure set_with_default_bounds
   end interface set_with_default


contains


   function CubedSphereGridFactory_from_parameters(unusable, grid_name, grid_type, &
        & im_world, lm, nx, ny, ims, jms, stretch_factor, target_lon, target_lat, &
        & rc) result(factory)
      type (CubedSphereGridFactory) :: factory
      class (KeywordEnforcer), optional, intent(in) :: unusable
      character(len=*), optional, intent(in) :: grid_name
      integer, optional, intent(in) :: grid_type

      ! grid details:
      integer, optional, intent(in) :: im_world
      integer, optional, intent(in) :: lm

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

      ! stretched grid
      real(REAL32), optional, intent(in) :: stretch_factor, target_lon, target_lat

      integer, optional, intent(out) :: rc

      integer :: status
      character(len=*), parameter :: Iam = MOD_NAME // 'CubedSphereGridFactory_from_parameters'

      if (present(unusable)) print*,shape(unusable)

      call set_with_default(factory%grid_name, grid_name, MAPL_GRID_NAME_DEFAULT)
      call set_with_default(factory%grid_type, grid_type, FV_GRID_TYPE_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%lm, lm, MAPL_UNDEFINED_INTEGER)

      call set_with_default(factory%stretch_factor,stretch_factor,MAPL_UNDEFINED_REAL)
      call set_with_default(factory%target_lon,target_lon,MAPL_UNDEFINED_REAL)
      call set_with_default(factory%target_lat,target_lat,MAPL_UNDEFINED_REAL)

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

      call factory%check_and_fill_consistency(rc=status)

      _VERIFY(status)

      _RETURN(_SUCCESS)

   end function CubedSphereGridFactory_from_parameters


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

      integer :: status
      character(len=*), parameter :: Iam = MOD_NAME // 'make_grid'

      _UNUSED_DUMMY(unusable)

      grid = this%create_basic_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 (CubedSphereGridFactory), intent(in) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      integer :: i,nTile
      integer, allocatable :: ims(:,:), jms(:,:)
      real(kind=ESMF_KIND_R8), pointer :: lats(:,:),lons(:,:)
      type(ESMF_CubedSphereTransform_Args) :: transformArgument
      integer :: status
      character(len=*), parameter :: Iam = MOD_NAME // 'create_basic_grid'

      _UNUSED_DUMMY(unusable)

      if (this%grid_type <=3) then
         nTile=6
      else
         nTile=1
      end if

      allocate(ims(this%nx,nTile))
      do i=1,nTile
         ims(:,i)=this%ims
      enddo

      if(allocated(this%jms_2d)) then
         _ASSERT(size(this%jms_2d,2) == 6,'incompatible shape')
         allocate(jms, source = this%jms_2d)
      else
         allocate(jms(this%ny,nTile))
         do i=1,nTile
            jms(:,i)=this%jms
         end do
      endif

      if (this%grid_type <= 3) then
         if (this%stretched_cube) then
            transformArgument%stretch_factor=this%stretch_factor
            transformArgument%target_lon=this%target_lon
            transformArgument%target_lat=this%target_lat
            grid = ESMF_GridCreateCubedSPhere(this%im_world,countsPerDEDim1PTile=ims, &
                      countsPerDEDim2PTile=jms ,name=this%grid_name, &
                      staggerLocList=[ESMF_STAGGERLOC_CENTER,ESMF_STAGGERLOC_CORNER], coordSys=ESMF_COORDSYS_SPH_RAD, &
                      transformArgs=transformArgument,rc=status)
            _VERIFY(status)
            if (this%stretch_factor/=MAPL_UNDEFINED_REAL .and. this%target_lon/=MAPL_UNDEFINED_REAL .and. &
                this%target_lat/=MAPL_UNDEFINED_REAL) then
               call ESMF_AttributeSet(grid, name='STRETCH_FACTOR', value=this%stretch_factor,rc=status)
               _VERIFY(status)
               call ESMF_AttributeSet(grid, name='TARGET_LON', value=this%target_lon_degrees,rc=status)
               _VERIFY(status)
               call ESMF_AttributeSet(grid, name='TARGET_LAT', value=this%target_lat_degrees,rc=status)
               _VERIFY(status)
            end if
         else
            grid = ESMF_GridCreateCubedSPhere(this%im_world,countsPerDEDim1PTile=ims, &
                      countsPerDEDim2PTile=jms ,name=this%grid_name, &
                      staggerLocList=[ESMF_STAGGERLOC_CENTER,ESMF_STAGGERLOC_CORNER], coordSys=ESMF_COORDSYS_SPH_RAD, rc=status)
            _VERIFY(status)
         end if
         call ESMF_AttributeSet(grid, name='GridType', value='Cubed-Sphere', rc=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)
         call ESMF_AttributeSet(grid, 'GridType', 'Doubly-Periodic', rc=status)
         _VERIFY(status)
         call ESMF_GridAddCoord(grid,rc=status)
         _VERIFY(status)
         call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, &
           staggerloc=ESMF_STAGGERLOC_CENTER, &
           farrayPtr=lons, rc=status)
         _VERIFY(status)
         lons=0.0
         call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, &
           staggerloc=ESMF_STAGGERLOC_CENTER, &
           farrayPtr=lats, rc=status)
         _VERIFY(status)
         lats=0.0

      end if

      deallocate(ims,jms)

      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, name='NEW_CUBE', value=1,rc=status)
      _VERIFY(status)

      _RETURN(_SUCCESS)
   end function create_basic_grid

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

      class (CubedSphereGridFactory), 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

      character(len=*), parameter :: Iam= MOD_NAME // 'initialize_from_file_metadata()'
      integer :: status
      logical :: hasLev,hasLevel,hasXdim,hasLon
      logical :: is_stretched
      character(:), allocatable :: lev_name
      type(Attribute), pointer :: attr
      class(*), pointer :: attr_val(:)

      associate(im => this%im_world)
         hasXdim = file_metadata%has_dimension('Xdim')
         hasLon = file_metadata%has_dimension('lon',rc=status)
         if (hasXdim .and. (.not.haslon)) then
            im = file_metadata%get_dimension('Xdim',rc=status)
            _VERIFY(status)
         else if (hasLon .and. (.not.hasXdim)) then
            im = file_metadata%get_dimension('lon',rc=status)
            _VERIFY(status)
         else
            _FAIL("can not identify dimenions of cubed-sphere file")
         end if
      end associate
      call this%make_arbitrary_decomposition(this%nx, this%ny, reduceFactor=6, rc=status)
      _VERIFY(status)

      is_stretched = file_metadata%has_attribute('STRETCH_FACTOR') .and. &
                     file_metadata%has_attribute('TARGET_LON') .and. &
                     file_metadata%has_attribute('TARGET_LAT')
      if (is_stretched) then
         attr => file_metadata%get_attribute('STRETCH_FACTOR')
         attr_val => attr%get_values()
         select type(q=>attr_val)
         type is (real(kind=REAL32))
            this%stretch_factor = q(1)
         class default
            _FAIL('unsupport subclass for stretch params')
         end select
         attr => file_metadata%get_attribute('TARGET_LAT')
         attr_val => attr%get_values()
         select type(q=>attr_val)
         type is (real(kind=REAL32))
            this%target_lat = q(1)
         type is (real(kind=REAL64))
            this%target_lat = q(1)
         class default
            _FAIL('unsupport subclass for stretch params')
         end select
         attr => file_metadata%get_attribute('TARGET_LON')
         attr_val => attr%get_values()
         select type(q=>attr_val)
         type is (real(kind=REAL32))
            this%target_lon = q(1)
         type is (real(kind=REAL64))
            this%target_lon = q(1)
         class default
            _FAIL('unsupport subclass for stretch params')
         end select
      end if


      hasLev=.false.
      hasLevel=.false.
      lev_name = 'lev'
      hasLev = file_metadata%has_dimension(lev_name)
      if (hasLev) then
         this%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
            this%lm = file_metadata%get_dimension(lev_name,rc=status)
            _VERIFY(status)
         end if
      end if

      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%im_world, this%jms, this%ny, min_DE_extent=2)
      call this%check_and_fill_consistency(rc=status)
      _VERIFY(status)

      _UNUSED_DUMMY(unusable)
      _UNUSED_DUMMY(force_file_coordinates)
   end subroutine initialize_from_file_metadata


   subroutine initialize_from_config_with_prefix(this, config, prefix, unusable, rc)
      use esmf
      class (CubedSphereGridFactory), 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=*), parameter :: Iam = MOD_NAME//'make_geos_grid_from_config'
      character(len=ESMF_MAXSTR) :: tmp
      type (ESMF_VM) :: vm
      integer :: vmcomm, ndes

      if (present(unusable)) print*,shape(unusable)

      call ESMF_VMGetCurrent(vm, rc=status)
      _VERIFY(status)

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

      call ESMF_ConfigGetAttribute(config, this%grid_type, label=prefix//'CS_GRID_TYPE:', default=FV_GRID_TYPE_DEFAULT)

      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%stretch_factor, label=prefix//'STRETCH_FACTOR:', default=MAPL_UNDEFINED_REAL)
      call ESMF_ConfigGetAttribute(config, this%target_lon, label=prefix//'TARGET_LON:', default=MAPL_UNDEFINED_REAL)
      call ESMF_ConfigGetAttribute(config, this%target_lat, label=prefix//'TARGET_LAT:', default=MAPL_UNDEFINED_REAL)

      call get_multi_integer(this%ims, 'IMS:', rc=status)
      _VERIFY(status)

      call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'JMS_FILE:', rc=status)
      if (status == _SUCCESS) then
         call get_jms_from_file(this%jms_2d, 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 this%check_and_fill_consistency(rc=status)
      _VERIFY(status)

      ! halo initialization

      call ESMF_VmGet(VM, mpicommunicator=vmcomm, petCount=ndes, 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, default=MAPL_UNDEFINED_INTEGER, 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,isPresent=isPresent,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_jms_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,k,face,total, unit, max_procs
         integer :: status, N_proc,NF
         integer, allocatable :: values_tmp(:), values_(:,:)


         N_proc = n*6 ! it has been devided by 6. get back the original NY
         allocate(values_tmp(N_proc), stat=status) ! no point in checking status
         _VERIFY(status)

         inquire(FILE = trim(file_name), EXIST=FileExists)
         if ( .not. FileExists) then
            print*, file_name //  " does not exist"
             _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, max_procs
            if (total /= N_proc) then
                print*, "n /= total"
                _RETURN(_FAILURE)
            endif
            do i = 1,total
                read(UNIT,*) values_tmp(i)
            enddo
            close(UNIT)
         endif

         call MAPL_CommsBcast(VM, max_procs,  n=1, ROOT=MAPL_Root, rc=status)
         call MAPL_CommsBcast(VM, values_tmp, n=N_proc, ROOT=MAPL_Root, rc=status)
         _VERIFY(STATUS)

         ! distributed to 6 faces
         allocate(values_(max_procs,6))
         values_ = 0
         k = 1
         do NF = 1, 6
            face = 0
            do i = 1, max_procs
               values_(i,NF) = values_tmp(k)
               face = face + values_tmp(k)
               k = k+1
               if (face == this%im_world) exit
            enddo
          enddo
          values = values_

         _RETURN(_SUCCESS)

      end subroutine get_jms_from_file

      subroutine get_bounds(bounds, label, rc)
         type(RealMinMax), intent(out) :: bounds
         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, bounds%min, rc=status)
         _VERIFY(status)
         call ESMF_ConfigGetAttribute(config, bounds%max, rc=status)
         _VERIFY(status)

         _RETURN(_SUCCESS)

      end subroutine get_bounds


   end subroutine initialize_from_config_with_prefix

   subroutine halo_init(this, halo_width,rc)
      class (CubedSphereGridFactory), intent(inout) :: this
      integer, optional, intent(in) :: halo_width
      integer, optional, intent(out) :: rc

      type(ESMF_Field) :: field
      type(ESMF_Grid), pointer :: grid
      integer :: useableHalo_width,status
      real, pointer :: ptr(:,:)
      character(len=*), parameter :: Iam = MOD_NAME // 'halo_init'

      if (present(halo_width)) then
         useableHalo_width=halo_width
      else
         useableHalo_width=1
      end if

      grid => this%get_grid(rc=status)
      _VERIFY(status)
      field = ESMF_FieldCreate(grid,ESMF_TYPEKIND_R4, &
              totalLWidth=[useableHalo_width,useableHalo_width], &
              totalUWidth=[useableHalo_width,useableHalo_width], &
              rc=status)
      _VERIFY(status)
      call ESMF_FieldGet(field,farrayPtr=ptr,rc=status)
      _VERIFY(status)
      ptr=0.0
      call ESMF_FieldHaloStore(field,this%rh,rc=status)
      _VERIFY(status)
      call ESMF_FieldDestroy(field,rc=status)
      _VERIFY(status)

   end subroutine halo_init

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

      _UNUSED_DUMMY(this)
      string = 'CubedSphereGridFactory'

   end function to_string

   subroutine check_and_fill_consistency(this, unusable, rc)
      use MAPL_BaseMod, only: MAPL_DecomposeDim
      use MAPL_Constants, only: PI => MAPL_PI_R8
      class (CubedSphereGridFactory), intent(inout) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      integer :: status
      character(len=*), parameter :: Iam = MOD_NAME // 'check_and_fill_consistency'

      _UNUSED_DUMMY(unusable)

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

      if (this%grid_type == MAPL_UNDEFINED_INTEGER) then
         this%grid_type = FV_GRID_TYPE_DEFAULT ! fv default
      end if

      if ( (this%target_lon /= MAPL_UNDEFINED_REAL) .and. &
           (this%target_lat /= MAPL_UNDEFINED_REAL) .and. &
           (this%stretch_factor /= MAPL_UNDEFINED_REAL) ) then
         _ASSERT(this%target_lat >= -90.0, 'Latitude should be greater than -90.0 degrees')
         _ASSERT(this%target_lat <= 90, 'Latitude should be less than 90.0 degrees')
         this%stretched_cube = .true.
         this%target_lon_degrees = this%target_lon
         this%target_lat_degrees = this%target_lat
         this%target_lon=this%target_lon*pi/180.d0
         this%target_lat=this%target_lat*pi/180.d0
      end if

      ! Check decomposition/bounds
      ! WY notes: not necessary for this assert
      !_ASSERT(allocated(this%ims) .eqv. allocated(this%jms),'inconsistent options')
      call verify(this%nx, this%im_world, this%ims, rc=status)
      if (allocated(this%jms_2d)) then
        _ASSERT(size(this%jms_2d,2)==6, 'incompatible shape')
        _ASSERT(sum(this%jms_2d) == 6*this%im_world, 'incompatible shape')
      else
         call verify(this%ny, this%im_world, this%jms, rc=status)
      endif

      _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, 'must be > 0 PEs in each dimension')

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

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

         else

            _ASSERT(n /= MAPL_UNDEFINED_INTEGER,'n not specified')
            _ASSERT(m_world /= MAPL_UNDEFINED_INTEGER,'m_world not specified')
            allocate(ms(n), stat=status)
            _VERIFY(status)

            call MAPL_DecomposeDim (m_world, ms, n, symmetric=.true.)

         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_real64(to, from, default)
      real(REAL64), intent(out) :: to
      real(REAL64), optional, intent(in) :: from
      real(REAL64), intent(in) :: default

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

   end subroutine set_with_default_real64

   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_bounds(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_bounds

   function decomps_are_equal(this, a) result(equal)
      class (CubedSphereGridFactory), intent(in) :: this
      class (AbstractGridFactory), intent(in) :: a
      integer :: a_nx,b_nx,a_ny,b_ny
      logical :: equal

      select type(a)
      class default
         equal = .false.
      class is (CubedSphereGridFactory)
         equal = .true.
         equal = size(a%ims) == size(this%ims)
         if (.not. equal) return
         equal = size(a%jms) == size(this%jms)
         if (.not. equal) return
         equal = all(a%ims == this%ims)
         if (.not. equal) return

         if ( allocated(a%jms) .and. allocated(this%jms)) then
            a_ny=size(a%jms)
            b_ny=size(this%ims)
            a_nx=size(a%ims)
            b_nx=size(this%ims)
            equal = a_nx*a_ny == b_nx*b_ny
            if (.not. equal) return
         else
            equal = all(a%jms_2d == this%jms_2d)
            if (.not. equal) return
         endif
      end select

   end function decomps_are_equal

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

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

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

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

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

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

      end select

   end function physical_params_are_equal

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

      select type (b)
      class default
         equals = .false.
         return
      class is (CubedSphereGridFactory)
         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

   subroutine initialize_from_esmf_distGrid(this, dist_grid, lon_array, lat_array, unusable, rc)
      class (CubedSphereGridFactory), 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

      character(len=*), parameter :: Iam = MOD_NAME // 'CubedSphereGridFactory_initialize_from_esmf_distGrid'

      _UNUSED_DUMMY(this)
      _UNUSED_DUMMY(dist_grid)
      _UNUSED_DUMMY(lon_array)
      _UNUSED_DUMMY(lat_array)
      _UNUSED_DUMMY(unusable)

      _FAIL('not implemented')

   end subroutine initialize_from_esmf_distGrid

   subroutine halo(this, array, unusable, halo_width, rc)
      use, intrinsic :: iso_fortran_env, only: REAL32
      class (CubedSphereGridFactory), 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
      character(len=*), parameter :: Iam = MOD_NAME // 'halo'
      type(ESMF_Field) :: field
      type(ESMF_Grid), pointer :: grid
      real, pointer :: ptr(:,:)
      integer :: useableHalo_width

      _UNUSED_DUMMY(unusable)

      if (.not. this%halo_initialized) then
         call this%halo_init(halo_width = halo_width)
         this%halo_initialized = .true.
      end if

      if (present(halo_width)) then
         useableHalo_width=halo_width
      else
         useableHalo_width=1
      end if
      grid => this%get_grid()
      field = ESMF_FieldCreate(grid,ESMF_TYPEKIND_R4, &
              totalLWidth=[useableHalo_width,useableHalo_width], &
              totalUWidth=[useableHalo_width,useableHalo_width], &
              rc=status)
      _VERIFY(status)
      call ESMF_FieldGet(field,farrayPtr=ptr,rc=status)
      _VERIFY(status)
      ptr = array
      call ESMF_FieldHalo(field,this%rh,rc=status)
      _VERIFY(status)
      array = ptr
      call ESMF_FieldDestroy(field,rc=status)
      _VERIFY(status)

      _RETURN(_SUCCESS)

   end subroutine halo

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

      character(len=4) :: im_string

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

      name = 'CF' // im_string //'x6C'

   end function generate_grid_name

   subroutine append_metadata(this, metadata)!, unusable, rc)
      class (CubedSphereGridFactory), intent(inout) :: this
      type (FileMetadata), intent(inout) :: metadata
!!$      class (KeywordEnforcer), optional, intent(in) :: unusable
!!$      integer, optional, intent(out) :: rc

      integer :: im
      type (Variable) :: v
      integer, parameter :: MAXLEN=80
      character(len=MAXLEN) :: gridspec_file_name
      !!! character(len=5), allocatable :: cvar(:,:)
      integer, allocatable :: ivar(:,:)
      integer, allocatable :: ivar2(:,:,:)

      real(REAL64), allocatable :: temp_coords(:)

      integer :: status
      integer, parameter :: ncontact = 4
      integer, parameter :: nf = 6

      ! Grid dimensions
      call metadata%add_dimension('Xdim', this%im_world, rc=status)
      call metadata%add_dimension('Ydim', this%im_world, rc=status)
      call metadata%add_dimension('XCdim', this%im_world+1, rc=status)
      call metadata%add_dimension('YCdim', this%im_world+1, rc=status)
      call metadata%add_dimension('nf', nf, rc=status)
      call metadata%add_dimension('ncontact', ncontact, rc=status)
      call metadata%add_dimension('orientationStrLen', 5, rc=status)

      ! Coordinate variables
      v = Variable(type=PFIO_REAL64, dimensions='Xdim')
      call v%add_attribute('long_name', 'Fake Longitude for GrADS Compatibility')
      call v%add_attribute('units', 'degrees_east')
      temp_coords = this%get_fake_longitudes()
      call metadata%add_variable('Xdim', CoordinateVariable(v, temp_coords))
      deallocate(temp_coords)

      v = Variable(type=PFIO_REAL64, dimensions='Ydim')
      call v%add_attribute('long_name', 'Fake Latitude for GrADS Compatibility')
      call v%add_attribute('units', 'degrees_north')
      temp_coords = this%get_fake_latitudes()
      call metadata%add_variable('Ydim', CoordinateVariable(v, temp_coords))
      deallocate(temp_coords)

      v = Variable(type=PFIO_INT32, dimensions='nf')
      call v%add_attribute('long_name','cubed-sphere face')
      call v%add_attribute('axis','e')
      call v%add_attribute('grads_dim','e')
      call v%add_const_value(UnlimitedEntity([1,2,3,4,5,6]))
      call metadata%add_variable('nf',v)

      v = Variable(type=PFIO_INT32, dimensions='ncontact')
      call v%add_attribute('long_name','number of contact points')
      call v%add_const_value(UnlimitedEntity([1,2,3,4]))
      call metadata%add_variable('ncontact',v)

      ! Other variables
      allocate(ivar(4,6))
      ivar = reshape( [5, 3, 2, 6, &
                       1, 3, 4, 6, &
                       1, 5, 4, 2, &
                       3, 5, 6, 2, &
                       3, 1, 6, 4, &
                       5, 1, 2, 4 ], [ncontact,nf])
      v = Variable(type=PFIO_INT32, dimensions='ncontact,nf')
      call v%add_attribute('long_name', 'adjacent face starting from left side going clockwise')
      call v%add_const_value(UnlimitedEntity(ivar))
      call metadata%add_variable('contacts', v)

      !!! At present pfio does not seem to work with string variables
      !!! allocate(cvar(4,6))
      !!! cvar =reshape([" Y:-X", " X:-Y", " Y:Y ", " X:X ", &
      !!!                " Y:Y ", " X:X ", " Y:-X", " X:-Y", &
      !!!                " Y:-X", " X:-Y", " Y:Y ", " X:X ", &
      !!!                " Y:Y ", " X:X ", " Y:-X", " X:-Y", &
      !!!                " Y:-X", " X:-Y", " Y:Y ", " X:X ", &
      !!!                " Y:Y ", " X:X ", " Y:-X", " X:-Y" ], [ncontact,nf])
      !!! v = Variable(type=PFIO_STRING, dimensions='orientationStrLen,ncontact,nf')
      !!! call v%add_attribute('long_name', 'orientation of boundary')
      !!! call v%add_const_value(UnlimitedEntity(cvar))
      !!! call metadata%add_variable('orientation', v)

      im = this%im_world
      allocate(ivar2(4,4,6))
      ivar2 = reshape(            &
               [[im, im,  1, im,  &
                  1, im,  1,  1,  &
                  1, im,  1,  1,  &
                 im, im,  1, im], &
                [im,  1, im, im,  &
                  1,  1, im,  1,  &
                  1,  1, im,  1,  &
                 im,  1, im, im], &
                [im, im,  1, im,  &
                  1, im,  1,  1,  &
                  1, im,  1,  1,  &
                 im, im,  1, im], &
                [im,  1, im, im,  &
                  1,  1, im,  1,  &
                  1,  1, im,  1,  &
                 im,  1, im, im], &
                [im, im,  1, im,  &
                  1, im,  1,  1,  &
                  1, im,  1,  1,  &
                 im, im,  1, im], &
                [im,  1, im, im,  &
                  1,  1, im,  1,  &
                  1,  1, im,  1,  &
                 im,  1, im, im] ], [ncontact,ncontact,nf])
      v = Variable(type=PFIO_INT32, dimensions='ncontact,ncontact,nf')
      call v%add_attribute('long_name', 'anchor point')
      call v%add_const_value(UnlimitedEntity(ivar2))
      call metadata%add_variable('anchor', v)

      call Metadata%add_attribute('grid_mapping_name', 'gnomonic cubed-sphere')
      call Metadata%add_attribute('file_format_version', '2.92')
      call Metadata%add_attribute('additional_vars', 'contacts,orientation,anchor')
      write(gridspec_file_name,'("C",i0,"_gridspec.nc4")') this%im_world
      call Metadata%add_attribute('gridspec_file', trim(gridspec_file_name))

      v = Variable(type=PFIO_REAL64, dimensions='Xdim,Ydim,nf')
      call v%add_attribute('long_name','longitude')
      call v%add_attribute('units','degrees_east')
      call metadata%add_variable('lons',v)

      v = Variable(type=PFIO_REAL64, dimensions='Xdim,Ydim,nf')
      call v%add_attribute('long_name','latitude')
      call v%add_attribute('units','degrees_north')
      call metadata%add_variable('lats',v)

      v = Variable(type=PFIO_REAL64, dimensions='XCdim,YCdim,nf')
      call v%add_attribute('long_name','longitude')
      call v%add_attribute('units','degrees_east')
      call metadata%add_variable('corner_lons',v)

      v = Variable(type=PFIO_REAL64, dimensions='XCdim,YCdim,nf')
      call v%add_attribute('long_name','latitude')
      call v%add_attribute('units','degrees_north')
      call metadata%add_variable('corner_lats',v)

      if (this%stretched_cube) then
         call metadata%add_attribute('STRETCH_FACTOR',this%stretch_factor)
         call metadata%add_attribute('TARGET_LON',this%target_lon_degrees)
         call metadata%add_attribute('TARGET_LAT',this%target_lat_degrees)
      end if

   end subroutine append_metadata

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

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

      vars = 'Xdim,Ydim,nf'

   end function get_grid_vars

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

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

      vars = 'Xdim,Ydim,nf,anchor,lons,lats,corner_lons,corner_lats,nf,ncontact,cubed_sphere,contacts,orientation'

   end function get_file_format_vars

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

      call var%add_attribute('coordinates','lons lats')
      call var%add_attribute('grid_mapping','cubed_sphere')

   end subroutine append_variable_metadata

   ! Routine to return fake coordinates for Grads+netCDF use cases.
   ! Result is the longitudes (in degrees) of cell centers along the
   ! "equator" of face 1.
   function get_fake_longitudes(this, unusable, rc) result(longitudes)
      use mpi
      use MAPL_BaseMod, only: MAPL_Grid_Interior
      real (kind=REAL64), allocatable :: longitudes(:)
      class (CubedSphereGridFactory), intent(inout) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      type (ESMF_Grid) :: grid
      real (kind=ESMF_KIND_R8), pointer :: centers(:,:)
      real (kind=REAL64), allocatable :: piece(:)
      integer :: n_loc
      integer, allocatable :: counts(:), displs(:)
      type (ESMF_VM) :: vm
      integer :: ierror
      integer :: npes, p, pet
      integer :: comm_grid
      integer :: i_1, i_n, j_1, j_n
      integer :: j_mid
      integer :: tile
      integer :: status

      character(len=*), parameter :: Iam = MOD_NAME // 'get_fake_longitudes()'

      _UNUSED_DUMMY(unusable)

      grid = this%make_grid()

      call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, &
           staggerloc=ESMF_STAGGERLOC_CENTER, &
           farrayPtr=centers, rc=status)
      _VERIFY(status)

      call ESMF_VMGetCurrent(vm, rc=status)
      _VERIFY(status)

      call ESMF_VMGet(vm, mpiCommunicator=comm_grid, petcount=npes, localpet=pet, rc=status)
      _VERIFY(status)

      call MAPL_grid_interior(grid, i_1, i_n, j_1, j_n)

      j_mid = 1 + this%im_world/2

      tile = 1 + (j_1-1)/this%im_world
      if (tile == 1 .and. (j_1 <= j_mid) .and. (j_mid <= j_n)) then
         allocate(piece(i_1:i_n))
         piece(:) = centers(:,j_mid-(j_1-1))
         n_loc = (i_n - i_1 + 1)
      else
         allocate(piece(1)) ! MPI does not like 0-sized arrays.
         piece(1) = 0
         n_loc = 0
      end if

      allocate(counts(0:npes-1), displs(0:npes-1))

      call MPI_Allgather(n_loc, 1, MPI_INTEGER, counts, 1, MPI_INTEGER, comm_grid, ierror)
      _VERIFY(ierror)

      displs(0) = 0
      do p = 1, npes-1
         displs(p) = displs(p-1) + counts(p-1)
      end do

      allocate(longitudes(this%im_world))
      call MPI_Allgatherv(piece, n_loc, MPI_REAL8, longitudes, counts, displs, MPI_REAL8, comm_grid, ierror)
      _VERIFY(ierror)

      longitudes = longitudes * MAPL_RADIANS_TO_DEGREES

   end function get_fake_longitudes

   function get_fake_latitudes(this, unusable, rc) result(latitudes)
      use mpi
      use MAPL_BaseMod, only: MAPL_Grid_Interior
      real (kind=REAL64), allocatable :: latitudes(:)
      class (CubedSphereGridFactory), intent(inout) :: this
      class (KeywordEnforcer), optional, intent(in) :: unusable
      integer, optional, intent(out) :: rc

      type (ESMF_Grid) :: grid
      real (kind=ESMF_KIND_R8), pointer :: centers(:,:)
      real (kind=REAL64), allocatable :: piece(:)
      integer :: n_loc
      integer, allocatable :: counts(:), displs(:)
      type (ESMF_VM) :: vm
      integer :: ierror
      integer :: npes, p, pet
      integer :: comm_grid
      integer :: i_1, i_n, j_1, j_n
      integer :: j_mid
      integer :: tile
      integer :: status

      character(len=*), parameter :: Iam = MOD_NAME // 'get_fake_latitudes()'

      _UNUSED_DUMMY(unusable)

      grid = this%make_grid()

      call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, &
           staggerloc=ESMF_STAGGERLOC_CENTER, &
           farrayPtr=centers, rc=status)
      _VERIFY(status)

      call ESMF_VMGetCurrent(vm, rc=status)
      _VERIFY(status)

      call ESMF_VMGet(vm, mpiCommunicator=comm_grid, petcount=npes, localpet=pet, rc=status)
      _VERIFY(status)

      call MAPL_grid_interior(grid, i_1, i_n, j_1, j_n)

      j_mid = 1 + this%im_world/2

      tile = 1 + (j_1-1)/this%im_world
      if (tile == 1 .and. (i_1 <= j_mid) .and. (j_mid <= i_n)) then
         allocate(piece(j_1:j_n))
         piece(:) = centers(j_mid-(i_1-1),:)
         n_loc = (j_n - j_1 + 1)
      else
         allocate(piece(1)) ! MPI does not like 0-sized arrays.
         piece(1) = 0
         n_loc = 0
      end if

      allocate(counts(0:npes-1), displs(0:npes-1))

      call MPI_Allgather(n_loc, 1, MPI_INTEGER, counts, 1, MPI_INTEGER, comm_grid, ierror)
      _VERIFY(ierror)

      displs(0) = 0
      do p = 1, npes-1
         displs(p) = displs(p-1) + counts(p-1)
      end do

      allocate(latitudes(this%im_world))
      call MPI_Allgatherv(piece, n_loc, MPI_REAL8, latitudes, counts, displs, MPI_REAL8, comm_grid, ierror)
      _VERIFY(ierror)

      latitudes = latitudes * MAPL_RADIANS_TO_DEGREES

   end function get_fake_latitudes

   subroutine generate_file_bounds(this,grid,local_start,global_start,global_count,metaData,rc)
      use MAPL_BaseMod
      class(CubedSphereGridFactory), 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,tile
      character(len=*), parameter :: Iam = MOD_NAME // 'generate_file_bounds'
      logical :: face_format


      if (present(metadata)) then
         face_format = metadata%has_dimension('nf')
      else
         face_format = .true.
      end if
      call MAPL_GridGet(grid,globalCellCountPerDim=global_dim,rc=status)
      _VERIFY(status)
      call MAPL_GridGetInterior(grid,i1,in,j1,jn)
      if (face_format) then
         tile = 1 + (j1-1)/global_dim(1)
         allocate(local_start,source=[i1,j1-(tile-1)*global_dim(1),tile])
         allocate(global_start,source=[1,1,1])
         allocate(global_count,source=[global_dim(1),global_dim(1),6])
      else
         allocate(local_start,source=[i1,j1])
         allocate(global_start,source=[1,1])
         allocate(global_count,source=[global_dim(1),global_dim(2)])
      end if

      _RETURN(_SUCCESS)
      _UNUSED_DUMMY(this)

   end subroutine generate_file_bounds

   subroutine generate_file_corner_bounds(this,grid,local_start,global_start,global_count,rc)
      use MAPL_BaseMod
      class(CubedSphereGridFactory), 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

      integer :: status
      integer :: global_dim(3),i1,j1,in,jn,tile
      integer :: face_j1, is, js
      character(len=*), parameter :: Iam = MOD_NAME // 'generate_file_bounds'
      _UNUSED_DUMMY(this)

      call MAPL_GridGet(grid,globalCellCountPerDim=global_dim,rc=status)
      _VERIFY(status)
      call MAPL_GridGetInterior(grid,i1,in,j1,jn)
      tile = 1 + (j1-1)/global_dim(1)
      face_j1 = j1-(tile-1)*global_dim(1)
      is = i1
      js = face_j1
      allocate(local_start,source=[is,js,tile])
      allocate(global_start,source=[1,1,1])
      allocate(global_count,source=[global_dim(1)+1,global_dim(1)+1,6])

      _RETURN(_SUCCESS)

   end subroutine generate_file_corner_bounds

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

   function generate_file_reference3D(this,fpointer,metadata) result(ref)
      use pFIO
      use, intrinsic :: ISO_C_BINDING
      type(ArrayReference) :: ref
      class(CubedSphereGridFactory), intent(inout) :: this
      real, pointer, intent(in) :: fpointer(:,:,:)
      type(FileMetadata), intent(in), optional :: metaData
      type(c_ptr) :: cptr
      real, pointer :: ptr_ref(:,:,:,:,:)
      logical :: face_format

      if (present(metadata)) then
         face_format = metadata%has_dimension('nf')
      else
         face_format = .true.
      end if

      if (face_format) then
         cptr = c_loc(fpointer)
         call C_F_pointer(cptr,ptr_ref,[size(fpointer,1),size(fpointer,2),1,size(fpointer,3),1])
         ref = ArrayReference(ptr_ref)
      else
         ref = ArrayReference(fpointer)
      end if

      _UNUSED_DUMMY(this)
   end function generate_file_reference3D

end module MAPL_CubedSphereGridFactoryMod