mpas_geom_mod.F90

Go to the documentation of this file.

! (C) Copyright 2017 UCAR
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 
module mpas_geom_mod
 
use atlas_module, only: atlas_functionspace, atlas_fieldset, atlas_field, atlas_real
use fckit_configuration_module, only: fckit_configuration, fckit_yamlconfiguration
use fckit_pathname_module, only: fckit_pathname
use fckit_log_module, only: fckit_log
use fckit_mpi_module, only: fckit_mpi_comm
use iso_c_binding
 
!oops
use oops_variables_mod, only: oops_variables
 
!ufo
use ufo_vars_mod, only: maxvarlen, ufo_vars_getindex
 
!MPAS-Model
use mpas_derived_types
use mpas_kind_types
use mpas_constants
use kinds, only : kind_real
use mpas_dmpar, only: mpas_dmpar_sum_int
use mpas_subdriver
use atm_core
use mpas_pool_routines
 
!mpas_jedi
use mpas_constants_mod
 
implicit none
private
public :: mpas_geom, &
          geo_setup, geo_clone, geo_delete, geo_info, geo_is_equal, &
          geo_set_atlas_lonlat, geo_fill_atlas_fieldset, pool_has_field, &
          getsolvedimsizes, getsolvedimnames, getvertlevels
 
public :: mpas_geom_registry
 
 
! ------------------------------------------------------------------------------
 
type :: templated_field
   character(len=MAXVARLEN) :: name
   character(len=MAXVARLEN) :: template
   character(len=MAXVARLEN) :: identity
end type templated_field
 
!> Fortran derived type to hold geometry definition
type :: mpas_geom
   integer :: ncellsglobal !Global count
   integer :: nedgesglobal !Global count
   integer :: nverticesglobal !Global count
   integer :: ncells !Memory count (Local + Halo)
   integer :: nedges !Memory count (Local + Halo)
   integer :: nvertices !Memory count (Local + Halo)
   integer :: ncellssolve !Local count
   integer :: nedgessolve !Local count
   integer :: nverticessolve !Local count
   integer :: nvertlevels
   integer :: nvertlevelsp1
   integer :: nsoillevels
   integer :: vertexdegree
   integer :: maxedges
   logical :: deallocate_nonda_fields
   logical :: use_bump_interpolation
   character(len=StrKIND) :: bump_vunit
   real(kind=kind_real), dimension(:),   allocatable :: latcell, loncell
   real(kind=kind_real), dimension(:),   allocatable :: areacell
   real(kind=kind_real), dimension(:),   allocatable :: latedge, lonedge
   real(kind=kind_real), dimension(:,:), allocatable :: edgenormalvectors
   real(kind=kind_real), dimension(:,:), allocatable :: zgrid
   integer, allocatable :: nedgesoncell(:)
   integer, allocatable :: cellsoncell(:,:)
   integer, allocatable :: edgesoncell(:,:)
   integer, allocatable :: cellsonvertex(:,:)
   integer, allocatable :: cellsonedge(:,:)
   integer, allocatable :: verticesonedge(:,:)
   real(kind=kind_real), DIMENSION(:), ALLOCATABLE :: dcedge, dvedge
   real(kind=kind_real), DIMENSION(:), ALLOCATABLE :: areatriangle, angleedge
   real(kind=kind_real), DIMENSION(:,:), ALLOCATABLE :: kiteareasonvertex, edgesoncell_sign
 
   type (domain_type), pointer :: domain => null()
   type (core_type), pointer :: corelist => null()
 
   type(fckit_mpi_comm) :: f_comm
 
   type(atlas_functionspace) :: afunctionspace
 
   type(templated_field), allocatable :: templated_fields(:)
 
   contains
 
   procedure, public :: is_templated => field_is_templated
   procedure, public :: template => template_fieldname
   procedure, public :: has_identity => field_has_identity
   procedure, public :: identity => identity_fieldname
   generic, public :: nlevels => variables_nlevels, var_nlevels
   procedure :: variables_nlevels
   procedure :: var_nlevels
 
end type mpas_geom
 
type :: idcounter
   integer :: id, counter
end type idcounter
 
type(idcounter), allocatable :: geom_count(:)
 
character(len=1024) :: message
 
character(len=MAXVARLEN), parameter :: &
   mpasverticalcoordinates(6) = &
      [character(len=maxvarlen) :: 'nVertLevels', 'nVertLevelsP1', &
         'nVertLevelsP2', 'nSoilLevels', 'nOznLevels', 'nAerLevels']
 
#define LISTED_TYPE mpas_geom
 
!> Linked list interface - defines registry_t type
#include <oops/util/linkedList_i.f>
 
!> Global registry
type(registry_t) :: mpas_geom_registry
 
! ------------------------------------------------------------------------------
contains
! ------------------------------------------------------------------------------
!> Linked list implementation
#include <oops/util/linkedList_c.f>
 
! ------------------------------------------------------------------------------
subroutine geo_setup(self, f_conf, f_comm)
 
   implicit none
 
   type(mpas_geom),           intent(inout) :: self
   type(fckit_configuration), intent(in)    :: f_conf
   type(fckit_mpi_comm),   intent(in)    :: f_comm
 
   character(len=StrKIND) :: string1
   type (mpas_pool_type), pointer :: meshpool, fg
   type (block_type), pointer :: block_ptr
   !character(len=120) :: fn
   character(len=512) :: nml_file, streams_file, fields_file
   character(len=:), allocatable :: str
   type(fckit_configuration) :: template_conf
   type(fckit_configuration), allocatable :: fields_conf(:)
 
   real (kind=kind_real), pointer :: r1d_ptr(:), r2d_ptr(:,:)
   integer, pointer :: i0d_ptr, i1d_ptr(:), i2d_ptr(:,:)
 
   type(idcounter), allocatable :: prev_count(:)
   integer :: ii, nprev
 
   logical :: deallocate_fields
   logical :: bump_interp
 
   call fckit_log%info('==> create geom')
 
   ! MPI communicator
   self%f_comm = f_comm
 
   ! "nml_file" and "streams_file" are mandatory for mpas_init.
   call f_conf%get_or_die("nml_file",str)
   nml_file = str
   call f_conf%get_or_die("streams_file",str)
   streams_file = str
 
   ! Domain decomposition and templates for state/increment variables
   call mpas_init( self % corelist, self % domain, mpi_comm = self%f_comm%communicator(), &
                 & namelistfileparam = trim(nml_file), streamsfileparam = trim(streams_file))
 
   ! The interpolation method specified here will be used to interpolate data between
   ! different geometries, but not in GetValues. GetValues has its own configuration
   ! for specifying its interpolation method.
   if (f_conf%get("interpolation type", str)) then
     select case (str)
       case ('bump')
         self%use_bump_interpolation = .true.
       case ('unstructured')
         self%use_bump_interpolation = .false.
       case default
         write(message,*) '--> geo_setup: interpolation type: ',str,' not implemented'
         call abor1_ftn(message)
     end select
   else
     self%use_bump_interpolation = .true. ! BUMP is default interpolation
   end if
 
   !Deallocate not-used fields for memory reduction
   if (f_conf%has("deallocate non-da fields")) then
      call f_conf%get_or_die("deallocate non-da fields",deallocate_fields)
      self % deallocate_nonda_fields = deallocate_fields
   else
      self % deallocate_nonda_fields = .false.
   end if
   if (self % deallocate_nonda_fields) call geo_deallocate_nonda_fields (self % domain)
 
   ! Set up the vertical coordinate for bump
   if (f_conf%has("bump vunit")) then
      call f_conf%get_or_die("bump vunit",str)
      self % bump_vunit = str
   else
      self % bump_vunit = 'modellevel'
   end if
 
   if (allocated(geom_count)) then
      nprev = size(geom_count)
      allocate(prev_count(nprev))
      do ii = 1, nprev
         prev_count(ii) = geom_count(ii)
         if (prev_count(ii)%id == self%domain%domainID) then
            call abor1_ftn("domainID already used")
         end if
      end do
      deallocate(geom_count)
   else
      nprev = 0
   end if
   allocate(geom_count(nprev+1))
   do ii = 1, nprev
      geom_count(ii) = prev_count(ii)
   end do
   geom_count(nprev+1)%id = self%domain%domainID
   geom_count(nprev+1)%counter = 1
 
   ! first read array of templated field configurations
   if (f_conf%get('template fields file',str)) then
     fields_file = str
   else
     fields_file = 'geovars.yaml'
   end if
   template_conf = fckit_yamlconfiguration(fckit_pathname(fields_file))
   call template_conf%get_or_die('fields',fields_conf)
 
   ! create templated field descriptors
   allocate(self%templated_fields(size(fields_conf)))
   do ii = 1, size(fields_conf)
      call fields_conf(ii)%get_or_die('field name',str)
      self%templated_fields(ii)%name = trim(str)
      call fields_conf(ii)%get_or_die('mpas template field',str)
      self%templated_fields(ii)%template = trim(str)
      if (fields_conf(ii)%get('mpas identity field',str)) then
        self%templated_fields(ii)%identity = trim(str)
      else
        self%templated_fields(ii)%identity = 'none'
      end if
   end do
   deallocate(fields_conf)
 
