ESMFL_Mod.F90 Source File


This file depends on

sourcefile~~esmfl_mod.f90~~EfferentGraph sourcefile~esmfl_mod.f90 ESMFL_Mod.F90 sourcefile~base_base.f90 Base_Base.F90 sourcefile~esmfl_mod.f90->sourcefile~base_base.f90 sourcefile~constants.f90 Constants.F90 sourcefile~esmfl_mod.f90->sourcefile~constants.f90 sourcefile~mapl_abstractgridfactory.f90 MAPL_AbstractGridFactory.F90 sourcefile~esmfl_mod.f90->sourcefile~mapl_abstractgridfactory.f90 sourcefile~mapl_comms.f90 MAPL_Comms.F90 sourcefile~esmfl_mod.f90->sourcefile~mapl_comms.f90 sourcefile~mapl_exceptionhandling.f90 MAPL_ExceptionHandling.F90 sourcefile~esmfl_mod.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~mapl_gridmanager.f90 MAPL_GridManager.F90 sourcefile~esmfl_mod.f90->sourcefile~mapl_gridmanager.f90 sourcefile~base_base.f90->sourcefile~constants.f90 sourcefile~mapl_keywordenforcer.f90 MAPL_KeywordEnforcer.F90 sourcefile~base_base.f90->sourcefile~mapl_keywordenforcer.f90 sourcefile~mapl_range.f90 MAPL_Range.F90 sourcefile~base_base.f90->sourcefile~mapl_range.f90 sourcefile~maplgrid.f90 MaplGrid.F90 sourcefile~base_base.f90->sourcefile~maplgrid.f90 sourcefile~internalconstants.f90 InternalConstants.F90 sourcefile~constants.f90->sourcefile~internalconstants.f90 sourcefile~mathconstants.f90 MathConstants.F90 sourcefile~constants.f90->sourcefile~mathconstants.f90 sourcefile~physicalconstants.f90 PhysicalConstants.F90 sourcefile~constants.f90->sourcefile~physicalconstants.f90 sourcefile~mapl_abstractgridfactory.f90->sourcefile~base_base.f90 sourcefile~mapl_abstractgridfactory.f90->sourcefile~constants.f90 sourcefile~mapl_abstractgridfactory.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~mapl_abstractgridfactory.f90->sourcefile~mapl_keywordenforcer.f90 sourcefile~pfio.f90 pFIO.F90 sourcefile~mapl_abstractgridfactory.f90->sourcefile~pfio.f90 sourcefile~mapl_comms.f90->sourcefile~base_base.f90 sourcefile~mapl_comms.f90->sourcefile~constants.f90 sourcefile~mapl_comms.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~shmem.f90 Shmem.F90 sourcefile~mapl_comms.f90->sourcefile~shmem.f90 sourcefile~mapl_errorhandling.f90 MAPL_ErrorHandling.F90 sourcefile~mapl_exceptionhandling.f90->sourcefile~mapl_errorhandling.f90 sourcefile~mapl_throw.f90 MAPL_Throw.F90 sourcefile~mapl_exceptionhandling.f90->sourcefile~mapl_throw.f90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_abstractgridfactory.f90 sourcefile~mapl_cubedspheregridfactory.f90 MAPL_CubedSphereGridFactory.F90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_cubedspheregridfactory.f90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_errorhandling.f90 sourcefile~mapl_externalgridfactory.f90 MAPL_ExternalGridFactory.F90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_externalgridfactory.f90 sourcefile~mapl_integer64gridfactorymap.f90 MAPL_Integer64GridFactoryMap.F90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_integer64gridfactorymap.f90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_keywordenforcer.f90 sourcefile~mapl_latlongridfactory.f90 MAPL_LatLonGridFactory.F90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_latlongridfactory.f90 sourcefile~mapl_stringgridfactorymap.f90 MAPL_StringGridFactoryMap.F90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_stringgridfactorymap.f90 sourcefile~mapl_swathgridfactory.f90 MAPL_SwathGridFactory.F90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_swathgridfactory.f90 sourcefile~mapl_tripolargridfactory.f90 MAPL_TripolarGridFactory.F90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_tripolargridfactory.f90 sourcefile~mapl_xygridfactory.f90 MAPL_XYGridFactory.F90 sourcefile~mapl_gridmanager.f90->sourcefile~mapl_xygridfactory.f90 sourcefile~mapl_gridmanager.f90->sourcefile~pfio.f90

Files dependent on this one

sourcefile~~esmfl_mod.f90~~AfferentGraph sourcefile~esmfl_mod.f90 ESMFL_Mod.F90 sourcefile~base.f90 Base.F90 sourcefile~base.f90->sourcefile~esmfl_mod.f90 sourcefile~base_base_implementation.f90 Base_Base_implementation.F90 sourcefile~base_base_implementation.f90->sourcefile~esmfl_mod.f90 sourcefile~comp_testing_driver.f90 Comp_Testing_Driver.F90 sourcefile~comp_testing_driver.f90->sourcefile~esmfl_mod.f90 sourcefile~extdatagridcompmod.f90 ExtDataGridCompMod.F90 sourcefile~extdatagridcompmod.f90->sourcefile~esmfl_mod.f90 sourcefile~extdatagridcompng.f90 ExtDataGridCompNG.F90 sourcefile~extdatagridcompng.f90->sourcefile~esmfl_mod.f90 sourcefile~extdatamasking.f90 ExtDataMasking.F90 sourcefile~extdatamasking.f90->sourcefile~esmfl_mod.f90 sourcefile~genericcplcomp.f90 GenericCplComp.F90 sourcefile~genericcplcomp.f90->sourcefile~esmfl_mod.f90 sourcefile~griddedio.f90 GriddedIO.F90 sourcefile~griddedio.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_bundleio_test.f90 mapl_bundleio_test.F90 sourcefile~mapl_bundleio_test.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_capgridcomp.f90 MAPL_CapGridComp.F90 sourcefile~mapl_capgridcomp.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_cfio.f90 MAPL_CFIO.F90 sourcefile~mapl_cfio.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_epochswathmod.f90 MAPL_EpochSwathMod.F90 sourcefile~mapl_epochswathmod.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_generic.f90 MAPL_Generic.F90 sourcefile~mapl_generic.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_historygridcomp.f90 MAPL_HistoryGridComp.F90 sourcefile~mapl_historygridcomp.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_locstreammod.f90 MAPL_LocStreamMod.F90 sourcefile~mapl_locstreammod.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_orbgridcompmod.f90 MAPL_OrbGridCompMod.F90 sourcefile~mapl_orbgridcompmod.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_resource.f90 MAPL_Resource.F90 sourcefile~mapl_resource.f90->sourcefile~esmfl_mod.f90 sourcefile~mapl_simplebundlemod.f90 MAPL_SimpleBundleMod.F90 sourcefile~mapl_simplebundlemod.f90->sourcefile~esmfl_mod.f90 sourcefile~regrid_util.f90 Regrid_Util.F90 sourcefile~regrid_util.f90->sourcefile~esmfl_mod.f90

Source Code

!------------------------------------------------------------------------------
!               Global Modeling and Assimilation Office (GMAO)                !
!                    Goddard Earth Observing System (GEOS)                    !
!                                 MAPL Component                              !
!------------------------------------------------------------------------------
!
#include "MAPL_ErrLog.h"
!
!>
!### MODULE: `ESMFL_MOD`
!
! Author: GMAO SI-Team
!

#if 0
stubbed routines
 ESMFL_RegridStore
 FieldRegrid1
 BundleRegrid1
 BundlePrep_ (no grid attributes)
 ESMFL_FieldGetDims
 ESMFL_GridDistBlockSet
#endif

#include "unused_dummy.H"

module ESMFL_MOD

  !BOP
  ! !MODULE: ESMFL_MOD

  !USES:
  use ESMF
  use MAPL_Constants
  use MAPL_BaseMod
  use MAPL_CommsMod
  use MAPL_ExceptionHandling
  use, intrinsic :: iso_fortran_env, only: REAL32, REAL64
  implicit none
  private

!
!ALT These need to be changed
!    values here are just to compile
!

  !PUBLIC MEMBER FUNCTIONS:
  public ESMFL_StateGetField
  public ESMFL_StateGetFieldArray
  public ESMFL_StateGetPointerToData
  public ESMFL_FieldGetPointerToData
  public ESMFL_BundleGetPointerToData
  public ESMFL_BundleCpyField
  public ESMFL_GridCoordGet
  ! public ESMFL_Connect2STATE
  public ESMFL_FCOLLECT
  public ESMFL_StateFreePointers
  public ESMFL_StateSetFieldNeeded
  public ESMFL_StateFieldIsNeeded
  public ESMFL_FieldGetDims
  public ESMFL_GridDistBlockSet

  public ESMFL_FieldRegrid ! alt: this should be MAPL_FieldRegrid
                           !      (topo_bin may need to be here)

  public ESMFL_RegridStore ! only used for regridding using ESMF_FieldRegrid
  public ESMFL_Regrid
  public ESMFL_State2Bundle
  public ESMFL_Bundle2State
  public ESMFL_Bundles2Bundle
  public ESMFL_Add2Bundle
  public ESMFL_HALO
  public ESMFL_BundleAddState
  public MAPL_AreaMean
  public ESMFL_Diff
  public ESMFL_statistics
  public ESMFL_field_is_undefined
  !EOP

  ! regridding
  interface ESMFL_Regrid

     module procedure BundleRegrid  ! Uses Larry's hinterp
     module procedure StateRegrid   ! Uses Larry's hinterp

     module procedure FieldRegrid1  ! Uses ESMF regrid
     module procedure BundleRegrid1 ! Uses ESMF regrid

  end interface

  ! compare two states or bundles
  interface ESMFL_Diff
    module procedure StateDiff
    module procedure BundleDiff
  end interface

  interface ESMFL_statistics
    module procedure StateStatistics
    module procedure BundleStatistics
  end interface

  ! Extract fields from a State and place them in a Bundle
  interface ESMFL_State2Bundle
    module procedure State2Bundle
  end interface
  ! Extract fields from a Bundle and place them in a State
  interface ESMFL_Bundle2State
    module procedure Bundle2State
  end interface
  ! Extract fields from Bundles and place them in a single Bundle
  interface ESMFL_Bundles2Bundle
    module procedure Bundles2Bundle
  end interface
  interface ESMFL_Add2Bundle
    module procedure Add2Bundle
  end interface

  interface ESMFL_BundleAddState
    module procedure BundleAddState_
  end interface

  interface ESMFL_StateGetPointerToData
     module procedure ESMFL_StateGetPtrToDataR4_1
     module procedure ESMFL_StateGetPtrToDataR4_2
     module procedure ESMFL_StateGetPtrToDataR4_3
     module procedure ESMFL_StateGetPtrToDataR4_4
     module procedure ESMFL_StateGetPtrToDataR8_1
     module procedure ESMFL_StateGetPtrToDataR8_2
     module procedure ESMFL_StateGetPtrToDataR8_3
     module procedure ESMFL_StateGetPtrToDataR8_4
  end interface

  interface ESMFL_FieldGetPointerToData
     module procedure ESMFL_FieldGetPtrToDataR4_1
     module procedure ESMFL_FieldGetPtrToDataR4_2
     module procedure ESMFL_FieldGetPtrToDataR4_3
     module procedure ESMFL_FieldGetPtrToDataR4_4
     module procedure ESMFL_FieldGetPtrToDataR8_1
     module procedure ESMFL_FieldGetPtrToDataR8_2
     module procedure ESMFL_FieldGetPtrToDataR8_3
     module procedure ESMFL_FieldGetPtrToDataR8_4
  end interface

  interface ESMFL_BundleGetPointerToData
     module procedure ESMFL_BundleGetPointerByIndex2
     module procedure ESMFL_BundleGetPointerByIndex3
     module procedure ESMFL_BundleGetPointerByName2
     module procedure ESMFL_BundleGetPointerByName3
  end interface

  interface ESMFL_FCOLLECT
     module procedure ESMFL_FCOLLECT_I4
     module procedure ESMFL_FCOLLECT_R4
     module procedure ESMFL_FCOLLECT_R8
  end interface

  interface ESMFL_HALO
     module procedure ESMFL_HALO_R4_2D
  end interface

  interface AdjustPtrBounds
     module procedure AdjustPtrBounds1dr4
     module procedure AdjustPtrBounds1dr8
     module procedure AdjustPtrBounds3dr4
     module procedure AdjustPtrBounds3dr8
  end interface

  interface MAPL_AreaMean
     module procedure MAPL_AreaMean_2d_r8_bitrep
     module procedure MAPL_AreaMean_2d_r8
  end interface

  Logical            :: verbose

!EOP
!
! !REVISION HISTORY:
!
!  17Aug2007  Todling  Implemented verbose
!
!-------------------------------------------------------------------------

contains

 subroutine ESMFL_StateGetFieldArray(state, NAME, ARRAY, RC)
   type(ESMF_State),  intent(IN)  :: state
   type(ESMF_Array),  intent(OUT) :: array
   character(len=*),   intent(IN) :: name
   integer, optional, intent(OUT) :: rc


   integer          :: status
   type(ESMF_Field) :: field
   type(ESMF_StateItem_Flag) :: itemType

   call ESMF_StateGet(state, name, itemType=itemType, rc=status)
   _VERIFY(STATUS)

   if (itemType /= ESMF_STATEITEM_FIELD) then
      if (present(RC)) then
         RC=STATUS
         return
      end if
   else
      call ESMF_StateGet(state, name, field, rc=status)
      _VERIFY(STATUS)

      call ESMF_FieldGet(field, Array=array, rc=status)
      _VERIFY(STATUS)
   end if

   _RETURN(ESMF_SUCCESS)
 end subroutine ESMFL_StateGetFieldArray


 subroutine ESMFL_StateGetField(State, FieldName, Bundle, FieldAlias, RC)
   type(ESMF_State),           intent(IN   ) :: State
   character(len=*),           intent(IN   ) :: FieldName(:)
   type(ESMF_FieldBundle),     intent(INOUT) :: Bundle
   character(len=*), optional, intent(IN   ) :: FieldAlias(:)
   integer, optional,  intent(OUT)   :: rc

   integer                    :: status

   character(len=ESMF_MAXSTR) :: NameInBundle
   character(len=ESMF_MAXSTR) :: NameInState
   type(ESMF_Field)           :: Bundlefield
   type(ESMF_Field)           :: Statefield
   integer                    :: I
   logical                    :: NotInState
   logical                    :: NeedNewField

   type(ESMF_StateItem_Flag)  :: itemType
   logical                    :: isPresent

!   integer                    :: gridRank

!  Adds fields from a state to a bundle. If a field
!  is not in the state, it adds a dummy field. It fails
!  if it tries to add a dummy field to an empty bundle,
!  since it does not know what grid to use in the dummy.
!  The name of the field in the bundle can be reset with
!  FieldAlias

   RequestedFields: do i=1,size(FIELDNAME)

      NameInState = FieldName (I)

      if(present(FieldAlias)) then
         NameInBundle = FieldAlias(I)
      else
         NameInBundle = NameInState
      end if

