module shr_dmodel_mod ! !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, & CX=>SHR_KIND_CX 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 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_translateAV_list public :: shr_dmodel_translate_list 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) 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, & 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 deallocate(lon) deallocate(lat) deallocate(area) deallocate(mask) deallocate(frac) 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, & 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_const_mod, only : shr_const_cDay 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 = shr_const_cday newData = .false. n_lb = -1 n_ub = -1 fn_lb = 'undefinedlb' fn_ub = 'undefinedub' oDateLB = mDateLB oSecLB = mSecLB oDateUB = mDateUB oSecUB = mSecUB 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 ) 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(stream%calendar,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) then write(logunit,F02) 'reading file: ',trim(path),trim(fn_lb),n_lb call shr_sys_flush(logunit) endif 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) then write(logunit,F02) 'reading file: ',trim(path),trim(fn_ub),n_ub call shr_sys_flush(logunit) endif 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) call shr_sys_flush(logunit) 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) call shr_sys_flush(logunit) 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, & 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) ! !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 !------------------------------------------ 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, & 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, & 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(CX) :: ilist ! input list for translation character(CX) :: olist ! output list for translation character(CX) :: 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) call shr_sys_flush(logunit) 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 !=============================================================================== subroutine shr_dmodel_translate_list( avi, avo, avifld, avofld, ilist, olist, cnt) 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 character(CL) ,intent(out) :: ilist ! input list for translation character(CL) ,intent(out) :: olist ! output list for translation integer(IN) ,intent(out) :: cnt ! indices !----- local ----- integer(IN) :: n,k,ka,kb,kc ! indices integer(IN) :: lsize ! lsize integer(IN) :: nflds ! number of fields in avi 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) call shr_sys_flush(logunit) 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 end subroutine shr_dmodel_translate_list !=============================================================================== subroutine shr_dmodel_translateAV_list( avi, avo, ilist, olist, rearr ) implicit none !----- arguments ----- type(mct_aVect), intent(in) :: avi ! input av type(mct_aVect), intent(inout) :: avo ! output av character(CL) ,intent(in) :: ilist ! input list for translation character(CL) ,intent(in) :: olist ! output list for translation type(mct_rearr), intent(in),optional :: rearr ! rearranger for diff decomp !----- local ----- integer(IN) :: lsize ! lsize type(mct_aVect) :: avtri,avtro ! translated av on input/output grid character(*), parameter :: subname = '(shr_dmodel_translateAV) ' !------------------------------------------------------------------------------- ! PURPOSE: Fill avo from avi !------------------------------------------------------------------------------- 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) end subroutine shr_dmodel_translateAV_list !=============================================================================== end module shr_dmodel_mod