!   if (associated(self % domain)) then
!       write(message,*) 'inside geom: geom % domain associated for domainID = ', self % domain % domainID
!       call fckit_log%debug(message)
!   end if
!   if (associated(self % corelist)) then
!       call fckit_log%debug('inside geom: geom % corelist associated')
!   else
!       call fckit_log%debug('inside geom: geom % corelist not associated')
!   end if
 
   !  These pool accesses refer to memory (local+halo) for a single MPAS block (standard)
   block_ptr => self % domain % blocklist
 
   call mpas_pool_get_subpool ( block_ptr % structs, 'mesh', meshpool )
 
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'nCells', i0d_ptr )
   self % nCells = i0d_ptr
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'nCellsSolve', i0d_ptr )
   self % nCellsSolve = i0d_ptr
   call mpas_dmpar_sum_int ( self % domain % dminfo, &
                             self % nCellsSolve, self % nCellsGlobal )
 
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'nEdges', i0d_ptr )
   self % nEdges = i0d_ptr
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'nEdgesSolve', i0d_ptr )
   self % nEdgesSolve = i0d_ptr
   call mpas_dmpar_sum_int ( self % domain % dminfo, &
                             self % nEdgesSolve, self % nEdgesGlobal )
 
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'nVertices', i0d_ptr )
   self % nVertices = i0d_ptr
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'nVerticesSolve', i0d_ptr )
   self % nVerticesSolve = i0d_ptr
   call mpas_dmpar_sum_int ( self % domain % dminfo, &
                             self % nVerticesSolve, self % nVerticesGlobal )
 
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'nVertLevels', i0d_ptr )
   self % nVertLevels = i0d_ptr
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'nVertLevelsP1', i0d_ptr )
   self % nVertLevelsP1 = i0d_ptr
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'nSoilLevels', i0d_ptr )
   self % nSoilLevels = i0d_ptr
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'vertexDegree', i0d_ptr )
   self % vertexDegree = i0d_ptr
   call mpas_pool_get_dimension ( block_ptr % dimensions, 'maxEdges', i0d_ptr )
   self % maxEdges = i0d_ptr
 
   allocate ( self % latCell ( self % nCells ) )
   allocate ( self % lonCell ( self % nCells ) )
   allocate ( self % latEdge ( self % nEdges ) )
   allocate ( self % lonEdge ( self % nEdges ) )
   allocate ( self % areaCell ( self % nCells ) )
   allocate ( self % edgeNormalVectors (3,  self % nEdges ) )
   allocate ( self % zgrid ( self % nVertLevelsP1, self % nCells ) )
   allocate ( self % nEdgesOnCell ( self % nCells ) )
   allocate ( self % edgesOnCell ( self % maxEdges, self % nCells ) )
   allocate ( self % cellsOnCell ( self % maxEdges, self % nCells ) )
   allocate ( self % cellsOnVertex ( self % vertexDegree, self % nVertices ) )
   allocate ( self % cellsOnEdge ( 2, self % nEdges ) )
   allocate ( self % verticesOnEdge ( 2, self % nEdges ) )
   allocate ( self % dcEdge ( self % nEdges ) )
   allocate ( self % dvEdge ( self % nEdges ) )
   allocate ( self % kiteAreasOnVertex ( self % vertexDegree, self % nVertices ) )
   allocate ( self % edgesOnCell_sign ( self % maxEdges, self % nCells ) )
   allocate ( self % areaTriangle ( self % nVertices ) )
   allocate ( self % angleEdge ( self % nEdges ) )
 
   call mpas_pool_get_array ( meshpool, 'latCell', r1d_ptr )
   self % latCell = r1d_ptr(1:self % nCells)
   where (self % latCell > mpas_jedi_piio2_kr)
       self % latCell = mpas_jedi_piio2_kr
   end where
   where (self % latCell < - mpas_jedi_piio2_kr)
       self % latCell = - mpas_jedi_piio2_kr
   end where
 
   call mpas_pool_get_array ( meshpool, 'lonCell', r1d_ptr )
   self % lonCell = r1d_ptr(1:self % nCells)
   call mpas_pool_get_array ( meshpool, 'areaCell', r1d_ptr )
   self % areaCell = r1d_ptr(1:self % nCells)
   call mpas_pool_get_array ( meshpool, 'latEdge', r1d_ptr )
   self % latEdge = r1d_ptr(1:self % nEdges)
   call mpas_pool_get_array ( meshpool, 'lonEdge', r1d_ptr )
   self % lonEdge = r1d_ptr(1:self % nEdges)
   call mpas_pool_get_array ( meshpool, 'edgeNormalVectors', r2d_ptr )
   self % edgeNormalVectors = r2d_ptr( 1:3, 1:self % nEdges )
   call mpas_pool_get_array ( meshpool, 'nEdgesOnCell', i1d_ptr )
   self % nEdgesOnCell = i1d_ptr(1:self % nCells)
   call mpas_pool_get_array ( meshpool, 'edgesOnCell', i2d_ptr )
   self % edgesOnCell = i2d_ptr( 1:self % maxEdges, 1:self % nCells )
   call mpas_pool_get_array ( meshpool, 'cellsOnCell', i2d_ptr )
   self % cellsOnCell = i2d_ptr( 1:self % maxEdges, 1:self % nCells )
   call mpas_pool_get_array ( meshpool, 'cellsOnVertex', i2d_ptr )
   self % cellsOnVertex = i2d_ptr( 1:self % vertexDegree, 1:self % nVertices )
   call mpas_pool_get_array ( meshpool, 'cellsOnEdge', i2d_ptr )
   self % cellsOnEdge = i2d_ptr( 1:2, 1:self % nEdges )
   call mpas_pool_get_array ( meshpool, 'verticesOnEdge', i2d_ptr )
   self % verticesOnEdge = i2d_ptr( 1:2, 1:self % nEdges )
   call mpas_pool_get_array ( meshpool, 'dcEdge', r1d_ptr )           
   self % dcEdge = r1d_ptr(1:self % nEdges)
   call mpas_pool_get_array ( meshpool, 'dvEdge', r1d_ptr )           
   self % dvEdge = r1d_ptr(1:self % nEdges)
   call mpas_pool_get_array ( meshpool, 'kiteAreasOnVertex', r2d_ptr )             
   self % kiteAreasOnVertex = r2d_ptr( 1:self % vertexDegree, 1:self % nVertices )
   call mpas_pool_get_array ( meshpool, 'edgesOnCell_sign', r2d_ptr )             
   self % edgesOnCell_sign = r2d_ptr( 1:self % maxEdges, 1:self % nCells )
   call mpas_pool_get_array ( meshpool, 'areaTriangle', r1d_ptr )           
   self % areaTriangle = r1d_ptr(1:self % nVertices)
   call mpas_pool_get_array ( meshpool, 'angleEdge', r1d_ptr )           
   self % angleEdge = r1d_ptr(1:self % nEdges)
 
   call mpas_pool_get_array ( meshpool, 'zgrid', r2d_ptr )
   self % zgrid = r2d_ptr( 1:self % nVertLevelsP1, 1:self % nCells )
 
   call fckit_log%debug('End of geo_setup')
   if (allocated(prev_count)) deallocate(prev_count)
   if (allocated(str)) deallocate(str)
 
end subroutine geo_setup
 
! --------------------------------------------------------------------------------------------------
 
