module seq_domain_mct

  use shr_kind_mod, only: R8=>shr_kind_r8, IN=>shr_kind_in
  use shr_kind_mod, only: CL=>shr_kind_cl
  use shr_sys_mod,  only: shr_sys_flush, shr_sys_abort
  use shr_mpi_mod,  only: shr_mpi_min, shr_mpi_max

  use mct_mod
  use seq_comm_mct
  use seq_infodata_mod
  use seq_map_mod     , only: seq_map_map
  use seq_map_type_mod, only: seq_map

  use component_type_mod

  implicit none
  private ! except
#include <mpif.h>
  save

  !--------------------------------------------------------------------------
  ! Public interfaces
  !--------------------------------------------------------------------------

  public :: seq_domain_check
  public :: seq_domain_compare
  public :: seq_domain_areafactinit

  !--------------------------------------------------------------------------
  ! Public variables
  !--------------------------------------------------------------------------

  real(R8), parameter :: eps_tiny   = 1.0e-16_R8 ! roundoff eps
  real(R8), parameter :: eps_big    = 1.0e+02_R8 ! big eps
  real(R8), parameter :: eps_frac_samegrid = 1.0e-9_R8 ! epsilon for fractions for samegrid

  !--------------------------------------------------------------------------
  ! Private interfaces
  !--------------------------------------------------------------------------

  private :: seq_domain_check_grid

  !================================================================================
