CESM_Driver: shr_dmodel_mod.F90 Source File

00001 module shr_dmodel_mod
00002 
00003 ! !USES:
00004 
00005   use shr_sys_mod
00006   use shr_kind_mod, only: IN=>SHR_KIND_IN, R8=>SHR_KIND_R8, &
00007                                CS=>SHR_KIND_CS, CL=>SHR_KIND_CL
00008   use shr_log_mod, only: loglev  => shr_log_Level
00009   use shr_log_mod, only: logunit => shr_log_Unit
00010   use shr_mpi_mod, only: shr_mpi_bcast
00011   use mct_mod
00012 !  use esmf_mod
00013   use perf_mod
00014   use pio
00015 
00016 ! !PUBLIC TYPES:
00017   implicit none
00018   private ! except
00019 
00020 !--------------------------------------------------------------------------
00021 ! Public interfaces
00022 !--------------------------------------------------------------------------
00023 
00024   public :: shr_dmodel_gsmapCreate
00025   public :: shr_dmodel_readLBUB
00026   public :: shr_dmodel_readgrid
00027   public :: shr_dmodel_gGridCompare
00028   public :: shr_dmodel_mapSet
00029   public :: shr_dmodel_translateAV
00030   public :: shr_dmodel_rearrGGrid
00031 
00032    interface shr_dmodel_mapSet; module procedure &
00033       shr_dmodel_mapSet_global
00034 !      shr_dmodel_mapSet_dest
00035    end interface
00036 
00037   integer(IN),parameter,public :: shr_dmodel_gGridCompareXYabs      = 1 ! X,Y  relative error
00038   integer(IN),parameter,public :: shr_dmodel_gGridCompareXYrel      = 2 ! X,Y  absolute error
00039   integer(IN),parameter,public :: shr_dmodel_gGridCompareAreaAbs    = 3 ! area relative error
00040   integer(IN),parameter,public :: shr_dmodel_gGridCompareAreaRel    = 4 ! area absolute error
00041   integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskIdent  = 5 ! masks are identical
00042   integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskZeros  = 6 ! masks have same zeros
00043   integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskSubset = 7 ! mask is subset of other
00044 
00045   ! masked methods
00046 
00047   integer(IN),parameter,public :: shr_dmodel_gGridCompareXYabsMask      = 101 ! X,Y  relative error
00048   integer(IN),parameter,public :: shr_dmodel_gGridCompareXYrelMask      = 102 ! X,Y  absolute error
00049   integer(IN),parameter,public :: shr_dmodel_gGridCompareAreaAbsMask    = 103 ! area relative error
00050   integer(IN),parameter,public :: shr_dmodel_gGridCompareAreaRelMask    = 104 ! area absolute error
00051   integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskIdentMask  = 105 ! masks are identical
00052   integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskZerosMask  = 106 ! masks have same zeros
00053   integer(IN),parameter,public :: shr_dmodel_gGridCompareMaskSubsetMask = 107 ! mask is subset of other
00054 
00055   integer(IN),parameter,public :: iotype_std_netcdf = -99 ! non pio option
00056 
00057 !--------------------------------------------------------------------------
00058 ! Private data
00059 !--------------------------------------------------------------------------
00060 
00061  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00062 CONTAINS
00063 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00064 
00065 !===============================================================================
00066 subroutine shr_dmodel_gsmapCreate(gsmap,gsize,compid,mpicom,decomp)
00067 
00068     implicit none
00069 
00070     type(mct_gsMap), intent(inout) :: gsmap
00071     integer(IN)    , intent(in)    :: gsize
00072     integer(IN)    , intent(in)    :: compid
00073     integer(IN)    , intent(in)    :: mpicom
00074     character(len=*),intent(in)    :: decomp
00075 
00076     ! local
00077 
00078     integer(IN) :: n,npes,ierr
00079     integer(IN), pointer :: start(:)     ! for gsmap initialization
00080     integer(IN), pointer :: length(:)    ! for gsmap initialization
00081     integer(IN), pointer :: pe_loc(:)    ! for gsmap initialization
00082     character(*), parameter :: subname = '(shr_dmodel_gsmapCreate) '
00083     character(*), parameter :: F00   = "('(shr_dmodel_gsmapCreate) ',8a)"
00084     character(*), parameter :: F01   = "('(shr_dmodel_gsmapCreate) ',a,5i8)"
00085 
00086     ! ---------------------------------------------
00087 
00088     if (gsize > 0) then
00089        call mpi_comm_size(mpicom,npes,ierr)
00090        allocate(start(npes),length(npes),pe_loc(npes))
00091 
00092        start = 0
00093        length = 0
00094        do n = 1,npes
00095           if (trim(decomp) == '1d') then
00096              length(n)  = gsize/npes
00097              if (n <= mod(gsize,npes)) length(n) = length(n) + 1
00098           elseif (trim(decomp) == 'root') then
00099              length = 0
00100              length(1) = gsize
00101           else
00102              write(logunit,F00) ' ERROR: decomp not allowed, ',trim(decomp)
00103              call shr_sys_abort(subname//' ERROR decomp')
00104           endif
00105           if (n == 1) then
00106              start(n) = 1
00107           else
00108              start(n) = start(n-1) + length(n-1)
00109           endif
00110           pe_loc(n) = n-1
00111        enddo
00112        call mct_gsMap_init( gsMap, COMPID, npes, gsize, start, length, pe_loc)
00113        deallocate(start,length,pe_loc)
00114     endif
00115 
00116 end subroutine shr_dmodel_gsmapCreate
00117 !===============================================================================
00118 
00119 subroutine shr_dmodel_readgrid( gGrid, gsMap, nxgo, nygo, filename, compid, mpicom, &
00120   decomp, lonname, latname, maskname, areaname, fracname, readfrac, &
00121   scmmode, scmlon, scmlat)
00122 
00123   use seq_flds_mod, only : seq_flds_dom_coord, seq_flds_dom_other
00124   use shr_file_mod, only : shr_file_noprefix, shr_file_queryprefix, shr_file_get
00125   use shr_string_mod, only : shr_string_lastindex
00126   use shr_ncread_mod, only : shr_ncread_domain, shr_ncread_vardimsizes, &
00127                              shr_ncread_varexists, shr_ncread_vardimnum
00128   implicit none
00129 
00130   !----- arguments -----
00131   type(mct_gGrid), intent(inout) :: gGrid
00132   type(mct_gsMap), intent(inout) :: gsMap
00133   integer(IN)    , intent(out)   :: nxgo
00134   integer(IN)    , intent(out)   :: nygo
00135   character(len=*),intent(in)    :: filename
00136   integer(IN)    , intent(in)    :: compid
00137   integer(IN)    , intent(in)    :: mpicom
00138   character(len=*),optional,intent(in)    ::   decomp ! decomp strategy for gsmap
00139   character(len=*),optional,intent(in)    ::  lonname ! name of  lon variable in file
00140   character(len=*),optional,intent(in)    ::  latname ! name of  lat variable in file
00141   character(len=*),optional,intent(in)    :: maskname ! name of mask variable in file
00142   character(len=*),optional,intent(in)    :: areaname ! name of area variable in file
00143   character(len=*),optional,intent(in)    :: fracname ! name of frac variable in file
00144   logical         ,optional,intent(in)    :: readfrac ! T <=> also read frac  in file
00145   logical         ,optional,intent(in)    :: scmmode  ! single column mode
00146   real(R8)        ,optional,intent(in)    :: scmlon   ! single column lon
00147   real(R8)        ,optional,intent(in)    :: scmlat   ! single column lat
00148 
00149   !----- local -----
00150   integer(IN)   :: n,k,j,i    ! indices  
00151   integer(IN)   :: lsize      ! lsize
00152   integer(IN)   :: gsize      ! gsize
00153   integer(IN)   :: my_task, master_task
00154   integer(IN)   :: ierr       ! error code
00155   logical       :: fileexists ! 
00156   integer(IN)   :: rCode      ! return code
00157   character(CL) :: remoteFn   ! input file name (possibly at an archival location)
00158   character(CL) :: localFn    ! file name to be opened (possibly a local copy)
00159   character(CS) :: prefix     ! file prefix
00160   character(CS) ::   ldecomp  ! decomp strategy
00161   character(CS) ::  llatname  ! name of  lat variable
00162   character(CS) ::  llonname  ! name of  lon variable
00163   character(CS) :: lmaskname  ! name of mask variable
00164   character(CS) :: lareaname  ! name of area variable
00165   character(CS) :: lfracname  ! name of area variable
00166   logical       :: lreadfrac  ! read fraction
00167   logical       :: maskexists ! is mask on dataset
00168   integer(IN)   :: nxg,nyg    ! size of input fields
00169   integer(IN)   :: ndims      ! number of dims
00170   integer(IN)   :: nlon,nlat,narea,nmask,nfrac
00171   logical       :: lscmmode   ! local scm mode
00172   real(R8)      :: dist,mind  ! scmmode point search
00173   integer(IN)   :: ni,nj      ! scmmode point search
00174   real(R8)      :: lscmlon    ! local copy of scmlon
00175 
00176   real   (R8),allocatable ::  lon(:,:) ! temp array for domain lon  info
00177   real   (R8),allocatable ::  lat(:,:) ! temp array for domain lat  info
00178   integer(IN),allocatable :: mask(:,:) ! temp array for domain mask info
00179   real   (R8),allocatable :: area(:,:) ! temp array for domain area info
00180   real   (R8),allocatable :: frac(:,:) ! temp array for domain frac info
00181 
00182   integer(IN), pointer :: idata(:)   ! temporary
00183   type(mct_ggrid)      :: gGridRoot       ! global mct ggrid
00184 
00185   character(*), parameter :: subname = '(shr_dmodel_readgrid) '
00186   character(*), parameter :: F00   = "('(shr_dmodel_readgrid) ',8a)"
00187   character(*), parameter :: F01   = "('(shr_dmodel_readgrid) ',a,5i8)"
00188 
00189 !-------------------------------------------------------------------------------
00190 ! PURPOSE:  Read MCT ggrid and set gsmap
00191 !----------------------------------------------------------------------------
00192 ! Notes:
00193 ! o as per shr_file_get(), the file name format is expected to be
00194 !   remoteFn = [location:][directory path]localFn
00195 !   eg. "foobar.nc"  "/home/user/foobar.nc"  "mss:/USER/fobar.nc"
00196 ! o assumes a very specific netCDF domain file format wrt var names, etc.
00197 !
00198 ! TO DO: have the calling routine select/input the domain's file name
00199 !----------------------------------------------------------------------------
00200 
00201   call MPI_COMM_RANK(mpicom,my_task,ierr)
00202   master_task = 0
00203 
00204   lscmmode = .false.
00205   if (present(scmmode)) then
00206      lscmmode = scmmode
00207      if (lscmmode) then
00208         if (.not.present(scmlon) .or. .not.present(scmlat)) then
00209            write(logunit,*) subname,' ERROR: scmmode must supply scmlon and scmlat'
00210            call shr_sys_abort(subname//' ERROR: scmmode1 lon lat')
00211         endif
00212         if (my_task > 0) then
00213            write(logunit,*) subname,' ERROR: scmmode must be run on one pe'
00214            call shr_sys_abort(subname//' ERROR: scmmode2 tasks')
00215         endif
00216      endif
00217   endif
00218 
00219   if (my_task == master_task) then
00220      if ( shr_file_queryPrefix(fileName,prefix=prefix) /= shr_file_noPrefix ) then
00221         n        = max(len_trim(prefix),shr_string_lastIndex(fileName,"/"))
00222         remoteFn = fileName
00223         localFn  = fileName(n+1: len_trim(fileName) )
00224         call shr_file_get(rCode,localFn,remoteFn)
00225      else
00226         remoteFn = "undefined" ! this isn't needed
00227         localFn  = fileName    ! file to open
00228      end if
00229      inquire(file=trim(localFn),exist=fileExists)
00230      if (.not. fileExists) then
00231          write(logunit,F00) "ERROR: file does not exist: ", trim(localFn)
00232          call shr_sys_abort(subName//"ERROR: file does not exist: "//trim(localFn))
00233      end if
00234   endif
00235 
00236   lreadfrac = .false.
00237   ldecomp   = "1d"
00238   llonname  = "xc"  ! default values / standard data model domain file format
00239   llatname  = "yc"
00240   lmaskname = "mask"
00241   lareaname = "area"
00242   lfracname = "frac"
00243   if (present(  decomp))   ldecomp =   decomp
00244   if (present(readfrac)) lreadfrac = readfrac
00245   if (present( lonname))  llonname =  lonname
00246   if (present( latname))  llatname =  latname
00247   if (present(maskname)) lmaskname = maskname
00248   if (present(areaname)) lareaname = areaname
00249   if (present(fracname)) lfracname = fracname
00250 
00251   ! Initialize mct domain type
00252   ! lat/lon in degrees,  area in radians^2, mask is 1 (ocean), 0 (non-ocean)
00253 
00254   if (my_task == master_task) then
00255      if (shr_ncread_varexists(localFn,lmaskname)) then
00256         maskexists = .true.
00257         call shr_ncread_varDimSizes(localFn,lmaskname,nxg,nyg)
00258      else
00259         maskexists = .false.
00260         call shr_ncread_varDimNum(localFn,llonName,ndims)
00261         if (ndims == 1) then
00262            call shr_ncread_varDimSizes(localFn,llonName,nxg)
00263            call shr_ncread_varDimSizes(localFn,llatName,nyg)
00264         else
00265            call shr_ncread_varDimSizes(localFn,llonName,nxg,nyg)
00266         endif
00267      endif
00268   endif
00269   call shr_mpi_bcast(nxg,mpicom)
00270   call shr_mpi_bcast(nyg,mpicom)
00271   if (lscmmode) then
00272      nxgo = 1
00273      nygo = 1
00274      gsize = 1
00275   else
00276      nxgo = nxg
00277      nygo = nyg
00278      gsize = nxg*nyg
00279      if (gsize < 1) return
00280   endif
00281 
00282   call shr_dmodel_gsmapCreate(gsMap,gsize,compid,mpicom,trim(ldecomp))
00283   lsize = mct_gsMap_lsize(gsMap, mpicom)
00284   call mct_gGrid_init( GGrid=gGrid, CoordChars=trim(seq_flds_dom_coord), &
00285      OtherChars=trim(seq_flds_dom_other), lsize=lsize )
00286 
00287   ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid
00288 
00289   call mct_gsMap_orderedPoints(gsMap, my_task, idata)
00290   call mct_gGrid_importIAttr(gGrid,'GlobGridNum',idata,lsize)
00291   deallocate(idata)
00292 
00293   ! Initialize attribute vector with special value
00294 
00295   gGrid%data%rAttr = -9999.0_R8
00296 
00297   ! Load file data into domG then scatter to ggrid
00298 
00299   if (my_task == master_task) then
00300 
00301      allocate(lon(nxg,nyg))
00302      allocate(lat(nxg,nyg))
00303      allocate(area(nxg,nyg))
00304      allocate(mask(nxg,nyg))
00305      allocate(frac(nxg,nyg))
00306 
00307      if (.not.maskexists) then
00308         call shr_ncread_domain(localFn,llonName,lon,llatName,lat)
00309         mask = 1
00310         frac = 1.0_R8
00311         area = 1.0e36_R8
00312      else
00313         if (lreadfrac) then
00314            call shr_ncread_domain(localFn,llonName,lon,llatName,lat, &
00315                 lmaskName,mask,lareaName,area,lfracName,frac)
00316         else ! assume frac = 1.0
00317            call shr_ncread_domain(localFn,llonName,lon,llatName,lat, &
00318                 lmaskName,mask,lareaName,area)
00319            where (mask == 0)
00320               frac = 0.0_R8
00321            elsewhere
00322               frac = 1.0_R8
00323            end where
00324         endif
00325      endif
00326 
00327      call mct_gGrid_init(gGridRoot,gGrid,gsize)
00328 
00329 ! initialize gGridRoot to avoid errors when using strict compiler checks
00330      gGridRoot%data%rAttr = -9999.0_R8
00331 
00332      nlon  = mct_aVect_indexRA(gGridRoot%data,'lon')
00333      nlat  = mct_aVect_indexRA(gGridRoot%data,'lat')
00334      narea = mct_aVect_indexRA(gGridRoot%data,'area')
00335      nmask = mct_aVect_indexRA(gGridRoot%data,'mask')
00336      nfrac = mct_aVect_indexRA(gGridRoot%data,'frac')
00337 
00338      if (lscmmode) then
00339         !--- assumes regular 2d grid for compatability with shr_scam_getCloseLatLon ---
00340         !--- want lon values between 0 and 360, assume 1440 is enough ---
00341         lscmlon = mod(scmlon+1440.0_r8,360.0_r8)
00342         lon     = mod(lon   +1440.0_r8,360.0_r8)
00343 
00344         !--- start with wraparound ---
00345         ni = 1
00346         mind = abs(lscmlon - (lon(1,1)+360.0_r8))
00347         do i=1,nxg
00348            dist = abs(lscmlon - lon(i,1))
00349            if (dist < mind) then
00350               mind = dist
00351               ni = i
00352            endif
00353         enddo        
00354 
00355         nj = -1
00356         mind = 1.0e20
00357         do j=1,nyg
00358            dist = abs(scmlat - lat(1,j))
00359            if (dist < mind) then
00360               mind = dist
00361               nj = j
00362            endif
00363         enddo
00364 
00365         n = 1
00366         i = ni
00367         j = nj
00368         gGridRoot%data%rAttr(nlat ,n) = lat(i,j)
00369         gGridRoot%data%rAttr(nlon ,n) = lon(i,j)
00370         gGridRoot%data%rAttr(narea,n) = area(i,j)
00371         gGridRoot%data%rAttr(nmask,n) = real(mask(i,j),R8)
00372         gGridRoot%data%rAttr(nfrac,n) = frac(i,j)
00373      else
00374         n=0                                                                                    
00375         do j=1,nyg
00376         do i=1,nxg
00377            n=n+1                                                                               
00378            gGridRoot%data%rAttr(nlat ,n) = lat(i,j)
00379            gGridRoot%data%rAttr(nlon ,n) = lon(i,j)
00380            gGridRoot%data%rAttr(narea,n) = area(i,j)
00381            gGridRoot%data%rAttr(nmask,n) = real(mask(i,j),R8)
00382            gGridRoot%data%rAttr(nfrac,n) = frac(i,j)
00383         enddo
00384         enddo
00385      endif
00386 
00387   endif
00388 
00389   call mct_gGrid_scatter(gGridRoot, gGrid, gsMap, master_task, mpicom)
00390   if (my_task == master_task) call mct_gGrid_clean(gGridRoot)
00391 
00392 end subroutine shr_dmodel_readgrid
00393 
00394 !===============================================================================
00395 
00396 subroutine shr_dmodel_readLBUB(stream,pio_subsystem,pio_iotype,pio_iodesc,mDate,mSec,mpicom,gsMap, &
00397                              avLB,mDateLB,mSecLB,avUB,mDateUB,mSecUB,newData,rmOldFile,istr)
00398 
00399   use shr_file_mod, only : shr_file_noprefix, shr_file_queryprefix, shr_file_get
00400   use shr_stream_mod
00401   implicit none
00402 
00403   !----- arguments -----
00404   type(shr_stream_streamType),intent(inout) :: stream
00405   type(iosystem_desc_t)      ,intent(inout) :: pio_subsystem
00406   integer(IN)                ,intent(in)    :: pio_iotype
00407   type(io_desc_t)            ,intent(inout) :: pio_iodesc
00408   integer(IN)                ,intent(in)    :: mDate  ,mSec
00409   integer(IN)                ,intent(in)    :: mpicom
00410   type(mct_gsMap)            ,intent(in)    :: gsMap
00411   type(mct_aVect)            ,intent(inout) :: avLB
00412   integer(IN)                ,intent(inout) :: mDateLB,mSecLB
00413   type(mct_aVect)            ,intent(inout) :: avUB
00414   integer(IN)                ,intent(inout) :: mDateUB,mSecUB
00415   logical                    ,intent(out)   :: newData
00416   logical,optional           ,intent(in)    :: rmOldFile
00417   character(len=*),optional  ,intent(in)    :: istr
00418 
00419   !----- local -----
00420   integer(IN)   :: n,k,j,i    ! indices  
00421   integer(IN)   :: lsize      ! lsize
00422   integer(IN)   :: gsize      ! gsize
00423   integer(IN)   :: my_task, master_task
00424   integer(IN)   :: ierr       ! error code
00425   integer(IN)   :: rCode      ! return code
00426   logical       :: localCopy,fileexists
00427   integer(IN)   :: ivals(6)   ! bcast buffer
00428 
00429   integer(IN)   :: oDateLB,oSecLB,dDateLB,oDateUB,oSecUB,dDateUB
00430   real(R8)      :: rDateM,rDateLB,rDateUB  ! model,LB,UB dates with fractional days
00431   integer(IN)   ::  n_lb, n_ub
00432   character(CL) :: fn_lb,fn_ub,fn_next,fn_prev
00433   character(CL) :: path
00434   character(len=32) :: lstr
00435 
00436   real(R8)      :: spd
00437 
00438   character(*), parameter :: subname = '(shr_dmodel_readLBUB) '
00439   character(*), parameter :: F00   = "('(shr_dmodel_readLBUB) ',8a)"
00440   character(*), parameter :: F01   = "('(shr_dmodel_readLBUB) ',a,5i8)"
00441   character(*), parameter :: F02   = "('(shr_dmodel_readLBUB) ',3a,i8)"
00442 
00443 !-------------------------------------------------------------------------------
00444 ! PURPOSE:  Read LB and UB stream data
00445 !----------------------------------------------------------------------------
00446 
00447   lstr = 'shr_dmodel_readLBUB'
00448   if (present(istr)) then
00449     lstr = trim(istr)
00450   endif
00451 
00452   call t_startf(trim(lstr)//'_setup')
00453   call MPI_COMM_RANK(mpicom,my_task,ierr)
00454   master_task = 0
00455   spd = 86400.0_R8
00456 
00457   newData = .false.
00458   n_lb = -1
00459   n_ub = -1
00460   fn_lb = 'undefinedlb'
00461   fn_ub = 'undefinedub'
00462 
00463   oDateLB = mDateLB
00464   oSecLB  = mSecLB
00465   oDateUB = mDateUB
00466   oSecUB  = mSecUB
00467 
00468   ! for backwards compatability (yikes)
00469   rDateM  = real(mDate  ,R8) + real(mSec  ,R8)/spd
00470   rDateLB = real(mDateLB,R8) + real(mSecLB,R8)/spd
00471   rDateUB = real(mDateUB,R8) + real(mSecUB,R8)/spd
00472   call t_stopf(trim(lstr)//'_setup')
00473 
00474   if (rDateM < rDateLB .or. rDateM > rDateUB) then
00475      call t_startf(trim(lstr)//'_fbound')
00476      if (my_task == master_task) then
00477 !        call shr_stream_findBounds(stream,mDate,mSec,                 &
00478 !                                   mDateLB,dDateLB,mSecLB,n_lb,fn_lb, &
00479 !                                   mDateUB,dDateUB,mSecUB,n_ub,fn_ub  )
00480         call shr_stream_findBounds(stream,mDate,mSec,                 &
00481                                    ivals(1),dDateLB,ivals(2),ivals(5),fn_lb, &
00482                                    ivals(3),dDateUB,ivals(4),ivals(6),fn_ub  )
00483         call shr_stream_getFilePath(stream,path)
00484         localCopy = (shr_file_queryPrefix(path) /= shr_file_noPrefix )
00485 !        write(logunit,*) 'tcx ivals ',ivals
00486      endif
00487      call t_stopf(trim(lstr)//'_fbound')
00488      call t_startf(trim(lstr)//'_bcast')
00489 
00490 !    --- change 4 bcasts to a single bcast and copy for performance ---
00491 !     call shr_mpi_bcast(mDateLB,mpicom)
00492 !     call shr_mpi_bcast(mSecLB,mpicom)
00493 !     call shr_mpi_bcast(mDateUB,mpicom)
00494 !     call shr_mpi_bcast(mSecUB,mpicom)
00495      call shr_mpi_bcast(ivals,mpicom)
00496      mDateLB = ivals(1)
00497      mSecLB  = ivals(2)
00498      mDateUB = ivals(3)
00499      mSecUB  = ivals(4)
00500      n_lb    = ivals(5)
00501      n_ub    = ivals(6)
00502      call t_stopf(trim(lstr)//'_bcast')
00503   endif
00504 
00505   if (mDateLB /= oDateLB .or. mSecLB /= oSecLB) then
00506      newdata = .true.
00507      if (mDateLB == oDateUB .and. mSecLB == oSecUB) then
00508         call t_startf(trim(lstr)//'_LB_copy')
00509         avLB%rAttr(:,:) = avUB%rAttr(:,:)
00510         call t_stopf(trim(lstr)//'_LB_copy')
00511      else
00512         if (my_task == master_task) write(logunit,F02) 'reading file: ',trim(path),trim(fn_lb),n_lb
00513         call shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gsMap, avLB, mpicom, &
00514              path, fn_lb, n_lb,istr=trim(lstr)//'_LB')
00515      endif
00516   endif
00517 
00518   if (mDateUB /= oDateUB .or. mSecUB /= oSecUB) then
00519      newdata = .true.
00520      if (my_task == master_task) write(logunit,F02) 'reading file: ',trim(path),trim(fn_ub),n_ub
00521      call shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gsMap, avUB, mpicom, &
00522           path, fn_ub, n_ub,istr=trim(lstr)//'_UB')
00523   endif
00524 
00525   call t_startf(trim(lstr)//'_filemgt')
00526   !--- determine previous & next data files in list of files ---
00527   if (my_task == master_task .and. newdata) then
00528   call shr_stream_getFilePath(stream,path)
00529   localCopy = (shr_file_queryPrefix(path) /= shr_file_noPrefix )
00530   if (localCopy) then
00531      call shr_stream_getPrevFileName(stream,fn_lb,fn_prev,path)
00532      call shr_stream_getNextFileName(stream,fn_ub,fn_next,path)
00533      inquire(file=trim(fn_next),exist=fileExists)
00534      if ( trim(fn_next) == "unknown" .or. fileExists) then
00535         ! do nothing
00536      else 
00537         call shr_file_get(rCode,fn_next,trim(path)//fn_next,async=.true.)
00538         write(logunit,F00) "get next file: ",trim(fn_next)
00539      end if
00540 
00541      !---  remove the old file? (only if acquiring local copies) ---
00542      if ( rmOldFile .and. &
00543           fn_prev/=fn_lb .and. fn_prev/=fn_ub .and. fn_prev/=fn_next ) then 
00544         !--- previous file is not in use and is not next in list ---
00545         inquire(file=trim(fn_prev),exist=fileExists)
00546         if ( fileExists ) then
00547            call shr_sys_system(" rm "//trim(fn_prev),rCode)
00548            write(logunit,F00) "rm  prev file: ",trim(fn_prev)
00549         end if
00550      end if
00551   endif
00552   endif
00553   call t_stopf(trim(lstr)//'_filemgt')
00554 
00555 end subroutine shr_dmodel_readLBUB
00556 
00557 !===============================================================================
00558 subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gsMap, av, mpicom, &
00559            path, fn, nt, istr)
00560 
00561   use shr_file_mod, only : shr_file_noprefix, shr_file_queryprefix, shr_file_get
00562   use shr_stream_mod
00563   use shr_ncread_mod
00564   implicit none
00565 
00566   !----- arguments -----
00567   type(shr_stream_streamType),intent(inout) :: stream
00568   type(iosystem_desc_t),intent(inout)  :: pio_subsystem
00569   integer(IN)     ,intent(in)    :: pio_iotype
00570   type(io_desc_t)      ,intent(inout)  :: pio_iodesc
00571   type(mct_gsMap) ,intent(in)    :: gsMap
00572   type(mct_aVect) ,intent(inout) :: av  
00573   integer(IN)     ,intent(in)    :: mpicom
00574   character(len=*),intent(in)    :: path
00575   character(len=*),intent(in)    :: fn
00576   integer(IN)     ,intent(in)    :: nt
00577   character(len=*),optional  ,intent(in)    :: istr
00578 
00579   !----- local -----
00580   integer(IN) :: my_task
00581   integer(IN) :: master_task
00582   integer(IN) :: ierr
00583   logical     :: localCopy,fileexists
00584   type(mct_avect) :: avG
00585   integer(IN) :: gsize,nx,ny
00586   integer(IN) :: k
00587   integer(IN) :: fid
00588   integer(IN) :: rCode      ! return code
00589   real(R8),allocatable :: data(:,:)
00590   character(CL) :: fileName
00591   character(CL) :: sfldName
00592   type(mct_avect) :: avtmp
00593   character(len=32) :: lstr
00594 
00595   integer(in) :: ndims
00596   integer(in),pointer :: dimid(:)
00597   type(file_desc_t) :: pioid
00598   type(var_desc_t) :: varid
00599   integer(kind=pio_offset) :: frame
00600 
00601   character(*), parameter :: subname = '(shr_dmodel_readstrm) '
00602   character(*), parameter :: F00   = "('(shr_dmodel_readstrm) ',8a)"
00603   character(*), parameter :: F01   = "('(shr_dmodel_readstrm) ',a,5i8)"
00604 
00605 !-------------------------------------------------------------------------------
00606 
00607   lstr = 'shr_dmodel_readstrm'
00608   if (present(istr)) then
00609     lstr = trim(istr)
00610   endif
00611 
00612   call t_barrierf(trim(lstr)//'_BARRIER',mpicom)
00613   call t_startf(trim(lstr)//'_setup')
00614   call MPI_COMM_RANK(mpicom,my_task,ierr)
00615   master_task = 0
00616 
00617   gsize = mct_gsmap_gsize(gsMap)
00618 
00619   if (my_task == master_task) then
00620      localCopy = (shr_file_queryPrefix(path) /= shr_file_noPrefix )
00621      if (localCopy) then
00622         call shr_file_get(rCode,fn,trim(path)//fn)
00623         fileName = fn
00624      else                 ! DON'T acquire a local copy of the data file
00625         fileName = trim(path)//fn
00626      end if
00627      inquire(file=trim(fileName),exist=fileExists)
00628      if (.not. fileExists) then
00629         write(logunit,F00) "ERROR: file does not exist: ", trim(fileName)
00630         call shr_sys_abort(subName//"ERROR: file does not exist: "//trim(fileName))
00631      end if
00632   endif
00633 
00634   if (my_task == master_task) then
00635      call shr_stream_getFileFieldName(stream,1,sfldName)
00636   endif
00637 
00638   call t_stopf(trim(lstr)//'_setup')
00639 
00640   if (pio_iotype == iotype_std_netcdf) then
00641 
00642      call t_startf(trim(lstr)//'_readcdf')
00643      if (my_task == master_task) then
00644         call shr_ncread_varDimSizes(trim(fileName),trim(sfldName),nx,ny)
00645         if (gsize /= nx*ny) then
00646            write(logunit,F01) "ERROR in data sizes ",nx,ny,gsize
00647            call shr_sys_abort(subname//"ERROR in data sizes")
00648         endif
00649         call mct_aVect_init(avG,av,gsize)
00650         allocate(data(nx,ny))
00651         call shr_ncread_open(trim(fileName),fid,rCode)
00652         do k = 1,mct_aVect_nRAttr(av)
00653            call shr_stream_getFileFieldName(stream,k,sfldName)
00654            call shr_ncread_tField(fileName,nt,sfldName,data,fidi=fid,rc=rCode)
00655            avG%rAttr(k,:) = reshape(data, (/gsize/))
00656         enddo
00657         call shr_ncread_close(fid,rCode)
00658         deallocate(data)
00659      endif
00660      call t_stopf(trim(lstr)//'_readcdf')
00661      call t_barrierf(trim(lstr)//'_scatter'//'_BARRIER',mpicom)
00662      call t_startf(trim(lstr)//'_scatter')
00663      call mct_aVect_scatter(avG,avtmp,gsMap,master_task,mpicom)
00664      call mct_aVect_copy(avtmp,av)
00665      if (my_task == master_task) call mct_aVect_clean(avG)
00666      call mct_aVect_clean(avtmp)
00667      call t_stopf(trim(lstr)//'_scatter')
00668 
00669   else
00670 
00671      call t_startf(trim(lstr)//'_readpio')
00672      call shr_mpi_bcast(sfldName,mpicom,'sfldName')
00673      call shr_mpi_bcast(filename,mpicom,'filename')
00674      rcode = pio_openfile(pio_subsystem, pioid, pio_iotype, trim(filename), pio_nowrite)
00675      call pio_seterrorhandling(pioid,PIO_INTERNAL_ERROR)
00676 
00677      rcode = pio_inq_varid(pioid,trim(sfldName),varid)
00678      rcode = pio_inq_varndims(pioid, varid, ndims)
00679      allocate(dimid(ndims))
00680      rcode = pio_inq_vardimid(pioid, varid, dimid(1:ndims))
00681      rcode = pio_inq_dimlen(pioid, dimid(1), nx)
00682      rcode = pio_inq_dimlen(pioid, dimid(2), ny)
00683      deallocate(dimid)
00684      if (gsize /= nx*ny) then
00685         write(logunit,F01) "ERROR in data sizes ",nx,ny,gsize
00686         call shr_sys_abort(subname//"ERROR in data sizes")
00687      endif
00688 
00689      do k = 1,mct_aVect_nRAttr(av)
00690         if (my_task == master_task) then
00691            call shr_stream_getFileFieldName(stream,k,sfldName)
00692         endif
00693         call shr_mpi_bcast(sfldName,mpicom,'sfldName')
00694         rcode = pio_inq_varid(pioid,trim(sfldName),varid)
00695         frame = nt
00696         call pio_setframe(varid,frame)
00697         call pio_read_darray(pioid,varid,pio_iodesc,av%rattr(k,:),rcode)
00698      enddo
00699 
00700      call pio_closefile(pioid)
00701      call t_stopf(trim(lstr)//'_readpio')
00702 
00703   endif
00704 
00705 end subroutine shr_dmodel_readstrm
00706 !===============================================================================
00707 !===============================================================================
00708 !BOP ===========================================================================
00709 !
00710 ! !IROUTINE: shr_dmodel_gGridCompare -- returns TRUE if two domains are the same.
00711 !
00712 ! !DESCRIPTION:
00713 !    Returns TRUE if two domains are the the same (within tolerance).
00714 !
00715 ! !REVISION HISTORY:
00716 !    2005-May-03 - B. Kauffman - added mulitiple methods
00717 !    2005-May-27 - T. Craig - first version
00718 !
00719 ! !INTERFACE: ------------------------------------------------------------------
00720 
00721 logical function shr_dmodel_gGridCompare(ggrid1,gsmap1,ggrid2,gsmap2,method,mpicom,eps)
00722 
00723 ! !USES:
00724 
00725    implicit none
00726 
00727 ! !INPUT/OUTPUT PARAMETERS:
00728 
00729    type(mct_gGrid)     ,intent(in)  :: ggrid1   ! 1st ggrid
00730    type(mct_gsmap)     ,intent(in)  :: gsmap1   ! 1st gsmap
00731    type(mct_gGrid)     ,intent(in)  :: ggrid2   ! 2nd ggrid
00732    type(mct_gsmap)     ,intent(in)  :: gsmap2   ! 2nd gsmap
00733    integer(IN)         ,intent(in)  :: method   ! selects what to compare
00734    integer(IN)         ,intent(in)  :: mpicom   ! mpicom
00735    real(R8)   ,optional,intent(in)  :: eps      ! epsilon compare value
00736 
00737 !EOP
00738 
00739    !--- local ---
00740    real(R8)    :: leps         ! local epsilon
00741    integer(IN) :: n            ! counters
00742    integer(IN) :: my_task,master_task
00743    integer(IN) :: gsize
00744    integer(IN) :: ierr
00745    integer(IN) :: nlon1, nlon2, nlat1, nlat2, nmask1, nmask2  ! av field indices
00746    logical     :: compare      ! local compare logical
00747    real(R8)    :: lon1,lon2    ! longitudes to compare
00748    real(R8)    :: lat1,lat2    ! latitudes to compare
00749    real(R8)    :: msk1,msk2    ! masks to compare
00750    integer(IN) :: nx,ni1,ni2   ! i grid size, i offset for 1 vs 2 and 2 vs 1
00751    integer(IN) :: n1,n2,i,j    ! local indices
00752    type(mct_aVect) :: avG1     ! global av
00753    type(mct_aVect) :: avG2     ! global av
00754    integer(IN) :: lmethod      ! local method
00755    logical     :: maskmethod, maskpoint ! masking on method
00756  
00757    !--- formats ---
00758    character(*),parameter :: subName = '(shr_dmodel_gGridCompare) '
00759    character(*),parameter :: F01     = "('(shr_dmodel_gGridCompare) ',4a)" 
00760 
00761 !-------------------------------------------------------------------------------
00762 !
00763 !-------------------------------------------------------------------------------
00764 
00765    call MPI_COMM_RANK(mpicom,my_task,ierr)
00766    master_task = 0
00767 
00768    leps = 1.0e-6_R8
00769    if (present(eps)) leps = eps
00770 
00771    lmethod = mod(method,100)
00772    if (method > 100) then
00773       maskmethod=.true.
00774    else
00775       maskmethod=.false.
00776    endif
00777 
00778    call mct_aVect_gather(gGrid1%data,avG1,gsmap1,master_task,mpicom)
00779    call mct_aVect_gather(gGrid2%data,avG2,gsmap2,master_task,mpicom)
00780 
00781    if (my_task == master_task) then
00782 
00783       compare = .true. 
00784       gsize = mct_aVect_lsize(avG1)
00785       if (gsize /= mct_aVect_lsize(avG2)) then
00786          compare = .false.
00787       endif
00788 
00789       if (.not. compare ) then
00790          !--- already failed the comparison test, check no futher ---
00791       else 
00792          nlon1 = mct_aVect_indexRA(avG1,'lon')
00793          nlat1 = mct_aVect_indexRA(avG1,'lat')
00794          nlon2 = mct_aVect_indexRA(avG2,'lon')
00795          nlat2 = mct_aVect_indexRA(avG2,'lat')
00796          nmask1 = mct_aVect_indexRA(avG1,'mask')
00797          nmask2 = mct_aVect_indexRA(avG2,'mask')
00798 
00799          ! To compare, want to be able to treat longitude wraparound generally.
00800          ! So we need to compute i index offset and we need to compute the size of the nx dimension
00801          ! First adjust the lon so it's in the range [0,360), add 1440 to lon to take into
00802          ! accounts lons less than 1440.  if any lon is less than -1440, abort.  1440 is arbitrary
00803          ! Next, comute ni1 and ni2.  These are the offsets of grid1 relative to grid2 and
00804          ! grid2 relative to grid1.  The sum of those offsets is nx.  Use ni1 to offset grid2
00805          ! in comparison and compute new grid2 index from ni1 and nx.  If ni1 is zero, then
00806          ! there is no offset, don't need to compute ni2, and nx can be anything > 0.
00807 
00808          !--- compute offset of grid2 compared to pt 1 of grid 1
00809          lon1 = minval(avG1%rAttr(nlon1,:))
00810          lon2 = minval(avG2%rAttr(nlon2,:))
00811          if ((lon1 < -1440.0_R8) .or. (lon2 < -1440.0_R8)) then
00812             write(logunit,*) subname,' ERROR: lon1 lon2 lt -1440 ',lon1,lon2
00813             call shr_sys_abort(subname//' ERROR: lon1 lon2 lt -1440')
00814          endif
00815 
00816          lon1 = mod(avG1%rAttr(nlon1,1)+1440.0_R8,360.0_R8)
00817          lat1 = avG1%rAttr(nlat1,1)
00818          ni1 = -1
00819          do n = 1,gsize
00820             lon2 = mod(avG2%rAttr(nlon2,n)+1440.0_R8,360.0_R8)
00821             lat2 = avG2%rAttr(nlat2,n)
00822             if ((ni1 < 0) .and. abs(lon1-lon2) <= leps .and. abs(lat1-lat2) <= leps) then
00823                ni1 = n - 1  ! offset, compare to first gridcell in grid 1
00824             endif
00825          enddo
00826 
00827          if (ni1 < 0) then        ! no match for grid point 1, so fails without going further
00828             compare = .false.
00829          elseif (ni1 == 0) then   ! no offset, set nx to anything > 0
00830             nx = 1
00831          else                     ! now compute ni2
00832             ! compute offset of grid1 compared to pt 1 of grid 2
00833             lon2 = mod(avG2%rAttr(nlon2,1)+1440.0_R8,360.0_R8)
00834             lat2 = avG2%rAttr(nlat2,1)
00835             ni2 = -1
00836             do n = 1,gsize
00837                lon1 = mod(avG1%rAttr(nlon1,n)+1440.0_R8,360.0_R8)
00838                lat1 = avG1%rAttr(nlat1,n)
00839                if ((ni2 < 0) .and. abs(lon1-lon2) <= leps .and. abs(lat1-lat2) <= leps) then
00840                   ni2 = n - 1  ! offset, compare to first gridcell in grid 1
00841                endif
00842             enddo
00843             if (ni2 < 0) then
00844                write(logunit,*) subname,' ERROR in ni2 ',ni1,ni2
00845                call shr_sys_abort(subname//' ERROR in ni2')
00846             endif
00847             nx = ni1 + ni2
00848          endif
00849 
00850          if (compare) then
00851             do n = 1,gsize
00852                j = ((n-1)/nx) + 1
00853                i = mod(n-1,nx) + 1
00854                n1 = (j-1)*nx + mod(n-1,nx) + 1
00855                n2 = (j-1)*nx + mod(n-1+ni1,nx) + 1
00856                if (n1 /= n) then    ! sanity check, could be commented out
00857                   write(logunit,*) subname,' ERROR in n1 n2 ',n,i,j,n1,n2
00858                   call shr_sys_abort(subname//' ERROR in n1 n2')
00859                endif
00860                lon1 = mod(avG1%rAttr(nlon1,n1)+1440.0_R8,360.0_R8)
00861                lat1 = avG1%rAttr(nlat1,n1)
00862                lon2 = mod(avG2%rAttr(nlon2,n2)+1440.0_R8,360.0_R8)
00863                lat2 = avG2%rAttr(nlat2,n2)
00864                msk1 = avG1%rAttr(nmask1,n1)
00865                msk2 = avG2%rAttr(nmask2,n2)
00866 
00867                maskpoint = .true.
00868                if (maskmethod .and. (msk1 == 0 .or. msk2 == 0)) then
00869                   maskpoint = .false.
00870                endif
00871 
00872                if (maskpoint) then
00873                if (lmethod == shr_dmodel_gGridCompareXYabs      ) then
00874                   if (abs(lon1 - lon2) > leps) compare = .false.
00875                   if (abs(lat1 - lat2) > leps) compare = .false.
00876                else if (lmethod == shr_dmodel_gGridCompareXYrel      ) then
00877                   if (rdiff(lon1,lon2) > leps) compare = .false.
00878                   if (rdiff(lat1,lat2) > leps) compare = .false.
00879                else if (lmethod == shr_dmodel_gGridCompareMaskIdent  ) then
00880                   if (msk1 /= msk2)compare = .false.
00881                else if (lmethod == shr_dmodel_gGridCompareMaskZeros  ) then
00882                   if (msk1 == 0 .and. msk2 /= 0) compare = .false.
00883                   if (msk1 /= 0 .and. msk2 == 0) compare = .false.
00884                else if (lmethod == shr_dmodel_gGridCompareMaskSubset ) then
00885                   if (msk1 /= 0 .and. msk2 == 0) compare = .false.
00886                else
00887                   write(logunit,F01) "ERROR: compare method not recognized, method = ",method
00888                   call shr_sys_abort(subName//"ERROR: compare method not recognized")
00889                endif  ! lmethod
00890                endif  ! maskpoint
00891             enddo ! gsize
00892          endif  ! compare
00893       endif   ! compare
00894    endif   ! master_task
00895 
00896    if (my_task == master_task) call mct_avect_clean(avG1)
00897    if (my_task == master_task) call mct_avect_clean(avG2)
00898 
00899    call shr_mpi_bcast(compare,mpicom)
00900    shr_dmodel_gGridCompare = compare
00901 
00902    return
00903 
00904 !-------------------------------------------------------------------------------
00905 contains   ! internal subprogram
00906 !-------------------------------------------------------------------------------
00907 
00908    real(R8) function rdiff(v1,v2) ! internal function
00909       !------------------------------------------
00910       real(R8),intent(in) :: v1,v2                 ! two values to compare
00911       real(R8),parameter  :: c0           = 0.0_R8 ! zero
00912       real(R8),parameter  :: large_number = 1.0e20_R8 ! infinity
00913       !------------------------------------------
00914       if (v1 == v2) then
00915          rdiff = c0
00916       elseif (v1 == c0 .and. v2 /= c0) then
00917          rdiff = large_number
00918       elseif (v2 == c0 .and. v1 /= c0) then
00919          rdiff = large_number
00920       else
00921 !        rdiff = abs((v2-v1)/v1)   ! old version, but rdiff(v1,v2) /= vdiff(v2,v1)
00922          rdiff = abs(2.0_R8*(v2-v1)/(v1+v2))
00923       endif
00924       !------------------------------------------
00925    end function rdiff
00926 
00927 end function shr_dmodel_gGridCompare
00928 
00929 !===============================================================================
00930 
00931 subroutine shr_dmodel_mapSet_global(smatp,ggridS,gsmapS,nxgS,nygS,ggridD,gsmapD,nxgD,nygD, &
00932            name,type,algo,mask,vect,compid,mpicom,strategy)
00933 
00934     use shr_map_mod
00935     implicit none
00936 
00937     !----- arguments -----
00938     type(mct_sMatP), intent(inout) :: smatp
00939     type(mct_gGrid), intent(in)    :: ggridS
00940     type(mct_gsmap), intent(in)    :: gsmapS
00941     integer(IN)    , intent(in)    :: nxgS
00942     integer(IN)    , intent(in)    :: nygS
00943     type(mct_gGrid), intent(in)    :: ggridD
00944     type(mct_gsmap), intent(in)    :: gsmapD
00945     integer(IN)    , intent(in)    :: nxgD
00946     integer(IN)    , intent(in)    :: nygD
00947     character(len=*),intent(in)    :: name
00948     character(len=*),intent(in)    :: type
00949     character(len=*),intent(in)    :: algo
00950     character(len=*),intent(in)    :: mask
00951     character(len=*),intent(in)    :: vect
00952     integer(IN)    , intent(in)    :: compid
00953     integer(IN)    , intent(in)    :: mpicom
00954     character(len=*),intent(in),optional :: strategy
00955 
00956     !----- local -----
00957 
00958     integer(IN) :: n,i,j
00959     integer(IN) :: lsizeS,gsizeS,lsizeD,gsizeD
00960     integer(IN) :: nlon,nlat,nmsk
00961     integer(IN) :: my_task,master_task,ierr
00962 
00963     real(R8)   , pointer :: Xsrc(:,:)
00964     real(R8)   , pointer :: Ysrc(:,:)
00965     integer(IN), pointer :: Msrc(:,:)
00966     real(R8)   , pointer :: Xdst(:,:)
00967     real(R8)   , pointer :: Ydst(:,:)
00968     integer(IN), pointer :: Mdst(:,:)
00969     type(shr_map_mapType) :: shrmap
00970     type(mct_aVect) :: AVl
00971     type(mct_aVect) :: AVg
00972 
00973     character(len=32) :: lstrategy
00974     integer(IN) :: nsrc,ndst,nwts
00975     integer(IN), pointer :: isrc(:)
00976     integer(IN), pointer :: idst(:)
00977     real(R8)   , pointer :: wgts(:)
00978     type(mct_sMat) :: sMat0
00979 
00980     character(*), parameter :: subname = '(shr_dmodel_mapSet_global) '
00981     character(*), parameter :: F00   = "('(shr_dmodel_mapSet_global) ',8a)"
00982     character(*), parameter :: F01   = "('(shr_dmodel_mapSet_global) ',a,5i8)"
00983 
00984 !-------------------------------------------------------------------------------
00985 ! PURPOSE:  Initialize sMatP from mct gGrid
00986 !-------------------------------------------------------------------------------
00987 
00988     call MPI_COMM_RANK(mpicom,my_task,ierr)
00989     master_task = 0
00990 
00991     !--- get sizes and allocate for SRC ---
00992 
00993     lsizeS = mct_aVect_lsize(ggridS%data)
00994     call mct_avect_init(AVl,rList='lon:lat:mask',lsize=lsizeS)
00995     call mct_avect_copy(ggridS%data,AVl,rList='lon:lat:mask')
00996     call mct_avect_gather(AVl,AVg,gsmapS,master_task,mpicom)
00997 
00998     if (my_task == master_task) then
00999        gsizeS = mct_aVect_lsize(AVg)
01000        if (gsizeS /= nxgS*nygS) then
01001           write(logunit,F01) ' ERROR: gsizeS ',gsizeS,nxgS,nygS
01002           call shr_sys_abort(subname//' ERROR gsizeS')
01003        endif
01004        allocate(Xsrc(nxgS,nygS),Ysrc(nxgS,nygS),Msrc(nxgS,nygS))
01005 
01006        nlon = mct_avect_indexRA(AVg,'lon')
01007        nlat = mct_avect_indexRA(AVg,'lat')
01008        nmsk = mct_avect_indexRA(AVg,'mask')
01009 
01010        n = 0
01011        Msrc = 1
01012        do j = 1,nygS
01013        do i = 1,nxgS
01014           n = n + 1
01015           Xsrc(i,j) = AVg%rAttr(nlon,n)
01016           Ysrc(i,j) = AVg%rAttr(nlat,n)
01017           if (abs(AVg%rAttr(nmsk,n)) < 0.5_R8) Msrc(i,j) = 0
01018        enddo
01019        enddo
01020     endif
01021 
01022     if (my_task == master_task) call mct_aVect_clean(AVg)
01023     call mct_aVect_clean(AVl)
01024 
01025     !--- get sizes and allocate for DST ---
01026 
01027     lsizeD = mct_aVect_lsize(ggridD%data)
01028     call mct_avect_init(AVl,rList='lon:lat:mask',lsize=lsizeD)
01029     call mct_avect_copy(ggridD%data,AVl,rList='lon:lat:mask')
01030     call mct_avect_gather(AVl,AVg,gsmapD,master_task,mpicom)
01031 
01032     if (my_task == master_task) then
01033        gsizeD = mct_aVect_lsize(AVg)
01034        if (gsizeD /= nxgD*nygD) then
01035           write(logunit,F01) ' ERROR: gsizeD ',gsizeD,nxgD,nygD
01036           call shr_sys_abort(subname//' ERROR gsizeD')
01037        endif
01038        allocate(Xdst(nxgD,nygD),Ydst(nxgD,nygD),Mdst(nxgD,nygD))
01039 
01040        nlon = mct_avect_indexRA(AVg,'lon')
01041        nlat = mct_avect_indexRA(AVg,'lat')
01042        nmsk = mct_avect_indexRA(AVg,'mask')
01043 
01044        n = 0
01045        Mdst = 1
01046        do j = 1,nygD
01047        do i = 1,nxgD
01048           n = n + 1
01049           Xdst(i,j) = AVg%rAttr(nlon,n)
01050           Ydst(i,j) = AVg%rAttr(nlat,n)
01051           if (abs(AVg%rAttr(nmsk,n)) < 0.5_R8) Mdst(i,j) = 0
01052        enddo
01053        enddo
01054     endif
01055 
01056     if (my_task == master_task) call mct_aVect_clean(AVg)
01057     call mct_aVect_clean(AVl)
01058 
01059     !--- set map ---
01060 
01061     if (my_task == master_task) then
01062        call shr_map_mapSet(shrmap,Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst, &
01063           trim(name),trim(type),trim(algo),trim(mask),trim(vect))
01064        deallocate(Xsrc,Ysrc,Msrc)
01065        deallocate(Xdst,Ydst,Mdst)
01066     endif
01067 
01068     !--- convert map to sMatP ---
01069 
01070     lstrategy = 'Xonly'
01071     if (present(strategy)) then
01072        lstrategy = trim(strategy)
01073     endif
01074 
01075     if (my_task == master_task) then
01076        call shr_map_get(shrmap,shr_map_fs_nsrc,nsrc)
01077        call shr_map_get(shrmap,shr_map_fs_ndst,ndst)
01078        call shr_map_get(shrmap,shr_map_fs_nwts,nwts)
01079        allocate(isrc(nwts),idst(nwts),wgts(nwts))
01080        call shr_map_get(shrmap,isrc,idst,wgts)
01081        call shr_map_clean(shrmap)
01082 
01083        call mct_sMat_init(sMat0,ndst,nsrc,nwts)
01084 
01085        call mct_sMat_ImpGColI (sMat0,isrc,nwts)
01086        call mct_sMat_ImpGRowI (sMat0,idst,nwts)
01087        call mct_sMat_ImpMatrix(sMat0,wgts,nwts)
01088        deallocate(isrc,idst,wgts)
01089     endif
01090 
01091     call mct_sMatP_Init(sMatP,sMat0,gsmapS,gsmapD,lstrategy,master_task,mpicom,compid)
01092 
01093     if (my_task == master_task) then
01094        call mct_sMat_clean(sMat0)
01095     endif
01096 
01097 end subroutine shr_dmodel_mapSet_global
01098 
01099 !===============================================================================
01100 
01101 subroutine shr_dmodel_mapSet_dest(smatp,ggridS,gsmapS,nxgS,nygS,ggridD,gsmapD,nxgD,nygD, &
01102            name,type,algo,mask,vect,compid,mpicom,strategy)
01103 
01104     use shr_map_mod
01105     implicit none
01106 
01107     !----- arguments -----
01108     type(mct_sMatP), intent(inout) :: smatp
01109     type(mct_gGrid), intent(in)    :: ggridS
01110     type(mct_gsmap), intent(in)    :: gsmapS
01111     integer(IN)    , intent(in)    :: nxgS
01112     integer(IN)    , intent(in)    :: nygS
01113     type(mct_gGrid), intent(in)    :: ggridD
01114     type(mct_gsmap), intent(in)    :: gsmapD
01115     integer(IN)    , intent(in)    :: nxgD
01116     integer(IN)    , intent(in)    :: nygD
01117     character(len=*),intent(in)    :: name
01118     character(len=*),intent(in)    :: type
01119     character(len=*),intent(in)    :: algo
01120     character(len=*),intent(in)    :: mask
01121     character(len=*),intent(in)    :: vect
01122     integer(IN)    , intent(in)    :: compid
01123     integer(IN)    , intent(in)    :: mpicom
01124     character(len=*),intent(in),optional :: strategy
01125 
01126     !----- local -----
01127 
01128     integer(IN) :: n,i,j
01129     integer(IN) :: lsizeS,gsizeS,lsizeD,gsizeD
01130     integer(IN) :: nlon,nlat,nmsk
01131     integer(IN) :: my_task,master_task,ierr
01132 
01133     real(R8)   , pointer :: Xsrc(:,:)
01134     real(R8)   , pointer :: Ysrc(:,:)
01135     integer(IN), pointer :: Msrc(:,:)
01136     real(R8)   , pointer :: Xdst(:)
01137     real(R8)   , pointer :: Ydst(:)
01138     integer(IN), pointer :: Mdst(:)
01139     type(shr_map_mapType) :: shrmap
01140     type(mct_aVect) :: AVl
01141     type(mct_aVect) :: AVg
01142 
01143     character(len=32) :: lstrategy
01144     integer(IN) :: nsrc,ndst,nwts
01145     integer(IN), pointer :: points(:)
01146     integer(IN), pointer :: isrc(:)
01147     integer(IN), pointer :: idst(:)
01148     real(R8)   , pointer :: wgts(:)
01149     type(mct_sMat) :: sMat0
01150 
01151     character(*), parameter :: subname = '(shr_dmodel_mapSet_dest) '
01152     character(*), parameter :: F00   = "('(shr_dmodel_mapSet_dest) ',8a)"
01153     character(*), parameter :: F01   = "('(shr_dmodel_mapSet_dest) ',a,5i8)"
01154 
01155 !-------------------------------------------------------------------------------
01156 ! PURPOSE:  Initialize sMatP from mct gGrid
01157 !-------------------------------------------------------------------------------
01158 
01159     call MPI_COMM_RANK(mpicom,my_task,ierr)
01160     master_task = 0
01161 
01162     !--- get sizes and allocate for SRC ---
01163 
01164     lsizeS = mct_aVect_lsize(ggridS%data)
01165     call mct_avect_init(AVl,rList='lon:lat:mask',lsize=lsizeS)
01166     call mct_avect_copy(ggridS%data,AVl,rList='lon:lat:mask')
01167 
01168     call mct_avect_gather(AVl,AVg,gsmapS,master_task,mpicom)
01169 
01170     allocate(Xsrc(nxgS,nygS),Ysrc(nxgS,nygS),Msrc(nxgS,nygS))
01171     if (my_task == master_task) then
01172        gsizeS = mct_aVect_lsize(AVg)
01173        if (gsizeS /= nxgS*nygS) then
01174           write(logunit,F01) ' ERROR: gsizeS ',gsizeS,nxgS,nygS
01175           call shr_sys_abort(subname//' ERROR gsizeS')
01176        endif
01177 
01178        nlon = mct_avect_indexRA(AVg,'lon')
01179        nlat = mct_avect_indexRA(AVg,'lat')
01180        nmsk = mct_avect_indexRA(AVg,'mask')
01181 
01182        n = 0
01183        Msrc = 1
01184        do j = 1,nygS
01185        do i = 1,nxgS
01186           n = n + 1
01187           Xsrc(i,j) = AVg%rAttr(nlon,n)
01188           Ysrc(i,j) = AVg%rAttr(nlat,n)
01189           if (abs(AVg%rAttr(nmsk,n)) < 0.5_R8) Msrc(i,j) = 0
01190        enddo
01191        enddo
01192     endif
01193     call shr_mpi_bcast(Xsrc,mpicom)
01194     call shr_mpi_bcast(Ysrc,mpicom)
01195     call shr_mpi_bcast(Msrc,mpicom)
01196 
01197     if (my_task == master_task) call mct_aVect_clean(AVg)
01198     call mct_aVect_clean(AVl)
01199 
01200     !--- get sizes and allocate for DST ---
01201 
01202     lsizeD = mct_aVect_lsize(ggridD%data)
01203     call mct_avect_init(AVl,rList='lon:lat:mask',lsize=lsizeD)
01204     call mct_avect_copy(ggridD%data,AVl,rList='lon:lat:mask')
01205 
01206 #if (1 == 0)
01207     call mct_avect_gather(AVl,AVg,gsmapD,master_task,mpicom)
01208 
01209     if (my_task == master_task) then
01210        gsizeD = mct_aVect_lsize(AVg)
01211        if (gsizeD /= nxgD*nygD) then
01212           write(logunit,F01) ' ERROR: gsizeD ',gsizeD,nxgD,nygD
01213           call shr_sys_abort(subname//' ERROR gsizeD')
01214        endif
01215        allocate(Xdst(nxgD,nygD),Ydst(nxgD,nygD),Mdst(nxgD,nygD))
01216 
01217        nlon = mct_avect_indexRA(AVg,'lon')
01218        nlat = mct_avect_indexRA(AVg,'lat')
01219        nmsk = mct_avect_indexRA(AVg,'mask')
01220 
01221        n = 0
01222        Mdst = 1
01223        do j = 1,nygD
01224        do i = 1,nxgD
01225           n = n + 1
01226           Xdst(i,j) = AVg%rAttr(nlon,n)
01227           Ydst(i,j) = AVg%rAttr(nlat,n)
01228           if (abs(AVg%rAttr(nmsk,n)) < 0.5_R8) Mdst(i,j) = 0
01229        enddo
01230        enddo
01231     endif
01232 
01233     if (my_task == master_task) call mct_aVect_clean(AVg)
01234 #endif
01235 
01236     allocate(Xdst(lsizeD),Ydst(lsizeD),Mdst(lsizeD))
01237 
01238     nlon = mct_avect_indexRA(AVl,'lon')
01239     nlat = mct_avect_indexRA(AVl,'lat')
01240     nmsk = mct_avect_indexRA(AVl,'mask')
01241 
01242     Mdst = 1
01243     do n = 1,lsizeD
01244        Xdst(n) = AVl%rAttr(nlon,n)
01245        Ydst(n) = AVl%rAttr(nlat,n)
01246        if (abs(AVl%rAttr(nmsk,n)) < 0.5_R8) Mdst(n) = 0
01247     enddo
01248 
01249     call mct_aVect_clean(AVl)
01250 
01251     !--- set map ---
01252 
01253     nsrc = nxgS*nygS
01254     ndst = nxgD*nygD
01255     call mct_gsmap_orderedPoints(gsmapD,my_task,points)
01256     if (size(points) /= size(Xdst)) then
01257        write(logunit,F01) ' ERROR: gsizeD ',size(points),size(Xdst)
01258        call shr_sys_abort(subname//' ERROR points size')
01259     endif
01260     call shr_map_mapSet(shrmap,Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst,ndst,points, &
01261        trim(name),trim(type),trim(algo),trim(mask),trim(vect))
01262     deallocate(points)
01263     deallocate(Xsrc,Ysrc,Msrc)
01264     deallocate(Xdst,Ydst,Mdst)
01265 
01266     !--- convert map to sMatP ---
01267 
01268     lstrategy = 'Xonly'
01269     if (present(strategy)) then
01270        lstrategy = trim(strategy)
01271     endif
01272 
01273     call shr_map_get(shrmap,shr_map_fs_nwts,nwts)
01274     allocate(isrc(nwts),idst(nwts),wgts(nwts))
01275     call shr_map_get(shrmap,isrc,idst,wgts)
01276     call shr_map_clean(shrmap)
01277 
01278     call mct_sMat_init(sMat0,ndst,nsrc,nwts)
01279 
01280     call mct_sMat_ImpLColI (sMat0,isrc,nwts)
01281     call mct_sMat_ImpLRowI (sMat0,idst,nwts)
01282     call mct_sMat_ImpMatrix(sMat0,wgts,nwts)
01283     deallocate(isrc,idst,wgts)
01284 
01285     call mct_sMatP_Init(sMatP,sMat0,gsmapS,gsmapD,master_task,mpicom,compid)
01286 
01287     call mct_sMat_clean(sMat0)
01288 
01289 end subroutine shr_dmodel_mapSet_dest
01290 
01291 !===============================================================================
01292 
01293 subroutine shr_dmodel_rearrGGrid( ggridi, ggrido, gsmap, rearr, mpicom )
01294 
01295   implicit none
01296 
01297   !----- arguments -----
01298   type(mct_ggrid), intent(in)    :: ggridi
01299   type(mct_ggrid), intent(inout) :: ggrido  
01300   type(mct_gsmap), intent(in)    :: gsmap
01301   type(mct_rearr), intent(in)    :: rearr
01302   integer(IN)    , intent(in)    :: mpicom
01303 
01304   !----- local -----
01305   integer(IN)          :: lsize      ! lsize
01306   real(R8)   , pointer :: data(:)    ! temporary
01307   integer(IN), pointer :: idata(:)   ! temporary
01308   integer(IN)          :: my_task    ! local pe number
01309   integer(IN)          :: ier        ! error code
01310   character(*), parameter :: subname = '(shr_dmodel_rearrGGrid) '
01311 
01312 !-------------------------------------------------------------------------------
01313 ! PURPOSE:  Determine MCT ggrid
01314 !-------------------------------------------------------------------------------
01315 
01316   ! Initialize mct ggrid type
01317   ! lat/lon in degrees,  area in radians^2, mask is 1 (ocean), 0 (non-ocean)
01318 
01319   call mpi_comm_rank(mpicom,my_task,ier)
01320 
01321   lsize = mct_gsMap_lsize(gsMap, mpicom)
01322   call mct_gGrid_init( ggrido, ggridi, lsize=lsize )
01323 
01324   ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid
01325 
01326   call mct_gsMap_orderedPoints(gsMap, my_task, idata)
01327   call mct_gGrid_importIAttr(ggrido,'GlobGridNum',idata,lsize)
01328   deallocate(idata)
01329 
01330   ! Initialize attribute vector with special value
01331 
01332   allocate(data(lsize))
01333 
01334   data(:) = -9999.0_R8 
01335   call mct_gGrid_importRAttr(ggrido,"lat"  ,data,lsize) 
01336   call mct_gGrid_importRAttr(ggrido,"lon"  ,data,lsize) 
01337   call mct_gGrid_importRAttr(ggrido,"area" ,data,lsize) 
01338   call mct_gGrid_importRAttr(ggrido,"aream",data,lsize) 
01339   data(:) = 0.0_R8     
01340   call mct_gGrid_importRAttr(ggrido,"mask",data,lsize) 
01341   call mct_gGrid_importRAttr(ggrido,"frac",data,lsize) 
01342 
01343   deallocate(data)
01344 
01345   call mct_rearr_rearrange(ggridi%data, ggrido%data, rearr)
01346 
01347 end subroutine shr_dmodel_rearrGGrid
01348 
01349 !===============================================================================
01350 
01351 subroutine shr_dmodel_translateAV( avi, avo, avifld, avofld, rearr )
01352 
01353   implicit none
01354 
01355   !----- arguments -----
01356   type(mct_aVect), intent(in)    :: avi       ! input av
01357   type(mct_aVect), intent(inout) :: avo       ! output av
01358   character(len=*),intent(in)    :: avifld(:) ! input field names for translation
01359   character(len=*),intent(in)    :: avofld(:) ! output field names for translation
01360   type(mct_rearr), intent(in),optional :: rearr     ! rearranger for diff decomp
01361 
01362   !----- local -----
01363   integer(IN)      :: n,k,ka,kb,kc,cnt  ! indices  
01364   integer(IN)      :: lsize      ! lsize
01365   integer(IN)      :: gsize      ! gsize
01366   integer(IN)      :: nflds      ! number of fields in avi
01367 
01368   type(mct_aVect)  :: avtri,avtro ! translated av on input/output grid
01369   character(CL)    :: ilist      ! input list for translation
01370   character(CL)    :: olist      ! output list for translation
01371   character(CL)    :: cfld       ! character field name
01372   type(mct_string) :: sfld       ! string field
01373   integer(IN)      :: ktrans
01374   character(*), parameter :: subname = '(shr_dmodel_translateAV) '
01375 
01376 !-------------------------------------------------------------------------------
01377 ! PURPOSE:  Fill avo from avi
01378 !-------------------------------------------------------------------------------
01379 
01380   if (size(avifld) /= size(avofld)) then
01381      write(logunit,*) subname,' ERROR": avi and avo fld list ',size(avifld),size(avofld)
01382   endif
01383   ktrans = size(avifld)
01384 
01385   ! generate fld lists
01386   nflds = mct_aVect_nRattr(avi)
01387   cnt = 0
01388   do ka = 1,nflds
01389      call mct_aVect_getRList(sfld,ka,avi)
01390      cfld = mct_string_toChar(sfld)
01391      call mct_string_clean(sfld)
01392 
01393      k = 0
01394      kb = 0
01395      kc = 0
01396      do while (kb == 0 .and. k < ktrans)
01397         k = k + 1
01398         if (trim(avifld(k)) == trim(cfld)) then
01399            kb = k
01400            kc = mct_aVect_indexRA(avo,trim(avofld(kb)),perrWith='quiet')
01401            if (ka > 0 .and. kc > 0) then
01402               cnt = cnt + 1
01403               if (cnt == 1) then
01404                  ilist = trim(avifld(kb))
01405                  olist = trim(avofld(kb))
01406               else
01407                  ilist = trim(ilist)//':'//trim(avifld(kb))
01408                  olist = trim(olist)//':'//trim(avofld(kb))
01409               endif
01410            endif
01411         endif
01412      enddo
01413   enddo
01414 
01415   if (cnt > 0) then
01416      lsize = mct_avect_lsize(avi)
01417      call mct_avect_init(avtri,rlist=olist,lsize=lsize)
01418      call mct_avect_Copy(avi,avtri,rList=ilist,TrList=olist)
01419 
01420      if (present(rearr)) then
01421         lsize = mct_avect_lsize(avo)
01422         call mct_avect_init(avtro,rlist=olist,lsize=lsize)
01423         call mct_avect_zero(avtro)
01424         call mct_rearr_rearrange(avtri, avtro, rearr)
01425         call mct_avect_Copy(avtro,avo)
01426         call mct_aVect_clean(avtro)
01427      else
01428         call mct_avect_Copy(avtri,avo)
01429      endif
01430 
01431      call mct_aVect_clean(avtri)
01432   endif
01433 
01434 end subroutine shr_dmodel_translateAV
01435 !===============================================================================
01436 
01437 end module shr_dmodel_mod