subroutine geo_set_atlas_lonlat(self, afieldset)
 
   implicit none
 
   type(mpas_geom),  intent(inout) :: self
   type(atlas_fieldset), intent(inout) :: afieldset
 
   real(kind_real), pointer :: real_ptr(:,:)
   type(atlas_field) :: afield
 
   ! Create lon/lat field
   afield = atlas_field(name="lonlat", kind=atlas_real(kind_real), shape=(/2,self%nCellsSolve/))
   call afield%data(real_ptr)
   real_ptr(1,:) = self%lonCell(1:self%nCellsSolve) * mpas_jedi_rad2deg_kr
   real_ptr(2,:) = self%latCell(1:self%nCellsSolve) * mpas_jedi_rad2deg_kr
   call afieldset%add(afield)
 
end subroutine geo_set_atlas_lonlat
 
! --------------------------------------------------------------------------------------------------
 
subroutine geo_fill_atlas_fieldset(self, afieldset)
 
   implicit none
 
   type(mpas_geom),  intent(inout) :: self
   type(atlas_fieldset), intent(inout) :: afieldset
 
   integer :: i, jz
   real(kind=kind_real), pointer :: real_ptr_1(:), real_ptr_2(:,:)
   type(atlas_field) :: afield
 
   ! Add area
   afield = self%afunctionspace%create_field(name='area', kind=atlas_real(kind_real), levels=0)
   call afield%data(real_ptr_1)
   real_ptr_1 = self%areaCell(1:self%nCellsSolve)
   call afieldset%add(afield)
   call afield%final()
 
   ! Add vertical unit
   afield = self%afunctionspace%create_field(name='vunit', kind=atlas_real(kind_real), levels=self%nVertLevels)
   call afield%data(real_ptr_2)
   do jz=1,self%nVertLevels
      if (trim(self % bump_vunit) .eq. 'modellevel') then
         real_ptr_2(jz,:) = real(jz, kind_real)
      else if (trim(self % bump_vunit) .eq. 'height') then
         real_ptr_2(jz,1:self%nCellsSolve) = mpas_jedi_half_kr * &
         ( self%zgrid(jz,1:self%nCellsSolve) + self%zgrid(jz+1,1:self%nCellsSolve) )
      end if
      !--Similarly for ln of pressure_base. I don't know if we can access to "total pressure" here.
      !real (kind=kind_real), dimension(:,:), pointer :: pressure_base
      !call mpas_pool_get_array(self%domain%blocklist%allFields,'pressure_base', pressure_base)
      !real_ptr_2(jz,1:self%nCellsSolve) = log( pressure_base(jz,1:self%nCellsSolve) )
   end do
   call afieldset%add(afield)
   call afield%final()
 
end subroutine geo_fill_atlas_fieldset
 
! ------------------------------------------------------------------------------
subroutine geo_deallocate_nonda_fields(domain)
 
   implicit none
   type (domain_type), pointer,    intent(inout) :: domain
   type (mpas_pool_type), pointer                :: pool_a, pool_b
   type (mpas_pool_data_type), pointer           :: mem
   type (mpas_pool_iterator_type)                :: poolItr_b
 
   type (field0DReal), pointer :: field0d
   type (field1DReal), pointer :: field1d
   type (field2DReal), pointer :: field2d
 
   integer, parameter :: num_da_fields = 36
   integer            :: i
   character (len=22), allocatable :: poolname_a(:), poolname_b(:)
   character (len=22), allocatable :: da_fieldnames(:)
 
   allocate(poolname_a(3))
   poolname_a(1)='tend'
   poolname_a(2)='tend_physics'
   poolname_a(3)='atm_input'
 
   allocate(poolname_b(3))
   poolname_b(1)='diag_physics'
   poolname_b(2)='diag'
   poolname_b(3)='sfc_input'
 
   !TODO: remove the following list of fieldnames and read from a stream_list.atmosphere file
   allocate(da_fieldnames(num_da_fields))
   da_fieldnames(1)='re_cloud'
   da_fieldnames(2)='re_ice'
   da_fieldnames(3)='re_snow'
   da_fieldnames(4)='cldfrac'
   da_fieldnames(5)='lai'
   da_fieldnames(6)='u10'
   da_fieldnames(7)='v10'
   da_fieldnames(8)='q2'
   da_fieldnames(9)='t2m'
   da_fieldnames(10)='pressure_p'
   da_fieldnames(11)='pressure'
   da_fieldnames(12)='pressure_base'
   da_fieldnames(13)='surface_pressure'
   da_fieldnames(14)='rho'
   da_fieldnames(15)='theta'
   da_fieldnames(16)='temperature'
   da_fieldnames(17)='relhum'
   da_fieldnames(18)='spechum'
   da_fieldnames(19)='u'
   da_fieldnames(20)='uReconstructZonal'
   da_fieldnames(21)='uReconstructMeridional'
   da_fieldnames(22)='stream_function'
   da_fieldnames(23)='velocity_potential'
   da_fieldnames(24)='landmask'
   da_fieldnames(25)='xice'
   da_fieldnames(26)='snowc'
   da_fieldnames(27)='skintemp'
   da_fieldnames(28)='ivgtyp'
   da_fieldnames(29)='isltyp'
   da_fieldnames(30)='snowh'
   da_fieldnames(31)='vegfra'
   da_fieldnames(32)='lai'
   da_fieldnames(33)='smois'
   da_fieldnames(34)='tslb'
   da_fieldnames(35)='vorticity'
   da_fieldnames(36)='w'
 
   do i=1,size(poolname_a)
      mem => pool_get_member(domain % blocklist % structs, poolname_a(i), mpas_pool_subpool)
      if (associated(mem)) then
         call mpas_pool_get_subpool(domain % blocklist % structs, poolname_a(i), pool_a)
         call mpas_pool_destroy_pool(pool_a)
         call mpas_pool_remove_subpool(domain % blocklist % structs, poolname_a(i))
      end if
   end do
 
   do i=1,size(poolname_b)
 
      mem => pool_get_member(domain % blocklist % structs, poolname_b(i), mpas_pool_subpool)
      if (associated(mem)) then
         call mpas_pool_get_subpool(domain % blocklist % structs, poolname_b(i), pool_b)
         call mpas_pool_begin_iteration(pool_b)
 
         do while ( mpas_pool_get_next_member(pool_b, poolitr_b) )
 
            if (poolitr_b % memberType == mpas_pool_field .and. poolitr_b % dataType == mpas_pool_real .and. &
               ufo_vars_getindex(da_fieldnames, poolitr_b % memberName) < 1) then
 
               if (poolitr_b % nDims == 0) then
                  call mpas_pool_get_field(pool_b,trim(poolitr_b % memberName),field0d)
                  if (associated(field0d)) then
                     call mpas_deallocate_field(field0d)
                     call mpas_pool_remove_field(pool_b, trim(poolitr_b % memberName))
                  end if
                  nullify(field0d)
               else if (poolitr_b % nDims == 1) then
                  call mpas_pool_get_field(pool_b,trim(poolitr_b % memberName),field1d)
                  if (associated(field1d)) then
                     call mpas_deallocate_field(field1d)
                     call mpas_pool_remove_field(pool_b, trim(poolitr_b % memberName))
                  end if
                  nullify(field1d)
               else if (poolitr_b % nDims == 2) then
                  call mpas_pool_get_field(pool_b,trim(poolitr_b % memberName),field2d)
                  if (associated(field2d)) then
                     call mpas_deallocate_field(field2d)
                     call mpas_pool_remove_field(pool_b, trim(poolitr_b % memberName))
                  end if
                  nullify(field2d)
               endif
 
            end if
 
         end do
 
      end if
 
   end do
 
   deallocate(poolname_a)
   deallocate(poolname_b)
   deallocate(da_fieldnames)
 