! Make sure field is not in bundle

      call ESMF_FieldBundleGet(BUNDLE, FIELDNAME=NameInBundle, isPresent=isPresent, rc=STATUS)
      _VERIFY(STATUS)
      _ASSERT(.not. isPresent, trim(NameInBundle) // ' found in bundle.')

! Get Field from State

      call ESMF_StateGet (STATE, NameInState, itemType=itemType, RC=STATUS)
      _VERIFY(STATUS)

      if (itemType /= ESMF_STATEITEM_FIELD) then
         NotInState = .TRUE.
      else
         NotInState = .FALSE.
         call ESMF_StateGet (STATE, NameInState, StateFIELD,  RC=STATUS)
         _VERIFY(STATUS)
      end if

      NeedNewField = NotInState .or. (NameInState/=NameInBundle)

      if (.not. NotInState) then
         call MAPL_AllocateCoupling(stateField, rc=status)
         _VERIFY(STATUS)
      end if

      if(NeedNewField) then

         ! Define Grid and Array for new field

         if(NotInState) then

            ! Create a new empty field
            BundleField = ESMF_FieldEmptyCreate(name=NameInBundle, rc=status)
            _VERIFY(STATUS)

         else

            ! Use the grid and array in the State Field, just rename it
            BundleField = MAPL_FieldCreate(stateField, name=NameInBundle, RC=STATUS)
            _VERIFY(STATUS)
         end if

      else

         BundleField = StateField

      end if

      call MAPL_FieldBundleAdd(BUNDLE, bundleField, rc=STATUS)
      _VERIFY(STATUS)

   end do RequestedFields

   _RETURN(ESMF_SUCCESS)
 end subroutine ESMFL_StateGetField



!BOP

! !IROUTINE: ESMFL_GridCoordGet - retrieves the coordinates of a particular axis in radians

! !INTERFACE:
subroutine ESMFL_GridCoordGet(GRID, coord, name, Location, Units, rc)

  !ARGUMENTS:
   type(ESMF_Grid),   intent(INout ) :: GRID
   real, dimension(:,:), pointer  :: coord
   character (len=*) , intent(IN) :: name
   type(ESMF_StaggerLoc)          :: location
   integer                        :: units
   integer, optional              :: rc
 !EOP

! local variables
  integer                   :: rank
  type(ESMF_CoordSys_Flag)  :: crdSys
  type(ESMF_TypeKind_Flag)  :: tk
  integer                   :: counts(ESMF_MAXDIM)
  integer                   :: crdOrder
  real(ESMF_KIND_R4), pointer :: r4d2(:,:)
  real(ESMF_KIND_R8), pointer :: r8d2(:,:)
  real(ESMF_KIND_R8)        :: conv2rad
  integer                   :: STATUS
  integer                   :: i
  integer                   :: j
!ALT  integer                   :: i1, in, j1, jn
  integer                   :: coordDimCount(ESMF_MAXDIM)
  character(len=ESMF_MAXSTR):: gridname

  _UNUSED_DUMMY(Units)

  call ESMF_GridGet(grid, coordSys=crdSys, coordTypeKind=tk, &
          dimCount=rank, coordDimCount=coordDimCount, rc=status)
  _VERIFY(STATUS)

  if (name == "Longitude") then
    crdOrder = 1
  else if (name == "Latitude") then
    crdOrder = 2
  else
   STATUS=ESMF_FAILURE
   _VERIFY(STATUS)
  endif

  call ESMF_GridGet(grid, name=gridname, rc=status)
  _VERIFY(STATUS)

  if (gridname(1:10) == 'tile_grid_') then

     call MAPL_GridGet(GRID, localCellCountPerDim=counts, rc=status)
     _VERIFY(STATUS)
     allocate(coord(counts(1), counts(2)), STAT=status)
     _VERIFY(STATUS)
     coord = 0.0 ! initialize just in case

!ALT     call MAPL_GRID_INTERIOR(GRID,I1,IN,J1,JN)

     do I=1, counts(1)
        do J=1, counts(2)
           if (crdOrder == 1) then
              coord(i,j) = i !ALT + I1 - 1
           else
              coord(i,j) = j !ALT + J1 - 1
           end if
        end do
     end do

     coord = coord * (MAPL_PI_R8 / 180.d+0)

     _RETURN(ESMF_SUCCESS)
  end if

  if (crdSys == ESMF_COORDSYS_SPH_DEG) then
     conv2rad = MAPL_PI_R8 / 180._ESMF_KIND_R8
  else if (crdSys == ESMF_COORDSYS_SPH_RAD) then
     conv2rad = 1._ESMF_KIND_R8
  else
     _FAIL('Unsupported coordinate system:  ESMF_COORDSYS_CART')
  end if

  if (tk == ESMF_TYPEKIND_R4) then
     if (coordDimCount(crdOrder)==2) then
        call ESMF_GridGetCoord(grid, localDE=0, coordDim=crdOrder, &
             staggerloc=location, &
             computationalCount=COUNTS,  &
             farrayPtr=R4D2, rc=status)
        _VERIFY(STATUS)
        allocate(coord(counts(1), counts(2)), STAT=status)
        _VERIFY(STATUS)
        coord = conv2rad * R4D2
     else
        _RETURN(ESMF_FAILURE)
     endif
  else if (tk == ESMF_TYPEKIND_R8) then
     if (coordDimCount(crdOrder)==2) then
        call ESMF_GridGetCoord(grid, localDE=0, coordDim=crdOrder, &
             staggerloc=location, &
             computationalCount=COUNTS,  &
             farrayPtr=R8D2, rc=status)
        _VERIFY(STATUS)
        allocate(coord(counts(1), counts(2)), STAT=status)
        _VERIFY(STATUS)
        coord = conv2rad * R8D2
     else
        _RETURN(ESMF_FAILURE)
     endif
  else
     _RETURN(ESMF_FAILURE)
  endif
  _RETURN(ESMF_SUCCESS)

 end subroutine ESMFL_GridCoordGet

 subroutine ESMFL_StateFreePointers(STATE, RC)
   type(ESMF_State),  intent(INOUT) :: STATE
   integer, optional, intent(  OUT) :: RC

   integer                          :: STATUS

   type(ESMF_Array)                 :: ARRAY
   type(ESMF_Field)                 :: FIELD
   integer                          :: RANK
   integer                          :: I
   integer                          :: ITEMCOUNT
   real, pointer                    :: PTR1(:)
   real, pointer                    :: PTR2(:,:)
   real, pointer                    :: PTR3(:,:,:)
   real, pointer                    :: PTR4(:,:,:,:)
   logical                          :: NEEDED

   character (len=ESMF_MAXSTR), pointer :: ITEMNAMELIST(:)
   type(ESMF_StateItem_Flag)  , pointer :: ITEMTYPELIST(:)

   type (ESMF_LocalArray), target  :: larrayList(1)
   type (ESMF_LocalArray), pointer :: larray
   integer        :: localDeCount

   logical :: isPresent

! Get information from state
!---------------------------

   call ESMF_StateGet(STATE,ITEMCOUNT=ITEMCOUNT,RC=STATUS)
   _VERIFY(STATUS)

   if(ITEMCOUNT==0) then
      _RETURN(ESMF_SUCCESS)
   end if

   allocate(ITEMNAMELIST(ITEMCOUNT),STAT=STATUS)
   _VERIFY(STATUS)
   allocate(ITEMTYPELIST(ITEMCOUNT),STAT=STATUS)
   _VERIFY(STATUS)

   call ESMF_StateGet(STATE,ITEMNAMELIST=ITEMNAMELIST,ITEMTYPELIST=ITEMTYPELIST,RC=STATUS)
   _VERIFY(STATUS)

   do I=1,ITEMCOUNT
      if(ITEMTYPELIST(I)==ESMF_STATEITEM_FIELD) then
         call ESMF_StateGet(STATE, trim(ITEMNAMELIST(I)), FIELD, RC=STATUS)
         _VERIFY(STATUS)

         call ESMF_AttributeGet  (FIELD, NAME="Needed", isPresent=isPresent, RC=STATUS)
         _VERIFY(STATUS)
         if(isPresent) then
            call ESMF_AttributeGet  (FIELD, NAME="Needed",VALUE=NEEDED, RC=STATUS)
            _VERIFY(STATUS)
         else
            NEEDED = .false.
         end if

         if( NEEDED .eqv. .false. ) then
            call ESMF_FieldGet(FIELD, Array=ARRAY, RC=STATUS)
            _VERIFY(STATUS)
            call ESMF_ArrayGet     (ARRAY, rank=RANK,   RC=STATUS)
            _VERIFY(STATUS)

            call ESMF_ArrayGet(array, localDeCount=localDeCount, rc=status)
            _VERIFY(STATUS)
            _ASSERT(localDeCount == 1, 'MAPL does not currently support multiple DEs per PET')
            call ESMF_ArrayGet(array, localarrayList=larrayList, rc=status)
            _VERIFY(STATUS)
            larray => lArrayList(1) ! alias

            select case (rank)
            case (1)
               call ESMF_LocalArrayGet(larray, PTR1, RC=status)
               _VERIFY(STATUS)
               if(associated(PTR1)) deallocate(PTR1)
            case (2)
               call ESMF_LocalArrayGet(larray, PTR2, RC=status)
               _VERIFY(STATUS)
               if(associated(PTR2)) deallocate(PTR2)
            case (3)
               call ESMF_LocalArrayGet(larray, PTR3, RC=status)
               _VERIFY(STATUS)
               if(associated(PTR3)) deallocate(PTR3)
            case (4)
               call ESMF_LocalArrayGet(larray, PTR4, RC=status)
               _VERIFY(STATUS)
               if(associated(PTR4)) deallocate(PTR4)
            end select
         end if
      end if
   end do

   deallocate(itemNameList,STAT=STATUS)
   _VERIFY(STATUS)
   deallocate(itemtypeList,STAT=STATUS)
   _VERIFY(STATUS)

   _RETURN(ESMF_SUCCESS)

 end subroutine ESMFL_StateFreePointers


 subroutine ESMFL_StateSetFieldNeeded(STATE, NAME, RC)
   type(ESMF_State),  intent(INOUT) :: STATE
   character(len=*),  intent(IN   ) :: NAME
   integer, optional, intent(  OUT) :: RC

   integer                          :: STATUS

   type(ESMF_Field)                 :: FIELD

   call ESMF_StateGet(STATE, trim(NAME), FIELD, RC=STATUS)
   _VERIFY(STATUS)

   call ESMF_AttributeSet  (FIELD, NAME="Needed",VALUE=.false., RC=STATUS)
   _VERIFY(STATUS)

   _RETURN(ESMF_SUCCESS)

 end subroutine ESMFL_StateSetFieldNeeded

function ESMFL_StateFieldIsNeeded(STATE, NAME, RC) result(NEEDED)
   type(ESMF_State),  intent(INOUT) :: STATE
   character(len=*),  intent(IN   ) :: NAME
   integer, optional, intent(  OUT) :: RC
   logical                          :: NEEDED

   integer                          :: STATUS

   type(ESMF_Field)                 :: FIELD

   call ESMF_StateGet(STATE, trim(NAME), FIELD, RC=STATUS)
   _VERIFY(STATUS)

   call ESMF_AttributeSet  (FIELD, NAME="Needed",VALUE=NEEDED, RC=STATUS)
   if(STATUS /= ESMF_SUCCESS) NEEDED = .false.

   _RETURN(ESMF_SUCCESS)

 end function ESMFL_StateFieldIsNeeded


#define VARTYPE_ 3

#define RANK_ 1
#include "GetFieldArray.H"
#define RANK_ 2
#include "GetFieldArray.H"
#define RANK_ 3
#include "GetFieldArray.H"
#define RANK_ 4
#include "GetFieldArray.H"

#undef VARTYPE_

#define VARTYPE_ 4

#define RANK_ 1
#include "GetFieldArray.H"
#define RANK_ 2
#include "GetFieldArray.H"
#define RANK_ 3
#include "GetFieldArray.H"
#define RANK_ 4
#include "GetFieldArray.H"

#undef VARTYPE_


#define VARTYPE_ 3

#define RANK_ 1
#include "GetPointer.H"
#define RANK_ 2
#include "GetPointer.H"
#define RANK_ 3
#include "GetPointer.H"
#define RANK_ 4
#include "GetPointer.H"

#undef VARTYPE_

#define VARTYPE_ 4

#define RANK_ 1
#include "GetPointer.H"
#define RANK_ 2
#include "GetPointer.H"
#define RANK_ 3
#include "GetPointer.H"
#define RANK_ 4
#include "GetPointer.H"

#undef VARTYPE_

 subroutine AdjustPtrBounds1dr4(PTR, A, I1, IN)
   real(KIND=ESMF_KIND_R4), pointer :: PTR(:)
   integer                          :: I1, IN
   real(KIND=ESMF_KIND_R4), target  :: A(I1:IN)

   ptr => A
 end subroutine AdjustPtrBounds1dr4

 subroutine AdjustPtrBounds1dr8(PTR, A, I1, IN)
   real(KIND=ESMF_KIND_R8), pointer :: PTR(:)
   integer                          :: I1, IN
   real(KIND=ESMF_KIND_R8), target  :: A(I1:IN)

   ptr => A
 end subroutine AdjustPtrBounds1dr8

 subroutine AdjustPtrBounds3dr4(PTR, A, I1, IN, J1, JN, L1, LN)
   real(KIND=ESMF_KIND_R4), pointer :: PTR(:,:,:)
   integer                          :: I1, IN, J1, JN, L1, LN
   real(KIND=ESMF_KIND_R4), target  :: A(I1:IN,J1:JN,L1:LN)


   ptr => A
 end subroutine AdjustPtrBounds3dr4

 subroutine AdjustPtrBounds3dr8(PTR, A, I1, IN, J1, JN, L1, LN)
   real(KIND=ESMF_KIND_R8), pointer :: PTR(:,:,:)
   integer                          :: I1, IN, J1, JN, L1, LN
   real(KIND=ESMF_KIND_R8), target  :: A(I1:IN,J1:JN,L1:LN)


   ptr => A
 end subroutine AdjustPtrBounds3dr8


   subroutine ESMFL_BundleGetPointerByIndex2(BUNDLE,INDEX,PTR,RC)
     type(ESMF_FieldBundle), intent(INOUT) :: BUNDLE !ALT: intent(in)
     integer,           intent(IN) :: INDEX
     real, pointer                 :: PTR(:,:)
     integer, optional, intent(OUT):: RC

     type(ESMF_FIELD)              :: FIELD
     integer                       :: status
     type(ESMF_FieldStatus_Flag)   :: fieldStatus

! ESMF 5 reorders items, be careful!
! ==================================

     call ESMF_FieldBundleGet (BUNDLE, INDEX, FIELD, RC=STATUS)
     _VERIFY(STATUS)

     call ESMF_FieldGet(field, status=fieldStatus, rc=status)
     _VERIFY(STATUS)
     if (fieldStatus == ESMF_FIELDSTATUS_COMPLETE) then
        call ESMF_FieldGet(field, 0, Ptr, rc=status)
        _VERIFY(STATUS)
     else
        NULLIFY(ptr)
     end if

     _RETURN(ESMF_SUCCESS)
   end subroutine ESMFL_BundleGetPointerByIndex2

   subroutine ESMFL_BundleGetPointerByIndex3(BUNDLE,INDEX,PTR,RC)
     type(ESMF_FieldBundle), intent(INOUT) :: BUNDLE !ALT: intent(in)
     integer,           intent(IN) :: INDEX
     real, pointer                 :: PTR(:,:,:)
     integer, optional, intent(OUT):: RC

     type(ESMF_FIELD)              :: FIELD
     type(ESMF_FieldStatus_Flag)   :: fieldStatus
     integer                       :: status

! ESMF 5 reorders items, be careful!
! ==================================

     call ESMF_FieldBundleGet (BUNDLE, INDEX, FIELD, RC=STATUS)
     _VERIFY(STATUS)
     call ESMF_FieldGet(field, status=fieldStatus, rc=status)
     _VERIFY(STATUS)
     if (fieldStatus == ESMF_FIELDSTATUS_COMPLETE) then
        call ESMF_FieldGet(field, 0, Ptr, rc=status)
        _VERIFY(STATUS)
     else
        NULLIFY(ptr)
     end if

     _RETURN(ESMF_SUCCESS)
   end subroutine ESMFL_BundleGetPointerByIndex3


   subroutine ESMFL_BundleGetPointerByName2(BUNDLE,NAME,PTR,RC)
     type(ESMF_FieldBundle), intent(INOUT) :: BUNDLE !ALT: intent(in)
     character(len=*), intent(IN) :: NAME
     real, pointer                 :: PTR(:,:)
     integer, optional, intent(OUT):: RC

     type(ESMF_FIELD)              :: FIELD
     type(ESMF_FieldStatus_Flag)   :: fieldStatus
     integer                       :: status

     call ESMF_FieldBundleGet (BUNDLE, NAME, FIELD=FIELD, RC=STATUS)
     _VERIFY(STATUS)
     call ESMF_FieldGet(field, status=fieldStatus, rc=status)
     _VERIFY(STATUS)
     if (fieldStatus == ESMF_FIELDSTATUS_COMPLETE) then
        call ESMF_FieldGet(field, 0, Ptr, rc=status)
        _VERIFY(STATUS)
     else
        NULLIFY(ptr)
     end if

     _RETURN(ESMF_SUCCESS)
   end subroutine ESMFL_BundleGetPointerByName2

   subroutine ESMFL_BundleGetPointerByName3(BUNDLE,NAME,PTR,RC)
     type(ESMF_FieldBundle), intent(INOUT) :: BUNDLE !ALT: intent(in)
     character(len=*), intent(IN) :: NAME
     real, pointer                 :: PTR(:,:,:)
     integer, optional, intent(OUT):: RC

     type(ESMF_FIELD)              :: FIELD
     type(ESMF_FieldStatus_Flag)   :: fieldStatus
     integer                       :: status

     call ESMF_FieldBundleGet (BUNDLE, NAME, FIELD=FIELD, RC=STATUS)
     _VERIFY(STATUS)
     call ESMF_FieldGet(field, status=fieldStatus, rc=status)
     _VERIFY(STATUS)
     if (fieldStatus == ESMF_FIELDSTATUS_COMPLETE) then
        call ESMF_FieldGet(field, 0, Ptr, rc=status)
        _VERIFY(STATUS)
     else
        NULLIFY(ptr)
     end if

     _RETURN(ESMF_SUCCESS)
   end subroutine ESMFL_BundleGetPointerByName3



    subroutine ESMFL_BundleCpyField (BUNDLE, FIELD, NAME, RC)
     type(ESMF_FieldBundle), intent(INOUT) :: BUNDLE
     type(ESMF_FIELD ), intent(INOUT) :: FIELD ! ALT: intent(in)
     character(len=ESMF_MAXSTR), optional, intent(IN) :: NAME
     integer, optional, intent(OUT)   :: RC

     type(ESMF_FIELD )       :: FIELD1


     integer          :: status

     FIELD1 = MAPL_FieldCreate(FIELD, name = name, RC=STATUS )
     _VERIFY(STATUS)

     call MAPL_FieldBundleAdd (BUNDLE, FIELD1, RC=STATUS )
     _VERIFY(STATUS)

     _RETURN(ESMF_SUCCESS)

   end subroutine  ESMFL_BundleCpyField


#ifdef DO_WE_NEED_THIS
 subroutine ESMFL_Connect2STATE(STATE, FIELD, RC)
   type(ESMF_State),  intent(INOUT) :: STATE
   type(ESMF_Field),  intent(INOUT) :: FIELD ! ALT: intent(in)
   integer, optional, intent(  OUT) :: rc

! local vars
   integer          :: status
   character(len=ESMF_MAXSTR)     :: NAME
   type (ESMF_Field)              :: f

   call ESMF_FieldGet(FIELD, NAME=NAME, RC=STATUS)
   _VERIFY(STATUS)

   call ESMF_StateGet(STATE, NAME, F, RC=STATUS)
   _VERIFY(STATUS)

   call MAPL_ConnectCoupling(F, FIELD, RC=STATUS)
   _VERIFY(STATUS)

   _RETURN(ESMF_SUCCESS)
 end subroutine ESMFL_Connect2STATE
#endif

 subroutine ESMFL_FCOLLECT_I4(GRID, FULLINPUT, INPUT, RC )
   type(ESMF_GRID), intent(IN   )     :: GRID
   integer,             intent(INOUT) :: FULLINPUT(:)
   integer,             intent(IN   ) :: INPUT(:)
   integer, optional,   intent(  OUT) :: rc

! local vars
   integer          :: status

   type (ESMF_DistGrid)                  :: distGrid
   type(ESMF_DELayout)                   :: LAYOUT
   type (ESMF_VM)                        :: vm
   integer,               allocatable    :: AL(:,:)
   integer,               allocatable    :: AU(:,:)

   integer, pointer, dimension(:)            :: recvcounts, displs
   integer                                       :: nDEs
   integer                                       :: sendcount

   integer                                       :: I, J, de, deId
   integer                                       :: I1, IN
   integer                                       :: gridRank


   call ESMF_GridGet(GRID, dimCount=gridRank, rc=STATUS)
   _VERIFY(STATUS)

   call ESMF_GridGet    (GRID,   distGrid=distGrid, rc=STATUS)
   _VERIFY(STATUS)
   call ESMF_DistGridGet(distGRID, delayout=layout, rc=STATUS)
   _VERIFY(STATUS)
   call ESMF_DELayoutGet(layout, vm=vm, rc=status)
   _VERIFY(STATUS)
   call ESMF_VmGet(vm, localPet=deId, petCount=nDEs, rc=status)
   _VERIFY(STATUS)

   allocate (AL(gridRank,0:nDEs-1),  stat=status)
   _VERIFY(STATUS)
   allocate (AU(gridRank,0:nDEs-1),  stat=status)
   _VERIFY(STATUS)

   call MAPL_DistGridGet(distgrid, &
        minIndex=AL, maxIndex=AU, rc=status)
   _VERIFY(STATUS)

   allocate (recvcounts(nDEs), displs(0:nDEs), stat=status)
   _VERIFY(STATUS)

   displs(0) = 0
   do I = 1,nDEs
      J = I - 1
      de = J
      I1 = AL(1,J)
      IN = AU(1,J)
      recvcounts(I) = (IN - I1 + 1)
      if (de == deId) then
         sendcount = recvcounts(I)
      endif
      displs(I) = displs(J) + recvcounts(I)
   enddo


   _ASSERT(sendcount == size(input), 'inconsistent sendcount')
   _ASSERT(sum(recvcounts) == size(fullinput), 'inconsistent recvcount')

   call MAPL_CommsAllGatherV(layout, input, sendcount, &
                             fullinput, recvcounts, displs, rc=status)
   _VERIFY(STATUS)

   deallocate(recvcounts, displs, AU, AL, stat=status)
   _VERIFY(STATUS)

   _RETURN(ESMF_SUCCESS)
 end subroutine ESMFL_FCOLLECT_I4

 subroutine ESMFL_FCOLLECT_R4(GRID, FULLINPUT, INPUT, RC )
   type(ESMF_GRID), intent(IN   )     :: GRID
   real,             intent(INOUT) :: FULLINPUT(:)
   real,             intent(IN   ) :: INPUT(:)
   integer, optional,   intent(  OUT) :: rc

! local vars
   integer          :: status

   type (ESMF_DistGrid)                  :: distGrid
   type(ESMF_DELayout)                   :: LAYOUT
   type (ESMF_VM)                        :: vm
   integer,               allocatable    :: AL(:,:)
   integer,               allocatable    :: AU(:,:)

   integer, pointer, dimension(:)            :: recvcounts, displs
   integer                                       :: nDEs
   integer                                       :: sendcount

   integer                                       :: I, J, de, deId
   integer                                       :: I1, IN
   integer                                       :: gridRank


   call ESMF_GridGet(GRID, dimCount=gridRank, rc=STATUS)
   _VERIFY(STATUS)

   call ESMF_GridGet    (GRID,   distGrid=distGrid, rc=STATUS)
   _VERIFY(STATUS)
   call ESMF_DistGridGet(distGRID, delayout=layout, rc=STATUS)
   _VERIFY(STATUS)
   call ESMF_DELayoutGet(layout, vm=vm, rc=status)
   _VERIFY(STATUS)
   call ESMF_VmGet(vm, localPet=deId, petCount=nDEs, rc=status)
   _VERIFY(STATUS)

   allocate (AL(gridRank,0:nDEs-1),  stat=status)
   _VERIFY(STATUS)
   allocate (AU(gridRank,0:nDEs-1),  stat=status)
   _VERIFY(STATUS)

   call MAPL_DistGridGet(distgrid, &
        minIndex=AL, maxIndex=AU, rc=status)
   _VERIFY(STATUS)

   allocate (recvcounts(nDEs), displs(0:nDEs), stat=status)
   _VERIFY(STATUS)

   displs(0) = 0
   do I = 1,nDEs
      J = I - 1
      de = J
      I1 = AL(1,J)
      IN = AU(1,J)
      recvcounts(I) = (IN - I1 + 1)
      if (de == deId) then
         sendcount = recvcounts(I)
      endif
      displs(I) = displs(J) + recvcounts(I)
   enddo

   _ASSERT(sendcount == size(input), 'inconsistent sendcount')
   _ASSERT(sum(recvcounts) == size(fullinput), 'inconsistent recvcount')

   call MAPL_CommsAllGatherV(layout, input, sendcount, &
                             fullinput, recvcounts, displs, rc=status)
   _VERIFY(STATUS)

   deallocate(recvcounts, displs, AU, AL, stat=status)
   _VERIFY(STATUS)

   _RETURN(ESMF_SUCCESS)
 end subroutine ESMFL_FCOLLECT_R4

 subroutine ESMFL_FCOLLECT_R8(GRID, FULLINPUT, INPUT, RC )
   type(ESMF_GRID), intent(IN   )     :: GRID
   real(kind= ESMF_KIND_R8), intent(INOUT) :: FULLINPUT(:)
   real(kind= ESMF_KIND_R8), intent(IN   ) :: INPUT(:)
   integer, optional,   intent(  OUT) :: rc

! local vars
   integer          :: status

   type (ESMF_DistGrid)                  :: distGrid
   type(ESMF_DELayout)                   :: LAYOUT
   type (ESMF_VM)                        :: vm
   integer,               allocatable    :: AL(:,:)
   integer,               allocatable    :: AU(:,:)

   integer, pointer, dimension(:)            :: recvcounts, displs
   integer                                       :: nDEs
   integer                                       :: sendcount

   integer                                       :: I, J, de, deId
   integer                                       :: I1, IN
   integer                                       :: gridRank


   call ESMF_GridGet(GRID, dimCount=gridRank, rc=STATUS)
   _VERIFY(STATUS)

   call ESMF_GridGet    (GRID,   distGrid=distGrid, rc=STATUS)
   _VERIFY(STATUS)
   call ESMF_DistGridGet(distGRID, delayout=layout, rc=STATUS)
   _VERIFY(STATUS)
   call ESMF_DELayoutGet(layout, vm=vm, rc=status)
   _VERIFY(STATUS)
   call ESMF_VmGet(vm, localPet=deId, petCount=nDEs, rc=status)
   _VERIFY(STATUS)

   allocate (AL(gridRank,0:nDEs-1),  stat=status)
   _VERIFY(STATUS)
   allocate (AU(gridRank,0:nDEs-1),  stat=status)
   _VERIFY(STATUS)

   call MAPL_DistGridGet(distgrid, &
        minIndex=AL, maxIndex=AU, rc=status)
   _VERIFY(STATUS)

   allocate (recvcounts(nDEs), displs(0:nDEs), stat=status)
   _VERIFY(STATUS)

   displs(0) = 0
   do I = 1,nDEs
      J = I - 1
      de = J
      I1 = AL(1,J)
      IN = AU(1,J)
      recvcounts(I) = (IN - I1 + 1)
      if (de == deId) then
         sendcount = recvcounts(I)
      endif
      displs(I) = displs(J) + recvcounts(I)
   enddo

   _ASSERT(sendcount == size(input), 'inconsistent sendcount')
   _ASSERT(sum(recvcounts) == size(fullinput), 'inconsistent recvcount')

   call MAPL_CommsAllGatherV(layout, input, sendcount, &
                             fullinput, recvcounts, displs, rc=status)
   _VERIFY(STATUS)

   deallocate(recvcounts, displs, AU, AL, stat=status)
   _VERIFY(STATUS)

   _RETURN(ESMF_SUCCESS)
 end subroutine ESMFL_FCOLLECT_R8

 subroutine ESMFL_FieldRegrid(src, dst, RC)
   type(ESMF_Field) :: src
   type(ESMF_Field) :: dst
   integer, optional, intent(OUT) :: rc


! local vars
#if 0
   type (ESMF_Grid) :: srcgrid, dstgrid
   type (ESMF_DELayout) :: layout
   type (ESMF_Array)    :: array
   integer              :: DIMS(3)
   integer              :: IM_SRC, JM_SRC
   integer              :: IM_DST, JM_DST
   integer              :: J
   real, dimension(:,:), pointer :: lats, lons, z0, z, h
   real, parameter               :: EPS=1.0E-3
   real                          :: ZPOLE
#endif

   _UNUSED_DUMMY(src)
   _UNUSED_DUMMY(dst)

#if 0
! begin
   call ESMF_FieldGetGrid(src, srcgrid, rc=status)
   _VERIFY(STATUS)

   call ESMF_GridGet(SRCGRID, global_cell_dim=dims, RC=STATUS)
   _VERIFY(STATUS)

   IM_SRC = DIMS(1)
   JM_SRC = DIMS(2)

   call ESMF_GridGetDELayout(srcgrid, layout=layout, rc=status)
   _VERIFY(STATUS)

   call ESMF_FieldGetGrid(dst, dstgrid, rc=status)
   _VERIFY(STATUS)

   call ESMF_GridGetDELocalInfo(DSTGRID, LOCAL_AXIS_LENGTH=DIMS, RC=STATUS)
   _VERIFY(STATUS)

   IM_DST = DIMS(1)
   JM_DST = DIMS(2)

   call ESMFL_GridCoordGet(   DSTGRID, LATS                   , &
                              Name     = "Latitude"              , &
                              Location = ESMF_CELL_CENTER        , &
                              Units    = MAPL_UnitsRadians      , &
                              RC       = STATUS                    )
   _VERIFY(STATUS)

   call ESMFL_GridCoordGet(   DSTGRID, LONS                   , &
                              Name     = "Longitude"             , &
                              Location = ESMF_CELL_CENTER        , &
                              Units    = MAPL_UnitsRadians      , &
                              RC       = STATUS                    )
   _VERIFY(STATUS)



   call ESMF_FieldAllGather(src, array=array, rc=status)
   _VERIFY(STATUS)

   call ESMF_ArrayGetData(array, z0, RC=status)
   _VERIFY(STATUS)

   call ESMF_FieldGetData(dst, array, rc=status)
   _VERIFY(STATUS)

   call ESMF_ArrayGetData(array, z, RC=status)
   _VERIFY(STATUS)

! binning

   call topo_bin (z0, im_src, jm_src, z, lons(:,1), lats(1,:), &
                  im_dst,jm_dst)
!   _VERIFY(STATUS)


! do I have the pole?
   call ESMF_FieldAllGather(dst, array=array, rc=status)
   _VERIFY(STATUS)

   call ESMF_ArrayGetData(array, h, RC=status)
   _VERIFY(STATUS)

   DO J = 1, JM_DST
      IF (ABS(ABS(LATS(1,J)) - 0.5*PI) .LT. EPS) then
! yes, I have the pole
         ZPOLE = SUM(h(:,J))/SIZE(H,1)
         Z(:,J) = ZPOLE
      end IF
   END DO
#endif

   _RETURN(ESMF_SUCCESS)
 end subroutine ESMFL_FieldRegrid

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! Given a `srcFLD` and its associated `3dGrid` and a `dstFLD` and its associated
! `3DGrid`, the subroutine `ESMFL_RegridStore` creates their corresponding
! `2DGrids` and a 2D routehandle.
!
!#### History
!- 17Oct2005  Cruz  Initial code.
!
   subroutine ESMFL_RegridStore (srcFLD, SRCgrid2D, dstFLD, DSTgrid2D, &
                                 vm, rh, rc)
!
   implicit NONE

! !ARGUMENTS:

   type(ESMF_Field), intent(inout)       :: srcFLD
   type(ESMF_Field), intent(inout)       :: dstFLD
   type(ESMF_Grid), intent(out)          :: SRCgrid2D
   type(ESMF_Grid), intent(out)          :: DSTgrid2D
   type(ESMF_RouteHandle), intent(inout) :: rh
   type(ESMF_VM), intent(in)             :: vm  ! should be intent IN
   integer, optional, intent(OUT)        :: rc
!
!-------------------------------------------------------------------------

#if 0
! local vars

   type(ESMF_DELayout)     :: layout
   type(ESMF_Field)        :: dstFld2D
   type(ESMF_Field)        :: srcFld2D
   type(ESMF_ArraySpec)    :: arrayspec
   type(ESMF_Array)        :: dstARR
   type(ESMF_Array)        :: srcARR
   type(ESMF_FieldDataMap) :: dmap
   type(ESMF_Grid)         :: grid3D
   real(kind=REAL32), pointer        :: Sptr2d(:,:)
   real(kind=REAL32), pointer        :: Dptr2d(:,:)
   real(ESMF_KIND_R8)      :: deltaX, deltaY
   real(ESMF_KIND_R8)      :: min(2), max(2)
   integer :: status, rank
   integer :: sCPD(3), dCPD(3) ! localCellCountPerDim
   integer :: gccpd(3)         ! globalCellCountPerDim

   type(ESMF_AxisIndex), dimension(:,:), pointer :: AI
   integer, allocatable, dimension(:)            :: ims, jms
   integer                                   :: nDEs, de
   integer                                   :: NX, NY
   integer                                   :: IX, JY
   integer                                   :: gridRank
   integer                                   :: decpd(7)
#endif

! start

   _UNUSED_DUMMY(srcFLD)
   _UNUSED_DUMMY(SRCgrid2D)
   _UNUSED_DUMMY(dstFLD)
   _UNUSED_DUMMY(DSTgrid2d)
   _UNUSED_DUMMY(rh)
   _UNUSED_DUMMY(vm)

#if 0
   ! can get the rank from either srcFLD or dstFLD
   call ESMF_FieldGetArray(srcFLD, srcARR, rc=status)
   _VERIFY(status)
   call ESMF_ArrayGet(srcARR, RANK=rank, rc=status)
   _VERIFY(status)
   ! datamap - with rank=2!
   call ESMF_FieldDataMapSetDefault(dmap, 2, rc=status)
   _VERIFY(STATUS)

   ! Create 2D grids and fields

   ! SOURCE FIELD

   call ESMF_FieldGet(srcFLD, grid=grid3D, rc=status)
   _VERIFY(status)

   call ESMF_GridGet(GRID3d, dimCount=gridRank, rc=STATUS)
   _VERIFY(STATUS)

   call ESMF_GridGetDELayout(grid3D, layout, rc=status)
   _VERIFY(STATUS)

   call ESMF_DELayoutGet(layout, deCount=nDEs, deCountPerDim = decpd, rc=status)

   NX = decpd(1)
   NY = decpd(2)

   allocate (ims(NX), jms(NY), stat=status)
   _VERIFY(STATUS)

   call ESMF_GridGet(grid3D, &
                     horzRelLoc=ESMF_CELL_CENTER, vertRelLoc=ESMF_CELL_CELL, &
                     globalCellCountPerDim=gccpd, &
                     minGlobalCoordPerDim=min, &
                     maxGlobalCoordPerDim=max, &
                     rc=status)
   _VERIFY(status)

#if 0
   ! kludge min,max DEGREE/RADIANS conversion
   if (min(1) < 0.0) then
      min(1) = min(1) * 180./MAPL_PI
      min(2) = min(2) * 180./MAPL_PI
      max(1) = max(1) * 180./MAPL_PI
      max(2) = max(2) * 180./MAPL_PI
   end if
   ! overide min,max - just create a simple 2-D grid
   min(1) = 0.0
   min(2) = 0.0
   max(1) = 360.0
   max(2) = 180.0
#endif

   ! this is probably incorrect unless parent grid is XYuni
   !SRCGrid2D = ESMF_GridCreateHorzXYUni( &
   !                           counts=gccpd(1:2), &
   !                           minGlobalCoordPerDim=min, &
   !                           maxGlobalCoordPerDim=max, &
   !                           horzStagger=ESMF_GRID_HORZ_STAGGER_A, &
   !                           periodic=(/ESMF_TRUE, ESMF_FALSE/),     &
   !                           name="SRC 2D grid", rc=status)
   !_VERIFY(status)
   ! instead use the following...
   deltaX = 2.0*MAPL_PI/gccpd(1)
   deltaY = MAPL_PI/(gccpd(2)-1)
   SRCGrid2D = ESMF_GridCreateHorzLatLonUni(         &
         counts = gccpd(1:2),        &
         minGlobalCoordPerDim=min(1:2),     &
         deltaPerDim=(/deltaX, deltaY /),        &
         horzStagger=ESMF_Grid_Horz_Stagger_A,   &
         periodic=(/ESMF_TRUE, ESMF_FALSE/),     &
         name='SRC 2D grid', rc=status)
   _VERIFY(STATUS)

   allocate (AI(nDEs,gridRank), stat=status)
   _VERIFY(STATUS)

   if (gridRank == 3) then
      call ESMF_GridGetAllAxisIndex(grid3d, &
                                    horzRelLoc=ESMF_CELL_CENTER, &
                                    vertRelLoc=ESMF_CELL_CELL, &
                                    globalAI=AI, rc=status)
   else
      call ESMF_GridGetAllAxisIndex(grid3d, &
                                    horzRelLoc=ESMF_CELL_CENTER, &
                                    globalAI=AI, rc=status)
   end if
   _VERIFY(STATUS)

   JY = 1
   DO IX = 1, NX

!ALT: this is a workaround to compute deId from Layout coords
      de = (jy-1)*NX + ix

      ims(IX) = AI(de,1)%max - AI(de,1)%min + 1
   END DO

   IX = 1
   DO JY = 1, NY

!ALT: same workaround as above
      de = (jy-1)*NX + ix

      jms(JY) = AI(de,2)%max - AI(de,2)%min + 1
   END DO

   deallocate(AI)

   if (verbose .and. MAPL_AM_I_Root()) then
      print *, 'ims=', ims
      print *, 'jms=', jms
   end if

   call ESMF_GridDistribute(SRCGrid2D, delayout=layout, &
                            countsPerDEDim1=ims,        &
                            countsPerDEDim2=jms,        &
                            rc=status)
   _VERIFY(status)

   deallocate(jms, ims)

   call ESMF_GridGetDELocalInfo(SRCGrid2D, &
         horzRelLoc=ESMF_CELL_CENTER, &
!         vertRelLoc=ESMF_CELL_CELL, &
         localCellCountPerDim=sCPD,RC=STATUS)
   _VERIFY(STATUS)
   allocate(Sptr2d(1:sCPD(1),1:sCPD(2)), STAT=STATUS)
   srcArr = ESMF_ArrayCreate(Sptr2d, ESMF_DATA_REF, RC=STATUS)
   _VERIFY(STATUS)
   srcFLD2D = ESMF_FieldCreate(SRCGrid2D, srcARR,    &
                 horzRelloc = ESMF_CELL_CENTER,         &
                 datamap=dmap,                       &
                 haloWidth=0, &
                 name   = "PS", rc=status )  ! can be any name, all we want
   _VERIFY(STATUS)                           ! is the route handle

   ! DESTINATION FIELD

   call ESMF_FieldGet(dstFLD, grid=grid3D, rc=status)
   _VERIFY(status)

   call ESMF_GridGetDELayout(grid3D, layout, rc=status)
   _VERIFY(STATUS)

   call ESMF_DELayoutGet(layout, deCount=nDEs, deCountPerDim = decpd, rc=status)

   NX = decpd(1)
   NY = decpd(2)

   allocate (ims(NX), jms(NY), stat=status)
   _VERIFY(STATUS)

   call ESMF_GridGet(grid3D, &
                     horzRelLoc=ESMF_CELL_CENTER, vertRelLoc=ESMF_CELL_CELL, &
                     globalCellCountPerDim=gccpd, &
                     minGlobalCoordPerDim=min, &
                     maxGlobalCoordPerDim=max, &
                     rc=status)
   _VERIFY(status)

   ! this is probably incorrect unless parent grid is XYuni
   !DSTGrid2D = ESMF_GridCreateHorzXYUni( &
   !                           counts=gccpd(1:2), &
   !                           minGlobalCoordPerDim=min, &
   !                           maxGlobalCoordPerDim=max, &
   !                           horzStagger=ESMF_GRID_HORZ_STAGGER_A, &
   !                           periodic=(/ESMF_TRUE, ESMF_FALSE/),     &
   !                           name="DST 2D grid", rc=status)
   !_VERIFY(status)
   ! instead use the following ...
   deltaX = 2.0*MAPL_PI/gccpd(1)
   deltaY = MAPL_PI/(gccpd(2)-1)
   DSTGrid2D = ESMF_GridCreateHorzLatLonUni(         &
         counts = gccpd(1:2),        &
         minGlobalCoordPerDim=min(1:2),     &
         deltaPerDim=(/deltaX, deltaY /),        &
         horzStagger=ESMF_Grid_Horz_Stagger_A,   &
         periodic=(/ESMF_TRUE, ESMF_FALSE/),     &
         name='DST 2D grid', rc=status)
   _VERIFY(STATUS)

   allocate (AI(nDEs,gridRank), stat=status)
   _VERIFY(STATUS)

   if (gridRank == 3) then
      call ESMF_GridGetAllAxisIndex(grid3d, &
                                    horzRelLoc=ESMF_CELL_CENTER, &
                                    vertRelLoc=ESMF_CELL_CELL, &
                                    globalAI=AI, rc=status)
   else
      call ESMF_GridGetAllAxisIndex(grid3d, &
                                    horzRelLoc=ESMF_CELL_CENTER, &
                                    globalAI=AI, rc=status)
   end if
   _VERIFY(STATUS)

   JY = 1
   DO IX = 1, NX

!ALT: this is a workaround to compute deId from Layout coords
      de = (jy-1)*NX + ix

      ims(IX) = AI(de,1)%max - AI(de,1)%min + 1
   END DO

   IX = 1
   DO JY = 1, NY

!ALT: same workaround as above
      de = (jy-1)*NX + ix

      jms(JY) = AI(de,2)%max - AI(de,2)%min + 1
   END DO

   deallocate(AI)

   if (verbose .and. MAPL_AM_I_Root()) then
      print *, 'ims=', ims
      print *, 'jms=', jms
   end if

   call ESMF_GridDistribute(DSTGrid2D, delayout=layout, &
                            countsPerDEDim1=ims,        &
                            countsPerDEDim2=jms,        &
                            rc=status)
   _VERIFY(status)

   deallocate(jms, ims)

   call ESMF_GridGetDELocalInfo(DSTGrid2D, &
         horzRelLoc=ESMF_CELL_CENTER, &
!         vertRelLoc=ESMF_CELL_CELL, &
         localCellCountPerDim=dCPD,RC=STATUS)
   _VERIFY(STATUS)
   allocate(Dptr2d(1:dCPD(1),1:dCPD(2)), STAT=STATUS)
   _VERIFY(STATUS)
   dstArr = ESMF_ArrayCreate(Dptr2d, ESMF_DATA_REF, RC=STATUS)
   _VERIFY(STATUS)
   dstFLD2D = ESMF_FieldCreate(DSTGrid2D, dstARR,    &
                 horzRelloc = ESMF_CELL_CENTER,         &
                 datamap=dmap,                       &
                 haloWidth=0, &
                 name   = "PS", rc=status )  ! can be any name, all we want
   _VERIFY(STATUS)                           ! is the correct route handle

   call ESMF_FIELDRegridStore(srcFLD2D, dstFLD2D, vm, rh,  &
                             regridmethod=ESMF_REGRID_METHOD_BILINEAR, &
                             rc=status)
   _VERIFY(status)

   deallocate(Sptr2d, stat=status)
   _VERIFY(STATUS)
   deallocate(Dptr2d, stat=status)
   _VERIFY(STATUS)
#endif

   if (present(rc)) then
      rc = 0
   end if

   end subroutine ESMFL_RegridStore

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `FieldRegrid1` regrids 3D fields using ESMF_FieldRegrid.
!
!#### History
!- 17Oct2005  Cruz  Initial code.
!
   subroutine FieldRegrid1 (srcFLD, Sgrid2D, dstFLD, Dgrid2D, &
                                   vm, rh, fname, rc)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_Field), intent(in)          :: srcFLD
   type(ESMF_Field), intent(inout)       :: dstFLD
   type(ESMF_Grid), intent(in)           :: Sgrid2D
   type(ESMF_Grid), intent(in)           :: Dgrid2D
   type(ESMF_RouteHandle), intent(inout) :: rh
   type(ESMF_VM), intent(inout)          :: vm
   ! assumes name of src and dst are the same!!
   character(len=*), intent(in)          :: fname
   integer, optional, intent(out)        :: rc
