mpas4da_mod.F90

Go to the documentation of this file.

! Copyright (c) 2018, National Atmospheric for Atmospheric Research (NCAR).
!
! Unless noted otherwise source code is licensed under the BSD license.
! Additional copyright and license information can be found in the LICENSE file
! distributed with this code, or at http://mpas-dev.github.com/license.html
 
module mpas4da_mod
 
   !***********************************************************************
   !
   !  Module mpas4da_mod to encapsulate operations needed for
   !  Data assimilation purpose.
   !  It can be used from /somewhere/MPAS/src/operators
   !  or from /somewhere/mpas-bundle/mpas/model (OOPS)
   !> \author  Gael Descombes/Mickael Duda NCAR/MMM
   !> \date    January 2018
   !
   !-----------------------------------------------------------------------
 
!fckit
use fckit_log_module, only: fckit_log
 
!oops
use kinds, only: kind_real
 
!saber?
use random_mod, only: normal_distribution
 
!ufo
use ufo_vars_mod
 
!MPAS-Model
use mpas_abort, only: mpas_dmpar_global_abort
use mpas_constants
use mpas_derived_types
use mpas_dmpar
use mpas_field_routines
use mpas_kind_types, only: shortstrkind
use mpas_pool_routines
 
!mpas-jedi
use mpas_constants_mod
use mpas_geom_mod, only: mpas_geom, pool_has_field, getsolvedimsizes
 
private
 
public :: &
   da_operator_addition, &
   da_copy_all2sub_fields, &
   da_copy_sub2all_fields, &
   da_template_pool, &
   !mpas_pool_template_field, &
   da_random, &
   da_operator, &
   da_self_mult, &
   da_constant, &
   da_posdef, &
   da_setval, &
   da_axpy, &
   da_gpnorm, &
   da_fldrms, &
   da_dot_product, &
   cvt_oopsmpas_date, &
   uv_cell_to_edges, &
   r3_normalize
 
character(len=1024) :: message
 