end subroutine geo_deallocate_nonda_fields
 
! ------------------------------------------------------------------------------
 
subroutine geo_clone(self, other)
 
   implicit none
 
   type(mpas_geom), intent(inout) :: self
   type(mpas_geom), intent(in) :: other
 
   integer :: ii
 
   ! Clone communicator
   self%f_comm = other%f_comm
 
   call fckit_log%debug('====> copy of geom array')
 
   self % nCellsGlobal  = other % nCellsGlobal
   self % nCells        = other % nCells
   self % nCellsSolve   = other % nCellsSolve
 
   self % nEdgesGlobal  = other % nEdgesGlobal
   self % nEdges        = other % nEdges
   self % nEdgesSolve   = other % nEdgesSolve
 
   self % nVerticesGlobal  = other % nVerticesGlobal
   self % nVertices        = other % nVertices
   self % nVerticesSolve   = other % nVerticesSolve
 
   self % nVertLevels   = other % nVertLevels
   self % nVertLevelsP1 = other % nVertLevelsP1
   self % nSoilLevels   = other % nSoilLevels
   self % vertexDegree  = other % vertexDegree
   self % maxEdges      = other % maxEdges
 
   if (.not.allocated(self % latCell)) allocate(self % latCell(other % nCells))
   if (.not.allocated(self % lonCell)) allocate(self % lonCell(other % nCells))
   if (.not.allocated(self % latEdge)) allocate(self % latEdge(other % nEdges))
   if (.not.allocated(self % lonEdge)) allocate(self % lonEdge(other % nEdges))
   if (.not.allocated(self % areaCell)) allocate(self % areaCell(other % nCells))
   if (.not.allocated(self % edgeNormalVectors)) allocate(self % edgeNormalVectors(3, other % nEdges))
   if (.not.allocated(self % zgrid)) allocate(self % zgrid(other % nVertLevelsP1, other % nCells))
   if (.not.allocated(self % nEdgesOnCell)) allocate(self % nEdgesOnCell(other % nCells))
   if (.not.allocated(self % edgesOnCell)) allocate(self % edgesOnCell(other % maxEdges, other % nCells))
   if (.not.allocated(self % cellsOnCell)) allocate(self % cellsOnCell (other % maxEdges, other % nCells))
   if (.not.allocated(self % cellsOnVertex)) allocate(self % cellsOnVertex(self % vertexDegree, self % nVertices))
   if (.not.allocated(self % cellsOnEdge)) allocate(self % cellsOnEdge (2, self % nEdges))
   if (.not.allocated(self % verticesOnEdge)) allocate(self % verticesOnEdge (2, self % nEdges))
   if (.not.allocated(self % dcEdge)) allocate(self % dcEdge (self % nEdges))
   if (.not.allocated(self % dvEdge)) allocate(self % dvEdge (self % nEdges))
   if (.not.allocated(self % kiteAreasOnVertex)) allocate(self % kiteAreasOnVertex(self % vertexDegree, self % nVertices))
   if (.not.allocated(self % edgesOnCell_sign)) allocate(self % edgesOnCell_sign(self % maxEdges, self % nCells))
   if (.not.allocated(self % areaTriangle)) allocate(self % areaTriangle(self % nVertices))
   if (.not.allocated(self % angleEdge)) allocate(self % angleEdge(self % nEdges))
 
   self % use_bump_interpolation = other % use_bump_interpolation
   self % templated_fields  = other % templated_fields
   self % latCell           = other % latCell
   self % lonCell           = other % lonCell
   self % areaCell          = other % areaCell
   self % latEdge           = other % latEdge
   self % lonEdge           = other % lonEdge
   self % edgeNormalVectors = other % edgeNormalVectors
   self % zgrid             = other % zgrid
   self % nEdgesOnCell      = other % nEdgesOnCell
   self % edgesOnCell       = other % edgesOnCell
   self % cellsOnCell       = other % cellsOnCell
   self % cellsOnVertex     = other % cellsOnVertex
   self % cellsOnEdge       = other % cellsOnEdge
   self % verticesOnEdge    = other % verticesOnEdge
   self % dcEdge            = other % dcEdge
   self % dvEdge            = other % dvEdge
   self % kiteAreasOnVertex = other % kiteAreasOnVertex
   self % edgesOnCell_sign  = other % edgesOnCell_sign
   self % areaTriangle      = other % areaTriangle
   self % angleEdge         = other % angleEdge
 
   call fckit_log%debug('====> copy of geom corelist and domain')
 
   self % corelist => other % corelist
   self % domain   => other % domain
   write(message,*) 'inside geo_clone: other % domain % domainID = ', other % domain % domainID
   call fckit_log%debug(message)
 
   do ii = 1, size(geom_count)
      if (geom_count(ii)%id == self%domain%domainID) then
         geom_count(ii)%counter = geom_count(ii)%counter + 1
      end if
   end do
 
   call fckit_log%debug('====> copy of geom done')
 