!
!-------------------------------------------------------------------------

! local vars
   _UNUSED_DUMMY(srcFLD)
   _UNUSED_DUMMY(dstFLD)
   _UNUSED_DUMMY(Sgrid2D)
   _UNUSED_DUMMY(Dgrid2D)
   _UNUSED_DUMMY(vm)
   _UNUSED_DUMMY(rh)
   _UNUSED_DUMMY(fname)

#if 0
   type(ESMF_Field)        :: dstFLD2D
   type(ESMF_Field)        :: srcFLD2D
   type(ESMF_Array)        :: srcARR
   type(ESMF_Array)        :: dstARR
   type(ESMF_Array)        :: newSrcARR
   type(ESMF_Array)        :: newDstARR
   type(ESMF_Grid)         :: grid3D
   type(ESMF_FieldDataMap) :: dmap
   real(kind=REAL32), dimension(:,:), pointer   :: sptr, dptr
   real(kind=REAL32), dimension(:,:), pointer   :: sptr2d, dptr2d
   real(kind=REAL32), dimension(:,:,:), pointer :: sptr3d, dptr3d
   integer :: sCPD(3), dCPD(3) ! localCellCountPerDim
   integer :: rank, k, deid, kmax, status

! begin

   call ESMF_VMGet(vm, localPet=deid, rc=status)
   _VERIFY(status)

   ! get rank (get from any FLD)
   call ESMF_FieldGetArray(dstFLD, dstARR, rc=status)
   _VERIFY(status)
   call ESMF_ArrayGet(dstARR, RANK=rank, rc=status)
   _VERIFY(status)

   ! datamap - with rank=2!
   call ESMF_FieldDataMapSetDefault(dmap, 2, rc=status)
   _VERIFY(STATUS)

   ! get sCPD from srcFLD, Sgrid2D does not have correct 3rd dim
   call ESMF_FieldGet(srcFLD, grid=grid3D, rc=status)
   _VERIFY(STATUS)
   call ESMF_GridGetDELocalInfo(grid3D, &
         horzRelLoc=ESMF_CELL_CENTER, &
         vertRelLoc=ESMF_CELL_CELL, &
         localCellCountPerDim=sCPD,RC=STATUS)
   _VERIFY(STATUS)
   ! get dCPD from Dgrid2D
   call ESMF_FieldGet(dstFLD, grid=grid3D, rc=status)
   _VERIFY(STATUS)
   call ESMF_GridGetDELocalInfo(grid3D, &
         horzRelLoc=ESMF_CELL_CENTER, &
         vertRelLoc=ESMF_CELL_CELL, &
         localCellCountPerDim=dCPD,RC=STATUS)
   _VERIFY(STATUS)

   ! for MAPL and GSI sCPD(3) = dCPD(3)

   ! assume 2D field
   kmax=1
   ! if rank is three then km = # levels in srcFLD
   if (rank==3) kmax=sCPD(3)

   ! get array from srcFLD and dstFLD
   call ESMF_FieldGetArray(srcFLD, srcARR, rc=status)
   _VERIFY(status)
   call ESMF_FieldGetArray(dstFLD, dstARR, rc=status)
   _VERIFY(status)

   ! allocate f90 pointer to hold data
   ! note halo width = 0
   if(rank==3) then
     call ESMF_ArrayGetData(srcARR, sptr3d, rc=status)
     _VERIFY(status)
     call ESMF_ArrayGetData(dstARR, dptr3d, rc=status)
     _VERIFY(status)
   else
     call ESMF_ArrayGetData(srcARR, sptr2d, rc=status)
     _VERIFY(status)
     call ESMF_ArrayGetData(dstARR, dptr2d, rc=status)
     _VERIFY(status)
   end if

   ! these are 2d pointers used to create 2d fields
   allocate(sptr(1:sCPD(1),1:sCPD(2)), STAT=STATUS)
   _VERIFY(status)
   allocate(dptr(1:dCPD(1),1:dCPD(2)), STAT=STATUS)
   _VERIFY(status)
   ! 2d ESMF arrays associated with pointers sptr and dptr
   newSrcARR = ESMF_ArrayCreate(sptr, ESMF_DATA_REF, RC=STATUS)
   _VERIFY(STATUS)
   newDstARR = ESMF_ArrayCreate(dptr, ESMF_DATA_REF, RC=STATUS)
   _VERIFY(STATUS)

   ! local fields built with 2d sized arrays and 2D grids
   srcFLD2D = ESMF_FieldCreate(Sgrid2D, newSrcARR,    &
                 horzRelloc = ESMF_CELL_CENTER,         &
!                 vertRelloc = ESMF_CELL_CELL,         &
                 datamap=dmap,                       &
                 haloWidth=0, &
                 name   = fname, rc=status )
   _VERIFY(STATUS)
   dstFLD2D = ESMF_FieldCreate(Dgrid2D, newDstARR,    &
                 horzRelloc = ESMF_CELL_CENTER,         &
!                 vertRelloc = ESMF_CELL_CELL,         &
                 datamap=dmap,                       &
                 haloWidth=0, &
                 name   = fname, rc=status )
   _VERIFY(STATUS)

   ! loop over field's vertical levels
   do k=1, kmax

     ! modify field's pointer by reference
     if (rank==3) then
       sptr(:,:) = sptr3d(:,:,k)
     else
       sptr(:,:) = sptr2d(:,:)
     end if

     call ESMF_FieldRegrid(srcFLD2D, dstFLD2D, rh, rc=status)
     _VERIFY(STATUS)

     !  Update contents (array pointer) of dstFLD
     call ESMF_FieldGetArray(dstFLD2D, dstARR, rc=status)
     _VERIFY(status)
     call ESMF_ArrayGetData(dstARR, dptr, rc=status)
     _VERIFY(status)

     if (rank==3) then
       dptr3d(:,:,k) = dptr(:,:)
     else
       dptr2d(:,:) = dptr(:,:)
     end if

   end do

   deallocate(sptr, dptr, stat=status)
   _VERIFY(STATUS)
#endif

   if (present(rc)) then
      rc = 0
   end if

   end subroutine FieldRegrid1

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `BundleRegrid1`
! regrids members of a bundle using ESMF_FieldRegrid.
!
!#### History
!- 17Apr2006  Cruz  Initial code.
!
   subroutine BundleRegrid1 (srcBUN, Sgrid2D, dstBUN, Dgrid2D, &
                             vm, rh, rc)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_FieldBundle), intent(inOUT)      :: srcBUN
   type(ESMF_FieldBundle), intent(inout)      :: dstBUN
   type(ESMF_Grid), intent(in)           :: Sgrid2D
   type(ESMF_Grid), intent(in)           :: Dgrid2D
   type(ESMF_RouteHandle), intent(inout) :: rh
   type(ESMF_VM), intent(inout)          :: vm
   integer, optional, intent(out)        :: rc
!
!-------------------------------------------------------------------------

! local vars
   type(ESMF_Field)        :: dstFLD
   type(ESMF_Field)        :: srcFLD
   type(ESMF_Array)        :: dstARR
   type(ESMF_Grid)         :: grid3D
   integer :: sCPD(3), dCPD(3)
   integer :: rank, deid, kmax, status, numVars

! begin

   _UNUSED_DUMMY(Sgrid2D)
   _UNUSED_DUMMY(Dgrid2D)
   _UNUSED_DUMMY(rh)

   call ESMF_VMGet(vm, localPet=deid, rc=status)
   _VERIFY(status)

   ! get number of fields in bundle
   ! number in srcBUN should be the same as in dstBUN (we can change later)
   call ESMF_FieldBundleGet      (srcBUN, FieldCount=NumVars,  RC=STATUS)
   _VERIFY(STATUS)

   call ESMF_FieldBundleGet (srcBUN,  1, srcFLD,            RC=STATUS)
   _VERIFY(STATUS)
   call ESMF_FieldBundleGet (dstBUN,  1, dstFLD,            RC=STATUS)
   _VERIFY(STATUS)

   ! get rank from any field
   call ESMF_FieldGet(dstFLD, array=dstARR, rc=status)
   _VERIFY(status)
   call ESMF_ArrayGet(dstARR, RANK=rank, rc=status)
   _VERIFY(status)

   ! get sCPD from srcFLD, Sgrid2D does not have correct 3rd dim
   call ESMF_FieldGet(srcFLD, grid=grid3D, rc=status)
   _VERIFY(STATUS)

   call MAPL_GridGet(grid3D, &
         localCellCountPerDim=sCPD, &
         RC=STATUS)
   _VERIFY(STATUS)
   ! get dCPD from Dgrid2D
   call ESMF_FieldGet(dstFLD, grid=grid3D, rc=status)
   _VERIFY(STATUS)
   call MAPL_GridGet(grid3D, &
         localCellCountPerDim=dCPD,RC=STATUS)
   _VERIFY(STATUS)

   ! for MAPL and GSI sCPD(3) = dCPD(3)

   ! assume 2D field
   kmax=1
   ! if rank is three then km = # levels in srcFLD
   if (rank==3) kmax=sCPD(3)


#if 0
   ! loop over fields in bundle

   do L=1,NumVars

     call ESMF_FieldBundleGet     (srcBUN, L, srcFLD,           RC=STATUS)
     _VERIFY(STATUS)
     call ESMF_FieldGet           (srcFLD, NAME=srcName,         RC=STATUS)
     _VERIFY(STATUS)
     call ESMF_FieldBundleGet     (dstBUN, L, dstFLD,           RC=STATUS)
     _VERIFY(STATUS)
     call ESMF_FieldGet           (dstFLD, NAME=dstName,         RC=STATUS)
     _VERIFY(STATUS)

     ! allocate f90 pointer to hold data
     ! note halo width = 0
     if(rank==3) then
       call ESMF_FieldGet (srcFLD, 0, sptr3d, rc = status)
       _VERIFY(status)
       call ESMF_FieldGet(dstFLD, dptr3d, rc = status)
       _VERIFY(status)
     else
       call ESMF_FieldGetDataPointer (srcFLD, sptr2d, rc = status)
       _VERIFY(status)
       call ESMF_FieldGetDataPointer (dstFLD, dptr2d, rc = status)
       _VERIFY(status)
     end if

     ! these are 2d pointers used to create 2d fields
     allocate(sptr(1:sCPD(1),1:sCPD(2)), STAT=STATUS)
     _VERIFY(status)
     allocate(dptr(1:dCPD(1),1:dCPD(2)), STAT=STATUS)
     _VERIFY(status)
     ! 2d ESMF arrays associated with pointers sptr and dptr
     newSrcARR = ESMF_ArrayCreate(sptr, ESMF_DATA_REF, RC=STATUS)
     _VERIFY(STATUS)
     newDstARR = ESMF_ArrayCreate(dptr, ESMF_DATA_REF, RC=STATUS)
     _VERIFY(STATUS)

     ! local fields built with 2d sized arrays and 2D grids
     srcFLD2D = ESMF_FieldCreate(Sgrid2D, newSrcARR,    &
                   horzRelloc = ESMF_CELL_CENTER,         &
                   datamap=dmap,                       &
                   haloWidth=0, &
                   name   = fname, rc=status )
     _VERIFY(STATUS)
     dstFLD2D = ESMF_FieldCreate(Dgrid2D, newDstARR,    &
                   horzRelloc = ESMF_CELL_CENTER,         &
                   datamap=dmap,                       &
                   haloWidth=0, &
                   name   = fname, rc=status )
     _VERIFY(STATUS)

     ! loop over field's vertical levels

     do k=1, kmax

     ! modify field's pointer by reference
     if (rank==3) then
       sptr(:,:) = sptr3d(:,:,k)
     else
       sptr(:,:) = sptr2d(:,:)
     end if

     call ESMF_FieldRegrid(srcFLD2D, dstFLD2D, rh, rc=status)
     _VERIFY(STATUS)

     !  Update contents (array pointer) of dstFLD
     call ESMF_FieldGetDataPointer (dstFLD2D, dptr, rc = status)
     _VERIFY(status)

     if (rank==3) then
       dptr3d(:,:,k) = dptr(:,:)
     else
       dptr2d(:,:) = dptr(:,:)
     end if

     end do ! k

     deallocate(sptr, dptr, stat=status)
     _VERIFY(STATUS)

   end do ! L

#else
   _RETURN(-1)
#endif
   end subroutine BundleRegrid1

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `BundleRegrid`
! regrids a source bundle (srcBUN) into a destination bundle (dstBUN)
! using hinterp. A bundle is thought of as being comprised of n 2D
! slices (nslices) distributed among the n PEs (ns_per_pe). The
! limits among each ns_per_pe region are given by n1 and n2 which
! are functions of mype (the local PE):
!
!```
!                                                     slice_pe
!           1 --- n1(pe=0)  -                     -->    0
!           2 ---            |                    -->    0
!              .             |_  ns_per_pe(pe=0)         .
!              .             |                           0
!              .             |                           0
!             --- n2(pe=0)  -                            0
!             --- n1(pe=1)                               1
!              .                                         .
!              .                                         .
!              .                                         .
!             --- n2(pe=1)                               1
!             --- n1(pe=2)                               2
!              .                                         .
!              .                                         .
!              .                                         .
!          ns ---                                  slice_pe(ns)
!              .                                         .
!              .                                         .
!              .                                         .
!     nslices --- n2(pe=n)                        -->   npe
!```
!
! Each slice is gathered, regridded (hinterp), and scattered on
! a PE determined by a slice-to-PE map (slice\_pe) to "load balance"
! the work of the serial hinterp function.
!
!#### History
!- 24Apr2006  Cruz  Initial code.
!
   subroutine BundleRegrid (srcBUN, dstBUN, rc)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_FieldBundle), intent(inout)      :: srcBUN !! source bundle
   type(ESMF_FieldBundle), intent(inout)      :: dstBUN !! destination bundle
   integer, optional, intent(out)             :: rc     !! return code

! local vars

   type(ESMF_VM)    :: vm
   type(ESMF_Grid)  :: srcGrid   ! grid associated with source bundle
   type(ESMF_Grid)  :: dstGrid   ! grid associated with destination bundle
   Logical          :: flip_poles
   Logical          :: flip_lons
   integer          :: numVars   ! number of fields in bundles
   integer          :: nslices   ! number of 2D slices in bundles
   integer          :: ns        ! counter for slices
   integer          :: mype      ! local PE
   integer          :: npe       ! number of PEs
   integer          :: nfirst    ! integer value of n1(mype)
   integer          :: nlast     ! integer value of n2(mype)
   integer          :: bufindex  ! local index of global buffer
   integer          :: ims_world ! x- global dimensions of src fields
   integer          :: jms_world ! y- global dimensions of src fields
   integer          :: imd_world ! x- global dimensions of dst fields
   integer          :: jmd_world ! y- global dimensions of dst fields
   integer          :: km_world  ! z- vertical dimension
   integer          :: ims       ! x- local dimensions of src fields
   integer          :: jms       ! y- local dimensions of src fields
   integer          :: imd       ! x- local dimensions of dst fields
   integer          :: jmd       ! y- local dimensions of dst fields
   integer          :: kfirst    ! lower bound of buffer's 3rd dim
   integer          :: status    ! return code status
   integer, allocatable, dimension(:    ) :: slice_pe
   integer, allocatable, dimension(:    ) :: ns_per_pe
   integer, allocatable, dimension(:    ) :: n1, n2
   real(kind=REAL32), allocatable, dimension(:,:,:) :: srcBUF  ! src buffer
   real(kind=REAL32), allocatable, dimension(:,:,:) :: dstBUF  ! dst buffer
   real(kind=REAL32),     pointer, dimension(:,:  ) :: llons, llats
   real(kind=REAL32), allocatable, dimension(:,:  ) :: glons, glats
   real(kind=REAL32), allocatable, dimension(:    ) :: srcLons, srcLats
   real(kind=REAL32), allocatable, dimension(:,:  ) :: srcWork
   real(kind=REAL32), allocatable, dimension(:,:  ) :: dstWork
   character(len=ESMF_MAXSTR), parameter  :: IAm = 'BundleRegrid'