contains
 
   !***********************************************************************
   !
   !  function field_is_scalar
   !
   !> \brief   Test for a form of water vapor or hydrometeor of interest
   !> \details
   !>  At various places in this module we wish to test for the case where
   !>  a string is one of several 'q?' values. Rather than repeat that
   !>  logic many times, it is encapsulated here.
   !
   !-----------------------------------------------------------------------
   pure function field_is_scalar(fieldName)
 
      implicit none
 
      character (len=*), intent(in) :: fieldname
      logical :: field_is_scalar
      field_is_scalar = any(trim(fieldname) == &
                    (/'qv', 'qc', 'qi', 'qr', 'qs', 'qg', 'qh', 'nr'/))
 
   end function
 
   !***********************************************************************
   !
   !  function match_scalar
   !
   !> \brief   Test for a form of water vapor or hydrometeor of interest
   !> \author  Steven Vahl
   !> \date    11 July 2019
   !> \details
   !>  Test for the case where one string is 'scalars' and another string
   !>   matches with available scalar variables.
   !
   !-----------------------------------------------------------------------
   pure function match_scalar(scalarName, fieldName)
 
      implicit none
 
      character (len=*), intent(in) :: scalarname
      character (len=*), intent(in) :: fieldname
      logical :: match_scalar
 
      match_scalar = ( &
         scalarname == 'scalars' .and. &
         field_is_scalar(fieldname) )
 
   end function
 
   !***********************************************************************
   !
   !  subroutine mpas_pool_demo
   !
   !> \brief   Demonstrate basic usage of MPAS pools
   !> \author  Michael Duda
   !> \date    20 December 2017
   !> \details
   !>  This routine provides a simple demonstration of how to construct a new
   !>  pool at runtime, add members (fields) to the pool, and to perform generic
   !>  operations on that pool.
   !
   !-----------------------------------------------------------------------
   subroutine mpas_pool_demo(block)
 
      implicit none
 
      type (block_type), pointer :: block
 
      type (mpas_pool_type), pointer :: structs
      type (mpas_pool_type), pointer :: allFields
      type (mpas_pool_type), pointer :: da_state
      type (mpas_pool_type), pointer :: da_state_incr
 
      type (field2DReal), pointer :: field
 
      write(0,*) '****** Begin pool demo routine ******'
 
      structs => block % structs
      allfields => block % allFields
 
      !
      ! Create a new pool
      !
      call mpas_pool_create_pool(da_state)
 
      !
      ! Get pointers to several fields from the allFields pool, and add
      ! those fields to the da_state pool as well
      !
      call mpas_pool_get_field(allfields, 'theta', field)
      call mpas_pool_add_field(da_state, 'theta', field)
      write(0,*) 'Now, max value of theta is ', maxval(field % array),minval(field % array)
      field % array(:,:) = mpas_jedi_one_kr
      write(0,*)'Dimensions Field: ',field % dimSizes(:)
 
      call mpas_pool_get_field(allfields, 'rho', field)
      call mpas_pool_add_field(da_state, 'rho', field)
      write(0,*) 'Now, max value of rho is ', maxval(field % array),minval(field % array)
      field % array(:,:) = mpas_jedi_one_kr
 
      !
      ! Create another pool
      !
      call mpas_pool_create_pool(da_state_incr)
 
      !
      ! Duplicate the members of da_state into da_state_incr, and do a deep
      ! copy of the fields from da_state to da_state_incr
      !
      call mpas_pool_clone_pool(da_state, da_state_incr)
 
      !
      ! Call example algebra routine to compute A = A + B for all fields in
      ! the da_state and da_state_inc pools
      !
      call da_operator_addition(da_state, da_state_incr)
 
      call mpas_pool_get_field(da_state_incr, 'rho', field)
      write(0,*) 'Now, max value of rho_incr is ', maxval(field % array)
 
      call mpas_pool_get_field(da_state, 'rho', field)
      write(0,*) 'Now, max value of rho is ', maxval(field % array)
 
      !
      ! Before destroying a pool, we should remove any fields that are
      ! still referenced by other active pools to avoid deallocating them
      !
      call mpas_pool_empty_pool(da_state)
 
      !
      ! Destroy the now-empty da_state pool
      !
      call mpas_pool_destroy_pool(da_state)
 
      !
      ! Destroy the da_state_incr pool, deallocating all of its
      ! fields in the process (because this pool was not emptied)
      !
      call mpas_pool_destroy_pool(da_state_incr)
 
      write(0,*) '****** End pool demo routine ******'
 
   end subroutine mpas_pool_demo
 
   !***********************************************************************
   !
   !  subroutine da_operator_addition
   !
   !> \brief   Performs A = A + B for pools A and B
   !> \author  Michael Duda
   !> \date    20 December 2017
   !> \details
   !>  Given two pools, A and B, where the fields in B are a subset of
   !>  the fields in A, this routine adds the fields in B to fields in A
   !>  with the same name. When A and B contain identical fields, this
   !>  is equivalent to A = A + B.
   !
   !-----------------------------------------------------------------------
   subroutine da_operator_addition(pool_a, pool_b)
 
      implicit none
 
      type (mpas_pool_type), pointer :: pool_a, pool_b
 
      type (mpas_pool_iterator_type) :: poolitr
      real (kind=kind_real), pointer :: r0d_ptr_a, r0d_ptr_b
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a, r1d_ptr_b
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a, r2d_ptr_b
      real (kind=kind_real), dimension(:,:,:), pointer :: r3d_ptr_a, r3d_ptr_b
 
      !
      ! Iterate over all fields in pool_b, adding them to fields of the same
      ! name in pool_a
      !
      call mpas_pool_begin_iteration(pool_b)
 
      do while ( mpas_pool_get_next_member(pool_b, poolitr) )
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr % dataType == mpas_pool_real) then
 
               ! Depending on the dimensionality of the field, we need to set pointers of
               ! the correct type
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r0d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r0d_ptr_b)
                  r0d_ptr_a = r0d_ptr_a + r0d_ptr_b
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r1d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r1d_ptr_b)
                  r1d_ptr_a = r1d_ptr_a + r1d_ptr_b
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r2d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r2d_ptr_b)
                  r2d_ptr_a = r2d_ptr_a + r2d_ptr_b
                  write(message,*) 'Operator add MIN/MAX: ',minval(r2d_ptr_a),maxval(r2d_ptr_a)
                  call fckit_log%debug(message)
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r3d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r3d_ptr_b)
                  r3d_ptr_a = r3d_ptr_a + r3d_ptr_b
               end if
 
            end if
         end if
      end do
 
   end subroutine da_operator_addition
 
 
   !***********************************************************************
   !
   !  subroutine da_copy_all2sub_fields
   !
   !> \brief   Performs a copy of allfield to a sub pool A
   !> \author  Gael Desccombes
   !> \date    5 February 2018
   !> \details
   !>  Given two pools, allfields and A, where the fields in A are a subset of
   !>  the fields in allfields, this routine copy the fields allfields to fields in A
   !>  with the same name.
   !
   !-----------------------------------------------------------------------
   subroutine da_copy_all2sub_fields(domain, pool_a)
 
      implicit none
 
      type (mpas_pool_type), pointer, intent(inout) :: pool_a
      type (mpas_pool_type), pointer :: pool_b, state
      type (domain_type), pointer, intent(in) :: domain
 
      type (mpas_pool_iterator_type) :: poolitr_a, poolitr_b
      real (kind=kind_real), pointer :: r0d_ptr_a, r0d_ptr_b
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a, r1d_ptr_b
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a, r2d_ptr_b
      integer, pointer :: index_scalar
 
      type (field2dreal), pointer :: field2d
      type (field3dreal), pointer :: field3d
 
 
      pool_b => domain % blocklist % allFields
      call mpas_pool_get_subpool(domain % blocklist % structs,'state',state)
      !
      ! Iterate over all fields in pool_b, adding them to fields of the same
      ! name in pool_a
      !
      call mpas_pool_begin_iteration(pool_b)
 
      do while ( mpas_pool_get_next_member(pool_b, poolitr_b) )
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr_b % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr_b % dataType == mpas_pool_real) then
 
             call mpas_pool_begin_iteration(pool_a)
             do while ( mpas_pool_get_next_member(pool_a, poolitr_a) )
 
               if (( trim(poolitr_b % memberName)).eq.(trim(poolitr_a % memberName)) ) then
 
                  ! Depending on the dimensionality of the field, we need to set pointers of
                  ! the correct type
                  if (poolitr_b % nDims == 0) then
                     call mpas_pool_get_array(pool_a, trim(poolitr_a % memberName), r0d_ptr_a)
                     call mpas_pool_get_array(pool_b, trim(poolitr_b % memberName), r0d_ptr_b)
                     r0d_ptr_a = r0d_ptr_b
                  else if (poolitr_b % nDims == 1) then
                     call mpas_pool_get_array(pool_a, trim(poolitr_a % memberName), r1d_ptr_a)
                     call mpas_pool_get_array(pool_b, trim(poolitr_b % memberName), r1d_ptr_b)
                     r1d_ptr_a = r1d_ptr_b
                  else if (poolitr_b % nDims == 2) then
                     write(message,*) 'poolItr_b % memberName=',trim(poolitr_b % memberName)
                     call fckit_log%debug(message)
                     call mpas_pool_get_array(pool_a, trim(poolitr_a % memberName), r2d_ptr_a)
                     call mpas_pool_get_array(pool_b, trim(poolitr_b % memberName), r2d_ptr_b)
                     r2d_ptr_a = r2d_ptr_b
                     write(message,*) 'Copy all2sub field MIN/MAX: ',trim(poolitr_b % memberName), &
                                      minval(r2d_ptr_a),maxval(r2d_ptr_a)
                     call fckit_log%debug(message)
                  end if
 
               else if ( match_scalar(trim(poolitr_b % memberName), trim(poolitr_a % memberName)) ) then
                  write(message,*) 'Copy all2sub field: Looking at SCALARS now',trim(poolitr_a % memberName)
                  call fckit_log%debug(message)
                  call mpas_pool_get_dimension(state, 'index_'//trim(poolitr_a % memberName), index_scalar)
                  if (index_scalar .gt. 0) then
                     call mpas_pool_get_field(pool_a, trim(poolitr_a % memberName), field2d)
                     call mpas_pool_get_field(pool_b, trim(poolitr_b % memberName), field3d)
                     field2d % array(:,:) = field3d % array(index_scalar,:,:)
                     write(message,*) 'Copy all2sub field MIN/MAX: ',trim(poolitr_a % memberName), &
                                      minval(field2d % array), maxval(field2d % array)
                     call fckit_log%debug(message)
                  else
                     write(message,*) 'WARNING in Copy all2sub field; ',trim(poolitr_a % memberName), &
                                      'not available from MPAS'
                     call fckit_log%debug(message)
                  end if
               end if
            end do
           end if
         end if
      end do
 
   end subroutine da_copy_all2sub_fields
 
 
   !***********************************************************************
   !
   !  subroutine da_copy_sub2all_fields
   !
   !> \brief   Performs a copy of a sub pool A to allfields
   !> \author  Gael Desccombes
   !> \date    5 February 2018
   !> \details
   !>  Given two pools, allfields and A, where the fields in A are a subset of
   !>  the fields in allfields, this routine copy the subfields to allfields
   !>  with the same name.
   !
   !-----------------------------------------------------------------------
   subroutine da_copy_sub2all_fields(domain, pool_a)
 
      implicit none
 
      type (mpas_pool_type), pointer, intent(in) :: pool_a
      type (mpas_pool_type), pointer :: pool_b, state
      type (domain_type), pointer, intent(inout) :: domain
 
      type (mpas_pool_iterator_type) :: poolitr_a, poolitr_b
      real (kind=kind_real), pointer :: r0d_ptr_a, r0d_ptr_b
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a, r1d_ptr_b
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a, r2d_ptr_b
      integer, pointer :: index_scalar
 
      type (field2dreal), pointer :: field2d
      type (field3dreal), pointer :: field3d
 
 
      pool_b => domain % blocklist % allFields
      call mpas_pool_get_subpool(domain % blocklist % structs,'state',state)
      !
      ! Iterate over all fields in pool_b, adding them to fields of the same
      ! name in pool_a
      !
      call mpas_pool_begin_iteration(pool_b)
 
      do while ( mpas_pool_get_next_member(pool_b, poolitr_b) )
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr_b % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr_b % dataType == mpas_pool_real) then
 
             call mpas_pool_begin_iteration(pool_a)
             do while ( mpas_pool_get_next_member(pool_a, poolitr_a) )
 
               if (( trim(poolitr_b % memberName)).eq.(trim(poolitr_a % memberName)) ) then
 
                  ! Depending on the dimensionality of the field, we need to set pointers of
                  ! the correct type
                  if (poolitr_b % nDims == 0) then
                     call mpas_pool_get_array(pool_a, trim(poolitr_a % memberName), r0d_ptr_a)
                     call mpas_pool_get_array(pool_b, trim(poolitr_b % memberName), r0d_ptr_b)
                     r0d_ptr_b = r0d_ptr_a
                  else if (poolitr_b % nDims == 1) then
                     call mpas_pool_get_array(pool_a, trim(poolitr_a % memberName), r1d_ptr_a)
                     call mpas_pool_get_array(pool_b, trim(poolitr_b % memberName), r1d_ptr_b)
                     r1d_ptr_b = r1d_ptr_a
                     write(message,*) 'Copy sub2all field MIN/MAX: ',trim(poolitr_b % memberName), &
                                      minval(r1d_ptr_a),maxval(r1d_ptr_a)
                     call fckit_log%debug(message)
                  else if (poolitr_b % nDims == 2) then
                     call mpas_pool_get_array(pool_a, trim(poolitr_a % memberName), r2d_ptr_a)
                     call mpas_pool_get_array(pool_b, trim(poolitr_b % memberName), r2d_ptr_b)
                     r2d_ptr_b = r2d_ptr_a
                     write(message,*) 'Copy sub2all field MIN/MAX: ',trim(poolitr_b % memberName), &
                                      minval(r2d_ptr_a),maxval(r2d_ptr_a)
                     call fckit_log%debug(message)
                  end if
 
               else if ( match_scalar(trim(poolitr_b % memberName), trim(poolitr_a % memberName)) ) then
                  write(message,*) 'Copy sub2all field: Looking at SCALARS now',trim(poolitr_a % memberName)
                  call fckit_log%debug(message)
                  call mpas_pool_get_dimension(state, 'index_'//trim(poolitr_a % memberName), index_scalar)
                  if (index_scalar .gt. 0) then
                     call mpas_pool_get_field(pool_a, trim(poolitr_a % memberName), field2d)
                     call mpas_pool_get_field(pool_b, trim(poolitr_b % memberName), field3d)
                     field3d % array(index_scalar,:,:) = field2d % array(:,:)
                     write(message,*) 'Copy sub2all field MIN/MAX: ',trim(poolitr_a % memberName), &
                                      minval(field2d % array), maxval(field2d % array)
                     call fckit_log%debug(message)
                  else
                     write(message,*) 'WARNING in Copy sub2all field; ',trim(poolitr_a % memberName), &
                                      'not available from MPAS'
                     call fckit_log%debug(message)
                  end if
!                       end if
               end if
            end do
           end if
         end if
      end do
 
   end subroutine da_copy_sub2all_fields
 
 
   !***********************************************************************
   !
   !  subroutine da_template_pool
   !
   !> \brief   Subset a pool from fields described in geom
   !> \details
   !>  Given an mpas_geom object, create templatePool containing a subset
   !>  of the fields in the geom % domain's allFields pool and the
   !>  geom % templated_fields
   !
   !-----------------------------------------------------------------------
   subroutine da_template_pool(geom, templatePool, nf, fieldnames)
 
      implicit none
 
      ! Arguments
      type (mpas_geom), intent(in) :: geom
      type (mpas_pool_type), pointer, intent(out) :: templatepool
      integer, intent(in) :: nf
      character (len=*), intent(in) :: fieldnames(nf)
 
      ! Local variables
      type (mpas_pool_type), pointer :: allfields, state
      character(len=MAXVARLEN) :: fieldname, template
      !type (field2DReal), pointer :: qField
      !type (field3DReal), pointer :: scalars
 
      integer, pointer :: index_scalar, dim0d
      integer :: ii
      integer, parameter :: ndims=10
      character(len=ShortStrKIND) :: dimnames(ndims)
 
      call mpas_pool_create_pool(templatepool)
 
      dimnames( 1) = 'nCellsSolve'
      dimnames( 2) = 'nEdgesSolve'
      dimnames( 3) = 'nVerticesSolve'
      dimnames( 4) = 'nVertLevels'
      dimnames( 5) = 'nVertLevelsP1'
      dimnames( 6) = 'nSoilLevels'
      dimnames( 7) = 'nCells'
      dimnames( 8) = 'nEdges'
      dimnames( 9) = 'nVertices'
      dimnames(10) = 'vertexDegree'
 
      do ii = 1, ndims
         call mpas_pool_get_dimension(geom % domain % blocklist % dimensions, trim(dimnames(ii)), dim0d)
         call mpas_pool_add_dimension(templatepool, trim(dimnames(ii)), dim0d)
      end do
 
      call mpas_pool_get_subpool(geom % domain % blocklist % structs, 'state',state)
      allfields => geom % domain % blocklist % allFields
 
      do ii=1, nf
         fieldname = fieldnames(ii)
 
         if (field_is_scalar(fieldname)) then
            ! Check if this scalar is activated
            call mpas_pool_get_dimension(state, 'index_'//trim(fieldname), index_scalar)
            if (index_scalar .gt. 0) then
               call mpas_pool_template_field('theta', allfields, fieldname, templatepool)
            else
               write(message,*)'--> da_template_pool: ',trim(fieldname), &
                           ' not available in MPAS domain'
               call abor1_ftn(message)
            end if
         else
            if (pool_has_field(allfields, fieldname)) then
               template = fieldname
            else if (geom % is_templated(fieldname)) then
               template = geom%template(fieldname)
               if (.not.pool_has_field(allfields, template)) then
                  write(message,*)'--> da_template_pool: ',trim(template), &
                              ' not available in MPAS domain'
                  call abor1_ftn(message)
               end if
            else
               write(message,*)'--> da_template_pool: ',trim(fieldname), &
                           ' not available in MPAS domain or geom % templated_fields'
               call abor1_ftn(message)
            end if
 
            call mpas_pool_template_field(template, allfields, fieldname, templatepool)
 
         end if
      end do
 
      call da_constant(templatepool, mpas_jedi_zero_kr)
 
   end subroutine da_template_pool
 
 
   !***********************************************************************
   !
   !  subroutine mpas_pool_template_field
   !
   !> \brief   Add a field to dstPool that is templated on a srcPool field
   !> \details
   !>  Duplicate srcFieldName from srcPool with the name dstFieldName
   !>  in dstPool
   !
   !-----------------------------------------------------------------------
   subroutine mpas_pool_template_field(srcFieldName, srcPool, dstFieldName, dstPool)
 
      ! Arguments
      type (mpas_pool_type), pointer, intent(in) :: srcPool
      character (len=*), intent(in) :: srcFieldName, dstFieldName
      type (mpas_pool_type), pointer, intent(inout) :: dstPool
 
      ! Local variables
      type (mpas_pool_iterator_type) :: poolItr
      type (field0DReal), pointer :: field0d_src, field0d_dst
      type (field1DReal), pointer :: field1d_src, field1d_dst
      type (field2DReal), pointer :: field2d_src, field2d_dst
      type (field3DReal), pointer :: field3d_src, field3d_dst, field3d
      type (field0DInteger), pointer :: ifield0d_src, ifield0d_dst
      type (field1DInteger), pointer :: ifield1d_src, ifield1d_dst
      type (field2DInteger), pointer :: ifield2d_src, ifield2d_dst
      type (field3DInteger), pointer :: ifield3d_src, ifield3d_dst
 
      !
      ! Iterate over srcPool, add the one that matches srcFieldName into dstPool
      !
      call mpas_pool_begin_iteration(srcpool)
 
      do while ( mpas_pool_get_next_member(srcpool, poolitr) )
 
         ! Only handle fields, ignore dimensions, namelist options, and other pools
         if ( poolitr % memberType == mpas_pool_field .and. &
              trim(srcfieldname) == trim(poolitr % memberName) ) then
 
            ! Handle real and integer, ignore logical and char
            if (poolitr % dataType == mpas_pool_real) then
 
               ! Use correctly dimensioned field pointers
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_field(srcpool, trim(srcfieldname), field0d_src)
                  call mpas_duplicate_field(field0d_src, field0d_dst)
                  field0d_dst % fieldName = dstfieldname
                  call mpas_pool_add_field(dstpool, trim(dstfieldname), field0d_dst)
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_field(srcpool, trim(srcfieldname), field1d_src)
                  call mpas_duplicate_field(field1d_src, field1d_dst)
                  field1d_dst % fieldName = dstfieldname
                  call mpas_pool_add_field(dstpool, trim(dstfieldname), field1d_dst)
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_field(srcpool, trim(srcfieldname), field2d_src)
                  call mpas_duplicate_field(field2d_src, field2d_dst)
                  field2d_dst % fieldName = dstfieldname
                  call mpas_pool_add_field(dstpool, trim(dstfieldname), field2d_dst)
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_field(srcpool, trim(srcfieldname), field3d_src)
                  call mpas_duplicate_field(field3d_src, field3d_dst)
                  field3d_dst % fieldName = dstfieldname
                  call mpas_pool_add_field(dstpool, trim(dstfieldname), field3d_dst)
               end if
 
            else if (poolitr % dataType == mpas_pool_integer) then
 
               ! Use correctly dimensioned field pointers
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_field(srcpool, trim(srcfieldname), ifield0d_src)
                  call mpas_duplicate_field(ifield0d_src, ifield0d_dst)
                  ifield0d_dst % fieldName = dstfieldname
                  call mpas_pool_add_field(dstpool, trim(dstfieldname), ifield0d_dst)
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_field(srcpool, trim(srcfieldname), ifield1d_src)
                  call mpas_duplicate_field(ifield1d_src, ifield1d_dst)
                  ifield1d_dst % fieldName = dstfieldname
                  call mpas_pool_add_field(dstpool, trim(dstfieldname), ifield1d_dst)
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_field(srcpool, trim(srcfieldname), ifield2d_src)
                  call mpas_duplicate_field(ifield2d_src, ifield2d_dst)
                  ifield2d_dst % fieldName = dstfieldname
                  call mpas_pool_add_field(dstpool, trim(dstfieldname), ifield2d_dst)
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_field(srcpool, trim(srcfieldname), ifield3d_src)
                  call mpas_duplicate_field(ifield3d_src, ifield3d_dst)
                  ifield3d_dst % fieldName = dstfieldname
                  call mpas_pool_add_field(dstpool, trim(dstfieldname), ifield3d_dst)
               end if
            end if
 
            return
 
         end if
      end do
 
      write(message,*)'--> mpas_pool_template_field: ',trim(srcfieldname), &
                  ' not available in srcPool'
      call abor1_ftn(message)
 
   end subroutine mpas_pool_template_field
 
 
   !***********************************************************************
   !
   !  function da_common_vars
   !
   !> \author  Gael Descombes
   !> \date    26 December 2017
   !> \details
   !>  Count the number of fields in a Pool related to a list of fields
   !
   !-----------------------------------------------------------------------
   function da_common_vars(pool_a, fieldname) result(nsize0)
 
      implicit none
      type (mpas_pool_type), pointer :: pool_a, pool_b
      type (mpas_pool_iterator_type) :: poolitr
      character (len=*) :: fieldname(:)
      integer :: ii, jj, nsize, nsize0
 
      nsize0 = 0
      nsize  = size(fieldname)
      call mpas_pool_begin_iteration(pool_a)
 
         do while ( mpas_pool_get_next_member(pool_a, poolitr) )
            ! Pools may in general contain dimensions, namelist options, fields, or other pools,
            ! so we select only those members of the pool that are fields
            if (poolitr % memberType == mpas_pool_field) then
               ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
               if (poolitr % dataType == mpas_pool_real) then
                  do ii=1, nsize
                     if ( trim(fieldname(ii)).eq.(trim(poolitr % memberName)) ) then
                        nsize0 = nsize0 + 1
                     else if (match_scalar(trim(poolitr % memberName), trim(fieldname(ii)))) then
                        nsize0 = nsize0 + 1
                     end if
                  end do
               end if
            end if
         end do
 
   end function da_common_vars
 
 
 
   !***********************************************************************
   !
   !  subroutine da_random
   !
   !> \brief   Performs random for pool A
   !> \author  Gael Descombes
   !> \date    January 2018
   !> \details
   !
   !-----------------------------------------------------------------------
   subroutine da_random(pool_a, fld_select)
 
      implicit none
 
      type (mpas_pool_type), pointer, intent(inout) :: pool_a
      character (len=*), optional,    intent(in)    :: fld_select(:)
 
      integer, parameter :: rseed = 7
      type (mpas_pool_iterator_type) :: poolitr
      real (kind=kind_real), pointer :: r0d_ptr_a
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a
      real (kind=kind_real), dimension(:,:,:), pointer :: r3d_ptr_a
 
      !
      ! Iterate over all fields in pool_b, adding them to fields of the same
      ! name in pool_a
      !
      call mpas_pool_begin_iteration(pool_a)
 
      do while ( mpas_pool_get_next_member(pool_a, poolitr) )
 
         if (present(fld_select)) then
            if (ufo_vars_getindex(fld_select,trim(poolitr % memberName)) < 0) cycle
         end if
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr % dataType == mpas_pool_real) then
 
               ! Depending on the dimensionality of the field, we need to set pointers of
               ! the correct type
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r0d_ptr_a)
                  !call normal_distribution(r0d_ptr_a, MPAS_JEDI_ZERO_kr, MPAS_JEDI_ONE_kr, rseed)
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r1d_ptr_a)
                  call normal_distribution(r1d_ptr_a, mpas_jedi_zero_kr, mpas_jedi_one_kr, rseed)
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r2d_ptr_a)
                  call normal_distribution(r2d_ptr_a, mpas_jedi_zero_kr, mpas_jedi_one_kr, rseed)
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r3d_ptr_a)
                  call normal_distribution(r3d_ptr_a, mpas_jedi_zero_kr, mpas_jedi_one_kr, rseed)
               end if
 
            end if
         end if
      end do
 
   end subroutine da_random
 
   !-----------------------------------------------------------------------
   !  subroutine da_operator
   !
   !> \brief   Performs A = A 'kind_op' B for pools A and B
   !> \author  Michael Duda
   !> \date    20 December 2017
   !> \details
   !>  Given two pools, A and B, where the fields in B are a subset of
   !>  the fields in A, this routine adds the fields in B to fields in A
   !>  with the same name. When A and B contain identical fields, this
   !>  is equivalent to A = A 'kind_op' B.
   !>  \modified by Gael DESCOMBES to apply diffferent operator
   !
   !-----------------------------------------------------------------------
   subroutine da_operator(kind_op, pool_a, pool_b, pool_c, fld_select)
 
      implicit none
 
      type (mpas_pool_type), pointer :: pool_a, pool_b
      type (mpas_pool_type), pointer, optional :: pool_c
      character (len=*) :: kind_op
      character (len=*), optional :: fld_select(:)
 
      type (mpas_pool_iterator_type) :: poolitr
      real (kind=kind_real), pointer :: r0d_ptr_a, r0d_ptr_b, r0d_ptr_c
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a, r1d_ptr_b, r1d_ptr_c
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a, r2d_ptr_b, r2d_ptr_c
      real (kind=kind_real), dimension(:,:,:), pointer :: r3d_ptr_a, r3d_ptr_b, r3d_ptr_c
      !
      ! Iterate over all fields in pool_b, adding them to fields of the same
      ! name in pool_a
      !
      call mpas_pool_begin_iteration(pool_b)
 
      write(message,*) 'Operator ',trim(kind_op)
      call fckit_log%debug(message)
 
      do while ( mpas_pool_get_next_member(pool_b, poolitr) )
 
         if (present(fld_select)) then
            if (ufo_vars_getindex(fld_select,trim(poolitr % memberName)) < 0) cycle
         end if
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr % dataType == mpas_pool_real) then
 
               ! Depending on the dimensionality of the field, we need to set pointers of
               ! the correct type
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r0d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r0d_ptr_b)
                  if (present(pool_c)) then
                     call mpas_pool_get_array(pool_c, trim(poolitr % memberName), r0d_ptr_c)
                     r0d_ptr_a = mpas_jedi_zero_kr
                  end if
                  if ( trim(kind_op).eq.'add' ) then
                     r0d_ptr_a = r0d_ptr_a + r0d_ptr_b
                     if (present(pool_c)) then
                       r0d_ptr_a = r0d_ptr_b + r0d_ptr_c
                     else
                        r0d_ptr_a = r0d_ptr_a + r0d_ptr_b
                     end if
                  else if ( trim(kind_op).eq.'schur' ) then
                     if (present(pool_c)) then
                        r0d_ptr_a = r0d_ptr_b * r0d_ptr_c
                     else
                        r0d_ptr_a = r0d_ptr_a * r0d_ptr_b
                     end if
                  else if ( trim(kind_op).eq.'sub' ) then
                     if (present(pool_c)) then
                        r0d_ptr_a = r0d_ptr_b - r0d_ptr_c
                     else
                        r0d_ptr_a = r0d_ptr_a - r0d_ptr_b
                     end if
                  end if
 
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r1d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r1d_ptr_b)
                  if (present(pool_c)) then
                     call mpas_pool_get_array(pool_c, trim(poolitr % memberName), r1d_ptr_c)
                     r1d_ptr_a = mpas_jedi_zero_kr
                  end if
                  if ( trim(kind_op).eq.'add' ) then
                     if (present(pool_c)) then
                        r1d_ptr_a = r1d_ptr_b + r1d_ptr_c
                     else
                        r1d_ptr_a = r1d_ptr_a + r1d_ptr_b
                     end if
                  else if ( trim(kind_op).eq.'schur' ) then
                     if (present(pool_c)) then
                        r1d_ptr_a = r1d_ptr_b * r1d_ptr_c
                     else
                        r1d_ptr_a = r1d_ptr_a * r1d_ptr_b
                     end if
                  else if ( trim(kind_op).eq.'sub' ) then
                     if (present(pool_c)) then
                        r1d_ptr_a = r1d_ptr_b - r1d_ptr_c
                     else
                        r1d_ptr_a = r1d_ptr_a - r1d_ptr_b
                     end if
                  end if
 
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r2d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r2d_ptr_b)
                  if (present(pool_c)) then
                     call mpas_pool_get_array(pool_c, trim(poolitr % memberName), r2d_ptr_c)
                     r2d_ptr_a = mpas_jedi_zero_kr
                  end if
                  if ( trim(kind_op).eq.'add' ) then
                     write(message,*) 'Operator_a add MIN/MAX: ',minval(r2d_ptr_a),maxval(r2d_ptr_a)
                     call fckit_log%debug(message)
                     write(message,*) 'Operator_b add MIN/MAX: ',minval(r2d_ptr_b),maxval(r2d_ptr_b)
                     call fckit_log%debug(message)
                     if (present(pool_c)) then
                        r2d_ptr_a = r2d_ptr_b + r2d_ptr_c
                     else
                        r2d_ptr_a = r2d_ptr_a + r2d_ptr_b
                     end if
                     write(message,*) 'Operator2 add MIN/MAX: ',minval(r2d_ptr_a),maxval(r2d_ptr_a)
                     call fckit_log%debug(message)
                  else if ( trim(kind_op).eq.'schur' ) then
                     if (present(pool_c)) then
                        r2d_ptr_a = r2d_ptr_b * r2d_ptr_c
                     else
                        r2d_ptr_a = r2d_ptr_a * r2d_ptr_b
                     end if
                  else if ( trim(kind_op).eq.'sub' ) then
                     if (present(pool_c)) then
                        r2d_ptr_a = r2d_ptr_b - r2d_ptr_c
                     else
                        r2d_ptr_a = r2d_ptr_a - r2d_ptr_b
                     end if
                  end if
 
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r3d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r3d_ptr_b)
                  if (present(pool_c)) then
                     call mpas_pool_get_array(pool_c, trim(poolitr % memberName), r3d_ptr_c)
                     r3d_ptr_a = mpas_jedi_zero_kr
                  end if
                  if ( trim(kind_op).eq.'add' ) then
                     if (present(pool_c)) then
                        r3d_ptr_a = r3d_ptr_b + r3d_ptr_c
                     else
                        r3d_ptr_a = r3d_ptr_a + r3d_ptr_b
                     end if
                  else if ( trim(kind_op).eq.'schur' ) then
                     if (present(pool_c)) then
                        r3d_ptr_a = r3d_ptr_b * r3d_ptr_c
                     else
                        r3d_ptr_a = r3d_ptr_a * r3d_ptr_b
                     end if
                  else if ( trim(kind_op).eq.'sub' ) then
                     if (present(pool_c)) then
                        r3d_ptr_a = r3d_ptr_b - r3d_ptr_c
                     else
                        r3d_ptr_a = r3d_ptr_a - r3d_ptr_b
                     end if
                  end if
               end if
 
            end if
         end if
      end do
 
   end subroutine da_operator
 
   !***********************************************************************
   !
   !  subroutine da_self_mult
   !
   !> \brief   Performs A = A * zz for pool A, zz a real number
   !> \author  Gael Descombes
   !> \date    22 December 2017
   !> \details
   !
   !-----------------------------------------------------------------------
   subroutine da_self_mult(pool_a, zz)
 
      implicit none
 
      type (mpas_pool_type), pointer :: pool_a
      real (kind=kind_real) :: zz
 
      type (mpas_pool_iterator_type) :: poolitr
      real (kind=kind_real), pointer :: r0d_ptr_a
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a
      real (kind=kind_real), dimension(:,:,:), pointer :: r3d_ptr_a
 
      !
      ! Iterate over all fields in pool_b, adding them to fields of the same
      ! name in pool_a
      !
      call mpas_pool_begin_iteration(pool_a)
 
      do while ( mpas_pool_get_next_member(pool_a, poolitr) )
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr % dataType == mpas_pool_real) then
 
               ! Depending on the dimensionality of the field, we need to set pointers of
               ! the correct type
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r0d_ptr_a)
                  r0d_ptr_a = r0d_ptr_a * zz
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r1d_ptr_a)
                  r1d_ptr_a = r1d_ptr_a * zz
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r2d_ptr_a)
                  r2d_ptr_a = r2d_ptr_a * zz
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r3d_ptr_a)
                  r3d_ptr_a = r3d_ptr_a * zz
               end if
 
            end if
         end if
      end do
 
   end subroutine da_self_mult
 
 
   !***********************************************************************
   !
   !  subroutine da_constant
   !
   !> \brief   Performs A = constant. for pool A
   !> \author  Gael Descombes
   !> \date    22 December 2017
   !> \details
   !
   !-----------------------------------------------------------------------
   subroutine da_constant(pool_a, realvalue, fld_select)
 
      implicit none
 
      type (mpas_pool_type), pointer, intent(inout) :: pool_a
      real (kind=kind_real),          intent(in)    :: realvalue
      character (len=*), optional,    intent(in)    :: fld_select(:)
 
      type (mpas_pool_iterator_type) :: poolitr
      real (kind=kind_real), pointer :: r0d_ptr_a
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a
      real (kind=kind_real), dimension(:,:,:), pointer :: r3d_ptr_a
 
      call mpas_pool_begin_iteration(pool_a)
 
      do while ( mpas_pool_get_next_member(pool_a, poolitr) )
 
         if (present(fld_select)) then
            if (ufo_vars_getindex(fld_select,trim(poolitr % memberName)) < 0) cycle
         end if
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr % dataType == mpas_pool_real) then
 
               ! Depending on the dimensionality of the field, we need to set pointers of
               ! the correct type
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r0d_ptr_a)
                  r0d_ptr_a = realvalue
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r1d_ptr_a)
                  r1d_ptr_a = realvalue
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r2d_ptr_a)
                  r2d_ptr_a = realvalue
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r3d_ptr_a)
                  r3d_ptr_a = realvalue
               end if
 
            end if
         end if
      end do
 
   end subroutine da_constant
 
   !***********************************************************************
   !
   !  subroutine da_posdef
   !
   !> \brief   Performs A = max(0.,A) for pool A
   !> \author  JJ Guerrette
   !> \date    12 July 2019
   !> \details
   !
   !-----------------------------------------------------------------------
   subroutine da_posdef(pool_a, fld_select)
 
      implicit none
 
      type (mpas_pool_type), pointer, intent(inout) :: pool_a
      character (len=*), optional,    intent(in)    :: fld_select(:)
 
      type (mpas_pool_iterator_type) :: poolitr
      real (kind=kind_real), pointer :: r0d_ptr_a
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a
      real (kind=kind_real), dimension(:,:,:), pointer :: r3d_ptr_a
 
      !
      ! Iterate over all fields in pool_b, adding them to fields of the same
      ! name in pool_a
      !
      call mpas_pool_begin_iteration(pool_a)
 
      do while ( mpas_pool_get_next_member(pool_a, poolitr) )
 
         if (present(fld_select)) then
            if (ufo_vars_getindex(fld_select,trim(poolitr % memberName)) < 0) cycle
         end if
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr % dataType == mpas_pool_real) then
 
               ! Depending on the dimensionality of the field, we need to set pointers of
               ! the correct type
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r0d_ptr_a)
                  r0d_ptr_a = max(mpas_jedi_zero_kr, r0d_ptr_a)
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r1d_ptr_a)
                  r1d_ptr_a = max(mpas_jedi_zero_kr, r1d_ptr_a)
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r2d_ptr_a)
                  r2d_ptr_a = max(mpas_jedi_zero_kr, r2d_ptr_a)
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r3d_ptr_a)
                  r3d_ptr_a = max(mpas_jedi_zero_kr, r3d_ptr_a)
               end if
 
            end if
         end if
      end do
 
   end subroutine da_posdef
 
   !***********************************************************************
   !
   !  subroutine da_setval
   !
   !> \brief   Performs A = Val_R. for pool A
   !> \author  Gael Descombes
   !> \date    22 December 2017
   !> \details
   !
   !-----------------------------------------------------------------------
   subroutine da_setval(pool_a,zz)
 
      implicit none
 
      type (mpas_pool_type), pointer :: pool_a
      real (kind=kind_real) :: zz
 
      type (mpas_pool_iterator_type) :: poolitr
      real (kind=kind_real), pointer :: r0d_ptr_a
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a
      real (kind=kind_real), dimension(:,:,:), pointer :: r3d_ptr_a
 
      !
      ! Iterate over all fields in pool_b, adding them to fields of the same
      ! name in pool_a
      !
      call mpas_pool_begin_iteration(pool_a)
 
      do while ( mpas_pool_get_next_member(pool_a, poolitr) )
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr % dataType == mpas_pool_real) then
 
               ! Depending on the dimensionality of the field, we need to set pointers of
               ! the correct type
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r0d_ptr_a)
                  r0d_ptr_a = zz
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r1d_ptr_a)
                  r1d_ptr_a = zz
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r2d_ptr_a)
                  r2d_ptr_a = zz
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r3d_ptr_a)
                  r3d_ptr_a = zz
               end if
 
            end if
         end if
      end do
 
   end subroutine da_setval
 
 
 
 
 
   !***********************************************************************
   !
   !  subroutine da_axpy
   !
   !> \brief   Performs A = A + B * zz for pools A and B
   !> \author  Gael Descombes
   !> \date    20 December 2017
   !> \details
   !>  Given two pools, A and B, where the fields in B are a subset of
   !>  the fields in A, this routine adds the fields in B to fields in A
   !>  with the same name. When A and B contain identical fields, this
   !>  is equivalent to A = A + B.
   !
   !-----------------------------------------------------------------------
   subroutine da_axpy(pool_a, pool_b, zz, fld_select)
 
      implicit none
      type (mpas_pool_type), pointer, intent(inout) :: pool_a
      type (mpas_pool_type), pointer, intent(in)    :: pool_b
      real (kind=kind_real), intent(in) :: zz
      character (len=*), optional, intent(in) :: fld_select(:)
 
 
 
      type (mpas_pool_iterator_type) :: poolitr
      real (kind=kind_real), pointer :: r0d_ptr_a, r0d_ptr_b
      real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a, r1d_ptr_b
      real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a, r2d_ptr_b
      real (kind=kind_real), dimension(:,:,:), pointer :: r3d_ptr_a, r3d_ptr_b
 
      !
      ! Iterate over all fields in pool_b, adding them to fields of the same
      ! name in pool_a
      !
      call mpas_pool_begin_iteration(pool_b)
 
      do while ( mpas_pool_get_next_member(pool_b, poolitr) )
 
         if (present(fld_select)) then
            if (ufo_vars_getindex(fld_select,trim(poolitr % memberName)) < 0) cycle
         end if
 
         ! Pools may in general contain dimensions, namelist options, fields, or other pools,
         ! so we select only those members of the pool that are fields
         if (poolitr % memberType == mpas_pool_field) then
 
            ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
            if (poolitr % dataType == mpas_pool_real) then
 
               ! Depending on the dimensionality of the field, we need to set pointers of
               ! the correct type
               if (poolitr % nDims == 0) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r0d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r0d_ptr_b)
                  r0d_ptr_a = r0d_ptr_a + r0d_ptr_b * zz
               else if (poolitr % nDims == 1) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r1d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r1d_ptr_b)
                  r1d_ptr_a = r1d_ptr_a + r1d_ptr_b * zz
               else if (poolitr % nDims == 2) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r2d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r2d_ptr_b)
                  r2d_ptr_a = r2d_ptr_a + r2d_ptr_b * zz
               else if (poolitr % nDims == 3) then
                  call mpas_pool_get_array(pool_a, trim(poolitr % memberName), r3d_ptr_a)
                  call mpas_pool_get_array(pool_b, trim(poolitr % memberName), r3d_ptr_b)
                  r3d_ptr_a = r3d_ptr_a + r3d_ptr_b * zz
               end if
 
            end if
         end if
      end do
 
   end subroutine da_axpy
 
 
   !***********************************************************************
   !
   !  subroutine da_gpnorm
   !
   !> \brief   Performs basic statistics min/max/norm given a pool
   !> \author  Gael Descombes
   !> \date    February 2018
   !> \details
   !>  Given a pool of fields, return min/max/norm array
   !
   !-----------------------------------------------------------------------
 
   subroutine da_gpnorm(pool_a, dminfo, nf, pstat, fld_select)
 
   implicit none
   type (mpas_pool_type), pointer, intent(in)  :: pool_a
   type (dm_info), pointer,        intent(in)  :: dminfo
   integer,                        intent(in)  :: nf
   character (len=*),              intent(in)  :: fld_select(nf)
   real(kind=kind_real),           intent(out) :: pstat(3, nf)
 
   type (mpas_pool_iterator_type) :: poolitr
   type (field1dreal), pointer :: field1d
   type (field2dreal), pointer :: field2d
   type (field3dreal), pointer :: field3d
   real(kind=kind_real) :: globalsum, globalmin, globalmax, dimtot, dimtot_global, prodtot
 
   integer :: jj, ndims
   integer :: dim1, dim2, dim3
   integer, allocatable :: dimsizes(:)
 
   pstat = mpas_jedi_zero_kr
 
   !
   ! Iterate over all fields in pool_a
   ! name in pool_a
   !
   call mpas_pool_begin_iteration(pool_a)
 
   do while ( mpas_pool_get_next_member(pool_a, poolitr) )
      jj = ufo_vars_getindex(fld_select,trim(poolitr % memberName))
      if ( jj < 0 .or. jj > nf ) cycle
 
      ! Pools may in general contain dimensions, namelist options, fields, or other pools,
      ! so we select only those members of the pool that are fields
      if (poolitr % memberType == mpas_pool_field) then
 
         ! Fields can be integer, logical, or real. Here, we operate only on real-valued fields
         if (poolitr % dataType == mpas_pool_real) then
 
            ! Depending on the dimensionality of the field, we need to set pointers of
            ! the correct type
            ndims = poolitr % nDims
            dimsizes = getsolvedimsizes(pool_a, poolitr%memberName)
            if (ndims == 1) then
               dim1 = dimsizes(1)
               call mpas_pool_get_field(pool_a, trim(poolitr % memberName), field1d)
               dimtot = real(dim1,kind_real)
               prodtot = sum(field1d % array(1:dim1)**2 )
               call mpas_dmpar_sum_real(dminfo, dimtot, dimtot_global)
               call mpas_dmpar_sum_real(dminfo, prodtot, globalsum)
               call mpas_dmpar_min_real(dminfo, minval(field1d % array(1:dim1)), globalmin)
               call mpas_dmpar_max_real(dminfo, maxval(field1d % array(1:dim1)), globalmax)
               pstat(1,jj) = globalmin
               pstat(2,jj) = globalmax
               pstat(3,jj) = sqrt( globalsum / dimtot_global )
            else if (ndims == 2) then
               dim1 = dimsizes(1)
               dim2 = dimsizes(2)
               call mpas_pool_get_field(pool_a, trim(poolitr % memberName), field2d)
               dimtot  = real(dim1*dim2,kind_real)
               prodtot = sum(field2d % array(1:dim1,1:dim2)**2 )
               call mpas_dmpar_sum_real(dminfo, dimtot, dimtot_global)
               call mpas_dmpar_sum_real(dminfo, prodtot, globalsum)
               call mpas_dmpar_min_real(dminfo, minval(field2d % array(1:dim1,1:dim2)), globalmin)
               call mpas_dmpar_max_real(dminfo, maxval(field2d % array(1:dim1,1:dim2)), globalmax)
               pstat(1,jj) = globalmin
               pstat(2,jj) = globalmax
               pstat(3,jj) = sqrt( globalsum / dimtot_global )
            else if (ndims == 3) then
               dim1 = dimsizes(1)
               dim2 = dimsizes(2)
               dim3 = dimsizes(3)
               call mpas_pool_get_field(pool_a, trim(poolitr % memberName), field3d)
               dimtot  = real(dim1*dim2*dim3,kind_real)
               prodtot = sum(field3d % array(1:dim1,1:dim2,1:dim3)**2 )
               call mpas_dmpar_sum_real(dminfo, dimtot, dimtot_global)
               call mpas_dmpar_sum_real(dminfo, prodtot, globalsum)
               call mpas_dmpar_min_real(dminfo, minval(field3d % array(1:dim1,1:dim2,1:dim3)), globalmin)
               call mpas_dmpar_max_real(dminfo, maxval(field3d % array(1:dim1,1:dim2,1:dim3)), globalmax)
               pstat(1,jj) = globalmin
               pstat(2,jj) = globalmax
               pstat(3,jj) = sqrt( globalsum / dimtot_global )
            end if
            deallocate(dimsizes)
         end if
      end if
   end do
 
   end subroutine da_gpnorm
 
 
   !***********************************************************************
   !
   !  subroutine da_fldrms
   !
   !> \brief   Performs basic statistics min/max/norm given a pool
   !> \author  Gael Descombes
   !> \date    February 2018
   !> \details
   !>  Given a pool of fields, return min/max/norm array
   !
   !-----------------------------------------------------------------------
 
   subroutine da_fldrms(pool_a, dminfo, fldrms, fld_select)
 
   implicit none
   type (mpas_pool_type), pointer, intent(in)  :: pool_a
   type (dm_info), pointer,        intent(in)  :: dminfo
   real(kind=kind_real),           intent(out) :: fldrms
   character (len=*), optional,    intent(in)  :: fld_select(:)
 
   type (mpas_pool_iterator_type) :: poolitr
   type (field1dreal), pointer :: field1d
   type (field2dreal), pointer :: field2d
   type (field3dreal), pointer :: field3d
   real(kind=kind_real) :: dimtot, dimtot_global, prodtot, prodtot_global
 
   integer :: ndims
   integer :: dim1, dim2, dim3
   integer, allocatable :: dimsizes(:)
 
   prodtot = mpas_jedi_zero_kr
   dimtot  = mpas_jedi_zero_kr
 
   !
   ! Iterate over all fields in pool_a
   ! named in pool_a
   !
   call mpas_pool_begin_iteration(pool_a)
 
   do while ( mpas_pool_get_next_member(pool_a, poolitr) )
      if (present(fld_select)) then
         if (ufo_vars_getindex(fld_select,trim(poolitr % memberName)) < 0) cycle
      end if
      if (poolitr % dataType == mpas_pool_real) then
         if (poolitr % memberType == mpas_pool_field) then
            ndims = poolitr % nDims
            dimsizes = getsolvedimsizes(pool_a, poolitr%memberName)
            if (ndims == 1) then
               dim1 = dimsizes(1)
               call mpas_pool_get_field(pool_a, trim(poolitr % memberName), field1d)
               dimtot  = dimtot + real(dim1,kind_real)
               prodtot = prodtot + sum( field1d % array(1:dim1)**2 )
            else if (ndims == 2) then
               dim1 = dimsizes(1)
               dim2 = dimsizes(2)
               call mpas_pool_get_field(pool_a, trim(poolitr % memberName), field2d)
               dimtot  = dimtot + real(dim1*dim2,kind_real)
               prodtot = prodtot + sum( field2d % array(1:dim1,1:dim2)**2 )
            else if (ndims == 3) then
               dim1 = dimsizes(1)
               dim2 = dimsizes(2)
               dim3 = dimsizes(3)
               call mpas_pool_get_field(pool_a, trim(poolitr % memberName), field3d)
               dimtot  = dimtot + real(dim1*dim2*dim3,kind_real)
               prodtot = prodtot + sum( field3d % array(1:dim1,1:dim2,1:dim3)**2 )
            end if
            deallocate(dimsizes)
         end if
      end if
   end do
 
   call mpas_dmpar_sum_real(dminfo, dimtot, dimtot_global)
   call mpas_dmpar_sum_real(dminfo, prodtot, prodtot_global)
   fldrms = sqrt(prodtot_global / dimtot_global)
 
   end subroutine da_fldrms
 
 
   !***********************************************************************
   !
   !  subroutine da_dot_product
   !
   !> \brief   Performs the dot_product given two pools of fields
   !> \author  Gael Descombes
   !> \date    February 2018
   !> \details
   !>  Given two pools of fields, compute the dot_product
   !
   !-----------------------------------------------------------------------
 
   subroutine da_dot_product(pool_a, pool_b, dminfo, zprod)
 
   implicit none
   type (mpas_pool_type), pointer, intent(in)  :: pool_a, pool_b
   type (dm_info), pointer,        intent(in)  :: dminfo
   real(kind=kind_real),           intent(out) :: zprod
 
   type (mpas_pool_iterator_type) :: poolitr
   type (field1dreal), pointer :: field1d_a, field1d_b
   type (field2dreal), pointer :: field2d_a, field2d_b
   type (field3dreal), pointer :: field3d_a, field3d_b
   real(kind=kind_real) :: fieldsum_local, zprod_local
 
   integer :: ndims
   integer :: dim1, dim2, dim3
   integer, allocatable :: dimsizes(:)
 
   !
   ! Iterate over all fields in pool_a
   ! named in pool_a
   !
   call mpas_pool_begin_iteration(pool_a)
 
   zprod_local = mpas_jedi_zero_kr
 
   do while ( mpas_pool_get_next_member(pool_a, poolitr) )
      if (poolitr % dataType == mpas_pool_real) then
         if (poolitr % memberType == mpas_pool_field) then
            ndims = poolitr % nDims
            dimsizes = getsolvedimsizes(pool_a, poolitr%memberName)
            !TODO: add check that dimSizes are the same between pool_a and pool_b for poolItr
            if (ndims == 1) then
               dim1 = dimsizes(1)
               call mpas_pool_get_field(pool_a, trim(poolitr % memberName), field1d_a)
               call mpas_pool_get_field(pool_b, trim(poolitr % memberName), field1d_b)
               fieldsum_local = sum(field1d_a % array(1:dim1) * field1d_b % array(1:dim1))
               zprod_local = zprod_local + fieldsum_local
            else if (ndims == 2) then
               dim1 = dimsizes(1)
               dim2 = dimsizes(2)
               call mpas_pool_get_field(pool_a, trim(poolitr % memberName), field2d_a)
               call mpas_pool_get_field(pool_b, trim(poolitr % memberName), field2d_b)
               fieldsum_local = sum(field2d_a % array(1:dim1,1:dim2) &
                                  * field2d_b % array(1:dim1,1:dim2))
               zprod_local = zprod_local + fieldsum_local
            else if (ndims == 3) then
               dim1 = dimsizes(1)
               dim2 = dimsizes(2)
               dim3 = dimsizes(3)
               call mpas_pool_get_field(pool_a, trim(poolitr % memberName), field3d_a)
               call mpas_pool_get_field(pool_b, trim(poolitr % memberName), field3d_b)
               fieldsum_local = sum(field3d_a % array(1:dim1,1:dim2,1:dim3) &
                                  * field3d_b % array(1:dim1,1:dim2,1:dim3))
               zprod_local = zprod_local + fieldsum_local
            end if
            deallocate(dimsizes)
         end if
      end if
   end do
 
   call mpas_dmpar_sum_real(dminfo, zprod_local, zprod)
 
   end subroutine da_dot_product
 
 
  subroutine cvt_oopsmpas_date(inString2,outString2,iconv)
 
     implicit none
 
     character (len=*), intent(in) :: instring2
     character (len=*), intent(inout) :: outstring2
     integer, intent(in) :: iconv
     integer :: i, curlen
     integer :: year, month, day, hour, minute, second
 
     character (len=ShortStrKIND) :: timepart
     character (len=ShortStrKIND) :: yearformat
     logical :: charexpand
     character (len=4) :: yyyy
     character (len=2) :: mm, dd, h, m, s
     character (len=21) :: outstring, instring
 
     ! 2017-08-08T00:00:00Z OOPS/YAML format
     ! 2010-10-24_02.00.00  MPAS format
     ! iconv=1: MPAS --> OOPS/YAML
     ! iconv=-1: OOPS/YAML --> MPAS
 
     if (iconv.eq.-1) then
        yyyy = instring2(1:4)
        mm   = instring2(6:7)
        dd   = instring2(9:10)
        h    = instring2(12:13)
        m    = instring2(15:16)
        s    = instring2(18:19)
     else
        yyyy = instring2(1:4)
        mm   = instring2(6:7)
        dd   = instring2(9:10)
        h    = instring2(12:13)
        m    = instring2(15:16)
        s    = instring2(18:19)
     end if
 
     write(message,*) 'cvt_oopsmpas_date instring: ',trim(yyyy),trim(mm),trim(dd),trim(h),trim(m),trim(s)
     call fckit_log%debug(message)
     write(message,*) 'cvt_oopsmpas_date input ',trim(instring2)
     call fckit_log%debug(message)
 
     write(outstring,*) ''
     instring = trim(outstring2)
 
     curlen = 0
     charexpand = .false.
     do i = 1, len_trim(instring)
           if (instring(i:i) == '$' ) then
               charexpand = .true.
           else if (instring(i:i) /= '$') then
               write(message,*) 'inString: ',trim(instring(i:i)),charexpand
               call fckit_log%debug(message)
               if (charexpand) then
                  select case (instring(i:i))
                     case ('Y')
                         outstring = trim(outstring) // trim(yyyy)
                     case ('M')
                         outstring = trim(outstring) // trim(mm)
                     case ('D')
                         outstring = trim(outstring) // trim(dd)
                     case ('h')
                         outstring = trim(outstring) // trim(h)
                     case ('m')
                         outstring = trim(outstring) // trim(m)
                     case ('s')
                         outstring = trim(outstring) // trim(s)
                     case default
                        call mpas_dmpar_global_abort('ERROR: mpas_timekeeping')
                  end select
                  curlen = len_trim(outstring)
                  charexpand = .false.
                  write(message,*) 'outString: ',trim(outstring)
                  call fckit_log%debug(message)
               else
                  outstring(curlen+1:curlen+1) = outstring2(i:i)
                  curlen = curlen+1
               end if
           end if
     end do
 
     outstring2 = trim(outstring)
     write(message,*) 'cvt_oopsmpas_date output ',trim(outstring2)
     call fckit_log%debug(message)
 
  end subroutine cvt_oopsmpas_date
 
 
 