end subroutine geo_clone
 
! ------------------------------------------------------------------------------
 
subroutine geo_delete(self)
 
   implicit none
 
   type(mpas_geom), intent(inout) :: self
   integer :: ii
 
   if (allocated(self%templated_fields)) deallocate(self%templated_fields)
   if (allocated(self%latCell)) deallocate(self%latCell)
   if (allocated(self%lonCell)) deallocate(self%lonCell)
   if (allocated(self%latEdge)) deallocate(self%latEdge)
   if (allocated(self%lonEdge)) deallocate(self%lonEdge)
   if (allocated(self%areaCell)) deallocate(self%areaCell)
   if (allocated(self%edgeNormalVectors)) deallocate(self%edgeNormalVectors)
   if (allocated(self%zgrid)) deallocate(self%zgrid)
   if (allocated(self%nEdgesOnCell)) deallocate(self%nEdgesOnCell)
   if (allocated(self%edgesOnCell)) deallocate(self%edgesOnCell)
   if (allocated(self%cellsOnCell)) deallocate(self%cellsOnCell)
   if (allocated(self%cellsOnVertex)) deallocate(self%cellsOnVertex)
   if (allocated(self%cellsOnEdge)) deallocate(self%cellsOnEdge)
   if (allocated(self%verticesOnEdge)) deallocate(self%verticesOnEdge)
   if (allocated(self%dcEdge)) deallocate(self%dcEdge)
   if (allocated(self%dvEdge)) deallocate(self%dvEdge)
   if (allocated(self%kiteAreasOnVertex)) deallocate(self%kiteAreasOnVertex)
   if (allocated(self%edgesOnCell_sign)) deallocate(self%edgesOnCell_sign)
   if (allocated(self%areaTriangle)) deallocate(self%areaTriangle)
   if (allocated(self%angleEdge)) deallocate(self%angleEdge)
 
   do ii = 1, size(geom_count)
      if (geom_count(ii)%id == self%domain%domainID) then
         geom_count(ii)%counter = geom_count(ii)%counter - 1
         if ((associated(self % corelist)).and.(associated(self % domain))) then
            if (geom_count(ii)%counter == 0) then
               ! Completely destroy corelist and domain
               !  + can only do this if they are not used by another copy of self
               !  + otherwise the persistance of the underlying domain
               !    and corelist is a known memory leak
               write(message,'(A,I3)') &
                  '==> destruct MPAS corelist and domain:', self%domain%domainID
               call fckit_log%info(message)
               call mpas_timer_set_context( self % domain )
               call mpas_finalize(self % corelist, self % domain)
            else
               nullify(self % corelist)
               nullify(self % domain)
            end if
            exit
         end if
      end if
   end do
 
end subroutine geo_delete
 
! ------------------------------------------------------------------------------
 
subroutine geo_is_equal(is_equal, self, other)
 
   implicit none
 
   logical*1,       intent(inout) :: is_equal
   type(mpas_geom), intent(in)    :: self
   type(mpas_geom), intent(in)    :: other
 
   ! Only compares two of many attributes of the two geometries, but
   ! sufficient for current needs.
   is_equal = (self%nCellsGlobal .eq. other%nCellsGlobal) .and. &
              (self%nVertLevels  .eq. other%nVertLevels)
 
end subroutine geo_is_equal
 
! ------------------------------------------------------------------------------
 
subroutine geo_info(self, nCellsGlobal, nCells, nCellsSolve, &
                          nEdgesGlobal, nEdges, nEdgesSolve, &
                          nVertLevels, nVertLevelsP1)
 
   implicit none
 
   type(mpas_geom), intent(in) :: self
   integer, intent(inout) :: ncellsglobal, ncells, ncellssolve
   integer, intent(inout) :: nedgesglobal, nedges, nedgessolve
   integer, intent(inout) :: nvertlevels
   integer, intent(inout) :: nvertlevelsp1
 
   ncellsglobal  = self%nCellsGlobal
   ncells        = self%nCells
   ncellssolve   = self%nCellsSolve
 
   nedgesglobal  = self%nEdgesGlobal
   nedges        = self%nEdges
   nedgessolve   = self%nEdgesSolve
 
   nvertlevels   = self%nVertLevels
   nvertlevelsp1 = self%nVertLevelsP1
 
end subroutine geo_info
 