contains
  !================================================================================

  !================================================================================

  subroutine seq_domain_check( infodata, &
       atm, ice, lnd, ocn, rof, glc, &
       samegrid_al, samegrid_ao, samegrid_ro, samegrid_lg)

    !-----------------------------------------------------------
    ! Uses
    !
    use prep_atm_mod, only: prep_atm_get_mapper_Fi2a
    use prep_atm_mod, only: prep_atm_get_mapper_Fl2a
    use prep_atm_mod, only: prep_atm_get_mapper_Fo2a
    use prep_lnd_mod, only: prep_lnd_get_mapper_Fa2l
    use prep_ocn_mod, only: prep_ocn_get_mapper_SFi2o
    use prep_glc_mod, only: prep_glc_get_mapper_Fl2g
    !
    ! Arguments
    !
    type (seq_infodata_type) , intent(inout) :: infodata
    type(component_type)     , intent(in)    :: atm
    type(component_type)     , intent(in)    :: ice
    type(component_type)     , intent(in)    :: lnd
    type(component_type)     , intent(in)    :: ocn
    type(component_type)     , intent(in)    :: rof
    type(component_type)     , intent(in)    :: glc
    logical                  , intent(in)    :: samegrid_al ! atm lnd grid same
    logical                  , intent(in)    :: samegrid_ao ! atm ocn grid same
    logical                  , intent(in)    :: samegrid_ro ! rof ocn grid same
    logical                  , intent(in)    :: samegrid_lg ! lnd glc grid same
    !
    ! Local variables
    !
    type(seq_map)   , pointer :: mapper_i2a ! inout needed for lower methods
    type(seq_map)   , pointer :: mapper_i2o ! inout needed for lower methods
    type(seq_map)   , pointer :: mapper_o2a !
    type(seq_map)   , pointer :: mapper_l2g !
    type(seq_map)   , pointer :: mapper_a2l !
    type(seq_map)   , pointer :: mapper_l2a !
    !
    type(mct_gGrid) , pointer :: atmdom_a   ! atm domain
    type(mct_gGrid) , pointer :: icedom_i   ! ice domain
    type(mct_gGrid) , pointer :: lnddom_l   ! lnd domain
    type(mct_gGrid) , pointer :: ocndom_o   ! ocn domain
    type(mct_gGrid) , pointer :: glcdom_g   ! glc domain
    !
    type(mct_gsMap) , pointer :: gsMap_a    ! atm global seg map
    type(mct_gsMap) , pointer :: gsMap_i    ! ice global seg map
    type(mct_gsMap) , pointer :: gsMap_l    ! lnd global seg map
    type(mct_gsMap) , pointer :: gsMap_o    ! ocn global seg map
    type(mct_gsMap) , pointer :: gsMap_r    ! ocn global seg map
    type(mct_gsMap) , pointer :: gsMap_g    ! glc global seg map
    !
    type(mct_gGrid) :: lnddom_a              ! lnd domain info on atm decomp
    type(mct_gGrid) :: lnddom_g              ! lnd domain info on glc decomp
    type(mct_gGrid) :: icedom_a              ! ice domain info on atm decomp (all grids same)
    type(mct_gGrid) :: ocndom_a              ! ocn domain info on atm decomp (all grids same)
    type(mct_gGrid) :: icedom_o              ! ocn domain info on ocn decomp (atm/ocn grid different)
    !
    real(R8), pointer :: fracl(:)            ! land fraction on atm decomp
    real(R8), pointer :: fraco(:)            ! ocn  fraction on atm decomp
    real(R8), pointer :: fraci(:)            ! ice  fraction on atm decomp
    real(R8), pointer :: maskl(:)            ! land mask on atm decomp (all grids same)
    real(R8), pointer :: maski(:)            ! ice  mask on atm decomp (all grids same)
    real(R8), pointer :: masko(:)            ! ocn  mask on atm decomp (all grids same)
    !
    integer(IN) :: n            ! indicies
    !
    integer(IN) :: mpicom_cplid
    !
    logical      :: atm_present              ! atm present flag
    logical      :: lnd_present              ! lnd present flag
    logical      :: ocn_present              ! ocn present flag
    logical      :: ice_present              ! ice present flag
    logical      :: glc_present              ! glc present flag
    logical      :: rof_present              ! rof present flag
    logical      :: ocnrof_prognostic        ! ocn rof prognostic flag
    integer(IN)  :: rcode                    ! error status
    integer(IN)  :: atmsize                  ! local  size of atm  grid
    integer(IN)  :: lndsize                  ! local  size of land grid
    integer(IN)  :: ocnsize                  ! local  size of ocn  grid
    integer(IN)  :: icesize                  ! local  size of ice  grid
    integer(IN)  :: glcsize                  ! local  size of glc  grid
    integer(IN)  :: gatmsize                 ! global size of atm  grid
    integer(IN)  :: glndsize                 ! global size of land grid
    integer(IN)  :: gocnsize                 ! global size of ocn  grid
    integer(IN)  :: grofsize                 ! global size of ocn  grid
    integer(IN)  :: gicesize                 ! global size of ice  grid
    integer(IN)  :: gglcsize                 ! global size of glc  grid
    integer(IN)  :: npts                     ! local size temporary
    real(R8)     :: diff,dmaxo,dmaxi         ! difference tracker
    logical      :: iamroot                  ! local masterproc
    real(R8)     :: eps_frac                 ! epsilon for fractions
    real(R8)     :: eps_axmask               ! epsilon for masks, atm/lnd
    real(R8)     :: eps_axgrid               ! epsilon for grid coords, atm/lnd
    real(R8)     :: eps_axarea               ! epsilon for areas, atm/lnd
    real(R8)     :: eps_oimask               ! epsilon for masks, ocn/ice
    real(R8)     :: eps_oigrid               ! epsilon for grid coords, ocn/ice
    real(R8)     :: eps_oiarea               ! epsilon for areas, ocn/ice
    real(R8)     :: my_eps_frac              ! local eps_frac value
    !
    real(R8),allocatable :: mask (:)         ! temporary real vector, domain mask
    !
    character(*),parameter :: F00 = "('(seq_domain_check) ',4a)"
    character(*),parameter :: F01 = "('(seq_domain_check) ',a,i6,a)"
    character(*),parameter :: F02 = "('(seq_domain_check) ',a,g23.15)"
    character(*),parameter :: F0R = "('(seq_domain_check) ',2A,2g23.15,A )"
    character(*),parameter :: subName = '(seq_domain_check) '
    !-----------------------------------------------------------

    mapper_i2a => prep_atm_get_mapper_Fi2a()
    mapper_i2o => prep_ocn_get_mapper_SFi2o()
    mapper_o2a => prep_atm_get_mapper_Fo2a()
    mapper_l2g => prep_glc_get_mapper_Fl2g()
    mapper_a2l => prep_lnd_get_mapper_Fa2l()
    mapper_l2a => prep_atm_get_mapper_Fl2a()

    call seq_comm_setptrs(CPLID,iamroot=iamroot, mpicom=mpicom_cplid)

    call seq_infodata_GetData( infodata,      &
         lnd_present=lnd_present,             &
         ocn_present=ocn_present,             &
         ice_present=ice_present,             &
         glc_present=glc_present,             &
         atm_present=atm_present,             &
         rof_present=rof_present,             &
         ocnrof_prognostic=ocnrof_prognostic, &
         eps_frac=eps_frac,                   &
         eps_amask=eps_axmask,                &
         eps_agrid=eps_axgrid,                &
         eps_aarea=eps_axarea,                &
         eps_omask=eps_oimask,                &
         eps_ogrid=eps_oigrid,                &
         eps_oarea=eps_oiarea )

    ! Get info

    if (atm_present) then
       gsmap_a  => component_get_gsmap_cx(atm) ! gsmap_ax
       atmdom_a => component_get_dom_cx(atm)   ! dom_ax
       atmsize  = mct_avect_lsize(atmdom_a%data)
       gatmsize = mct_gsMap_gsize(gsMap_a)
    end if

    if (atm_present .and. lnd_present) then
       gsmap_l  => component_get_gsmap_cx(lnd) ! gsmap_lx
       lnddom_l => component_get_dom_cx(lnd)   ! dom_lx
       lndsize  = mct_avect_lsize(lnddom_l%data)
       glndsize = mct_gsMap_gsize(gsMap_l)

       if (samegrid_al .and. gatmsize /= glndsize) then
          write(logunit,*) subname,' error: global atmsize = ',&
               gatmsize,' global lndsize= ',glndsize
          call shr_sys_flush(logunit)
          call shr_sys_abort(subname//' atm and lnd grid must have the same global size')
       end if
       if (iamroot) write(logunit,F00) ' --- checking land maskfrac ---'
       call seq_domain_check_fracmask(lnddom_l%data)
       call mct_gGrid_init(oGGrid=lnddom_a, iGGrid=lnddom_l, lsize=atmsize)
       call mct_aVect_zero(lnddom_a%data)
       call seq_map_map(mapper_l2a, lnddom_l%data, lnddom_a%data, norm=.false., &
            omit_nonlinear=.true.)
       allocate(maskl(atmsize),stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' allocate maskl')
       allocate(fracl(atmsize),stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' allocate fracl')
       call mct_aVect_exportRAttr(lnddom_a%data, 'mask', maskl, atmsize)
       call mct_aVect_exportRAttr(lnddom_a%data, 'frac', fracl, atmsize)
    endif

    if (atm_present .and. ocn_present) then
       gsmap_o  => component_get_gsmap_cx(ocn) ! gsmap_ox
       ocndom_o => component_get_dom_cx(ocn)   ! dom_ox
       ocnsize  = mct_avect_lsize(ocndom_o%data)
       gocnsize = mct_gsMap_gsize(gsMap_o)

       if (samegrid_ao .and. gatmsize /= gocnsize) then
          write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global ocnsize= ',gocnsize
          call shr_sys_flush(logunit)
          call shr_sys_abort(subname//' atm and ocn grid must have the same global size')
       end if
       if (iamroot) write(logunit,F00) ' --- checking ocean maskfrac ---'
       call seq_domain_check_fracmask(ocndom_o%data)
       call mct_gGrid_init(oGGrid=ocndom_a, iGGrid=ocndom_o, lsize=atmsize)
       call mct_aVect_zero(ocndom_a%data)
       call seq_map_map(mapper_o2a, ocndom_o%data, ocndom_a%data, norm=.false., &
            omit_nonlinear=.true.)
       allocate(masko(atmsize),stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' allocate masko')
       allocate(fraco(atmsize),stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' allocate fraco')
       call mct_aVect_exportRAttr(ocndom_a%data, 'mask', masko, atmsize)
       if (samegrid_ao) then
          call mct_aVect_exportRattr(ocndom_a%data, 'frac', fraco, atmsize)
       else
          call mct_aVect_exportRattr(ocndom_a%data, 'mask', fraco, atmsize)
       endif
    endif

    if (atm_present .and. ice_present) then
       gsmap_i  => component_get_gsmap_cx(ice) ! gsmap_ix
       icedom_i => component_get_dom_cx(ice)   ! dom_ix
       icesize  = mct_avect_lsize(icedom_i%data)
       gicesize = mct_gsMap_gsize(gsMap_i)

       if (samegrid_ao .and. gatmsize /= gicesize) then
          write(logunit,*) subname,' error: global atmsize = ',&
               gatmsize,' global icesize= ',gicesize
          call shr_sys_flush(logunit)
          call shr_sys_abort(subname//' atm and ice grid must have the same global size')
       end if
       if (iamroot) write(logunit,F00) ' --- checking ice maskfrac ---'
       call seq_domain_check_fracmask(icedom_i%data)
       call mct_gGrid_init(oGGrid=icedom_a, iGGrid=icedom_i, lsize=atmsize)
       call mct_aVect_zero(icedom_a%data)
       call seq_map_map(mapper_i2a, icedom_i%data, icedom_a%data, norm=.false., &
            omit_nonlinear=.true.)
       allocate(maski(atmsize),stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' allocate maski')
       allocate(fraci(atmsize),stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' allocate fraci')
       call mct_aVect_exportRAttr(icedom_a%data, 'mask', maski, atmsize)
       if (samegrid_ao) then
          call mct_aVect_exportRattr(icedom_a%data, 'frac', fraci, atmsize)
       else
          call mct_aVect_exportRattr(icedom_a%data, 'mask', fraci, atmsize)
       endif
    endif

    if (lnd_present .and. glc_present) then
       gsmap_l  => component_get_gsmap_cx(lnd) ! gsmap_lx
       lnddom_l => component_get_dom_cx(lnd)   ! dom_lx
       lndsize  = mct_avect_lsize(lnddom_l%data)
       glndsize = mct_gsMap_gsize(gsMap_l)

       gsmap_g  => component_get_gsmap_cx(glc) ! gsmap_gx
       glcdom_g => component_get_dom_cx(glc)   ! dom_gx
       glcsize  = mct_avect_lsize(glcdom_g%data)
       gglcsize = mct_gsMap_gsize(gsMap_g)

       if (samegrid_lg .and. gglcsize /= glndsize) then
          write(logunit,*) subname,' error: global glcsize = ',gglcsize,' global lndsize= ',glndsize
          call shr_sys_flush(logunit)
          call shr_sys_abort(subname//' glc and lnd grid must have the same global size')
       end if

       if (iamroot) write(logunit,F00) ' --- checking glc maskfrac ---'
       call seq_domain_check_fracmask(glcdom_g%data)
       if (iamroot) write(logunit,F00) ' --- checking lnd maskfrac ---'
       call seq_domain_check_fracmask(lnddom_l%data)

       if (samegrid_lg) then
          call mct_gGrid_init(oGGrid=lnddom_g, iGGrid=lnddom_l, lsize=glcsize)
          call mct_aVect_zero(lnddom_g%data)
          call seq_map_map(mapper_l2g, lnddom_l%data, lnddom_g%data, norm=.false., &
               omit_nonlinear=.true.)
          if (iamroot) write(logunit,F00) ' --- checking glc/lnd domains ---'
          npts = glcsize
          allocate(mask(npts),stat=rcode)
          if(rcode /= 0) call shr_sys_abort(subname//' allocate mask')
          call mct_aVect_getRAttr(lnddom_g%data,"mask",mask,rcode)
          where (mask < eps_axmask) mask = 0.0_R8
          call seq_domain_check_grid(glcdom_g%data, lnddom_g%data, 'mask', eps=eps_axmask, mpicom=mpicom_cplid, mask=mask)
          call seq_domain_check_grid(glcdom_g%data, lnddom_g%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_cplid, mask=mask)
          call seq_domain_check_grid(glcdom_g%data, lnddom_g%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_cplid, mask=mask)
          call seq_domain_check_grid(glcdom_g%data, lnddom_g%data, 'area', eps=eps_axarea, mpicom=mpicom_cplid, mask=mask)
          deallocate(mask,stat=rcode)
          if(rcode /= 0) call shr_sys_abort(subname//' deallocate mask')
       end if

    endif

    if (ice_present .and. ocn_present) then
       gsmap_i  => component_get_gsmap_cx(ice) ! gsmap_ix
       icedom_i => component_get_dom_cx(ice)   ! dom_ix
       icesize  = mct_avect_lsize(icedom_i%data)
       gicesize = mct_gsMap_gsize(gsMap_i)

       gsmap_o  => component_get_gsmap_cx(ocn) ! gsmap_ox
       ocndom_o => component_get_dom_cx(ocn)   ! dom_ox
       ocnsize  = mct_avect_lsize(ocndom_o%data)
       gocnsize = mct_gsMap_gsize(gsMap_o)

       if (gocnsize /= gicesize) then
          write(logunit,*) subname,' error: global ocnsize = ',gocnsize,' global icesize= ',gicesize
          call shr_sys_flush(logunit)
          call shr_sys_abort(subname//' ocean and ice grid must have the same global size')
       endif
       call mct_gGrid_init(oGGrid=icedom_o, iGGrid=icedom_i, lsize=ocnsize)
       call mct_aVect_zero(icedom_o%data)
       call seq_map_map(mapper_i2o, icedom_i%data, icedom_o%data, norm=.false., &
            omit_nonlinear=.true.)
    end if

    if (rof_present .and. ocnrof_prognostic .and. samegrid_ro) then
       gsmap_r  => component_get_gsmap_cx(glc) ! gsmap_gx
       grofsize = mct_gsMap_gsize(gsMap_r)

       if (gocnsize /= grofsize) then
          write(logunit,*) subname,' error: global ocnsize = ',gocnsize,' global rofsize= ',grofsize
          call shr_sys_flush(logunit)
          call shr_sys_abort(subname//' ocean and rof grid must have the same global size')
       endif
    end if

    !------------------------------------------------------------------------------
    ! Check ice/ocean grid consistency
    !------------------------------------------------------------------------------

    if (ocn_present .and. ice_present) then
       !    if (samegrid_oi) then       ! doesn't yet exist

       npts = ocnsize
       allocate(mask(npts),stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' allocate mask')

       if (iamroot) write(logunit,F00) ' --- checking ocn/ice domains ---'
       call seq_domain_check_grid(ocndom_o%data, icedom_o%data,'mask', eps=eps_oigrid, mpicom=mpicom_cplid)
       call mct_aVect_getRAttr(ocndom_o%data,"mask",mask,rcode)
       where (mask < eps_oimask) mask = 0.0_R8

       call seq_domain_check_grid(ocndom_o%data, icedom_o%data,'lat' , eps=eps_oigrid, mpicom=mpicom_cplid, mask=mask)
       call seq_domain_check_grid(ocndom_o%data, icedom_o%data,'lon' , eps=eps_oigrid, mpicom=mpicom_cplid, mask=mask)
       call seq_domain_check_grid(ocndom_o%data, icedom_o%data,'area', eps=eps_oiarea, mpicom=mpicom_cplid, mask=mask)

       deallocate(mask,stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' deallocate mask')

       !    endif
    endif

    !------------------------------------------------------------------------------
    ! Check atm/lnd grid consistency
    !------------------------------------------------------------------------------

    if (atm_present .and. lnd_present .and. samegrid_al) then
       if (iamroot) write(logunit,F00) ' --- checking atm/land domains ---'
       call seq_domain_check_grid(atmdom_a%data, lnddom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_cplid, mask=maskl)
       call seq_domain_check_grid(atmdom_a%data, lnddom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_cplid, mask=maskl)
       call seq_domain_check_grid(atmdom_a%data, lnddom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_cplid, mask=maskl)
    endif

    !------------------------------------------------------------------------------
    ! Check atm/ocn and atm/ice grid consistency (if samegrid)
    !------------------------------------------------------------------------------

    if (atm_present .and. ice_present .and. samegrid_ao) then
       if (iamroot) write(logunit,F00) ' --- checking atm/ice domains ---'
       call seq_domain_check_grid(atmdom_a%data, icedom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_cplid, mask=maski)
       call seq_domain_check_grid(atmdom_a%data, icedom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_cplid, mask=maski)
       call seq_domain_check_grid(atmdom_a%data, icedom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_cplid, mask=maski)
    endif

    if (atm_present .and. ocn_present .and. samegrid_ao) then
       if (iamroot) write(logunit,F00) ' --- checking atm/ocn domains ---'
       call seq_domain_check_grid(atmdom_a%data, ocndom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_cplid, mask=masko)
       call seq_domain_check_grid(atmdom_a%data, ocndom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_cplid, mask=masko)
       call seq_domain_check_grid(atmdom_a%data, ocndom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_cplid, mask=masko)
    endif

    !------------------------------------------------------------------------------
    ! Check consistency of land fraction with ocean mask on grid
    !------------------------------------------------------------------------------

    if (atm_present) then
       my_eps_frac = eps_frac
       if (samegrid_ao) my_eps_frac = eps_frac_samegrid
       if (.not. samegrid_al) my_eps_frac = eps_big

       if (iamroot) write(logunit,F00) ' --- checking fractions in domains ---'
       dmaxi = 0.0_R8
       dmaxo = 0.0_R8
       do n = 1,atmsize
          if (lnd_present .and. ice_present) then
             diff = abs(1._R8 - fracl(n) - fraci(n))
             dmaxi = max(diff,dmaxi)
             if (diff > my_eps_frac) then
                write(logunit,*)'inconsistency between land fraction and sea ice fraction'
                write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraci= ',fraci(n),' sum= ',fracl(n)+fraci(n)
                call shr_sys_flush(logunit)
                call shr_sys_abort(subname//' inconsistency between land fraction and sea ice fraction')
             end if
             if ((1._R8-fraci(n)) > eps_frac .and. fracl(n) < eps_tiny) then
                write(logunit,*)'inconsistency between land mask and sea ice mask'
                write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraci= ',fraci(n)
                call shr_sys_flush(logunit)
                call shr_sys_abort(subname//'  inconsistency between land mask and sea ice mask')
             end if
          endif
          if (lnd_present .and. ocn_present) then
             diff = abs(1._R8 - fracl(n) - fraco(n))
             dmaxo = max(diff,dmaxo)
             if (diff > my_eps_frac) then
                write(logunit,*)'inconsistency between land fraction and ocn land fraction'
                write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraco= ',fraco(n),' sum= ',fracl(n)+fraco(n)
                call shr_sys_flush(logunit)
                call shr_sys_abort(subname//'  inconsistency between land fraction and ocn land fraction')
             end if
             if ((1._R8-fraco(n)) > eps_frac .and. fracl(n) < eps_tiny) then
                write(logunit,*)'inconsistency between land mask and ocn land mask'
                write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraco= ',fraco(n)
                call shr_sys_flush(logunit)
                call shr_sys_abort(subname//'  inconsistency between land mask and ocn land mask')
             end if
          endif
       end do
       if (iamroot) then
          write(logunit,F02) ' maximum           difference for ofrac sum ',dmaxo
          write(logunit,F02) ' maximum           difference for ifrac sum ',dmaxi
          write(logunit,F02) ' maximum allowable difference for  frac sum ',my_eps_frac
          write(logunit,F02) ' maximum allowable tolerance for valid frac ',eps_frac
          call shr_sys_flush(logunit)
       end if
    end if

    !------------------------------------------------------------------------------
    ! Clean up allocated memory
    !------------------------------------------------------------------------------

    if (atm_present .and. lnd_present) then
       deallocate(fracl,stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' deallocate fracl')
       deallocate(maskl,stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' deallocate maskl')
       call mct_gGrid_clean(lnddom_a, rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' clean lnddom_a')
    endif

    if (atm_present .and. ocn_present) then
       deallocate(fraco,stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' deallocate fraco')
       deallocate(masko,stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' deallocate masko')
       call mct_gGrid_clean(ocndom_a, rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' clean ocndom_a')
    endif

    if (atm_present .and. ice_present) then
       deallocate(fraci,stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' deallocate fraci')
       deallocate(maski,stat=rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' deallocate maski')
       call mct_gGrid_clean(icedom_a, rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' clean icedom_o')
    endif

    if (ocn_present .and. ice_present) then
       call mct_gGrid_clean(icedom_o, rcode)
       if(rcode /= 0) call shr_sys_abort(subname//' clean icedom_o')
    endif

    call shr_sys_flush(logunit)

  end subroutine seq_domain_check

  !===============================================================================

  subroutine seq_domain_compare(dom1, dom2, mpicom, eps)

    !-----------------------------------------------------------

    ! Arguments

    type(mct_gGrid)  , intent(in) :: dom1
    type(mct_gGrid)  , intent(in) :: dom2
    integer(IN)      , intent(in) :: mpicom
    real(R8),optional, intent(in) :: eps    ! error condition for compare

    ! Local variables
    real(R8) :: leps
    character(*),parameter :: F00 = "('(seq_domain_compare) ',4a)"
    character(*),parameter :: F01 = "('(seq_domain_compare) ',a,i12,a)"
    character(*),parameter :: F02 = "('(seq_domain_compare) ',2a,g23.15)"
    character(*),parameter :: F0R = "('(seq_domain_compare) ',2A,2g23.15,A )"
    character(*),parameter :: subName = '(seq_domain_compare) '

    leps = eps_tiny
    if (present(eps)) then
       leps = eps
    endif

    call seq_domain_check_grid(dom1%data, dom2%data, 'mask', eps=leps, mpicom=mpicom)
    call seq_domain_check_grid(dom1%data, dom2%data, 'lat' , eps=leps, mpicom=mpicom)
    call seq_domain_check_grid(dom1%data, dom2%data, 'lon' , eps=leps, mpicom=mpicom)
    call seq_domain_check_grid(dom1%data, dom2%data, 'area', eps=leps, mpicom=mpicom)

  end subroutine seq_domain_compare

  !===============================================================================

  subroutine seq_domain_check_fracmask(dom1)

    !-----------------------------------------------------------

    ! Arguments

    type(mct_aVect) , intent(in) :: dom1

    ! Local variables
    integer(in) :: n,npts,ndiff
    integer(in) :: rcode
    real(R8), pointer :: dmask(:)           ! temporaries
    real(R8), pointer :: dfrac(:)           ! temporaries

    character(*),parameter :: F00 = "('(seq_domain_check_fracmask) ',4a)"
    character(*),parameter :: F01 = "('(seq_domain_check_fracmask) ',a,i12,a)"
    character(*),parameter :: F02 = "('(seq_domain_check_fracmask) ',2a,g23.15)"
    character(*),parameter :: F0R = "('(seq_domain_check_fracmask) ',2A,2g23.15,A )"
    character(*),parameter :: subName = '(seq_domain_check_fracmask) '
    !-----------------------------------------------------------

    npts = mct_aVect_lsize(dom1)

    allocate(dmask(npts),stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' allocate dmask')
    allocate(dfrac(npts),stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' allocate dfrac')

    call mct_aVect_exportRAttr(dom1, 'mask', dmask, npts)
    call mct_aVect_exportRAttr(dom1, 'frac', dfrac, npts)

    ndiff = 0
    do n = 1,npts
       if (abs(dfrac(n)) > eps_tiny .and. abs(dmask(n)) < eps_tiny) then
          !debug            write(logunit,*)'n= ',n,' dfrac= ',dfrac(n),' dmask= ',dmask(n)
          ndiff = ndiff + 1
       endif
    enddo

    if (ndiff > 0) then
       write(logunit,*) trim(subname)," ERROR: incompatible domain mask and frac values"
       call shr_sys_flush(logunit)
       call shr_sys_abort(subName//" incompatible domain mask and frac values")
    endif

    deallocate(dmask,stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' deallocate dmask')
    deallocate(dfrac,stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' deallocate dfrac')

  end subroutine seq_domain_check_fracmask

  !===============================================================================

  subroutine seq_domain_check_grid(dom1, dom2, attr, eps, mpicom, mask)

    !-----------------------------------------------------------

    ! Arguments

    type(mct_aVect) , intent(in) :: dom1
    type(mct_aVect) , intent(in) :: dom2
    character(len=*), intent(in) :: attr   ! grid attribute to compare
    real(R8)        , intent(in) :: eps    ! error condition for compare
    integer(IN)     , intent(in) :: mpicom
    real(R8)        , intent(in), optional :: mask(:)

    ! Local variables

    integer(in)       :: n,ndiff            ! indices
    integer(in)       :: npts1,npts2,npts   ! counters
    integer(in)       :: rcode              ! error code
    real(R8)          :: diff,max_diff      ! temporaries
    real(R8)          :: tot_diff           ! maximum diff across all pes
    integer(IN)       :: ier                ! error code
    real(R8), pointer :: data1(:)           ! temporaries
    real(R8), pointer :: data2(:)           ! temporaries
    real(R8), pointer :: lmask(:)           ! temporaries
    logical           :: iamroot            ! local masterproc

    character(*),parameter :: F00 = "('(seq_domain_check_grid) ',4a)"
    character(*),parameter :: F01 = "('(seq_domain_check_grid) ',a,i12,a)"
    character(*),parameter :: F02 = "('(seq_domain_check_grid) ',2a,g23.15)"
    character(*),parameter :: F0R = "('(seq_domain_check_grid) ',2A,2g23.15,A )"
    character(*),parameter :: subName = '(seq_domain_check_grid) '
    !-----------------------------------------------------------

    call seq_comm_setptrs(CPLID,iamroot=iamroot)

    npts1 = mct_aVect_lsize(dom1)
    npts2 = mct_aVect_lsize(dom2)
    npts  = npts1

    if (npts1 == npts2) then
       if (iamroot) write(logunit,F01) " the domain size is = ", npts
    else
       write(logunit,*) trim(subname)," domain size #1 = ", npts1
       write(logunit,*) trim(subname)," domain size #2 = ", npts2
       write(logunit,*) trim(subname)," ERROR: domain size mis-match"
       call shr_sys_abort(subName//" ERROR: domain size mis-match")
    end if

    allocate(data1(npts),stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' allocate data1')
    allocate(data2(npts),stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' allocate data2')
    allocate(lmask(npts),stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' allocate lmask')

    call mct_aVect_exportRAttr(dom1, trim(attr), data1, npts)
    call mct_aVect_exportRAttr(dom2, trim(attr), data2, npts)
    lmask = 1.0_R8
    if (present(mask)) then
       if (size(mask) /= npts) then
          call shr_sys_abort(subName//" ERROR: mask size mis-match")
       endif
       lmask = mask
    endif

    ! --- adjust lons to address wraparound issues, we're assuming degree here! ---

    if (trim(attr) == "lon") then
       do n = 1,npts
          if (data2(n) > data1(n)) then
             do while ( (data1(n)+360.0_R8) < (data2(n)+180.0_R8) ) ! longitude is periodic
                data1(n) = data1(n) + 360.0_R8
             end do
          else
             do while ( (data2(n)+360.0_R8) < (data1(n)+180.0_R8) ) ! longitude is periodic
                data2(n) = data2(n) + 360.0_R8
             end do
          endif
       enddo
    endif

    ! Only check consistency where mask is greater than zero, if mask is present

    max_diff = 0.0_R8
    ndiff = 0
    do n=1,npts
       if (lmask(n) > eps_tiny) then
          diff = abs(data1(n)-data2(n))
          max_diff = max(max_diff,diff)
          if (diff > eps) then
             write(logunit,150) n,data1(n),data2(n),diff,eps
             ndiff = ndiff + 1
          endif
       end if
    end do
150 format('seq_domain_check_grid - n:',I3,' d1:',F12.6,' d2:',F12.6,' diff:',F18.14,' eps:',F18.14)

    call mpi_reduce(max_diff,tot_diff,1,MPI_REAL8,MPI_MAX,0,mpicom,ier)
    if (iamroot) then
       write(logunit,F02) " maximum           difference for ",trim(attr),tot_diff
       write(logunit,F02) " maximum allowable difference for ",trim(attr),eps
       call shr_sys_flush(logunit)
    endif
    call mpi_barrier(mpicom,ier)

    if (ndiff > 0) then
       write(logunit,*) trim(subname)," ERROR: incompatible domain grid coordinates"
       call shr_sys_flush(logunit)
       call shr_sys_abort(subName//" incompatible domain grid coordinates")
    endif

    deallocate(data1,stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' deallocate data1')
    deallocate(data2,stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' deallocate data2')
    deallocate(lmask,stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' deallocate lmask')

  end subroutine seq_domain_check_grid

  !===============================================================================

  subroutine seq_domain_areafactinit(domain, mdl2drv, drv2mdl, &
       samegrid, mpicom, iamroot, comment)
    !-----------------------------------------------------------
    !
    ! Arguments
    !
    type(mct_gGrid)  , pointer             :: domain     ! component domain on component pes
    real(R8)         , pointer             :: mdl2drv(:) ! comp->cpl factor on component pes
    real(R8)         , pointer             :: drv2mdl(:) ! cpl->comp factor on component pes
    logical          , intent(in)          :: samegrid   ! true => two grids are same
    integer          , intent(in)          :: mpicom     ! mpi communicator on component pes
    logical          , intent(in)          :: iamroot
    character(len=*) , optional,intent(in) :: comment
    !
    ! Local variables
    !
    integer                :: j1,j2,m1,n,rcode
    integer                :: gridsize
    real(R8)               :: rmin1,rmax1,rmin,rmax
    real(R8)               :: rmask,rarea,raream
    character(cl)          :: lcomment
    character(len=*),parameter :: subName = '(seq_domain_areafactinit) '
    character(len=*),parameter :: F0R = "(2A,2g23.15,A )"
    !
    !-----------------------------------------------------------

    lcomment = ''
    if (present(comment)) lcomment = comment

    ! get sizes

    gridsize = mct_gGrid_lsize(domain)
    allocate(drv2mdl(gridsize),mdl2drv(gridsize),stat=rcode)
    if(rcode /= 0) call shr_sys_abort(subname//' allocate area correction factors')

    j1 = mct_gGrid_indexRA(domain,"area"    ,dieWith=subName)
    j2 = mct_gGrid_indexRA(domain,"aream"   ,dieWith=subName)
    m1 = mct_gGrid_indexRA(domain,"mask"    ,dieWith=subName)

    mdl2drv(:)=1.0_R8
    drv2mdl(:)=1.0_R8

    if (samegrid) then
       ! default 1.0
    else
       do n=1,gridsize
          rmask  = domain%data%rAttr(m1,n)
          rarea  = domain%data%rAttr(j1,n)
          raream = domain%data%rAttr(j2,n)
          if ( abs(rmask) >= 1.0e-06) then
             if (rarea * raream /= 0.0_R8) then
                mdl2drv(n) = rarea/raream
                drv2mdl(n) = 1.0_R8/mdl2drv(n)
                !if (mdl2drv(n) > 10.0 .or. mdl2drv(n) < 0.1) then
                !   write(logunit,*) trim(subname),' WARNING area,aream= ', &
                !      domain%data%rAttr(j1,n),domain%data%rAttr(j2,n),' in ',n,gridsize
                !endif
             else
                write(logunit,*) trim(subname),' ERROR area,aream= ', &
                     rarea,raream,' in ',n,gridsize
                call shr_sys_flush(logunit)
                call shr_sys_abort()
             endif
          endif
       enddo
    end if

    rmin1 = minval(mdl2drv)
    rmax1 = maxval(mdl2drv)
    call shr_mpi_min(rmin1,rmin,mpicom)
    call shr_mpi_max(rmax1,rmax,mpicom)
    if (iamroot) write(logunit,F0R) trim(subname),' : min/max mdl2drv ',rmin,rmax,trim(lcomment)

    rmin1 = minval(drv2mdl)
    rmax1 = maxval(drv2mdl)
    call shr_mpi_min(rmin1,rmin,mpicom)
    call shr_mpi_max(rmax1,rmax,mpicom)
    if (iamroot) write(logunit,F0R) trim(subname),' : min/max drv2mdl ',rmin,rmax,trim(lcomment)
    if (iamroot) call shr_sys_flush(logunit)

  end subroutine seq_domain_areafactinit

  !===============================================================================

end module seq_domain_mct
