CESM_Driver: shr_mct_mod.F90 Source File

00001 !===============================================================================
00002 ! SVN $Id: shr_mct_mod.F90 18548 2009-09-26 23:55:51Z tcraig $
00003 ! SVN $URL: https://svn-ccsm-models.cgd.ucar.edu/csm_share/trunk_tags/share3_091114/shr/shr_mct_mod.F90 $
00004 !===============================================================================
00005 !BOP ===========================================================================
00006 !
00007 ! !MODULE: shr_mct_mod -- higher level mct type routines
00008 !     needed to prevent some circular dependencies
00009 !
00010 ! !REVISION HISTORY:
00011 !     2009-Dec-16 - T. Craig - first prototype
00012 !
00013 ! !INTERFACE: ------------------------------------------------------------------
00014 module shr_mct_mod
00015 
00016 ! !USES:
00017 
00018    use shr_kind_mod         ! shared kinds
00019    use shr_sys_mod          ! share system routines
00020    use shr_mpi_mod          ! mpi layer
00021    use shr_const_mod        ! constants
00022    use mct_mod
00023 
00024    use shr_log_mod          ,only: s_loglev               => shr_log_Level
00025    use shr_log_mod          ,only: s_logunit              => shr_log_Unit
00026 
00027    implicit none
00028    private
00029 
00030 ! PUBLIC: Public interfaces
00031 
00032    public :: shr_mct_sMatReadnc
00033    public :: shr_mct_sMatPInitnc
00034    public :: shr_mct_sMatReaddnc
00035    public :: shr_mct_sMatWritednc
00036 
00037 !EOP
00038 
00039    !--- local kinds ---
00040    integer,parameter,private :: R8 = SHR_KIND_R8
00041    integer,parameter,private :: IN = SHR_KIND_IN
00042    integer,parameter,private :: CL = SHR_KIND_CL
00043 
00044 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00045 contains
00046 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00047 
00048 !===============================================================================
00049 !===============================================================================
00050 !BOP ===========================================================================
00051 !
00052 ! !IROUTINE: shr_mct_sMatReadnc - read all mapping data from a NetCDF SCRIP file
00053 !                              in to a full SparseMatrix
00054 !
00055 ! !DESCRIPTION:
00056 !   Read in mapping matrix data from a SCRIP netCDF data file so a sMat.
00057 !
00058 ! !REMARKS:
00059 !   Based on cpl_map_read
00060 !
00061 ! !REVISION HISTORY:
00062 !     2006 Nov 27: R. Jacob
00063 !
00064 ! !INTERFACE: ------------------------------------------------------------------
00065 
00066 subroutine shr_mct_sMatReadnc(sMat,fileName)
00067 
00068 #include <netcdf.inc>
00069 
00070 ! !INPUT/OUTPUT PARAMETERS:
00071 
00072    type(mct_sMat),intent(inout)  :: sMat
00073    character(*),intent(in)  :: filename  ! netCDF file to read
00074 
00075 !EOP
00076 
00077    !--- local ---
00078    integer(IN)           :: n       ! generic loop indicies
00079    integer(IN)           :: na      ! size of source domain
00080    integer(IN)           :: nb      ! size of destination domain
00081    integer(IN)           :: ns      ! number of non-zero elements in matrix
00082    integer(IN)           :: ni,nj   ! number of row and col in the matrix
00083    integer(IN)           :: igrow   ! aVect index for matrix row
00084    integer(IN)           :: igcol   ! aVect index for matrix column
00085    integer(IN)           :: iwgt    ! aVect index for matrix element
00086 
00087    real(R8)   ,allocatable :: rtemp(:)  ! reals
00088    integer(IN),allocatable :: itemp(:)  ! ints
00089 
00090    integer(IN)           :: rcode   ! netCDF routine return code
00091    integer(IN)           :: fid     ! netCDF file      ID
00092    integer(IN)           :: vid     ! netCDF variable  ID
00093    integer(IN)           :: did     ! netCDF dimension ID
00094 
00095    character(*),parameter :: subName = '(shr_mct_sMatReadnc) '
00096    character(*),parameter :: F00 = "('(shr_mct_sMatReadnc) ',4a)"
00097    character(*),parameter :: F01 = '("(shr_mct_sMatReadnc) ",2(a,i9))'
00098 
00099    if (s_loglev > 0) write(s_logunit,F00) "reading mapping matrix data..."
00100 
00101    !----------------------------------------------------------------------------
00102    ! open & read the file
00103    !----------------------------------------------------------------------------
00104    if (s_loglev > 0) write(s_logunit,F00) "* file name                  : ",trim(fileName)
00105    rcode = nf_open(filename,NF_NOWRITE,fid)
00106    if (rcode /= NF_NOERR) then
00107       write(s_logunit,F00) nf_strerror(rcode)
00108       call mct_die(subName,"error opening Netcdf file")
00109    endif
00110 
00111    !--- allocate memory & get matrix data ----------
00112    rcode = nf_inq_dimid (fid, 'n_s', did)  ! size of sparse matrix
00113    rcode = nf_inq_dimlen(fid, did  , ns)
00114    rcode = nf_inq_dimid (fid, 'n_a', did)  ! size of  input vector
00115    rcode = nf_inq_dimlen(fid, did  , na)
00116    rcode = nf_inq_dimid (fid, 'n_b', did)  ! size of output vector
00117    rcode = nf_inq_dimlen(fid, did  , nb)
00118 
00119    if (s_loglev > 0) write(s_logunit,F01) "* matrix dimensions src x dst: ",na,' x',nb
00120    if (s_loglev > 0) write(s_logunit,F01) "* number of non-zero elements: ",ns
00121 
00122    !----------------------------------------------------------------------------
00123    ! init the mct sMat data type
00124    !----------------------------------------------------------------------------
00125    ! mct_sMat_init must be given the number of rows and columns that
00126    ! would be in the full matrix.  Nrows= size of output vector=nb.
00127    ! Ncols = size of input vector = na.
00128    call mct_sMat_init(sMat, nb, na, ns)
00129 
00130    igrow = mct_sMat_indexIA(sMat,'grow')
00131    igcol = mct_sMat_indexIA(sMat,'gcol')
00132    iwgt  = mct_sMat_indexRA(sMat,'weight')
00133 
00134    !!!!!!!!!!!!!!!!!!!!!!!!!!
00135    ! read and load matrix weights
00136    allocate(rtemp(ns),stat=rcode)
00137    if (rcode /= 0) &
00138      call mct_die(subName,':: allocate weights',rcode)
00139 
00140    rcode = nf_inq_varid     (fid,'S'  ,vid)
00141    rcode = nf_get_var_double(fid,vid  ,rtemp  )
00142    if (rcode /= NF_NOERR .and. s_loglev > 0) write(s_logunit,F00) nf_strerror(rcode)
00143 
00144    sMat%data%rAttr(iwgt ,:) =   rtemp(:)
00145 
00146    deallocate(rtemp, stat=rcode)
00147    if (rcode /= 0) call mct_perr_die(subName,':: deallocate weights',rcode)
00148 
00149    !!!!!!!!!!!!!!!!!!!!!!!!!!
00150    ! read and load rows
00151    allocate(itemp(ns),stat=rcode)
00152    if (rcode /= 0) call mct_perr_die(subName,':: allocate rows',rcode)
00153 
00154    rcode = nf_inq_varid     (fid,'row',vid)
00155    rcode = nf_get_var_int   (fid,vid  ,itemp)
00156    if (rcode /= NF_NOERR .and. s_loglev > 0) write(s_logunit,F00) nf_strerror(rcode)
00157 
00158    sMat%data%iAttr(igrow,:) = itemp(:)
00159 
00160 
00161    !!!!!!!!!!!!!!!!!!!!!!!!!!
00162    ! read and load columns
00163    itemp(:) = 0
00164 
00165    rcode = nf_inq_varid     (fid,'col',vid)
00166    rcode = nf_get_var_int   (fid,vid  ,itemp)
00167    if (rcode /= NF_NOERR .and. s_loglev > 0) write(s_logunit,F00) nf_strerror(rcode)
00168 
00169    sMat%data%iAttr(igcol,:) = itemp(:)
00170 
00171    deallocate(itemp, stat=rcode)
00172    if (rcode /= 0) call mct_perr_die(subName,':: deallocate cols',rcode)
00173 
00174    rcode = nf_close(fid)
00175 
00176    if (s_loglev > 0) write(s_logunit,F00) "... done reading file"
00177    call shr_sys_flush(s_logunit)
00178 
00179 end subroutine shr_mct_sMatReadnc
00180 !===============================================================================
00181 !BOP ===========================================================================
00182 !
00183 ! !IROUTINE: shr_mct_sMatPInitnc - initialize a SparseMatrixPlus.
00184 !
00185 ! !DESCRIPTION:
00186 !   Read in mapping matrix data from a SCRIP netCDF data file in first an
00187 !   Smat and then an SMatPlus
00188 !
00189 ! !REMARKS:
00190 !
00191 ! !REVISION HISTORY:
00192 !     2006 Nov 27: R. Jacob
00193 !
00194 ! !INTERFACE: ------------------------------------------------------------------
00195 
00196 subroutine shr_mct_sMatPInitnc(sMatP, gsMapX, gsMapY, ConfigFileName, &
00197                             MapLabel, MapTypeLabel, mpicom,&
00198                             ni_i, nj_i, ni_o, nj_o, &
00199                             areasrc, areadst)
00200 
00201 ! !INPUT/OUTPUT PARAMETERS:
00202 
00203    type(mct_sMatP),intent(inout)         :: sMatP
00204    type(mct_gsMap),intent(in)            :: gsMapX
00205    type(mct_gsMap),intent(in)            :: gsMapY
00206    character(*)   ,intent(in)            :: ConfigFileName  ! config file to read
00207    character(*)   ,intent(in)            :: MapLabel        ! map name
00208    character(*)   ,intent(in)            :: MapTypeLabel    ! map type
00209    integer        ,intent(in)            :: mpicom
00210    integer        ,intent(out), optional :: ni_i            ! number of longitudes on input grid
00211    integer        ,intent(out), optional :: nj_i            ! number of latitudes  on input grid
00212    integer        ,intent(out), optional :: ni_o            ! number of longitudes on output grid
00213    integer        ,intent(out), optional :: nj_o            ! number of latitudes  on output grid
00214    type(mct_Avect),intent(out), optional :: areasrc         ! area of src grid from mapping file
00215    type(mct_Avect),intent(out), optional :: areadst         ! area of src grid from mapping file
00216 
00217 !EOP
00218    type(mct_sMat ) :: sMati    ! initial sMat from read (either root or decomp)
00219    type(mct_Avect) :: areasrc_map ! area of src grid from mapping file
00220    type(mct_Avect) :: areadst_map ! area of dst grid from mapping file
00221 
00222    integer          :: lsize
00223    integer          :: iret
00224    integer          :: pe_loc
00225    character(256)   :: fileName
00226    character(1)     :: maptype
00227    logical          :: usevector 
00228    character(len=3) :: Smaptype
00229    character(*),parameter :: areaAV_field = 'aream'
00230    character(*),parameter :: F00 = "('(shr_mct_sMatPInitnc) ',4a)"
00231    character(*),parameter :: F01 = "('(shr_mct_sMatPInitnc) ',a,i10)"
00232 
00233    call shr_mpi_commrank(mpicom,pe_loc)
00234 
00235    if (s_loglev > 0) write(s_logunit,*) " "
00236    if (s_loglev > 0) write(s_logunit,F00) "Initializing SparseMatrixPlus"
00237    call I90_allLoadF(ConfigFileName,0,mpicom,iret)
00238       if(iret /= 0) then
00239         write(s_logunit,*)"Cant find config file ",ConfigFileName
00240         call mct_die("mct_sMapP_initnc","File Not Found")
00241       endif
00242 
00243    call i90_label(trim(MapLabel),iret)
00244       if(iret /= 0) then
00245         write(s_logunit,*)"Cant find label ",MapLabel
00246         call mct_die("mct_sMapP_initnc","Label Not Found")
00247       endif
00248 
00249    call i90_gtoken(fileName,iret)
00250       if(iret /= 0) then
00251         write(s_logunit,*)"Error reading token ",fileName
00252         call mct_die("mct_sMapP_initnc","Error on read")
00253       endif
00254    if (s_loglev > 0) write(s_logunit,F00) "map weight filename ",trim(fileName)
00255 
00256    call i90_label(trim(MapTypeLabel),iret)
00257    call i90_gtoken(maptype,iret)
00258 
00259    if (s_loglev > 0) write(s_logunit,F00) "SmatP maptype ",maptype
00260 
00261    if (maptype == "X") then
00262       Smaptype = "src"
00263    else if(maptype == "Y") then
00264       Smaptype = "dst"
00265    end if
00266 
00267    call shr_mpi_commrank(mpicom, pe_loc)
00268 
00269    lsize = mct_gsMap_lsize(gsMapX, mpicom)
00270    call mct_aVect_init(areasrc_map, rList=areaAV_field, lsize=lsize)
00271 
00272    lsize = mct_gsMap_lsize(gsMapY, mpicom)
00273    call mct_aVect_init(areadst_map, rList=areaAV_field, lsize=lsize)
00274 
00275    if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
00276       call shr_mct_sMatReaddnc(sMati, gsMapX, gsMapY, Smaptype, areasrc_map, areadst_map, &
00277            fileName, pe_loc, mpicom, ni_i, nj_i, ni_o, nj_o)
00278    else
00279       call shr_mct_sMatReaddnc(sMati, gsMapX, gsMapY, Smaptype, areasrc_map, areadst_map, &
00280            fileName, pe_loc, mpicom)
00281    end if
00282    call mct_sMatP_Init(sMatP, sMati, gsMapX, gsMapY, 0, mpicom, gsMapX%comp_id)
00283 
00284 #ifdef CPP_VECTOR
00285    !--- initialize the vector parts of the sMat ---
00286    call mct_sMatP_Vecinit(sMatP)
00287 #endif
00288 
00289    lsize = mct_smat_gNumEl(sMatP%Matrix,mpicom)
00290    if (s_loglev > 0) write(s_logunit,F01) "Done initializing SmatP, nElements = ",lsize
00291 
00292 #ifdef CPP_VECTOR
00293    usevector = .true.
00294 #else
00295    usevector = .false.
00296 #endif
00297    if (present(areasrc)) then
00298       call mct_aVect_copy(aVin=areasrc_map, aVout=areasrc, vector=usevector)
00299    end if
00300    if (present(areadst)) then
00301       call mct_aVect_copy(aVin=areadst_map, aVout=areadst, vector=usevector)
00302    end if
00303 
00304    call mct_aVect_clean(areasrc_map)
00305    call mct_aVect_clean(areadst_map)
00306 
00307    call mct_sMat_Clean(sMati)
00308    call I90_Release(iret)
00309 
00310 end subroutine shr_mct_sMatPInitnc
00311 
00312 !BOP ===========================================================================
00313 !
00314 ! !IROUTINE:  shr_mct_sMatReaddnc - Do a distributed read of a NetCDF SCRIP file and
00315 !                                return weights in a distributed SparseMatrix
00316 !
00317 ! !DESCRIPTION: 
00318 !     Read in mapping matrix data from a SCRIP netCDF data file using
00319 !     a low memory method and then scatter to all pes.
00320 !
00321 ! !REMARKS:
00322 !   This routine leverages gsmaps to determine scatter pattern
00323 !   The scatter is implemented as a bcast of all weights then a local
00324 !     computation on each pe to determine with weights to keep based
00325 !     on gsmap information.
00326 !   The algorithm to determine whether a weight belongs on a pe involves
00327 !     creating a couple local arrays (lsstart and lscount) which are
00328 !     the local values of start and length from the gsmap.  these are
00329 !     sorted via a bubble sort and then searched via a binary search
00330 !     to check whether a global index is on the local pe.
00331 !   The local buffer sizes are estimated up front based on ngridcell/npes
00332 !     plus 20% (see 1.2 below).  If the local buffer size fills up, then
00333 !     the buffer is reallocated 50% large (see 1.5 below) and the fill
00334 !     continues.  The idea is to trade off memory reallocation and copy
00335 !     with memory usage.  1.2 and 1.5 are arbitary, other values may
00336 !     result in better performance.
00337 !   Once all the matrix weights have been read, the sMat is initialized,
00338 !     the values from the buffers are copied in, and everything is deallocated.
00339 
00340 ! !SEE ALSO: 
00341 !    mct/m_SparseMatrix.F90 (MCT source code)
00342 !
00343 ! !REVISION HISTORY: 
00344 !    2007-Jan-18 - T. Craig -- first version
00345 !    2007-Mar-20 - R. Jacob -- rename to shr_mct_sMatReaddnc.  Remove use of cpl_
00346 !                  variables and move to shr_mct_mod
00347 ! 
00348 ! !INTERFACE:  -----------------------------------------------------------------
00349 
00350 subroutine shr_mct_sMatReaddnc(sMat,SgsMap,DgsMap,newdom,areasrc,areadst, &
00351                             fileName,mytask, mpicom, ni_i,nj_i,ni_o,nj_o )
00352 
00353 ! !USES:
00354 
00355 #include <netcdf.inc>
00356 
00357 ! !INPUT/OUTPUT PARAMETERS:
00358 
00359    type(mct_sMat)  ,intent(out)           :: sMat    ! mapping data
00360    type(mct_gsMap) ,intent(in) ,target    :: SgsMap  ! src gsmap
00361    type(mct_gSMap) ,intent(in) ,target    :: DgsMap  ! dst gsmap
00362    character(*)    ,intent(in)            :: newdom  ! type of sMat (src or dst)
00363    type(mct_Avect) ,intent(out), optional :: areasrc ! area of src grid from mapping file
00364    type(mct_Avect) ,intent(out), optional :: areadst ! area of dst grid from mapping file
00365    character(*)    ,intent(in)            :: filename! netCDF file to read
00366    integer(IN)     ,intent(in)            :: mytask   ! processor id
00367    integer(IN)     ,intent(in)            :: mpicom  ! communicator
00368    integer(IN)     ,intent(out), optional :: ni_i    ! number of lons on input grid   
00369    integer(IN)     ,intent(out), optional :: nj_i    ! number of lats on input grid   
00370    integer(IN)     ,intent(out), optional :: ni_o    ! number of lons on output grid   
00371    integer(IN)     ,intent(out), optional :: nj_o    ! number of lats on output grid   
00372 
00373 ! !EOP
00374 
00375    !--- local ---
00376    integer(IN)           :: n,m     ! generic loop indicies
00377    integer(IN)           :: na      ! size of source domain
00378    integer(IN)           :: nb      ! size of destination domain
00379    integer(IN)           :: ns      ! number of non-zero elements in matrix
00380    integer(IN)           :: ni,nj   ! number of row and col in the matrix
00381    integer(IN)           :: igrow   ! aVect index for matrix row
00382    integer(IN)           :: igcol   ! aVect index for matrix column
00383    integer(IN)           :: iwgt    ! aVect index for matrix element
00384    integer(IN)           :: iarea   ! aVect index for area
00385    integer(IN)           :: rsize   ! size of read buffer
00386    integer(IN)           :: cnt     ! local num of wgts
00387    integer(IN)           :: cntold  ! local num of wgts, previous read
00388    integer(IN)           :: start(1)! netcdf read
00389    integer(IN)           :: count(1)! netcdf read
00390    integer(IN)           :: bsize   ! buffer size
00391    integer(IN)           :: nread   ! number of reads 
00392    logical               :: mywt    ! does this weight belong on my pe
00393 
00394    !--- buffers for i/o ---
00395    real(R8)   ,allocatable :: rtemp(:) ! real temporary
00396    real(R8)   ,allocatable :: Sbuf(:)  ! real weights
00397    integer(IN),allocatable :: Rbuf(:)  ! ints rows
00398    integer(IN),allocatable :: Cbuf(:)  ! ints cols
00399 
00400    !--- variables associated with local computation of global indices
00401    integer(IN)             :: lsize     ! size of local seg map
00402    integer(IN)             :: commsize  ! size of local communicator
00403    integer(IN),allocatable :: lsstart(:) ! local seg map info
00404    integer(IN),allocatable :: lscount(:) ! local seg map info
00405    type(mct_gsMap),pointer :: mygsmap ! pointer to one of the gsmaps
00406    integer(IN)             :: l1,l2     ! generice indices for sort
00407    logical                 :: found     ! for sort
00408 
00409    !--- variable assocaited with local data buffers and reallocation
00410    real(R8)   ,allocatable :: Snew(:),Sold(:)  ! reals
00411    integer(IN),allocatable :: Rnew(:),Rold(:)  ! ints
00412    integer(IN),allocatable :: Cnew(:),Cold(:)  ! ints
00413 
00414    character,allocatable :: str(:)  ! variable length char string
00415    character(CL)         :: attstr  ! netCDF attribute name string
00416    integer(IN)           :: rcode   ! netCDF routine return code
00417    integer(IN)           :: fid     ! netCDF file      ID
00418    integer(IN)           :: vid     ! netCDF variable  ID
00419    integer(IN)           :: did     ! netCDF dimension ID
00420    !--- arbitrary size of read buffer, this is the chunk size weights reading
00421    integer(IN),parameter :: rbuf_size = 100000
00422 
00423    !--- global source and destination areas ---
00424    type(mct_Avect) :: areasrc0   ! area of src grid from mapping file
00425    type(mct_Avect) :: areadst0   ! area of src grid from mapping file
00426 
00427    character(*),parameter :: areaAV_field = 'aream'
00428 
00429    !--- formats ---
00430    character(*),parameter :: subName = '(shr_mct_sMatReaddnc) '
00431    character(*),parameter :: F00 = '("(shr_mct_sMatReaddnc) ",4a)'
00432    character(*),parameter :: F01 = '("(shr_mct_sMatReaddnc) ",2(a,i10))'
00433 
00434 !-------------------------------------------------------------------------------
00435 !
00436 !-------------------------------------------------------------------------------
00437 
00438  call shr_mpi_commsize(mpicom,commsize)
00439  if (mytask == 0) then
00440    if (s_loglev > 0) write(s_logunit,F00) "reading mapping matrix data decomposed..."
00441 
00442    !----------------------------------------------------------------------------
00443    ! open & read the file
00444    !----------------------------------------------------------------------------
00445    if (s_loglev > 0) write(s_logunit,F00) "* file name                  : ",trim(fileName)
00446    rcode = nf_open(filename,NF_NOWRITE,fid)
00447    if (rcode /= NF_NOERR .and. s_loglev > 0) write(s_logunit,F00) nf_strerror(rcode)
00448 
00449    !--- get matrix dimensions ----------
00450    rcode = nf_inq_dimid (fid, 'n_s', did)  ! size of sparse matrix
00451    rcode = nf_inq_dimlen(fid, did  , ns)
00452    rcode = nf_inq_dimid (fid, 'n_a', did)  ! size of  input vector
00453    rcode = nf_inq_dimlen(fid, did  , na)
00454    rcode = nf_inq_dimid (fid, 'n_b', did)  ! size of output vector
00455    rcode = nf_inq_dimlen(fid, did  , nb)
00456    
00457    if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
00458       rcode = nf_inq_dimid (fid, 'ni_a', did)  ! number of lons in input grid
00459       rcode = nf_inq_dimlen(fid, did  , ni_i)
00460       rcode = nf_inq_dimid (fid, 'nj_a', did)  ! number of lats in input grid
00461       rcode = nf_inq_dimlen(fid, did  , nj_i)
00462       rcode = nf_inq_dimid (fid, 'ni_b', did)  ! number of lons in output grid
00463       rcode = nf_inq_dimlen(fid, did  , ni_o)
00464       rcode = nf_inq_dimid (fid, 'nj_b', did)  ! number of lats in output grid
00465       rcode = nf_inq_dimlen(fid, did  , nj_o)
00466    end if
00467 
00468    if (s_loglev > 0) write(s_logunit,F01) "* matrix dims src x dst      : ",na,' x',nb
00469    if (s_loglev > 0) write(s_logunit,F01) "* number of non-zero elements: ",ns
00470 
00471  endif
00472  
00473    !--- read and load area_a ---
00474    if (present(areasrc)) then
00475    if (mytask == 0) then
00476       call mct_aVect_init(areasrc0,' ',areaAV_field,na)
00477       rcode = nf_inq_varid     (fid,'area_a',vid)
00478       if (rcode /= NF_NOERR) write(6,F00) nf_strerror(rcode)
00479       rcode = nf_get_var_double(fid, vid, areasrc0%rAttr)
00480       if (rcode /= NF_NOERR) write(6,F00) nf_strerror(rcode)
00481    endif
00482    call mct_aVect_scatter(areasrc0, areasrc, SgsMap, 0, mpicom, rcode)
00483    if (rcode /= 0) call mct_die("shr_mct_sMatReaddnc","Error on scatter of areasrc0")
00484    if (mytask == 0) then
00485 !      if (present(dbug)) then
00486 !         if (dbug > 2) then
00487 !            write(6,*) subName,'Size of src ',mct_aVect_lSize(areasrc0)
00488 !            write(6,*) subName,'min/max src ',minval(areasrc0%rAttr(1,:)),maxval(areasrc0%rAttr(1,:))
00489 !         endif
00490 !      end if
00491       call mct_aVect_clean(areasrc0)
00492    end if
00493    end if
00494 
00495    !--- read and load area_b ---
00496    if (present(areadst)) then
00497    if (mytask == 0) then
00498       call mct_aVect_init(areadst0,' ',areaAV_field,nb)
00499       rcode = nf_inq_varid     (fid,'area_b',vid)
00500       if (rcode /= NF_NOERR) write(6,F00) nf_strerror(rcode)
00501       rcode = nf_get_var_double(fid, vid, areadst0%rAttr)
00502       if (rcode /= NF_NOERR) write(6,F00) nf_strerror(rcode)
00503    endif
00504    call mct_aVect_scatter(areadst0, areadst, DgsMap, 0, mpicom, rcode)
00505    if (rcode /= 0) call mct_die("shr_mct_sMatReaddnc","Error on scatter of areadst0")
00506    if (mytask == 0) then
00507 !      if (present(dbug)) then
00508 !         if (dbug > 2) then
00509 !            write(6,*) subName,'Size of dst ',mct_aVect_lSize(areadst0)
00510 !            write(6,*) subName,'min/max dst ',minval(areadst0%rAttr(1,:)),maxval(areadst0%rAttr(1,:))
00511 !         endif
00512 !      end if
00513       call mct_aVect_clean(areadst0)
00514    endif
00515    endif
00516 
00517    if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
00518       call shr_mpi_bcast(ni_i,mpicom,subName//" MPI in ni_i bcast")
00519       call shr_mpi_bcast(nj_i,mpicom,subName//" MPI in nj_i bcast")
00520       call shr_mpi_bcast(ni_o,mpicom,subName//" MPI in ni_o bcast")
00521       call shr_mpi_bcast(nj_o,mpicom,subName//" MPI in nj_o bcast")
00522    end if
00523 
00524    call shr_mpi_bcast(ns,mpicom,subName//" MPI in ns bcast")
00525    call shr_mpi_bcast(na,mpicom,subName//" MPI in na bcast")
00526    call shr_mpi_bcast(nb,mpicom,subName//" MPI in nb bcast")
00527 
00528    !--- setup local seg map, sorted
00529    if (newdom == 'src') then
00530       mygsmap => DgsMap
00531    elseif (newdom == 'dst') then
00532       mygsmap => SgsMap
00533    else
00534       write(s_logunit,F00) 'ERROR: invalid newdom value = ',newdom
00535       call shr_sys_abort(trim(subName)//" invalid newdom value")
00536    endif
00537    lsize = 0
00538    do n = 1,size(mygsmap%start)
00539       if (mygsmap%pe_loc(n) == mytask) then
00540          lsize=lsize+1
00541       endif
00542    enddo
00543    allocate(lsstart(lsize),lscount(lsize),stat=rcode)
00544    if (rcode /= 0) call mct_perr_die(subName,':: allocate Lsstart',rcode)
00545 
00546    lsize = 0
00547    do n = 1,size(mygsmap%start)
00548       if (mygsmap%pe_loc(n) == mytask) then  ! on my pe
00549          lsize=lsize+1
00550          found = .false.
00551          l1 = 1
00552          do while (.not.found .and. l1 < lsize)         ! bubble sort copy
00553             if (mygsmap%start(n) < lsstart(l1)) then
00554                do l2 = lsize, l1+1, -1
00555                   lsstart(l2) = lsstart(l2-1)
00556                   lscount(l2) = lscount(l2-1)
00557                enddo
00558                found = .true.
00559             else
00560                l1 = l1 + 1
00561             endif
00562          enddo
00563          lsstart(l1) = mygsmap%start(n)
00564          lscount(l1) = mygsmap%length(n)
00565       endif
00566    enddo
00567    do n = 1,lsize-1
00568       if (lsstart(n) > lsstart(n+1)) then
00569          write(s_logunit,F00) ' ERROR: lsstart not properly sorted'
00570          call shr_sys_abort()
00571       endif
00572    enddo
00573 
00574    rsize = min(rbuf_size,ns)                     ! size of i/o chunks
00575    bsize = ((ns/commsize) + 1 ) * 1.2   ! local temporary buffer size
00576    if (ns == 0) then
00577       nread = 0
00578    else
00579       nread = (ns-1)/rsize + 1                      ! num of reads to do
00580    endif
00581 
00582    allocate(Sbuf(rsize),Rbuf(rsize),Cbuf(rsize),stat=rcode)
00583    if (rcode /= 0) call mct_perr_die(subName,':: allocate Sbuf',rcode)
00584    allocate(Snew(bsize),Cnew(bsize),Rnew(bsize),stat=rcode)
00585    if (rcode /= 0) call mct_perr_die(subName,':: allocate Snew1',rcode)
00586 
00587    cnt = 0
00588    do n = 1,nread
00589       start(1) = (n-1)*rsize + 1
00590       count(1) = min(rsize,ns-start(1)+1)
00591 
00592       !--- read data on root pe
00593       if (mytask== 0) then
00594          rcode = nf_inq_varid      (fid,'S'  ,vid)
00595          rcode = nf_get_vara_double(fid,vid,start,count,Sbuf)
00596          if (rcode /= NF_NOERR .and. s_loglev > 0) write(s_logunit,F00) nf_strerror(rcode)
00597 
00598          rcode = nf_inq_varid      (fid,'row',vid)
00599          rcode = nf_get_vara_int   (fid,vid,start,count,Rbuf)
00600          if (rcode /= NF_NOERR .and. s_loglev > 0) write(s_logunit,F00) nf_strerror(rcode)
00601 
00602          rcode = nf_inq_varid      (fid,'col',vid)
00603          rcode = nf_get_vara_int   (fid,vid,start,count,Cbuf)
00604          if (rcode /= NF_NOERR .and. s_loglev > 0) write(s_logunit,F00) nf_strerror(rcode)
00605       endif
00606 
00607       !--- send S, row, col to all pes
00608       call shr_mpi_bcast(Sbuf,mpicom,subName//" MPI in Sbuf bcast")
00609       call shr_mpi_bcast(Rbuf,mpicom,subName//" MPI in Rbuf bcast")
00610       call shr_mpi_bcast(Cbuf,mpicom,subName//" MPI in Cbuf bcast")
00611 
00612       !--- now each pe keeps what it should
00613       do m = 1,count(1)
00614          !--- should this weight be on my pe
00615          if (newdom == 'src') then
00616             mywt = mct_myindex(Rbuf(m),lsstart,lscount)
00617          elseif (newdom == 'dst') then
00618             mywt = mct_myindex(Cbuf(m),lsstart,lscount)
00619          endif
00620 
00621          if (mywt) then
00622             cntold = cnt
00623             cnt = cnt + 1
00624 
00625             !--- new arrays need to be bigger
00626             if (cnt > bsize) then
00627                !--- allocate old arrays and copy new into old
00628                allocate(Sold(cntold),Rold(cntold),Cold(cntold),stat=rcode)
00629                if (rcode /= 0) call mct_perr_die(subName,':: allocate old',rcode)
00630                Sold(1:cntold) = Snew(1:cntold)
00631                Rold(1:cntold) = Rnew(1:cntold)
00632                Cold(1:cntold) = Cnew(1:cntold)
00633 
00634                !--- reallocate new to bigger size, increase buffer by 50% (arbitrary)
00635                deallocate(Snew,Rnew,Cnew,stat=rcode)
00636                if (rcode /= 0) call mct_perr_die(subName,':: allocate new',rcode)
00637                bsize = 1.5 * bsize
00638                if (s_loglev > 1) write(s_logunit,F01) ' reallocate bsize to ',bsize
00639                allocate(Snew(bsize),Rnew(bsize),Cnew(bsize),stat=rcode)
00640                if (rcode /= 0) call mct_perr_die(subName,':: allocate old',rcode)
00641 
00642                !--- copy data back into new
00643                Snew(1:cntold) = Sold(1:cntold)
00644                Rnew(1:cntold) = Rold(1:cntold)
00645                Cnew(1:cntold) = Cold(1:cntold)
00646                deallocate(Sold,Rold,Cold,stat=rcode)
00647                if (rcode /= 0) call mct_perr_die(subName,':: deallocate old',rcode)
00648             endif
00649 
00650             Snew(cnt) = Sbuf(m)
00651             Rnew(cnt) = Rbuf(m)
00652             Cnew(cnt) = Cbuf(m)
00653          endif
00654       enddo  ! count
00655    enddo   ! nread
00656 
00657    deallocate(Sbuf,Rbuf,Cbuf, stat=rcode)
00658    if (rcode /= 0) call mct_perr_die(subName,':: deallocate Sbuf',rcode)
00659 
00660    !----------------------------------------------------------------------------
00661    ! init the mct sMat data type
00662    !----------------------------------------------------------------------------
00663    ! mct_sMat_init must be given the number of rows and columns that
00664    ! would be in the full matrix.  Nrows= size of output vector=nb.
00665    ! Ncols = size of input vector = na.
00666    call mct_sMat_init(sMat, nb, na, cnt)
00667 
00668    igrow = mct_sMat_indexIA(sMat,'grow')
00669    igcol = mct_sMat_indexIA(sMat,'gcol')
00670    iwgt  = mct_sMat_indexRA(sMat,'weight')
00671 
00672    if (cnt /= 0) then
00673       sMat%data%rAttr(iwgt ,1:cnt) = Snew(1:cnt)
00674       sMat%data%iAttr(igrow,1:cnt) = Rnew(1:cnt)
00675       sMat%data%iAttr(igcol,1:cnt) = Cnew(1:cnt)
00676    endif
00677    deallocate(Snew,Rnew,Cnew, stat=rcode)
00678    deallocate(lsstart,lscount,stat=rcode)
00679    if (rcode /= 0) call mct_perr_die(subName,':: deallocate new',rcode)
00680 
00681    if (mytask == 0) then
00682       rcode = nf_close(fid)
00683       if (s_loglev > 0) write(s_logunit,F00) "... done reading file"
00684       call shr_sys_flush(s_logunit)
00685    endif
00686 
00687 end subroutine shr_mct_sMatReaddnc
00688 
00689 !BOP ===========================================================================
00690 !
00691 ! !IROUTINE:  shr_mct_sMatWritednc - Do a distributed write of a NetCDF SCRIP file
00692 !                                 based on a distributed SparseMatrix
00693 !
00694 ! !DESCRIPTION: 
00695 !     Write out mapping matrix data from a SCRIP netCDF data file using
00696 !     a low memory method.
00697 !
00698 ! !SEE ALSO: 
00699 !    mct/m_SparseMatrix.F90 (MCT source code)
00700 !
00701 ! !REVISION HISTORY: 
00702 !    2009-Dec-15 - T. Craig -- first version
00703 ! 
00704 ! !INTERFACE:  -----------------------------------------------------------------
00705 
00706 subroutine shr_mct_sMatWritednc(sMat,fileName,compid, mpicom)
00707 
00708 ! !USES:
00709 
00710    use shr_pcdf_mod, only : shr_pcdf_readwrite
00711    implicit none
00712 #include <netcdf.inc>
00713 #include <mpif.h>
00714 
00715 ! !INPUT/OUTPUT PARAMETERS:
00716 
00717    type(mct_sMat)  ,intent(in)   :: sMat    ! mapping data
00718    character(*)    ,intent(in)   :: filename! netCDF file to read
00719    integer(IN)     ,intent(in)   :: compid   ! processor id
00720    integer(IN)     ,intent(in)   :: mpicom  ! communicator
00721 
00722  ! !local
00723    integer(IN) :: na,nb,ns,lsize,npes,ierr,my_task,n
00724    integer(IN), pointer :: start(:),count(:),ssize(:),pe_loc(:)
00725    integer(IN), pointer :: expvari(:)
00726    real(R8)   , pointer :: expvarr(:)
00727    type(mct_gsmap) :: gsmap
00728    type(mct_avect) :: AV
00729    character(*),parameter :: subName = '(shr_mct_sMatWritednc) '
00730 
00731 !----------------------------------------
00732 
00733    call MPI_COMM_SIZE(mpicom,npes,ierr)
00734    call MPI_COMM_RANK(mpicom,my_task,ierr)
00735    allocate(start(npes),count(npes),ssize(npes),pe_loc(npes))
00736 
00737    na = mct_sMat_ncols(sMat)
00738    nb = mct_sMat_nrows(sMat)
00739    ns = mct_sMat_gNumEl(sMat,mpicom)
00740    lsize = mct_sMat_lsize(sMat)
00741 
00742    count(:) = -999
00743    pe_loc(:) = -999
00744    ssize(:) = 1
00745    call MPI_GATHER(lsize,1,MPI_INTEGER,count,ssize,MPI_INTEGER,0,mpicom,ierr)
00746 
00747    if (my_task == 0) then
00748       if (minval(count) < 0) then
00749          call shr_sys_abort(subname//' ERROR: count invalid')
00750       endif
00751 
00752       start(1) = 1
00753       pe_loc(1) = 0
00754       do n = 2,npes
00755          start(n) = start(n-1)+count(n-1)
00756          pe_loc(n) = n-1
00757       enddo
00758 
00759    endif
00760 
00761    call mct_gsmap_init(gsmap,npes,start,count,pe_loc,0,mpicom,compid,ns)
00762    deallocate(start,count,ssize,pe_loc)
00763 
00764    call mct_aVect_init(AV,iList='row:col',rList='S',lsize=lsize)
00765    allocate(expvari(lsize))
00766    call mct_sMat_ExpGRowI(sMat,expvari)
00767    AV%iAttr(1,:) = expvari(:)
00768    call mct_sMat_ExpGColI(sMat,expvari)
00769    AV%iAttr(2,:) = expvari(:)
00770    deallocate(expvari)
00771    allocate(expvarr(lsize))
00772    call mct_sMat_ExpMatrix(sMat,expvarr)
00773    AV%rAttr(1,:) = expvarr(:)
00774    deallocate(expvarr)
00775 
00776    call shr_pcdf_readwrite('write',trim(filename),mpicom,gsmap,clobber=.false.,cdf64=.true., &
00777       id1=na,id1n='n_a',id2=nb,id2n='n_b',id3=ns,id3n='n_s',av1=AV,av1n='')
00778 
00779    call mct_gsmap_clean(gsmap)
00780    call mct_avect_clean(AV)
00781 
00782 end subroutine shr_mct_sMatWritednc
00783 !===============================================================================
00784 
00785 end module shr_mct_mod
00786