! ------------------------------------------------------------------------------
 
function field_is_templated(self, fieldname) result(is_templated)
   class(mpas_geom), intent(in) :: self
   character(len=*), intent(in) :: fieldname
   integer :: ii
   logical :: is_templated
   is_templated = .false.
   if (allocated(self%templated_fields)) then
      do ii = 1, size(self%templated_fields)
         if (trim(self%templated_fields(ii)%name) == trim(fieldname)) then
            is_templated = .true.
            return
         end if
      end do
   end if
end function field_is_templated
 
! ------------------------------------------------------------------------------
 
function template_fieldname(self, fieldname) result(template)
   class(mpas_geom), intent(in) :: self
   character(len=*), intent(in) :: fieldname
   integer :: ii
   character(len=MAXVARLEN) :: template
   do ii = 1, size(self%templated_fields)
      if (trim(self%templated_fields(ii)%name) == trim(fieldname)) then
         template = self%templated_fields(ii)%template
         return
      end if
   end do
   write(message,*) '--> mpas_geom % template_fieldname: fieldname is not templated, ',fieldname
   call abor1_ftn(message)
end function template_fieldname
 
! ------------------------------------------------------------------------------
 
function field_has_identity(self, fieldname) result(has_identity)
   class(mpas_geom), intent(in) :: self
   character(len=*), intent(in) :: fieldname
   integer :: ii
   logical :: has_identity
   has_identity = .false.
   if (allocated(self%templated_fields)) then
      do ii = 1, size(self%templated_fields)
         if (trim(self%templated_fields(ii)%name) == trim(fieldname) .and. &
             trim(self%templated_fields(ii)%identity) /= 'none') then
            has_identity = .true.
            return
         end if
      end do
   end if
end function field_has_identity
 
! ------------------------------------------------------------------------------
 
function identity_fieldname(self, fieldname) result(identity)
   class(mpas_geom), intent(in) :: self
   character(len=*), intent(in) :: fieldname
   integer :: ii
   character(len=MAXVARLEN) :: identity
   do ii = 1, size(self%templated_fields)
      if (trim(self%templated_fields(ii)%name) == trim(fieldname)) then
         identity = self%templated_fields(ii)%identity
         return
      end if
   end do
   write(message,*) '--> mpas_geom % identity_fieldname: fieldname does not have identity, ',fieldname
   call abor1_ftn(message)
end function identity_fieldname
 
! ------------------------------------------------------------------------------
 
!***********************************************************************
!
!  subroutine pool_has_field
!
!> \brief Check for presence of field in pool
!
!-----------------------------------------------------------------------
function pool_has_field(pool, field) result(has)
 
   implicit none
 
   type(mpas_pool_type), pointer, intent(in) :: pool
   character(len=*), intent(in) :: field
   type(mpas_pool_iterator_type) :: poolitr
   logical :: has
 
   has = .false.
 
   call mpas_pool_begin_iteration(pool)
 
   do while ( mpas_pool_get_next_member(pool, poolitr) )
      if (poolitr % memberType == mpas_pool_field .and. &
          trim(field) == trim(poolitr % memberName) ) then
         has = .true.
         exit
      endif
   end do
 
end function pool_has_field
 
 
!*******************************************************************************
!
! function isDecomposed
!
!> \brief whether horizontal dimension name is decomposed across MPI tasks
!
!*******************************************************************************
 
logical function isdecomposed(dimName)
   implicit none
   character(len=*), intent(in) :: dimname
   if (trim(dimname) == 'nCells') then
      isdecomposed = .true.
      return
   else if (trim(dimname) == 'nEdges') then
      isdecomposed = .true.
      return
   else if (trim(dimname) == 'nVertices') then
      isdecomposed = .true.
      return
   else
      isdecomposed = .false.
      return
   end if
end function isdecomposed
 
 
!***********************************************************************
!
!  function getSolveDimNames
!
!> \brief   Returns an array with the dimension names of a pool field member
!
!-----------------------------------------------------------------------
 
function getsolvedimnames(pool, key) result(solveDimNames)
 
   implicit none
 
   ! Arguments
   type(mpas_pool_type), pointer, intent(in) :: pool
   character(len=*), intent(in) :: key
 
   character (len=StrKIND), allocatable :: solvedimnames(:)
 
   type(mpas_pool_data_type), pointer :: dptr
   character (len=StrKIND) :: dimname
   integer :: id
 
   dptr => pool_get_member(pool, key, mpas_pool_field)
 
   if (.not. associated(dptr)) then
      write(message,*)'--> getSolveDimNames: problem finding ',trim(key), &
                  ' MPAS_POOL_FIELD in pool'
      call abor1_ftn(message)
   end if
 
   allocate(solvedimnames(dptr%contentsDims))
 
   do id = 1, dptr%contentsDims
 
      ! fields with only one time level
      !real
      if (associated(dptr%r1)) then
         dimname = trim(dptr%r1%dimNames(id))
      else if (associated(dptr%r2)) then
         dimname = trim(dptr%r2%dimNames(id))
      else if (associated(dptr%r3)) then
         dimname = trim(dptr%r3%dimNames(id))
      !integer
      else if (associated(dptr%i1)) then
         dimname = trim(dptr%i1%dimNames(id))
      else if (associated(dptr%i2)) then
         dimname = trim(dptr%i2%dimNames(id))
      else if (associated(dptr%i3)) then
         dimname = trim(dptr%i3%dimNames(id))
 
      ! fields with multiple time levels (use first time level)
      !real
      else if (associated(dptr%r1a)) then
         dimname = trim(dptr%r1a(1)%dimNames(id))
      else if (associated(dptr%r2a)) then
         dimname = trim(dptr%r2a(1)%dimNames(id))
      else if (associated(dptr%r3a)) then
         dimname = trim(dptr%r3a(1)%dimNames(id))
      !integer
      else if (associated(dptr%i1a)) then
         dimname = trim(dptr%i1a(1)%dimNames(id))
      else if (associated(dptr%i2a)) then
         dimname = trim(dptr%i2a(1)%dimNames(id))
      else if (associated(dptr%i3a)) then
         dimname = trim(dptr%i3a(1)%dimNames(id))
 
      else
         write(message,*) '--> getSolveDimNames: unsupported field type for key,&
            & contentsDims, contentsType = ',trim(key),',',dptr%contentsDims, &
            ',',dptr%contentsType
         call abor1_ftn(message)
      end if
 
      if (isdecomposed(dimname)) then
         solvedimnames(id) = trim(dimname)//'Solve'
      else
         solvedimnames(id) = trim(dimname)
      end if
   end do
