Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

mod_oasis_map.F90 @ 6331

Last change on this file since 6331 was 6331, checked in by aclsce, 17 months ago
Moved oasis-mct_5.0 in oasis3-mct/branches directory.
File size: 83.7 KB

Line
1
2	!> OASIS map (interpolation) data and methods
3
4
5	MODULE mod_oasis_map
6
7	USE mod_oasis_kinds
8	USE mod_oasis_data
9	USE mod_oasis_parameters
10	USE mod_oasis_namcouple
11	USE mod_oasis_sys
12	USE mod_oasis_part
13	USE mod_oasis_var
14	USE mod_oasis_mpi
15	USE mod_oasis_string
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
27	public oasis_map_sMatReaddnc_orig
28	public oasis_map_sMatReaddnc_ceg
29
30	!--- Type of data
31
32	public prism_mapper_type
33
34	!> Mapper data for interpolating data between grids
35	type prism_mapper_type
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_long):: file2 !< alternate file to read/write
41	character(len=ic_med) :: loc !< location setting, src or dst model
42	character(len=ic_med) :: opt !< optimization setting, bfb, sum, or opt
43	character(len=ic_med) :: optval !< mct map setting, src or dst, derived from opt
44	logical :: init !< flag indicating initialization complete
45	integer(kind=ip_i4_p) :: spart !< src partition
46	integer(kind=ip_i4_p) :: dpart !< dst partition
47	character(len=ic_med) :: srcgrid !< source grid name
48	character(len=ic_med) :: dstgrid !< target grid name
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	logical, parameter :: local_timers_on = .false.
67
68	!------------------------------------------------------------
69	CONTAINS
70	!------------------------------------------------------------
71
72	!------------------------------------------------------------
73
74	!> Routine to generate mapping weights data via a direct SCRIP call
75
76	!> This routine reads in grid data from files and passes that data
77	!> to SCRIP. Mapping weights are generated and written to a file.
78	!> This entire operation is done on a single task.
79
80	SUBROUTINE oasis_map_genmap(mapid,namid)
81
82	IMPLICIT NONE
83
84	integer(ip_i4_p), intent(in) :: mapid !< map id
85	integer(ip_i4_p), intent(in) :: namid !< namcouple id
86	!----------------------------------------------------------
87
88	integer(ip_i4_p) :: src_size,src_rank, ncrn_src
89	integer(ip_i4_p) ,allocatable :: src_dims(:),src_mask(:)
90	real(ip_double_p),allocatable :: src_lon(:),src_lat(:)
91	real(ip_double_p),allocatable :: src_area(:)
92	real(ip_double_p),allocatable :: src_corner_lon(:,:),src_corner_lat(:,:)
93	integer(ip_i4_p) :: dst_size,dst_rank, ncrn_dst
94	integer(ip_i4_p) ,allocatable :: dst_dims(:),dst_mask(:)
95	real(ip_double_p),allocatable :: dst_lon(:),dst_lat(:)
96	real(ip_double_p),allocatable :: dst_area(:)
97	real(ip_double_p),allocatable :: dst_corner_lon(:,:),dst_corner_lat(:,:)
98	integer(ip_i4_p) ,allocatable :: ifld2(:,:)
99	real(ip_double_p),allocatable :: fld2(:,:),fld3(:,:,:)
100	integer(ip_i4_p) :: i,j,k,icnt,nx,ny,nc
101	logical :: lextrapdone
102	logical :: do_corners
103	logical :: ll_src_area_in
104	logical :: ll_dst_area_in
105	character(len=ic_med) :: filename
106	character(len=ic_med) :: fldname
107	character(len=*),parameter :: subname = '(oasis_map_genmap)'
108
109	call oasis_debug_enter(subname)
110
111	lextrapdone = .false.
112	nx = -1 ! must be read
113	ny = -1 ! must be read
114	nc = 1 ! might not be read, set to something reasonable
115
116	!--- checks first ---
117
118	if (trim(namscrtyp(namID)) /= 'SCALAR') then
119	write(nulprt,*) subname,estr,'only scrip type SCALAR mapping supported'
120	call oasis_abort(file=__FILE__,line=__LINE__)
121	endif
122
123	do_corners = .false.
124	if (trim(namscrmet(namID)) == 'CONSERV') then
125	do_corners = .true.
126	endif
127
128	!--- src data ---
129
130	filename = 'grids.nc'
131	if (do_corners) then
132	fldname = trim(namsrcgrd(namID))//'.clo'
133	call oasis_io_read_field_fromroot(filename,fldname,nx=nx,ny=ny,nz=nc)
134	else
135	fldname = trim(namsrcgrd(namID))//'.lon'
136	call oasis_io_read_field_fromroot(filename,fldname,nx=nx,ny=ny)
137	endif
138	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',&
139	trim(fldname),nx,ny,nc,do_corners
140	src_rank = 2
141	src_size = nx*ny
142	allocate(src_dims(src_rank))
143	src_dims(1) = nx
144	src_dims(2) = ny
145	ncrn_src = nc
146	allocate(src_mask(src_size))
147	allocate(src_lon (src_size))
148	allocate(src_lat (src_size))
149	allocate(src_area (src_size))
150	allocate(src_corner_lon(ncrn_src,src_size))
151	allocate(src_corner_lat(ncrn_src,src_size))
152
153	filename = 'masks.nc'
154	fldname = trim(namsrcgrd(namID))//'.msk'
155	allocate(ifld2(nx,ny))
156	call oasis_io_read_field_fromroot(filename,fldname,ifld2=ifld2)
157	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
158	src_mask(icnt) = ifld2(i,j)
159	enddo; enddo
160	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
161	minval(src_mask),maxval(src_mask)
162	deallocate(ifld2)
163
164	allocate(fld2(nx,ny))
165
166	filename = 'areas.nc'
167	fldname = trim(namsrcgrd(namID))//'.srf'
168	if (oasis_io_varexists(filename,fldname)) then
169	ll_src_area_in = .true.
170	call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
171	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
172	src_area(icnt) = fld2(i,j)
173	enddo; enddo
174	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
175	minval(src_area),maxval(src_area)
176	else
177	ll_src_area_in = .false.
178	src_area = -9999.
179	endif
180
181	filename = 'grids.nc'
182	fldname = trim(namsrcgrd(namID))//'.lon'
183	call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
184	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
185	src_lon(icnt) = fld2(i,j)
186	enddo; enddo
187	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
188	minval(src_lon),maxval(src_lon)
189
190	fldname = trim(namsrcgrd(namID))//'.lat'
191	call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
192	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
193	src_lat(icnt) = fld2(i,j)
194	enddo; enddo
195	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
196	minval(src_lat),maxval(src_lat)
197
198	deallocate(fld2)
199
200	if (do_corners) then
201	allocate(fld3(nx,ny,nc))
202	filename = 'grids.nc'
203	fldname = trim(namsrcgrd(namID))//'.clo'
204	call oasis_io_read_field_fromroot(filename,fldname,fld3=fld3)
205	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
206	do k = 1,nc
207	src_corner_lon(k,icnt) = fld3(i,j,k)
208	enddo
209	enddo; enddo
210	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
211	minval(src_corner_lon),maxval(src_corner_lon)
212	fldname = trim(namsrcgrd(namID))//'.cla'
213	call oasis_io_read_field_fromroot(filename,fldname,fld3=fld3)
214	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
215	do k = 1,nc
216	src_corner_lat(k,icnt) = fld3(i,j,k)
217	enddo
218	enddo; enddo
219	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
220	minval(src_corner_lat),maxval(src_corner_lat)
221	deallocate(fld3)
222	else
223	src_corner_lon = -9999.
224	src_corner_lat = -9999.
225	endif
226
227	!--- dst data ---
228
229	filename = 'grids.nc'
230	if (do_corners) then
231	fldname = trim(namdstgrd(namID))//'.clo'
232	call oasis_io_read_field_fromroot(filename,fldname,nx=nx,ny=ny,nz=nc)
233	else
234	fldname = trim(namdstgrd(namID))//'.lon'
235	call oasis_io_read_field_fromroot(filename,fldname,nx=nx,ny=ny)
236	endif
237	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname),nx,ny,nc
238	dst_rank = 2
239	dst_size = nx*ny
240	allocate(dst_dims(dst_rank))
241	dst_dims(1) = nx
242	dst_dims(2) = ny
243	ncrn_dst = nc
244	allocate(dst_mask(dst_size))
245	allocate(dst_lon (dst_size))
246	allocate(dst_lat (dst_size))
247	allocate(dst_area (dst_size))
248	allocate(dst_corner_lon(ncrn_dst,dst_size))
249	allocate(dst_corner_lat(ncrn_dst,dst_size))
250
251	filename = 'masks.nc'
252	fldname = trim(namdstgrd(namID))//'.msk'
253	allocate(ifld2(nx,ny))
254	call oasis_io_read_field_fromroot(filename,fldname,ifld2=ifld2)
255	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
256	dst_mask(icnt) = ifld2(i,j)
257	enddo; enddo
258	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
259	minval(dst_mask),maxval(dst_mask)
260	deallocate(ifld2)
261
262	allocate(fld2(nx,ny))
263
264	filename = 'areas.nc'
265	fldname = trim(namdstgrd(namID))//'.srf'
266	if (oasis_io_varexists(filename,fldname)) then
267	ll_dst_area_in = .true.
268	call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
269	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
270	dst_area(icnt) = fld2(i,j)
271	enddo; enddo
272	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
273	minval(dst_area),maxval(dst_area)
274	else
275	ll_dst_area_in = .false.
276	dst_area = -9999.
277	endif
278
279	filename = 'grids.nc'
280	fldname = trim(namdstgrd(namID))//'.lon'
281	call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
282	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
283	dst_lon(icnt) = fld2(i,j)
284	enddo; enddo
285	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
286	minval(dst_lon),maxval(dst_lon)
287
288	fldname = trim(namdstgrd(namID))//'.lat'
289	call oasis_io_read_field_fromroot(filename,fldname,fld2=fld2)
290	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
291	dst_lat(icnt) = fld2(i,j)
292	enddo; enddo
293	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
294	minval(dst_lat),maxval(dst_lat)
295
296	deallocate(fld2)
297
298	if (do_corners) then
299	allocate(fld3(nx,ny,nc))
300	filename = 'grids.nc'
301	fldname = trim(namdstgrd(namID))//'.clo'
302	call oasis_io_read_field_fromroot(filename,fldname,fld3=fld3)
303	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
304	do k = 1,nc
305	dst_corner_lon(k,icnt) = fld3(i,j,k)
306	enddo
307	enddo; enddo
308	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
309	minval(dst_corner_lon),maxval(dst_corner_lon)
310	fldname = trim(namdstgrd(namID))//'.cla'
311	call oasis_io_read_field_fromroot(filename,fldname,fld3=fld3)
312	icnt = 0; do j = 1,ny; do i = 1,nx; icnt = icnt + 1
313	do k = 1,nc
314	dst_corner_lat(k,icnt) = fld3(i,j,k)
315	enddo
316	enddo; enddo
317	if (OASIS_debug >= 15) write(nulprt,*) subname,' read ',trim(filename),' ',trim(fldname), &
318	minval(dst_corner_lat),maxval(dst_corner_lat)
319	deallocate(fld3)
320	else
321	dst_corner_lon = -9999.
322	dst_corner_lat = -9999.
323	endif
324
325	if (OASIS_debug >= 15) THEN
326	write(nulprt,*) subname,' call grid_init '
327	call oasis_flush(nulprt)
328	endif
329
330	!--- 0/1 mask convention opposite in scrip vs oasis
331	src_mask = 1 - src_mask
332	dst_mask = 1 - dst_mask
333	call grid_init(namscrmet(namID),namscrres(namID),namscrbin(namID), &
334	src_size, dst_size, src_dims, dst_dims, &
335	src_rank, dst_rank, ncrn_src, ncrn_dst, &
336	src_mask, dst_mask, namsrcgrd(namID), namdstgrd(namID), &
337	src_lat, src_lon, dst_lat, dst_lon, &
338	src_corner_lat, src_corner_lon, &
339	dst_corner_lat, dst_corner_lon, &
340	ll_src_area_in, src_area, &
341	ll_dst_area_in, dst_area, &
342	ilogunit=nulprt,ilogprt=OASIS_debug)
343	if (OASIS_debug >= 15) then
344	write(nulprt,*) subname,' done grid_init '
345	call oasis_flush(nulprt)
346	endif
347
348	if (OASIS_debug >= 15) THEN
349	write(nulprt,*) subname,' call scrip '
350	call oasis_flush(nulprt)
351	endif
352	if (local_timers_on) call oasis_timer_start('cpl_genmap_scrip')
353	call scrip(prism_mapper(mapid)%file,prism_mapper(mapid)%file,namscrmet(namID), &
354	namscrnor(namID),lextrapdone,namscrvam(namID),namscrnbr(namID),namscrord(namID), &
355	namscrnth(namID),namscrsth(namID), &
356	mpi_comm_map, mpi_size_map, mpi_rank_map, mpi_root_map, namcdftyp)
357	if (local_timers_on) call oasis_timer_stop('cpl_genmap_scrip')
358	if (OASIS_debug >= 15) THEN
359	write(nulprt,*) subname,' done scrip '
360	call oasis_flush(nulprt)
361	endif
362
363	deallocate(src_dims, dst_dims)
364	deallocate(src_mask)
365	deallocate(src_lon)
366	deallocate(src_lat)
367	deallocate(src_area)
368	deallocate(src_corner_lon)
369	deallocate(src_corner_lat)
370	deallocate(dst_mask)
371	deallocate(dst_lon)
372	deallocate(dst_lat)
373	deallocate(dst_area)
374	deallocate(dst_corner_lon)
375	deallocate(dst_corner_lat)
376
377
378	call oasis_debug_exit(subname)
379
380	END SUBROUTINE oasis_map_genmap
381
382	!------------------------------------------------------------
383	!BOP ===========================================================================
384	!> Read in mapping matrix data from a SCRIP netCDF weights file
385	!
386	! !IROUTINE: oasis_map_sMatReaddnc_orig - Do a distributed read of a NetCDF SCRIP file and
387	! return weights in a distributed SparseMatrix
388	!
389	! !DESCRIPTION:
390	!> Read in mapping matrix data from a SCRIP netCDF data file using
391	!> a low memory method and then scatter to all pes. Based on
392	!> the sMatReaddnc method from CESM1.0.3.
393	!>
394	! !REMARKS:
395	!> This routine leverages gsmaps to determine scatter pattern.
396	!>
397	!> The scatter is implemented as a broadcast of all weights then a local
398	!> computation on each pe to determine with weights to keep based
399	!> on gsmap information.
400	!>
401	!> The algorithm to determine whether a weight belongs on a pe involves
402	!> creating a couple local arrays (lsstart and lscount) which are
403	!> the local values of start and length from the gsmap. These are
404	!> sorted via a bubble sort and then searched via a binary search
405	!> to check whether a global index is on the local pe.
406	!>
407	!> The local buffer sizes are estimated up front based on ngridcell/npes
408	!> plus 20% (search for 1.2 below). If the local buffer size fills up, then
409	!> the buffer is reallocated 50% larger (search for 1.5 below) and the fill
410	!> continues. The idea is to trade off memory reallocation and copy
411	!> with memory usage. 1.2 and 1.5 are arbitary, other values may
412	!> result in better performance.
413	!>
414	!> Once all the matrix weights have been read, the sMat is initialized,
415	!> the values from the buffers are copied in, and everything is deallocated.
416	!
417	! !INTERFACE: -----------------------------------------------------------------
418
419	subroutine oasis_map_sMatReaddnc_orig(sMat,SgsMap,DgsMap,newdom, &
420	fileName,mytask,mpicom,nwgts, &
421	areasrc,areadst,ni_i,nj_i,ni_o,nj_o )
422
423	! !USES:
424
425	!--- local kinds ---
426	integer,parameter :: R8 = ip_double_p
427	integer,parameter :: IN = ip_i4_p
428
429	! !INPUT/OUTPUT PARAMETERS:
430
431	type(mct_sMat) ,intent(out),pointer :: sMat(:) !< mapping data
432	type(mct_gsMap) ,intent(in) ,target :: SgsMap !< src gsmap
433	type(mct_gSMap) ,intent(in) ,target :: DgsMap !< dst gsmap
434	character(*) ,intent(in) :: newdom !< type of sMat (src or dst)
435	!< src = rearrange and map (bfb), dst = map and rearrange (partial sums)
436	character(*) ,intent(in) :: filename!< netCDF file to read
437	integer(IN) ,intent(in) :: mytask !< processor id
438	integer(IN) ,intent(in) :: mpicom !< mpi communicator
439	integer(IN) ,intent(out) :: nwgts !< number of weights
440	type(mct_Avect) ,intent(out), optional :: areasrc !< area of src grid from mapping file
441	type(mct_Avect) ,intent(out), optional :: areadst !< area of dst grid from mapping file
442	integer(IN) ,intent(out), optional :: ni_i !< number of lons on input grid
443	integer(IN) ,intent(out), optional :: nj_i !< number of lats on input grid
444	integer(IN) ,intent(out), optional :: ni_o !< number of lons on output grid
445	integer(IN) ,intent(out), optional :: nj_o !< number of lats on output grid
446
447	! !EOP
448
449	!--- local ---
450	integer :: n,m ! generic loop indicies
451	integer :: na ! size of source domain
452	integer :: nb ! size of destination domain
453	integer :: ns ! number of non-zero elements in matrix
454	integer :: ni,nj ! number of row and col in the matrix
455	integer :: igrow ! aVect index for matrix row
456	integer :: igcol ! aVect index for matrix column
457	integer :: iwgt ! aVect index for matrix element
458	integer :: iarea ! aVect index for area
459	integer :: rsize ! size of read buffer
460	integer :: cnt ! local num of wgts
461	integer :: cntold ! local num of wgts, previous read
462	integer :: start(1)! netcdf read
463	integer :: count(1)! netcdf read
464	integer :: start2(2)! netcdf read
465	integer :: count2(2)! netcdf read
466	integer :: bsize ! buffer size
467	integer :: nread ! number of reads
468	logical :: mywt ! does this weight belong on my pe
469	integer :: dims(2)
470	logical :: abort_weight ! flag if bad weight found
471
472	!--- buffers for i/o ---
473	real(R8) ,allocatable :: rtemp(:) ! real temporary
474	real(R8) ,allocatable :: Sbuf(:,:) ! real weights
475	real(R8) ,allocatable :: remaps(:,:) ! real weights with num_wgts dim
476	integer,allocatable :: Rbuf(:) ! ints rows
477	integer,allocatable :: Cbuf(:) ! ints cols
478
479	!--- variables associated with local computation of global indices
480	integer :: lsize ! size of local seg map
481	integer :: commsize ! size of local communicator
482	integer,allocatable :: lsstart(:) ! local seg map info
483	integer,allocatable :: lscount(:) ! local seg map info
484	type(mct_gsMap),pointer :: mygsmap ! pointer to one of the gsmaps
485	integer :: l1,l2 ! generice indices for sort
486	logical :: found ! for sort
487
488	!--- variable assocaited with local data buffers and reallocation
489	real(R8) ,allocatable :: Snew(:,:),Sold(:,:) ! reals
490	integer,allocatable :: Rnew(:),Rold(:) ! ints
491	integer,allocatable :: Cnew(:),Cold(:) ! ints
492
493	character,allocatable :: str(:) ! variable length char string
494	character(len=ic_long):: attstr ! netCDF attribute name string
495	integer :: status ! netCDF routine return code
496	integer :: fid ! netCDF file ID
497	integer :: vid ! netCDF variable ID
498	integer :: did ! netCDF dimension ID
499	!--- arbitrary size of read buffer, this is the chunk size weights reading
500	integer,parameter :: rbuf_size = 100000
501
502	!--- global source and destination areas ---
503	type(mct_Avect) :: areasrc0 ! area of src grid from mapping file
504	type(mct_Avect) :: areadst0 ! area of src grid from mapping file
505
506	character(*),parameter :: areaAV_field = 'aream'
507
508	!--- formats ---
509	character(*),parameter :: subName = '(oasis_map_sMatReaddnc_orig)'
510
511	!-------------------------------------------------------------------------------
512	!
513	!-------------------------------------------------------------------------------
514
515	call oasis_debug_enter(subname)
516	call oasis_mpi_commsize(mpicom,commsize)
517	nwgts = -1
518	if (local_timers_on) call oasis_timer_start('map_read_orig_read')
519
520	if (mytask == 0) then
521	if (OASIS_debug >= 2) write(nulprt,*) subname," reading mapping matrix data decomposed..."
522
523	!----------------------------------------------------------------------------
524	!> * Open and read the file SCRIP weights size on the root task
525	!----------------------------------------------------------------------------
526
527	if (OASIS_debug >=2 ) write(nulprt,) subname," file name : ",trim(fileName)
528	status = nf90_open(trim(filename),NF90_NOWRITE,fid)
529	if (status /= NF90_NOERR) then
530	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
531	write(nulprt,*) subname,estr,'mapping file not found = ',trim(filename)
532	call oasis_abort(file=__FILE__,line=__LINE__)
533	endif
534
535	!--- get matrix dimensions ----------
536	! status = nf90_inq_dimid (fid, 'n_s', did) ! size of sparse matrix
537	status = nf90_inq_dimid (fid, 'num_links', did) ! size of sparse matrix
538	if (status /= NF90_NOERR) then
539	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
540	write(nulprt,*) subname,estr,'mapping file not found = ',trim(filename)
541	call oasis_abort(file=__FILE__,line=__LINE__)
542	endif
543	status = nf90_inquire_dimension(fid, did , len = ns)
544	if (status /= NF90_NOERR) then
545	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
546	call oasis_abort(file=__FILE__,line=__LINE__)
547	endif
548	! status = nf90_inq_dimid (fid, 'n_a', did) ! size of input vector
549	status = nf90_inq_dimid (fid, 'src_grid_size', did) ! size of input vector
550	if (status /= NF90_NOERR) then
551	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
552	write(nulprt,*) subname,estr,'dim not found = ','src_grid_size'
553	call oasis_abort(file=__FILE__,line=__LINE__)
554	endif
555	status = nf90_inquire_dimension(fid, did , len = na)
556	if (status /= NF90_NOERR) then
557	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
558	call oasis_abort(file=__FILE__,line=__LINE__)
559	endif
560	! status = nf90_inq_dimid (fid, 'n_b', did) ! size of output vector
561	status = nf90_inq_dimid (fid, 'dst_grid_size', did) ! size of output vector
562	if (status /= NF90_NOERR) then
563	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
564	write(nulprt,*) subname,estr,'dim not found = ','dst_grid_size'
565	call oasis_abort(file=__FILE__,line=__LINE__)
566	endif
567	status = nf90_inquire_dimension(fid, did , len = nb)
568	if (status /= NF90_NOERR) then
569	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
570	call oasis_abort(file=__FILE__,line=__LINE__)
571	endif
572	status = nf90_inq_dimid (fid, 'num_wgts', did) ! size of output vector
573	if (status /= NF90_NOERR) then
574	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
575	write(nulprt,*) subname,estr,'dim not found = ','num_wgts'
576	call oasis_abort(file=__FILE__,line=__LINE__)
577	endif
578	status = nf90_inquire_dimension(fid, did , len = nwgts)
579	if (status /= NF90_NOERR) then
580	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
581	call oasis_abort(file=__FILE__,line=__LINE__)
582	endif
583
584	if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
585	! status = nf90_inq_dimid (fid, 'ni_a', did) ! number of lons in input grid
586	! status = nf90_inquire_dimension(fid, did , len = ni_i)
587	! status = nf90_inq_dimid (fid, 'nj_a', did) ! number of lats in input grid
588	! status = nf90_inquire_dimension(fid, did , len = nj_i)
589	! status = nf90_inq_dimid (fid, 'ni_b', did) ! number of lons in output grid
590	! status = nf90_inquire_dimension(fid, did , len = ni_o)
591	! status = nf90_inq_dimid (fid, 'nj_b', did) ! number of lats in output grid
592	! status = nf90_inquire_dimension(fid, did , len = nj_o)
593	status = nf90_inq_varid(fid, 'src_grid_dims', vid)
594	if (status /= NF90_NOERR) then
595	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
596	write(nulprt,*) subname,estr,'var not found = ','src_grid_dims'
597	call oasis_abort(file=__FILE__,line=__LINE__)
598	endif
599	status = nf90_get_var(fid, vid, dims)
600	if (status /= NF90_NOERR) then
601	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
602	call oasis_abort(file=__FILE__,line=__LINE__)
603	endif
604	ni_i = dims(1)
605	nj_i = dims(2)
606	status = nf90_inq_varid(fid, 'dst_grid_dims', vid)
607	if (status /= NF90_NOERR) then
608	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
609	write(nulprt,*) subname,estr,'var not found = ','dst_grid_dims'
610	call oasis_abort(file=__FILE__,line=__LINE__)
611	endif
612	status = nf90_get_var(fid, vid, dims)
613	if (status /= NF90_NOERR) then
614	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
615	call oasis_abort(file=__FILE__,line=__LINE__)
616	endif
617	ni_o = dims(1)
618	nj_o = dims(2)
619	end if
620
621	if (OASIS_debug >= 2) write(nulprt,) subname," matrix dims src x dst : ",na,' x',nb
622	if (OASIS_debug >= 2) write(nulprt,) subname," number of non-zero elements: ",ns
623
624	endif
625
626	!----------------------------------------------------------------------------
627	!> * Read and load area_a on root task
628	!----------------------------------------------------------------------------
629
630	if (present(areasrc)) then
631	if (mytask == 0) then
632	call mct_aVect_init(areasrc0,' ',areaAV_field,na)
633	! status = nf90_inq_varid (fid,'area_a',vid)
634	status = nf90_inq_varid (fid,'src_grid_area',vid)
635	if (status /= NF90_NOERR) THEN
636	write(nulprt,*) subname,' nf90_strerrr = ',TRIM(nf90_strerror(status))
637	write(nulprt,*) subname,estr,'var not found = ','src_grid_area'
638	call oasis_flush(nulprt)
639	endif
640	status = nf90_get_var(fid, vid, areasrc0%rAttr)
641	if (status /= NF90_NOERR) THEN
642	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
643	call oasis_flush(nulprt)
644	endif
645	endif
646	call mct_aVect_scatter(areasrc0, areasrc, SgsMap, 0, mpicom, status)
647	if (status /= 0) call mct_die(subname,"Error on scatter of areasrc0")
648	if (mytask == 0) then
649	! if (present(dbug)) then
650	! if (dbug > 2) then
651	! write(nulprt,*) subName,'Size of src ',mct_aVect_lSize(areasrc0)
652	! write(nulprt,*) subName,'min/max src ',minval(areasrc0%rAttr(1,:)),maxval(areasrc0%rAttr(1,:))
653	! endif
654	! end if
655	call mct_aVect_clean(areasrc0)
656	end if
657	end if
658
659	!----------------------------------------------------------------------------
660	!> * Read and load area_b on root task
661	!----------------------------------------------------------------------------
662
663	if (present(areadst)) then
664	if (mytask == 0) then
665	call mct_aVect_init(areadst0,' ',areaAV_field,nb)
666	! status = nf90_inq_varid (fid,'area_b',vid)
667	status = nf90_inq_varid (fid,'dst_grid_area',vid)
668	if (status /= NF90_NOERR) THEN
669	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
670	write(nulprt,*) subname,estr,'var not found = ','dst_grid_area'
671	call oasis_flush(nulprt)
672	endif
673	status = nf90_get_var(fid, vid, areadst0%rAttr)
674	if (status /= NF90_NOERR) THEN
675	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
676	call oasis_flush(nulprt)
677	endif
678	endif
679	call mct_aVect_scatter(areadst0, areadst, DgsMap, 0, mpicom, status)
680	if (status /= 0) call mct_die(subname,"Error on scatter of areadst0")
681	if (mytask == 0) then
682	! if (present(dbug)) then
683	! if (dbug > 2) then
684	! write(nulprt,*) subName,'Size of dst ',mct_aVect_lSize(areadst0)
685	! write(nulprt,*) subName,'min/max dst ',minval(areadst0%rAttr(1,:)),maxval(areadst0%rAttr(1,:))
686	! endif
687	! end if
688	call mct_aVect_clean(areadst0)
689	endif
690	endif
691
692	!----------------------------------------------------------------------------
693	!> * Broadcast ni and nj if requested
694	!----------------------------------------------------------------------------
695
696	if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
697	call oasis_mpi_bcast(ni_i,mpicom,subName//" MPI in ni_i bcast")
698	call oasis_mpi_bcast(nj_i,mpicom,subName//" MPI in nj_i bcast")
699	call oasis_mpi_bcast(ni_o,mpicom,subName//" MPI in ni_o bcast")
700	call oasis_mpi_bcast(nj_o,mpicom,subName//" MPI in nj_o bcast")
701	end if
702
703	!----------------------------------------------------------------------------
704	!> * Broadcast array sizes and allocate arrays for local storage
705	!----------------------------------------------------------------------------
706
707	call oasis_mpi_bcast(ns,mpicom,subName//" MPI in ns bcast")
708	call oasis_mpi_bcast(na,mpicom,subName//" MPI in na bcast")
709	call oasis_mpi_bcast(nb,mpicom,subName//" MPI in nb bcast")
710	call oasis_mpi_bcast(nwgts,mpicom,subName//" MPI in nwgts bcast")
711	if (local_timers_on) call oasis_timer_stop('map_read_orig_read')
712
713	if (local_timers_on) call oasis_timer_start('map_read_orig_sort')
714	!--- setup local seg map, sorted
715	if (newdom == 'src') then
716	mygsmap => DgsMap
717	elseif (newdom == 'dst') then
718	mygsmap => SgsMap
719	else
720	write(nulprt,*) subname,estr,'invalid newdom value, expect src or dst = ',newdom
721	call oasis_abort(file=__FILE__,line=__LINE__)
722	endif
723	lsize = 0
724	do n = 1,size(mygsmap%start)
725	if (mygsmap%pe_loc(n) == mytask) then
726	lsize=lsize+1
727	endif
728	enddo
729	allocate(lsstart(lsize),lscount(lsize),stat=status)
730	if (status /= 0) call mct_perr_die(subName,':: allocate Lsstart',status)
731
732	!----------------------------------------------------------------------------
733	!> * Initialize lsstart and lscount, the sorted list of local indices
734	!----------------------------------------------------------------------------
735
736	lsize = 0
737	do n = 1,size(mygsmap%start)
738	if (mygsmap%pe_loc(n) == mytask) then ! on my pe
739	lsize=lsize+1
740	found = .false.
741	l1 = 1
742	do while (.not.found .and. l1 < lsize) ! bubble sort copy
743	if (mygsmap%start(n) < lsstart(l1)) then
744	do l2 = lsize, l1+1, -1
745	lsstart(l2) = lsstart(l2-1)
746	lscount(l2) = lscount(l2-1)
747	enddo
748	found = .true.
749	else
750	l1 = l1 + 1
751	endif
752	enddo
753	lsstart(l1) = mygsmap%start(n)
754	lscount(l1) = mygsmap%length(n)
755	endif
756	enddo
757	do n = 1,lsize-1
758	if (lsstart(n) > lsstart(n+1)) then
759	write(nulprt,*) subname,estr,'lsstart not properly sorted'
760	call oasis_abort(file=__FILE__,line=__LINE__)
761	endif
762	enddo
763	if (local_timers_on) call oasis_timer_stop('map_read_orig_sort')
764
765	!----------------------------------------------------------------------------
766	!> * Compute the number of chunks to read, read size is rbuf_size
767	!----------------------------------------------------------------------------
768
769	if (local_timers_on) call oasis_timer_start('map_read_orig_dist')
770	rsize = min(rbuf_size,ns) ! size of i/o chunks
771	bsize = ((ns/commsize) + 1 ) * 1.2 ! local temporary buffer size
772	if (ns == 0) then
773	nread = 0
774	else
775	nread = (ns-1)/rsize + 1 ! num of reads to do
776	endif
777
778	!----------------------------------------------------------------------------
779	!> * Allocate arrays for local weights plus row and column indices
780	!----------------------------------------------------------------------------
781
782	if (mytask == 0) then
783	allocate(remaps(nwgts,rsize),stat=status)
784	if (status /= 0) call mct_perr_die(subName,':: allocate remaps',status)
785	endif
786	allocate(Smat(nwgts),stat=status)
787	if (status /= 0) call mct_perr_die(subName,':: allocate Smat',status)
788	allocate(Sbuf(nwgts,rsize),Rbuf(rsize),Cbuf(rsize),stat=status)
789	if (status /= 0) call mct_perr_die(subName,':: allocate Sbuf',status)
790	allocate(Snew(nwgts,bsize),Cnew(bsize),Rnew(bsize),stat=status)
791	if (status /= 0) call mct_perr_die(subName,':: allocate Snew1',status)
792
793	!----------------------------------------------------------------------------
794	!> * Loop over the chunks of weights data
795	!----------------------------------------------------------------------------
796
797	cnt = 0
798	do n = 1,nread
799	start(1) = (n-1)*rsize + 1
800	count(1) = min(rsize,ns-start(1)+1)
801	start2(1) = 1
802	count2(1) = nwgts
803	start2(2) = start(1)
804	count2(2) = count(1)
805
806	!----------------------------------------------------------------------------
807	!> * Read chunk of data on root pe
808	!----------------------------------------------------------------------------
809
810	if (mytask== 0) then
811	! status = nf90_inq_varid (fid,'S' ,vid)
812	status = nf90_inq_varid (fid,'remap_matrix' ,vid)
813	if (status /= NF90_NOERR) then
814	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
815	write(nulprt,*) subname,estr,'var not found = ','remap_matrix'
816	call oasis_abort(file=__FILE__,line=__LINE__)
817	endif
818	! status = nf90_get_var(fid,vid,start,count,Sbuf)
819	status = nf90_get_var(fid,vid,remaps,start2,count2)
820	if (status /= NF90_NOERR) THEN
821	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
822	call oasis_flush(nulprt)
823	endif
824	Sbuf(:,:) = remaps(:,:)
825
826	! status = nf90_inq_varid (fid,'row',vid)
827	status = nf90_inq_varid (fid,'dst_address',vid)
828	if (status /= NF90_NOERR) then
829	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
830	write(nulprt,*) subname,estr,'var not found = ','dst_address'
831	call oasis_abort(file=__FILE__,line=__LINE__)
832	endif
833	status = nf90_get_var (fid,vid,Rbuf,start,count)
834	if (status /= NF90_NOERR) THEN
835	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
836	call oasis_flush(nulprt)
837	endif
838
839	! status = nf90_inq_varid (fid,'col',vid)
840	status = nf90_inq_varid (fid,'src_address',vid)
841	if (status /= NF90_NOERR) then
842	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
843	write(nulprt,*) subname,estr,'var not found = ','src_address'
844	call oasis_abort(file=__FILE__,line=__LINE__)
845	endif
846	status = nf90_get_var (fid,vid,Cbuf,start,count)
847	if (status /= NF90_NOERR) THEN
848	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
849	call oasis_flush(nulprt)
850	endif
851	endif
852
853	!----------------------------------------------------------------------------
854	!> * Broadcast S, row, col to all tasks
855	!----------------------------------------------------------------------------
856
857	call oasis_mpi_bcast(Sbuf,mpicom,subName//" MPI in Sbuf bcast")
858	call oasis_mpi_bcast(Rbuf,mpicom,subName//" MPI in Rbuf bcast")
859	call oasis_mpi_bcast(Cbuf,mpicom,subName//" MPI in Cbuf bcast")
860
861	!----------------------------------------------------------------------------
862	!> * Each task keeps only the data required
863	!----------------------------------------------------------------------------
864
865	if (namwgtopt == "abort_on_bad_index") then
866	abort_weight = .false.
867	do m = 1,count(1)
868	!--- check for bad weights
869	if ((Rbuf(m) <= 0 .or. Rbuf(m) > nb .or. &
870	Cbuf(m) <= 0 .or. Cbuf(m) > na) &
871	! tcx weight = 0
872	.and. (minval(Sbuf(:,m)) /= 0._R8 .or. maxval(Sbuf(:,m)) /= 0._R8) &
873	) then
874	abort_weight = .true.
875	write(nulprt,'(3A,I12,A,I12,A,I12,A,G14.7,A,G14.7,A)') &
876	subname,wstr,'BAD weight found in '//trim(filename), &
877	m,'=id',Cbuf(m),'=src',Rbuf(m),'=dst',minval(Sbuf(:,m)),'=minS',maxval(Sbuf(:,m)),'=maxS'
878	endif
879	enddo
880	if (abort_weight) then
881	write(nulprt,*) subname,wstr,'BAD weight found, aborting'
882	call oasis_abort(file=__FILE__,line=__LINE__)
883	endif
884	endif
885
886	do m = 1,count(1)
887	!--- check for bad weights
888	if ((namwgtopt(1:16) == "ignore_bad_index") .and. &
889	(Rbuf(m) <= 0 .or. Rbuf(m) > nb .or. &
890	Cbuf(m) <= 0 .or. Cbuf(m) > na)) then
891	mywt = .false.
892	! tcx weight = 0
893	if (minval(Sbuf(:,m)) /= 0._R8 .or. maxval(Sbuf(:,m)) /= 0._R8) then
894	if (OASIS_debug >= 2 .and. namwgtopt /= "ignore_bad_index_silently") then
895	write(nulprt,'(3A,I12,A,I12,A,I12,A,G14.7,A,G14.7,A)') &
896	subname,wstr,'BAD weight found in '//trim(filename), &
897	m,'=id',Cbuf(m),'=src',Rbuf(m),'=dst',minval(Sbuf(:,m)),'=minS',maxval(Sbuf(:,m)),'=maxS'
898	endif
899	endif
900	elseif (newdom == 'src') then
901	mywt = check_myindex(Rbuf(m),lsstart,lscount)
902	elseif (newdom == 'dst') then
903	mywt = check_myindex(Cbuf(m),lsstart,lscount)
904	endif
905
906	if (mywt) then
907	cntold = cnt
908	cnt = cnt + 1
909
910	!----------------------------------------------------------------------------
911	!> * Reallocate local weights arrays if they need to be bigger
912	!----------------------------------------------------------------------------
913	if (cnt > bsize) then
914	!--- allocate old arrays and copy new into old
915	allocate(Sold(1:nwgts,cntold),Rold(cntold),Cold(cntold),stat=status)
916	if (status /= 0) call mct_perr_die(subName,':: allocate old',status)
917	Sold(1:nwgts,1:cntold) = Snew(1:nwgts,1:cntold)
918	Rold(1:cntold) = Rnew(1:cntold)
919	Cold(1:cntold) = Cnew(1:cntold)
920
921	!--- reallocate new to bigger size, increase buffer by 50% (arbitrary)
922	deallocate(Snew,Rnew,Cnew,stat=status)
923	if (status /= 0) call mct_perr_die(subName,':: allocate new',status)
924	bsize = 1.5 * bsize
925	if (OASIS_debug > 15) write(nulprt,*) subname,' reallocate bsize to ',bsize
926	allocate(Snew(nwgts,bsize),Rnew(bsize),Cnew(bsize),stat=status)
927	if (status /= 0) call mct_perr_die(subName,':: allocate old',status)
928
929	!--- copy data back into new
930	Snew(1:nwgts,1:cntold) = Sold(1:nwgts,1:cntold)
931	Rnew(1:cntold) = Rold(1:cntold)
932	Cnew(1:cntold) = Cold(1:cntold)
933	deallocate(Sold,Rold,Cold,stat=status)
934	if (status /= 0) call mct_perr_die(subName,':: deallocate old',status)
935	endif
936
937	Snew(1:nwgts,cnt) = Sbuf(1:nwgts,m)
938	Rnew(cnt) = Rbuf(m)
939	Cnew(cnt) = Cbuf(m)
940	endif
941	enddo ! count
942	enddo ! nread
943
944	!----------------------------------------------------------------------------
945	!> * Clean up arrays
946	!----------------------------------------------------------------------------
947
948	if (mytask == 0) then
949	deallocate(remaps, stat=status)
950	if (status /= 0) call mct_perr_die(subName,':: deallocate remaps',status)
951	endif
952	deallocate(Sbuf,Rbuf,Cbuf, stat=status)
953	if (status /= 0) call mct_perr_die(subName,':: deallocate Sbuf',status)
954	if (local_timers_on) call oasis_timer_stop('map_read_orig_dist')
955
956	!----------------------------------------------------------------------------
957	!> * Initialize the mct sMat data type
958	!----------------------------------------------------------------------------
959
960	if (local_timers_on) call oasis_timer_start('map_read_orig_smati')
961	! mct_sMat_init must be given the number of rows and columns that
962	! would be in the full matrix. Nrows= size of output vector=nb.
963	! Ncols = size of input vector = na.
964	do n = 1,nwgts
965	call mct_sMat_init(sMat(n), nb, na, cnt)
966	enddo
967
968	igrow = mct_sMat_indexIA(sMat(1),'grow')
969	igcol = mct_sMat_indexIA(sMat(1),'gcol')
970	iwgt = mct_sMat_indexRA(sMat(1),'weight')
971
972	if (cnt /= 0) then
973	do n = 1,nwgts
974	sMat(n)%data%rAttr(iwgt ,1:cnt) = Snew(n,1:cnt)
975	sMat(n)%data%iAttr(igrow,1:cnt) = Rnew(1:cnt)
976	sMat(n)%data%iAttr(igcol,1:cnt) = Cnew(1:cnt)
977	enddo
978	endif
979	if (local_timers_on) call oasis_timer_stop('map_read_orig_smati')
980
981	!----------------------------------------------------------------------------
982	!> * More clean up
983	!----------------------------------------------------------------------------
984
985	if (local_timers_on) call oasis_timer_start('map_read_orig_clean')
986	deallocate(Snew,Rnew,Cnew, stat=status)
987	deallocate(lsstart,lscount,stat=status)
988	if (status /= 0) call mct_perr_die(subName,':: deallocate new',status)
989
990	if (mytask == 0) then
991	status = nf90_close(fid)
992	if (status /= NF90_NOERR) then
993	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
994	call oasis_abort(file=__FILE__,line=__LINE__)
995	endif
996	if (OASIS_debug >= 2) THEN
997	write(nulprt,*) subname," ... done reading file"
998	call oasis_flush(nulprt)
999	endif
1000	endif
1001	if (local_timers_on) call oasis_timer_stop('map_read_orig_clean')
1002
1003	call oasis_debug_exit(subname)
1004
1005	end subroutine oasis_map_sMatReaddnc_orig
1006
1007	!===============================================================================
1008	!BOP ===========================================================================
1009	!> Read in mapping matrix data from a SCRIP netCDF file using smart scatter (ceg)
1010	!
1011	! !IROUTINE: oasis_map_sMatReaddnc_ceg - Do a distributed read of a NetCDF SCRIP file and
1012	! return weights in a distributed SparseMatrix
1013	!
1014	! !DESCRIPTION:
1015	!> Read in mapping matrix data from a SCRIP netCDF data file using
1016	!> a low memory method and then scatter to all pes using a smart method
1017	!> where only select data is sent to tasks. Based on the sMatReaddnc method
1018	!> from CESM1.0.3
1019	! !REMARKS:
1020	!> This routine leverages gsmaps to determine scatter pattern
1021	!>
1022	!> The scatter is implemented via the root task reading the data and
1023	!> then determining which task gets which weights from the gsmap.
1024	!> The root the sends specific data to each task.
1025	!>
1026	!> The algorithm to determine which task a weigth belongs to involves
1027	!> checking the task ownership for a given global index.
1028	!>
1029	!> The local buffer sizes are estimated up front based on ngridcell/npes
1030	!> plus 20% (see 1.2 below). If the local buffer size fills up, then
1031	!> the buffer is reallocated 50% larger (see 1.5 below) and the fill
1032	!> continues. The idea is to trade off memory reallocation and copy
1033	!> with memory usage. 1.2 and 1.5 are arbitary, other values may
1034	!> result in better performance.
1035	!>
1036	!> Once all the matrix weights have been read, the sMat is initialized,
1037	!> the values from the buffers are copied in, and everything is deallocated.
1038	!
1039	! !INTERFACE: -----------------------------------------------------------------
1040
1041	subroutine oasis_map_sMatReaddnc_ceg(sMat,SgsMap,DgsMap,newdom, &
1042	fileName,mytask,mpicom,nwgts, &
1043	areasrc,areadst,ni_i,nj_i,ni_o,nj_o )
1044
1045	! !USES:
1046
1047	! !INPUT/OUTPUT PARAMETERS:
1048
1049	type(mct_sMat) ,intent(out),pointer :: sMat(:) !< mapping data
1050	type(mct_gsMap) ,intent(in) ,target :: SgsMap !< src gsmap
1051	type(mct_gSMap) ,intent(in) ,target :: DgsMap !< dst gsmap
1052	character(*) ,intent(in) :: newdom !< type of sMat (src or dst)
1053	!< src = rearrange and map (bfb), dst = map and rearrange (partial sums)
1054	character(*) ,intent(in) :: filename!< netCDF file to read
1055	integer(IN) ,intent(in) :: mytask !< processor id
1056	integer(IN) ,intent(in) :: mpicom !< mpi communicator
1057	integer(IN) ,intent(out) :: nwgts !< number of weights
1058	type(mct_Avect) ,intent(out), optional :: areasrc !< area of src grid from mapping file
1059	type(mct_Avect) ,intent(out), optional :: areadst !< area of dst grid from mapping file
1060	integer(IN) ,intent(out), optional :: ni_i !< number of lons on input grid
1061	integer(IN) ,intent(out), optional :: nj_i !< number of lats on input grid
1062	integer(IN) ,intent(out), optional :: ni_o !< number of lons on output grid
1063	integer(IN) ,intent(out), optional :: nj_o !< number of lats on output grid
1064
1065	! !EOP
1066
1067	!--- local ---
1068	integer :: n,m ! generic loop indicies
1069	integer :: na ! size of source domain
1070	integer :: nb ! size of destination domain
1071	integer :: ns ! number of non-zero elements in matrix
1072	integer :: ni,nj ! number of row and col in the matrix
1073	integer :: igrow ! aVect index for matrix row
1074	integer :: igcol ! aVect index for matrix column
1075	integer :: iwgt ! aVect index for matrix element
1076	integer :: iarea ! aVect index for area
1077	integer :: rsize ! size of read buffer
1078	integer :: cnt ! local num of wgts
1079	integer :: start(1)! netcdf read
1080	integer :: count(1)! netcdf read
1081	integer :: start2(2)! netcdf read
1082	integer :: count2(2)! netcdf read
1083	integer :: bsize ! buffer size
1084	integer :: nread ! number of reads
1085	logical :: mywt ! does this weight belong on my pe
1086	integer :: dims(2)
1087	logical :: abort_weight ! flag if bad weight found
1088
1089	!--- buffers for i/o ---
1090	real(R8) ,allocatable :: rtemp(:) ! real temporary
1091	real(R8) ,allocatable :: Sbuf(:,:) ! real weights
1092	real(R8) ,allocatable :: remaps(:,:) ! real weights with num_wgts dim
1093	integer,allocatable :: Rbuf(:) ! ints rows
1094	integer,allocatable :: Cbuf(:) ! ints cols
1095	real(R8),allocatable :: SReadData(:,:) !
1096	integer,allocatable :: RReadData(:) !
1097	integer,allocatable :: CReadData(:) !
1098	integer,allocatable :: pesave(:) !
1099	real(R8),allocatable :: SDistData(:,:) !
1100	integer,allocatable :: RDistData(:) !
1101	integer,allocatable :: CDistData(:) !
1102
1103	!--- variables associated with local computation of global indices
1104	integer :: commsize ! size of local communicator
1105	type(mct_gsMap),pointer :: mygsmap ! pointer to one of the gsmaps
1106	integer :: l1,l2,lsize ! generice indices for sort
1107	integer,allocatable :: lsstart(:) ! local seg map info
1108	integer,allocatable :: lscount(:) ! local seg map info
1109	integer,allocatable :: lspeloc(:) ! local seg map info
1110	integer,allocatable :: sortkey(:) ! local sort key info
1111	logical :: found ! for sort
1112	integer :: pe ! Process ID of owning process
1113	integer :: reclen ! Length of Rbuf/Cbuf to be received
1114	integer :: ierr ! MPI error check
1115	integer,allocatable :: cntrs(:) ! Counters of number of elements stored for each processor
1116	integer :: mpistatus(MPI_STATUS_SIZE) ! MPI_Status object
1117
1118	character,allocatable :: str(:) ! variable length char string
1119	character(len=ic_long):: attstr ! netCDF attribute name string
1120	integer :: status ! netCDF routine return code
1121	integer :: fid ! netCDF file ID
1122	integer :: vid ! netCDF variable ID
1123	integer :: did ! netCDF dimension ID
1124	!--- arbitrary size of read buffer, this is the chunk size weights reading
1125	integer,parameter :: rbuf_size = 100000
1126	integer :: dimbuffer(8)
1127
1128	!--- global source and destination areas ---
1129	type(mct_Avect) :: areasrc0 ! area of src grid from mapping file
1130	type(mct_Avect) :: areadst0 ! area of src grid from mapping file
1131
1132	character(*),parameter :: areaAV_field = 'aream'
1133
1134	!--- formats ---
1135	character(*),parameter :: subName = '(oasis_map_sMatReaddnc_ceg)'
1136	character*80 :: fname
1137
1138	!-------------------------------------------------------------------------------
1139	!
1140	!-------------------------------------------------------------------------------
1141
1142	call oasis_debug_enter(subname)
1143	call oasis_mpi_commsize(mpicom,commsize)
1144	nwgts = -1
1145	if (local_timers_on) call oasis_timer_start('map_read_ceg_read')
1146
1147	if (mytask == 0) then
1148	if (OASIS_debug >= 2) write(nulprt,*) subname," reading mapping matrix data decomposed..."
1149
1150	!----------------------------------------------------------------------------
1151	!> * Open and read the file SCRIP weights size on the root task
1152	!----------------------------------------------------------------------------
1153
1154	if (OASIS_debug >=2 ) write(nulprt,) subname," file name : ",trim(fileName)
1155	status = nf90_open(trim(filename),NF90_NOWRITE,fid)
1156	if (status /= NF90_NOERR) then
1157	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1158	write(nulprt,*) subname,estr,'mapping file not found = ',trim(filename)
1159	call oasis_abort(file=__FILE__,line=__LINE__)
1160	endif
1161
1162	!--- get matrix dimensions ----------
1163	! status = nf90_inq_dimid (fid, 'n_s', did) ! size of sparse matrix
1164	status = nf90_inq_dimid (fid, 'num_links', did) ! size of sparse matrix
1165	if (status /= NF90_NOERR) then
1166	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1167	write(nulprt,*) subname,estr,'dim not found = ','num_links'
1168	call oasis_abort(file=__FILE__,line=__LINE__)
1169	endif
1170	status = nf90_inquire_dimension(fid, did , len = ns)
1171	if (status /= NF90_NOERR) then
1172	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1173	call oasis_abort(file=__FILE__,line=__LINE__)
1174	endif
1175	! status = nf90_inq_dimid (fid, 'n_a', did) ! size of input vector
1176	status = nf90_inq_dimid (fid, 'src_grid_size', did) ! size of input vector
1177	if (status /= NF90_NOERR) then
1178	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1179	write(nulprt,*) subname,estr,'dim not found = ','src_grid_size'
1180	call oasis_abort(file=__FILE__,line=__LINE__)
1181	endif
1182	status = nf90_inquire_dimension(fid, did , len = na)
1183	if (status /= NF90_NOERR) then
1184	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1185	call oasis_abort(file=__FILE__,line=__LINE__)
1186	endif
1187	! status = nf90_inq_dimid (fid, 'n_b', did) ! size of output vector
1188	status = nf90_inq_dimid (fid, 'dst_grid_size', did) ! size of output vector
1189	if (status /= NF90_NOERR) then
1190	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1191	write(nulprt,*) subname,estr,'dim not found = ','dst_grid_size'
1192	call oasis_abort(file=__FILE__,line=__LINE__)
1193	endif
1194	status = nf90_inquire_dimension(fid, did , len = nb)
1195	if (status /= NF90_NOERR) then
1196	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1197	call oasis_abort(file=__FILE__,line=__LINE__)
1198	endif
1199	status = nf90_inq_dimid (fid, 'num_wgts', did) ! size of output vector
1200	if (status /= NF90_NOERR) then
1201	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1202	write(nulprt,*) subname,estr,'dim not found = ','num_wgts'
1203	call oasis_abort(file=__FILE__,line=__LINE__)
1204	endif
1205	status = nf90_inquire_dimension(fid, did , len = nwgts)
1206	if (status /= NF90_NOERR) then
1207	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1208	call oasis_abort(file=__FILE__,line=__LINE__)
1209	endif
1210
1211	if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
1212	! status = nf90_inq_dimid (fid, 'ni_a', did) ! number of lons in input grid
1213	! status = nf90_inquire_dimension(fid, did , len = ni_i)
1214	! status = nf90_inq_dimid (fid, 'nj_a', did) ! number of lats in input grid
1215	! status = nf90_inquire_dimension(fid, did , len = nj_i)
1216	! status = nf90_inq_dimid (fid, 'ni_b', did) ! number of lons in output grid
1217	! status = nf90_inquire_dimension(fid, did , len = ni_o)
1218	! status = nf90_inq_dimid (fid, 'nj_b', did) ! number of lats in output grid
1219	! status = nf90_inquire_dimension(fid, did , len = nj_o)
1220	status = nf90_inq_varid(fid, 'src_grid_dims', vid)
1221	if (status /= NF90_NOERR) then
1222	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1223	write(nulprt,*) subname,estr,'var not found = ','src_grid_dims'
1224	call oasis_abort(file=__FILE__,line=__LINE__)
1225	endif
1226	status = nf90_get_var(fid, vid, dims)
1227	if (status /= NF90_NOERR) then
1228	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1229	call oasis_abort(file=__FILE__,line=__LINE__)
1230	endif
1231	ni_i = dims(1)
1232	nj_i = dims(2)
1233	status = nf90_inq_varid(fid, 'dst_grid_dims', vid)
1234	if (status /= NF90_NOERR) then
1235	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1236	write(nulprt,*) subname,estr,'var not found = ','dst_grid_dims'
1237	call oasis_abort(file=__FILE__,line=__LINE__)
1238	endif
1239	status = nf90_get_var(fid, vid, dims)
1240	if (status /= NF90_NOERR) then
1241	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1242	call oasis_abort(file=__FILE__,line=__LINE__)
1243	endif
1244	ni_o = dims(1)
1245	nj_o = dims(2)
1246	end if
1247
1248	if (OASIS_debug >= 2) write(nulprt,) subname," matrix dims src x dst : ",na,' x',nb
1249	if (OASIS_debug >= 2) write(nulprt,) subname," number of non-zero elements: ",ns
1250
1251	endif
1252
1253	!----------------------------------------------------------------------------
1254	!> * Read and load area_a on root task
1255	!----------------------------------------------------------------------------
1256
1257	if (present(areasrc)) then
1258	if (mytask == 0) then
1259	call mct_aVect_init(areasrc0,' ',areaAV_field,na)
1260	! status = nf90_inq_varid (fid,'area_a',vid)
1261	status = nf90_inq_varid (fid,'src_grid_area',vid)
1262	if (status /= NF90_NOERR) THEN
1263	write(nulprt,*) subname,' nf90_strerrr = ',TRIM(nf90_strerror(status))
1264	write(nulprt,*) subname,estr,'var not found = ','src_grid_area'
1265	call oasis_flush(nulprt)
1266	endif
1267	status = nf90_get_var(fid, vid, areasrc0%rAttr)
1268	if (status /= NF90_NOERR) THEN
1269	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
1270	call oasis_flush(nulprt)
1271	endif
1272	endif
1273	call mct_aVect_scatter(areasrc0, areasrc, SgsMap, 0, mpicom, status)
1274	if (status /= 0) call mct_die(subname,"Error on scatter of areasrc0")
1275	if (mytask == 0) then
1276	! if (present(dbug)) then
1277	! if (dbug > 2) then
1278	! write(nulprt,*) subName,'min/max src ',minval(areasrc0%rAttr(1,:)),maxval(areasrc0%rAttr(1,:))
1279	! endif
1280	! end if
1281	call mct_aVect_clean(areasrc0)
1282	end if
1283	end if
1284
1285	!----------------------------------------------------------------------------
1286	!> * Read and load area_b on root task
1287	!----------------------------------------------------------------------------
1288
1289	if (present(areadst)) then
1290	if (mytask == 0) then
1291	call mct_aVect_init(areadst0,' ',areaAV_field,nb)
1292	! status = nf90_inq_varid (fid,'area_b',vid)
1293	status = nf90_inq_varid (fid,'dst_grid_area',vid)
1294	if (status /= NF90_NOERR) THEN
1295	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
1296	write(nulprt,*) subname,estr,'var not found = ','dst_grid_area'
1297	call oasis_flush(nulprt)
1298	endif
1299	status = nf90_get_var(fid, vid, areadst0%rAttr)
1300	if (status /= NF90_NOERR) THEN
1301	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
1302	call oasis_flush(nulprt)
1303	endif
1304	endif
1305	call mct_aVect_scatter(areadst0, areadst, DgsMap, 0, mpicom, status)
1306	if (status /= 0) call mct_die(subname,"Error on scatter of areadst0")
1307	if (mytask == 0) then
1308	! if (present(dbug)) then
1309	! if (dbug > 2) then
1310	! write(nulprt,*) subName,'min/max dst ',minval(areadst0%rAttr(1,:)),maxval(areadst0%rAttr(1,:))
1311	! endif
1312	! end if
1313	call mct_aVect_clean(areadst0)
1314	endif
1315	endif
1316
1317	!----------------------------------------------------------------------------
1318	!> * Broadcast ni and nj if requested
1319	!> * Broadcast array sizes and allocate arrays for local storage
1320	! Replace 8 MPI_Bcasts with just one
1321	!----------------------------------------------------------------------------
1322
1323	if (mpi_rank_local.eq.0) then
1324	dimbuffer(1) = ns
1325	dimbuffer(2) = na
1326	dimbuffer(3) = nb
1327	dimbuffer(4) = nwgts
1328	if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
1329	dimbuffer(5) = ni_i
1330	dimbuffer(6) = nj_i
1331	dimbuffer(7) = ni_o
1332	dimbuffer(8) = nj_o
1333	end if
1334	endif
1335	call oasis_mpi_bcast(dimbuffer,mpicom,subName//" MPI of dimbuffer")
1336	if (mpi_rank_local.ne.0) then
1337	ns = dimbuffer(1)
1338	na = dimbuffer(2)
1339	nb = dimbuffer(3)
1340	nwgts = dimbuffer(4)
1341	if (present(ni_i) .and. present(nj_i) .and. present(ni_o) .and. present(nj_o)) then
1342	ni_i = dimbuffer(5)
1343	nj_i = dimbuffer(6)
1344	ni_o = dimbuffer(7)
1345	nj_o = dimbuffer(8)
1346	end if
1347	endif
1348
1349	! call oasis_mpi_bcast(ni_i,mpicom,subName//" MPI in ni_i bcast")
1350	! call oasis_mpi_bcast(nj_i,mpicom,subName//" MPI in nj_i bcast")
1351	! call oasis_mpi_bcast(ni_o,mpicom,subName//" MPI in ni_o bcast")
1352	! call oasis_mpi_bcast(nj_o,mpicom,subName//" MPI in nj_o bcast")
1353
1354	!--- setup local seg map, sorted
1355	if (newdom == 'src') then
1356	mygsmap => DgsMap
1357	elseif (newdom == 'dst') then
1358	mygsmap => SgsMap
1359	else
1360	write(nulprt,*) subname,estr,'invalid newdom value, expect src or dst = ',newdom
1361	call oasis_abort(file=__FILE__,line=__LINE__)
1362	endif
1363
1364	!----------------------------------------------------------------------------
1365	!> * Compute the number of chunks to read, read size is rbuf_size
1366	!----------------------------------------------------------------------------
1367
1368	rsize = min(rbuf_size,ns) ! size of i/o chunks
1369	bsize = ((ns/commsize) + 1 ) * 1.2 ! local temporary buffer size
1370	if (ns == 0) then
1371	nread = 0
1372	else
1373	nread = (ns-1)/rsize + 1 ! num of reads to do
1374	endif
1375
1376	!----------------------------------------------------------------------------
1377	!> * Allocate arrays for local weights plus row and column indices
1378	!----------------------------------------------------------------------------
1379
1380	allocate(Smat(nwgts),stat=status)
1381	if (status /= 0) call mct_perr_die(subName,':: allocate Smat',status)
1382	allocate(Sbuf(nwgts,rsize),Rbuf(rsize),Cbuf(rsize),stat=status)
1383	if (status /= 0) call mct_perr_die(subName,':: allocate Sbuf',status)
1384	allocate(Snew(nwgts,bsize),Cnew(bsize),Rnew(bsize),stat=status)
1385	if (status /= 0) call mct_perr_die(subName,':: allocate Snew1',status)
1386	! write(nulprt,) subname,mpi_rank_local,rsize,rsizenwgts
1387
1388	cnt = 1
1389	if (local_timers_on) call oasis_timer_stop('map_read_ceg_read')
1390
1391	!----------------------------------------------------------------------------
1392	!> * On the root task
1393	!----------------------------------------------------------------------------
1394
1395	if (mytask == 0) then
1396
1397	!----------------------------------------------------------------------------
1398	!> * Initialize lsstart, lscount, and lspeloc, the sorted list of local indices
1399	!> * Sort via bubble sort or merge sort depending on size of array
1400	!----------------------------------------------------------------------------
1401
1402	if (local_timers_on) call oasis_timer_start('map_read_ceg_sort')
1403	lsize = size(mygsmap%start)
1404	allocate(lsstart(lsize),lscount(lsize),lspeloc(lsize),stat=status)
1405	if (status /= 0) call mct_perr_die(subName,':: allocate Lsstart',status)
1406
1407	if (lsize > 10) then
1408	! merge sort
1409	allocate(sortkey(lsize),stat=status)
1410	if (status /= 0) call mct_perr_die(subName,':: allocate sortkey',status)
1411	do n = 1,lsize
1412	lsstart(n) = mygsmap%start(n)
1413	lscount(n) = mygsmap%length(n)
1414	lspeloc(n) = mygsmap%pe_loc(n)
1415	sortkey(n) = n
1416	enddo
1417	call oasis_sys_sortI(lsize,lsstart,sortkey)
1418	call oasis_sys_sortIkey(lsize,lscount,sortkey)
1419	call oasis_sys_sortIkey(lsize,lspeloc,sortkey)
1420	deallocate(sortkey)
1421	else
1422	! bubble sort
1423	do n = 1,lsize
1424	found = .false.
1425	l1 = 1
1426	do while (.not.found .and. l1 < n)
1427	if (mygsmap%start(n) < lsstart(l1)) then
1428	do l2 = n, l1+1, -1
1429	lsstart(l2) = lsstart(l2-1)
1430	lscount(l2) = lscount(l2-1)
1431	lspeloc(l2) = lspeloc(l2-1)
1432	enddo
1433	found = .true.
1434	else
1435	l1 = l1 + 1
1436	endif
1437	enddo
1438	lsstart(l1) = mygsmap%start(n)
1439	lscount(l1) = mygsmap%length(n)
1440	lspeloc(l1) = mygsmap%pe_loc(n)
1441	enddo
1442	endif
1443
1444	do n = 1,size(mygsmap%start)-1
1445	if (lsstart(n) > lsstart(n+1)) then
1446	write(nulprt,*) subname,estr,'lsstart not properly sorted'
1447	call oasis_abort(file=__FILE__,line=__LINE__)
1448	endif
1449	enddo
1450	if (local_timers_on) call oasis_timer_stop('map_read_ceg_sort')
1451	if (local_timers_on) call oasis_timer_start('map_read_ceg_dist')
1452
1453	!----------------------------------------------------------------------------
1454	!> * Allocate arrays for reading data on the root
1455	!----------------------------------------------------------------------------
1456
1457	allocate(remaps(nwgts,rsize),stat=status)
1458	if (status /= 0) call mct_perr_die(subName,':: allocate remaps',status)
1459	allocate(SReadData(nwgts,rsize),RReadData(rsize),CReadData(rsize),stat=status)
1460	if (status /= 0) call mct_perr_die(subName,':: allocate SReadData',status)
1461	allocate(pesave(rsize),stat=status)
1462	if (status /= 0) call mct_perr_die(subName,':: allocate pesave',status)
1463	allocate(SDistData(nwgts,rsize),RDistData(rsize),CDistData(rsize), stat=status)
1464	if (status /= 0) call mct_perr_die(subName,':: allocate SDistData',status)
1465	allocate(cntrs(0:mpi_size_local), stat=status) ! Purposefully allocating one extra
1466	if (status /= 0) call mct_perr_die(subName,':: allocate cntrs',status)
1467
1468	!----------------------------------------------------------------------------
1469	!> * Loop over the chunks of weights data
1470	!----------------------------------------------------------------------------
1471
1472	do n = 1,nread
1473	start(1) = (n-1)*rsize + 1
1474	count(1) = min(rsize,ns-start(1)+1)
1475	start2(1) = 1
1476	count2(1) = nwgts
1477	start2(2) = start(1)
1478	count2(2) = count(1)
1479
1480	!----------------------------------------------------------------------------
1481	!> * Read chunk of data
1482	!----------------------------------------------------------------------------
1483
1484	! status = nf90_inq_varid (fid,'S' ,vid)
1485	status = nf90_inq_varid (fid,'remap_matrix' ,vid)
1486	if (status /= NF90_NOERR) then
1487	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1488	write(nulprt,*) subname,estr,'var not found = ','remap_matrix'
1489	call oasis_abort(file=__FILE__,line=__LINE__)
1490	endif
1491	! status = nf90_get_var(fid,vid,start,count,Sbuf)
1492	status = nf90_get_var(fid,vid,remaps,start2,count2)
1493	if (status /= NF90_NOERR) THEN
1494	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
1495	call oasis_flush(nulprt)
1496	endif
1497	SReadData(:,:) = remaps(:,:)
1498	! status = nf90_inq_varid (fid,'row',vid)
1499	status = nf90_inq_varid (fid,'dst_address',vid)
1500	if (status /= NF90_NOERR) then
1501	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1502	write(nulprt,*) subname,estr,'var not found = ','dst_address'
1503	call oasis_abort(file=__FILE__,line=__LINE__)
1504	endif
1505	status = nf90_get_var (fid,vid,RReadData,start,count)
1506	if (status /= NF90_NOERR) THEN
1507	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
1508	call oasis_flush(nulprt)
1509	endif
1510
1511	! status = nf90_inq_varid (fid,'col',vid)
1512	status = nf90_inq_varid (fid,'src_address',vid)
1513	if (status /= NF90_NOERR) then
1514	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1515	write(nulprt,*) subname,estr,'var not found = ','src_address'
1516	call oasis_abort(file=__FILE__,line=__LINE__)
1517	endif
1518	status = nf90_get_var (fid,vid,CReadData,start,count)
1519	if (status /= NF90_NOERR) THEN
1520	write(nulprt,*) subname,' nf90_strerror = ',TRIM(nf90_strerror(status))
1521	call oasis_flush(nulprt)
1522	endif
1523
1524	! Two stage process
1525	! 1. Count how many to send to each process
1526	! 2. Build the sending array up in process receiving order
1527
1528	! Initialize cntrs array for number of accumulations per process
1529	cntrs =0
1530	pesave = -99
1531
1532	!----------------------------------------------------------------------------
1533	!> * Determine which process owns each weight and count them
1534	!> * Determine offsets in the array
1535	!----------------------------------------------------------------------------
1536
1537	if (namwgtopt == "abort_on_bad_index") then
1538	abort_weight = .false.
1539	do m = 1,count(1)
1540	!--- check for bad weights
1541	if ((RReadData(m) <= 0 .or. RReadData(m) > nb .or. &
1542	CReadData(m) <= 0 .or. CReadData(m) > na) &
1543	! tcx weight = 0
1544	.and. (minval(SReadData(:,m)) /= 0._R8 .or. maxval(SReadData(:,m)) /= 0._R8) &
1545	) then
1546	abort_weight = .true.
1547	write(nulprt,'(3A,I12,A,I12,A,I12,A,G14.7,A,G14.7,A)') &
1548	subname,wstr,'BAD weight found in '//trim(filename), &
1549	m,'=id',CReadData(m),'=src',RReadData(m),'=dst',minval(SReadData(:,m)),'=minS',maxval(SReadData(:,m)),'=maxS'
1550	endif
1551	enddo
1552	if (abort_weight) then
1553	write(nulprt,*) subname,wstr,'BAD weight found, aborting'
1554	call oasis_abort(file=__FILE__,line=__LINE__)
1555	endif
1556	endif
1557
1558	do m = 1,count(1)
1559	!--- check for bad weights
1560	if ((namwgtopt(1:16) == "ignore_bad_index") .and. &
1561	(RReadData(m) <= 0 .or. RReadData(m) > nb .or. &
1562	CReadData(m) <= 0 .or. CReadData(m) > na)) then
1563	pe = -11
1564	! tcx weight = 0
1565	if (minval(SReadData(:,m)) /= 0._R8 .or. maxval(SReadData(:,m)) /= 0._R8) then
1566	if (OASIS_debug >= 2 .and. namwgtopt /= "ignore_bad_index_silently") then
1567	write(nulprt,'(3A,I12,A,I12,A,I12,A,G14.7,A,G14.7,A)') &
1568	subname,wstr,'BAD weight found in '//trim(filename), &
1569	m,'=id',CReadData(m),'=src',RReadData(m),'=dst',minval(SReadData(:,m)),'=minS',maxval(SReadData(:,m)),'=maxS'
1570	endif
1571	endif
1572	else if (newdom == 'src') then
1573	pe = get_cegindex(RReadData(m),lsstart,lscount,lspeloc)
1574	else if (newdom == 'dst') then
1575	pe = get_cegindex(CReadData(m),lsstart,lscount,lspeloc)
1576	endif
1577	pesave(m) = pe
1578
1579	! Copy data, incrementing index array, pe = 1 to mpi_size_local
1580	if (pe == -11) then
1581	! skip it, understood, get_cegindex will error with -99
1582	elseif (pe+1 < 1 .or. pe+1 > mpi_size_local) then
1583	! skip this case too
1584	! write(nulprt,*) subname,wstr,'get_cegindex search error', m,count(1),CReadData(m),RReadData(m),SReadData(:,m)
1585	else
1586	cntrs(pe+1) = cntrs(pe+1) + 1 ! Note incrementing 1->noprocs rather than 0->noprocs-1
1587	endif
1588
1589	end do
1590
1591	cntrs(0) = 1
1592	do pe = 1, mpi_size_local
1593	cntrs(pe) = cntrs(pe-1) + cntrs(pe)
1594	end do
1595
1596	!----------------------------------------------------------------------------
1597	!> * Determine which process owns each weight and fill arrays
1598	!----------------------------------------------------------------------------
1599
1600	do m = 1,count(1)
1601
1602	pe = pesave(m)
1603
1604	! Copy data, incrementing index array, pe = 1 to mpi_size_local
1605	if (pe+1 < 1 .or. pe+1 > mpi_size_local) then
1606	! skip it
1607	! write(nulprt,*) subname,'get_cegindex search error2', m,count(1),CReadData(m),RReadData(m),SReadData(:,m)
1608	else
1609	! Copy data, incrementing index array
1610	SDistData(:,cntrs(pe)) = SReadData(:,m)
1611	RDistData(cntrs(pe)) = RReadData(m)
1612	CDistData(cntrs(pe)) = CReadData(m)
1613
1614	cntrs(pe) = cntrs(pe) + 1
1615	endif
1616
1617	enddo ! count
1618
1619	!----------------------------------------------------------------------------
1620	!> * Send select data from root to other processes
1621	! Copy root process data into local array
1622	! Send data from root to other processes
1623	!----------------------------------------------------------------------------
1624
1625	if (cntrs(0).gt.1) then ! Need to do it differently if it is for me!
1626	reclen = cntrs(0)-1
1627	! Reallocate local weights arrays if they need to be bigger
1628	call augment_arrays(cnt, reclen, bsize, nwgts)
1629
1630	!--- now p0 copies straight into the ?new arrays
1631	Snew(1:nwgts,cnt:cnt+reclen-1) = SDistData(1:nwgts,1:reclen)
1632	Rnew(cnt:cnt+reclen-1) = RDistData(1:reclen)
1633	Cnew(cnt:cnt+reclen-1) = CDistData(1:reclen)
1634	cnt = cnt + reclen
1635
1636	endif
1637
1638	do pe = 1, mpi_size_local
1639	if (cntrs(pe).gt.cntrs(pe-1)) then
1640	! Send data out
1641	m = cntrs(pe)-cntrs(pe-1)
1642	! write(nulprt,*) subname, 'Sending [to, datalen] ', pe, m, cntrs(pe-1), cntrs(pe)
1643	call MPI_Send(m, 1, MPI_INTEGER, pe, 4000, mpicom, ierr)
1644	call MPI_Send(SDistData(1,cntrs(pe-1)), nwgts*m, MPI_REAL8, pe, 1000, mpicom, ierr)
1645	call MPI_Send(RDistData(cntrs(pe-1)), m, MPI_INTEGER, pe, 2000, mpicom, ierr)
1646	call MPI_Send(CDistData(cntrs(pe-1)), m, MPI_INTEGER, pe, 3000, mpicom, ierr)
1647	endif
1648	enddo
1649
1650	enddo ! nread
1651
1652	!----------------------------------------------------------------------------
1653	!> * Deallocate memory on root process
1654	!----------------------------------------------------------------------------
1655
1656	m = -1
1657	do pe = 1, mpi_size_local-1
1658	!write(nulprt,*) subname, 'Final sending [to, datalen] ', m, cntrs(m)
1659	call MPI_Send(m, 1, MPI_INTEGER, pe, 4000, mpicom, ierr)
1660	end do
1661
1662	! Free memory used during the distribution
1663	deallocate(lsstart,lscount,lspeloc, stat=status)
1664	if (status /= 0) call mct_perr_die(subName,':: deallocate lsstart',status)
1665	deallocate(pesave, stat=status)
1666	if (status /= 0) call mct_perr_die(subName,':: deallocate pesave',status)
1667	deallocate(SReadData,RReadData,CReadData, stat=status)
1668	if (status /= 0) call mct_perr_die(subName,':: deallocate SReadData',status)
1669	deallocate(SDistData,RDistData,CDistData, stat=status)
1670	if (status /= 0) call mct_perr_die(subName,':: deallocate SDistData',status)
1671	deallocate(cntrs, stat=status)
1672	if (status /= 0) call mct_perr_die(subName,':: deallocate cntrs',status)
1673	deallocate(remaps, stat=status)
1674	if (status /= 0) call mct_perr_die(subName,':: deallocate remaps',status)
1675
1676	if (local_timers_on) call oasis_timer_stop('map_read_ceg_dist')
1677
1678	!----------------------------------------------------------------------------
1679	!> * On non-root processes
1680	!----------------------------------------------------------------------------
1681
1682	else ! mytask == 0
1683
1684	if (local_timers_on) call oasis_timer_start('map_read_ceg_dist')
1685	!----------------------------------------------------------------------------
1686	!> * Receive data from root
1687	!----------------------------------------------------------------------------
1688
1689	call oasis_mpi_recv(reclen, 0, 4000, mpicom, subName//" MPI in reclen recv")
1690	! Check we haven't now had the last data
1691	do while (reclen.ne.-1)
1692
1693	!--- recv S, row, col on all other pes
1694	!write(nulprt,*) subname, 'Receiving [global-id, datalen] ', mpi_rank_global, reclen
1695	call MPI_Recv(Sbuf, reclen*nwgts, MPI_REAL8, 0, 1000, mpicom, MPI_STATUS_IGNORE, ierr)
1696	call MPI_Recv(Rbuf, reclen, MPI_INTEGER, 0, 2000, mpicom, MPI_STATUS_IGNORE, ierr)
1697	call MPI_Recv(Cbuf, reclen, MPI_INTEGER, 0, MPI_ANY_TAG, mpicom, mpistatus, ierr)
1698
1699	!----------------------------------------------------------------------------
1700	!> * Reallocate local weights arrays if they need to be bigger
1701	!----------------------------------------------------------------------------
1702
1703	call augment_arrays(cnt, reclen, bsize, nwgts)
1704
1705	!--- now each pe keeps what it has been sent
1706	Snew(1:nwgts,cnt:cnt+reclen-1) = Sbuf(1:nwgts,1:reclen)
1707	Rnew(cnt:cnt+reclen-1) = Rbuf(1:reclen)
1708	Cnew(cnt:cnt+reclen-1) = Cbuf(1:reclen)
1709	cnt = cnt + reclen
1710
1711	call oasis_mpi_recv(reclen, 0, 4000, mpicom, subName//" MPI in reclen recv")
1712	end do
1713	if (local_timers_on) call oasis_timer_stop('map_read_ceg_dist')
1714
1715	endif ! mytask == 0
1716
1717	if (local_timers_on) call oasis_timer_start('map_read_ceg_smati')
1718	! Fix cnt to be the length of the array
1719	cnt = cnt-1
1720
1721	!----------------------------------------------------------------------------
1722	!> * Clean up arrays
1723	!----------------------------------------------------------------------------
1724
1725	deallocate(Sbuf,Rbuf,Cbuf, stat=status)
1726	if (status /= 0) call mct_perr_die(subName,':: deallocate Sbuf',status)
1727
1728	!----------------------------------------------------------------------------
1729	!> * Initialize the mct sMat data type
1730	!----------------------------------------------------------------------------
1731
1732	! mct_sMat_init must be given the number of rows and columns that
1733	! would be in the full matrix. Nrows= size of output vector=nb.
1734	! Ncols = size of input vector = na.
1735	do n = 1,nwgts
1736	call mct_sMat_init(sMat(n), nb, na, cnt)
1737	enddo
1738
1739	igrow = mct_sMat_indexIA(sMat(1),'grow')
1740	igcol = mct_sMat_indexIA(sMat(1),'gcol')
1741	iwgt = mct_sMat_indexRA(sMat(1),'weight')
1742
1743	if (cnt /= 0) then
1744	do n = 1,nwgts
1745	sMat(n)%data%rAttr(iwgt ,1:cnt) = Snew(n,1:cnt)
1746	sMat(n)%data%iAttr(igrow,1:cnt) = Rnew(1:cnt)
1747	sMat(n)%data%iAttr(igcol,1:cnt) = Cnew(1:cnt)
1748	enddo
1749	endif
1750	if (local_timers_on) call oasis_timer_stop('map_read_ceg_smati')
1751
1752	!----------------------------------------------------------------------------
1753	!> * More clean up
1754	!----------------------------------------------------------------------------
1755
1756	if (local_timers_on) call oasis_timer_start('map_read_ceg_clean')
1757	deallocate(Snew,Rnew,Cnew, stat=status)
1758	if (status /= 0) call mct_perr_die(subName,':: deallocate new',status)
1759
1760	if (mytask == 0) then
1761	status = nf90_close(fid)
1762	if (status /= NF90_NOERR) then
1763	write(nulprt,*) subname,' nf90_strerror = ',trim(nf90_strerror(status))
1764	call oasis_abort(file=__FILE__,line=__LINE__)
1765	endif
1766	if (OASIS_debug >= 2) THEN
1767	write(nulprt,*) subname," ... done reading file"
1768	call oasis_flush(nulprt)
1769	endif
1770	endif
1771
1772	if (local_timers_on) call oasis_timer_stop('map_read_ceg_clean')
1773	call oasis_debug_exit(subname)
1774
1775	end subroutine oasis_map_sMatReaddnc_ceg
1776
1777	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1778	!!!!!!!!!!!!!!!!!!!!!!!!! Extend array and store data !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1779	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1780
1781	!> Function that increases temporary work array size of Snew, Rnew, Cnew
1782
1783	subroutine augment_arrays(cnt, reclen, bsize, nwgts)
1784
1785	! !USES:
1786
1787	integer, intent(inout) :: cnt !< elements in current array
1788	integer, intent(in) :: reclen !< elements of new data
1789	integer, intent(inout) :: bsize !< max size of current array
1790	integer, intent(in) :: nwgts !< number of weights in S
1791
1792	integer :: status ! netCDF routine return code
1793	!--- formats ---
1794	character(*),parameter :: subName = '(ceg_coupler_augment_arrays)'
1795
1796
1797	!--- ?new arrays need to be bigger?
1798	if (cnt+reclen > bsize) then
1799
1800	! write(nulprt,*) subname,mpi_rank_global, 'EXTEND'
1801
1802	!--- allocate old arrays and copy new into old
1803	allocate(Sold(1:nwgts,cnt),Rold(cnt),Cold(cnt),stat=status)
1804	if (status /= 0) call mct_perr_die(subName,':: allocate old',status)
1805	Sold(1:nwgts,1:cnt-1) = Snew(1:nwgts,1:cnt-1)
1806	Rold(1:cnt-1) = Rnew(1:cnt-1)
1807	Cold(1:cnt-1) = Cnew(1:cnt-1)
1808
1809	!--- reallocate new to bigger size, increase buffer by 50% (arbitrary)
1810	deallocate(Snew,Rnew,Cnew,stat=status)
1811	if (status /= 0) call mct_perr_die(subName,':: allocate new',status)
1812	bsize = 1.5 * (cnt+reclen)
1813	if (OASIS_debug > 15) write(nulprt,*) subname,' reallocate bsize to ',bsize
1814	allocate(Snew(nwgts,bsize),Rnew(bsize),Cnew(bsize),stat=status)
1815	if (status /= 0) call mct_perr_die(subName,':: allocate old',status)
1816
1817	!--- copy data back into new
1818	Snew(1:nwgts,1:cnt-1) = Sold(1:nwgts,1:cnt-1)
1819	Rnew(1:cnt-1) = Rold(1:cnt-1)
1820	Cnew(1:cnt-1) = Cold(1:cnt-1)
1821	deallocate(Sold,Rold,Cold,stat=status)
1822	if (status /= 0) call mct_perr_die(subName,':: deallocate old',status)
1823	endif
1824
1825	end subroutine augment_arrays
1826
1827	!------------------------------------------------------------
1828	! !BOP ===========================================================================
1829	!
1830	!> Function that checks whether an index is part of a start and count list
1831	!
1832	!> Does a binary search on a sorted start and count list to determine
1833	!> whether index is a value in the list. The list values consist of
1834	!> the values start(i):start(i)+count(i)-1 for all i.
1835	!
1836	! !DESCRIPTION:
1837	! Do a binary search to see if a value is contained in a list of
1838	! values. return true or false. starti must be monotonically
1839	! increasing, function does NOT check this.
1840	!
1841	! !INTERFACE: -----------------------------------------------------------------
1842
1843	logical function check_myindex(index,starti,counti)
1844
1845	! !USES:
1846
1847	!--- local kinds ---
1848	integer,parameter :: R8 = ip_double_p
1849	integer,parameter :: IN = ip_i4_p
1850
1851	! !INPUT/OUTPUT PARAMETERS:
1852
1853	integer(IN) :: index !< index to search
1854	integer(IN) :: starti(:) !< start list
1855	integer(IN) :: counti(:) !< count list
1856
1857	! !EOP
1858
1859	!--- local ---
1860	integer(IN) :: nl,nc,nr,ncprev
1861	integer(IN) :: lsize
1862	logical :: stopnow
1863
1864	!--- formats ---
1865	character(*),parameter :: subName = '(check_myindex) '
1866
1867	!-------------------------------------------------------------------------------
1868	!
1869	!-------------------------------------------------------------------------------
1870
1871	! call oasis_debug_enter(subname)
1872	check_myindex = .false.
1873
1874	lsize = size(starti)
1875	if (lsize < 1) then
1876	! call oasis_debug_exit(subname)
1877	return
1878	endif
1879
1880	nl = 0
1881	nr = lsize + 1
1882	nc = (nl+nr)/2
1883	stopnow = .false.
1884	do while (.not.stopnow)
1885	if (index < starti(nc)) then
1886	nr = nc
1887	elseif (index > (starti(nc) + counti(nc) - 1)) then
1888	nl = nc
1889	else
1890	check_myindex = .true.
1891	! call oasis_debug_exit(subname)
1892	return
1893	endif
1894	ncprev = nc
1895	nc = (nl + nr)/2
1896	if (nc == ncprev .or. nc < 1 .or. nc > lsize) stopnow = .true.
1897	enddo
1898
1899	check_myindex = .false.
1900
1901	! call oasis_debug_exit(subname)
1902
1903	end function check_myindex
1904
1905	!------------------------------------------------------------
1906	! !BOP ===========================================================================
1907	!
1908	!> Function that carrys out a binary search for index in list
1909	!
1910	! !DESCRIPTION:
1911	! Do a binary search to see if a value is contained in a list of
1912	! values. return value of processor that owns it
1913	! starti must be monotonically
1914	! increasing, function does NOT check this.
1915	!
1916	! !INTERFACE: -----------------------------------------------------------------
1917
1918	integer function get_cegindex(index,starti,counti,peloci)
1919
1920	! !USES:
1921
1922
1923	!--- local kinds ---
1924	integer,parameter :: R8 = ip_double_p
1925	integer,parameter :: IN = ip_i4_p
1926
1927	! !INPUT/OUTPUT PARAMETERS:
1928
1929	integer(IN) :: index !< index to search
1930	integer(IN) :: starti(:) !< start list
1931	integer(IN) :: counti(:) !< count list
1932	integer(IN) :: peloci(:) !< pe list
1933
1934	! !EOP
1935
1936
1937	!--- local ---
1938	integer(IN) :: nl,nc,nr,ncprev
1939	integer(IN) :: lsize
1940	logical :: stopnow
1941
1942	!--- formats ---
1943	character(*),parameter :: subName = '(get_cegindex) '
1944
1945	!-------------------------------------------------------------------------------
1946	!
1947	!-------------------------------------------------------------------------------
1948	! call oasis_debug_enter(subname)
1949
1950	get_cegindex = -99
1951	lsize = size(starti)
1952	if (lsize < 1) then
1953	! call oasis_debug_exit(subname)
1954	return
1955	endif
1956
1957	nl = 0
1958	nr = lsize + 1
1959	nc = (nl+nr)/2
1960	stopnow = .false.
1961	do while (.not.stopnow)
1962	if (index < starti(nc)) then
1963	nr = nc
1964	elseif (index > (starti(nc) + counti(nc) - 1)) then
1965	nl = nc
1966	else
1967	get_cegindex = peloci(nc)
1968	! call oasis_debug_exit(subname)
1969	return
1970	endif
1971	ncprev = nc
1972	nc = (nl + nr)/2
1973	if (nc == ncprev .or. nc < 1 .or. nc > lsize) stopnow = .true.
1974	enddo
1975
1976
1977	! call oasis_debug_exit(subname)
1978
1979
1980	end function get_cegindex
1981
1982	!===============================================================================
1983	END MODULE mod_oasis_map
1984
1985

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/lib/psmile/src/mod_oasis_map.F90 @ 6331

Download in other formats: