module shr_dmodel_mod 3,6

! !USES:

  use shr_sys_mod
  use shr_kind_mod, only: IN=>SHR_KIND_IN, R8=>SHR_KIND_R8, &
                               CS=>SHR_KIND_CS, CL=>SHR_KIND_CL
  use shr_log_mod, only: loglev  => shr_log_Level
  use shr_log_mod, only: logunit => shr_log_Unit
  use shr_mpi_mod, only: shr_mpi_bcast
  use mct_mod
!  use esmf_mod
  use perf_mod
  use pio

! !PUBLIC TYPES:
  implicit none
  private ! except

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

  public :: shr_dmodel_gsmapCreate
  public :: shr_dmodel_readLBUB
  public :: shr_dmodel_readgrid
  public :: shr_dmodel_gGridCompare
  public :: shr_dmodel_mapSet
  public :: shr_dmodel_translateAV
  public :: shr_dmodel_rearrGGrid

   interface shr_dmodel_mapSet; module procedure &
      shr_dmodel_mapSet_global
!      shr_dmodel_mapSet_dest
   end interface

  integer(IN),parameter,public :: shr_dmodel_gGridCompareXYabs      = 1 ! X,Y  relative error
  integer(IN),parameter,public :: shr_dmodel_gGridCompareXYrel      = 2 ! X,Y  absolute error
  integer(IN),parameter,public :: shr_dmodel_gGridCompareAreaAbs    = 3 ! area relative error
  integer(IN),parameter,public :: shr_dmodel_gGridCompareAreaRel    = 4 ! area absolute error
  integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskIdent  = 5 ! masks are identical
  integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskZeros  = 6 ! masks have same zeros
  integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskSubset = 7 ! mask is subset of other

  ! masked methods

  integer(IN),parameter,public :: shr_dmodel_gGridCompareXYabsMask      = 101 ! X,Y  relative error
  integer(IN),parameter,public :: shr_dmodel_gGridCompareXYrelMask      = 102 ! X,Y  absolute error
  integer(IN),parameter,public :: shr_dmodel_gGridCompareAreaAbsMask    = 103 ! area relative error
  integer(IN),parameter,public :: shr_dmodel_gGridCompareAreaRelMask    = 104 ! area absolute error
  integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskIdentMask  = 105 ! masks are identical
  integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskZerosMask  = 106 ! masks have same zeros
  integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskSubsetMask = 107 ! mask is subset of other

  integer(IN),parameter,public :: iotype_std_netcdf = -99 ! non pio option

!--------------------------------------------------------------------------
! Private data
!--------------------------------------------------------------------------

 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
CONTAINS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