! begin

   call ESMF_VMGetCurrent(vm)
   call ESMF_VMGet(vm, petCount=npe, localPet=mype, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

! get bundle grid information:
! global and local counts perdimension, km_world, nslices, numVars

   call Bundle_Prep_ (srcBUN, dstBUN)

! assign slices to PEs - create slice_pe array

   allocate(slice_pe(nslices), stat = status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(ns_per_pe(npe), stat = status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(n1(npe), n2(npe), stat = status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call assign_slices_ (nslices, mype, npe, slice_pe, nfirst, nlast)

! allocate buffers srcBUF, dstBUF

   call alloc_ (ims_world, jms_world, &
                imd_world, jmd_world, &
                nfirst, nlast)

! gather srcFLDs on all PEs into a srcBUF

   if(verbose .and. mype==MAPL_Root) print *,' - Gather...'
   call Do_Gathers_ (srcBUN, srcBUF)

! loop through slices and interpolate, use local addressing

   if(verbose .and. mype==MAPL_Root) then
      print *,' - Regrid ',nslices,' slices '
      print *,' SRC res ',ims_world,' x ', jms_world,&
           ' DST res ',imd_world,' x ', jmd_world
   end if

   bufindex = 0
   kfirst = lbound(srcBUF,3)
   do ns = 1, nslices
     bufindex = ns - n1(mype+1)
     bufindex = kfirst + bufindex
     if ( mype == slice_pe(ns) ) &
     call Do_Regrid_ (ns, srcBuf(:,:,bufindex), dstBuf(:,:,bufindex))
   end do

! scatter dstBUF to dstFLDs on all PEs

   if(verbose .and. mype==MAPL_Root) print *,' - Scatter...'
   call Do_Scatters_ (dstBUN, dstBUF)

   deallocate(slice_pe, ns_per_pe, n1, n2, stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call dealloc_

   _RETURN(ESMF_SUCCESS)

   contains

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `Bundle_Prep_` prepares for regridding.
!
!#### History
!- 24Apr2006  Cruz  Initial code.
!
   subroutine Bundle_Prep_ (srcBUN, dstBUN, only_vars)
!
! !ARGUMENTS:

   type(ESMF_FieldBundle), intent(inout)      :: srcBUN !! ALT: intent(in)
   type(ESMF_FieldBundle), intent(inout)      :: dstBUN
   character(len=*), optional, intent(in):: only_vars   !! comma separated, no spaces
!-------------------------------------------------------------------------

! locals

   type(ESMF_Array) :: srcArr
   type(ESMF_Field) :: srcFld, dstFld
   integer :: rank
   integer :: sCPD(3), dCPD(3)   ! src and dst counts per dimension (local)
   integer :: gsCPD(3), gdCPD(3) ! src and dst counts per dimension (global)
   integer :: srcLM              ! src vertical levels
   integer :: dstLM              ! dst
   integer :: dstNvars           ! numVars can be redefined if only_vars
                                 ! is specified
   integer :: L
   logical :: isPresent
   character(len=ESMF_MAXSTR)            :: name
   character(len=ESMF_MAXSTR), parameter :: IAm = 'Bundle_Prep_'


! get number of fields in bundles

   call ESMF_FieldBundleGet (srcBUN, FieldCount=NumVars, RC=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   nslices = 0
   dstNvars = 0
   do L = 1, NumVars

     call ESMF_FieldBundleGet(srcBUN, L, srcFLD, RC=STATUS)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     call ESMF_FieldBundleGet(dstBUN, L, dstFLD, RC=STATUS)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

     if ( present(ONLY_VARS) ) then
       if ( index(','//trim(ONLY_VARS)  //',', &
       ','//trim(name)//',') < 1 ) cycle
     endif
     dstNvars = dstNvars + 1

     call ESMFL_FieldGetDims(srcFLD, lm=srcLM)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     call ESMFL_FieldGetDims(dstFLD, lm=dstLM)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

! we can only do horizontal interpolation

     if(srcLM /= dstLM) then
       if(verbose .and. mype==MAPL_Root) then
         print *, 'Vertical interpolation is not implemented'
       end if
       _RETURN(ESMF_FAILURE)
     end if
     km_world = srcLM

     call ESMF_FieldGet(srcFLD, Array=srcARR, rc=status)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     call ESMF_ArrayGet(srcARR, RANK=rank, rc=status)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

     if ( rank == 3 ) then
       nslices = nslices + km_world
     else
       nslices = nslices + 1
     end if

   end do
   if ( present(ONLY_VARS) ) NumVars = dstNvars

! prepare buffers

   call ESMF_FieldBundleGet(srcBUN, 1, srcFLD, RC=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call ESMFL_FieldGetDims(srcFLD, gCPD=gsCPD)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   ims_world = gsCPD(1); jms_world = gsCPD(2)
   call ESMF_FieldBundleGet(dstBUN, 1, dstFLD, RC=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call ESMFL_FieldGetDims(dstFLD, gCPD=gdCPD)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   imd_world = gdCPD(1); jmd_world = gdCPD(2)

! Get local dimensions

   call ESMFL_FieldGetDims(srcFLD, lCPD=sCPD)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   ims = sCPD(1); jms = sCPD(2)
   call ESMFL_FieldGetDims(dstFLD, lCPD=dCPD)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   imd = dCPD(1); jmd = dCPD(2)

! get bundle grids

   call ESMF_FieldBundleGet (srcBUN, 1, srcFLD, RC=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call ESMF_FieldGet(srcFLD, grid=srcGRID, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call ESMF_FieldBundleGet (dstBUN, 1, dstFLD, RC=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call ESMF_FieldGet(dstFLD, grid=dstGRID, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

! get lons,lats used by hhinterp
   call ESMFL_GridCoordGet(   srcGrid, llats                     , &
                              Name     = "Latitude"              , &
                              Location = ESMF_STAGGERLOC_CENTER  , &
                              Units    = MAPL_UnitsRadians      , &
                              RC       = STATUS                    )
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(gLats(ims_world,jms_world), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call ArrayGather(llats, gLats, srcGrid,  RC=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(srcLats(jms_world), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   if(mype==MAPL_Root) then
      srcLats = gLats(1,:)
   end if
   call ESMF_VMBroadcast(vm, srcLats, jms_world, MAPL_Root, rc=status)

   call ESMFL_GridCoordGet(   srcGrid, llons                     , &
                              Name     = "Longitude"             , &
                              Location = ESMF_STAGGERLOC_CENTER  , &
                              Units    = MAPL_UnitsRadians      , &
                              RC       = STATUS                    )
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(gLons(ims_world,jms_world), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call ArrayGather(llons, gLons, srcGrid,  RC=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(srcLons(ims_world), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   if(mype==MAPL_Root) then
      srcLons = gLons(:,1)
   end if
   call ESMF_VMBroadcast(vm, srcLons, ims_world, MAPL_Root, rc=status)

   call ESMF_AttributeGet(dstGrid, 'VERBOSE', isPresent=isPresent, rc=status)
   if (isPresent) then
      call ESMF_AttributeGet(dstGrid, 'VERBOSE', verbose, rc=status)
      _VERIFY(STATUS)
   else
      verbose =.FALSE.
   end if

   call ESMF_AttributeGet(dstGrid, 'FLIP_LONS', isPresent=isPresent, rc=status)
   _VERIFY(STATUS)
   if (isPresent) then
      call ESMF_AttributeGet(dstGrid, 'FLIP_LONS', flip_lons, rc=status)
      _VERIFY(STATUS)
   else
      flip_lons = .FALSE.
   end if

   call ESMF_AttributeGet(dstGrid, 'FLIP_POLES', isPresent=isPresent, rc=status)
   _VERIFY(STATUS)
   if (isPresent) then
      call ESMF_AttributeGet(dstGrid, 'FLIP_POLES', flip_poles, rc=status)
      _VERIFY(STATUS)
   else
      flip_poles = .FALSE.
   end if

   if(mype==MAPL_Root.and.verbose) then
      if(flip_lons) print *, trim(Iam)//': We will flip lons'
      if(flip_poles) print *, trim(Iam)//': We will flip poles'
   end if
   _RETURN(ESMF_SUCCESS)

   end subroutine Bundle_Prep_

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `assign_slices_`
! determines number of bundle slices per PE and "load balanced"
! map of slices-to-pes (slice_pe).
!
!#### History
!- 24Apr2006  Cruz  Initial code.
!
   subroutine assign_slices_ (nslices, mype, npe, slice_pe, nfirst, nlast)
!
  implicit NONE

! !ARGUMENTS:

   integer, intent(in)    :: nslices       !! number of slices
   integer, intent(in)    :: mype          !! local PE
   integer, intent(in)    :: npe           !! number of PEs
   integer, intent(inout) :: slice_pe(:)   !! slice-to-pe map
   integer, intent(out)   :: nfirst
   integer, intent(out)   :: nlast
!
!-------------------------------------------------------------------------

! locals

   integer                            :: ns_rem, i, ipe, peidx
   integer, allocatable, dimension(:) :: displ

! start

   peidx = mype + 1

! calculate number of slices per PE

   ns_per_pe = nslices/npe
   ns_rem  = mod(nslices,npe)

! redistribute remaining slices

   allocate(displ(npe), stat = status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   do ipe = 1, npe
     if(ipe-1 < ns_rem) ns_per_pe(ipe) = ns_per_pe(ipe) + 1
     displ(ipe) = ns_per_pe(ipe) * ipe
     if(ipe-1 >= ns_rem) displ(ipe) = displ(ipe) + ns_rem
!?     call ESMF_VMBroadcast(vm, ns_per_pe, npe, ipe-1, rc=status)
!?     call ESMF_VMBroadcast(vm, displ, npe, ipe-1, rc=status)
   end do

! limits for slice_pe map

   n2 = displ
   n1 = n2 - ns_per_pe + 1
   nfirst = n1(peidx)
   nlast  = n2(peidx)

! determine slice_pe - this determines which PEs do the hinterp

   slice_pe = -999
   do ipe = 1, npe
     slice_pe (n1(ipe):n2(ipe)) = ipe-1
!?     call ESMF_VMBroadcast(vm, slice_pe, nslices, ipe-1, rc=status)
   end do

   deallocate(displ, stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   end subroutine assign_slices_

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `Do_Gathers_` gathers FLDs in a BUNdle on all PEs into a BUFfer.
!
! Note: local adressing is used.
!
!#### History
!- 28Apr2006  Cruz  Initial code.
!
   subroutine Do_Gathers_ (BUN, BUF)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_FieldBundle), intent(inout)         :: BUN
   real(kind=REAL32), intent(inout), dimension(:,:,:) :: BUF
!
!-------------------------------------------------------------------------

! locals

   type(ESMF_Field) :: FLD   ! ESMF field
   integer          :: n     ! number of vars in a bundle counter
   integer          :: L     ! vertical dim counter
   integer          :: rank  ! field rank
   integer          :: lmFLD ! field vertical dimension (1 or km_world)
   integer          :: bufi, kf
   integer          :: halowidth, hw
   real(kind=REAL32), dimension(:,:  ), pointer    :: ptr2d
   real(kind=REAL32), dimension(:,:,:), pointer    :: ptr3d
   character(len=ESMF_MAXSTR)            :: name
   character(len=ESMF_MAXSTR), parameter :: IAm = 'Do_Gathers_'

! start

   ns = 0
   bufi = 0
   kf = lbound(BUF,3)
   do n = 1, NumVars

     call ESMF_FieldBundleGet (BUN, n, FLD, RC=STATUS)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     ! get field name
     call ESMF_FieldGet(FLD, name=name, rc=status)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     call ESMFL_FieldGetDims (FLD, lm=lmFLD, ar=rank)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     ! check if field has halo, initialize to no halo
     hw = 0
     call ESMF_AttributeGet(FLD, "HALOWIDTH", halowidth, &
                            rc=status)
     if (status == ESMF_SUCCESS) hw = halowidth
     if (verbose .and. mype==MAPL_Root .and. n==1) print *, ' halowidth = ',hw

     if (rank==2) then
       call ESMF_FieldGet (FLD, localDE=0, farrayPtr=ptr2d, rc = status)
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     else
       call ESMF_FieldGet (FLD, localDE=0, farrayPtr=ptr3d, rc = status)
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     end if

     do L = 1, lmFLD

       ns = ns + 1

       ! gather from distributed pointer to global buffer BUF
       if (rank==2) then
          if(hw>0) then
             call ArrayGather(ptr2d       , srcWork, srcGRID, &
                  depe = slice_pe(ns), hw=hw, RC=STATUS)
          else
             call ArrayGather(ptr2d       , srcWork, srcGRID, &
                  depe = slice_pe(ns), RC=STATUS)
          end if
       else
          if(hw>0) then
             call ArrayGather(ptr3d(:,:,L), srcWork, srcGRID, &
                  depe = slice_pe(ns), hw=hw, RC=STATUS)
          else
             call ArrayGather(ptr3d(:,:,L), srcWork, srcGRID, &
                  depe = slice_pe(ns), RC=STATUS)
          end if
       end if
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

       call ESMF_VMBarrier(vm, rc=status)
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

       bufi = ns - n1(mype+1)
       bufi = kf + bufi

       if(mype==slice_pe(ns)) then
         BUF(:,:,bufi) = srcWork
       else
         srcWork = 0.0
       end if

       if (L==1) then
       if(verbose .and. mype==slice_pe(ns)) &
         write(*,'(1x,2(a,i4),2(a,f12.4))') &
         ' *** gathered buf index ',bufi,' on PE ',slice_pe(ns), &
         ' SRC min = ',minval(BUF(:,:,bufi)), &
             ' max = ',maxval(BUF(:,:,bufi))
       end if
     end do

   end do

   _RETURN(ESMF_SUCCESS)

   end subroutine Do_Gathers_

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `Do_Regrid_` calls `hinterp` on local PE.
!
!#### History
!- 28Apr2006  Cruz  Initial code.
!
   subroutine Do_Regrid_ (n, inBuf, outBuf)
!
  implicit NONE

! !ARGUMENTS:

   integer, intent(in)                    :: n      ! slice index
   real(kind=REAL32), intent(inout), dimension(:,:) :: inBuf  ! source buffer
   real(kind=REAL32), intent(inout), dimension(:,:) :: outBuf ! destination buffer
!
!-------------------------------------------------------------------------

! locals

   character(len=ESMF_MAXSTR), parameter :: IAm = 'Do_Regrid_'

 ! start

   call hhinterp ( inBuf , ims_world, jms_world   , &
                   outBuf, imd_world, jmd_world, 1, &
                   MAPL_Undef, srcLons, srcLats)

   if (flip_lons ) call FlipLons_ (outBuf)
   if (flip_poles) call FlipPoles_(outBuf)

   end subroutine Do_Regrid_

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `Do_Scatters_` scatters from BUffer onto FLDs in a BUNdle.
!
! Note: local adressing is used.
!
!#### History
!- 28Apr2006  Cruz  Initial code.
!
   subroutine Do_Scatters_ (BUN, BUF)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_FieldBundle), intent(inout)         :: BUN
   real(kind=REAL32), intent(inout), dimension(:,:,:) :: BUF
!
!-------------------------------------------------------------------------

! locals

   type(ESMF_Field) :: FLD
   integer          :: n     ! number of vars in a bundle counter
   integer          :: L     ! vertical dim counter
   integer          :: rank  ! field rank
   integer          :: lmFLD ! field vertical dimension (1 or km_world)
   integer          :: bufi, kf
   integer          :: halowidth, hw
   real(kind=REAL32), dimension(:,:  ), pointer    :: ptr2d
   real(kind=REAL32), dimension(:,:,:), pointer    :: ptr3d
   character(len=ESMF_MAXSTR)            :: name
   character(len=ESMF_MAXSTR), parameter :: IAm = 'Do_Scatters_'

! start

   ns = 0
   bufi = 0
   kf = lbound(BUF,3)
   do n = 1, NumVars

     call ESMF_FieldBundleGet (BUN, n, FLD, RC=STATUS)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     call ESMF_FieldGet(FLD, name=name, rc=status)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     call ESMFL_FieldGetDims (FLD, lm=lmFLD, ar=rank)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     ! check if field has halo, initialize to no halo
     hw = 0
     call ESMF_AttributeGet(FLD, "HALOWIDTH", halowidth, &
                                 rc=status)
     if (status == ESMF_SUCCESS) hw = halowidth
     if (verbose .and. mype==MAPL_Root .and. n==1) print *, ' halowidth = ',hw

     if (rank==2) then
       call ESMF_FieldGet (FLD, localDE=0, farrayPtr=ptr2d, rc = status)
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     else
       call ESMF_FieldGet (FLD, localDE=0, farrayPtr=ptr3d, rc = status)
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     end if

     do L = 1, lmFLD

       ns = ns + 1

       bufi = ns - n1(mype+1)
       bufi = kf + bufi

       if (L==1) then
       if(verbose .and. mype==slice_pe(ns)) &
         write(*,'(1x,2(a,i4),2(a,f12.4))') &
         ' *** scatter buf index ',bufi,' from PE ',mype, &
         ' DST min = ',minval(BUF(:,:,bufi)), &
             ' max = ',maxval(BUF(:,:,bufi))
       end if

       if(mype==slice_pe(ns)) then
         dstWork = BUF(:,:,bufi)
       else
         dstWork = 0.0
       end if

       ! scatter from buffer BUF to distributed pointer
       if (rank==2) then
          if(hw>0) then
             call ArrayScatter(ptr2d       , dstWork, dstGRID, &
                  depe=slice_pe(ns), hw=hw, RC=STATUS)
          else
             call ArrayScatter(ptr2d       , dstWork, dstGRID, &
                  depe=slice_pe(ns), RC=STATUS)
          endif
       else
          if(hw>0) then
             call ArrayScatter(ptr3d(:,:,L), dstWork, dstGRID, &
                  depe=slice_pe(ns), hw=hw, RC=STATUS)
          else
             call ArrayScatter(ptr3d(:,:,L), dstWork, dstGRID, &
                  depe=slice_pe(ns), RC=STATUS)
          endif
       end if
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

       ! CC: This barrier does not seem to be necessary
       !call ESMF_VMBarrier(vm, rc=status)
       !if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

     end do

   end do

   end subroutine  Do_Scatters_

!-------------------------------------------------------------------------
   subroutine alloc_ (ims,jms,imd,jmd,nfirst,nlast)
!-------------------------------------------------------------------------

   integer, intent(in) :: ims,jms,imd,jmd,nfirst,nlast
!
   integer :: status
   character(len=ESMF_MAXSTR), parameter :: IAm = 'alloc_'
!
   allocate(srcBUF(1:ims,1:jms,nfirst:nlast), STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   srcBUF = 0.0
   allocate(dstBUF(1:imd,1:jmd,nfirst:nlast), STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   dstBUF = 0.0
   allocate(srcWork(1:ims,1:jms), STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(dstWork(1:imd,1:jmd), STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   end subroutine alloc_

!-------------------------------------------------------------------------
   subroutine dealloc_
!-------------------------------------------------------------------------

   integer :: status
   character(len=ESMF_MAXSTR), parameter :: IAm = 'dealloc_'
!
   deallocate(srcBUF, STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   deallocate(dstBUF, STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   deallocate(srcWork, STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   deallocate(dstWork, STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   ! these are allocated in Bundle_Prep_
   deallocate(srcLons, srcLats, STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   deallocate(gLons, gLats, STAT=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   end subroutine dealloc_

!-------------------------------------------------------------------------
   subroutine FlipLons_(q)
!-------------------------------------------------------------------------
     implicit none
     real,dimension(:,:),intent(inout) :: q
     integer :: j
     do j=1,size(q,2)
        call FlipLonsi_(q(:,j))
     end do
   end subroutine FlipLons_

!-------------------------------------------------------------------------
   subroutine FlipLonsi_(q)
!-------------------------------------------------------------------------
     implicit none
     real,dimension(:),intent(inout) :: q
     integer :: im
     real,dimension(size(q,1)/2) :: d
     im=size(q,1)
     d(     1:im/2) = q(     1:im/2)
     q(     1:im/2) = q(im/2+1:im  )
     q(im/2+1:im  ) = d(     1:im/2)
   end subroutine FlipLonsi_

!-------------------------------------------------------------------------
   subroutine FlipPoles_(rbufr)
!-------------------------------------------------------------------------
     implicit none
     real(kind=REAL32),dimension(:,:),intent(inout) :: rbufr
     real(kind=REAL32),allocatable, dimension(:,:) :: sfcin
     integer :: im,jm

     im=size(rbufr,1)
     jm=size(rbufr,2)
     allocate(sfcin(im,jm))
     sfcin = rbufr
     rbufr(:,1:jm) = sfcin(:,jm:1:-1)
     deallocate(sfcin)
   end subroutine FlipPoles_

   end subroutine BundleRegrid

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `StateRegrid` regrids a state.
!
!#### History
!- 19Apr2006  Cruz  Initial code.
!
   subroutine StateRegrid (srcSTA, dstSTA, rc)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_State), intent(inout)  :: srcSTA
   type(ESMF_State), intent(inout)  :: dstSTA
   integer, optional, intent(out)   :: rc      !! return code
!
!-------------------------------------------------------------------------

! locals

   type(ESMF_VM)              :: vm
   type(ESMF_Field)           :: FLD
   type(ESMF_Grid)            :: GRID
   type(ESMF_FieldBundle)          :: srcBUN
   type(ESMF_FieldBundle)          :: dstBUN
   integer :: status, mype, npe, i, itemcount, srcFCount,dstFCount
   character (len=ESMF_MAXSTR), pointer  :: itemnamelist(:)
   character (len=ESMF_MAXSTR), allocatable  :: dstNames(:)
   character (len=ESMF_MAXSTR), allocatable  :: srcNames(:)
   type(ESMF_StateItem_Flag)  , pointer  :: itemtypelist(:)
   character(len=ESMF_MAXSTR)            :: name
   character(len=ESMF_MAXSTR), parameter :: IAm = 'StateRegrid'

! start
   call ESMF_VMGetCurrent(vm)
   call ESMF_VMGet(vm, localPet=mype, petCount=npe, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

! Get information from src state
!-------------------------------
! query the states for item number and types. If item type is other
! than a field then return(1)

   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call ESMF_StateGet(srcSTA, itemcount=itemcount, rc=status)

   if(itemcount==0) then
      call ESMFL_FailedRC(mype,Iam//': ESMF_state is empty')
   end if

   allocate(itemnamelist(itemcount), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(itemtypelist(itemcount), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   call ESMF_StateGet(srcSTA, itemnamelist=itemnamelist, &
                      itemtypelist=itemtypelist, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   srcFCount=0
   do i=1,itemcount
      if(itemtypelist(i)==ESMF_STATEITEM_FIELD) then
        srcFCount=srcFCount+1
        !call ESMFL_FailedRC(mype,trim(Iam)//': State item '//trim(itemNameList(i))//' is not a field.')
      end if
   end do
   allocate(srcNames(srcFCount))
   srcFCount=0
   do i=1,itemcount
      if(itemtypelist(i)==ESMF_STATEITEM_FIELD) then
        srcFCount=srcFCount+1
        srcNames(srcFCount)=trim(itemNameList(i))
      end if
   end do

! get grid from first field in state, and create bundle

   if(verbose .and. mype==MAPL_Root) print *, ' Get grid from ',trim(srcNames(1))
   call ESMF_StateGet(srcSTA, trim(srcNames(1)), FLD, &
                      rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' ESMF_StateGetField failed')
   call ESMF_StateGet(srcSTA, name=name, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' ESMF_StateGet failed')
   name = name//'to bundle'
   call ESMF_FieldGet (FLD, grid=GRID, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' ESMF_FieldGet failed')
   ! create bundle with input grid
   srcBUN = ESMF_FieldBundleCreate ( name=name, rc=status )
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' src ESMF_FieldBundleCreate failed' )
   call ESMF_FieldBundleSet(srcBUN, grid, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' src ESMF_FieldBundleSet failed' )

! populate source bundle

   call ESMFL_StateGetField(srcSTA, srcNames , srcBUN, RC=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' SRC ESMFL_StateGetField failed')

! Get information from dst state
!-------------------------------

   call ESMF_StateGet(dstSTA, itemcount=itemcount, rc=status)
   _VERIFY(status)
   deallocate(itemNameList)
   deallocate(itemTypeList)
   allocate(itemNameList(itemCount))
   allocate(itemTypeList(itemCount))
   call ESMF_StateGet(dstSTA, itemnamelist=itemnamelist, &
                      itemtypelist=itemtypelist, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   dstFCount=0
   do i=1,itemcount
      if(itemtypelist(i)==ESMF_STATEITEM_FIELD) then
        dstFCount=dstFCount+1
      end if
   end do
   allocate(dstNames(dstFCount))
   dstFCount=0
   do i=1,itemcount
      if(itemtypelist(i)==ESMF_STATEITEM_FIELD) then
        dstFCount=dstFCount+1
        dstNames(dstFCount)=trim(itemNameList(i))
      end if
   end do
   _ASSERT(dstFCount==srcFCount,"source and destination bundle sizes do not match")
! get grid from first field in state, and create bundle

   if(verbose .and. mype==MAPL_Root) print *, ' Get grid from ',trim(dstNames(1))
   call ESMF_StateGet(dstSTA, trim(dstNames(1)), FLD, &
                      rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' ESMF_StateGetField failed')
   call ESMF_StateGet(dstSTA, name=name, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' ESMF_StateGet failed')
   name = name//'to bundle'
   call ESMF_FieldGet (FLD, grid=GRID, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' ESMF_FieldGet failed')
   ! create bundle with input grid
   dstBUN = ESMF_FieldBundleCreate ( name=name, rc=status )
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' ESMF_FieldBundleCreate failed' )
   call ESMF_FieldBundleSet(dstBUN, grid, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' src ESMF_FieldBundleSet failed' )

! populate destination bundle

   call ESMFL_StateGetField(dstSTA, dstNames, dstBUN, RC=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam//' DST ESMFL_StateGetField failed')

! now call BundleRegrid

   call BundleRegrid (srcBUN, dstBUN, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

! clean up

   deallocate(itemnamelist, stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   deallocate(itemtypelist, stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   _RETURN(ESMF_SUCCESS)

   end subroutine StateRegrid

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `ESMFL_FieldGetDims`
! returns some grid information associated from an ESMF field.
!
!#### History
!- 24Apr2006  Cruz  Initial code.
!
   subroutine ESMFL_FieldGetDims(FLD, gCPD, lCPD, lm, ar)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_Field), intent(inout)         :: FLD !ALT: intent(in)
   integer, optional, intent(out)       :: gCPD(3)
   integer, optional, intent(out)       :: lCPD(3)
   integer, optional, intent(out)       :: lm
   integer, optional, intent(out)       :: ar
!
!-------------------------------------------------------------------------

! locals

   type(ESMF_VM)    :: vm
   type(ESMF_Grid)  :: GRID
   type(ESMF_Array) :: ARR
   integer          :: globalCPD(3)
   integer          :: localCPD(3)
   integer          :: status, thisrank, mype
   character(len=ESMF_MAXSTR), parameter :: IAm = 'ESMFL_FieldGetDims'

! start

   call ESMF_VMGetCurrent(vm)
   call ESMF_VMGet(vm, localPet=mype, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   call ESMF_FieldGet(FLD, grid=GRID, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call MAPL_GridGet(GRID, &
        globalCellCountPerDim=globalCPD, &
        localCellCountPerDim=localCPD, &
        rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   if (present(gCPD)) gCPD = globalCPD
   if (present(lCPD)) lCPD = localCPD

   call ESMF_FieldGet(FLD, Array=ARR, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   call ESMF_ArrayGet(ARR, RANK=thisrank, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   if (present(lm)) then
     if ( thisrank >= 3 ) then
       lm = globalCPD(3)
     else
       lm = 1
     end if
   end if

   if (present(ar)) ar = thisrank

   end subroutine ESMFL_FieldGetDims

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! The subroutine `BundleDiff` determines the diff of two bundles.
!
!#### History
!- 24Apr2006  Cruz  Initial code.
!
   subroutine BundleDiff (srcBUN, dstBUN, rc)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_FieldBundle),           intent(inout) :: srcBUN
   type(ESMF_FieldBundle), optional, intent(inout) :: dstBUN
   integer,                optional, intent(out)   :: rc     !! return code
!
!-------------------------------------------------------------------------

! local vars

   type(ESMF_VM)           :: vm
   type(ESMF_Field)        :: dstFLD
   type(ESMF_Field)        :: srcFLD
   type(ESMF_Grid)         :: grid3D
   real(kind=REAL32), dimension(:,:    ), pointer   :: sptr2d, dptr2d
   real(kind=REAL32), dimension(:,:,:  ), pointer   :: sptr3d, dptr3d
   real(kind=REAL32), dimension(:,:,:,:), pointer   :: sptr4d, dptr4d
   real(kind=REAL32), allocatable, dimension(:,:) :: srcBuf, dstBuf, diff
   integer :: rank, k, mype, kmax, status, L, numVars, u, umax ! umax is the 4th dimension
   integer :: globalCPD(3)
   character(len=ESMF_MAXSTR)            :: srcName
   character(len=ESMF_MAXSTR), parameter :: fname = 'aleph'
   character(len=ESMF_MAXSTR), parameter :: IAm = 'BundleDiff'

! begin

   call ESMF_VMGetCurrent(vm)
   call ESMF_VMGet(vm, localPet=mype, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   ! get number of fields in bundles
   ! number in srcBUN should be the same as in dstBUN (we can change later)
   call ESMF_FieldBundleGet      (srcBUN, FieldCount=NumVars,  RC=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   ! use only srcBUN; bundles are conformant
   call ESMF_FieldBundleGet     (srcBUN, 1, srcFLD, RC=STATUS)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   ! get grid
   call ESMF_FieldGet(srcFLD, grid=grid3D, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   ! get global counts per dim
   call ESMFL_FieldGetDims(srcFLD, gCPD=globalCPD)
   allocate(srcbuf(globalCPD(1),globalCPD(2)), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(dstbuf(globalCPD(1),globalCPD(2)), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(diff(globalCPD(1),globalCPD(2)), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   ! loop over fields in bundle

   if(verbose .and. mype==MAPL_Root) print *,' Loop through bundles'

   do L=1,NumVars

     call ESMF_FieldBundleGet     (srcBUN, L, srcFLD, RC=STATUS)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     call ESMF_FieldGet           (srcFLD, NAME=srcName, RC=STATUS)
     if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     if(present(dstBUN)) then
        call ESMF_FieldBundleGet     (dstBUN, L, dstFLD, RC=STATUS)
        if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
     endif
     ! get rank from any field
     call ESMFL_FieldGetDims(srcFLD, ar=rank)

     ! allocate f90 pointer to hold data
     ! note halo width = 0
     umax = 1 ! by default, no 4th dimension
     if(rank==4) then
       call ESMF_FieldGet (srcFLD, localDE=0, farrayPtr=sptr4d, rc = status)
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
       if (present(dstBUN)) then
         call ESMF_FieldGet (dstFLD, localDE=0, farrayPtr=dptr4d, rc = status)
         if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
       endif
       umax = size(sptr4d,4)
     else if(rank==3) then
       call ESMF_FieldGet (srcFLD, localDE=0, farrayPtr=sptr3d, rc = status)
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
       if (present(dstBUN)) then
         call ESMF_FieldGet (dstFLD, localDE=0, farrayPtr=dptr3d, rc = status)
         if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
       endif
     else if(rank==2) then
       call ESMF_FieldGet (srcFLD, localDE=0, farrayPtr=sptr2d, rc = status)
       if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
       if (present(dstBUN)) then
          call ESMF_FieldGet (dstFLD, localDE=0, farrayPtr=dptr2d, rc = status)
          if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
       endif
     else
       if (MAPL_AM_I_Root()) print*, "WARNING: rank=1 field not implmented yet"
       cycle
     end if

    ! loop over field's vertical levels

     call ESMFL_FieldGetDims(srcFLD, lm=kmax)

     do u = 1, umax ! loop 4th ungrid dimension

        if ( umax > 1) then
          srcName = trim(srcName)//'_'//trim(ESMF_UtilStringInt2String (u))
        endif

        do k=1, kmax

          if (rank==4) then
            call ArrayGather(sptr4d(:,:,k,u), srcBuf, grid3D, RC=STATUS)
            if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
            if (present(dstBUN)) then
               call ArrayGather(dptr4d(:,:,k,u), dstBuf, grid3D, RC=STATUS)
               if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
            endif
          elseif (rank==3) then
            call ArrayGather(sptr3d(:,:,k), srcBuf, grid3D, RC=STATUS)
            if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
            if (present(dstBUN)) then
               call ArrayGather(dptr3d(:,:,k), dstBuf, grid3D, RC=STATUS)
               if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
            endif
          else
            call ArrayGather(sptr2d, srcBuf, grid3D, RC=STATUS)
            if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
            if (present(dstBUN)) then
              call ArrayGather(dptr2d, dstBuf, grid3D, RC=STATUS)
              if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
            endif
          end if
          call ESMF_VMBarrier(vm, rc=status)
          if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

          if(verbose .and. mype==MAPL_Root) then
            ! stats on level-by-level basis
            if(present(dstBUN)) then
               call stats_ ( 6, globalCPD(1), globalCPD(2), k, &
                         srcBuf,                    &
                         trim(srcName), 'lev', MAPL_undef , &
                         trim(srcName)//' diff statistics', 1, &
                         dstBuf )
            else
               call stats_ ( 6, globalCPD(1), globalCPD(2), k, &
                         srcBuf,                    &
                         trim(srcName), 'lev', MAPL_undef , &
                         trim(srcName)//' statistics', 1 &
                         )

            endif

          end if

        end do ! k
     end do ! u, ungrid fourth
   end do ! L

   if(present(dstBUN)) deallocate(dstBuf)

   deallocate(srcbuf, diff, stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   _RETURN(ESMF_SUCCESS)

CONTAINS

         logical function spv(aspv, amiss)
            real, intent(in) :: aspv
            real, intent(in) :: amiss

            real, parameter :: rfrcval = 1.0e-6

            spv = abs((aspv-amiss)/amiss).le.rfrcval
         end function spv

!>
! Print statistics of one 3-d variable. This is from the PSAS library.
! with some simplifications.
!
        subroutine stats_ (lu,mx,my,k,a1,&
                           atype,htype,amiss,header,inc, a2)

        implicit none

        integer lu              ! Output unit
        integer mx,my           ! Array sizes
        integer k               ! level
        real a1(mx,my)          ! The array1
        character(*), intent(in) :: atype       ! Type of the variable
        character(*), intent(in) :: htype       ! Typf of the levels
        real amiss              ! missing value flag of a
        character(len=*) header    ! A header message
        integer inc             ! order of the listing
        real,optional :: a2(mx,my)        ! The array2
!
        integer i,j
        integer imx,imn,jmx,jmn
        integer cnt
        real a(mx,my)        ! diff
        real amx,amn
        real avg,dev,d,rms,sumsq,relerr
        logical first

!       ..A practical value for the magnitude of the fraction of a real
!       number.

        character(len=255) dash

!       ..function

!       compute diff
        if (present(a2)) then
           a = a1 - a2
        else
           a = a1
        endif
        ! only print header when k==1

        if (k==1) then

           do i = 1, 255
              dash(i:i) = '-'
           end do
           i = len(trim(header))

           write(lu,'(//a)') trim(header)
           write(lu,'(a/)')  dash(1:i)
           if(present(a2)) then
              write(lu,'(a,1x,a,5x,a,6x,a,9x,a,9x,a,15x,a)') 'lvl',&
               'count','mean','rmse','rele','maxi','mini'
           else
              write(lu,'(a,1x,a,5x,a,6x,a,9x,a,15x,a)') 'lvl',&
               'count','mean','stdev','maxi','mini'

           endif ! present
        end if ! k==1

! compute some statistics on diff

          cnt=0
          avg=0.
          rms=0.
          do j=1,my
            do i=1,mx
              if(.not.spv(a(i,j),amiss)) then
                cnt=cnt+1
                avg=avg+a(i,j)
                if(present(a2)) then
                   sumsq = sumsq + a(i,j)*a(i,j)
                   relerr = relerr + (a1(i,j) - a2(i,j)) !/a1(i,j)
                endif
              endif
            end do
          end do

          avg=avg/max(1,cnt)
          if(present(a2)) then
             sumsq=sumsq/max(1,cnt)
             rms = sqrt(sumsq)
             relerr=relerr/max(1,cnt)
          endif

          dev=0.
          do j=1,my
            do i=1,mx
              if(.not.spv(a(i,j),amiss)) then
                d=a(i,j)-avg
                dev=dev+d*d
              endif
            end do
          end do
          dev=sqrt(dev/max(1,cnt-1))

          amx=a(1,1)
          amn=a(1,1)
          first=.true.
          do j=1,my
            do i=1,mx
              if(.not.spv(a(i,j),amiss)) then
                if(first) then
                  imx=i
                  imn=i
                  jmx=j
                  jmn=j
                  amx=a(imx,jmx)
                  amn=a(imn,jmn)
                  first=.false.
                else
                  if(a(i,j).gt.amx) then
                    amx=a(i,j)
                    imx=i
                    jmx=j
                  endif
                  if(a(i,j).lt.amn) then
                    amn=a(i,j)
                    imn=i
                    jmn=j
                  endif
                endif
              endif
            end do
          end do

          if(present(a2)) then
             write(lu,'(i3,1p,i7,1p,3e10.3e1,0p,'// &
                '2(1p,e10.3e1,0p,a,i3,a,i3,a))') &
                k,cnt,avg,rms,relerr, &
                amx,'(',imx,',',jmx,')', &
                amn,'(',imn,',',jmn,')'
           else
             write(lu,'(i3,1p,i7,1p,2e10.3e1,0p,'// &
                '2(1p,e10.3e1,0p,a,i3,a,i3,a))') &
                k,cnt,avg,dev, &
                amx,'(',imx,',',jmx,')', &
                amn,'(',imn,',',jmn,')'

           endif

           _UNUSED_DUMMY(ATYPE)
           _UNUSED_DUMMY(HTYPE)
           _UNUSED_DUMMY(INC)
      end  subroutine stats_

   end subroutine BundleDiff

!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! Determine the diff of two state.
!
!#### History
!- 19Apr2006  Cruz  Initial code.
!
   subroutine StateDiff (srcSTA, dstSTA, rc)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_State), intent(inout)  :: srcSTA
   type(ESMF_State),optional, intent(inout)  :: dstSTA
   integer, optional, intent(out)   :: rc     ! return code
!
!-------------------------------------------------------------------------

! locals

   type(ESMF_VM)              :: vm
   type(ESMF_FieldBundle)          :: srcBUN
   type(ESMF_FieldBundle)          :: dstBUN
   integer :: status, mype, i, itemcount
   character (len=ESMF_MAXSTR), pointer  :: itemnamelist(:)
   type(ESMF_StateItem_Flag)  , pointer  :: itemtypelist(:)
   character(len=ESMF_MAXSTR), parameter :: IAm = 'StateDiff'

! start

   call ESMF_VMGetCurrent(vm)
   call ESMF_VMGet(vm, localPet=mype, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

! Get information from state
!---------------------------
! query the states for item number and types. If item type is other
! than a field then return(1)

   call ESMF_StateGet(srcSTA, itemcount=itemcount, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   if(itemcount==0) then
      !call ESMFL_FailedRC(mype,Iam//': ESMF_state is empty')
      if (MAPL_AM_I_Root()) print*, "WARNING: "//Iam//': ESMF_state is empty'
      _RETURN(_SUCCESS)
   end if

   allocate(itemnamelist(itemcount), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
   allocate(itemtypelist(itemcount), stat=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   call ESMF_StateGet(srcSTA, itemnamelist=itemnamelist, &
                      itemtypelist=itemtypelist, rc=status)
   if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

   do i=1,itemcount
      if(itemtypelist(i)/=ESMF_STATEITEM_FIELD) then
        !call ESMFL_FailedRC(mype,Iam//': State item is not a field.')
        if (MAPL_AM_I_Root()) print*, "WARNING: "//Iam//':  State item is not a field'
      end if
   end do

! populate source bundle with fields from source state
! ----------------------------------------------------
   call State2Bundle (srcSTA, srcBun)

   if(present(dstSTA)) then

      call ESMF_StateGet(dstSTA, itemnamelist=itemnamelist, &
                      itemtypelist=itemtypelist, rc=status)
      if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

      do i=1,itemcount
         if(itemtypelist(i)/=ESMF_STATEITEM_FIELD) then
          ! call ESMFL_FailedRC(mype,Iam//': State item is not a field.')
            if (MAPL_AM_I_Root()) print*, "WARNING: "//Iam//':  State item is not a field'
         end if
      end do
      call State2Bundle (dstSTA, dstBun)

   ! now call BundleDiff
      call BundleDiff (srcBUN, dstBUN, rc=status)
      if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)
    else

      call BundleDiff (srcBUN, rc=status)
      if (status /= ESMF_SUCCESS) call ESMFL_FailedRC(mype,Iam)

    endif

   _RETURN(ESMF_SUCCESS)

   end subroutine StateDiff

   subroutine StateStatistics(srcSTA, rc)
     implicit NONE
     type(ESMF_State), intent(inout)  :: srcSTA
     integer, optional, intent(out)   :: rc     ! return code
     integer :: status

     call StateDiff(srcSTA, rc=status)
     _VERIFY(status)
     _RETURN(ESMF_SUCCESS)

   end subroutine StateStatistics

   subroutine BundleStatistics(srcBUN, rc)
     implicit NONE
     type(ESMF_FieldBundle), intent(inout) :: srcBUN
     integer, optional, intent(out)   :: rc     ! return code

     integer :: status

     call Bundlediff(srcBUN, rc=status)
     _VERIFY(status)
     _RETURN(ESMF_SUCCESS)

   end subroutine BundleStatistics
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!  NASA/GSFC, Global Modeling and Assimilation Office, Code 610.3, GMAO  !
!-------------------------------------------------------------------------
!>
! `ESMFL_GridDistBlockSet` ...
!
!#### History
!-16Jun2006  Cruz  Initial code.
!
   subroutine ESMFL_GridDistBlockSet (Egrid, ist, jst, il, jl, &
                                      rlons, rlats, rc)
!
  implicit NONE

! !ARGUMENTS:

   type(ESMF_Grid), intent(inout)            :: Egrid
   integer,         intent(in), dimension(:) :: ist, jst, il, jl
   real(kind=REAL64), optional, dimension(:) :: rlats
   real(kind=REAL64), optional, dimension(:) :: rlons
   integer,           optional, intent(out)  :: rc     !! return code

   _UNUSED_DUMMY(Egrid)
   _UNUSED_DUMMY(ist)
   _UNUSED_DUMMY(jst)
   _UNUSED_DUMMY(il)
   _UNUSED_DUMMY(jl)
   _UNUSED_DUMMY(rlats)
   _UNUSED_DUMMY(rlons)

#if 0
!
!-------------------------------------------------------------------------

   integer, allocatable, dimension(:) :: iAI
   integer :: ncnt,i,j,k,count,n,status,nDE,mype
   integer, allocatable, dimension(:,:) :: myindices
   type(ESMF_Logical) :: HasBlockSet
   type(ESMF_DELayout):: layout
   type(ESMF_VM):: vm
   type(ESMF_AxisIndex), allocatable, dimension(:,:) :: AI
   character(len=ESMF_MAXSTR), parameter :: IAm = 'ESMFL_GridDistBlockSet'

! start

   call ESMF_VMGetCurrent(vm)
   call ESMF_VMGet(vm, petCount=nDE, localPet=mype, rc=status)
   if (status /= ESMF_SUCCESS) &
   call ESMFL_FailedRC(mype,Iam//': ESMF_VMGet failed')

   allocate(myindices(il(mype+1)*jl(mype+1),2), stat=status)
   if (status /= ESMF_SUCCESS) &
   call ESMFL_FailedRC(mype,Iam//': allocate failed')

   ncnt=0
   do i=jst(mype+1),jst(mype+1)+jl(mype+1)-1
     do j=ist(mype+1),ist(mype+1)+il(mype+1)-1
       ncnt=ncnt+1
       myindices(ncnt,1) = i
       myindices(ncnt,2) = j
     end do
   end do
   !print *,mype,' ncnt = ',ncnt,shape(myindices),il(mype+1)*jl(mype+1)

   layout = ESMF_DELayoutCreate(vm, deCountList=(/1, 8/), rc=status)
   if (status /= ESMF_SUCCESS) &
   call ESMFL_FailedRC(mype,Iam//': ESMF_DELayoutCreate failed')

   call ESMF_GridDistribute(Egrid,               &
         deLayout=layout,                        &
         mycount=ncnt,                           &
         myindices=myindices,                    &
         rc=status)
   if (status /= ESMF_SUCCESS) &
   call ESMFL_FailedRC(mype,Iam//': ESMF_GridDistribute failed')

! Encode the AI information. Let's use the grid attributes to do that.
! First, serialize the AI object: create an integer 1D array to hold the
! AI information: 3D grids will be larger

   count = nDE*6
   allocate(iAI(count))
   iAI = 0
   k = 1
   do n = 1, nDE
     iAI(k) = ist(n); k=k+1
     iAI(k) = jst(n); k=k+1
     iAI(k) = il(k) ; k=k+1
     iAI(k) = jl(k) ; k=k+1
     iAI(k) = 0     ; k=k+1
     iAI(k) = 0     ; k=k+1
   end do

! Attach attributes to the grid:

   call ESMF_GridSetAttribute(Egrid, 'GlobalAxisIndex', count, iAI, rc=status)

   HasBlockSet = ESMF_TRUE
   call ESMF_GridSetAttribute(Egrid, 'HasBlockSet', HasBlockSet, rc)

! if rlons and rlats were present attach them to grid

   if(present(rlats)) &
   call ESMF_GridSetAttribute(Egrid, 'GSI rlats', size(rlats), rlats, rc)
   if(present(rlons)) &
   call ESMF_GridSetAttribute(Egrid, 'GSI rlons', size(rlons), rlons, rc)
!
   deallocate(iAI, myindices)
#endif

   if (present(rc)) then
      rc = 0
   end if
   end subroutine ESMFL_GridDistBlockSet

! THE FOLLOWING ROUTINE MAY BE TEMPORARY

!-------------------------------------------------------------------------
   subroutine State2Bundle (STA, BUN, usrfldlist, rc)
!-------------------------------------------------------------------------
! extract fields from state and add to a bundle
   type(ESMF_FieldBundle), intent(inout) :: BUN
   type(ESMF_State) , intent(inout) :: STA
   character(len=*),intent(in),optional :: usrfldlist
   integer, intent(out), optional   :: rc

! locals

   integer          :: status
   integer          :: mype
   integer          :: i, itemcount
   type(ESMF_VM)    :: VM
   type(ESMF_Field) :: FLD
   type(ESMF_Grid)  :: GRID
   character (len=ESMF_MAXSTR)           :: fldname
   character (len=ESMF_MAXSTR), pointer  :: namelist(:)
   type(ESMF_StateItem_Flag)  , pointer  :: typelist(:)
   character(len=ESMF_MAXSTR)            :: name
   character(len=ESMF_MAXSTR), parameter :: IAm='State2Bundle'
   logical                               :: field_found_already

! start

   call ESMF_VMGetCurrent(VM)
   call ESMF_VMGet(vm, localPet=mype, rc=status)
   _VERIFY(STATUS)
   call ESMF_StateGet(STA, itemcount=itemcount, rc=status)
   _VERIFY(STATUS)
   if(itemcount==0) then
      if (verbose) print *,' state is empty'
      _RETURN(ESMF_FAILURE)
   end if

   allocate(namelist(itemcount), stat=status)
   _VERIFY(STATUS)
   allocate(typelist(itemcount), stat=status)
   _VERIFY(STATUS)

   call ESMF_StateGet(STA, itemnamelist=namelist, &
                      itemtypelist=typelist, rc=status)
   _VERIFY(STATUS)

   field_found_already = .false.
   do i=1,itemcount
      if(typelist(i)/=ESMF_STATEITEM_FIELD) cycle

      call ESMF_StateGet(STA, trim(namelist(i)), FLD, &
           rc=status)
      _VERIFY(STATUS)
      if(verbose .and. MAPL_AM_I_ROOT()) print *, ' Got ',trim(namelist(i)),' from field'
      ! get grid from first field in state, and create bundle
      if (.not. field_found_already) then
         field_found_already = .true.

         call ESMF_StateGet(STA, name=name, rc=status)
         _VERIFY(STATUS)
         name = name//'to bundle'
         call ESMF_FieldGet (FLD, grid=GRID, rc=status)
         _VERIFY(STATUS)
         ! create bundle with input grid
         BUN = ESMF_FieldBundleCreate ( name=name, rc=status )
         _VERIFY(STATUS)
         call ESMF_FieldBundleSet(BUN, grid=grid, rc=status)
         _VERIFY(STATUS)

      end if

     ! populate bundle

     if (present(usrfldlist) ) then
         call ESMF_FieldGet(FLD, NAME=fldname, rc=status )
         _VERIFY(STATUS)
         if(.not. check_list_(fldname,usrfldlist)) cycle
     endif
     call MAPL_FieldBundleAdd (BUN, FLD, rc=status)
     _VERIFY(STATUS)

   end do

   deallocate(typelist)
   deallocate(namelist)

   if (.not. field_found_already) then
      _RETURN(ESMF_FAILURE)
   end if

   _RETURN(ESMF_SUCCESS)

   end subroutine State2Bundle

!-------------------------------------------------------------------
   subroutine Bundle2State (BUN, STA, rc)
!-------------------------------------------------------------------
! extract fields from bundle and add to a state
   type(ESMF_FieldBundle), intent(inout) :: BUN !ALT: intent(in)
   type(ESMF_State) , intent(inout):: STA
  integer, intent(out), optional   :: rc
! locals
   integer          :: L, status,NumVars
   type(ESMF_Field) :: FLD
   type(ESMF_Field) :: stateFIELD
   character(len=ESMF_MAXSTR)              :: name
   character(len=ESMF_MAXSTR), parameter   :: IAm='Bundle2State'

   real, pointer :: src_ptr1d(:)
   real, pointer :: dst_ptr1d(:)
   real, pointer :: src_ptr2d(:,:)
   real, pointer :: dst_ptr2d(:,:)
   real, pointer :: src_ptr3d(:,:,:)
   real, pointer :: dst_ptr3d(:,:,:)
   real(kind=ESMF_KIND_R8), pointer :: src_pr81d(:)
   real(kind=ESMF_KIND_R8), pointer :: dst_pr81d(:)
   real(kind=ESMF_KIND_R8), pointer :: src_pr82d(:,:)
   real(kind=ESMF_KIND_R8), pointer :: dst_pr82d(:,:)
   real(kind=ESMF_KIND_R8), pointer :: src_pr83d(:,:,:)
   real(kind=ESMF_KIND_R8), pointer :: dst_pr83d(:,:,:)
   type(ESMF_TypeKind_Flag) :: src_tk, dst_tk
   integer                  :: src_fieldRank, dst_fieldRank
   logical                  :: NotInState
   character(len=ESMF_MAXSTR) :: NameInBundle
   type(ESMF_StateItem_Flag) :: itemType


!
   call ESMF_FieldBundleGet (BUN, name=name, FieldCount=NumVars, RC=STATUS)
   if(verbose .and. MAPL_AM_I_ROOT()) print *,' Unwrapping ',trim(name)
   do L=1,NumVars
      call ESMF_FieldBundleGet (BUN, L, FLD, RC=STATUS)
      _VERIFY(STATUS)

      call ESMF_FieldGet(FLD, Name=NameInBundle, rc=STATUS)
      _VERIFY(STATUS)

      call ESMF_StateGet (STA, nameInBundle, stateField,  RC=STATUS)
      notInState=(status/=ESMF_SUCCESS)
      if (.not.notINState) then
         call ESMF_StateGet (STA, NameInBundle, itemType=itemType,  RC=STATUS)
         _VERIFY(STATUS)
      end if

      if (itemType /= ESMF_STATEITEM_FIELD) then
         cycle
      end if
      if (NotInState) then
         call MAPL_StateAdd (STA, FLD, rc=status)
         _VERIFY(STATUS)
      else
         ! copy content to StateField
         call ESMF_FieldGet(stateFIELD, &
              dimCount=dst_fieldRank, &
              typeKind=dst_tk, rc=status)
         _VERIFY(STATUS)
         call ESMF_FieldGet(FLD, dimCount=src_fieldRank, &
              typeKind=src_tk, rc=status)
         _VERIFY(STATUS)
         ! make sure that type kind rank agrees between
         ! fld and stateField
         _ASSERT(dst_fieldRank == src_fieldRank, 'Rank mismatch between fld and stateField')
         _ASSERT(dst_tk == src_tk, 'Type Kind mismatch between fld and stateField')
         ! for appropriate TKR:
         if (dst_tk == ESMF_TypeKind_R4) then
            select case (dst_fieldRank)
            case (1)
               call ESMF_FieldGet(fld, farrayPtr=src_ptr1d, rc=status)
               _VERIFY(STATUS)
               call ESMF_FieldGet(stateField, farrayPtr=dst_ptr1d, rc=status)
               _VERIFY(STATUS)
               dst_ptr1d = src_ptr1d
            case (2)
               call ESMF_FieldGet(fld, farrayPtr=src_ptr2d, rc=status)
               _VERIFY(STATUS)
               call ESMF_FieldGet(stateField, farrayPtr=dst_ptr2d, rc=status)
               _VERIFY(STATUS)
               dst_ptr2d = src_ptr2d
            case (3)
               call ESMF_FieldGet(fld, farrayPtr=src_ptr3d, rc=status)
               _VERIFY(STATUS)
               call ESMF_FieldGet(stateField, farrayPtr=dst_ptr3d, rc=status)
               _VERIFY(STATUS)
               dst_ptr3d = src_ptr3d
            case default
            end select
         else if (dst_tk == ESMF_TypeKind_R8) then
            select case (dst_fieldRank)
            case (1)
              call ESMF_FieldGet(fld, farrayPtr=src_pr81d, rc=status)
              _VERIFY(STATUS)
              call ESMF_FieldGet(stateField, farrayPtr=dst_pr81d, rc=status)
              _VERIFY(STATUS)
              dst_pr81d = src_pr81d
           case (2)
              call ESMF_FieldGet(fld, farrayPtr=src_pr82d, rc=status)
              _VERIFY(STATUS)
              call ESMF_FieldGet(stateField, farrayPtr=dst_pr82d, rc=status)
              _VERIFY(STATUS)
              dst_pr82d = src_pr82d
           case (3)
              call ESMF_FieldGet(fld, farrayPtr=src_pr83d, rc=status)
              _VERIFY(STATUS)
              call ESMF_FieldGet(stateField, farrayPtr=dst_pr83d, rc=status)
              _VERIFY(STATUS)
              dst_pr83d = src_pr83d
            case default
               _FAIL( 'unsupported rank (>= 4)')
            end select
         end if
      end if

   end do
   _RETURN(ESMF_SUCCESS)

   end subroutine Bundle2State

!-------------------------------------------------------------------
   subroutine Bundles2Bundle (BUN1, BUN2, mergedBUN, rc)
!-------------------------------------------------------------------
   type(ESMF_FieldBundle), intent(inout) :: BUN1
   type(ESMF_FieldBundle), intent(inout) :: BUN2
   type(ESMF_FieldBundle), intent(inout) :: mergedBUN
   integer, intent(out), optional :: rc
!
   type(ESMF_VM)           :: vm
   type(ESMF_Field)        :: FLD
   integer :: status, mype, L, numVars1, numVars2
   character(len=ESMF_MAXSTR)            :: name1, name2, name3
   character(len=ESMF_MAXSTR), parameter :: IAm = 'Bundles2Bundle'

! begin

   call ESMF_VMGetCurrent(vm)
   call ESMF_VMGet(vm, localPet=mype, rc=status)
   _VERIFY(STATUS)

   ! get number of fields in bundles
   call ESMF_FieldBundleGet (BUN1, name=name1, FieldCount=NumVars1, RC=STATUS)
   _VERIFY(STATUS)
   call ESMF_FieldBundleGet (BUN2, name=name2, FieldCount=NumVars2, RC=STATUS)
   _VERIFY(STATUS)
   ! mergedBUN is an empty bundle created by the parent routine
   call ESMF_FieldBundleGet (mergedBUN, name=name3,  RC=STATUS)
   _VERIFY(STATUS)

   ! loop over fields in bundle and add to mergedBUN

   if(verbose .and. MAPL_AM_I_ROOT()) &
   print *,' Add ',NumVars1,' fields from ',trim(name1), &
   ' into ', trim(name3)
   do L=1,NumVars1
     call ESMF_FieldBundleGet (BUN1, L, FLD, RC=STATUS)
     _VERIFY(STATUS)
     call MAPL_FieldBundleAdd (mergedBUN, FLD, rc=status)
     _VERIFY(STATUS)
   end do ! L
   if(verbose .and. MAPL_AM_I_ROOT()) &
   print *,' Add ',NumVars2,' fields from ',trim(name2), &
   ' into ', trim(name3)
   do L=1,NumVars2
     call ESMF_FieldBundleGet (BUN2, L, FLD, RC=STATUS)
     _VERIFY(STATUS)
     call MAPL_FieldBundleAdd (mergedBUN, FLD, rc=status)
     _VERIFY(STATUS)
   end do ! L

   _RETURN(ESMF_SUCCESS)

   end subroutine Bundles2Bundle

!-------------------------------------------------------------------
   subroutine Add2Bundle (BUN, mergedBUN, rc)
!-------------------------------------------------------------------
   type(ESMF_FieldBundle), intent(inout) :: BUN
   type(ESMF_FieldBundle), intent(inout) :: mergedBUN
   integer, intent(out), optional :: rc
!
   type(ESMF_VM)           :: vm
   type(ESMF_Field)        :: FLD
   integer :: status, mype, L, numVars1, numVars2
   character(len=ESMF_MAXSTR)            :: name1, name2
   character(len=ESMF_MAXSTR), parameter :: IAm = 'Add2Bundle'

! begin

   call ESMF_VMGetCurrent(vm)
   call ESMF_VMGet(vm, localPet=mype, rc=status)
   _VERIFY(STATUS)

   ! get number of fields in bundles
   call ESMF_FieldBundleGet (BUN, name=name1, FieldCount=NumVars1, RC=STATUS)
   _VERIFY(STATUS)

   ! loop over fields in bundle and add to mergedBUN

   do L=1,NumVars1
     call ESMF_FieldBundleGet (BUN, L, FLD, RC=STATUS)
     _VERIFY(STATUS)
     call MAPL_FieldBundleAdd (mergedBUN, FLD, rc=status)
     _VERIFY(STATUS)
   end do ! L
   call ESMF_FieldBundleGet (mergedBUN, name=name2, FieldCount=NumVars2, RC=STATUS)
   _VERIFY(STATUS)
   if(verbose .and. MAPL_AM_I_ROOT()) &
   print *,' Added ',NumVars1,' fields from bundle ',trim(name1), &
   ' into ', trim(name2), ' for a total number of fields of ', NumVars2

   _RETURN(ESMF_SUCCESS)

   end subroutine Add2Bundle

! THE FOLLOWING ROUTINE IS TEMPORARY

!-------------------------------------------------------------------------
   subroutine ESMFL_FailedRC(mype, name)
!-------------------------------------------------------------------------
!
! A local error-exiting routine for ESMFL
!
   integer, intent(in) :: mype
   character(len=*), intent(in) :: name

   _UNUSED_DUMMY(mype)

   write(*,*) ' *** ',trim(name),' failed ***'
   call ESMF_Finalize()
   end subroutine ESMFL_FailedRC

!-------------------------------------------------------------------------
  SUBROUTINE ESMFL_HALO_R4_2D(GRID, INPUT, RC)
!-------------------------------------------------------------------------
     use MAPL_GridManagerMod, only: get_factory
     use MAPL_AbstractGridFactoryMod
     TYPE(ESMF_Grid),         INTENT(INout) :: GRID
     REAL,                    INTENT(INOUT) :: INPUT(:,:)
     integer, optional,       intent(OUT)   :: RC

! Local variables

     integer                    :: STATUS



     INTEGER :: LAST, LEN1, LEN2

     integer :: nn_north, nn_south, nn_east, nn_west
     integer :: j_north, j_south, i_east, i_west
     integer, save                            :: NX=-1
     integer, save                            :: NY=-1
     integer                                  :: NX0, NY0
     type (ESMF_DELayout)                     :: layout
     integer                                  :: myid
     integer                                  :: ndes
     integer                                  :: dimCount
     integer, allocatable                     :: minindex(:,:)
     integer, allocatable                     :: maxindex(:,:)
     integer, pointer                         :: ims(:) => null()
     integer, pointer                         :: jms(:) => null()
     type (ESMF_VM)                           :: VM
     type (ESMF_DistGrid)                     :: distGrid
     integer                                  :: IM, JM, DIMS(3)

#define MAX_HALOTYPES 8
     type HaloType
        character(len=ESMF_MAXSTR)            :: gridname = ""
        type(ESMF_DELayout)                   :: layout
        integer                               :: NX = -1
        integer                               :: NY = -1
        integer                               :: domainIdx = -1
        integer                               :: myId
        logical                               :: isCube = .false.
        logical                               :: inUse = .false.
     end type HaloType

     type(HaloType), save                     :: myHaloTypes(MAX_HALOTYPES)
     type(HaloType)                           :: thisHaloType
     logical                                  :: found
     integer                                  :: I
     integer                                  :: isd, ied
     integer                                  :: jsd, jed
     character(len=ESMF_MAXSTR)               :: gridname
     class (AbstractGridFactory), pointer :: factory
     call ESMF_GridGet    (GRID,   name=gridname, RC=STATUS)
     _VERIFY(STATUS)

     found = .false.
     DO I = 1, MAX_HALOTYPES
        if (gridname == myHaloTypes(I)%gridname) then
           found = .true.
           thisHaloType = myHaloTypes(I)
           exit
        end if
     END DO

     if (.not. found) then

! try to find unused slot to save this Halo for future use
        DO I = 1, MAX_HALOTYPES
           if(.not.myHaloTypes(I)%inUse) then
              found = .true.
              thisHaloType%domainIdx = I
              exit
           end if
        END DO

        if (.not.found) then
           print *, "Error: need bigger MAX_HALOTYPES value"
           _FAIL( 'no unused slot for halo types')
        end if

        call ESMF_GridGet(GRID,   distGrid=distGrid, dimCount=dimCount, RC=STATUS)
        _VERIFY(STATUS)
        call ESMF_DistGridGet(distGRID, delayout=layout, rc=STATUS)
        _VERIFY(STATUS)

        call ESMF_DELayoutGet (layout, VM=vm, RC=STATUS)
        _VERIFY(STATUS)

        call ESMF_VmGet(VM, localPet=MYID, petCount=ndes, rc=status)
        _VERIFY(STATUS)

        thisHaloType%myId = MyId

!    Layout sizes
!    ------------

        allocate(minindex(dimCount,ndes), maxindex(dimCount,ndes), stat=status)
        _VERIFY(STATUS)

! Processors in each direction
!-----------------------------

        call MAPL_DistGridGet(distgrid, &
             minIndex=minindex, &
             maxIndex=maxindex, rc=status)
        _VERIFY(STATUS)

        call MAPL_GetImsJms(Imins=minindex(1,:),Imaxs=maxindex(1,:),&
             Jmins=minindex(2,:),Jmaxs=maxindex(2,:),Ims=ims,Jms=jms,rc=status)
        _VERIFY(STATUS)

        NX = size(ims)
        NY = size(jms)

        deallocate(jms, ims)
        deallocate(maxindex, minindex)

        thisHaloType%gridname = gridname
        thisHaloType%inUse    = .true.
        thisHaloType%NX = NX
        thisHaloType%NY = NY
        thisHaloType%layout = layout

        call MAPL_GridGet(grid, globalCellCountPerDim=DIMS, RC=status)
        _VERIFY(STATUS)
        IM = DIMS(1)
        JM = DIMS(2)
        if (JM == 6*IM) then
           thisHaloType%isCube = .true.
        else
           thisHaloType%isCube = .false.
        end if

        myHaloTypes(thisHaloType%domainIdx) = thisHaloType

     end if

     if (thisHaloType%isCube) then
        factory => get_factory(grid, rc=status)
        _VERIFY(status)
        isd = lbound(input,1)
        ied = ubound(input,1)
        jsd = lbound(input,2)
        jed = ubound(input,2)
        ! fill the corners with MAPL_UNDEF, they will be overwritten
        ! by the CubeHalo, unless these are the actual cube corners
        input(isd,jsd) = MAPL_UNDEF
        input(ied,jsd) = MAPL_UNDEF
        input(isd,jed) = MAPL_UNDEF
        input(ied,jed) = MAPL_UNDEF
        call factory%halo(input, rc=status)
        _VERIFY(status)
        _RETURN(ESMF_SUCCESS)
     end if

! This is section of the code is valid ONLY for rectilinear grids!

     NX = thisHaloType%NX
     NY = thisHaloType%NY
     myId = thisHaloType%myId
     layout = thisHaloType%layout

!    My processor coordinates
!    -------------------------

!ALT: we will use 0-based NX0, NY0 here, unlike in MAPL_Generic !

     NX0 = mod(MYID,NX)
     NY0 = MYID/NX

     j_north = mod(NY0   +1,NY)
     j_south = mod(NY0+NY-1,NY)
     i_east  = MOD(NX0   +1,NX)
     i_west  = MOD(NX0+NX-1,NX)

     ! Nearest neighbors processor' ids

     nn_north = NX0 + NX*j_north
     nn_south = NX0 + NX*j_south
     nn_west  = i_west + NX*NY0
     nn_east  = i_east + NX*NY0

!  Fill NORTHERN ghost region

     LAST = SIZE(INPUT,2)-1
     LEN1 = SIZE(INPUT,1)

     call MAPL_CommsSendRecv(LAYOUT, &
          INPUT(:,2        ),  LEN1,  NN_SOUTH,  &
          INPUT(:,last+1   ),  LEN1,  NN_NORTH,  &
          RC=STATUS)
     _VERIFY(STATUS)
     if(NY0==NY-1) then
        INPUT(:,last+1   ) = INPUT(:,last )
     end if

!  Fill SOUTHERN ghost region

     CALL MAPL_CommsSendRecv(layout,             &
          INPUT(:,last     ),  LEN1,  NN_NORTH,  &
          INPUT(:,1        ),  LEN1,  NN_SOUTH,  &
          RC=STATUS)
     _VERIFY(STATUS)
     if(NY0==0) then
        INPUT(:,1   ) = INPUT(:,2 )
     endif

     LAST = SIZE(INPUT,1)-1
     LEN2 = SIZE(INPUT,2)

!  Fill EASTERN ghost region

     CALL MAPL_CommsSendRecv(layout,           &
          INPUT(2     , : ),  LEN2,  NN_WEST,  &
          INPUT(last+1, : ),  LEN2,  NN_EAST,  &
          RC=STATUS)
     _VERIFY(STATUS)

!  Fill WESTERN ghost region

     CALL MAPL_CommsSendRecv(layout,           &
          INPUT(last  , : ),  LEN2,  NN_EAST,  &
          INPUT(1     , : ),  LEN2,  NN_WEST,  &
          RC=STATUS)
     _VERIFY(STATUS)

     _RETURN(ESMF_SUCCESS)

   end SUBROUTINE ESMFL_HALO_R4_2D

!---------------------------------------------------------------------
!>
! The recursive subrountine `BundleAddState_` adds contents of State to Bundle.
! Extracts fields from an ESMF State and adds them to a
! ESMF Bundle. In essesence, it serializes an ESMF state
! in a flat Bundle. The BUNDLE must have been created prior to calling
! this routine.
!
    RECURSIVE subroutine BundleAddState_ ( BUNDLE, STATE, rc, &
                                           GRID, VALIDATE )
!
! !ARGUMENTS:
!
    implicit NONE
    type(ESMF_FieldBundle), intent(inout)         :: BUNDLE
    type(ESMF_State),  intent(INout)            :: STATE
    integer, optional                        :: rc
    type(ESMF_Grid),  optional, intent(in)   :: GRID
    logical, optional, intent(in)            :: VALIDATE
!
    character(len=*), parameter          :: Iam="ESMFL_StateSerialize"
    integer                              :: STATUS

    type(ESMF_State)                     :: tSTATE
    type(ESMF_FieldBundle)                    :: tBUNDLE
    type(ESMF_Field)                     :: tFIELD
    type(ESMF_Grid)                      :: tGRID

    integer                              :: I, J
    integer                              :: ItemCount, FieldCount
    type (ESMF_StateItem_Flag), pointer  :: ItemTypes(:)
    character(len=ESMF_MAXSTR ), pointer :: ItemNames(:), FieldNames(:)
    logical                              :: needGrid = .true.
    logical                              :: validate_ = .false.
    integer, allocatable :: orderlist(:)
    integer :: jj
    character(len=ESMF_MAXSTR)           :: attrName
    character(len=ESMF_MAXSTR), allocatable :: currList(:)
    integer                                 :: natt

!                           ---


    if ( present(validate) ) validate_ = validate

!   Query state for number of items and allocate space for them
!   -----------------------------------------------------------
    call ESMF_StateGet(STATE,ItemCount=ItemCount,RC=STATUS)
    _VERIFY(STATUS)
    _ASSERT(ItemCount>0, 'itemCount should be > 0')
    allocate ( ItemNames(ItemCount), stat=STATUS)
    _VERIFY(STATUS)
    allocate ( ItemTypes(ItemCount), stat=STATUS)
    _VERIFY(STATUS)

!   Next, retrieve the names and types of each item in the state
!   ------------------------------------------------------------
    call ESMF_StateGet ( STATE,      ItemNameList = ItemNames, &
                                ItemtypeList = ItemTypes, &
                                rc=STATUS)
    _VERIFY(STATUS)

!   Loop over each item on STATE
!   ----------------------------
    attrName = MAPL_StateItemOrderList
    call ESMF_AttributeGet(state, NAME=attrName, itemcount=natt, RC=STATUS)
    _VERIFY(STATUS)

    _ASSERT(natt > 0, 'natt should be > 0')
    allocate(orderlist(natt), stat=status)
    _VERIFY(STATUS)
    allocate(currList(natt), stat=status)
    _VERIFY(STATUS)

    ! get the current list
    call ESMF_AttributeGet(state, NAME=attrName, VALUELIST=currList, rc=status)
    _VERIFY(STATUS)

    orderList = -1 ! not found
    do i = 1, natt
       ! search loop
       do jj = 1, ITEMCOUNT
          if (itemNames(jj) == currList(i)) then
             orderList(i) = jj
             exit
          end if
       end do
    end do

    deallocate(currList)

    do JJ = 1, natt

       I = ORDERLIST(JJ)
!@@@    do I = 1, ItemCount

!         Got a field
!         -----------
          if (ItemTypes(I) == ESMF_StateItem_Field) THEN

             call ESMF_StateGet ( STATE, ItemNames(i), tFIELD, rc=status)
             _VERIFY(STATUS)
             call AddThisField_()

!         Got a Bundle
!         ------------
          else if (ItemTypes(I) == ESMF_StateItem_FieldBundle) then

             call ESMF_StateGet(STATE, ItemNames(i), tBUNDLE, rc=STATUS)
             _VERIFY(STATUS)
             call ESMF_FieldBundleGet ( tBUNDLE, FieldCount = FieldCount, rc=STATUS)
             _VERIFY(STATUS)
             _ASSERT(FieldCount>0, 'FieldCount should be > 0')
             do j = 1, FieldCount
                call ESMF_FieldBundleGet ( tBUNDLE, j, tFIELD, rc=STATUS)
                _VERIFY(STATUS)
                call AddThisField_()
             end do

!         Got another state
!         -----------------
          else if (ItemTypes(I) == ESMF_StateItem_State) then

             call ESMF_StateGet(STATE, ItemNames(i), tSTATE, rc=STATUS)
             _VERIFY(STATUS)
             call BundleAddState_ ( BUNDLE, tSTATE, rc=STATUS )
             _VERIFY(STATUS)
             if (needGrid) then
                call ESMF_FieldBundleGet ( BUNDLE, GRID=tGRID, rc=STATUS )
                _VERIFY(STATUS)
                needGrid = .false.
             end if

!         What is this?
!         ------------
          else
             cycle
          end IF

    end do

    deallocate(orderlist)

!   Make sure the Bundle is not empty
!   ---------------------------------
    call ESMF_FieldBundleGet ( BUNDLE, FieldCount = FieldCount, rc=STATUS)
    _VERIFY(STATUS)
    _ASSERT(FieldCount>0, 'FieldCount should be > 0')

!   Set the grid
!   ------------
    if ( present(GRID) ) then
       call ESMF_FieldBundleSet ( BUNDLE, grid=GRID, rc=STATUS )
       _VERIFY(STATUS)
       needGrid = .false.
    else
       _ASSERT(.not. needGrid, 'could not find a grid')
    end if

!   Make sure field names are unique
!   --------------------------------
    allocate ( FieldNames(FieldCount), stat=STATUS )
    _VERIFY(STATUS)
    call ESMF_FieldBundleGet ( BUNDLE, FieldNameList=FieldNames, rc=STATUS )
    _VERIFY(STATUS)

!   Make sure field names are unique
!   --------------------------------
    if ( validate_ ) then
       do j = 1, FieldCount
          do i = j+1, FieldCount
             if ( trim(FieldNames(i)) == trim(FieldNames(j)) ) then
                STATUS = -1              ! same name
                _VERIFY(STATUS)
             end if
          end do
       end do
    end if

!   All done
!   --------
    deallocate(ItemNames)
    deallocate(ItemTypes)
    deallocate(FieldNames)

    _RETURN(ESMF_SUCCESS)

CONTAINS

    subroutine AddThisField_()

      call MAPL_FieldBundleAdd ( BUNDLE, tField, rc=STATUS )
      _VERIFY(STATUS)
      if ( needGrid ) then
         call ESMF_FieldGet ( tFIELD, grid=tGRID, rc=STATUS )
         _VERIFY(STATUS)
         needGrid = .false.
      end if
    end subroutine AddThisField_

  end subroutine BundleAddState_

  subroutine MAPL_AreaMean_2d_r8_bitrep ( qave, q, area, grid, bitreproducible, rc )

    implicit none

    ! Arguments
    real(kind=REAL64), intent(  OUT) :: qave
    real,              intent(IN   ) :: q(:,:)
    real,              intent(IN   ) :: area(:,:)
    type(ESMF_Grid),   intent(INout) :: grid
    logical,           intent(IN   ) :: bitreproducible
    integer, optional, intent(  OUT) :: rc

    ! Log err vars
    integer :: status
    character(len=ESMF_MAXSTR), parameter :: Iam='MAPL_AreaMeanBR'

    ! Local vars
    real(kind=REAL64)  :: qdum(2)
    integer :: im,jm
    integer :: DIMS(3)

    integer :: i,j
    logical :: amIRoot

    real, allocatable :: qglobal(:,:)
    real, allocatable :: aglobal(:,:)

    type(ESMF_VM) :: vm

    if (.not. bitreproducible) then
       call MAPL_AreaMean(qave, q, area, grid, rc=status )
       _RETURN(STATUS)
    end if

    ! get VM (should get from the grid, but this is quicker)
    call ESMF_VmGetCurrent(vm, rc=status)
    _VERIFY(STATUS)

    amIRoot = MAPL_AM_I_Root(vm)

    if (amIRoot) then
       call MAPL_GridGet(GRID, globalCellCountPerDim=DIMS, RC=STATUS)
       im = DIMS(1) ! global grid dim
       jm = DIMS(2) ! global grid dim
    else
       im = 0 ! dummy sizes
       jm = 0 ! dummy sizes
    end if

    allocate(qglobal(im,jm), stat=status)
    _VERIFY(STATUS)
    allocate(aglobal(im,jm), stat=status)
    _VERIFY(STATUS)

    call ArrayGather(local_array=area, global_array=aglobal, grid=grid, rc=status)
    _VERIFY(STATUS)
    call ArrayGather(local_array=q,    global_array=qglobal, grid=grid, rc=status)
    _VERIFY(STATUS)
    qdum = 0.0d+0
    ! do calculation on root
    if (amIRoot) then
       do j=1,jm
          do i=1,im
             if (qglobal(i,j) == MAPL_Undef) cycle ! exclude any undefs
             qdum(1) = qdum(1) + qglobal(i,j)*aglobal(i,j)
             qdum(2) = qdum(2) + aglobal(i,j)
          enddo
       end do

       if (qdum(2) /= 0.0d+0) then

          qave = qdum(1) / qdum(2)

          !ALT: convert the the result to single precision
          !     This technically is not needed here, but is is done to be in sync
          !     with the parallel method below
          !     qave = real(qave, kind=4)
       else
          qave = MAPL_Undef
       end if
    end if
    deallocate(aglobal)
    deallocate(qglobal)

    call MAPL_CommsBcast(vm, DATA=qave, N=1, Root=MAPL_Root, RC=status)
    _VERIFY(STATUS)

    _RETURN(ESMF_SUCCESS)
  end subroutine MAPL_AreaMean_2d_r8_bitrep

  subroutine MAPL_AreaMean_2d_r8 ( qave, q, area, grid, rc )

    implicit none

    ! Arguments
    real(kind=REAL64), intent(  OUT) :: qave
    real,              intent(IN   ) :: q(:,:)
    real,              intent(IN   ) :: area(:,:)
    type(ESMF_Grid),   intent(INout) :: grid
    integer, optional, intent(  OUT) :: rc

    ! Log err vars
    integer :: status
    character(len=ESMF_MAXSTR), parameter :: Iam='MAPL_AreaMean'

    ! Local vars
    real(kind=REAL64)  :: qdum(2)
    real(kind=REAL64)  :: qdumloc(2)
    integer :: im,jm

    integer :: i,j

    type(ESMF_VM) :: vm

    _UNUSED_DUMMY(grid)

    ! get VM (should get from the grid, but this is quicker)
    call ESMF_VmGetCurrent(vm, rc=status)
    _VERIFY(STATUS)

    im = size(area,1) ! local grid dim
    jm = size(area,2) ! local grid dim

    ! do calculation on every PE

    ! compute local sum
    qdumloc = 0.0d+0
    do j=1,jm
       do i=1,im
          if (q(i,j) == MAPL_Undef) cycle ! exclude any undefs
          qdumloc(1) = qdumloc(1) + q(i,j)*area(i,j)
          qdumloc(2) = qdumloc(2) + area(i,j)
       enddo
    end do

    call MAPL_CommsAllReduceSum(vm, sendbuf=qdumloc, recvbuf=qdum, &
         cnt=2, RC=status)
    _VERIFY(STATUS)

    if (qdum(2) /= 0.0d+0) then

       qave = qdum(1) / qdum(2)

       !ALT: convert the the result to single precision to get rid of
       !     numerical non-associativity in floating point numbers
       !     qave = real(qave, kind=4)
    else
       qave = MAPL_Undef
    end if

    _RETURN(ESMF_SUCCESS)
  end subroutine MAPL_AreaMean_2d_r8

  logical function check_list_ (this,list)
  implicit none
  character(len=*), intent(in) :: this
  character(len=*), intent(in) :: list
  integer i,ll,ls,le
  check_list_=.false.
  if (trim(list)=="NONE") then ! this is when does not specify variables
      check_list_=.true.
  endif
  ll  = len_trim(list)
  ls  = 1
  do i=1,ll
     le=i
     if(list(i:i)==",") then
        if (trim(this)==list(ls:le-1)) then
           check_list_=.true.
           exit
        endif
        ls=le+1
     endif
  enddo
  if (trim(this)==list(ls:ll)) then
      check_list_=.true.
  endif
  end function check_list_

  function ESMFL_field_is_undefined(field,rc) result(field_is_undefined)
     use mpi
     logical :: field_is_undefined
     type(ESMF_Field), intent(in) :: field
     integer, optional, intent(out) :: rc

     integer :: status
     type(ESMF_Grid) :: grid
     type(ESMF_VM) :: vm
     integer :: my_mpi_comm, local_undef, global_undef,grid_size(3),rank
     real, pointer :: ptr2d(:,:), ptr3d(:,:,:)

     call ESMF_VMGetCurrent(VM,_RC)
     call ESMF_VMGet(VM,mpiCommunicator=my_mpi_comm,_RC)
     call ESMF_FieldGet(field,rank=rank,grid=grid,_RC)
     call MAPL_GridGet(grid,globalcellcountperdim=grid_size,_RC)
     if (rank ==2) then
        call ESMF_FieldGet(field,0,farrayPtr=ptr2d,_RC)
        local_undef = count(ptr2d==MAPL_UNDEF)
        call MPI_AllReduce(local_undef,global_undef,1,MPI_INTEGER,MPI_SUM,my_mpi_comm,status)
        _VERIFY(status)
        field_is_undefined = global_undef == (grid_size(1)*grid_size(2))
     else if (rank ==3) then
        call ESMF_FieldGet(field,0,farrayPtr=ptr3d,_RC)
        local_undef = count(ptr3d==MAPL_UNDEF)
        call MPI_AllReduce(local_undef,global_undef,1,MPI_INTEGER,MPI_SUM,my_mpi_comm,status)
        _VERIFY(status)
        field_is_undefined = global_undef == (grid_size(1)*grid_size(2)*grid_size(3))
     else
        _FAIL("Unsupported rank when checking for undef")
     end if

     _RETURN(_SUCCESS)
  end function

end module ESMFL_MOD