end function getsolvedimnames
 
 
!***********************************************************************
!
!  function getSolveDimSizes
!
!> \brief   Returns an array with the dimension sizes of a pool field member
!
!> \details
!    In some cases, the pool containing the field will not contain any
!    dimensions.  Then the optional dimPool argument can be used.
!
!-----------------------------------------------------------------------
 
function getsolvedimsizes(pool, key, dimPool) result(solveDimSizes)
 
   implicit none
 
   ! Arguments
   type (mpas_pool_type), pointer, intent(in) :: pool
   character(len=*), intent(in) :: key
   type (mpas_pool_type), pointer, optional, intent(in) :: dimpool
 
   integer, allocatable :: solvedimsizes(:)
 
   character (len=StrKIND), allocatable :: dimnames(:)
   integer :: id
   integer, pointer :: dimsize
 
   dimnames = getsolvedimnames(pool, key)
 
   allocate(solvedimsizes(size(dimnames)))
 
   do id = 1, size(dimnames)
      if (present(dimpool)) then
         call mpas_pool_get_dimension(dimpool, dimnames(id), dimsize)
      else
         call mpas_pool_get_dimension(pool, dimnames(id), dimsize)
      end if
 
      if (associated(dimsize)) then
         solvedimsizes(id) = dimsize
      else
         write(message,*)'--> getSolveDimSizes: ',trim(dimnames(id)), &
                     ' not available'
         call abor1_ftn(message)
      end if
   end do
 
   deallocate(dimnames)
 
end function getsolvedimsizes
 
 
!***********************************************************************
!
!  function getSolveDimSizes
!
!> \brief   Returns the number of vertical levels for a pool field member
!
!> \details
!    In some cases, the pool containing the field will not contain any
!    dimensions.  Then the optional dimPool argument can be used.
!
!-----------------------------------------------------------------------
 
function getvertlevels(pool, key, dimPool) result(nlevels)
 
   ! Arguments
   type (mpas_pool_type), pointer, intent(in) :: pool
   character(len=*), intent(in) :: key
   type (mpas_pool_type), pointer, optional, intent(in) :: dimpool
 
   integer :: nlevels, id
   integer, allocatable :: dimsizes(:)
   character (len=StrKIND), allocatable :: dimnames(:)
 
   ! assume default value of 1, unless stated otherwise in field attributes
   nlevels = 1
 
   dimnames = getsolvedimnames(pool, key)
   if (present(dimpool)) then
      dimsizes = getsolvedimsizes(pool, key, dimpool)
   else
      dimsizes = getsolvedimsizes(pool, key)
   end if
   do id = 1, size(dimsizes)
      ! check if dimNames(id) is one of the available vertical coordinates
      if (ufo_vars_getindex(mpasverticalcoordinates, dimnames(id)) > 0) then
         nlevels = dimsizes(id)
      end if
   end do
 
   deallocate(dimnames, dimsizes)
 
end function getvertlevels
 
! ------------------------------------------------------------------------------
 
subroutine var_nlevels(self, var, nlevels)
 
   implicit none
 
   class(mpas_geom), intent(in) :: self
   character(len=*), intent(in) :: var
   integer, intent(out) :: nlevels
 
   character(len=StrKIND) :: poolVar
   type(mpas_pool_type), pointer :: allFields
 
   allfields => self % domain % blocklist % allFields
 
   if (pool_has_field(allfields, var)) then
      poolvar = var
   else if (self % is_templated(var)) then
      poolvar = self%template(var)
      if (.not.pool_has_field(allfields, poolvar)) then
         write(message,*)'--> vars_nlevels: ',trim(poolvar), &
                     ' not available in MPAS domain'
         call abor1_ftn(message)
      end if
   else
      write(message,*)'--> vars_nlevels: ',trim(var), &
                  ' not available in MPAS domain or self % templated_fields'
      call abor1_ftn(message)
   end if
 
   nlevels = getvertlevels(allfields, poolvar, self % domain % blocklist % dimensions)
 
end subroutine var_nlevels
 
! ------------------------------------------------------------------------------
 
subroutine variables_nlevels(self, vars, nv, nlevels)
 
   implicit none
 
   class(mpas_geom), intent(in) :: self
   type(oops_variables), intent(in) :: vars
   integer(c_size_t), intent(in) :: nv
   integer(c_size_t), intent(inout) :: nlevels(nv)
 
   integer :: ivar, nn
 
   do ivar = 1, nv
      call self%nlevels(vars%variable(ivar), nn)
      nlevels(ivar) = int(nn, c_size_t)
   end do
 
end subroutine variables_nlevels
 
! ------------------------------------------------------------------------------
 
end module mpas_geom_mod