! ------------------------------------------------------------------------
!  chunk of code from DART
! ------------------------------------------------------------------------
 
subroutine uv_cell_to_edges(domain, u_field, v_field, du, lonCell, latCell, &
                            &  nCells, edgeNormalVectors, nEdgesOnCell, edgesOnCell, nVertLevels)
 
   ! Project u, v wind increments at cell centers onto the edges.
   ! FIXME:
   !        we can hard-code R3 here since it comes from the (3d) x/y/z cartesian coordinate.
   !        We define nEdgesOnCell in get_grid_dims, and read edgesOnCell in get_grid.
   !        We read edgeNormalVectors in get_grid to use this subroutine.
   !        Here "U" is the prognostic variable in MPAS, and we update it with the wind
   !        increments at cell centers.
 
   implicit none
 
   type (domain_type), pointer, intent(inout) :: domain
   type (field2dreal), pointer, intent(in) :: u_field    ! u wind updated from filter
   type (field2dreal), pointer, intent(in) :: v_field    ! v wind updated from filter
   type (field2dreal), pointer, intent(inout) :: du       ! normal velocity increment on the edges
   real(kind_real), intent(in) :: loncell(1:ncells), latcell(1:ncells) ! lon, lat at cell centers in radians
   real(kind_real), intent(in) :: edgenormalvectors(:,:)
   integer, intent(in) :: nedgesoncell(:), edgesoncell(:,:)
   integer, intent(in) :: ncells, nvertlevels
 
   ! Local variables
   integer, parameter :: r3 = 3
   real(kind_real), dimension(:,:), allocatable :: east, north
   integer  :: icell, iedge, jedge, k
 
   ! allocation
   allocate(east(r3,ncells))
   allocate(north(r3,ncells))
 
   ! Initialization
   du%array(:,:) = mpas_jedi_zero_kr
 
   ! Compute unit vectors in east and north directions for each cell:
   do icell = 1, ncells
       east(1,icell) = -sin(loncell(icell))
       east(2,icell) =  cos(loncell(icell))
       east(3,icell) =  mpas_jedi_zero_kr
       call r3_normalize(east(1,icell), east(2,icell), east(3,icell))
       north(1,icell) = -cos(loncell(icell))*sin(latcell(icell))
       north(2,icell) = -sin(loncell(icell))*sin(latcell(icell))
       north(3,icell) =  cos(latcell(icell))
       call r3_normalize(north(1,icell), north(2,icell), north(3,icell))
   end do
 
   ! Project analysis increments from the cell centers to the edges
 
   do icell = 1, ncells
      do jedge = 1, nedgesoncell(icell)
         iedge = edgesoncell(jedge, icell)
            do k = 1, nvertlevels
               du%array(k,iedge) = du%array(k,iedge) + mpas_jedi_half_kr * u_field%array(k,icell)   &
                     * (edgenormalvectors(1,iedge) * east(1,icell)  &
                     +  edgenormalvectors(2,iedge) * east(2,icell)  &
                     +  edgenormalvectors(3,iedge) * east(3,icell)) &
                     + mpas_jedi_half_kr * v_field%array(k,icell)            &
                     * (edgenormalvectors(1,iedge) * north(1,icell) &
                     +  edgenormalvectors(2,iedge) * north(2,icell) &
                     +  edgenormalvectors(3,iedge) * north(3,icell))
            end do
      end do
   end do
 
   ! deallocation
   deallocate(east)
   deallocate(north)
 
end subroutine uv_cell_to_edges
 
!----------------------------------------------------------------------
 
subroutine r3_normalize(ax, ay, az)
 
   implicit none
 
   real(kind_real), intent(inout) :: ax, ay, az
   real(kind_real) :: mi
 
   mi = mpas_jedi_one_kr / sqrt(ax**2 + ay**2 + az**2)
 
   ax = ax * mi
   ay = ay * mi
   az = az * mi
 
end subroutine r3_normalize
 
!===============================================================================================================
 
end module mpas4da_mod