subroutine shr_dmodel_gsmapCreate(gsmap,gsize,compid,mpicom,decomp) 1,1

    implicit none

    type(mct_gsMap), intent(inout) :: gsmap
    integer(IN)    , intent(in)    :: gsize
    integer(IN)    , intent(in)    :: compid
    integer(IN)    , intent(in)    :: mpicom
    character(len=*),intent(in)    :: decomp

    ! local

    integer(IN) :: n,npes,ierr
    integer(IN), pointer :: start(:)     ! for gsmap initialization
    integer(IN), pointer :: length(:)    ! for gsmap initialization
    integer(IN), pointer :: pe_loc(:)    ! for gsmap initialization
    character(*), parameter :: subname = '(shr_dmodel_gsmapCreate) '
    character(*), parameter :: F00   = "('(shr_dmodel_gsmapCreate) ',8a)"
    character(*), parameter :: F01   = "('(shr_dmodel_gsmapCreate) ',a,5i8)"

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

    if (gsize > 0) then
       call mpi_comm_size(mpicom,npes,ierr)
       allocate(start(npes),length(npes),pe_loc(npes))

       start = 0
       length = 0
       do n = 1,npes
          if (trim(decomp) == '1d') then
             length(n)  = gsize/npes
             if (n <= mod(gsize,npes)) length(n) = length(n) + 1
          elseif (trim(decomp) == 'root') then
             length = 0
             length(1) = gsize
          else
             write(logunit,F00) ' ERROR: decomp not allowed, ',trim(decomp)
             call shr_sys_abort(subname//' ERROR decomp')
          endif
          if (n == 1) then
             start(n) = 1
          else
             start(n) = start(n-1) + length(n-1)
          endif
          pe_loc(n) = n-1
       enddo
       call mct_gsMap_init( gsMap, COMPID, npes, gsize, start, length, pe_loc)
       deallocate(start,length,pe_loc)
    endif

end subroutine shr_dmodel_gsmapCreate
!===============================================================================


subroutine shr_dmodel_readgrid( gGrid, gsMap, nxgo, nygo, filename, compid, mpicom, & 6,20
  decomp, lonname, latname, maskname, areaname, fracname, readfrac, &
  scmmode, scmlon, scmlat)

  use seq_flds_mod, only : seq_flds_dom_coord, seq_flds_dom_other
  use shr_file_mod, only : shr_file_noprefix, shr_file_queryprefix, shr_file_get
  use shr_string_mod, only : shr_string_lastindex
  use shr_ncread_mod, only : shr_ncread_domain, shr_ncread_vardimsizes, &
                             shr_ncread_varexists, shr_ncread_vardimnum
  implicit none

  !----- arguments -----
  type(mct_gGrid), intent(inout) :: gGrid
  type(mct_gsMap), intent(inout) :: gsMap
  integer(IN)    , intent(out)   :: nxgo
  integer(IN)    , intent(out)   :: nygo
  character(len=*),intent(in)    :: filename
  integer(IN)    , intent(in)    :: compid
  integer(IN)    , intent(in)    :: mpicom
  character(len=*),optional,intent(in)    ::   decomp ! decomp strategy for gsmap
  character(len=*),optional,intent(in)    ::  lonname ! name of  lon variable in file
  character(len=*),optional,intent(in)    ::  latname ! name of  lat variable in file
  character(len=*),optional,intent(in)    :: maskname ! name of mask variable in file
  character(len=*),optional,intent(in)    :: areaname ! name of area variable in file
  character(len=*),optional,intent(in)    :: fracname ! name of frac variable in file
  logical         ,optional,intent(in)    :: readfrac ! T <=> also read frac  in file
  logical         ,optional,intent(in)    :: scmmode  ! single column mode
  real(R8)        ,optional,intent(in)    :: scmlon   ! single column lon
  real(R8)        ,optional,intent(in)    :: scmlat   ! single column lat

  !----- local -----
  integer(IN)   :: n,k,j,i    ! indices  
  integer(IN)   :: lsize      ! lsize
  integer(IN)   :: gsize      ! gsize
  integer(IN)   :: my_task, master_task
  integer(IN)   :: ierr       ! error code
  logical       :: fileexists ! 
  integer(IN)   :: rCode      ! return code
  character(CL) :: remoteFn   ! input file name (possibly at an archival location)
  character(CL) :: localFn    ! file name to be opened (possibly a local copy)
  character(CS) :: prefix     ! file prefix
  character(CS) ::   ldecomp  ! decomp strategy
  character(CS) ::  llatname  ! name of  lat variable
  character(CS) ::  llonname  ! name of  lon variable
  character(CS) :: lmaskname  ! name of mask variable
  character(CS) :: lareaname  ! name of area variable
  character(CS) :: lfracname  ! name of area variable
  logical       :: lreadfrac  ! read fraction
  logical       :: maskexists ! is mask on dataset
  integer(IN)   :: nxg,nyg    ! size of input fields
  integer(IN)   :: ndims      ! number of dims
  integer(IN)   :: nlon,nlat,narea,nmask,nfrac
  logical       :: lscmmode   ! local scm mode
  real(R8)      :: dist,mind  ! scmmode point search
  integer(IN)   :: ni,nj      ! scmmode point search
  real(R8)      :: lscmlon    ! local copy of scmlon

  real   (R8),allocatable ::  lon(:,:) ! temp array for domain lon  info
  real   (R8),allocatable ::  lat(:,:) ! temp array for domain lat  info
  integer(IN),allocatable :: mask(:,:) ! temp array for domain mask info
  real   (R8),allocatable :: area(:,:) ! temp array for domain area info
  real   (R8),allocatable :: frac(:,:) ! temp array for domain frac info

  integer(IN), pointer :: idata(:)   ! temporary
  type(mct_ggrid)      :: gGridRoot       ! global mct ggrid

  character(*), parameter :: subname = '(shr_dmodel_readgrid) '
  character(*), parameter :: F00   = "('(shr_dmodel_readgrid) ',8a)"
  character(*), parameter :: F01   = "('(shr_dmodel_readgrid) ',a,5i8)"

!-------------------------------------------------------------------------------
! PURPOSE:  Read MCT ggrid and set gsmap
!----------------------------------------------------------------------------
! Notes:
! o as per shr_file_get(), the file name format is expected to be
!   remoteFn = [location:][directory path]localFn
!   eg. "foobar.nc"  "/home/user/foobar.nc"  "mss:/USER/fobar.nc"
! o assumes a very specific netCDF domain file format wrt var names, etc.
!
! TO DO: have the calling routine select/input the domain's file name
!----------------------------------------------------------------------------

  call MPI_COMM_RANK(mpicom,my_task,ierr)
  master_task = 0

  lscmmode = .false.
  if (present(scmmode)) then
     lscmmode = scmmode
     if (lscmmode) then
        if (.not.present(scmlon) .or. .not.present(scmlat)) then
           write(logunit,*) subname,' ERROR: scmmode must supply scmlon and scmlat'
           call shr_sys_abort(subname//' ERROR: scmmode1 lon lat')
        endif
        if (my_task > 0) then
           write(logunit,*) subname,' ERROR: scmmode must be run on one pe'
           call shr_sys_abort(subname//' ERROR: scmmode2 tasks')
        endif
     endif
  endif

  if (my_task == master_task) then
     if ( shr_file_queryPrefix(fileName,prefix=prefix) /= shr_file_noPrefix ) then
        n        = max(len_trim(prefix),shr_string_lastIndex(fileName,"/"))
        remoteFn = fileName
        localFn  = fileName(n+1: len_trim(fileName) )
        call shr_file_get(rCode,localFn,remoteFn)
     else
        remoteFn = "undefined" ! this isn't needed
        localFn  = fileName    ! file to open
     end if
     inquire(file=trim(localFn),exist=fileExists)
     if (.not. fileExists) then
         write(logunit,F00) "ERROR: file does not exist: ", trim(localFn)
         call shr_sys_abort(subName//"ERROR: file does not exist: "//trim(localFn))
     end if
  endif

  lreadfrac = .false.
  ldecomp   = "1d"
  llonname  = "xc"  ! default values / standard data model domain file format
  llatname  = "yc"
  lmaskname = "mask"
  lareaname = "area"
  lfracname = "frac"
  if (present(  decomp))   ldecomp =   decomp
  if (present(readfrac)) lreadfrac = readfrac
  if (present( lonname))  llonname =  lonname
  if (present( latname))  llatname =  latname
  if (present(maskname)) lmaskname = maskname
  if (present(areaname)) lareaname = areaname
  if (present(fracname)) lfracname = fracname

  ! Initialize mct domain type
  ! lat/lon in degrees,  area in radians^2, mask is 1 (ocean), 0 (non-ocean)

  if (my_task == master_task) then
     if (shr_ncread_varexists(localFn,lmaskname)) then
        maskexists = .true.
        call shr_ncread_varDimSizes(localFn,lmaskname,nxg,nyg)
     else
        maskexists = .false.
        call shr_ncread_varDimNum(localFn,llonName,ndims)
        if (ndims == 1) then
           call shr_ncread_varDimSizes(localFn,llonName,nxg)
           call shr_ncread_varDimSizes(localFn,llatName,nyg)
        else
           call shr_ncread_varDimSizes(localFn,llonName,nxg,nyg)
        endif
     endif
  endif
  call shr_mpi_bcast(nxg,mpicom)
  call shr_mpi_bcast(nyg,mpicom)
  if (lscmmode) then
     nxgo = 1
     nygo = 1
     gsize = 1
  else
     nxgo = nxg
     nygo = nyg
     gsize = nxg*nyg
     if (gsize < 1) return
  endif

  call shr_dmodel_gsmapCreate(gsMap,gsize,compid,mpicom,trim(ldecomp))
  lsize = mct_gsMap_lsize(gsMap, mpicom)
  call mct_gGrid_init( GGrid=gGrid, CoordChars=trim(seq_flds_dom_coord), &
     OtherChars=trim(seq_flds_dom_other), lsize=lsize )

  ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid

  call mct_gsMap_orderedPoints(gsMap, my_task, idata)
  call mct_gGrid_importIAttr(gGrid,'GlobGridNum',idata,lsize)
  deallocate(idata)

  ! Initialize attribute vector with special value

  gGrid%data%rAttr = -9999.0_R8

  ! Load file data into domG then scatter to ggrid

  if (my_task == master_task) then

     allocate(lon(nxg,nyg))
     allocate(lat(nxg,nyg))
     allocate(area(nxg,nyg))
     allocate(mask(nxg,nyg))
     allocate(frac(nxg,nyg))

     if (.not.maskexists) then
        call shr_ncread_domain(localFn,llonName,lon,llatName,lat)
        mask = 1
        frac = 1.0_R8
        area = 1.0e36_R8
     else
        if (lreadfrac) then
           call shr_ncread_domain(localFn,llonName,lon,llatName,lat, &
                lmaskName,mask,lareaName,area,lfracName,frac)
        else ! assume frac = 1.0
           call shr_ncread_domain(localFn,llonName,lon,llatName,lat, &
                lmaskName,mask,lareaName,area)
           where (mask == 0)
              frac = 0.0_R8
           elsewhere
              frac = 1.0_R8
           end where
        endif
     endif

     call mct_gGrid_init(gGridRoot,gGrid,gsize)

! initialize gGridRoot to avoid errors when using strict compiler checks
     gGridRoot%data%rAttr = -9999.0_R8

     nlon  = mct_aVect_indexRA(gGridRoot%data,'lon')
     nlat  = mct_aVect_indexRA(gGridRoot%data,'lat')
     narea = mct_aVect_indexRA(gGridRoot%data,'area')
     nmask = mct_aVect_indexRA(gGridRoot%data,'mask')
     nfrac = mct_aVect_indexRA(gGridRoot%data,'frac')

     if (lscmmode) then
        !--- assumes regular 2d grid for compatability with shr_scam_getCloseLatLon ---
        !--- want lon values between 0 and 360, assume 1440 is enough ---
        lscmlon = mod(scmlon+1440.0_r8,360.0_r8)
        lon     = mod(lon   +1440.0_r8,360.0_r8)

        !--- start with wraparound ---
        ni = 1
        mind = abs(lscmlon - (lon(1,1)+360.0_r8))
        do i=1,nxg
           dist = abs(lscmlon - lon(i,1))
           if (dist < mind) then
              mind = dist
              ni = i
           endif
        enddo        

        nj = -1
        mind = 1.0e20
        do j=1,nyg
           dist = abs(scmlat - lat(1,j))
           if (dist < mind) then
              mind = dist
              nj = j
           endif
        enddo

        n = 1
        i = ni
        j = nj
        gGridRoot%data%rAttr(nlat ,n) = lat(i,j)
        gGridRoot%data%rAttr(nlon ,n) = lon(i,j)
        gGridRoot%data%rAttr(narea,n) = area(i,j)
        gGridRoot%data%rAttr(nmask,n) = real(mask(i,j),R8)
        gGridRoot%data%rAttr(nfrac,n) = frac(i,j)
     else
        n=0                                                                                    
        do j=1,nyg
        do i=1,nxg
           n=n+1                                                                               
           gGridRoot%data%rAttr(nlat ,n) = lat(i,j)
           gGridRoot%data%rAttr(nlon ,n) = lon(i,j)
           gGridRoot%data%rAttr(narea,n) = area(i,j)
           gGridRoot%data%rAttr(nmask,n) = real(mask(i,j),R8)
           gGridRoot%data%rAttr(nfrac,n) = frac(i,j)
        enddo
        enddo
     endif

  endif

  call mct_gGrid_scatter(gGridRoot, gGrid, gsMap, master_task, mpicom)
  if (my_task == master_task) call mct_gGrid_clean(gGridRoot)

end subroutine shr_dmodel_readgrid

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


subroutine shr_dmodel_readLBUB(stream,pio_subsystem,pio_iotype,pio_iodesc,mDate,mSec,mpicom,gsMap, & 1,13
                             avLB,mDateLB,mSecLB,avUB,mDateUB,mSecUB,newData,rmOldFile,istr)

  use shr_file_mod, only : shr_file_noprefix, shr_file_queryprefix, shr_file_get
  use shr_stream_mod
  implicit none

  !----- arguments -----
  type(shr_stream_streamType),intent(inout) :: stream
  type(iosystem_desc_t)      ,intent(inout) :: pio_subsystem
  integer(IN)                ,intent(in)    :: pio_iotype
  type(io_desc_t)            ,intent(inout) :: pio_iodesc
  integer(IN)                ,intent(in)    :: mDate  ,mSec
  integer(IN)                ,intent(in)    :: mpicom
  type(mct_gsMap)            ,intent(in)    :: gsMap
  type(mct_aVect)            ,intent(inout) :: avLB
  integer(IN)                ,intent(inout) :: mDateLB,mSecLB
  type(mct_aVect)            ,intent(inout) :: avUB
  integer(IN)                ,intent(inout) :: mDateUB,mSecUB
  logical                    ,intent(out)   :: newData
  logical,optional           ,intent(in)    :: rmOldFile
  character(len=*),optional  ,intent(in)    :: istr

  !----- local -----
  integer(IN)   :: n,k,j,i    ! indices  
  integer(IN)   :: lsize      ! lsize
  integer(IN)   :: gsize      ! gsize
  integer(IN)   :: my_task, master_task
  integer(IN)   :: ierr       ! error code
  integer(IN)   :: rCode      ! return code
  logical       :: localCopy,fileexists
  integer(IN)   :: ivals(6)   ! bcast buffer

  integer(IN)   :: oDateLB,oSecLB,dDateLB,oDateUB,oSecUB,dDateUB
  real(R8)      :: rDateM,rDateLB,rDateUB  ! model,LB,UB dates with fractional days
  integer(IN)   ::  n_lb, n_ub
  character(CL) :: fn_lb,fn_ub,fn_next,fn_prev
  character(CL) :: path
  character(len=32) :: lstr

  real(R8)      :: spd

  character(*), parameter :: subname = '(shr_dmodel_readLBUB) '
  character(*), parameter :: F00   = "('(shr_dmodel_readLBUB) ',8a)"
  character(*), parameter :: F01   = "('(shr_dmodel_readLBUB) ',a,5i8)"
  character(*), parameter :: F02   = "('(shr_dmodel_readLBUB) ',3a,i8)"

!-------------------------------------------------------------------------------
! PURPOSE:  Read LB and UB stream data
!----------------------------------------------------------------------------

  lstr = 'shr_dmodel_readLBUB'
  if (present(istr)) then
    lstr = trim(istr)
  endif

  call t_startf(trim(lstr)//'_setup')
  call MPI_COMM_RANK(mpicom,my_task,ierr)
  master_task = 0
  spd = 86400.0_R8

  newData = .false.
  n_lb = -1
  n_ub = -1
  fn_lb = 'undefinedlb'
  fn_ub = 'undefinedub'

  oDateLB = mDateLB
  oSecLB  = mSecLB
  oDateUB = mDateUB
  oSecUB  = mSecUB

  ! for backwards compatability (yikes)
  rDateM  = real(mDate  ,R8) + real(mSec  ,R8)/spd
  rDateLB = real(mDateLB,R8) + real(mSecLB,R8)/spd
  rDateUB = real(mDateUB,R8) + real(mSecUB,R8)/spd
  call t_stopf(trim(lstr)//'_setup')

  if (rDateM < rDateLB .or. rDateM > rDateUB) then
     call t_startf(trim(lstr)//'_fbound')
     if (my_task == master_task) then
!        call shr_stream_findBounds(stream,mDate,mSec,                 &
!                                   mDateLB,dDateLB,mSecLB,n_lb,fn_lb, &
!                                   mDateUB,dDateUB,mSecUB,n_ub,fn_ub  )
        call shr_stream_findBounds(stream,mDate,mSec,                 &
                                   ivals(1),dDateLB,ivals(2),ivals(5),fn_lb, &
                                   ivals(3),dDateUB,ivals(4),ivals(6),fn_ub  )
        call shr_stream_getFilePath(stream,path)
        localCopy = (shr_file_queryPrefix(path) /= shr_file_noPrefix )
!        write(logunit,*) 'tcx ivals ',ivals
     endif
     call t_stopf(trim(lstr)//'_fbound')
     call t_startf(trim(lstr)//'_bcast')

!    --- change 4 bcasts to a single bcast and copy for performance ---
!     call shr_mpi_bcast(mDateLB,mpicom)
!     call shr_mpi_bcast(mSecLB,mpicom)
!     call shr_mpi_bcast(mDateUB,mpicom)
!     call shr_mpi_bcast(mSecUB,mpicom)
     call shr_mpi_bcast(ivals,mpicom)
     mDateLB = ivals(1)
     mSecLB  = ivals(2)
     mDateUB = ivals(3)
     mSecUB  = ivals(4)
     n_lb    = ivals(5)
     n_ub    = ivals(6)
     call t_stopf(trim(lstr)//'_bcast')
  endif

  if (mDateLB /= oDateLB .or. mSecLB /= oSecLB) then
     newdata = .true.
     if (mDateLB == oDateUB .and. mSecLB == oSecUB) then
        call t_startf(trim(lstr)//'_LB_copy')
        avLB%rAttr(:,:) = avUB%rAttr(:,:)
        call t_stopf(trim(lstr)//'_LB_copy')
     else
        if (my_task == master_task) write(logunit,F02) 'reading file: ',trim(path),trim(fn_lb),n_lb
        call shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gsMap, avLB, mpicom, &
             path, fn_lb, n_lb,istr=trim(lstr)//'_LB')
     endif
  endif

  if (mDateUB /= oDateUB .or. mSecUB /= oSecUB) then
     newdata = .true.
     if (my_task == master_task) write(logunit,F02) 'reading file: ',trim(path),trim(fn_ub),n_ub
     call shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gsMap, avUB, mpicom, &
          path, fn_ub, n_ub,istr=trim(lstr)//'_UB')
  endif

  call t_startf(trim(lstr)//'_filemgt')
  !--- determine previous & next data files in list of files ---
  if (my_task == master_task .and. newdata) then
  call shr_stream_getFilePath(stream,path)
  localCopy = (shr_file_queryPrefix(path) /= shr_file_noPrefix )
  if (localCopy) then
     call shr_stream_getPrevFileName(stream,fn_lb,fn_prev,path)
     call shr_stream_getNextFileName(stream,fn_ub,fn_next,path)
     inquire(file=trim(fn_next),exist=fileExists)
     if ( trim(fn_next) == "unknown" .or. fileExists) then
        ! do nothing
     else 
        call shr_file_get(rCode,fn_next,trim(path)//fn_next,async=.true.)
        write(logunit,F00) "get next file: ",trim(fn_next)
     end if

     !---  remove the old file? (only if acquiring local copies) ---
     if ( rmOldFile .and. &
          fn_prev/=fn_lb .and. fn_prev/=fn_ub .and. fn_prev/=fn_next ) then 
        !--- previous file is not in use and is not next in list ---
        inquire(file=trim(fn_prev),exist=fileExists)
        if ( fileExists ) then
           call shr_sys_system(" rm "//trim(fn_prev),rCode)
           write(logunit,F00) "rm  prev file: ",trim(fn_prev)
        end if
     end if
  endif
  endif
  call t_stopf(trim(lstr)//'_filemgt')

end subroutine shr_dmodel_readLBUB

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

subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gsMap, av, mpicom, & 2,14
           path, fn, nt, istr)

  use shr_file_mod, only : shr_file_noprefix, shr_file_queryprefix, shr_file_get
  use shr_stream_mod
  use shr_ncread_mod
  implicit none

  !----- arguments -----
  type(shr_stream_streamType),intent(inout) :: stream
  type(iosystem_desc_t),intent(inout)  :: pio_subsystem
  integer(IN)     ,intent(in)    :: pio_iotype
  type(io_desc_t)      ,intent(inout)  :: pio_iodesc
  type(mct_gsMap) ,intent(in)    :: gsMap
  type(mct_aVect) ,intent(inout) :: av  
  integer(IN)     ,intent(in)    :: mpicom
  character(len=*),intent(in)    :: path
  character(len=*),intent(in)    :: fn
  integer(IN)     ,intent(in)    :: nt
  character(len=*),optional  ,intent(in)    :: istr

  !----- local -----
  integer(IN) :: my_task
  integer(IN) :: master_task
  integer(IN) :: ierr
  logical     :: localCopy,fileexists
  type(mct_avect) :: avG
  integer(IN) :: gsize,nx,ny
  integer(IN) :: k
  integer(IN) :: fid
  integer(IN) :: rCode      ! return code
  real(R8),allocatable :: data(:,:)
  character(CL) :: fileName
  character(CL) :: sfldName
  type(mct_avect) :: avtmp
  character(len=32) :: lstr

  integer(in) :: ndims
  integer(in),pointer :: dimid(:)
  type(file_desc_t) :: pioid
  type(var_desc_t) :: varid
  integer(kind=pio_offset) :: frame

  character(*), parameter :: subname = '(shr_dmodel_readstrm) '
  character(*), parameter :: F00   = "('(shr_dmodel_readstrm) ',8a)"
  character(*), parameter :: F01   = "('(shr_dmodel_readstrm) ',a,5i8)"

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

  lstr = 'shr_dmodel_readstrm'
  if (present(istr)) then
    lstr = trim(istr)
  endif

  call t_barrierf(trim(lstr)//'_BARRIER',mpicom)
  call t_startf(trim(lstr)//'_setup')
  call MPI_COMM_RANK(mpicom,my_task,ierr)
  master_task = 0

  gsize = mct_gsmap_gsize(gsMap)

  if (my_task == master_task) then
     localCopy = (shr_file_queryPrefix(path) /= shr_file_noPrefix )
     if (localCopy) then
        call shr_file_get(rCode,fn,trim(path)//fn)
        fileName = fn
     else                 ! DON'T acquire a local copy of the data file
        fileName = trim(path)//fn
     end if
     inquire(file=trim(fileName),exist=fileExists)
     if (.not. fileExists) then
        write(logunit,F00) "ERROR: file does not exist: ", trim(fileName)
        call shr_sys_abort(subName//"ERROR: file does not exist: "//trim(fileName))
     end if
  endif

  if (my_task == master_task) then
     call shr_stream_getFileFieldName(stream,1,sfldName)
  endif

  call t_stopf(trim(lstr)//'_setup')

  if (pio_iotype == iotype_std_netcdf) then

     call t_startf(trim(lstr)//'_readcdf')
     if (my_task == master_task) then
        call shr_ncread_varDimSizes(trim(fileName),trim(sfldName),nx,ny)
        if (gsize /= nx*ny) then
           write(logunit,F01) "ERROR in data sizes ",nx,ny,gsize
           call shr_sys_abort(subname//"ERROR in data sizes")
        endif
        call mct_aVect_init(avG,av,gsize)
        allocate(data(nx,ny))
        call shr_ncread_open(trim(fileName),fid,rCode)
        do k = 1,mct_aVect_nRAttr(av)
           call shr_stream_getFileFieldName(stream,k,sfldName)
           call shr_ncread_tField(fileName,nt,sfldName,data,fidi=fid,rc=rCode)
           avG%rAttr(k,:) = reshape(data, (/gsize/))
        enddo
        call shr_ncread_close(fid,rCode)
        deallocate(data)
     endif
     call t_stopf(trim(lstr)//'_readcdf')
     call t_barrierf(trim(lstr)//'_scatter'//'_BARRIER',mpicom)
     call t_startf(trim(lstr)//'_scatter')
     call mct_aVect_scatter(avG,avtmp,gsMap,master_task,mpicom)
     call mct_aVect_copy(avtmp,av)
     if (my_task == master_task) call mct_aVect_clean(avG)
     call mct_aVect_clean(avtmp)
     call t_stopf(trim(lstr)//'_scatter')

  else

     call t_startf(trim(lstr)//'_readpio')
     call shr_mpi_bcast(sfldName,mpicom,'sfldName')
     call shr_mpi_bcast(filename,mpicom,'filename')
     rcode = pio_openfile(pio_subsystem, pioid, pio_iotype, trim(filename), pio_nowrite)
     call pio_seterrorhandling(pioid,PIO_INTERNAL_ERROR)

     rcode = pio_inq_varid(pioid,trim(sfldName),varid)
     rcode = pio_inq_varndims(pioid, varid, ndims)
     allocate(dimid(ndims))
     rcode = pio_inq_vardimid(pioid, varid, dimid(1:ndims))
     rcode = pio_inq_dimlen(pioid, dimid(1), nx)
     rcode = pio_inq_dimlen(pioid, dimid(2), ny)
     deallocate(dimid)
     if (gsize /= nx*ny) then
        write(logunit,F01) "ERROR in data sizes ",nx,ny,gsize
        call shr_sys_abort(subname//"ERROR in data sizes")
     endif

     do k = 1,mct_aVect_nRAttr(av)
        if (my_task == master_task) then
           call shr_stream_getFileFieldName(stream,k,sfldName)
        endif
        call shr_mpi_bcast(sfldName,mpicom,'sfldName')
        rcode = pio_inq_varid(pioid,trim(sfldName),varid)
        frame = nt
        call pio_setframe(varid,frame)
        call pio_read_darray(pioid,varid,pio_iodesc,av%rattr(k,:),rcode)
     enddo

     call pio_closefile(pioid)
     call t_stopf(trim(lstr)//'_readpio')

  endif

end subroutine shr_dmodel_readstrm
!===============================================================================
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_dmodel_gGridCompare -- returns TRUE if two domains are the same.
!
! !DESCRIPTION:
!    Returns TRUE if two domains are the the same (within tolerance).
!
! !REVISION HISTORY:
!    2005-May-03 - B. Kauffman - added mulitiple methods
!    2005-May-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------


logical function shr_dmodel_gGridCompare(ggrid1,gsmap1,ggrid2,gsmap2,method,mpicom,eps) 3,6

! !USES:

   implicit none

! !INPUT/OUTPUT PARAMETERS:

   type(mct_gGrid)     ,intent(in)  :: ggrid1   ! 1st ggrid
   type(mct_gsmap)     ,intent(in)  :: gsmap1   ! 1st gsmap
   type(mct_gGrid)     ,intent(in)  :: ggrid2   ! 2nd ggrid
   type(mct_gsmap)     ,intent(in)  :: gsmap2   ! 2nd gsmap
   integer(IN)         ,intent(in)  :: method   ! selects what to compare
   integer(IN)         ,intent(in)  :: mpicom   ! mpicom
   real(R8)   ,optional,intent(in)  :: eps      ! epsilon compare value

!EOP

   !--- local ---
   real(R8)    :: leps         ! local epsilon
   integer(IN) :: n            ! counters
   integer(IN) :: my_task,master_task
   integer(IN) :: gsize
   integer(IN) :: ierr
   integer(IN) :: nlon1, nlon2, nlat1, nlat2, nmask1, nmask2  ! av field indices
   logical     :: compare      ! local compare logical
   real(R8)    :: lon1,lon2    ! longitudes to compare
   real(R8)    :: lat1,lat2    ! latitudes to compare
   real(R8)    :: msk1,msk2    ! masks to compare
   integer(IN) :: nx,ni1,ni2   ! i grid size, i offset for 1 vs 2 and 2 vs 1
   integer(IN) :: n1,n2,i,j    ! local indices
   type(mct_aVect) :: avG1     ! global av
   type(mct_aVect) :: avG2     ! global av
   integer(IN) :: lmethod      ! local method
   logical     :: maskmethod, maskpoint ! masking on method
 
   !--- formats ---
   character(*),parameter :: subName = '(shr_dmodel_gGridCompare) '
   character(*),parameter :: F01     = "('(shr_dmodel_gGridCompare) ',4a)" 

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

   call MPI_COMM_RANK(mpicom,my_task,ierr)
   master_task = 0

   leps = 1.0e-6_R8
   if (present(eps)) leps = eps

   lmethod = mod(method,100)
   if (method > 100) then
      maskmethod=.true.
   else
      maskmethod=.false.
   endif

   call mct_aVect_gather(gGrid1%data,avG1,gsmap1,master_task,mpicom)
   call mct_aVect_gather(gGrid2%data,avG2,gsmap2,master_task,mpicom)

   if (my_task == master_task) then

      compare = .true. 
      gsize = mct_aVect_lsize(avG1)
      if (gsize /= mct_aVect_lsize(avG2)) then
         compare = .false.
      endif

      if (.not. compare ) then
         !--- already failed the comparison test, check no futher ---
      else 
         nlon1 = mct_aVect_indexRA(avG1,'lon')
         nlat1 = mct_aVect_indexRA(avG1,'lat')
         nlon2 = mct_aVect_indexRA(avG2,'lon')
         nlat2 = mct_aVect_indexRA(avG2,'lat')
         nmask1 = mct_aVect_indexRA(avG1,'mask')
         nmask2 = mct_aVect_indexRA(avG2,'mask')

         ! To compare, want to be able to treat longitude wraparound generally.
         ! So we need to compute i index offset and we need to compute the size of the nx dimension
         ! First adjust the lon so it's in the range [0,360), add 1440 to lon to take into
         ! accounts lons less than 1440.  if any lon is less than -1440, abort.  1440 is arbitrary
         ! Next, comute ni1 and ni2.  These are the offsets of grid1 relative to grid2 and
         ! grid2 relative to grid1.  The sum of those offsets is nx.  Use ni1 to offset grid2
         ! in comparison and compute new grid2 index from ni1 and nx.  If ni1 is zero, then
         ! there is no offset, don't need to compute ni2, and nx can be anything > 0.

         !--- compute offset of grid2 compared to pt 1 of grid 1
         lon1 = minval(avG1%rAttr(nlon1,:))
         lon2 = minval(avG2%rAttr(nlon2,:))
         if ((lon1 < -1440.0_R8) .or. (lon2 < -1440.0_R8)) then
            write(logunit,*) subname,' ERROR: lon1 lon2 lt -1440 ',lon1,lon2
            call shr_sys_abort(subname//' ERROR: lon1 lon2 lt -1440')
         endif

         lon1 = mod(avG1%rAttr(nlon1,1)+1440.0_R8,360.0_R8)
         lat1 = avG1%rAttr(nlat1,1)
         ni1 = -1
         do n = 1,gsize
            lon2 = mod(avG2%rAttr(nlon2,n)+1440.0_R8,360.0_R8)
            lat2 = avG2%rAttr(nlat2,n)
            if ((ni1 < 0) .and. abs(lon1-lon2) <= leps .and. abs(lat1-lat2) <= leps) then
               ni1 = n - 1  ! offset, compare to first gridcell in grid 1
            endif
         enddo

         if (ni1 < 0) then        ! no match for grid point 1, so fails without going further
            compare = .false.
         elseif (ni1 == 0) then   ! no offset, set nx to anything > 0
            nx = 1
         else                     ! now compute ni2
            ! compute offset of grid1 compared to pt 1 of grid 2
            lon2 = mod(avG2%rAttr(nlon2,1)+1440.0_R8,360.0_R8)
            lat2 = avG2%rAttr(nlat2,1)
            ni2 = -1
            do n = 1,gsize
               lon1 = mod(avG1%rAttr(nlon1,n)+1440.0_R8,360.0_R8)
               lat1 = avG1%rAttr(nlat1,n)
               if ((ni2 < 0) .and. abs(lon1-lon2) <= leps .and. abs(lat1-lat2) <= leps) then
                  ni2 = n - 1  ! offset, compare to first gridcell in grid 1
               endif
            enddo
            if (ni2 < 0) then
               write(logunit,*) subname,' ERROR in ni2 ',ni1,ni2
               call shr_sys_abort(subname//' ERROR in ni2')
            endif
            nx = ni1 + ni2
         endif

         if (compare) then
            do n = 1,gsize
               j = ((n-1)/nx) + 1
               i = mod(n-1,nx) + 1
               n1 = (j-1)*nx + mod(n-1,nx) + 1
               n2 = (j-1)*nx + mod(n-1+ni1,nx) + 1
               if (n1 /= n) then    ! sanity check, could be commented out
                  write(logunit,*) subname,' ERROR in n1 n2 ',n,i,j,n1,n2
                  call shr_sys_abort(subname//' ERROR in n1 n2')
               endif
               lon1 = mod(avG1%rAttr(nlon1,n1)+1440.0_R8,360.0_R8)
               lat1 = avG1%rAttr(nlat1,n1)
               lon2 = mod(avG2%rAttr(nlon2,n2)+1440.0_R8,360.0_R8)
               lat2 = avG2%rAttr(nlat2,n2)
               msk1 = avG1%rAttr(nmask1,n1)
               msk2 = avG2%rAttr(nmask2,n2)

               maskpoint = .true.
               if (maskmethod .and. (msk1 == 0 .or. msk2 == 0)) then
                  maskpoint = .false.
               endif

               if (maskpoint) then
               if (lmethod == shr_dmodel_gGridCompareXYabs      ) then
                  if (abs(lon1 - lon2) > leps) compare = .false.
                  if (abs(lat1 - lat2) > leps) compare = .false.
               else if (lmethod == shr_dmodel_gGridCompareXYrel      ) then
                  if (rdiff(lon1,lon2) > leps) compare = .false.
                  if (rdiff(lat1,lat2) > leps) compare = .false.
               else if (lmethod == shr_dmodel_gGridCompareMaskIdent  ) then
                  if (msk1 /= msk2)compare = .false.
               else if (lmethod == shr_dmodel_gGridCompareMaskZeros  ) then
                  if (msk1 == 0 .and. msk2 /= 0) compare = .false.
                  if (msk1 /= 0 .and. msk2 == 0) compare = .false.
               else if (lmethod == shr_dmodel_gGridCompareMaskSubset ) then
                  if (msk1 /= 0 .and. msk2 == 0) compare = .false.
               else
                  write(logunit,F01) "ERROR: compare method not recognized, method = ",method
                  call shr_sys_abort(subName//"ERROR: compare method not recognized")
               endif  ! lmethod
               endif  ! maskpoint
            enddo ! gsize
         endif  ! compare
      endif   ! compare
   endif   ! master_task

   if (my_task == master_task) call mct_avect_clean(avG1)
   if (my_task == master_task) call mct_avect_clean(avG2)

   call shr_mpi_bcast(compare,mpicom)
   shr_dmodel_gGridCompare = compare

   return

!-------------------------------------------------------------------------------
contains   ! internal subprogram
!-------------------------------------------------------------------------------


   real(R8) function rdiff(v1,v2) ! internal function 3
      !------------------------------------------
      real(R8),intent(in) :: v1,v2                 ! two values to compare
      real(R8),parameter  :: c0           = 0.0_R8 ! zero
      real(R8),parameter  :: large_number = 1.0e20_R8 ! infinity
      !------------------------------------------
      if (v1 == v2) then
         rdiff = c0
      elseif (v1 == c0 .and. v2 /= c0) then
         rdiff = large_number
      elseif (v2 == c0 .and. v1 /= c0) then
         rdiff = large_number
      else
!        rdiff = abs((v2-v1)/v1)   ! old version, but rdiff(v1,v2) /= vdiff(v2,v1)
         rdiff = abs(2.0_R8*(v2-v1)/(v1+v2))
      endif
      !------------------------------------------
   end function rdiff

end function shr_dmodel_gGridCompare

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


subroutine shr_dmodel_mapSet_global(smatp,ggridS,gsmapS,nxgS,nygS,ggridD,gsmapD,nxgD,nygD, &,4
           name,type,algo,mask,vect,compid,mpicom,strategy)

    use shr_map_mod
    implicit none

    !----- arguments -----
    type(mct_sMatP), intent(inout) :: smatp
    type(mct_gGrid), intent(in)    :: ggridS
    type(mct_gsmap), intent(in)    :: gsmapS
    integer(IN)    , intent(in)    :: nxgS
    integer(IN)    , intent(in)    :: nygS
    type(mct_gGrid), intent(in)    :: ggridD
    type(mct_gsmap), intent(in)    :: gsmapD
    integer(IN)    , intent(in)    :: nxgD
    integer(IN)    , intent(in)    :: nygD
    character(len=*),intent(in)    :: name
    character(len=*),intent(in)    :: type
    character(len=*),intent(in)    :: algo
    character(len=*),intent(in)    :: mask
    character(len=*),intent(in)    :: vect
    integer(IN)    , intent(in)    :: compid
    integer(IN)    , intent(in)    :: mpicom
    character(len=*),intent(in),optional :: strategy

    !----- local -----

    integer(IN) :: n,i,j
    integer(IN) :: lsizeS,gsizeS,lsizeD,gsizeD
    integer(IN) :: nlon,nlat,nmsk
    integer(IN) :: my_task,master_task,ierr

    real(R8)   , pointer :: Xsrc(:,:)
    real(R8)   , pointer :: Ysrc(:,:)
    integer(IN), pointer :: Msrc(:,:)
    real(R8)   , pointer :: Xdst(:,:)
    real(R8)   , pointer :: Ydst(:,:)
    integer(IN), pointer :: Mdst(:,:)
    type(shr_map_mapType) :: shrmap
    type(mct_aVect) :: AVl
    type(mct_aVect) :: AVg

    character(len=32) :: lstrategy
    integer(IN) :: nsrc,ndst,nwts
    integer(IN), pointer :: isrc(:)
    integer(IN), pointer :: idst(:)
    real(R8)   , pointer :: wgts(:)
    type(mct_sMat) :: sMat0

    character(*), parameter :: subname = '(shr_dmodel_mapSet_global) '
    character(*), parameter :: F00   = "('(shr_dmodel_mapSet_global) ',8a)"
    character(*), parameter :: F01   = "('(shr_dmodel_mapSet_global) ',a,5i8)"

!-------------------------------------------------------------------------------
! PURPOSE:  Initialize sMatP from mct gGrid
!-------------------------------------------------------------------------------

    call MPI_COMM_RANK(mpicom,my_task,ierr)
    master_task = 0

    !--- get sizes and allocate for SRC ---

    lsizeS = mct_aVect_lsize(ggridS%data)
    call mct_avect_init(AVl,rList='lon:lat:mask',lsize=lsizeS)
    call mct_avect_copy(ggridS%data,AVl,rList='lon:lat:mask')
    call mct_avect_gather(AVl,AVg,gsmapS,master_task,mpicom)

    if (my_task == master_task) then
       gsizeS = mct_aVect_lsize(AVg)
       if (gsizeS /= nxgS*nygS) then
          write(logunit,F01) ' ERROR: gsizeS ',gsizeS,nxgS,nygS
          call shr_sys_abort(subname//' ERROR gsizeS')
       endif
       allocate(Xsrc(nxgS,nygS),Ysrc(nxgS,nygS),Msrc(nxgS,nygS))

       nlon = mct_avect_indexRA(AVg,'lon')
       nlat = mct_avect_indexRA(AVg,'lat')
       nmsk = mct_avect_indexRA(AVg,'mask')

       n = 0
       Msrc = 1
       do j = 1,nygS
       do i = 1,nxgS
          n = n + 1
          Xsrc(i,j) = AVg%rAttr(nlon,n)
          Ysrc(i,j) = AVg%rAttr(nlat,n)
          if (abs(AVg%rAttr(nmsk,n)) < 0.5_R8) Msrc(i,j) = 0
       enddo
       enddo
    endif

    if (my_task == master_task) call mct_aVect_clean(AVg)
    call mct_aVect_clean(AVl)

    !--- get sizes and allocate for DST ---

    lsizeD = mct_aVect_lsize(ggridD%data)
    call mct_avect_init(AVl,rList='lon:lat:mask',lsize=lsizeD)
    call mct_avect_copy(ggridD%data,AVl,rList='lon:lat:mask')
    call mct_avect_gather(AVl,AVg,gsmapD,master_task,mpicom)

    if (my_task == master_task) then
       gsizeD = mct_aVect_lsize(AVg)
       if (gsizeD /= nxgD*nygD) then
          write(logunit,F01) ' ERROR: gsizeD ',gsizeD,nxgD,nygD
          call shr_sys_abort(subname//' ERROR gsizeD')
       endif
       allocate(Xdst(nxgD,nygD),Ydst(nxgD,nygD),Mdst(nxgD,nygD))

       nlon = mct_avect_indexRA(AVg,'lon')
       nlat = mct_avect_indexRA(AVg,'lat')
       nmsk = mct_avect_indexRA(AVg,'mask')

       n = 0
       Mdst = 1
       do j = 1,nygD
       do i = 1,nxgD
          n = n + 1
          Xdst(i,j) = AVg%rAttr(nlon,n)
          Ydst(i,j) = AVg%rAttr(nlat,n)
          if (abs(AVg%rAttr(nmsk,n)) < 0.5_R8) Mdst(i,j) = 0
       enddo
       enddo
    endif

    if (my_task == master_task) call mct_aVect_clean(AVg)
    call mct_aVect_clean(AVl)

    !--- set map ---

    if (my_task == master_task) then
       call shr_map_mapSet(shrmap,Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst, &
          trim(name),trim(type),trim(algo),trim(mask),trim(vect))
       deallocate(Xsrc,Ysrc,Msrc)
       deallocate(Xdst,Ydst,Mdst)
    endif

    !--- convert map to sMatP ---

    lstrategy = 'Xonly'
    if (present(strategy)) then
       lstrategy = trim(strategy)
    endif

    if (my_task == master_task) then
       call shr_map_get(shrmap,shr_map_fs_nsrc,nsrc)
       call shr_map_get(shrmap,shr_map_fs_ndst,ndst)
       call shr_map_get(shrmap,shr_map_fs_nwts,nwts)
       allocate(isrc(nwts),idst(nwts),wgts(nwts))
       call shr_map_get(shrmap,isrc,idst,wgts)
       call shr_map_clean(shrmap)

       call mct_sMat_init(sMat0,ndst,nsrc,nwts)

       call mct_sMat_ImpGColI (sMat0,isrc,nwts)
       call mct_sMat_ImpGRowI (sMat0,idst,nwts)
       call mct_sMat_ImpMatrix(sMat0,wgts,nwts)
       deallocate(isrc,idst,wgts)
    endif

    call mct_sMatP_Init(sMatP,sMat0,gsmapS,gsmapD,lstrategy,master_task,mpicom,compid)

    if (my_task == master_task) then
       call mct_sMat_clean(sMat0)
    endif

end subroutine shr_dmodel_mapSet_global

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


subroutine shr_dmodel_mapSet_dest(smatp,ggridS,gsmapS,nxgS,nygS,ggridD,gsmapD,nxgD,nygD, &,5
           name,type,algo,mask,vect,compid,mpicom,strategy)

    use shr_map_mod
    implicit none

    !----- arguments -----
    type(mct_sMatP), intent(inout) :: smatp
    type(mct_gGrid), intent(in)    :: ggridS
    type(mct_gsmap), intent(in)    :: gsmapS
    integer(IN)    , intent(in)    :: nxgS
    integer(IN)    , intent(in)    :: nygS
    type(mct_gGrid), intent(in)    :: ggridD
    type(mct_gsmap), intent(in)    :: gsmapD
    integer(IN)    , intent(in)    :: nxgD
    integer(IN)    , intent(in)    :: nygD
    character(len=*),intent(in)    :: name
    character(len=*),intent(in)    :: type
    character(len=*),intent(in)    :: algo
    character(len=*),intent(in)    :: mask
    character(len=*),intent(in)    :: vect
    integer(IN)    , intent(in)    :: compid
    integer(IN)    , intent(in)    :: mpicom
    character(len=*),intent(in),optional :: strategy

    !----- local -----

    integer(IN) :: n,i,j
    integer(IN) :: lsizeS,gsizeS,lsizeD,gsizeD
    integer(IN) :: nlon,nlat,nmsk
    integer(IN) :: my_task,master_task,ierr

    real(R8)   , pointer :: Xsrc(:,:)
    real(R8)   , pointer :: Ysrc(:,:)
    integer(IN), pointer :: Msrc(:,:)
    real(R8)   , pointer :: Xdst(:)
    real(R8)   , pointer :: Ydst(:)
    integer(IN), pointer :: Mdst(:)
    type(shr_map_mapType) :: shrmap
    type(mct_aVect) :: AVl
    type(mct_aVect) :: AVg

    character(len=32) :: lstrategy
    integer(IN) :: nsrc,ndst,nwts
    integer(IN), pointer :: points(:)
    integer(IN), pointer :: isrc(:)
    integer(IN), pointer :: idst(:)
    real(R8)   , pointer :: wgts(:)
    type(mct_sMat) :: sMat0

    character(*), parameter :: subname = '(shr_dmodel_mapSet_dest) '
    character(*), parameter :: F00   = "('(shr_dmodel_mapSet_dest) ',8a)"
    character(*), parameter :: F01   = "('(shr_dmodel_mapSet_dest) ',a,5i8)"

!-------------------------------------------------------------------------------
! PURPOSE:  Initialize sMatP from mct gGrid
!-------------------------------------------------------------------------------

    call MPI_COMM_RANK(mpicom,my_task,ierr)
    master_task = 0

    !--- get sizes and allocate for SRC ---

    lsizeS = mct_aVect_lsize(ggridS%data)
    call mct_avect_init(AVl,rList='lon:lat:mask',lsize=lsizeS)
    call mct_avect_copy(ggridS%data,AVl,rList='lon:lat:mask')

    call mct_avect_gather(AVl,AVg,gsmapS,master_task,mpicom)

    allocate(Xsrc(nxgS,nygS),Ysrc(nxgS,nygS),Msrc(nxgS,nygS))
    if (my_task == master_task) then
       gsizeS = mct_aVect_lsize(AVg)
       if (gsizeS /= nxgS*nygS) then
          write(logunit,F01) ' ERROR: gsizeS ',gsizeS,nxgS,nygS
          call shr_sys_abort(subname//' ERROR gsizeS')
       endif

       nlon = mct_avect_indexRA(AVg,'lon')
       nlat = mct_avect_indexRA(AVg,'lat')
       nmsk = mct_avect_indexRA(AVg,'mask')

       n = 0
       Msrc = 1
       do j = 1,nygS
       do i = 1,nxgS
          n = n + 1
          Xsrc(i,j) = AVg%rAttr(nlon,n)
          Ysrc(i,j) = AVg%rAttr(nlat,n)
          if (abs(AVg%rAttr(nmsk,n)) < 0.5_R8) Msrc(i,j) = 0
       enddo
       enddo
    endif
    call shr_mpi_bcast(Xsrc,mpicom)
    call shr_mpi_bcast(Ysrc,mpicom)
    call shr_mpi_bcast(Msrc,mpicom)

    if (my_task == master_task) call mct_aVect_clean(AVg)
    call mct_aVect_clean(AVl)

    !--- get sizes and allocate for DST ---

    lsizeD = mct_aVect_lsize(ggridD%data)
    call mct_avect_init(AVl,rList='lon:lat:mask',lsize=lsizeD)
    call mct_avect_copy(ggridD%data,AVl,rList='lon:lat:mask')

#if (1 == 0)
    call mct_avect_gather(AVl,AVg,gsmapD,master_task,mpicom)

    if (my_task == master_task) then
       gsizeD = mct_aVect_lsize(AVg)
       if (gsizeD /= nxgD*nygD) then
          write(logunit,F01) ' ERROR: gsizeD ',gsizeD,nxgD,nygD
          call shr_sys_abort(subname//' ERROR gsizeD')
       endif
       allocate(Xdst(nxgD,nygD),Ydst(nxgD,nygD),Mdst(nxgD,nygD))

       nlon = mct_avect_indexRA(AVg,'lon')
       nlat = mct_avect_indexRA(AVg,'lat')
       nmsk = mct_avect_indexRA(AVg,'mask')

       n = 0
       Mdst = 1
       do j = 1,nygD
       do i = 1,nxgD
          n = n + 1
          Xdst(i,j) = AVg%rAttr(nlon,n)
          Ydst(i,j) = AVg%rAttr(nlat,n)
          if (abs(AVg%rAttr(nmsk,n)) < 0.5_R8) Mdst(i,j) = 0
       enddo
       enddo
    endif

    if (my_task == master_task) call mct_aVect_clean(AVg)
#endif

    allocate(Xdst(lsizeD),Ydst(lsizeD),Mdst(lsizeD))

    nlon = mct_avect_indexRA(AVl,'lon')
    nlat = mct_avect_indexRA(AVl,'lat')
    nmsk = mct_avect_indexRA(AVl,'mask')

    Mdst = 1
    do n = 1,lsizeD
       Xdst(n) = AVl%rAttr(nlon,n)
       Ydst(n) = AVl%rAttr(nlat,n)
       if (abs(AVl%rAttr(nmsk,n)) < 0.5_R8) Mdst(n) = 0
    enddo

    call mct_aVect_clean(AVl)

    !--- set map ---

    nsrc = nxgS*nygS
    ndst = nxgD*nygD
    call mct_gsmap_orderedPoints(gsmapD,my_task,points)
    if (size(points) /= size(Xdst)) then
       write(logunit,F01) ' ERROR: gsizeD ',size(points),size(Xdst)
       call shr_sys_abort(subname//' ERROR points size')
    endif
    call shr_map_mapSet(shrmap,Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst,ndst,points, &
       trim(name),trim(type),trim(algo),trim(mask),trim(vect))
    deallocate(points)
    deallocate(Xsrc,Ysrc,Msrc)
    deallocate(Xdst,Ydst,Mdst)

    !--- convert map to sMatP ---

    lstrategy = 'Xonly'
    if (present(strategy)) then
       lstrategy = trim(strategy)
    endif

    call shr_map_get(shrmap,shr_map_fs_nwts,nwts)
    allocate(isrc(nwts),idst(nwts),wgts(nwts))
    call shr_map_get(shrmap,isrc,idst,wgts)
    call shr_map_clean(shrmap)

    call mct_sMat_init(sMat0,ndst,nsrc,nwts)

    call mct_sMat_ImpLColI (sMat0,isrc,nwts)
    call mct_sMat_ImpLRowI (sMat0,idst,nwts)
    call mct_sMat_ImpMatrix(sMat0,wgts,nwts)
    deallocate(isrc,idst,wgts)

    call mct_sMatP_Init(sMatP,sMat0,gsmapS,gsmapD,master_task,mpicom,compid)

    call mct_sMat_clean(sMat0)

end subroutine shr_dmodel_mapSet_dest

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


subroutine shr_dmodel_rearrGGrid( ggridi, ggrido, gsmap, rearr, mpicom )

  implicit none

  !----- arguments -----
  type(mct_ggrid), intent(in)    :: ggridi
  type(mct_ggrid), intent(inout) :: ggrido  
  type(mct_gsmap), intent(in)    :: gsmap
  type(mct_rearr), intent(in)    :: rearr
  integer(IN)    , intent(in)    :: mpicom

  !----- local -----
  integer(IN)          :: lsize      ! lsize
  real(R8)   , pointer :: data(:)    ! temporary
  integer(IN), pointer :: idata(:)   ! temporary
  integer(IN)          :: my_task    ! local pe number
  integer(IN)          :: ier        ! error code
  character(*), parameter :: subname = '(shr_dmodel_rearrGGrid) '

!-------------------------------------------------------------------------------
! PURPOSE:  Determine MCT ggrid
!-------------------------------------------------------------------------------

  ! Initialize mct ggrid type
  ! lat/lon in degrees,  area in radians^2, mask is 1 (ocean), 0 (non-ocean)

  call mpi_comm_rank(mpicom,my_task,ier)

  lsize = mct_gsMap_lsize(gsMap, mpicom)
  call mct_gGrid_init( ggrido, ggridi, lsize=lsize )

  ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid

  call mct_gsMap_orderedPoints(gsMap, my_task, idata)
  call mct_gGrid_importIAttr(ggrido,'GlobGridNum',idata,lsize)
  deallocate(idata)

  ! Initialize attribute vector with special value

  allocate(data(lsize))

  data(:) = -9999.0_R8 
  call mct_gGrid_importRAttr(ggrido,"lat"  ,data,lsize) 
  call mct_gGrid_importRAttr(ggrido,"lon"  ,data,lsize) 
  call mct_gGrid_importRAttr(ggrido,"area" ,data,lsize) 
  call mct_gGrid_importRAttr(ggrido,"aream",data,lsize) 
  data(:) = 0.0_R8     
  call mct_gGrid_importRAttr(ggrido,"mask",data,lsize) 
  call mct_gGrid_importRAttr(ggrido,"frac",data,lsize) 

  deallocate(data)

  call mct_rearr_rearrange(ggridi%data, ggrido%data, rearr)

end subroutine shr_dmodel_rearrGGrid

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


subroutine shr_dmodel_translateAV( avi, avo, avifld, avofld, rearr )

  implicit none

  !----- arguments -----
  type(mct_aVect), intent(in)    :: avi       ! input av
  type(mct_aVect), intent(inout) :: avo       ! output av
  character(len=*),intent(in)    :: avifld(:) ! input field names for translation
  character(len=*),intent(in)    :: avofld(:) ! output field names for translation
  type(mct_rearr), intent(in),optional :: rearr     ! rearranger for diff decomp

  !----- local -----
  integer(IN)      :: n,k,ka,kb,kc,cnt  ! indices  
  integer(IN)      :: lsize      ! lsize
  integer(IN)      :: gsize      ! gsize
  integer(IN)      :: nflds      ! number of fields in avi

  type(mct_aVect)  :: avtri,avtro ! translated av on input/output grid
  character(CL)    :: ilist      ! input list for translation
  character(CL)    :: olist      ! output list for translation
  character(CL)    :: cfld       ! character field name
  type(mct_string) :: sfld       ! string field
  integer(IN)      :: ktrans
  character(*), parameter :: subname = '(shr_dmodel_translateAV) '

!-------------------------------------------------------------------------------
! PURPOSE:  Fill avo from avi
!-------------------------------------------------------------------------------

  if (size(avifld) /= size(avofld)) then
     write(logunit,*) subname,' ERROR": avi and avo fld list ',size(avifld),size(avofld)
  endif
  ktrans = size(avifld)

  ! generate fld lists
  nflds = mct_aVect_nRattr(avi)
  cnt = 0
  do ka = 1,nflds
     call mct_aVect_getRList(sfld,ka,avi)
     cfld = mct_string_toChar(sfld)
     call mct_string_clean(sfld)

     k = 0
     kb = 0
     kc = 0
     do while (kb == 0 .and. k < ktrans)
        k = k + 1
        if (trim(avifld(k)) == trim(cfld)) then
           kb = k
           kc = mct_aVect_indexRA(avo,trim(avofld(kb)),perrWith='quiet')
           if (ka > 0 .and. kc > 0) then
              cnt = cnt + 1
              if (cnt == 1) then
                 ilist = trim(avifld(kb))
                 olist = trim(avofld(kb))
              else
                 ilist = trim(ilist)//':'//trim(avifld(kb))
                 olist = trim(olist)//':'//trim(avofld(kb))
              endif
           endif
        endif
     enddo
  enddo

  if (cnt > 0) then
     lsize = mct_avect_lsize(avi)
     call mct_avect_init(avtri,rlist=olist,lsize=lsize)
     call mct_avect_Copy(avi,avtri,rList=ilist,TrList=olist)

     if (present(rearr)) then
        lsize = mct_avect_lsize(avo)
        call mct_avect_init(avtro,rlist=olist,lsize=lsize)
        call mct_avect_zero(avtro)
        call mct_rearr_rearrange(avtri, avtro, rearr)
        call mct_avect_Copy(avtro,avo)
        call mct_aVect_clean(avtro)
     else
        call mct_avect_Copy(avtri,avo)
     endif

     call mct_aVect_clean(avtri)
  endif

end subroutine shr_dmodel_translateAV
!===============================================================================

end module shr_dmodel_mod