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