Oasis3-MCT
 All Classes Files Functions Variables Macros Pages
mod_oasis_map.F90
Go to the documentation of this file.
1 
2 !> OASIS map (interpolation) data and methods
3 
4 
6 
11  USE mod_oasis_sys
12  USE mod_oasis_part
13  USE mod_oasis_var
14  USE mod_oasis_mpi
16  USE mod_oasis_io
17  USE mod_oasis_timer
18  USE mct_mod
19  USE grids ! scrip
20  USE netcdf
21 
22  IMPLICIT NONE
23 
24  private
25 
26  public oasis_map_genmap
29 
30  !--- Type of data
31 
32  public prism_mapper_type
33 
34  !> Mapper data for interpolating data between grids
36  !--- fixed at initialization ---
37  type(mct_smatp),pointer :: smatp(:) !< stores mapping data such as weights
38  integer(kind=ip_i4_p) :: nwgts !< number of weights in weights file
39  character(len=ic_long):: file !< file to read/write
40  character(len=ic_med) :: loc !< location setting, src or dst model
41  character(len=ic_med) :: opt !< optimization setting, bfb, sum, or opt
42  character(len=ic_med) :: optval !< mct map setting, src or dst, derived from opt
43  logical :: init !< flag indicating initialization complete
44  integer(kind=ip_i4_p) :: spart !< src partition
45  integer(kind=ip_i4_p) :: dpart !< dst partition
46  logical :: avred !< flag indicating AV_ms, AV_md data has been read
47  type(mct_avect) :: av_ms !< stores data for CONSERV for src such as mask and area
48  type(mct_avect) :: av_md !< stores data for CONSERV for dst such as mask and area
49  end type prism_mapper_type
50 
51  integer(kind=ip_i4_p) ,public :: prism_mmapper !< max mappers
52  integer(kind=ip_i4_p) ,public :: prism_nmapper = 0 !< mapper counter
53  type(prism_mapper_type) ,public, pointer :: prism_mapper(:) !< list of defined mappers
54 
55  !--- local kinds ---
56 
57  integer,private,parameter :: r8 = ip_double_p
58  integer,private,parameter :: in = ip_i4_p
59 
60  !--- variable assocaited with local data buffers and reallocation
61 
62  real(R8),private,allocatable :: snew(:,:),sold(:,:) ! reals
63  integer,private,allocatable :: rnew(:),rold(:) ! ints
64  integer,private,allocatable :: cnew(:),cold(:) ! ints
65 
66 !------------------------------------------------------------
67 CONTAINS
68 !------------------------------------------------------------
69 
70 !------------------------------------------------------------
71 
72 !> Routine to generate mapping weights data via a direct SCRIP call
73 
74 !> This routine reads in grid data from files and passes that data
75 !> to SCRIP. Mapping weights are generated and written to a file.
76 !> This entire operation is done on a single task.
77 
78  SUBROUTINE oasis_map_genmap(mapid,namid)
79 
80  IMPLICIT NONE
81 
82  integer(ip_i4_p), intent(in) :: mapid !< map id
83  integer(ip_i4_p), intent(in) :: namid !< namcouple id
84  !----------------------------------------------------------
85 
86  integer(ip_i4_p) :: src_size,src_rank, ncrn_src
87  integer(ip_i4_p) ,allocatable :: src_dims(:),src_mask(:)
88  real(ip_double_p),allocatable :: src_lon(:),src_lat(:)
89  real(ip_double_p),allocatable :: src_corner_lon(:,:),src_corner_lat(:,:)
90  integer(ip_i4_p) :: dst_size,dst_rank, ncrn_dst
91  integer(ip_i4_p) ,allocatable :: dst_dims(:),dst_mask(:)
92  real(ip_double_p),allocatable :: dst_lon(:),dst_lat(:)
93  real(ip_double_p),allocatable :: dst_corner_lon(:,:),dst_corner_lat(:,:)
94  integer(ip_i4_p) ,allocatable :: ifld2(:,:)
95  real(ip_double_p),allocatable :: fld2(:,:),fld3(:,:,:)
96  integer(ip_i4_p) :: i,j,k,icnt,nx,ny,nc
97  logical :: lextrapdone
98  logical :: do_corners
99  character(len=ic_med) :: filename
100  character(len=ic_med) :: fldname
101  character(len=*),parameter :: subname = '(oasis_map_genmap)'
102 
103  call oasis_debug_enter(subname)
104 
105  lextrapdone = .false.
106  nx = -1 ! must be read
107  ny = -1 ! must be read
108  nc = 1 ! might not be read, set to something reasonable
109 
110  !--- checks first ---
111 
112  if (trim(namscrtyp(namid)) /= 'SCALAR') then
113  write(nulprt,*) subname,estr,'only scrip type SCALAR mapping supported'
114  CALL oasis_abort()
115  endif
116 
117 
118  do_corners = .false.
119  if (trim(namscrmet(namid)) == 'CONSERV') then
120  do_corners = .true.
121  endif
122 
123  !--- src data ---
124 
125  filename = 'grids.nc'
126  if (do_corners) then
127  fldname = trim(namsrcgrd(namid))//'.clo'
128  call oasis_io_read_field_fromroot(filename,fldname,nx=nx,ny=ny,nz=nc)
129  else
130  fldname = trim(namsrcgrd(namid))//'.lon'
131  call oasis_io_read_field_fromroot(filename,fldname,nx=nx,ny=ny)
132  endif
133  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',&
134  trim(fldname),nx,ny,nc,do_corners
135  src_rank = 2
136  src_size = nx*ny
137  allocate(src_dims(src_rank))
138  src_dims(1) = nx
139  src_dims(2) = ny
140  ncrn_src = nc
141  allocate(src_mask(src_size))
142  allocate(src_lon(src_size))
143  allocate(src_lat(src_size))
144  allocate(src_corner_lon(ncrn_src,src_size))
145  allocate(src_corner_lat(ncrn_src,src_size))
146 
147  allocate(ifld2(nx,ny))
148  filename = 'masks.nc'
149  fldname = trim(namsrcgrd(namid))//'.msk'
150  call oasis_io_read_field_fromroot(filename,fldname,ifld2=ifld2)
151  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
152  src_mask(icnt) = ifld2(i,j)
153  enddo; enddo
154  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
155  minval(src_mask),maxval(src_mask)
156  deallocate(ifld2)
157 
158  allocate(fld2(nx,ny))
159  filename = 'grids.nc'
160  fldname = trim(namsrcgrd(namid))//'.lon'
161  call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
162  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
163  src_lon(icnt) = fld2(i,j)
164  enddo; enddo
165  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
166  minval(src_lon),maxval(src_lon)
167  fldname = trim(namsrcgrd(namid))//'.lat'
168  call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
169  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
170  src_lat(icnt) = fld2(i,j)
171  enddo; enddo
172  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
173  minval(src_lat),maxval(src_lat)
174  deallocate(fld2)
175 
176  if (do_corners) then
177  allocate(fld3(nx,ny,nc))
178  filename = 'grids.nc'
179  fldname = trim(namsrcgrd(namid))//'.clo'
180  call oasis_io_read_field_fromroot(filename,fldname,fld3=fld3)
181  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
182  do k = 1,nc
183  src_corner_lon(k,icnt) = fld3(i,j,k)
184  enddo
185  enddo; enddo
186  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
187  minval(src_corner_lon),maxval(src_corner_lon)
188  fldname = trim(namsrcgrd(namid))//'.cla'
189  call oasis_io_read_field_fromroot(filename,fldname,fld3=fld3)
190  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
191  do k = 1,nc
192  src_corner_lat(k,icnt) = fld3(i,j,k)
193  enddo
194  enddo; enddo
195  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
196  minval(src_corner_lat),maxval(src_corner_lat)
197  deallocate(fld3)
198  else
199  src_corner_lon = -9999.
200  src_corner_lat = -9999.
201  endif
202 
203  !--- dst data ---
204 
205  filename = 'grids.nc'
206  if (do_corners) then
207  fldname = trim(namdstgrd(namid))//'.clo'
208  call oasis_io_read_field_fromroot(filename,fldname,nx=nx,ny=ny,nz=nc)
209  else
210  fldname = trim(namdstgrd(namid))//'.lon'
211  call oasis_io_read_field_fromroot(filename,fldname,nx=nx,ny=ny)
212  endif
213  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname),nx,ny,nc
214  dst_rank = 2
215  dst_size = nx*ny
216  allocate(dst_dims(dst_rank))
217  dst_dims(1) = nx
218  dst_dims(2) = ny
219  ncrn_dst = nc
220  allocate(dst_mask(dst_size))
221  allocate(dst_lon(dst_size))
222  allocate(dst_lat(dst_size))
223  allocate(dst_corner_lon(ncrn_dst,dst_size))
224  allocate(dst_corner_lat(ncrn_dst,dst_size))
225 
226  allocate(ifld2(nx,ny))
227  filename = 'masks.nc'
228  fldname = trim(namdstgrd(namid))//'.msk'
229  call oasis_io_read_field_fromroot(filename,fldname,ifld2=ifld2)
230  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
231  dst_mask(icnt) = ifld2(i,j)
232  enddo; enddo
233  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
234  minval(dst_mask),maxval(dst_mask)
235  deallocate(ifld2)
236 
237  allocate(fld2(nx,ny))
238  filename = 'grids.nc'
239  fldname = trim(namdstgrd(namid))//'.lon'
240  call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
241  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
242  dst_lon(icnt) = fld2(i,j)
243  enddo; enddo
244  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
245  minval(dst_lon),maxval(dst_lon)
246  fldname = trim(namdstgrd(namid))//'.lat'
247  call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
248  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
249  dst_lat(icnt) = fld2(i,j)
250  enddo; enddo
251  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
252  minval(dst_lat),maxval(dst_lat)
253  deallocate(fld2)
254 
255  if (do_corners) then
256  allocate(fld3(nx,ny,nc))
257  filename = 'grids.nc'
258  fldname = trim(namdstgrd(namid))//'.clo'
259  call oasis_io_read_field_fromroot(filename,fldname,fld3=fld3)
260  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
261  do k = 1,nc
262  dst_corner_lon(k,icnt) = fld3(i,j,k)
263  enddo
264  enddo; enddo
265  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
266  minval(dst_corner_lon),maxval(dst_corner_lon)
267  fldname = trim(namdstgrd(namid))//'.cla'
268  call oasis_io_read_field_fromroot(filename,fldname,fld3=fld3)
269  icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
270  do k = 1,nc
271  dst_corner_lat(k,icnt) = fld3(i,j,k)
272  enddo
273  enddo; enddo
274  if (oasis_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
275  minval(dst_corner_lat),maxval(dst_corner_lat)
276  deallocate(fld3)
277  else
278  dst_corner_lon = -9999.
279  dst_corner_lat = -9999.
280  endif
281 
282  IF (oasis_debug >= 15) THEN
283  WRITE(nulprt,*) subname,' call grid_init '
284  CALL oasis_flush(nulprt)
285  ENDIF
286 
287  !--- 0/1 mask convention opposite in scrip vs oasis
288  src_mask = 1 - src_mask
289  dst_mask = 1 - dst_mask
290  call grid_init(namscrmet(namid),namscrres(namid),namscrbin(namid), &
291  src_size, dst_size, src_dims, dst_dims, &
292  src_rank, dst_rank, ncrn_src, ncrn_dst, &
293  src_mask, dst_mask, namsrcgrd(namid), namdstgrd(namid), &
294  src_lat, src_lon, dst_lat, dst_lon, &
295  src_corner_lat, src_corner_lon, &
296  dst_corner_lat, dst_corner_lon, &
297  logunit=nulprt)
298  if (oasis_debug >= 15) then
299  WRITE(nulprt,*) subname,' done grid_init '
300  CALL oasis_flush(nulprt)
301  ENDIF
302 
303  IF (oasis_debug >= 15) THEN
304  WRITE(nulprt,*) subname,' call scrip '
305  CALL oasis_flush(nulprt)
306  ENDIF
307  call scrip(prism_mapper(mapid)%file,prism_mapper(mapid)%file,namscrmet(namid), &
308  namscrnor(namid),lextrapdone,namscrvam(namid),namscrnbr(namid))
309  IF (oasis_debug >= 15) THEN
310  WRITE(nulprt,*) subname,' done scrip '
311  CALL oasis_flush(nulprt)
312  ENDIF
313 
314  deallocate(src_dims, dst_dims)
315  deallocate(src_mask)
316  deallocate(src_lon)
317  deallocate(src_lat)
318  deallocate(src_corner_lon)
319  deallocate(src_corner_lat)
320  deallocate(dst_mask)
321  deallocate(dst_lon)
322  deallocate(dst_lat)
323  deallocate(dst_corner_lon)
324  deallocate(dst_corner_lat)
325 
326 
327  call oasis_debug_exit(subname)
328 
329  END SUBROUTINE oasis_map_genmap
330 
331 !------------------------------------------------------------
332 !BOP ===========================================================================
333 !> Read in mapping matrix data from a SCRIP netCDF weights file
334 !
335 ! !IROUTINE: oasis_map_sMatReaddnc_orig - Do a distributed read of a NetCDF SCRIP file and
336 ! return weights in a distributed SparseMatrix
337 !
338 ! !DESCRIPTION:
339 !> Read in mapping matrix data from a SCRIP netCDF data file using
340 !> a low memory method and then scatter to all pes. Based on
341 !> the sMatReaddnc method from CESM1.0.3.
342 !>
343 ! !REMARKS:
344 !> This routine leverages gsmaps to determine scatter pattern.
345 !>
346 !> The scatter is implemented as a broadcast of all weights then a local
347 !> computation on each pe to determine with weights to keep based
348 !> on gsmap information.
349 !>
350 !> The algorithm to determine whether a weight belongs on a pe involves
351 !> creating a couple local arrays (lsstart and lscount) which are
352 !> the local values of start and length from the gsmap. These are
353 !> sorted via a bubble sort and then searched via a binary search
354 !> to check whether a global index is on the local pe.
355 !>
356 !> The local buffer sizes are estimated up front based on ngridcell/npes
357 !> plus 20% (search for 1.2 below). If the local buffer size fills up, then
358 !> the buffer is reallocated 50% larger (search for 1.5 below) and the fill
359 !> continues. The idea is to trade off memory reallocation and copy
360 !> with memory usage. 1.2 and 1.5 are arbitary, other values may
361 !> result in better performance.
362 !>
363 !> Once all the matrix weights have been read, the sMat is initialized,
364 !> the values from the buffers are copied in, and everything is deallocated.
365 !
366 ! !INTERFACE: -----------------------------------------------------------------
367 
368 subroutine oasis_map_smatreaddnc_orig(sMat,SgsMap,DgsMap,newdom, &
369  filename,mytask,mpicom,nwgts, &
370  areasrc,areadst,ni_i,nj_i,ni_o,nj_o )
371 
372 ! !USES:
373 
374  !--- local kinds ---
375  integer,parameter :: R8 = ip_double_p
376  integer,parameter :: IN = ip_i4_p
377 
378 ! !INPUT/OUTPUT PARAMETERS:
379 
380  type(mct_smat) ,intent(out),pointer :: sMat(:) !< mapping data
381  type(mct_gsmap) ,intent(in) ,target :: SgsMap !< src gsmap
382  type(mct_gsmap) ,intent(in) ,target :: DgsMap !< dst gsmap
383  character(*) ,intent(in) :: newdom !< type of sMat (src or dst)
384  !< src = rearrange and map (bfb), dst = map and rearrange (partial sums)
385  character(*) ,intent(in) :: filename!< netCDF file to read
386  integer(IN) ,intent(in) :: mytask !< processor id
387  integer(IN) ,intent(in) :: mpicom !< mpi communicator
388  integer(IN) ,intent(out) :: nwgts !< number of weights
389  type(mct_avect) ,intent(out), optional :: areasrc !< area of src grid from mapping file
390  type(mct_avect) ,intent(out), optional :: areadst !< area of dst grid from mapping file
391  integer(IN) ,intent(out), optional :: ni_i !< number of lons on input grid
392  integer(IN) ,intent(out), optional :: nj_i !< number of lats on input grid
393  integer(IN) ,intent(out), optional :: ni_o !< number of lons on output grid
394  integer(IN) ,intent(out), optional :: nj_o !< number of lats on output grid
395 
396 ! !EOP
397 
398  !--- local ---
399  integer :: n,m ! generic loop indicies
400  integer :: na ! size of source domain
401  integer :: nb ! size of destination domain
402  integer :: ns ! number of non-zero elements in matrix
403  integer :: ni,nj ! number of row and col in the matrix
404  integer :: igrow ! aVect index for matrix row
405  integer :: igcol ! aVect index for matrix column
406  integer :: iwgt ! aVect index for matrix element
407  integer :: iarea ! aVect index for area
408  integer :: rsize ! size of read buffer
409  integer :: cnt ! local num of wgts
410  integer :: cntold ! local num of wgts, previous read
411  integer :: start(1)! netcdf read
412  integer :: count(1)! netcdf read
413  integer :: start2(2)! netcdf read
414  integer :: count2(2)! netcdf read
415  integer :: bsize ! buffer size
416  integer :: nread ! number of reads
417  logical :: mywt ! does this weight belong on my pe
418  integer :: dims(2)
419 
420  !--- buffers for i/o ---
421  real(R8) ,allocatable :: rtemp(:) ! real temporary
422  real(R8) ,allocatable :: Sbuf(:,:) ! real weights
423  real(R8) ,allocatable :: remaps(:,:) ! real weights with num_wgts dim
424  integer,allocatable :: Rbuf(:) ! ints rows
425  integer,allocatable :: Cbuf(:) ! ints cols
426 
427  !--- variables associated with local computation of global indices
428  integer :: lsize ! size of local seg map
429  integer :: commsize ! size of local communicator
430  integer,allocatable :: lsstart(:) ! local seg map info
431  integer,allocatable :: lscount(:) ! local seg map info
432  type(mct_gsmap),pointer :: mygsmap ! pointer to one of the gsmaps
433  integer :: l1,l2 ! generice indices for sort
434  logical :: found ! for sort
435 
436  !--- variable assocaited with local data buffers and reallocation
437  real(R8) ,allocatable :: Snew(:,:),Sold(:,:) ! reals
438  integer,allocatable :: Rnew(:),Rold(:) ! ints
439  integer,allocatable :: Cnew(:),Cold(:) ! ints
440 
441  character,allocatable :: str(:) ! variable length char string
442  character(len=ic_long):: attstr ! netCDF attribute name string
443  integer :: status ! netCDF routine return code
444  integer :: fid ! netCDF file ID
445  integer :: vid ! netCDF variable ID
446  integer :: did ! netCDF dimension ID
447  !--- arbitrary size of read buffer, this is the chunk size weights reading
448  integer,parameter :: rbuf_size = 100000
449 
450  !--- global source and destination areas ---
451  type(mct_avect) :: areasrc0 ! area of src grid from mapping file
452  type(mct_avect) :: areadst0 ! area of src grid from mapping file
453 
454  character(*),parameter :: areaAV_field = 'aream'
455 
456  !--- formats ---
457  character(*),parameter :: subName = '(oasis_map_sMatReaddnc_orig)'
458 
459 !-------------------------------------------------------------------------------
460 !
461 !-------------------------------------------------------------------------------
462 
463  call oasis_debug_enter(subname)
464  call oasis_mpi_commsize(mpicom,commsize)
465  nwgts = -1
466  if (mytask == 0) then
467  if (oasis_debug >= 2) write(nulprt,*) subname," reading mapping matrix data decomposed..."
468 
469  !----------------------------------------------------------------------------
470  !> * Open and read the file SCRIP weights size on the root task
471  !----------------------------------------------------------------------------
472  if (oasis_debug >=2 ) write(nulprt,*) subname," * file name : ",trim(filename)
473  status = nf90_open(trim(filename),nf90_nowrite,fid)
474  if (status /= nf90_noerr) then
475  write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
476  WRITE(nulprt,*) subname,estr,'mapping file not found = ',trim(filename)
477  call oasis_abort()
478  endif
479 
480  !--- get matrix dimensions ----------
481 ! status = nf90_inq_dimid (fid, 'n_s', did) ! size of sparse matrix
482  status = nf90_inq_dimid(fid, 'num_links', did) ! size of sparse matrix
483  status = nf90_inquire_dimension(fid, did , len = ns)
484 ! status = nf90_inq_dimid (fid, 'n_a', did) ! size of input vector
485  status = nf90_inq_dimid(fid, 'src_grid_size', did) ! size of input vector
486  status = nf90_inquire_dimension(fid, did , len = na)
487 ! status = nf90_inq_dimid (fid, 'n_b', did) ! size of output vector
488  status = nf90_inq_dimid(fid, 'dst_grid_size', did) ! size of output vector
489  status = nf90_inquire_dimension(fid, did , len = nb)
490  status = nf90_inq_dimid(fid, 'num_wgts', did) ! size of output vector
491  status = nf90_inquire_dimension(fid, did , len = nwgts)
492 
493  if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
494 ! status = nf90_inq_dimid (fid, 'ni_a', did) ! number of lons in input grid
495 ! status = nf90_inquire_dimension(fid, did , len = ni_i)
496 ! status = nf90_inq_dimid (fid, 'nj_a', did) ! number of lats in input grid
497 ! status = nf90_inquire_dimension(fid, did , len = nj_i)
498 ! status = nf90_inq_dimid (fid, 'ni_b', did) ! number of lons in output grid
499 ! status = nf90_inquire_dimension(fid, did , len = ni_o)
500 ! status = nf90_inq_dimid (fid, 'nj_b', did) ! number of lats in output grid
501 ! status = nf90_inquire_dimension(fid, did , len = nj_o)
502  status = nf90_inq_varid(fid, 'src_grid_dims', vid)
503  status = nf90_get_var(fid, vid, dims)
504  ni_i = dims(1)
505  nj_i = dims(2)
506  status = nf90_inq_varid(fid, 'dst_grid_dims', vid)
507  status = nf90_get_var(fid, vid, dims)
508  ni_o = dims(1)
509  nj_o = dims(2)
510  end if
511 
512  if (oasis_debug >= 2) write(nulprt,*) subname," * matrix dims src x dst : ",na,' x',nb
513  if (oasis_debug >= 2) write(nulprt,*) subname," * number of non-zero elements: ",ns
514 
515  endif
516 
517  !----------------------------------------------------------------------------
518  !> * Read and load area_a on root task
519  !----------------------------------------------------------------------------
520  if (present(areasrc)) then
521  if (mytask == 0) then
522  call mct_avect_init(areasrc0,' ',areaav_field,na)
523 ! status = nf90_inq_varid (fid,'area_a',vid)
524  status = nf90_inq_varid(fid,'src_grid_area',vid)
525  IF (status /= nf90_noerr) THEN
526  WRITE(nulprt,*) subname,' nf90_strerrr = ',trim(nf90_strerror(status))
527  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
528  CALL oasis_flush(nulprt)
529  ENDIF
530  status = nf90_get_var(fid, vid, areasrc0%rAttr)
531  IF (status /= nf90_noerr) THEN
532  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
533  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
534  CALL oasis_flush(nulprt)
535  ENDIF
536  endif
537  call mct_avect_scatter(areasrc0, areasrc, sgsmap, 0, mpicom, status)
538  if (status /= 0) call mct_die(subname,"Error on scatter of areasrc0")
539  if (mytask == 0) then
540 ! if (present(dbug)) then
541 ! if (dbug > 2) then
542 ! write(nulprt,*) subName,'Size of src ',mct_aVect_lSize(areasrc0)
543 ! write(nulprt,*) subName,'min/max src ',minval(areasrc0%rAttr(1,:)),maxval(areasrc0%rAttr(1,:))
544 ! endif
545 ! end if
546  call mct_avect_clean(areasrc0)
547  end if
548  end if
549 
550  !----------------------------------------------------------------------------
551  !> * Read and load area_b on root task
552  !----------------------------------------------------------------------------
553  if (present(areadst)) then
554  if (mytask == 0) then
555  call mct_avect_init(areadst0,' ',areaav_field,nb)
556 ! status = nf90_inq_varid (fid,'area_b',vid)
557  status = nf90_inq_varid(fid,'dst_grid_area',vid)
558  IF (status /= nf90_noerr) THEN
559  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
560  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
561  CALL oasis_flush(nulprt)
562  ENDIF
563  status = nf90_get_var(fid, vid, areadst0%rAttr)
564  IF (status /= nf90_noerr) THEN
565  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
566  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
567  CALL oasis_flush(nulprt)
568  ENDIF
569  endif
570  call mct_avect_scatter(areadst0, areadst, dgsmap, 0, mpicom, status)
571  if (status /= 0) call mct_die(subname,"Error on scatter of areadst0")
572  if (mytask == 0) then
573 ! if (present(dbug)) then
574 ! if (dbug > 2) then
575 ! write(nulprt,*) subName,'Size of dst ',mct_aVect_lSize(areadst0)
576 ! write(nulprt,*) subName,'min/max dst ',minval(areadst0%rAttr(1,:)),maxval(areadst0%rAttr(1,:))
577 ! endif
578 ! end if
579  call mct_avect_clean(areadst0)
580  endif
581  endif
582 
583  !----------------------------------------------------------------------------
584  !> * Broadcast ni and nj if requested
585  !----------------------------------------------------------------------------
586  if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
587  call oasis_mpi_bcast(ni_i,mpicom,subname//" MPI in ni_i bcast")
588  call oasis_mpi_bcast(nj_i,mpicom,subname//" MPI in nj_i bcast")
589  call oasis_mpi_bcast(ni_o,mpicom,subname//" MPI in ni_o bcast")
590  call oasis_mpi_bcast(nj_o,mpicom,subname//" MPI in nj_o bcast")
591  end if
592 
593  !----------------------------------------------------------------------------
594  !> * Broadcast array sizes and allocate arrays for local storage
595  !----------------------------------------------------------------------------
596  call oasis_mpi_bcast(ns,mpicom,subname//" MPI in ns bcast")
597  call oasis_mpi_bcast(na,mpicom,subname//" MPI in na bcast")
598  call oasis_mpi_bcast(nb,mpicom,subname//" MPI in nb bcast")
599  call oasis_mpi_bcast(nwgts,mpicom,subname//" MPI in nwgts bcast")
600 
601  !--- setup local seg map, sorted
602  if (newdom == 'src') then
603  mygsmap => dgsmap
604  elseif (newdom == 'dst') then
605  mygsmap => sgsmap
606  else
607  write(nulprt,*) subname,estr,'invalid newdom value, expect src or dst = ',newdom
608  call oasis_abort()
609  endif
610  lsize = 0
611  do n = 1,size(mygsmap%start)
612  if (mygsmap%pe_loc(n) == mytask) then
613  lsize=lsize+1
614  endif
615  enddo
616  allocate(lsstart(lsize),lscount(lsize),stat=status)
617  if (status /= 0) call mct_perr_die(subname,':: allocate Lsstart',status)
618 
619  !----------------------------------------------------------------------------
620  !> * Initialize lsstart and lscount, the sorted list of local indices
621  !----------------------------------------------------------------------------
622  lsize = 0
623  do n = 1,size(mygsmap%start)
624  if (mygsmap%pe_loc(n) == mytask) then ! on my pe
625  lsize=lsize+1
626  found = .false.
627  l1 = 1
628  do while (.not.found .and. l1 < lsize) ! bubble sort copy
629  if (mygsmap%start(n) < lsstart(l1)) then
630  do l2 = lsize, l1+1, -1
631  lsstart(l2) = lsstart(l2-1)
632  lscount(l2) = lscount(l2-1)
633  enddo
634  found = .true.
635  else
636  l1 = l1 + 1
637  endif
638  enddo
639  lsstart(l1) = mygsmap%start(n)
640  lscount(l1) = mygsmap%length(n)
641  endif
642  enddo
643  do n = 1,lsize-1
644  if (lsstart(n) > lsstart(n+1)) then
645  write(nulprt,*) subname,estr,'lsstart not properly sorted'
646  call oasis_abort()
647  endif
648  enddo
649 
650  !----------------------------------------------------------------------------
651  !> * Compute the number of chunks to read, read size is rbuf_size
652  !----------------------------------------------------------------------------
653  rsize = min(rbuf_size,ns) ! size of i/o chunks
654  bsize = ((ns/commsize) + 1 ) * 1.2 ! local temporary buffer size
655  if (ns == 0) then
656  nread = 0
657  else
658  nread = (ns-1)/rsize + 1 ! num of reads to do
659  endif
660 
661  !----------------------------------------------------------------------------
662  !> * Allocate arrays for local weights plus row and column indices
663  !----------------------------------------------------------------------------
664  if (mytask == 0) then
665  allocate(remaps(nwgts,rsize),stat=status)
666  if (status /= 0) call mct_perr_die(subname,':: allocate remaps',status)
667  endif
668  allocate(smat(nwgts),stat=status)
669  if (status /= 0) call mct_perr_die(subname,':: allocate Smat',status)
670  allocate(sbuf(nwgts,rsize),rbuf(rsize),cbuf(rsize),stat=status)
671  if (status /= 0) call mct_perr_die(subname,':: allocate Sbuf',status)
672  allocate(snew(nwgts,bsize),cnew(bsize),rnew(bsize),stat=status)
673  if (status /= 0) call mct_perr_die(subname,':: allocate Snew1',status)
674 
675  !----------------------------------------------------------------------------
676  !> * Loop over the chunks of weights data
677  !----------------------------------------------------------------------------
678  cnt = 0
679  do n = 1,nread
680  start(1) = (n-1)*rsize + 1
681  count(1) = min(rsize,ns-start(1)+1)
682  start2(1) = 1
683  count2(1) = nwgts
684  start2(2) = start(1)
685  count2(2) = count(1)
686 
687  !----------------------------------------------------------------------------
688  !> * Read chunk of data on root pe
689  !----------------------------------------------------------------------------
690  if (mytask== 0) then
691 ! status = nf90_inq_varid (fid,'S' ,vid)
692  status = nf90_inq_varid(fid,'remap_matrix' ,vid)
693 ! status = nf90_get_var(fid,vid,start,count,Sbuf)
694  status = nf90_get_var(fid,vid,remaps,start2,count2)
695  sbuf(:,:) = remaps(:,:)
696  IF (status /= nf90_noerr) THEN
697  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
698  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
699  CALL oasis_flush(nulprt)
700  ENDIF
701 
702 ! status = nf90_inq_varid (fid,'row',vid)
703  status = nf90_inq_varid(fid,'dst_address',vid)
704  status = nf90_get_var(fid,vid,rbuf,start,count)
705  IF (status /= nf90_noerr) THEN
706  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
707  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
708  CALL oasis_flush(nulprt)
709  ENDIF
710 
711 ! status = nf90_inq_varid (fid,'col',vid)
712  status = nf90_inq_varid(fid,'src_address',vid)
713  status = nf90_get_var(fid,vid,cbuf,start,count)
714  IF (status /= nf90_noerr) THEN
715  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
716  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
717  CALL oasis_flush(nulprt)
718  ENDIF
719  endif
720 
721  !----------------------------------------------------------------------------
722  !> * Broadcast S, row, col to all tasks
723  !----------------------------------------------------------------------------
724 
725  call oasis_mpi_bcast(sbuf,mpicom,subname//" MPI in Sbuf bcast")
726  call oasis_mpi_bcast(rbuf,mpicom,subname//" MPI in Rbuf bcast")
727  call oasis_mpi_bcast(cbuf,mpicom,subname//" MPI in Cbuf bcast")
728 
729  !----------------------------------------------------------------------------
730  !> * Each task keeps only the data required
731  !----------------------------------------------------------------------------
732  do m = 1,count(1)
733  !--- should this weight be on my pe
734  if (newdom == 'src') then
735  mywt = check_myindex(rbuf(m),lsstart,lscount)
736  elseif (newdom == 'dst') then
737  mywt = check_myindex(cbuf(m),lsstart,lscount)
738  endif
739 
740  if (mywt) then
741  cntold = cnt
742  cnt = cnt + 1
743 
744  !----------------------------------------------------------------------------
745  !> * Reallocate local weights arrays if they need to be bigger
746  !----------------------------------------------------------------------------
747  if (cnt > bsize) then
748  !--- allocate old arrays and copy new into old
749  allocate(sold(1:nwgts,cntold),rold(cntold),cold(cntold),stat=status)
750  if (status /= 0) call mct_perr_die(subname,':: allocate old',status)
751  sold(1:nwgts,1:cntold) = snew(1:nwgts,1:cntold)
752  rold(1:cntold) = rnew(1:cntold)
753  cold(1:cntold) = cnew(1:cntold)
754 
755  !--- reallocate new to bigger size, increase buffer by 50% (arbitrary)
756  deallocate(snew,rnew,cnew,stat=status)
757  if (status /= 0) call mct_perr_die(subname,':: allocate new',status)
758  bsize = 1.5 * bsize
759  if (oasis_debug > 15) write(nulprt,*) subname,' reallocate bsize to ',bsize
760  allocate(snew(nwgts,bsize),rnew(bsize),cnew(bsize),stat=status)
761  if (status /= 0) call mct_perr_die(subname,':: allocate old',status)
762 
763  !--- copy data back into new
764  snew(1:nwgts,1:cntold) = sold(1:nwgts,1:cntold)
765  rnew(1:cntold) = rold(1:cntold)
766  cnew(1:cntold) = cold(1:cntold)
767  deallocate(sold,rold,cold,stat=status)
768  if (status /= 0) call mct_perr_die(subname,':: deallocate old',status)
769  endif
770 
771  snew(1:nwgts,cnt) = sbuf(1:nwgts,m)
772  rnew(cnt) = rbuf(m)
773  cnew(cnt) = cbuf(m)
774  endif
775  enddo ! count
776  enddo ! nread
777 
778  !----------------------------------------------------------------------------
779  !> * Clean up arrays
780  !----------------------------------------------------------------------------
781  if (mytask == 0) then
782  deallocate(remaps, stat=status)
783  if (status /= 0) call mct_perr_die(subname,':: deallocate remaps',status)
784  endif
785  deallocate(sbuf,rbuf,cbuf, stat=status)
786  if (status /= 0) call mct_perr_die(subname,':: deallocate Sbuf',status)
787 
788  !----------------------------------------------------------------------------
789  !> * Initialize the mct sMat data type
790  !----------------------------------------------------------------------------
791  ! mct_sMat_init must be given the number of rows and columns that
792  ! would be in the full matrix. Nrows= size of output vector=nb.
793  ! Ncols = size of input vector = na.
794  do n = 1,nwgts
795  call mct_smat_init(smat(n), nb, na, cnt)
796  enddo
797 
798  igrow = mct_smat_indexia(smat(1),'grow')
799  igcol = mct_smat_indexia(smat(1),'gcol')
800  iwgt = mct_smat_indexra(smat(1),'weight')
801 
802  if (cnt /= 0) then
803  do n = 1,nwgts
804  smat(n)%data%rAttr(iwgt ,1:cnt) = snew(n,1:cnt)
805  smat(n)%data%iAttr(igrow,1:cnt) = rnew(1:cnt)
806  smat(n)%data%iAttr(igcol,1:cnt) = cnew(1:cnt)
807  enddo
808  endif
809 
810  !----------------------------------------------------------------------------
811  !> * More clean up
812  !----------------------------------------------------------------------------
813  deallocate(snew,rnew,cnew, stat=status)
814  deallocate(lsstart,lscount,stat=status)
815  if (status /= 0) call mct_perr_die(subname,':: deallocate new',status)
816 
817  if (mytask == 0) then
818  status = nf90_close(fid)
819  IF (oasis_debug >= 2) THEN
820  WRITE(nulprt,*) subname," ... done reading file"
821  CALL oasis_flush(nulprt)
822  ENDIF
823  endif
824 
825  call oasis_debug_exit(subname)
826 
827 end subroutine oasis_map_smatreaddnc_orig
828 
829 !------------------------------------------------------------
830 ! !BOP ===========================================================================
831 !
832 !> Function that checks whether an index is part of a start and count list
833 !
834 !> Does a binary search on a sorted start and count list to determine
835 !> whether index is a value in the list. The list values consist of
836 !> the values start(i):start(i)+count(i)-1 for all i.
837 !
838 ! !DESCRIPTION:
839 ! Do a binary search to see if a value is contained in a list of
840 ! values. return true or false. starti must be monotonically
841 ! increasing, function does NOT check this.
842 !
843 ! !INTERFACE: -----------------------------------------------------------------
844 
845 logical function check_myindex(index,starti,counti)
846 
847 ! !USES:
848 
849  !--- local kinds ---
850  integer,parameter :: R8 = ip_double_p
851  integer,parameter :: IN = ip_i4_p
852 
853 ! !INPUT/OUTPUT PARAMETERS:
854 
855  integer(IN) :: index !< index to search
856  integer(IN) :: starti(:) !< start list
857  integer(IN) :: counti(:) !< count list
858 
859 ! !EOP
860 
861  !--- local ---
862  integer(IN) :: nl,nc,nr,ncprev
863  integer(IN) :: lsize
864  logical :: stopnow
865 
866  !--- formats ---
867  character(*),parameter :: subName = '(check_myindex) '
868 
869 !-------------------------------------------------------------------------------
870 !
871 !-------------------------------------------------------------------------------
872 
873 ! call oasis_debug_enter(subname)
874  check_myindex = .false.
875 
876  lsize = size(starti)
877  if (lsize < 1) then
878 ! call oasis_debug_exit(subname)
879  return
880  endif
881 
882  nl = 0
883  nr = lsize + 1
884  nc = (nl+nr)/2
885  stopnow = .false.
886  do while (.not.stopnow)
887  if (index < starti(nc)) then
888  nr = nc
889  elseif (index > (starti(nc) + counti(nc) - 1)) then
890  nl = nc
891  else
892  check_myindex = .true.
893 ! call oasis_debug_exit(subname)
894  return
895  endif
896  ncprev = nc
897  nc = (nl + nr)/2
898  if (nc == ncprev .or. nc < 1 .or. nc > lsize) stopnow = .true.
899  enddo
900 
901  check_myindex = .false.
902 
903 ! call oasis_debug_exit(subname)
904 
905 end function check_myindex
906 
907 !===============================================================================
908 !BOP ===========================================================================
909 !> Read in mapping matrix data from a SCRIP netCDF file using smart scatter (ceg)
910 !
911 ! !IROUTINE: oasis_map_sMatReaddnc_ceg - Do a distributed read of a NetCDF SCRIP file and
912 ! return weights in a distributed SparseMatrix
913 !
914 ! !DESCRIPTION:
915 !> Read in mapping matrix data from a SCRIP netCDF data file using
916 !> a low memory method and then scatter to all pes using a smart method
917 !> where only select data is sent to tasks. Based on the sMatReaddnc method
918 !> from CESM1.0.3
919 ! !REMARKS:
920 !> This routine leverages gsmaps to determine scatter pattern
921 !>
922 !> The scatter is implemented via the root task reading the data and
923 !> then determining which task gets which weights from the gsmap.
924 !> The root the sends specific data to each task.
925 !>
926 !> The algorithm to determine which task a weigth belongs to involves
927 !> checking the task ownership for a given global index.
928 !>
929 !> The local buffer sizes are estimated up front based on ngridcell/npes
930 !> plus 20% (see 1.2 below). If the local buffer size fills up, then
931 !> the buffer is reallocated 50% larger (see 1.5 below) and the fill
932 !> continues. The idea is to trade off memory reallocation and copy
933 !> with memory usage. 1.2 and 1.5 are arbitary, other values may
934 !> result in better performance.
935 !>
936 !> Once all the matrix weights have been read, the sMat is initialized,
937 !> the values from the buffers are copied in, and everything is deallocated.
938 !
939 ! !INTERFACE: -----------------------------------------------------------------
940 
941 subroutine oasis_map_smatreaddnc_ceg(sMat,SgsMap,DgsMap,newdom, &
942  filename,mytask,mpicom,nwgts, &
943  areasrc,areadst,ni_i,nj_i,ni_o,nj_o )
944 
945 ! !USES:
946 
947 ! !INPUT/OUTPUT PARAMETERS:
948 
949  type(mct_smat) ,intent(out),pointer :: sMat(:) !< mapping data
950  type(mct_gsmap) ,intent(in) ,target :: SgsMap !< src gsmap
951  type(mct_gsmap) ,intent(in) ,target :: DgsMap !< dst gsmap
952  character(*) ,intent(in) :: newdom !< type of sMat (src or dst)
953  !< src = rearrange and map (bfb), dst = map and rearrange (partial sums)
954  character(*) ,intent(in) :: filename!< netCDF file to read
955  integer(IN) ,intent(in) :: mytask !< processor id
956  integer(IN) ,intent(in) :: mpicom !< mpi communicator
957  integer(IN) ,intent(out) :: nwgts !< number of weights
958  type(mct_avect) ,intent(out), optional :: areasrc !< area of src grid from mapping file
959  type(mct_avect) ,intent(out), optional :: areadst !< area of dst grid from mapping file
960  integer(IN) ,intent(out), optional :: ni_i !< number of lons on input grid
961  integer(IN) ,intent(out), optional :: nj_i !< number of lats on input grid
962  integer(IN) ,intent(out), optional :: ni_o !< number of lons on output grid
963  integer(IN) ,intent(out), optional :: nj_o !< number of lats on output grid
964 
965 ! !EOP
966 
967  !--- local ---
968  integer :: n,m ! generic loop indicies
969  integer :: na ! size of source domain
970  integer :: nb ! size of destination domain
971  integer :: ns ! number of non-zero elements in matrix
972  integer :: ni,nj ! number of row and col in the matrix
973  integer :: igrow ! aVect index for matrix row
974  integer :: igcol ! aVect index for matrix column
975  integer :: iwgt ! aVect index for matrix element
976  integer :: iarea ! aVect index for area
977  integer :: rsize ! size of read buffer
978  integer :: cnt ! local num of wgts
979  integer :: start(1)! netcdf read
980  integer :: count(1)! netcdf read
981  integer :: start2(2)! netcdf read
982  integer :: count2(2)! netcdf read
983  integer :: bsize ! buffer size
984  integer :: nread ! number of reads
985  logical :: mywt ! does this weight belong on my pe
986  integer :: dims(2)
987 
988  !--- buffers for i/o ---
989  real(R8) ,allocatable :: rtemp(:) ! real temporary
990  real(R8) ,allocatable :: Sbuf(:,:) ! real weights
991  real(R8) ,allocatable :: remaps(:,:) ! real weights with num_wgts dim
992  integer,allocatable :: Rbuf(:) ! ints rows
993  integer,allocatable :: Cbuf(:) ! ints cols
994  real(R8),allocatable :: SReadData(:,:) !
995  integer,allocatable :: RReadData(:) !
996  integer,allocatable :: CReadData(:) !
997  integer,allocatable :: pesave(:) !
998  real(R8),allocatable :: SDistData(:,:) !
999  integer,allocatable :: RDistData(:) !
1000  integer,allocatable :: CDistData(:) !
1001 
1002  !--- variables associated with local computation of global indices
1003  integer :: commsize ! size of local communicator
1004  type(mct_gsmap),pointer :: mygsmap ! pointer to one of the gsmaps
1005  integer :: l1,l2,lsize ! generice indices for sort
1006  integer,allocatable :: lsstart(:) ! local seg map info
1007  integer,allocatable :: lscount(:) ! local seg map info
1008  integer,allocatable :: lspeloc(:) ! local seg map info
1009  logical :: found ! for sort
1010  integer :: pe ! Process ID of owning process
1011  integer :: reclen ! Length of Rbuf/Cbuf to be received
1012  integer :: ierr ! MPI error check
1013  integer,allocatable :: cntrs(:) ! Counters of number of elements stored for each processor
1014  integer :: mpistatus(mpi_status_size) ! MPI_Status object
1015 
1016  character,allocatable :: str(:) ! variable length char string
1017  character(len=ic_long):: attstr ! netCDF attribute name string
1018  integer :: status ! netCDF routine return code
1019  integer :: fid ! netCDF file ID
1020  integer :: vid ! netCDF variable ID
1021  integer :: did ! netCDF dimension ID
1022  !--- arbitrary size of read buffer, this is the chunk size weights reading
1023  integer,parameter :: rbuf_size = 100000
1024  integer :: dimbuffer(8)
1025 
1026  !--- global source and destination areas ---
1027  type(mct_avect) :: areasrc0 ! area of src grid from mapping file
1028  type(mct_avect) :: areadst0 ! area of src grid from mapping file
1029 
1030  character(*),parameter :: areaAV_field = 'aream'
1031 
1032  !--- formats ---
1033  character(*),parameter :: subName = '(oasis_map_sMatReaddnc_ceg)'
1034  character*80 :: fname
1035 
1036 !-------------------------------------------------------------------------------
1037 !
1038 !-------------------------------------------------------------------------------
1039 
1040  call oasis_debug_enter(subname)
1041  call oasis_mpi_commsize(mpicom,commsize)
1042  nwgts = -1
1043  if (mytask == 0) then
1044  if (oasis_debug >= 2) write(nulprt,*) subname," reading mapping matrix data decomposed..."
1045 
1046  !----------------------------------------------------------------------------
1047  !> * Open and read the file SCRIP weights size on the root task
1048  !----------------------------------------------------------------------------
1049  if (oasis_debug >=2 ) write(nulprt,*) subname," * file name : ",trim(filename)
1050  status = nf90_open(trim(filename),nf90_nowrite,fid)
1051  if (status /= nf90_noerr) then
1052  write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1053  WRITE(nulprt,*) subname,estr,'mapping file not found = ',trim(filename)
1054  call oasis_abort()
1055  endif
1056 
1057  !--- get matrix dimensions ----------
1058 ! status = nf90_inq_dimid (fid, 'n_s', did) ! size of sparse matrix
1059  status = nf90_inq_dimid(fid, 'num_links', did) ! size of sparse matrix
1060  status = nf90_inquire_dimension(fid, did , len = ns)
1061 ! status = nf90_inq_dimid (fid, 'n_a', did) ! size of input vector
1062  status = nf90_inq_dimid(fid, 'src_grid_size', did) ! size of input vector
1063  status = nf90_inquire_dimension(fid, did , len = na)
1064 ! status = nf90_inq_dimid (fid, 'n_b', did) ! size of output vector
1065  status = nf90_inq_dimid(fid, 'dst_grid_size', did) ! size of output vector
1066  status = nf90_inquire_dimension(fid, did , len = nb)
1067  status = nf90_inq_dimid(fid, 'num_wgts', did) ! size of output vector
1068  status = nf90_inquire_dimension(fid, did , len = nwgts)
1069 
1070  if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
1071 ! status = nf90_inq_dimid (fid, 'ni_a', did) ! number of lons in input grid
1072 ! status = nf90_inquire_dimension(fid, did , len = ni_i)
1073 ! status = nf90_inq_dimid (fid, 'nj_a', did) ! number of lats in input grid
1074 ! status = nf90_inquire_dimension(fid, did , len = nj_i)
1075 ! status = nf90_inq_dimid (fid, 'ni_b', did) ! number of lons in output grid
1076 ! status = nf90_inquire_dimension(fid, did , len = ni_o)
1077 ! status = nf90_inq_dimid (fid, 'nj_b', did) ! number of lats in output grid
1078 ! status = nf90_inquire_dimension(fid, did , len = nj_o)
1079  status = nf90_inq_varid(fid, 'src_grid_dims', vid)
1080  status = nf90_get_var(fid, vid, dims)
1081  ni_i = dims(1)
1082  nj_i = dims(2)
1083  status = nf90_inq_varid(fid, 'dst_grid_dims', vid)
1084  status = nf90_get_var(fid, vid, dims)
1085  ni_o = dims(1)
1086  nj_o = dims(2)
1087  end if
1088 
1089  if (oasis_debug >= 2) write(nulprt,*) subname," * matrix dims src x dst : ",na,' x',nb
1090  if (oasis_debug >= 2) write(nulprt,*) subname," * number of non-zero elements: ",ns
1091 
1092  endif
1093 
1094  !----------------------------------------------------------------------------
1095  !> * Read and load area_a on root task
1096  !----------------------------------------------------------------------------
1097  if (present(areasrc)) then
1098  if (mytask == 0) then
1099  call mct_avect_init(areasrc0,' ',areaav_field,na)
1100 ! status = nf90_inq_varid (fid,'area_a',vid)
1101  status = nf90_inq_varid(fid,'src_grid_area',vid)
1102  IF (status /= nf90_noerr) THEN
1103  WRITE(nulprt,*) subname,' nf90_strerrr = ',trim(nf90_strerror(status))
1104  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
1105  CALL oasis_flush(nulprt)
1106  ENDIF
1107  status = nf90_get_var(fid, vid, areasrc0%rAttr)
1108  IF (status /= nf90_noerr) THEN
1109  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1110  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
1111  CALL oasis_flush(nulprt)
1112  ENDIF
1113  endif
1114  call mct_avect_scatter(areasrc0, areasrc, sgsmap, 0, mpicom, status)
1115  if (status /= 0) call mct_die(subname,"Error on scatter of areasrc0")
1116  if (mytask == 0) then
1117 ! if (present(dbug)) then
1118 ! if (dbug > 2) then
1119 ! write(nulprt,*) subName,'min/max src ',minval(areasrc0%rAttr(1,:)),maxval(areasrc0%rAttr(1,:))
1120 ! endif
1121 ! end if
1122  call mct_avect_clean(areasrc0)
1123  end if
1124  end if
1125 
1126  !----------------------------------------------------------------------------
1127  !> * Read and load area_b on root task
1128  !----------------------------------------------------------------------------
1129  if (present(areadst)) then
1130  if (mytask == 0) then
1131  call mct_avect_init(areadst0,' ',areaav_field,nb)
1132 ! status = nf90_inq_varid (fid,'area_b',vid)
1133  status = nf90_inq_varid(fid,'dst_grid_area',vid)
1134  IF (status /= nf90_noerr) THEN
1135  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1136  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
1137  CALL oasis_flush(nulprt)
1138  ENDIF
1139  status = nf90_get_var(fid, vid, areadst0%rAttr)
1140  IF (status /= nf90_noerr) THEN
1141  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1142  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
1143  CALL oasis_flush(nulprt)
1144  ENDIF
1145  endif
1146  call mct_avect_scatter(areadst0, areadst, dgsmap, 0, mpicom, status)
1147  if (status /= 0) call mct_die(subname,"Error on scatter of areadst0")
1148  if (mytask == 0) then
1149 ! if (present(dbug)) then
1150 ! if (dbug > 2) then
1151 ! write(nulprt,*) subName,'min/max dst ',minval(areadst0%rAttr(1,:)),maxval(areadst0%rAttr(1,:))
1152 ! endif
1153 ! end if
1154  call mct_avect_clean(areadst0)
1155  endif
1156  endif
1157 
1158  !----------------------------------------------------------------------------
1159  !> * Broadcast ni and nj if requested
1160  !> * Broadcast array sizes and allocate arrays for local storage
1161  ! Replace 8 MPI_Bcasts with just one
1162  !----------------------------------------------------------------------------
1163  if (mpi_rank_local.eq.0) then
1164  dimbuffer(1) = ns
1165  dimbuffer(2) = na
1166  dimbuffer(3) = nb
1167  dimbuffer(4) = nwgts
1168  if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
1169  dimbuffer(5) = ni_i
1170  dimbuffer(6) = nj_i
1171  dimbuffer(7) = ni_o
1172  dimbuffer(8) = nj_o
1173  end if
1174  endif
1175  call oasis_mpi_bcast(dimbuffer,mpicom,subname//" MPI of dimbuffer")
1176  if (mpi_rank_local.ne.0) then
1177  ns = dimbuffer(1)
1178  na = dimbuffer(2)
1179  nb = dimbuffer(3)
1180  nwgts = dimbuffer(4)
1181  if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
1182  ni_i = dimbuffer(5)
1183  nj_i = dimbuffer(6)
1184  ni_o = dimbuffer(7)
1185  nj_o = dimbuffer(8)
1186  end if
1187  endif
1188 
1189 ! call oasis_mpi_bcast(ni_i,mpicom,subName//" MPI in ni_i bcast")
1190 ! call oasis_mpi_bcast(nj_i,mpicom,subName//" MPI in nj_i bcast")
1191 ! call oasis_mpi_bcast(ni_o,mpicom,subName//" MPI in ni_o bcast")
1192 ! call oasis_mpi_bcast(nj_o,mpicom,subName//" MPI in nj_o bcast")
1193 
1194  !--- setup local seg map, sorted
1195  if (newdom == 'src') then
1196  mygsmap => dgsmap
1197  elseif (newdom == 'dst') then
1198  mygsmap => sgsmap
1199  else
1200  write(nulprt,*) subname,estr,'invalid newdom value, expect src or dst = ',newdom
1201  call oasis_abort()
1202  endif
1203 
1204  !----------------------------------------------------------------------------
1205  !> * Compute the number of chunks to read, read size is rbuf_size
1206  !----------------------------------------------------------------------------
1207  rsize = min(rbuf_size,ns) ! size of i/o chunks
1208  bsize = ((ns/commsize) + 1 ) * 1.2 ! local temporary buffer size
1209  if (ns == 0) then
1210  nread = 0
1211  else
1212  nread = (ns-1)/rsize + 1 ! num of reads to do
1213  endif
1214 
1215  !----------------------------------------------------------------------------
1216  !> * Allocate arrays for local weights plus row and column indices
1217  !----------------------------------------------------------------------------
1218  allocate(smat(nwgts),stat=status)
1219  if (status /= 0) call mct_perr_die(subname,':: allocate Smat',status)
1220  allocate(sbuf(nwgts,rsize),rbuf(rsize),cbuf(rsize),stat=status)
1221  if (status /= 0) call mct_perr_die(subname,':: allocate Sbuf',status)
1222  allocate(snew(nwgts,bsize),cnew(bsize),rnew(bsize),stat=status)
1223  if (status /= 0) call mct_perr_die(subname,':: allocate Snew1',status)
1224 ! write(nulprt,*) subname,mpi_rank_local,rsize,rsize*nwgts
1225 
1226  cnt = 1
1227 
1228  !----------------------------------------------------------------------------
1229  !> * On the root task
1230  !----------------------------------------------------------------------------
1231  if (mytask == 0) then
1232 
1233  !----------------------------------------------------------------------------
1234  !> * Initialize lsstart and lscount, the sorted list of local indices
1235  !----------------------------------------------------------------------------
1236  lsize = size(mygsmap%start)
1237  allocate(lsstart(lsize),lscount(lsize),lspeloc(lsize),stat=status)
1238  if (status /= 0) call mct_perr_die(subname,':: allocate Lsstart',status)
1239 
1240  do n = 1,size(mygsmap%start)
1241  found = .false.
1242  l1 = 1
1243  do while (.not.found .and. l1 < n) ! bubble sort copy
1244  if (mygsmap%start(n) < lsstart(l1)) then
1245  do l2 = n, l1+1, -1
1246  lsstart(l2) = lsstart(l2-1)
1247  lscount(l2) = lscount(l2-1)
1248  lspeloc(l2) = lspeloc(l2-1)
1249  enddo
1250  found = .true.
1251  else
1252  l1 = l1 + 1
1253  endif
1254  enddo
1255  lsstart(l1) = mygsmap%start(n)
1256  lscount(l1) = mygsmap%length(n)
1257  lspeloc(l1) = mygsmap%pe_loc(n)
1258  enddo
1259  do n = 1,size(mygsmap%start)-1
1260  if (lsstart(n) > lsstart(n+1)) then
1261  write(nulprt,*) subname,estr,'lsstart not properly sorted'
1262  call oasis_abort()
1263  endif
1264  enddo
1265 
1266  !----------------------------------------------------------------------------
1267  !> * Allocate arrays for reading data on the root
1268  !----------------------------------------------------------------------------
1269  allocate(remaps(nwgts,rsize),stat=status)
1270  if (status /= 0) call mct_perr_die(subname,':: allocate remaps',status)
1271  allocate(sreaddata(nwgts,rsize),rreaddata(rsize),creaddata(rsize),stat=status)
1272  if (status /= 0) call mct_perr_die(subname,':: allocate SReadData',status)
1273  allocate(pesave(rsize),stat=status)
1274  if (status /= 0) call mct_perr_die(subname,':: allocate pesave',status)
1275  allocate(sdistdata(nwgts,rsize),rdistdata(rsize),cdistdata(rsize), stat=status)
1276  if (status /= 0) call mct_perr_die(subname,':: allocate SDistData',status)
1277  allocate(cntrs(0:mpi_size_local), stat=status) ! Purposefully allocating one extra
1278  if (status /= 0) call mct_perr_die(subname,':: allocate cntrs',status)
1279 
1280  !----------------------------------------------------------------------------
1281  !> * Loop over the chunks of weights data
1282  !----------------------------------------------------------------------------
1283  do n = 1,nread
1284  start(1) = (n-1)*rsize + 1
1285  count(1) = min(rsize,ns-start(1)+1)
1286  start2(1) = 1
1287  count2(1) = nwgts
1288  start2(2) = start(1)
1289  count2(2) = count(1)
1290 
1291  !----------------------------------------------------------------------------
1292  !> * Read chunk of data
1293  !----------------------------------------------------------------------------
1294 ! status = nf90_inq_varid (fid,'S' ,vid)
1295  status = nf90_inq_varid(fid,'remap_matrix' ,vid)
1296 ! status = nf90_get_var(fid,vid,start,count,Sbuf)
1297  status = nf90_get_var(fid,vid,remaps,start2,count2)
1298  sreaddata(:,:) = remaps(:,:)
1299  IF (status /= nf90_noerr) THEN
1300  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1301  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
1302  CALL oasis_flush(nulprt)
1303  ENDIF
1304 ! status = nf90_inq_varid (fid,'row',vid)
1305  status = nf90_inq_varid(fid,'dst_address',vid)
1306  status = nf90_get_var(fid,vid,rreaddata,start,count)
1307  IF (status /= nf90_noerr) THEN
1308  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1309  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
1310  CALL oasis_flush(nulprt)
1311  ENDIF
1312 
1313 ! status = nf90_inq_varid (fid,'col',vid)
1314  status = nf90_inq_varid(fid,'src_address',vid)
1315  status = nf90_get_var(fid,vid,creaddata,start,count)
1316  IF (status /= nf90_noerr) THEN
1317  WRITE(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1318  WRITE(nulprt,*) subname,'model :',compid,' proc :',mpi_rank_local
1319  CALL oasis_flush(nulprt)
1320  ENDIF
1321 
1322  ! Two stage process
1323  ! 1. Count how many to send to each process
1324  ! 2. Build the sending array up in process receiving order
1325 
1326  ! Initialize cntrs array for number of accumulations per process
1327  cntrs =0
1328  pesave = -99
1329 
1330  !----------------------------------------------------------------------------
1331  !> * Determine which process owns each weight and count them
1332  !> * Determine offsets in the array
1333  !----------------------------------------------------------------------------
1334  do m = 1,count(1)
1335  !--- which process owns this point?
1336  if (newdom == 'src') then
1337  pe = get_cegindex(rreaddata(m),lsstart,lscount,lspeloc)
1338  else if (newdom == 'dst') then
1339  pe = get_cegindex(creaddata(m),lsstart,lscount,lspeloc)
1340  endif
1341  pesave(m) = pe
1342 
1343  ! Copy data, incrementing index array, pe = 1 to mpi_size_local
1344  if (pe+1 < 1 .or. pe+1 > mpi_size_local) then
1345  write(nulprt,*) subname,wstr,'get_cegindex search error', m,count(1),creaddata(m),rreaddata(m),sreaddata(:,m)
1346  else
1347  cntrs(pe+1) = cntrs(pe+1) + 1 ! Note incrementing 1->noprocs rather than 0->noprocs-1
1348  endif
1349 
1350  end do
1351 
1352  cntrs(0) = 1
1353  do pe = 1, mpi_size_local
1354  cntrs(pe) = cntrs(pe-1) + cntrs(pe)
1355  end do
1356 
1357  !----------------------------------------------------------------------------
1358  !> * Determine which process owns each weight and fill arrays
1359  !----------------------------------------------------------------------------
1360  do m = 1,count(1)
1361  !--- which process owns this point?
1362  if (newdom == 'src') then
1363  pe = get_cegindex(rreaddata(m),lsstart,lscount,lspeloc)
1364  else if (newdom == 'dst') then
1365  pe = get_cegindex(creaddata(m),lsstart,lscount,lspeloc)
1366  endif
1367 
1368  pe = pesave(m)
1369 
1370  ! Copy data, incrementing index array, pe = 1 to mpi_size_local
1371  if (pe+1 < 1 .or. pe+1 > mpi_size_local) then
1372  ! skip it
1373  ! write(nulprt,*) subname,'get_cegindex search error2', m,count(1),CReadData(m),RReadData(m),SReadData(:,m)
1374  else
1375  ! Copy data, incrementing index array
1376  sdistdata(:,cntrs(pe)) = sreaddata(:,m)
1377  rdistdata(cntrs(pe)) = rreaddata(m)
1378  cdistdata(cntrs(pe)) = creaddata(m)
1379 
1380  cntrs(pe) = cntrs(pe) + 1
1381  endif
1382 
1383  enddo ! count
1384 
1385  !----------------------------------------------------------------------------
1386  !> * Send select data from root to other processes
1387  ! Copy root process data into local array
1388  ! Send data from root to other processes
1389  !----------------------------------------------------------------------------
1390  if (cntrs(0).gt.1) then ! Need to do it differently if it is for me!
1391  reclen = cntrs(0)-1
1392  ! Reallocate local weights arrays if they need to be bigger
1393  call augment_arrays(cnt, reclen, bsize, nwgts)
1394 
1395  !--- now p0 copies straight into the ?new arrays
1396  snew(1:nwgts,cnt:cnt+reclen-1) = sdistdata(1:nwgts,1:reclen)
1397  rnew(cnt:cnt+reclen-1) = rdistdata(1:reclen)
1398  cnew(cnt:cnt+reclen-1) = cdistdata(1:reclen)
1399  cnt = cnt + reclen
1400 
1401  endif
1402 
1403  do pe = 1, mpi_size_local
1404  if (cntrs(pe).gt.cntrs(pe-1)) then
1405  ! Send data out
1406  m = cntrs(pe)-cntrs(pe-1)
1407 ! write(nulprt,*) subname, 'Sending [to, datalen] ', pe, m, cntrs(pe-1), cntrs(pe)
1408  call mpi_send(m, 1, mpi_integer, pe, 4000, mpicom, ierr)
1409  call mpi_send(sdistdata(1,cntrs(pe-1)), nwgts*m, mpi_real8, pe, 1000, mpicom, ierr)
1410  call mpi_send(rdistdata(cntrs(pe-1)), m, mpi_integer, pe, 2000, mpicom, ierr)
1411  call mpi_send(cdistdata(cntrs(pe-1)), m, mpi_integer, pe, 3000, mpicom, ierr)
1412  endif
1413  enddo
1414 
1415  enddo ! nread
1416 
1417  !----------------------------------------------------------------------------
1418  !> * Deallocate memory on root process
1419  !----------------------------------------------------------------------------
1420  m = -1
1421  do pe = 1, mpi_size_local-1
1422  !write(nulprt,*) subname, 'Final sending [to, datalen] ', m, cntrs(m)
1423  call mpi_send(m, 1, mpi_integer, pe, 4000, mpicom, ierr)
1424  end do
1425 
1426  ! Free memory used during the distribution
1427  deallocate(lsstart,lscount,lspeloc, stat=status)
1428  if (status /= 0) call mct_perr_die(subname,':: deallocate lsstart',status)
1429  deallocate(pesave, stat=status)
1430  if (status /= 0) call mct_perr_die(subname,':: deallocate pesave',status)
1431  deallocate(sreaddata,rreaddata,creaddata, stat=status)
1432  if (status /= 0) call mct_perr_die(subname,':: deallocate SReadData',status)
1433  deallocate(sdistdata,rdistdata,cdistdata, stat=status)
1434  if (status /= 0) call mct_perr_die(subname,':: deallocate SDistData',status)
1435  deallocate(cntrs, stat=status)
1436  if (status /= 0) call mct_perr_die(subname,':: deallocate cntrs',status)
1437  deallocate(remaps, stat=status)
1438  if (status /= 0) call mct_perr_die(subname,':: deallocate remaps',status)
1439 
1440  !----------------------------------------------------------------------------
1441  !> * On non-root processes
1442  !----------------------------------------------------------------------------
1443 
1444  else ! mytask == 0
1445 
1446  !----------------------------------------------------------------------------
1447  !> * Receive data from root
1448  !----------------------------------------------------------------------------
1449  call oasis_mpi_recv(reclen, 0, 4000, mpicom, subname//" MPI in reclen recv")
1450  ! Check we haven't now had the last data
1451  do while (reclen.ne.-1)
1452 
1453  !--- recv S, row, col on all other pes
1454  !write(nulprt,*) subname, 'Receiving [global-id, datalen] ', mpi_rank_global, reclen
1455  call mpi_recv(sbuf, reclen*nwgts, mpi_real8, 0, 1000, mpicom, mpi_status_ignore, ierr)
1456  call mpi_recv(rbuf, reclen, mpi_integer, 0, 2000, mpicom, mpi_status_ignore, ierr)
1457  call mpi_recv(cbuf, reclen, mpi_integer, 0, mpi_any_tag, mpicom, mpistatus, ierr)
1458 
1459  !----------------------------------------------------------------------------
1460  !> * Reallocate local weights arrays if they need to be bigger
1461  !----------------------------------------------------------------------------
1462  call augment_arrays(cnt, reclen, bsize, nwgts)
1463 
1464  !--- now each pe keeps what it has been sent
1465  snew(1:nwgts,cnt:cnt+reclen-1) = sbuf(1:nwgts,1:reclen)
1466  rnew(cnt:cnt+reclen-1) = rbuf(1:reclen)
1467  cnew(cnt:cnt+reclen-1) = cbuf(1:reclen)
1468  cnt = cnt + reclen
1469 
1470  call oasis_mpi_recv(reclen, 0, 4000, mpicom, subname//" MPI in reclen recv")
1471  end do
1472 
1473  endif ! mytask == 0
1474 
1475  ! Fix cnt to be the length of the array
1476  cnt = cnt-1
1477 
1478  !----------------------------------------------------------------------------
1479  !> * Clean up arrays
1480  !----------------------------------------------------------------------------
1481  deallocate(sbuf,rbuf,cbuf, stat=status)
1482  if (status /= 0) call mct_perr_die(subname,':: deallocate Sbuf',status)
1483 
1484  !----------------------------------------------------------------------------
1485  !> * Initialize the mct sMat data type
1486  !----------------------------------------------------------------------------
1487  ! mct_sMat_init must be given the number of rows and columns that
1488  ! would be in the full matrix. Nrows= size of output vector=nb.
1489  ! Ncols = size of input vector = na.
1490  do n = 1,nwgts
1491  call mct_smat_init(smat(n), nb, na, cnt)
1492  enddo
1493 
1494  igrow = mct_smat_indexia(smat(1),'grow')
1495  igcol = mct_smat_indexia(smat(1),'gcol')
1496  iwgt = mct_smat_indexra(smat(1),'weight')
1497 
1498  if (cnt /= 0) then
1499  do n = 1,nwgts
1500  smat(n)%data%rAttr(iwgt ,1:cnt) = snew(n,1:cnt)
1501  smat(n)%data%iAttr(igrow,1:cnt) = rnew(1:cnt)
1502  smat(n)%data%iAttr(igcol,1:cnt) = cnew(1:cnt)
1503  enddo
1504  endif
1505 
1506  !----------------------------------------------------------------------------
1507  !> * More clean up
1508  !----------------------------------------------------------------------------
1509  deallocate(snew,rnew,cnew, stat=status)
1510  if (status /= 0) call mct_perr_die(subname,':: deallocate new',status)
1511 
1512  if (mytask == 0) then
1513  status = nf90_close(fid)
1514  IF (oasis_debug >= 2) THEN
1515  WRITE(nulprt,*) subname," ... done reading file"
1516  CALL oasis_flush(nulprt)
1517  ENDIF
1518  endif
1519 
1520  call oasis_debug_exit(subname)
1521 
1522 end subroutine oasis_map_smatreaddnc_ceg
1523 
1524 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1525 !!!!!!!!!!!!!!!!!!!!!!!!! Extend array and store data !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1526 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1527 
1528 subroutine augment_arrays(cnt, reclen, bsize, nwgts)
1529 
1530 ! !USES:
1531 
1532  integer, intent(inout) :: cnt
1533  integer, intent(in) :: reclen
1534  integer, intent(inout) :: bsize
1535  integer, intent(in) :: nwgts
1536 
1537  integer :: status ! netCDF routine return code
1538  !--- formats ---
1539  character(*),parameter :: subName = '(ceg_coupler_augment_arrays)'
1540 
1541 
1542  !--- ?new arrays need to be bigger?
1543  if (cnt+reclen > bsize) then
1544 
1545 ! write(nulprt,*) subname,mpi_rank_global, 'EXTEND'
1546 
1547  !--- allocate old arrays and copy new into old
1548  allocate(sold(1:nwgts,cnt),rold(cnt),cold(cnt),stat=status)
1549  if (status /= 0) call mct_perr_die(subname,':: allocate old',status)
1550  sold(1:nwgts,1:cnt-1) = snew(1:nwgts,1:cnt-1)
1551  rold(1:cnt-1) = rnew(1:cnt-1)
1552  cold(1:cnt-1) = cnew(1:cnt-1)
1553 
1554  !--- reallocate new to bigger size, increase buffer by 50% (arbitrary)
1555  deallocate(snew,rnew,cnew,stat=status)
1556  if (status /= 0) call mct_perr_die(subname,':: allocate new',status)
1557  bsize = 1.5 * (cnt+reclen)
1558  if (oasis_debug > 15) write(nulprt,*) subname,' reallocate bsize to ',bsize
1559  allocate(snew(nwgts,bsize),rnew(bsize),cnew(bsize),stat=status)
1560  if (status /= 0) call mct_perr_die(subname,':: allocate old',status)
1561 
1562  !--- copy data back into new
1563  snew(1:nwgts,1:cnt-1) = sold(1:nwgts,1:cnt-1)
1564  rnew(1:cnt-1) = rold(1:cnt-1)
1565  cnew(1:cnt-1) = cold(1:cnt-1)
1566  deallocate(sold,rold,cold,stat=status)
1567  if (status /= 0) call mct_perr_die(subname,':: deallocate old',status)
1568  endif
1569 
1570 end subroutine augment_arrays
1571 
1572 !------------------------------------------------------------
1573 ! !BOP ===========================================================================
1574 !
1575 ! !IROUTINE: get_cegindex - binary search for index in list
1576 !
1577 ! !DESCRIPTION:
1578 ! Do a binary search to see if a value is contained in a list of
1579 ! values. return value of processor that owns it
1580 ! starti must be monotonically
1581 ! increasing, function does NOT check this.
1582 !
1583 ! !INTERFACE: -----------------------------------------------------------------
1584 
1585 integer function get_cegindex(index,starti,counti,peloci)
1586 
1587 ! !USES:
1588 
1589 
1590  !--- local kinds ---
1591  integer,parameter :: R8 = ip_double_p
1592  integer,parameter :: IN = ip_i4_p
1593 
1594 ! !INPUT/OUTPUT PARAMETERS:
1595 
1596  integer(IN) :: index !< index to search
1597  integer(IN) :: starti(:) !< start list
1598  integer(IN) :: counti(:) !< count list
1599  integer(IN) :: peloci(:) !< pe list
1600 
1601 ! !EOP
1602 
1603 
1604  !--- local ---
1605  integer(IN) :: nl,nc,nr,ncprev
1606  integer(IN) :: lsize
1607  logical :: stopnow
1608 
1609  !--- formats ---
1610  character(*),parameter :: subName = '(get_cegindex) '
1611 
1612 !-------------------------------------------------------------------------------
1613 !
1614 !-------------------------------------------------------------------------------
1615 ! call oasis_debug_enter(subname)
1616 
1617  get_cegindex = -99
1618  lsize = size(starti)
1619  if (lsize < 1) then
1620 ! call oasis_debug_exit(subname)
1621  return
1622  endif
1623 
1624  nl = 0
1625  nr = lsize + 1
1626  nc = (nl+nr)/2
1627  stopnow = .false.
1628  do while (.not.stopnow)
1629  if (index < starti(nc)) then
1630  nr = nc
1631  elseif (index > (starti(nc) + counti(nc) - 1)) then
1632  nl = nc
1633  else
1634  get_cegindex = peloci(nc)
1635 ! call oasis_debug_exit(subname)
1636  return
1637  endif
1638  ncprev = nc
1639  nc = (nl + nr)/2
1640  if (nc == ncprev .or. nc < 1 .or. nc > lsize) stopnow = .true.
1641  enddo
1642 
1643 
1644 ! call oasis_debug_exit(subname)
1645 
1646 
1647 end function get_cegindex
1648 
1649 !===============================================================================
1650 END MODULE mod_oasis_map
1651 
1652 
subroutine, public oasis_debug_exit(string)
Used when a subroutine is exited, write info to log file at some debug level.
OASIS partition data and methods.
System type methods.
Reads the namcouple file for use in OASIS.
Provides a generic and simpler interface into MPI calls for OASIS.
Generic overloaded interface into MPI broadcast.
subroutine, public oasis_map_genmap(mapid, namid)
Routine to generate mapping weights data via a direct SCRIP call.
subroutine, public oasis_mpi_commsize(comm, size, string)
Get the total number of tasks associated with a communicator.
Provides a common location for several OASIS variables.
OASIS map (interpolation) data and methods.
subroutine, public oasis_debug_enter(string)
Used when a subroutine is entered, write info to log file at some debug level.
Defines kinds for OASIS.
subroutine, public oasis_flush(nu)
Flushes output to file.
subroutine augment_arrays(cnt, reclen, bsize, nwgts)
Character string manipulation methods.
Performance timer methods.
subroutine, public oasis_abort(id_compid, cd_routine, cd_message)
OASIS abort method, publically available to users.
logical function check_myindex(index, starti, counti)
Function that checks whether an index is part of a start and count list.
subroutine, public oasis_map_smatreaddnc_ceg(sMat, SgsMap, DgsMap, newdom, fileName, mytask, mpicom, nwgts, areasrc, areadst, ni_i, nj_i, ni_o, nj_o)
Read in mapping matrix data from a SCRIP netCDF file using smart scatter (ceg)
subroutine, public oasis_map_smatreaddnc_orig(sMat, SgsMap, DgsMap, newdom, fileName, mytask, mpicom, nwgts, areasrc, areadst, ni_i, nj_i, ni_o, nj_o)
Read in mapping matrix data from a SCRIP netCDF weights file.
subroutine, public oasis_io_read_field_fromroot(filename, fldname, ifld2, fld2, fld3, nx, ny, nz)
Read a field on the root task from a file into an array.
Provides reusable IO routines for OASIS.
Definition: mod_oasis_io.F90:4
integer function get_cegindex(index, starti, counti, peloci)
Mapper data for interpolating data between grids.
OASIS variable data and methods.
Generic overloaded interface into MPI receive.
Defines parameters for OASIS.