Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

remap_conserv.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: 129.6 KB

Line
1	!****
2	! ************************
3	! * OASIS MODULE *
4	! * ------------ *
5	! ************************
6	!****
7	!***********************************************************************
8	! This module belongs to the SCRIP library. It is modified to run
9	! within OASIS.
10	! Main modifications:
11	! - Some allocated array will be freed in the end to allow multiple
12	! calls of SCRIP
13	! - Introduction of a logical flag to distinguish between the first
14	! and following calls of scrip conservative remapping
15	! - Masking of overlapping grid points for source and target grid
16	! For these points, links and weights of overlapped point are used
17	!
18	! Modified by V. Gayler, M&D 20.09.2001
19	! Modified by D. Declat, CERFACS 27.06.2002
20	! Modified by A. Piacentini, CERFACS 24.01.2018
21	!***********************************************************************
22	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23	!
24	! this module contains necessary routines for computing addresses
25	! and weights for a conservative interpolation between any two
26	! grids on a sphere. the weights are computed by performing line
27	! integrals around all overlap regions of the two grids. see
28	! Dukowicz and Kodis, SIAM J. Sci. Stat. Comput. 8, 305 (1987) and
29	! Jones, P.W. Monthly Weather Review (submitted).
30	!
31	!-----------------------------------------------------------------------
32	!
33	! CVS:$Id: remap_conserv.f 1811 2008-12-19 08:41:41Z valcke $
34	!
35	! Copyright (c) 1997, 1998 the Regents of the University of
36	! California.
37	!
38	! This software and ancillary information (herein called software)
39	! called SCRIP is made available under the terms described here.
40	! The software has been approved for release with associated
41	! LA-CC Number 98-45.
42	!
43	! Unless otherwise indicated, this software has been authored
44	! by an employee or employees of the University of California,
45	! operator of the Los Alamos National Laboratory under Contract
46	! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S.
47	! Government has rights to use, reproduce, and distribute this
48	! software. The public may copy and use this software without
49	! charge, provided that this Notice and any statement of authorship
50	! are reproduced on all copies. Neither the Government nor the
51	! University makes any warranty, express or implied, or assumes
52	! any liability or responsibility for the use of this software.
53	!
54	! If software is modified to produce derivative works, such modified
55	! software should be clearly marked, so as not to confuse it with
56	! the version available from Los Alamos National Laboratory.
57	!
58	!***********************************************************************
59
60	module remap_conservative
61
62	!-----------------------------------------------------------------------
63
64	use kinds_mod ! defines common data types
65	use constants ! defines common constants
66	use timers ! module for timing
67	use grids ! module containing grid information
68	use remap_vars ! module containing remap information
69	use mod_oasis_flush
70	!$ use omp_lib
71
72	implicit none
73
74	!-----------------------------------------------------------------------
75	!
76	! module variables
77	!
78	!-----------------------------------------------------------------------
79
80	real (kind=dbl_kind) :: north_thresh = 2.00_dbl_kind, & ! threshold for coord transf.
81	south_thresh =-2.00_dbl_kind ! threshold for coord transf.
82
83	integer (kind=int_kind) :: il_nbthreads = 1
84
85	!***********************************************************************
86
87	contains
88
89	!***********************************************************************
90
91	subroutine remap_conserv(mpi_comm_map, mpi_size_map, mpi_rank_map, mpi_root_map)
92
93	!-----------------------------------------------------------------------
94	!
95	! this routine traces the perimeters of every grid cell on each
96	! grid checking for intersections with the other grid and computing
97	! line integrals for each subsegment.
98	!
99	!-----------------------------------------------------------------------
100
101	!-----------------------------------------------------------------------
102	!
103	! input variables
104	!
105	!-----------------------------------------------------------------------
106
107	integer (kind=int_kind) :: mpi_comm_map, mpi_rank_map, mpi_size_map, mpi_root_map
108
109	!-----------------------------------------------------------------------
110	!
111	! local variables
112	!
113	!-----------------------------------------------------------------------
114
115	integer (kind=int_kind), parameter :: max_subseg = 10000 ! max number of subsegments per segment
116	! to prevent infinite loop
117
118	integer (kind=int_kind) :: num_srch_cells ! num cells in restricted search arrays
119
120	integer (kind=int_kind), dimension(:), allocatable :: srch_add ! global address of cells in srch arrays
121
122	integer (kind=int_kind) :: grid1_add, & ! current linear address for grid1 cell
123	grid2_add, & ! current linear address for grid2 cell
124	min_add, & ! addresses for restricting search of
125	max_add, & ! destination grid
126	min_cell, & ! addresses for restricting search of
127	max_cell, & ! destination grid
128	n, nwgt, & ! generic counters
129	corner, & ! corner of cell that segment starts from
130	next_corn, & ! corner of cell that segment ends on
131	num_subseg ! number of subsegments
132
133	logical (kind=log_kind) :: lcoinc, & ! flag for coincident segments
134	lrevers, & ! flag for reversing direction of segment
135	lbegin, & ! flag for first integration of a segment
136	ltake, & ! flag for longitude intersection
137	full, ll_debug ! for debug outputs
138
139
140	real (kind=dbl_kind) :: &
141	intrsct_lat, intrsct_lon, & ! lat/lon of next intersect
142	beglat, endlat, beglon, endlon, & ! endpoints of current seg.
143	norm_factor, & ! factor for normalizing wts
144	norm_correc, & ! true area correction for norm wts
145	delta, & ! precision
146	r2d
147
148	real (kind=dbl_kind) :: dgbb1, dgbb2, dgbb3, dgbb4
149	integer (kind=int_kind) :: dgbbp
150
151	integer (kind=int_kind) :: num_links_map1_sweep1
152
153	real (kind=dbl_kind), dimension(:), allocatable :: grid2_centroid_lat, grid2_centroid_lon, & ! centroid coords
154	grid1_centroid_lat, grid1_centroid_lon ! on each grid
155
156	real (kind=dbl_kind), dimension(2) :: begseg ! begin lat/lon for
157	! full segment
158
159	real (kind=dbl_kind), dimension(6) :: weights ! local wgt array
160	real (kind=dbl_kind) :: intrsct_lat_off, intrsct_lon_off ! lat/lon coords offset
161
162	logical :: ll_timing=.true.
163
164	integer (kind=int_kind) :: il_splitsize
165	integer (kind=int_kind) :: ib_proc
166	integer (kind=int_kind) :: ib_thread
167	integer (kind=int_kind) :: buff_base
168	integer (kind=int_kind), dimension(:), allocatable :: ila_mpi_sz
169	integer (kind=int_kind), dimension(:), allocatable :: ila_mpi_mn
170	integer (kind=int_kind), dimension(:), allocatable :: ila_mpi_mx
171	integer (kind=int_kind), dimension(:), allocatable :: ila_thr_sz
172	integer (kind=int_kind), dimension(:), allocatable :: ila_thr_mn
173	integer (kind=int_kind), dimension(:), allocatable :: ila_thr_mx
174	integer (kind=int_kind), dimension(:), allocatable :: ila_num_links_mpi
175	integer (kind=int_kind), dimension(:,:), allocatable :: ila_req_mpi
176	integer (kind=int_kind), dimension(:,:,:), allocatable :: ila_sta_mpi
177	character (LEN=14) :: cl_envvar
178	integer (kind=int_kind) :: il_envthreads, il_err, est_num_neighbors
179	integer (kind=int_kind) :: grid1_np, grid1_sp, grid2_np, grid2_sp
180
181	#ifdef TREAT_OVERLAY
182	integer (kind=int_kind) :: cell
183	integer (kind=int_kind), dimension(:), allocatable :: ila_idx
184	integer (kind=int_kind), dimension(:), allocatable :: grid2_overlap ! overlapping points
185	#endif
186
187	!-----------------------------------------------------------------------
188	!
189	if (nlogprt .ge. 2) then
190	write (UNIT = nulou,FMT = *)' '
191	write (UNIT = nulou,FMT = *)'Entering routine remap_conserv'
192	call OASIS_FLUSH_SCRIP(nulou)
193	endif
194
195	weights = 0.
196
197	est_num_neighbors = 4
198	!
199	!-----------------------------------------------------------------------
200	!
201	! initialize centroid arrays
202	!
203	!-----------------------------------------------------------------------
204
205	allocate( grid1_centroid_lat(grid1_size), grid1_centroid_lon(grid1_size), &
206	grid2_centroid_lat(grid2_size), grid2_centroid_lon(grid2_size))
207
208	grid1_centroid_lat = zero
209	grid1_centroid_lon = zero
210	grid2_centroid_lat = zero
211	grid2_centroid_lon = zero
212
213	!-----------------------------------------------------------------------
214	!
215	! integrate around each cell on grid1
216	!
217	!-----------------------------------------------------------------------
218
219	#ifdef TREAT_OVERLAY
220
221	if (ll_timing) call timer_start(2,'remap_conserv overlay')
222
223	! Coordinate equality tolerance
224
225	delta = epsilon(1.)
226
227	! Check overlapping point of the source grid
228
229	if (nlogprt .ge. 2) then
230	write(nulou,*) 'Check overlapping point of the source grid'
231	call OASIS_FLUSH_SCRIP(nulou)
232	endif
233
234	! Initialise array that contains addresses of overlap grid point
235
236	do grid1_add = 1,grid1_size
237	grid1_add_repl1(grid1_add)=grid1_add
238	end do
239
240	allocate(ila_idx(grid1_size))
241
242	ila_idx(:) = [(n,n=1,grid1_size)]
243
244	! Sort the source grid indexes by lon (varying first) and lat
245
246	call hpsort_eps (grid1_size, grid1_center_lon, grid1_center_lat, ila_idx, delta)
247
248	! Check neighbours equality (with tolerance)
249
250	do cell = 1, grid1_size-1
251	n = max(ila_idx(cell),ila_idx(cell+1))
252	grid1_add = min(ila_idx(cell),ila_idx(cell+1))
253	if ( (abs(grid1_center_lon(grid1_add)- grid1_center_lon(n))<delta).and. &
254	(abs(grid1_center_lat(grid1_add)- grid1_center_lat(n))<delta)) then
255	if (grid1_mask(n) .or. grid1_mask(grid1_add)) then
256	if (grid1_mask(n)) then
257	grid1_mask(n) = .false.
258	grid1_mask(grid1_add) = .true.
259	grid1_add_repl1(grid1_add) = n
260	endif
261	endif
262	end if
263	end do
264
265	deallocate(ila_idx)
266
267	! Check overlapping point of the target grid
268
269	if (nlogprt .ge. 2) then
270	write(nulou,*) 'Check overlapping point of the target grid'
271	call OASIS_FLUSH_SCRIP(nulou)
272	endif
273
274	! Initialise array that contains addresses of overlap grid point
275
276	allocate(grid2_overlap(grid2_size))
277	grid2_overlap = -1
278
279	allocate(ila_idx(grid2_size))
280
281	ila_idx(:) = [(n,n=1,grid2_size)]
282
283	! Sort the destination grid indexes by lon (varying first) and lat
284
285	call hpsort_eps (grid2_size, grid2_center_lon, grid2_center_lat, ila_idx, delta)
286
287	! Check neighbours equality (with tolerance)
288
289	do cell = 1, grid2_size-1
290	n = max(ila_idx(cell),ila_idx(cell+1))
291	grid2_add = min(ila_idx(cell),ila_idx(cell+1))
292	if ( (abs(grid2_center_lon(grid2_add)- grid2_center_lon(n))<delta).and. &
293	(abs(grid2_center_lat(grid2_add)- grid2_center_lat(n))<delta)) then
294	grid2_overlap(grid2_add) = n
295	end if
296	end do
297
298	deallocate(ila_idx)
299
300	if (ll_timing) call timer_stop(2)
301
302	#endif
303
304	allocate(ila_mpi_mn(mpi_size_map), ila_mpi_mx(mpi_size_map) )
305
306	if (mpi_size_map .gt. 1) then
307
308	allocate(ila_mpi_sz(mpi_size_map))
309	il_splitsize = grid1_size
310	ila_mpi_sz(:) = floor(real(il_splitsize)/mpi_size_map)
311	ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) = ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) + 1
312
313	ila_mpi_mn(1) = 1
314	ila_mpi_mx(1) = ila_mpi_sz(1)
315
316	do ib_proc = 2, mpi_size_map
317	ila_mpi_mn(ib_proc) = sum(ila_mpi_sz(1:ib_proc-1)) + 1
318	ila_mpi_mx(ib_proc) = sum(ila_mpi_sz(1:ib_proc))
319	end do
320
321	deallocate(ila_mpi_sz)
322
323	else
324
325	ila_mpi_mn(1) = 1
326	ila_mpi_mx(1) = grid1_size
327
328	endif
329
330
331	call get_environment_variable(name='OASIS_OMP_NUM_THREADS', value=cl_envvar, status=il_err)
332
333	if ( il_err .ne. 0) then
334	il_envthreads = 0
335	else
336	read(cl_envvar,*) il_envthreads
337	end if
338
339	if (ll_timing) call timer_start(3,'remap_conserv grid1 sweep')
340
341	!$OMP PARALLEL NUM_THREADS(il_envthreads) DEFAULT(NONE) &
342	!$OMP SHARED(il_envthreads) &
343	!$OMP SHARED(num_wts) &
344	#ifdef TREAT_OVERLAY
345	!$OMP SHARED(grid2_overlap) &
346	#endif
347
348	!$OMP SHARED(grid1_corners,grid2_corners) &
349	!$OMP SHARED(grid1_bound_box,grid2_bound_box) &
350	!$OMP SHARED(grid1_bbox_per,grid2_bbox_per) &
351	!$OMP SHARED(grid2_size,bin_addr1,bin_addr2) &
352	!$OMP SHARED(est_num_neighbors,num_srch_bins) &
353	!$OMP SHARED(il_nbthreads) &
354	!$OMP SHARED(nulou,sga_remap) &
355	!$OMP SHARED(grid2_center_lat) &
356	!$OMP SHARED(grid1_center_lat,grid1_center_lon) &
357	!$OMP SHARED(grid2_center_lon) &
358	!$OMP SHARED(grid1_corner_lat,grid1_corner_lon) &
359	!$OMP SHARED(grid2_corner_lat,grid2_corner_lon) &
360	!$OMP SHARED(grid1_mask) &
361	!$OMP SHARED(mpi_rank_map,mpi_root_map,ila_mpi_mn,ila_mpi_mx) &
362	!$OMP SHARED(ila_thr_sz,ila_thr_mn,ila_thr_mx) &
363	!$OMP REDUCTION(+:grid2_frac) &
364	!!$OMP SHARED(grid2_frac) &
365	!$OMP SHARED(grid1_area) &
366	!$OMP SHARED(grid1_centroid_lat) &
367	!$OMP SHARED(grid1_centroid_lon) &
368	!$OMP PRIVATE(nlogprt) &
369	!$OMP PRIVATE(grid1_add,n,srch_add,corner) &
370	!$OMP PRIVATE(dgbb1,dgbb2,dgbb3,dgbb4,dgbbp) &
371	!$OMP PRIVATE(min_add,max_add) &
372	!$OMP PRIVATE(min_cell,max_cell,grid2_add,num_srch_cells) &
373	!$OMP PRIVATE(next_corn,beglat,beglon,endlat,endlon,lrevers) &
374	!$OMP PRIVATE(begseg,lbegin,num_subseg) &
375	!$OMP PRIVATE(intrsct_lat,intrsct_lon,lcoinc,weights,ll_debug) &
376	!$OMP PRIVATE(intrsct_lat_off,intrsct_lon_off) &
377	!$OMP PRIVATE(ib_thread,il_splitsize)
378
379	!$OMP SINGLE
380
381	il_nbthreads = 1
382	!$ il_nbthreads = OMP_GET_NUM_THREADS ()
383
384	allocate(ila_thr_mn(il_nbthreads))
385	allocate(ila_thr_mx(il_nbthreads))
386
387	if (il_nbthreads .gt. 1) then
388
389	nlogprt = 0
390
391	allocate(ila_thr_sz(il_nbthreads))
392	il_splitsize = ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1
393	ila_thr_sz(:) = floor(real(il_splitsize)/il_nbthreads)
394	ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) = ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) + 1
395
396	ila_thr_mn(1) = ila_mpi_mn(mpi_rank_map+1)
397	ila_thr_mx(1) = ila_thr_mn(1) + ila_thr_sz(1) - 1
398
399	do ib_thread = 2, il_nbthreads
400	ila_thr_mn(ib_thread) = ila_mpi_mn(mpi_rank_map+1) + sum(ila_thr_sz(1:ib_thread-1))
401	ila_thr_mx(ib_thread) = ila_thr_mn(ib_thread) + ila_thr_sz(ib_thread) - 1
402	end do
403
404	allocate(sga_remap(il_nbthreads))
405
406	do ib_thread = 1, il_nbthreads
407	il_splitsize = est_num_neighbors*ila_thr_sz(ib_thread)
408	sga_remap(ib_thread)%max_links = il_splitsize
409	sga_remap(ib_thread)%num_links = 0
410	allocate(sga_remap(ib_thread)%grid1_add(il_splitsize))
411	allocate(sga_remap(ib_thread)%grid2_add(il_splitsize))
412	allocate(sga_remap(ib_thread)%wts(num_wts,il_splitsize))
413	end do
414
415	deallocate(ila_thr_sz)
416
417	else
418
419	ila_thr_mn(1) = ila_mpi_mn(mpi_rank_map+1)
420	ila_thr_mx(1) = ila_mpi_mx(mpi_rank_map+1)
421
422	end if
423	!$OMP END SINGLE
424
425	! Sweeps
426
427
428	if (nlogprt .ge. 2) then
429	write (nulou,*) 'grid1 sweep '
430	call OASIS_FLUSH_SCRIP(nulou)
431	endif
432
433	!-----------------------------------------------------------------------
434	!
435	! integrate around each cell on grid1
436	!
437	!-----------------------------------------------------------------------
438	!$OMP DO SCHEDULE(STATIC,1)
439	thread_loop1: do ib_thread = 1, il_nbthreads
440
441	ll_debug = .false.
442
443	grid_loop1: do grid1_add = ila_thr_mn(ib_thread), ila_thr_mx(ib_thread)
444
445	if (ll_debug) then
446	write(nulou,*) 'grid1_add=', grid1_add
447	call OASIS_FLUSH_SCRIP(nulou)
448	endif
449
450	!***
451	!*** destination grid bounding box (preloaded)
452	!***
453
454	dgbb1 = grid1_bound_box(1,grid1_add)
455	dgbb2 = grid1_bound_box(2,grid1_add)
456	dgbb3 = grid1_bound_box(3,grid1_add)
457	dgbb4 = grid1_bound_box(4,grid1_add)
458	dgbbp = grid1_bbox_per(grid1_add)
459
460	!***
461	!*** restrict searches first using search bins
462	!***
463
464	min_add = grid2_size
465	max_add = 1
466	do n=1,num_srch_bins
467	if (grid1_add >= bin_addr1(1,n) .and. grid1_add <= bin_addr1(2,n)) then
468	min_add = min(min_add, bin_addr2(1,n))
469	max_add = max(max_add, bin_addr2(2,n))
470	endif
471	end do
472
473	!***
474	!*** select the search cells
475	!***
476
477	num_srch_cells = 0
478	min_cell = max_add
479	max_cell = min_add
480
481	gather1: DO grid2_add = min_add, max_add
482
483	!***
484	!*** further restrict searches using bounding boxes
485	!***
486
487	if ( grid2_bound_box(1,grid2_add) <= dgbb2 ) then
488	if ( grid2_bound_box(2,grid2_add) >= dgbb1 ) then
489	if ( lf_loncross(dgbbp,dgbb3,dgbb4,grid2_add,&
490	grid2_bbox_per,grid2_bound_box) ) then
491
492	num_srch_cells = num_srch_cells+1
493	min_cell = MIN(min_cell, grid2_add)
494	max_cell = MAX(max_cell, grid2_add)
495
496	end if
497	end if
498	end if
499
500	end do gather1
501
502	!***
503	!*** store the search cells array
504	!***
505
506	if (num_srch_cells .ne. 0) then
507
508	allocate(srch_add(num_srch_cells))
509
510	n=0
511	pick1: DO grid2_add = min_cell,max_cell
512
513	if ( grid2_bound_box(1,grid2_add) <= dgbb2 ) then
514	if ( grid2_bound_box(2,grid2_add) >= dgbb1 ) then
515	if ( lf_loncross(dgbbp,dgbb3,dgbb4,grid2_add,&
516	grid2_bbox_per,grid2_bound_box) ) then
517
518	n = n+1
519	srch_add(n) = grid2_add
520
521	end if
522	end if
523	end if
524
525	end do pick1
526
527	end if
528
529	if (ll_debug) then
530	write(nulou,*)' '
531	write(nulou,)' * Grid1 cell and associated search cells **'
532	write(nulou,1110) grid1_add, grid1_center_lon(grid1_add), grid1_center_lat(grid1_add)
533	DO corner = 1,grid1_corners
534	write(nulou,1111) corner, grid1_corner_lon(corner,grid1_add), grid1_corner_lat(corner,grid1_add)
535	enddo
536	write(nulou,1969) grid1_add, grid1_bound_box(1,grid1_add), grid1_bound_box(2,grid1_add)
537	write(nulou,1971) grid1_add, grid1_bound_box(3,grid1_add), grid1_bound_box(4,grid1_add)
538	write(nulou,*) ' num_srch_cells=', num_srch_cells
539	write(nulou,*) ' '
540	if (num_srch_cells .ne. 0) then
541	write(nulou,*) ' srch_add(:)=', srch_add(1:num_srch_cells)
542	do n=1, num_srch_cells
543	do corner = 1,grid2_corners
544	write(nulou,1112) n, grid2_corner_lon(corner,srch_add(n)), grid2_corner_lat(corner,srch_add(n))
545	end do
546	end do
547	end if
548	write(nulou,)' **************************************'
549	write(nulou,*)' '
550	call OASIS_FLUSH_SCRIP(nulou)
551	endif
552	1110 format (' grid1 center, lon, lat = ', I8, 2X, F12.4, 2X, F12.4)
553	1111 format (' grid1 corner, lon, lat = ', I8, 2X, F12.4, 2X, F12.4)
554	1112 format (' srch cell, lon, lat = ', I8, 2X, F12.4, 2X, F12.4)
555	1969 format (' grid1 index, bbox lat = ', I8, 2X, F12.4, 2X, F12.4)
556	1971 format (' grid1 index, bbox lon = ', I8, 2X, F12.4, 2X, F12.4)
557
558	if (num_srch_cells .ne. 0) then
559
560	!***
561	!*** integrate around this cell
562	!***
563
564	do corner = 1,grid1_corners
565	next_corn = mod(corner,grid1_corners) + 1
566
567	!***
568	!*** define endpoints of the current segment
569	!***
570
571	beglat = grid1_corner_lat(corner,grid1_add)
572	beglon = grid1_corner_lon(corner,grid1_add)
573	endlat = grid1_corner_lat(next_corn,grid1_add)
574	endlon = grid1_corner_lon(next_corn,grid1_add)
575	lrevers = .false.
576
577	!***
578	!*** to ensure exact path taken during both
579	!*** sweeps, always integrate segments in the same
580	!*** direction (SW to NE).
581	!***
582
583	if ((endlat < beglat) .or. &
584	(endlat == beglat .and. endlon < beglon)) then
585	beglat = grid1_corner_lat(next_corn,grid1_add)
586	beglon = grid1_corner_lon(next_corn,grid1_add)
587	endlat = grid1_corner_lat(corner,grid1_add)
588	endlon = grid1_corner_lon(corner,grid1_add)
589	lrevers = .true.
590	if (ll_debug) write(nulou, *) ' sweep1 LREVERS TRUE'
591	endif
592	begseg(1) = beglat
593	begseg(2) = beglon
594	lbegin = .true.
595	num_subseg = 0
596
597	!***
598	!*** if this is a constant-longitude segment, skip the rest
599	!*** since the line integral contribution will be zero.
600	!***
601	if (ll_debug) then
602	if (endlon .eq. beglon) then
603	write(nulou,1113) beglon, beglat
604	write(nulou,1114) endlon, endlat
605	write(nulou, *)' sweep1 endlon == beglon; skip segment'
606	write(nulou,*) ' '
607	endif
608	endif
609	1113 format (' endlon == beglon; beglon, beglat = ', 2X, F12.4, 2X, F12.4)
610	1114 format (' endlon == beglon; endlon, endlat = ', 2X, F12.4, 2X, F12.4)
611
612	if (endlon /= beglon) then
613	!***
614	!*** integrate along this segment, detecting intersections
615	!*** and computing the line integral for each sub-segment
616	!***
617
618	do while (beglat /= endlat .or. beglon /= endlon)
619	!***
620	!*** prevent infinite loops if integration gets stuck
621	!*** near cell or threshold boundary
622	!***
623
624	num_subseg = num_subseg + 1
625	if (num_subseg > max_subseg) then
626	write(nulou,*) 'ERROR=>integration stalled:'
627	write(nulou,*) 'num_subseg exceeded limit'
628	write(nulou,*) '=>Verify corners in grids.nc, especially'
629	write(nulou,*) 'if calculated by OASIS routine corners'
630	write(nulou,*) 'integration stalled: num_subseg exceeded limit'
631	call OASIS_FLUSH_SCRIP(nulou)
632	stop
633	endif
634
635	!***
636	!*** find next intersection of this segment with a grid
637	!*** line on grid 2.
638	!***
639
640	if (ll_debug) then
641	write(nulou,*) ' '
642	write(nulou,1115) beglon, beglat
643	write(nulou,1116) endlon, endlat
644	write(nulou,*) ' '
645	call OASIS_FLUSH_SCRIP(nulou)
646	endif
647	1115 format (' avant intersection; beglon, beglat = ', 2X, F12.4, 2X, F12.4)
648	1116 format (' avant intersection; endlon, endlat = ', 2X, F12.4, 2X, F12.4)
649
650	call intersection(grid2_add,intrsct_lat,intrsct_lon, &
651	lcoinc, beglat, beglon, endlat, endlon, begseg, &
652	lbegin, lrevers, &
653	grid2_corners, num_srch_cells, &
654	srch_add, grid2_corner_lon, grid2_corner_lat, &
655	intrsct_lat_off,intrsct_lon_off)
656
657	if (ll_debug) then
658	write(nulou,*) ' After call intersection, grid2_add', grid2_add
659	write(nulou,1117) beglon, beglat
660	write(nulou,1118) endlon, endlat
661	write(nulou,1119) intrsct_lon, intrsct_lat
662	write(nulou,*) ' '
663	call OASIS_FLUSH_SCRIP(nulou)
664	endif
665	1117 format(' after intersection; beglon, beglat = ', 2X, F12.4, 2X, F12.4, 2X, F12.4, 2X, F12.4)
666	1118 format (' after intersection; endlon, endlat = ', 2X, F12.4, 2X, F12.4, 2X, F12.4, 2X, F12.4)
667	1119 format (' after intersection; intrsct_lon, intrsct_lat = ', 2X, F12.4, 2X, F12.4, 2X, F12.4, 2X, F12.4)
668
669	lbegin = .false.
670
671	!***
672	!*** compute line integral for this subsegment.
673	!***
674
675	if (grid2_add /= 0) then
676	call line_integral(weights, num_wts, &
677	beglon, intrsct_lon, beglat, intrsct_lat, &
678	grid1_center_lon(grid1_add), &
679	grid2_center_lon(grid2_add))
680
681	if (ll_debug) then
682	write(nulou,*) ' A1) WEIGHTS for this subsegment =', weights(1)
683	write(nulou,*) ' '
684	call OASIS_FLUSH_SCRIP(nulou)
685	endif
686
687	else
688	call line_integral(weights, num_wts, &
689	beglon, intrsct_lon, beglat, intrsct_lat, &
690	grid1_center_lon(grid1_add), &
691	grid1_center_lon(grid1_add))
692
693	if (ll_debug) then
694	write(nulou,*) ' B1) WEIGHTS for this subsegment =', weights(1)
695	write(nulou,*) ' '
696	call OASIS_FLUSH_SCRIP(nulou)
697	endif
698
699	endif
700
701	!***
702	!*** if integrating in reverse order, change
703	!*** sign of weights
704	!***
705
706	if (lrevers) then
707	weights = -weights
708	if (ll_debug) then
709	write(nulou,*) ' LREVERS; WEIGHTS for this subsegment =', weights(1)
710	write(nulou,*) ' '
711	endif
712	endif
713
714
715	!***
716	!*** store the appropriate addresses and weights.
717	!*** also add contributions to cell areas and centroids.
718	!***
719	1120 FORMAT (' STORE add1,add2,blon,blat,ilon,ilat,WEIGHTS=', 1X,I8,1X,I8,1X,F12.8,1X,F12.8,1X,F12.8,1X,F12.8,1X, E16.8)
720	1121 FORMAT (' overlap STORE grid1_add, grid2_add, WEIGHTS=', 1X,I8,1X,I8,1X,F12.8,1X,F12.8,1X,F12.8,1X,F12.8,1X, E16.8)
721	1122 FORMAT (' lfracnei STORE grid1_add, grid2_add, WEIGHTS=', 1X,I8,1X,I8,1X,F12.8,1X,F12.8,1X,F12.8,1X,F12.8,1X, E16.8)
722	if (grid2_add /= 0) then
723	if (grid1_mask(grid1_add)) then
724	call store_link_cnsrv(grid1_add, grid2_add, weights, ib_thread)
725
726	if (ll_debug) then
727	write(nulou,*) ' after store_link_cnsrv norm1'
728	write(nulou,1120) grid1_add, grid2_add,beglon, beglat, intrsct_lon,intrsct_lat,weights(1)
729	endif
730	#ifdef TREAT_OVERLAY
731	if (grid2_overlap(grid2_add)/=-1) then
732	call store_link_cnsrv(grid1_add, grid2_overlap(grid2_add), weights, ib_thread)
733
734	if (ll_debug) then
735	write(nulou,*) ' after store_link_cnsrv overlap1'
736	write(nulou,1121) grid1_add, grid2_add,beglon, beglat, intrsct_lon,intrsct_lat,weights(1)
737	endif
738	endif
739	#endif
740
741	!EM grid1_frac(grid1_add) = grid1_frac(grid1_add) + weights(1)
742	grid2_frac(grid2_add) = grid2_frac(grid2_add) + weights(num_wts+1)
743	#ifdef TREAT_OVERLAY
744	if (grid2_overlap(grid2_add)/=-1) &
745	grid2_frac(grid2_overlap(grid2_add)) = &
746	grid2_frac(grid2_overlap(grid2_add)) + weights(num_wts+1)
747	#endif
748	endif
749	endif
750
751	grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1)
752	grid1_centroid_lat(grid1_add) = grid1_centroid_lat(grid1_add) + weights(2)
753	grid1_centroid_lon(grid1_add) = grid1_centroid_lon(grid1_add) + weights(3)
754
755	!***
756	!*** reset beglat and beglon for next subsegment.
757	!***
758
759	beglat = intrsct_lat
760	beglon = intrsct_lon
761
762	end do
763
764	endif
765
766	!***
767	!*** end of segment
768	!***
769
770	end do
771
772	!***
773	!*** finished with this cell: deallocate search array and
774	!*** start on next cell
775
776	deallocate(srch_add)
777
778	end if ! num_srch_cells .NE. 0
779
780	end do grid_loop1
781
782	end do thread_loop1
783	!$OMP END DO
784
785	!$OMP END PARALLEL
786
787	if (il_nbthreads .gt. 1) then
788
789	sga_remap(1)%start_pos = 1
790	il_splitsize = sga_remap(1)%num_links
791	do ib_thread = 2, il_nbthreads
792	il_splitsize = il_splitsize + sga_remap(ib_thread)%num_links
793	sga_remap(ib_thread)%start_pos = sga_remap(ib_thread-1)%start_pos + &
794	sga_remap(ib_thread-1)%num_links
795	end do
796
797	num_links_map1 = il_splitsize
798	if (num_links_map1 > max_links_map1) &
799	call resize_remap_vars(1,num_links_map1-max_links_map1)
800
801	do ib_thread = 1, il_nbthreads
802	grid1_add_map1(sga_remap(ib_thread)%start_pos: &
803	sga_remap(ib_thread)%start_pos+ &
804	sga_remap(ib_thread)%num_links-1) = &
805	sga_remap(ib_thread)%grid1_add(1:sga_remap(ib_thread)%num_links)
806	grid2_add_map1(sga_remap(ib_thread)%start_pos: &
807	sga_remap(ib_thread)%start_pos+ &
808	sga_remap(ib_thread)%num_links-1) = &
809	sga_remap(ib_thread)%grid2_add(1:sga_remap(ib_thread)%num_links)
810	wts_map1 (:,sga_remap(ib_thread)%start_pos: &
811	sga_remap(ib_thread)%start_pos+ &
812	sga_remap(ib_thread)%num_links-1) = &
813	sga_remap(ib_thread)%wts(:,1:sga_remap(ib_thread)%num_links)
814
815	end do
816
817	if (nlogprt.ge.2) then
818
819	do ib_thread = 1, il_nbthreads
820	if (sga_remap(ib_thread)%nb_resize.gt.0) then
821	write(nulou,*) ' Number of thread_resize_remap_vars on thread ',&
822	ib_thread, ' = ', sga_remap(ib_thread)%nb_resize
823	end if
824	end do
825
826	end if
827
828	end if
829
830	! Gather the complete results on master proc
831
832	if (mpi_size_map .gt. 1) then
833
834	IF (mpi_rank_map == mpi_root_map) THEN
835	ALLOCATE(ila_num_links_mpi(mpi_size_map))
836	ELSE
837	ALLOCATE(ila_num_links_mpi(1))
838	END IF
839
840	CALL MPI_Gather (num_links_map1, 1,MPI_INT,&
841	& ila_num_links_mpi,1,MPI_INT,&
842	& mpi_root_map,mpi_comm_map,il_err)
843
844	CALL MPI_Allreduce(MPI_IN_PLACE, grid2_frac, grid2_size, MPI_DOUBLE, MPI_SUM, &
845	mpi_comm_map,il_err)
846
847	IF (mpi_rank_map == mpi_root_map) THEN
848	num_links_map1 = SUM(ila_num_links_mpi)
849	if (num_links_map1 > max_links_map1) &
850	call resize_remap_vars(1,num_links_map1-max_links_map1)
851
852	ALLOCATE(ila_req_mpi(6,mpi_size_map-1))
853	ALLOCATE(ila_sta_mpi(MPI_STATUS_SIZE,6,mpi_size_map-1))
854
855	DO n = 1, mpi_size_map-1
856	buff_base = SUM(ila_num_links_mpi(1:n))+1
857	CALL MPI_IRecv(grid1_add_map1(buff_base),&
858	& ila_num_links_mpi(n+1),MPI_INT,n,1,mpi_comm_map,&
859	& ila_req_mpi(1,n),il_err)
860
861	CALL MPI_IRecv(grid2_add_map1(buff_base),&
862	& ila_num_links_mpi(n+1),MPI_INT,n,2,mpi_comm_map,&
863	& ila_req_mpi(2,n),il_err)
864
865	CALL MPI_IRecv(wts_map1(:,buff_base),&
866	& num_wts*ila_num_links_mpi(n+1),MPI_DOUBLE,n,3,mpi_comm_map,&
867	& ila_req_mpi(3,n),il_err)
868
869	CALL MPI_IRecv(grid1_area(ila_mpi_mn(n+1)),&
870	& ila_mpi_mx(n+1)-ila_mpi_mn(n+1)+1,MPI_DOUBLE,n,4,mpi_comm_map,&
871	& ila_req_mpi(4,n),il_err)
872
873	CALL MPI_IRecv(grid1_centroid_lat(ila_mpi_mn(n+1)),&
874	& ila_mpi_mx(n+1)-ila_mpi_mn(n+1)+1,MPI_DOUBLE,n,5,mpi_comm_map,&
875	& ila_req_mpi(5,n),il_err)
876
877	CALL MPI_IRecv(grid1_centroid_lon(ila_mpi_mn(n+1)),&
878	& ila_mpi_mx(n+1)-ila_mpi_mn(n+1)+1,MPI_DOUBLE,n,6,mpi_comm_map,&
879	& ila_req_mpi(6,n),il_err)
880
881	END DO
882
883	DO n=1,6
884	CALL MPI_Waitall(mpi_size_map-1,ila_req_mpi(n,:),ila_sta_mpi(1,n,1),il_err)
885	END DO
886
887	DEALLOCATE(ila_req_mpi)
888	DEALLOCATE(ila_sta_mpi)
889
890	ELSE
891
892	CALL MPI_Send(grid1_add_map1,num_links_map1,MPI_INT,&
893	& mpi_root_map,1,mpi_comm_map,il_err)
894
895	CALL MPI_Send(grid2_add_map1,num_links_map1,MPI_INT,&
896	& mpi_root_map,2,mpi_comm_map,il_err)
897
898	CALL MPI_Send(wts_map1,num_wts*num_links_map1,MPI_DOUBLE,&
899	& mpi_root_map,3,mpi_comm_map,il_err)
900
901	CALL MPI_Send(grid1_area(ila_mpi_mn(mpi_rank_map+1)),&
902	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
903	& mpi_root_map,4,mpi_comm_map,il_err)
904
905	CALL MPI_Send(grid1_centroid_lat(ila_mpi_mn(mpi_rank_map+1)),&
906	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
907	& mpi_root_map,5,mpi_comm_map,il_err)
908
909	CALL MPI_Send(grid1_centroid_lon(ila_mpi_mn(mpi_rank_map+1)),&
910	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
911	& mpi_root_map,6,mpi_comm_map,il_err)
912
913	END IF
914
915	deallocate(ila_num_links_mpi)
916
917	end if
918
919	num_links_map1_sweep1 = num_links_map1
920
921	if (nlogprt .ge. 2) &
922	write(6,*) ' num_links_map1 after sweep 1', num_links_map1_sweep1
923
924	deallocate(ila_thr_mn)
925	deallocate(ila_thr_mx)
926	if (il_nbthreads .gt. 1) then
927	do ib_thread = 1, il_nbthreads
928	deallocate(sga_remap(ib_thread)%grid1_add)
929	deallocate(sga_remap(ib_thread)%grid2_add)
930	deallocate(sga_remap(ib_thread)%wts)
931	end do
932	deallocate(sga_remap)
933	end if
934	!
935	#ifdef TREAT_OVERLAY
936	deallocate(grid2_overlap)
937	#endif
938
939	if (nlogprt .ge. 2) then
940	write(nulou,*) 'grid1 end sweep '
941	call OASIS_FLUSH_SCRIP(nulou)
942	endif
943
944	if (ll_timing) call timer_stop(3)
945
946	!-----------------------------------------------------------------------
947	!
948	! integrate around each cell on grid2
949	!
950	!-----------------------------------------------------------------------
951
952	if (ll_timing) call timer_start(4,'remap_conserv grid2 sweep')
953
954	if (nlogprt .ge. 2) then
955	write(nulou,*) 'grid2 sweep '
956	call OASIS_FLUSH_SCRIP(nulou)
957	endif
958
959	if (mpi_size_map .gt. 1) then
960
961	allocate(ila_mpi_sz(mpi_size_map))
962	il_splitsize = grid2_size
963	ila_mpi_sz(:) = floor(real(il_splitsize)/mpi_size_map)
964	ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) = ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) + 1
965
966	ila_mpi_mn(1) = 1
967	ila_mpi_mx(1) = ila_mpi_sz(1)
968
969	do ib_proc = 2, mpi_size_map
970	ila_mpi_mn(ib_proc) = sum(ila_mpi_sz(1:ib_proc-1)) + 1
971	ila_mpi_mx(ib_proc) = sum(ila_mpi_sz(1:ib_proc))
972	end do
973
974	deallocate(ila_mpi_sz)
975
976	else
977
978	ila_mpi_mn(1) = 1
979	ila_mpi_mx(1) = grid2_size
980
981	endif
982
983	!$OMP PARALLEL NUM_THREADS(il_envthreads) DEFAULT(NONE) &
984	!$OMP SHARED(num_wts) &
985	!$OMP SHARED(il_envthreads) &
986	!$OMP SHARED(grid1_corners,grid2_corners) &
987	!$OMP SHARED(grid1_bound_box,grid2_bound_box) &
988	!$OMP SHARED(grid1_bbox_per,grid2_bbox_per) &
989	!$OMP SHARED(grid1_size,bin_addr1,bin_addr2) &
990	!$OMP SHARED(est_num_neighbors,num_srch_bins) &
991	!$OMP SHARED(il_nbthreads) &
992	!$OMP SHARED(nulou,sga_remap) &
993	!$OMP SHARED(grid2_center_lat) &
994	!$OMP SHARED(grid1_center_lon) &
995	!$OMP SHARED(grid2_center_lon) &
996	!$OMP SHARED(grid1_corner_lat,grid1_corner_lon) &
997	!$OMP SHARED(grid2_corner_lat,grid2_corner_lon) &
998	!$OMP SHARED(grid1_mask) &
999	!$OMP SHARED(mpi_rank_map,mpi_root_map,ila_mpi_mn,ila_mpi_mx) &
1000	!$OMP SHARED(ila_thr_sz,ila_thr_mn,ila_thr_mx) &
1001	!$OMP SHARED(grid2_frac) &
1002	!$OMP SHARED(grid2_area) &
1003	!$OMP SHARED(grid2_centroid_lat) &
1004	!$OMP SHARED(grid2_centroid_lon) &
1005	!$OMP PRIVATE(nlogprt) &
1006	!$OMP PRIVATE(grid1_add,n,srch_add,corner) &
1007	!$OMP PRIVATE(dgbb1,dgbb2,dgbb3,dgbb4,dgbbp) &
1008	!$OMP PRIVATE(min_add,max_add) &
1009	!$OMP PRIVATE(min_cell,max_cell,grid2_add,num_srch_cells) &
1010	!$OMP PRIVATE(next_corn,beglat,beglon,endlat,endlon,lrevers) &
1011	!$OMP PRIVATE(begseg,lbegin,num_subseg) &
1012	!$OMP PRIVATE(intrsct_lat,intrsct_lon,lcoinc,weights,ll_debug) &
1013	!$OMP PRIVATE(intrsct_lat_off,intrsct_lon_off) &
1014	!$OMP PRIVATE(ib_thread,il_splitsize)
1015
1016	!$OMP SINGLE
1017
1018	il_nbthreads = 1
1019	!$ il_nbthreads = OMP_GET_NUM_THREADS ()
1020
1021	allocate(ila_thr_mn(il_nbthreads))
1022	allocate(ila_thr_mx(il_nbthreads))
1023
1024	if (il_nbthreads .gt. 1) then
1025
1026	nlogprt = 0
1027
1028	allocate(ila_thr_sz(il_nbthreads))
1029	il_splitsize = ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1
1030	ila_thr_sz(:) = floor(real(il_splitsize)/il_nbthreads)
1031	ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) = ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) + 1
1032
1033	ila_thr_mn(1) = ila_mpi_mn(mpi_rank_map+1)
1034	ila_thr_mx(1) = ila_thr_mn(1) + ila_thr_sz(1) - 1
1035
1036	do ib_thread = 2, il_nbthreads
1037	ila_thr_mn(ib_thread) = ila_mpi_mn(mpi_rank_map+1) + sum(ila_thr_sz(1:ib_thread-1))
1038	ila_thr_mx(ib_thread) = ila_thr_mn(ib_thread) + ila_thr_sz(ib_thread) - 1
1039	end do
1040
1041	allocate(sga_remap(il_nbthreads))
1042
1043	do ib_thread = 1, il_nbthreads
1044	il_splitsize = est_num_neighbors*ila_thr_sz(ib_thread)
1045	sga_remap(ib_thread)%max_links = il_splitsize
1046	sga_remap(ib_thread)%num_links = 0
1047	allocate(sga_remap(ib_thread)%grid1_add(il_splitsize))
1048	allocate(sga_remap(ib_thread)%grid2_add(il_splitsize))
1049	allocate(sga_remap(ib_thread)%wts(num_wts,il_splitsize))
1050	end do
1051
1052	deallocate(ila_thr_sz)
1053
1054	else
1055
1056	ila_thr_mn(1) = ila_mpi_mn(mpi_rank_map+1)
1057	ila_thr_mx(1) = ila_mpi_mx(mpi_rank_map+1)
1058
1059	end if
1060	!$OMP END SINGLE
1061
1062	!$OMP DO SCHEDULE(STATIC,1)
1063	thread_loop2: do ib_thread = 1, il_nbthreads
1064
1065	ll_debug = .false.
1066
1067	grid_loop2: do grid2_add = ila_thr_mn(ib_thread), ila_thr_mx(ib_thread)
1068
1069	if (ll_debug) then
1070	write(nulou,*) 'grid2_add=', grid2_add
1071	call OASIS_FLUSH_SCRIP(nulou)
1072	endif
1073
1074	!***
1075	!*** destination grid bounding box (preloaded)
1076	!***
1077
1078	dgbb1 = grid2_bound_box(1,grid2_add)
1079	dgbb2 = grid2_bound_box(2,grid2_add)
1080	dgbb3 = grid2_bound_box(3,grid2_add)
1081	dgbb4 = grid2_bound_box(4,grid2_add)
1082	dgbbp = grid2_bbox_per(grid2_add)
1083
1084	!***
1085	!*** restrict searches first using search bins
1086	!***
1087
1088	min_add = grid1_size
1089	max_add = 1
1090	do n=1,num_srch_bins
1091	if (grid2_add >= bin_addr2(1,n) .and. grid2_add <= bin_addr2(2,n)) then
1092	min_add = min(min_add, bin_addr1(1,n))
1093	max_add = max(max_add, bin_addr1(2,n))
1094	endif
1095	end do
1096
1097	!***
1098	!*** select the search cells
1099	!***
1100
1101	num_srch_cells = 0
1102	min_cell = max_add
1103	max_cell = min_add
1104
1105	gather2: do grid1_add = min_add, max_add
1106
1107	!***
1108	!*** further restrict searches using bounding boxes
1109	!***
1110
1111	if ( grid1_bound_box(1,grid1_add) <= dgbb2 ) then
1112	if ( grid1_bound_box(2,grid1_add) >= dgbb1 ) then
1113	if ( lf_loncross(dgbbp,dgbb3,dgbb4,grid1_add,&
1114	grid1_bbox_per,grid1_bound_box) ) then
1115
1116	num_srch_cells = num_srch_cells+1
1117	min_cell = min(min_cell, grid1_add)
1118	max_cell = max(max_cell, grid1_add)
1119
1120	end if
1121	end if
1122	end if
1123
1124	end do gather2
1125
1126	!***
1127	!*** store the search cells array
1128	!***
1129
1130	if (num_srch_cells .ne. 0) then
1131
1132	allocate(srch_add(num_srch_cells))
1133
1134	n=0
1135	pick2: do grid1_add = min_cell,max_cell
1136
1137	if ( grid1_bound_box(1,grid1_add) <= dgbb2 ) then
1138	if ( grid1_bound_box(2,grid1_add) >= dgbb1 ) then
1139	if ( lf_loncross(dgbbp,dgbb3,dgbb4,grid1_add, &
1140	grid1_bbox_per,grid1_bound_box) ) then
1141
1142	n = n+1
1143	srch_add(n) = grid1_add
1144
1145	end if
1146	end if
1147	end if
1148
1149	end do pick2
1150
1151	end if
1152
1153	if (ll_debug) then
1154	write(nulou,*)' '
1155	write(nulou,)' * Grid2 cell and associated search cells **'
1156	write(nulou,1130) grid2_add, grid2_center_lon(grid2_add), grid2_center_lat(grid2_add)
1157	do corner = 1,grid2_corners
1158	write(nulou,1131) corner, grid2_corner_lon(corner,grid2_add), grid2_corner_lat(corner,grid2_add)
1159	enddo
1160	write(nulou,1804) grid2_add, grid2_bound_box(1,grid2_add), grid2_bound_box(2,grid2_add)
1161	write(nulou,1408) grid2_add, grid2_bound_box(3,grid2_add), grid2_bound_box(4,grid2_add)
1162	write(nulou,*) ' num_srch_cells=', num_srch_cells
1163	write(nulou,*) ' '
1164	if (num_srch_cells .ne. 0) then
1165	write(nulou,*) ' srch_add(:)=', srch_add(1:num_srch_cells)
1166	write(nulou,*) ' msk_srch_add(:)=', (.true.,n=1,num_srch_cells)
1167	do n=1, num_srch_cells
1168	do corner = 1,grid2_corners
1169	write(nulou,1132) n, grid1_corner_lon(corner,srch_add(n)), grid1_corner_lat(corner,srch_add(n))
1170	end do
1171	end do
1172	end if
1173	write(nulou,)' **************************************'
1174	write(nulou,*)' '
1175	endif
1176	1130 format (' grid2 center, lon, lat = ', I8, 2X, F12.4, 2X, F12.4)
1177	1131 format (' grid2 corner, lon, lat = ', I8, 2X, F12.4, 2X, F12.4)
1178	1132 format (' srch cell, lon, lat = ', I8, 2X, F12.4, 2X, F12.4)
1179	1804 format (' grid2 index, bbox lat = ', I8, 2X, F12.4, 2X, F12.4)
1180	1408 format (' grid2 index, bbox lon = ', I8, 2X, F12.4, 2X, F12.4)
1181
1182	if (num_srch_cells .ne. 0) then
1183
1184
1185	!***
1186	!*** integrate around this cell
1187	!***
1188
1189	! full = .false.
1190	! do grid1_add = min_add,max_add
1191	! if (grid1_mask(grid1_add)) full = .true.
1192	! end do
1193	! if (full) then
1194
1195	do corner = 1,grid2_corners
1196	next_corn = mod(corner,grid2_corners) + 1
1197
1198	beglat = grid2_corner_lat(corner,grid2_add)
1199	beglon = grid2_corner_lon(corner,grid2_add)
1200	endlat = grid2_corner_lat(next_corn,grid2_add)
1201	endlon = grid2_corner_lon(next_corn,grid2_add)
1202	lrevers = .false.
1203
1204	!***
1205	!*** to ensure exact path taken during both
1206	!*** sweeps, always integrate in the same direction
1207	!***
1208
1209	if ((endlat < beglat) .or. &
1210	(endlat == beglat .and. endlon < beglon)) then
1211	beglat = grid2_corner_lat(next_corn,grid2_add)
1212	beglon = grid2_corner_lon(next_corn,grid2_add)
1213	endlat = grid2_corner_lat(corner,grid2_add)
1214	endlon = grid2_corner_lon(corner,grid2_add)
1215	lrevers = .true.
1216	if (ll_debug) &
1217	write(nulou, *) ' sweep2 LREVERS TRUE'
1218	endif
1219	begseg(1) = beglat
1220	begseg(2) = beglon
1221	lbegin = .true.
1222
1223	!***
1224	!*** if this is a constant-longitude segment, skip the rest
1225	!*** since the line integral contribution will be zero.
1226	!***
1227
1228	if (ll_debug) then
1229	if (endlon .eq. beglon) then
1230	write(nulou,1113) beglon, beglat
1231	write(nulou,1114) endlon, endlat
1232	write(nulou, *)' sweep2 endlon == beglon; skip segment'
1233	write(nulou,*) ' '
1234	endif
1235	endif
1236	if (endlon /= beglon) then
1237	num_subseg = 0
1238
1239	!***
1240	!*** integrate along this segment, detecting intersections
1241	!*** and computing the line integral for each sub-segment
1242	!***
1243
1244	do while (beglat /= endlat .or. beglon /= endlon)
1245
1246	!***
1247	!*** prevent infinite loops if integration gets stuck
1248	!*** near cell or threshold boundary
1249	!***
1250
1251	num_subseg = num_subseg + 1
1252	if (num_subseg > max_subseg) then
1253	write(nulou,*) 'ERROR=>integration stalled:'
1254	write(nulou,*) 'num_subseg exceeded limit'
1255	write(nulou,*) 'Verify corners in grids.nc, especially'
1256	write(nulou,*) 'if calculated by OASIS routine corners'
1257	write(nulou,*) 'integration stalled: num_subseg exceeded limit'
1258	call OASIS_FLUSH_SCRIP(nulou)
1259	stop
1260	endif
1261
1262	!***
1263	!*** find next intersection of this segment with a line
1264	!*** on grid 1.
1265	!***
1266
1267	if (ll_debug) then
1268	write(nulou,1115) beglon, beglat
1269	write(nulou,1116) endlon, endlat
1270	write(nulou,*) ' '
1271	endif
1272
1273	call intersection(grid1_add,intrsct_lat,intrsct_lon,lcoinc, &
1274	beglat, beglon, endlat, endlon, begseg, &
1275	lbegin, lrevers, &
1276	grid1_corners, num_srch_cells, &
1277	srch_add, grid1_corner_lon, grid1_corner_lat, &
1278	intrsct_lat_off,intrsct_lon_off)
1279
1280	if (ll_debug) then
1281	write(nulou,*) ' After call intersection, grid1_add', grid1_add
1282	write(nulou,1117) beglon, beglat
1283	write(nulou,1118) endlon, endlat
1284	write(nulou,1119) intrsct_lon, intrsct_lat
1285	write(nulou,*) ' '
1286	endif
1287
1288	lbegin = .false.
1289
1290	!***
1291	!*** compute line integral for this subsegment.
1292	!***
1293
1294	if (grid1_add /= 0) then
1295	call line_integral(weights, num_wts, &
1296	beglon, intrsct_lon, beglat, intrsct_lat,&
1297	grid1_center_lon(grid1_add), &
1298	grid2_center_lon(grid2_add))
1299
1300	if (ll_debug) then
1301	write(nulou,*) ' A2) WEIGHTS for this subsegment =', weights(1)
1302	write(nulou,*) ' '
1303	endif
1304	else
1305	call line_integral(weights, num_wts, &
1306	beglon, intrsct_lon, beglat, intrsct_lat,&
1307	grid2_center_lon(grid2_add), &
1308	grid2_center_lon(grid2_add))
1309
1310	if (ll_debug) then
1311	write(nulou,*) ' B2) WEIGHTS for this subsegment =', weights(1)
1312	write(nulou,*) ' '
1313	endif
1314	endif
1315
1316	if (lrevers) then
1317	weights = -weights
1318	if (ll_debug) then
1319	write(nulou,*) ' LREVERS; WEIGHTS for this subsegment =', weights(1)
1320	write(nulou,*) ' '
1321	endif
1322	endif
1323	!***
1324	!*** store the appropriate addresses and weights.
1325	!*** also add contributions to cell areas and centroids.
1326	!*** if there is a coincidence, do not store weights
1327	!*** because they have been captured in the previous loop.
1328	!*** the grid1 mask is the master mask
1329	!***
1330
1331	if (ll_debug .and. lcoinc) &
1332	write(nulou,*) ' LCOINC is TRUE; weight not stored'
1333
1334	if (.not. lcoinc .and. grid1_add /= 0) then
1335	if (grid1_mask(grid1_add)) then
1336	call store_link_cnsrv(grid1_add, grid2_add, weights, ib_thread)
1337
1338	if (ll_debug) then
1339	write(nulou,*) ' after store_link_cnsrv norm2'
1340	write(nulou,1120) grid1_add, grid2_add,beglon, beglat, intrsct_lon,intrsct_lat,weights(1)
1341	endif
1342	!EM grid1_frac(grid1_add) = grid1_frac(grid1_add) + weights(1)
1343	grid2_frac(grid2_add) = grid2_frac(grid2_add) + weights(num_wts+1)
1344	endif
1345
1346	endif
1347
1348	grid2_area(grid2_add) = grid2_area(grid2_add) + weights(num_wts+1)
1349	grid2_centroid_lat(grid2_add) = grid2_centroid_lat(grid2_add) + weights(num_wts+2)
1350	grid2_centroid_lon(grid2_add) = grid2_centroid_lon(grid2_add) + weights(num_wts+3)
1351
1352	!***
1353	!*** reset beglat and beglon for next subsegment.
1354	!***
1355
1356	beglat = intrsct_lat
1357	beglon = intrsct_lon
1358
1359	end do
1360
1361	end if
1362
1363	!***
1364	!*** end of segment
1365	!***
1366
1367	end do
1368
1369	!***
1370	!*** finished with this cell: deallocate search array and
1371	!*** start on next cell
1372
1373	deallocate(srch_add)
1374
1375	end if ! num_srch_cells .NE. 0
1376
1377	end do grid_loop2
1378
1379	end do thread_loop2
1380	!$OMP END DO
1381
1382	!$OMP END PARALLEL
1383
1384	if (il_nbthreads .gt. 1) then
1385
1386	sga_remap(1)%start_pos = num_links_map1_sweep1 + 1
1387
1388	il_splitsize = sga_remap(1)%num_links
1389	do ib_thread = 2, il_nbthreads
1390	il_splitsize = il_splitsize + sga_remap(ib_thread)%num_links
1391	sga_remap(ib_thread)%start_pos = sga_remap(ib_thread-1)%start_pos + &
1392	sga_remap(ib_thread-1)%num_links
1393	end do
1394
1395	num_links_map1 = num_links_map1_sweep1 + il_splitsize
1396	if (num_links_map1 > max_links_map1) &
1397	call resize_remap_vars(1,num_links_map1-max_links_map1)
1398
1399	do ib_thread = 1, il_nbthreads
1400	grid1_add_map1(sga_remap(ib_thread)%start_pos: &
1401	sga_remap(ib_thread)%start_pos+ &
1402	sga_remap(ib_thread)%num_links-1) = &
1403	sga_remap(ib_thread)%grid1_add(1:sga_remap(ib_thread)%num_links)
1404	grid2_add_map1(sga_remap(ib_thread)%start_pos: &
1405	sga_remap(ib_thread)%start_pos+ &
1406	sga_remap(ib_thread)%num_links-1) = &
1407	sga_remap(ib_thread)%grid2_add(1:sga_remap(ib_thread)%num_links)
1408	wts_map1 (:,sga_remap(ib_thread)%start_pos: &
1409	sga_remap(ib_thread)%start_pos+ &
1410	sga_remap(ib_thread)%num_links-1) = &
1411	sga_remap(ib_thread)%wts(:,1:sga_remap(ib_thread)%num_links)
1412
1413	end do
1414
1415	if (nlogprt.ge.2) then
1416
1417	do ib_thread = 1, il_nbthreads
1418	if (sga_remap(ib_thread)%nb_resize.gt.0) then
1419	write(nulou,*) ' Number of thread_resize_remap_vars on thread ',&
1420	ib_thread, ' = ', sga_remap(ib_thread)%nb_resize
1421	end if
1422	end do
1423
1424	end if
1425
1426	end if
1427
1428	! Gather the complete results on master proc
1429
1430	if (mpi_size_map .gt. 1) then
1431
1432	IF (mpi_rank_map == mpi_root_map) THEN
1433	ALLOCATE(ila_num_links_mpi(mpi_size_map))
1434	ELSE
1435	ALLOCATE(ila_num_links_mpi(1))
1436	END IF
1437
1438	CALL MPI_Gather (num_links_map1-num_links_map1_sweep1, 1,MPI_INT,&
1439	& ila_num_links_mpi,1,MPI_INT,&
1440	& mpi_root_map,mpi_comm_map,il_err)
1441
1442	IF (mpi_rank_map == mpi_root_map) THEN
1443	num_links_map1 = num_links_map1_sweep1 + SUM(ila_num_links_mpi)
1444	if (num_links_map1 > max_links_map1) &
1445	call resize_remap_vars(1,num_links_map1-max_links_map1)
1446
1447	ALLOCATE(ila_req_mpi(7,mpi_size_map-1))
1448	ALLOCATE(ila_sta_mpi(MPI_STATUS_SIZE,7,mpi_size_map-1))
1449
1450	DO n = 1, mpi_size_map-1
1451	buff_base = num_links_map1_sweep1+SUM(ila_num_links_mpi(1:n))+1
1452
1453	CALL MPI_IRecv(grid1_add_map1(buff_base),&
1454	& ila_num_links_mpi(n+1),MPI_INT,n,1,mpi_comm_map,&
1455	& ila_req_mpi(1,n),il_err)
1456
1457	CALL MPI_IRecv(grid2_add_map1(buff_base),&
1458	& ila_num_links_mpi(n+1),MPI_INT,n,2,mpi_comm_map,&
1459	& ila_req_mpi(2,n),il_err)
1460
1461	CALL MPI_IRecv(wts_map1(:,buff_base),&
1462	& num_wts*ila_num_links_mpi(n+1),MPI_DOUBLE,n,3,mpi_comm_map,&
1463	& ila_req_mpi(3,n),il_err)
1464
1465	CALL MPI_IRecv(grid2_area(ila_mpi_mn(n+1)),&
1466	& ila_mpi_mx(n+1)-ila_mpi_mn(n+1)+1,MPI_DOUBLE,n,4,mpi_comm_map,&
1467	& ila_req_mpi(4,n),il_err)
1468
1469	CALL MPI_IRecv(grid2_centroid_lat(ila_mpi_mn(n+1)),&
1470	& ila_mpi_mx(n+1)-ila_mpi_mn(n+1)+1,MPI_DOUBLE,n,5,mpi_comm_map,&
1471	& ila_req_mpi(5,n),il_err)
1472
1473	CALL MPI_IRecv(grid2_centroid_lon(ila_mpi_mn(n+1)),&
1474	& ila_mpi_mx(n+1)-ila_mpi_mn(n+1)+1,MPI_DOUBLE,n,6,mpi_comm_map,&
1475	& ila_req_mpi(6,n),il_err)
1476
1477	CALL MPI_IRecv(grid2_frac(ila_mpi_mn(n+1)),&
1478	& ila_mpi_mx(n+1)-ila_mpi_mn(n+1)+1,MPI_DOUBLE,n,7,mpi_comm_map,&
1479	& ila_req_mpi(7,n),il_err)
1480
1481	END DO
1482
1483	DO n=1,7
1484	CALL MPI_Waitall(mpi_size_map-1,ila_req_mpi(n,:),ila_sta_mpi(1,n,1),il_err)
1485	END DO
1486
1487	DEALLOCATE(ila_req_mpi)
1488	DEALLOCATE(ila_sta_mpi)
1489
1490	ELSE
1491
1492	CALL MPI_Send(grid1_add_map1(num_links_map1_sweep1+1),num_links_map1-num_links_map1_sweep1,MPI_INT,&
1493	& mpi_root_map,1,mpi_comm_map,il_err)
1494
1495	CALL MPI_Send(grid2_add_map1(num_links_map1_sweep1+1),num_links_map1-num_links_map1_sweep1,MPI_INT,&
1496	& mpi_root_map,2,mpi_comm_map,il_err)
1497
1498	CALL MPI_Send(wts_map1(1,num_links_map1_sweep1+1),num_wts*(num_links_map1-num_links_map1_sweep1),MPI_DOUBLE,&
1499	& mpi_root_map,3,mpi_comm_map,il_err)
1500
1501	CALL MPI_Send(grid2_area(ila_mpi_mn(mpi_rank_map+1)),&
1502	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
1503	& mpi_root_map,4,mpi_comm_map,il_err)
1504
1505	CALL MPI_Send(grid2_centroid_lat(ila_mpi_mn(mpi_rank_map+1)),&
1506	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
1507	& mpi_root_map,5,mpi_comm_map,il_err)
1508
1509	CALL MPI_Send(grid2_centroid_lon(ila_mpi_mn(mpi_rank_map+1)),&
1510	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
1511	& mpi_root_map,6,mpi_comm_map,il_err)
1512
1513	CALL MPI_Send(grid2_frac(ila_mpi_mn(mpi_rank_map+1)),&
1514	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
1515	& mpi_root_map,7,mpi_comm_map,il_err)
1516
1517	END IF
1518
1519	deallocate(ila_num_links_mpi)
1520
1521	end if
1522
1523	if (nlogprt .ge. 2) &
1524	write(6,*) ' num_links_map1 after sweep 2', num_links_map1
1525
1526	deallocate(ila_thr_mn)
1527	deallocate(ila_thr_mx)
1528
1529	deallocate(ila_mpi_mn, ila_mpi_mx)
1530
1531	if (il_nbthreads .gt. 1) then
1532	do ib_thread = 1, il_nbthreads
1533	deallocate(sga_remap(ib_thread)%grid1_add)
1534	deallocate(sga_remap(ib_thread)%grid2_add)
1535	deallocate(sga_remap(ib_thread)%wts)
1536	end do
1537	deallocate(sga_remap)
1538	end if
1539
1540	if (nlogprt .ge. 2) then
1541	write(nulou,*) 'grid2 end sweep '
1542	call OASIS_FLUSH_SCRIP(nulou)
1543	endif
1544
1545	ll_debug = .false.
1546
1547	if (ll_debug) then
1548	do n=1,num_links_map1
1549	write(nulou,*) 'grid1, grid2, weight= ', grid1_add_map1(n), grid2_add_map1(n), wts_map1(1,n)
1550	end do
1551	endif
1552
1553	if (ll_timing) call timer_stop(4)
1554
1555	!-----------------------------------------------------------------------
1556	! WARNING
1557	! After this line, calculations MUST be single threaded
1558	!-----------------------------------------------------------------------
1559	il_nbthreads = 1
1560
1561	!-----------------------------------------------------------------------
1562	!
1563	! correct for situations where N/S pole not explicitly included in
1564	! grid (i.e. as a grid corner point). if pole is missing from only
1565	! one grid, need to correct only the area and centroid of that
1566	! grid. if missing from both, do complete weight calculation.
1567	!
1568	!-----------------------------------------------------------------------
1569
1570	!*** North Pole
1571	weights(1) = pi2
1572	weights(2) = pi*pi
1573	weights(3) = zero
1574	weights(4) = pi2
1575	weights(5) = pi*pi
1576	weights(6) = zero
1577
1578	grid1_np = -1
1579	pole_loop1: do n=1,grid1_size
1580	if (grid1_area(n) < -three*pih .and. grid1_center_lat(n) > zero) then
1581	grid1_np = n
1582	exit pole_loop1
1583	endif
1584	end do pole_loop1
1585
1586	grid2_np = -1
1587	pole_loop2: do n=1,grid2_size
1588	if (grid2_area(n) < -three*pih .and. grid2_center_lat(n) > zero) then
1589	grid2_np = n
1590	exit pole_loop2
1591	endif
1592	end do pole_loop2
1593
1594	if (grid1_np > 0) then
1595	grid1_area(grid1_np) = grid1_area(grid1_np) + weights(1)
1596	grid1_centroid_lat(grid1_np) = grid1_centroid_lat(grid1_np) + weights(2)
1597	grid1_centroid_lon(grid1_np) = grid1_centroid_lon(grid1_np) + weights(3)
1598	endif
1599
1600	if (grid2_np > 0) then
1601	grid2_area(grid2_np) = grid2_area(grid2_np) + weights(num_wts+1)
1602	grid2_centroid_lat(grid2_np) = grid2_centroid_lat(grid2_np) + weights(num_wts+2)
1603	grid2_centroid_lon(grid2_np) = grid2_centroid_lon(grid2_np) + weights(num_wts+3)
1604	endif
1605
1606	if (grid1_np > 0 .and. grid2_np > 0) then
1607	call store_link_cnsrv(grid1_np, grid2_np, weights, 1)
1608	!EM grid1_frac(grid1_np) = grid1_frac(grid1_np) + weights(1)
1609	grid2_frac(grid2_np) = grid2_frac(grid2_np) + weights(num_wts+1)
1610	endif
1611
1612	!*** South Pole
1613	weights(1) = pi2
1614	weights(2) = -pi*pi
1615	weights(3) = zero
1616	weights(4) = pi2
1617	weights(5) = -pi*pi
1618	weights(6) = zero
1619
1620	grid1_sp = -1
1621	pole_loop3: do n=1,grid1_size
1622	if (grid1_area(n) < -three*pih .and. grid1_center_lat(n) < zero) then
1623	grid1_sp = n
1624	exit pole_loop3
1625	endif
1626	end do pole_loop3
1627
1628	grid2_sp = -1
1629	pole_loop4: do n=1,grid2_size
1630	if (grid2_area(n) < -three*pih .and. grid2_center_lat(n) < zero) then
1631	grid2_sp = n
1632	exit pole_loop4
1633	endif
1634	end do pole_loop4
1635
1636	if (grid1_sp > 0) then
1637	grid1_area(grid1_sp) = grid1_area(grid1_sp) + weights(1)
1638	grid1_centroid_lat(grid1_sp) = grid1_centroid_lat(grid1_sp) + weights(2)
1639	grid1_centroid_lon(grid1_sp) = grid1_centroid_lon(grid1_sp) + weights(3)
1640	endif
1641
1642	if (grid2_sp > 0) then
1643	grid2_area(grid2_sp) = grid2_area(grid2_sp) + weights(num_wts+1)
1644	grid2_centroid_lat(grid2_sp) = grid2_centroid_lat(grid2_sp) + weights(num_wts+2)
1645	grid2_centroid_lon(grid2_sp) = grid2_centroid_lon(grid2_sp) + weights(num_wts+3)
1646	endif
1647
1648	if (grid1_sp > 0 .and. grid2_sp > 0) then
1649	call store_link_cnsrv(grid1_sp, grid2_sp, weights, 1)
1650
1651	!EM grid1_frac(grid1_sp) = grid1_frac(grid1_sp) + weights(1)
1652	grid2_frac(grid2_sp) = grid2_frac(grid2_sp) + weights(num_wts+1)
1653	endif
1654
1655	!-----------------------------------------------------------------------
1656	!
1657	! finish centroid computation
1658	!
1659	!-----------------------------------------------------------------------
1660
1661	where (grid1_area /= zero)
1662	grid1_centroid_lat = grid1_centroid_lat/grid1_area
1663	grid1_centroid_lon = grid1_centroid_lon/grid1_area
1664	end where
1665
1666	where (grid2_area /= zero)
1667	grid2_centroid_lat = grid2_centroid_lat/grid2_area
1668	grid2_centroid_lon = grid2_centroid_lon/grid2_area
1669	end where
1670
1671	!-----------------------------------------------------------------------
1672	!
1673	! include centroids in weights and normalize using destination
1674	! area if requested
1675	!
1676	!-----------------------------------------------------------------------
1677
1678
1679	do n=1,num_links_map1
1680
1681	grid1_add = grid1_add_map1(n)
1682	grid2_add = grid2_add_map1(n)
1683	do nwgt=1,num_wts
1684	weights( nwgt) = wts_map1(nwgt,n)
1685	! if (num_maps > 1) then
1686	! weights(num_wts+nwgt) = wts_map2(nwgt,n)
1687	! endif
1688	end do
1689
1690	select case(norm_opt)
1691	case (norm_opt_dstarea)
1692	if (grid2_area(grid2_add) /= zero) then
1693	if (luse_grid2_area) then
1694	if (.not. lstore_grid2_area) then
1695	write(nulou,*) 'ERROR: remap_conserv norm_opt_dstarea failed with missing grid2_area_in'
1696	call OASIS_FLUSH_SCRIP(nulou)
1697	stop
1698	endif
1699	norm_factor = one/grid2_area_in(grid2_add)
1700	else
1701	norm_factor = one/grid2_area(grid2_add)
1702	endif
1703	else
1704	norm_factor = zero
1705	endif
1706	case (norm_opt_frcarea)
1707	if (grid2_frac(grid2_add) /= zero) then
1708	if (luse_grid2_area) then
1709	if (.not. lstore_grid2_area) then
1710	write(nulou,*) 'ERROR: remap_conserv norm_opt_frcarea failed with missing grid2_area_in'
1711	call OASIS_FLUSH_SCRIP(nulou)
1712	stop
1713	endif
1714	norm_factor = grid2_area(grid2_add)/ &
1715	(grid2_frac(grid2_add)*grid2_area_in(grid2_add))
1716	else
1717	norm_factor = one/grid2_frac(grid2_add)
1718	endif
1719	else
1720	norm_factor = zero
1721	endif
1722	case (norm_opt_dstartr)
1723	if (grid2_area(grid2_add) /= zero) then
1724	norm_factor = one/grid2_area(grid2_add)
1725	if (grid1_add .ne. grid1_np .and. grid1_add .ne. grid1_sp &
1726	.and. grid1_area(grid1_add) /= zero) then
1727	if (.not. lstore_grid1_area) then
1728	write(nulou,*) 'ERROR: remap_conserv norm_opt_dstartr failed with missing grid1_area_in'
1729	call OASIS_FLUSH_SCRIP(nulou)
1730	stop
1731	endif
1732	norm_correc = grid1_area_in(grid1_add)/grid1_area(grid1_add)
1733	if (norm_correc.gt.0.8.and.norm_correc.lt.1.2) &
1734	norm_factor = norm_factor * norm_correc
1735	end if
1736	if (grid2_add .ne. grid2_np .and. grid2_add .ne. grid2_sp) then
1737	if (.not. lstore_grid2_area) then
1738	write(nulou,*) 'ERROR: remap_conserv norm_opt_dstartr failed with missing grid2_area_in'
1739	call OASIS_FLUSH_SCRIP(nulou)
1740	stop
1741	endif
1742	norm_correc = grid2_area(grid2_add)/grid2_area_in(grid2_add)
1743	if (norm_correc.gt.0.8.and.norm_correc.lt.1.2) &
1744	norm_factor = norm_factor * norm_correc
1745	end if
1746	else
1747	norm_factor = zero
1748	endif
1749	case (norm_opt_frcartr)
1750	if (grid2_frac(grid2_add) /= zero) then
1751	norm_factor = one/grid2_frac(grid2_add)
1752	if (grid1_add .ne. grid1_np .and. grid1_add .ne. grid1_sp &
1753	.and. grid1_area(grid1_add) /= zero) then
1754	if (.not. lstore_grid1_area) then
1755	write(nulou,*) 'ERROR: remap_conserv norm_opt_frcartr failed with missing grid1_area_in'
1756	call OASIS_FLUSH_SCRIP(nulou)
1757	stop
1758	endif
1759	norm_correc = grid1_area_in(grid1_add)/grid1_area(grid1_add)
1760	if (norm_correc.gt.0.8.and.norm_correc.lt.1.2) &
1761	norm_factor = norm_factor * norm_correc
1762	end if
1763	if (grid2_add .ne. grid2_np .and. grid2_add .ne. grid2_sp) then
1764	if (.not. lstore_grid2_area) then
1765	write(nulou,*) 'ERROR: remap_conserv norm_opt_frcartr failed with missing grid2_area_in'
1766	call OASIS_FLUSH_SCRIP(nulou)
1767	stop
1768	endif
1769	norm_correc = grid2_area(grid2_add)/grid2_area_in(grid2_add)
1770	if (norm_correc.gt.0.8.and.norm_correc.lt.1.2) &
1771	norm_factor = norm_factor * norm_correc
1772	end if
1773	else
1774	norm_factor = zero
1775	endif
1776	case (norm_opt_none)
1777	norm_factor = one
1778	end select
1779
1780	wts_map1(1,n) = weights(1)*norm_factor
1781	wts_map1(2,n) = (weights(2) - weights(1)grid1_centroid_lat(grid1_add))norm_factor
1782	wts_map1(3,n) = (weights(3) - weights(1)grid1_centroid_lon(grid1_add))norm_factor
1783
1784	! if (num_maps > 1) then
1785	! select case(norm_opt)
1786	! case (norm_opt_dstarea)
1787	! if (grid1_area(grid1_add) /= zero) then
1788	! if (luse_grid1_area) then
1789	! norm_factor = one/grid1_area_in(grid1_add)
1790	! else
1791	! norm_factor = one/grid1_area(grid1_add)
1792	! endif
1793	! else
1794	! norm_factor = zero
1795	! endif
1796	! case (norm_opt_frcarea)
1797	! if (grid1_frac(grid1_add) /= zero) then
1798	! if (luse_grid1_area) then
1799	! norm_factor = grid1_area(grid1_add)/
1800	! & (grid1_frac(grid1_add)*
1801	! & grid1_area_in(grid1_add))
1802	! else
1803	! norm_factor = one/grid1_frac(grid1_add)
1804	! endif
1805	! else
1806	! norm_factor = zero
1807	! endif
1808	! case (norm_opt_none)
1809	! norm_factor = one
1810	! end select
1811	!
1812	! wts_map2(1,n) = weights(num_wts+1)*norm_factor
1813	! wts_map2(2,n) = (weights(num_wts+2) - weights(num_wts+1)*
1814	! & grid2_centroid_lat(grid2_add))*
1815	! & norm_factor
1816	! wts_map2(3,n) = (weights(num_wts+3) - weights(num_wts+1)*
1817	! & grid2_centroid_lon(grid2_add))*
1818	! & norm_factor
1819	! endif
1820
1821	end do
1822
1823	if (nlogprt .ge. 2) then
1824	write(nulou,*) 'Total number of links = ',num_links_map1
1825	call OASIS_FLUSH_SCRIP(nulou)
1826	endif
1827
1828	!EM where (grid1_area /= zero) grid1_frac = grid1_frac/grid1_area
1829	where (grid2_area /= zero) grid2_frac = grid2_frac/grid2_area
1830
1831	!-----------------------------------------------------------------------
1832	!
1833	! perform some error checking on final weights
1834	!
1835	!-----------------------------------------------------------------------
1836
1837	grid2_centroid_lat = zero
1838	grid2_centroid_lon = zero
1839
1840	do n=1,grid1_size
1841	if (nlogprt .ge. 3) then
1842	if (grid1_area(n) < -.01) then
1843	write(nulou,*) 'Grid 1 area error: n, area, mask =',n,grid1_area(n), grid1_mask(n)
1844	endif
1845	if (grid1_centroid_lat(n) < -pih-.01 .or.grid1_centroid_lat(n) > pih+.01) then
1846	write(nulou,*)'Grid1 centroid lat error: n, centroid_lat, mask=' &
1847	,n,grid1_centroid_lat(n), grid1_mask(n)
1848	endif
1849	call OASIS_FLUSH_SCRIP(nulou)
1850	endif
1851	grid1_centroid_lat(n) = zero
1852	grid1_centroid_lon(n) = zero
1853	end do
1854
1855	do n=1,grid2_size
1856	if (nlogprt .ge. 3) then
1857	if (grid2_area(n) < -.01) then
1858	write(nulou,*) 'Grid 2 area error: n, area, mask =' ,n,grid2_area(n), grid2_mask(n)
1859	endif
1860	if (grid2_centroid_lat(n) < -pih-.01 .or. grid2_centroid_lat(n) > pih+.01) then
1861	write(nulou,*) 'Grid 2 centroid lat error: n, centroid_lat, mask=' &
1862	,n,grid2_centroid_lat(n), grid2_mask(n)
1863	endif
1864	call OASIS_FLUSH_SCRIP(nulou)
1865	endif
1866	grid2_centroid_lat(n) = zero
1867	grid2_centroid_lon(n) = zero
1868	end do
1869
1870	do n=1,num_links_map1
1871	grid2_add = grid2_add_map1(n)
1872	grid2_centroid_lat(grid2_add) = grid2_centroid_lat(grid2_add) + wts_map1(1,n)
1873
1874	! if (num_maps > 1) then
1875	! if (wts_map2(1,n) < -.01) then
1876	! print *,'Map 2 weight < 0 '
1877	! PRINT *,
1878	! & 'grid1_add,grid2_add, wts_map2, grid1_mask, grid2_mask',
1879	! & grid1_add, grid2_add, wts_map2(1,n),
1880	! & grid1_mask(grid1_add), grid2_mask(grid2_add)
1881	! endif
1882	! if (norm_opt /= norm_opt_none .and. wts_map2(1,n) > 1.01) then
1883	! print *,'Map 2 weight < 0 '
1884	! PRINT *,
1885	! & 'grid1_add,grid2_add,wts_map2, grid1_mask,grid2_mask',
1886	! & grid1_add, grid2_add, wts_map2(1,n),
1887	! & grid1_mask(grid1_add), grid2_mask(grid2_add)
1888	! endif
1889	! grid1_centroid_lat(grid1_add) =
1890	! & grid1_centroid_lat(grid1_add) + wts_map2(1,n)
1891	! endif
1892	end do
1893
1894	do n=1,grid2_size
1895	select case(norm_opt)
1896	case (norm_opt_dstarea)
1897	norm_factor = grid2_frac(n)
1898	case (norm_opt_frcarea)
1899	norm_factor = one
1900	case (norm_opt_dstartr)
1901	norm_factor = grid2_frac(n)
1902	case (norm_opt_frcartr)
1903	norm_factor = one
1904	case (norm_opt_none)
1905	if (luse_grid2_area) then
1906	if (.not. lstore_grid2_area) then
1907	write(nulou,*) 'ERROR: remap_conserv norm_opt_none failed with missing grid2_area_in'
1908	call OASIS_FLUSH_SCRIP(nulou)
1909	stop
1910	endif
1911	norm_factor = grid2_area_in(n)
1912	else
1913	norm_factor = grid2_area(n)
1914	endif
1915	end select
1916	if (nlogprt .ge. 3) then
1917	if (abs(grid2_centroid_lat(n)-norm_factor) > .01) then
1918	print *, 'Error sum wts map1:grid2_add,grid2_centroid_lat,norm_factor=', &
1919	n,grid2_centroid_lat(n),norm_factor,grid2_mask(n)
1920	endif
1921	endif
1922	end do
1923
1924	! if (num_maps > 1) then
1925	! do n=1,grid1_size
1926	! select case(norm_opt)
1927	! case (norm_opt_dstarea)
1928	! norm_factor = grid1_frac(n)
1929	! case (norm_opt_frcarea)
1930	! norm_factor = one
1931	! case (norm_opt_none)
1932	! if (luse_grid1_area) then
1933	! norm_factor = grid1_area_in(n)
1934	! else
1935	! norm_factor = grid1_area(n)
1936	! endif
1937	! end select
1938	! if (abs(grid1_centroid_lat(n)-norm_factor) > .01) then
1939	! print *,
1940	! &'Error sum wts map2:grid1_add,grid1_centroid_lat,norm_factor='
1941	! &,n,grid1_centroid_lat(n),
1942	! &norm_factor,grid1_mask(n)
1943	! endif
1944	! end do
1945	! endif
1946
1947	!-----------------------------------------------------------------------
1948	!
1949	! deallocate allocated arrays
1950	!
1951	!-----------------------------------------------------------------------
1952
1953	deallocate (grid1_centroid_lat, grid1_centroid_lon, grid2_centroid_lat, grid2_centroid_lon)
1954
1955	if (nlogprt .ge. 2) then
1956	write (UNIT = nulou,FMT = *)' '
1957	write (UNIT = nulou,FMT = *)'Leaving routine remap_conserv'
1958	call OASIS_FLUSH_SCRIP(nulou)
1959	endif
1960
1961	!-----------------------------------------------------------------------
1962
1963	CONTAINS
1964
1965	logical (kind=log_kind) function lf_loncross(dgbbp,dgbb3,dgbb4,grid_add,&
1966	bbox_per, bound_box)
1967
1968	!-----------------------------------------------------------------------
1969	!
1970	! this function tests wether the bounding box
1971	! of the currently checked cell at address grid2_add in grid2
1972	! (grid1_add in grid1) intersects in longitude the destination cell
1973	! on grid1 (grid2).
1974	! Periodicity is taken into account
1975	!
1976	!-----------------------------------------------------------------------
1977
1978	integer (kind=int_kind), INTENT(IN) :: dgbbp, grid_add
1979	real (kind=dbl_kind) , INTENT(IN) :: dgbb3, dgbb4
1980	integer (kind=int_kind), dimension(:) , INTENT(IN) :: bbox_per
1981	real (kind=dbl_kind), dimension(:,:) , INTENT(IN) :: bound_box
1982
1983	!--------------------------------------------------------------------
1984	!
1985	! Work variables (argument)
1986	!
1987	! dgbpp : prefetched periodicity class of the destination grid cell
1988	! -1 : the west boundary is < 0
1989	! 0 : both boundaries are in [0,2pi]
1990	! 1 : the east boundary is > 2pi
1991	!
1992	! dgbb3, dgbb4 : prefetched west and east longitudes of destination
1993	! cell
1994	!
1995	! grid_add : address of the source grid
1996	!
1997	! grid_bbox_per : periodicity class of the source grid cells
1998	!
1999	! grid_bound_box : coordinates of the source grid bounding boxes
2000	!
2001	!--------------------------------------------------------------------
2002
2003	!--------------------------------------------------------------------
2004	!
2005	! Principle of the test
2006	!
2007	! A source cell with longitude boundaries [w_s, e_s] intersects
2008	! a destination cell with boundaries [w_d, e_d] if simultaneously
2009	! w_s is west of e_d and e_s is east of w_d
2010	! If one (and only one) of the cells crosses the periodicity
2011	! threshold the test checks if the other cell crosses its boundary
2012	! inside [0,2pi] (with the standard test) or the boundary stored
2013	! with periodicity (by a reverted test on the modulo image of the
2014	! boundary).
2015	!
2016	!--------------------------------------------------------------------
2017
2018	lf_loncross = .true.
2019
2020	select case(dgbbp)
2021	case(-1)
2022
2023	! The destination cell crosses 0
2024	! It necessarily intersects source cells crossing 0 or 2pi.
2025	! It is west of any cell in [0,2pi]: only the w_s is west of e_d
2026	! test is needed (periodicity accounted for)
2027
2028	if (bbox_per(grid_add).eq.0) then
2029	lf_loncross = bound_box(3,grid_add) <= dgbb4 .or. &
2030	& bound_box(4,grid_add) >= dgbb3 + pi2
2031	end if
2032	case(0)
2033
2034	! The destination cell is in [0,2pi]
2035
2036	select case(bbox_per(grid_add))
2037	case(-1)
2038
2039	! The source cell crosses 0
2040	! Only the e_s is east of w_d test is needed
2041	! (periodicity accounted for)
2042
2043	lf_loncross = bound_box(4,grid_add) >= dgbb3 .or. &
2044	& bound_box(3,grid_add) + pi2 <= dgbb4
2045	case(0)
2046
2047	! The source cell is in [0,2pi]
2048	! Standard test
2049
2050	lf_loncross = bound_box(3,grid_add) <= dgbb4 ! w_s is west of e_d
2051	if (lf_loncross) &
2052	& lf_loncross = bound_box(4,grid_add) >= dgbb3 ! e_s is east of w_d
2053	case(1)
2054
2055	! The source cell crosses 2pi
2056	! Only the w_s is west of e_d test is needed
2057	! (periodicity accounted for)
2058
2059	lf_loncross = bound_box(3,grid_add) <= dgbb4 .or. &
2060	& bound_box(4,grid_add) - pi2 >= dgbb3
2061	end select
2062	case(1)
2063
2064	! The destination cell crosses 2pi
2065	! It necessarily intersects source cells crossing 0 or 2pi.
2066	! It is east of any cell in [0,2pi]: only the e_s is east of w_d
2067	! test is needed (periodicity accounted for)
2068
2069	if (bbox_per(grid_add).eq.0) then
2070	lf_loncross = bound_box(4,grid_add) >= dgbb3 .or. &
2071	& bound_box(3,grid_add) <= dgbb4 - pi2
2072	end if
2073	end select
2074
2075	end function lf_loncross
2076
2077	end subroutine remap_conserv
2078
2079	!***********************************************************************
2080
2081	subroutine intersection(location,intrsct_lat,intrsct_lon,lcoinc, &
2082	beglat, beglon, endlat, endlon, begseg, &
2083	lbegin, lrevers, &
2084	srch_corners, num_srch_cells, &
2085	srch_add, grid_corner_lon, grid_corner_lat, &
2086	intrsct_lat_off, intrsct_lon_off)
2087
2088	!-----------------------------------------------------------------------
2089	!
2090	! this routine finds the next intersection of a destination grid
2091	! line with the line segment given by beglon, endlon, etc.
2092	! a coincidence flag is returned if the segment is entirely
2093	! coincident with an ocean grid line. the cells in which to search
2094	! for an intersection must have already been restricted in the
2095	! calling routine.
2096	!
2097	!-----------------------------------------------------------------------
2098
2099	!-----------------------------------------------------------------------
2100	!
2101	! intent(in):
2102	!
2103	!-----------------------------------------------------------------------
2104
2105	logical (kind=log_kind), intent(in) :: lbegin, & ! flag for first integration along this segment
2106	lrevers ! flag whether segment integrated in reverse
2107
2108	real (kind=dbl_kind), intent(in) :: beglat, beglon, & ! beginning lat/lon endpoints for segment
2109	endlat, endlon ! ending lat/lon endpoints for segment
2110
2111	real (kind=dbl_kind), dimension(2), intent(inout) :: begseg ! begin lat/lon of full segment
2112	integer (kind=int_kind), intent(in) :: srch_corners ! Nb of corners in sarched grid
2113
2114	integer (kind=int_kind), intent(in) :: num_srch_cells ! Nb of preselected search cells
2115
2116	integer (kind=int_kind), dimension(:), intent(in) :: srch_add ! Address of search cells
2117
2118	real (kind=dbl_kind), dimension(:,:), intent(in) :: grid_corner_lat,&
2119	& grid_corner_lon ! lat and lon of each corner of srch grids
2120
2121	!-----------------------------------------------------------------------
2122	!
2123	! intent(out):
2124	!
2125	!-----------------------------------------------------------------------
2126
2127	integer (kind=int_kind), intent(out) :: location ! address in destination array containing this
2128	! segment
2129
2130	logical (kind=log_kind), intent(out) :: lcoinc ! flag segments which are entirely coincident
2131	! with a grid line
2132
2133	real (kind=dbl_kind), intent(out) :: intrsct_lat, intrsct_lon ! lat/lon coords of next intersect.
2134	real (kind=dbl_kind), intent(inout) :: intrsct_lat_off, intrsct_lon_off ! lat/lon coords offset
2135
2136	!-----------------------------------------------------------------------
2137	!
2138	! local variables
2139	!
2140	!-----------------------------------------------------------------------
2141
2142	integer (kind=int_kind) :: n, next_n, cell, gcell
2143	! for test of non-convexe cell
2144	integer (kind=int_kind) :: next2_n, neg, pos
2145
2146	integer (kind=int_kind), save :: last_loc ! save location when crossing threshold
2147	!$OMP THREADPRIVATE(last_loc)
2148
2149	logical (kind=log_kind) :: loutside ! flags points outside grid
2150
2151	logical (kind=log_kind), save :: lthresh = .false. ! flags segments crossing threshold bndy
2152	!$OMP THREADPRIVATE(lthresh)
2153
2154	real (kind=dbl_kind) :: lon1, lon2, & ! local longitude variables for segment
2155	lat1, lat2, & ! local latitude variables for segment
2156	grdlon1, grdlon2, & ! local longitude variables for grid cell
2157	grdlat1, grdlat2, & ! local latitude variables for grid cell
2158	vec1_lat, vec1_lon, & ! vectors and cross products used
2159	vec2_lat, vec2_lon, & ! during grid search
2160	cross_product, &
2161	eps, offset, & ! small offset away from intersect
2162	s1, s2, determ, & ! variables used for linear solve to
2163	mat1, mat2, mat3, mat4, rhs1, rhs2, & ! find intersection
2164	rl_halfpi, rl_v2lonmpi2, rl_v2lonppi2
2165
2166	! for next search
2167
2168	!-----------------------------------------------------------------------
2169	!
2170	! initialize defaults, flags, etc.
2171	!
2172	!-----------------------------------------------------------------------
2173
2174	location = 0
2175	lcoinc = .false.
2176	intrsct_lat = endlat
2177	intrsct_lon = endlon
2178
2179	if (num_srch_cells == 0) return
2180
2181	if (beglat > north_thresh .or. beglat < south_thresh) then
2182
2183	if (lthresh) location = last_loc
2184	call pole_intersection(location, &
2185	intrsct_lat,intrsct_lon,lcoinc,lthresh, &
2186	beglat, beglon, endlat, endlon, begseg, lrevers, &
2187	srch_corners, num_srch_cells, &
2188	srch_add, grid_corner_lon, grid_corner_lat)
2189	if (lthresh) then
2190	last_loc = location
2191	intrsct_lat_off = intrsct_lat
2192	intrsct_lon_off = intrsct_lon
2193	endif
2194	return
2195
2196	endif
2197
2198	loutside = .false.
2199	if (lbegin) then
2200	lat1 = beglat
2201	lon1 = beglon
2202	else
2203	lat1 = intrsct_lat_off
2204	lon1 = intrsct_lon_off
2205	endif
2206	lat2 = endlat
2207	lon2 = endlon
2208	if ((lon2-lon1) > three*pih) then
2209	lon2 = lon2 - pi2
2210	else if ((lon2-lon1) < -three*pih) then
2211	lon2 = lon2 + pi2
2212	endif
2213	s1 = zero
2214
2215	!-----------------------------------------------------------------------
2216	!
2217	! search for location of this segment in ocean grid using cross
2218	! product method to determine whether a point is enclosed by a cell
2219	!
2220	!-----------------------------------------------------------------------
2221
2222
2223	srch_loop: do
2224
2225	!***
2226	!*** if last segment crossed threshold, use that location
2227	!***
2228
2229	if (lthresh) then
2230
2231	do cell=1,num_srch_cells
2232
2233	if (srch_add(cell) == last_loc) then
2234	location = last_loc
2235	eps = tiny
2236	exit srch_loop
2237	endif
2238	end do
2239	endif
2240
2241	!***
2242	!*** otherwise normal search algorithm
2243	!***
2244
2245	cell_loop: do cell=1,num_srch_cells
2246	gcell = srch_add(cell)
2247	lcoinc = .false.
2248	corner_loop: do n=1,srch_corners
2249	next_n = mod(n,srch_corners) + 1
2250
2251	!***
2252	!*** here we take the cross product of the vector making
2253	!*** up each cell side with the vector formed by the vertex
2254	!*** and search point. if all the cross products are
2255	!*** positive, the point is contained in the cell.
2256	!***
2257
2258	vec1_lat = grid_corner_lat(next_n,gcell) - grid_corner_lat(n ,gcell)
2259	vec1_lon = grid_corner_lon(next_n,gcell) - grid_corner_lon(n ,gcell)
2260	vec2_lat = lat1 - grid_corner_lat(n,gcell)
2261	vec2_lon = lon1 - grid_corner_lon(n,gcell)
2262
2263	!***
2264	!*** if endpoint coincident with vertex, offset
2265	!*** the endpoint
2266	!***
2267
2268	if (vec2_lat == 0 .and. vec2_lon == 0) then
2269	lat1 = lat1 + 1.d-10*(lat2-lat1)
2270	lon1 = lon1 + 1.d-10*(lon2-lon1)
2271	vec2_lat = lat1 - grid_corner_lat(n,gcell)
2272	vec2_lon = lon1 - grid_corner_lon(n,gcell)
2273	endif
2274
2275	!***
2276	!*** check for 0,2pi crossings
2277	!***
2278
2279	if (vec1_lon > pi) then
2280	vec1_lon = vec1_lon - pi2
2281	else if (vec1_lon < -pi) then
2282	vec1_lon = vec1_lon + pi2
2283	endif
2284	if (vec2_lon > pi) then
2285	vec2_lon = vec2_lon - pi2
2286	else if (vec2_lon < -pi) then
2287	vec2_lon = vec2_lon + pi2
2288	endif
2289
2290	cross_product = vec1_lonvec2_lat - vec2_lonvec1_lat
2291	!***
2292	!*** if the cross product for a side is zero, the point
2293	!*** lies exactly on the side or the side is degenerate
2294	!*** (zero length). if degenerate, set the cross
2295	!*** product to a positive number. otherwise perform
2296	!*** another cross product between the side and the
2297	!*** segment itself.
2298	!*** if this cross product is also zero, the line is
2299	!*** coincident with the cell boundary - perform the
2300	!*** dot product and only choose the cell if the dot
2301	!*** product is positive (parallel vs anti-parallel).
2302	!***
2303
2304	if (cross_product == zero) then
2305	if (vec1_lat /= zero .or. vec1_lon /= zero) then
2306	vec2_lat = lat2 - lat1
2307	vec2_lon = lon2 - lon1
2308
2309	if (vec2_lon > pi) then
2310	vec2_lon = vec2_lon - pi2
2311	else if (vec2_lon < -pi) then
2312	vec2_lon = vec2_lon + pi2
2313	endif
2314
2315	cross_product = vec1_lonvec2_lat - vec2_lonvec1_lat
2316	else
2317	cross_product = one
2318	endif
2319
2320	if (cross_product == zero) then
2321	lcoinc = .true.
2322	cross_product = vec1_lonvec2_lon + vec1_latvec2_lat
2323	if (lrevers) cross_product = -cross_product
2324	endif
2325	endif
2326
2327	!***
2328	!*** if cross product is less than zero, this cell
2329	!*** doesn't work
2330	!***
2331
2332	if (cross_product < zero) exit corner_loop
2333
2334	end do corner_loop
2335
2336	!***
2337	!*** if cross products all positive, we found the location
2338	!***
2339
2340	if (n > srch_corners) then
2341	location = srch_add(cell)
2342
2343	!***
2344	!*** if the beginning of this segment was outside the
2345	!*** grid, invert the segment so the intersection found
2346	!*** will be the first intersection with the grid
2347	!***
2348	if (loutside) then
2349	!*** do a test to see if the cell really is outside
2350	!*** or if it is a non-convexe cell
2351	neg=0
2352	pos=0
2353	do n=1,srch_corners
2354	next_n = mod(n,srch_corners) + 1
2355	next2_n = mod(next_n,srch_corners) + 1
2356
2357	vec1_lat = grid_corner_lat(next_n,gcell) - grid_corner_lat(n,gcell)
2358	vec1_lon = grid_corner_lon(next_n,gcell) - grid_corner_lon(n,gcell)
2359	vec2_lat = grid_corner_lat(next2_n,gcell) - grid_corner_lat(next_n,gcell)
2360	vec2_lon = grid_corner_lon(next2_n,gcell) - grid_corner_lon(next_n,gcell)
2361
2362	if (vec1_lon > three*pih) then
2363	vec1_lon = vec1_lon - pi2
2364	else if (vec1_lon < -three*pih) then
2365	vec1_lon = vec1_lon + pi2
2366	endif
2367
2368	if (vec2_lon > three*pih) then
2369	vec2_lon = vec2_lon - pi2
2370	else if (vec2_lon < -three*pih) then
2371	vec2_lon = vec2_lon + pi2
2372	endif
2373
2374	cross_product = vec1_latvec2_lon - vec2_latvec1_lon
2375
2376	if (cross_product < zero) then
2377	neg=neg+1
2378	else if (cross_product > zero) then
2379	pos=pos+1
2380	endif
2381	enddo
2382	!*** the cell is non-convexe if not all cross_products
2383	!*** have the same signe
2384	if (neg/=0 .and. pos/=0) then
2385	loutside=.false.
2386	if (nlogprt .ge. 2) then
2387	write(nulou,*) 'The mesh ',gcell,' is non-convex'
2388	write(nulou,*) 'srch_corner_lat=',grid_corner_lat(:,gcell)
2389	write(nulou,*) 'srch_corner_lon=',grid_corner_lon(:,gcell)
2390	call OASIS_FLUSH_SCRIP(nulou)
2391	endif
2392	endif
2393	endif
2394	if (loutside) then
2395	lat2 = beglat
2396	lon2 = beglon
2397	location = 0
2398	eps = -tiny
2399	else
2400	eps = tiny
2401	endif
2402
2403	exit srch_loop
2404	endif
2405
2406	!***
2407	!*** otherwise move on to next cell
2408	!***
2409
2410	end do cell_loop
2411
2412	!***
2413	!*** if still no cell found, the point lies outside the grid.
2414	!*** take some baby steps along the segment to see if any
2415	!*** part of the segment lies inside the grid.
2416	!***
2417
2418	loutside = .true.
2419	s1 = s1 + 0.001_dbl_kind
2420	lat1 = beglat + s1*(endlat - beglat)
2421	lon1 = beglon + s1*(lon2 - beglon)
2422
2423	!***
2424	!*** reached the end of the segment and still outside the grid
2425	!*** return no intersection
2426	!***
2427
2428	if (s1 >= one) return
2429
2430	end do srch_loop
2431
2432	!-----------------------------------------------------------------------
2433	!
2434	! now that a cell is found, search for the next intersection.
2435	! loop over sides of the cell to find intersection with side
2436	! must check all sides for coincidences or intersections
2437	!
2438	!-----------------------------------------------------------------------
2439
2440	intrsct_loop: do n=1,srch_corners
2441	next_n = mod(n,srch_corners) + 1
2442
2443	grdlon1 = grid_corner_lon(n ,gcell)
2444	grdlon2 = grid_corner_lon(next_n,gcell)
2445	grdlat1 = grid_corner_lat(n ,gcell)
2446	grdlat2 = grid_corner_lat(next_n,gcell)
2447
2448	!***
2449	!*** set up linear system to solve for intersection
2450	!***
2451
2452	mat1 = lat2 - lat1
2453	mat2 = grdlat1 - grdlat2
2454	mat3 = lon2 - lon1
2455	mat4 = grdlon1 - grdlon2
2456	rhs1 = grdlat1 - lat1
2457	rhs2 = grdlon1 - lon1
2458
2459	if (mat3 > pi) then
2460	mat3 = mat3 - pi2
2461	else if (mat3 < -pi) then
2462	mat3 = mat3 + pi2
2463	endif
2464	if (mat4 > pi) then
2465	mat4 = mat4 - pi2
2466	else if (mat4 < -pi) then
2467	mat4 = mat4 + pi2
2468	endif
2469	if (rhs2 > pi) then
2470	rhs2 = rhs2 - pi2
2471	else if (rhs2 < -pi) then
2472	rhs2 = rhs2 + pi2
2473	endif
2474
2475	determ = mat1mat4 - mat2mat3
2476
2477	!***
2478	!*** if the determinant is zero, the segments are either
2479	!*** parallel or coincident. coincidences were detected
2480	!*** above so do nothing.
2481	!*** if the determinant is non-zero, solve for the linear
2482	!*** parameters s for the intersection point on each line
2483	!*** segment.
2484	!*** if 0<s1,s2<1 then the segment intersects with this side.
2485	!*** return the point of intersection (adding a small
2486	!*** number so the intersection is off the grid line).
2487	!***
2488
2489	if (abs(determ) > 1.e-30) then
2490
2491	s1 = (rhs1mat4 - mat2rhs2)/determ
2492	s2 = (mat1rhs2 - rhs1mat3)/determ
2493
2494	!EM bug F77:
2495	! if (s2 >= zero .and. s2 <= one .and.
2496	! & s1 > zero. and. s1 <= one) then
2497
2498	if (s2 >= zero .and. s2 <= one .and. &
2499	s1 > zero .and. s1 <= one) then
2500
2501	!***
2502	!*** recompute intersection based on full segment
2503	!*** so intersections are consistent for both sweeps
2504	!***
2505
2506	if (.not. loutside) then
2507	mat1 = lat2 - begseg(1)
2508	mat3 = lon2 - begseg(2)
2509	rhs1 = grdlat1 - begseg(1)
2510	rhs2 = grdlon1 - begseg(2)
2511	else
2512	mat1 = begseg(1) - endlat
2513	mat3 = begseg(2) - endlon
2514	rhs1 = grdlat1 - endlat
2515	rhs2 = grdlon1 - endlon
2516	endif
2517
2518	if (mat3 > pi) then
2519	mat3 = mat3 - pi2
2520	else if (mat3 < -pi) then
2521	mat3 = mat3 + pi2
2522	endif
2523	if (rhs2 > pi) then
2524	rhs2 = rhs2 - pi2
2525	else if (rhs2 < -pi) then
2526	rhs2 = rhs2 + pi2
2527	endif
2528
2529	determ = mat1mat4 - mat2mat3
2530
2531	!***
2532	!*** sometimes due to roundoff, the previous
2533	!*** determinant is non-zero, but the lines
2534	!*** are actually coincident. if this is the
2535	!*** case, skip the rest.
2536	!***
2537
2538	if (determ /= zero) then
2539	s1 = (rhs1mat4 - mat2rhs2)/determ
2540	s2 = (mat1rhs2 - rhs1mat3)/determ
2541
2542	offset = s1 + eps/determ
2543	if (offset > one) offset = one
2544
2545	if (.not. loutside) then
2546	intrsct_lat = begseg(1) + mat1*s1
2547	intrsct_lon = begseg(2) + mat3*s1
2548	intrsct_lat_off = begseg(1) + mat1*offset
2549	intrsct_lon_off = begseg(2) + mat3*offset
2550	else
2551	intrsct_lat = endlat + mat1*s1
2552	intrsct_lon = endlon + mat3*s1
2553	intrsct_lat_off = endlat + mat1*offset
2554	intrsct_lon_off = endlon + mat3*offset
2555	endif
2556	exit intrsct_loop
2557	endif
2558
2559	endif
2560	endif
2561
2562	!***
2563	!*** no intersection this side, move on to next side
2564	!***
2565
2566	end do intrsct_loop
2567
2568	!-----------------------------------------------------------------------
2569	!
2570	! if the segment crosses a pole threshold, reset the intersection
2571	! to be the threshold latitude. only check if this was not a
2572	! threshold segment since sometimes coordinate transform can end
2573	! up on other side of threshold again.
2574	!
2575	!-----------------------------------------------------------------------
2576
2577	if (lthresh) then
2578	if (intrsct_lat < north_thresh .or. intrsct_lat > south_thresh) &
2579	lthresh = .false.
2580	else if (lat1 > zero .and. intrsct_lat > north_thresh) then
2581	intrsct_lat = north_thresh + tiny
2582	intrsct_lat_off = north_thresh + eps*mat1
2583	s1 = (intrsct_lat - begseg(1))/mat1
2584	intrsct_lon = begseg(2) + s1*mat3
2585	intrsct_lon_off = begseg(2) + (s1+eps)*mat3
2586	last_loc = location
2587	lthresh = .true.
2588	else if (lat1 < zero .and. intrsct_lat < south_thresh) then
2589	intrsct_lat = south_thresh - tiny
2590	intrsct_lat_off = south_thresh + eps*mat1
2591	s1 = (intrsct_lat - begseg(1))/mat1
2592	intrsct_lon = begseg(2) + s1*mat3
2593	intrsct_lon_off = begseg(2) + (s1+eps)*mat3
2594	last_loc = location
2595	lthresh = .true.
2596	endif
2597
2598	!-----------------------------------------------------------------------
2599
2600	end subroutine intersection
2601
2602	!***********************************************************************
2603
2604	subroutine pole_intersection(location, &
2605	intrsct_lat,intrsct_lon,lcoinc,lthresh, &
2606	beglat, beglon, endlat, endlon, begseg, lrevers, &
2607	srch_corners, num_srch_cells, &
2608	srch_add, grid_corner_lon, grid_corner_lat)
2609
2610	!-----------------------------------------------------------------------
2611	!
2612	! this routine is identical to the intersection routine except
2613	! that a coordinate transformation (using a Lambert azimuthal
2614	! equivalent projection) is performed to treat polar cells more
2615	! accurately.
2616	!
2617	!-----------------------------------------------------------------------
2618
2619	!-----------------------------------------------------------------------
2620	!
2621	! intent(in):
2622	!
2623	!-----------------------------------------------------------------------
2624
2625	real (kind=dbl_kind), intent(in) :: beglat, beglon, & ! beginning lat/lon endpoints for segment
2626	endlat, endlon ! ending lat/lon endpoints for segment
2627
2628	real (kind=dbl_kind), dimension(2), intent(inout) :: begseg ! begin lat/lon of full segment
2629
2630	logical (kind=log_kind), intent(in) :: lrevers ! flag true if segment integrated in reverse
2631
2632	integer (kind=int_kind), intent(in) :: srch_corners ! Nb of corners in sarched grid
2633
2634	integer (kind=int_kind), intent(in) :: num_srch_cells ! Nb of preselected search cells
2635
2636	integer (kind=int_kind), dimension(:), intent(in) :: srch_add ! Address of search cells
2637
2638	real (kind=dbl_kind), dimension(:,:), intent(in) :: grid_corner_lat,&
2639	& grid_corner_lon ! lat and lon of each corner of srch grids
2640
2641	!-----------------------------------------------------------------------
2642	!
2643	! intent(out):
2644	!
2645	!-----------------------------------------------------------------------
2646
2647	integer (kind=int_kind), intent(inout) :: location ! address in destination array containing this
2648	! segment -- also may contain last location on
2649	! entry
2650
2651	logical (kind=log_kind), intent(out) :: lcoinc ! flag segment coincident with grid line
2652
2653	logical (kind=log_kind), intent(inout) :: lthresh ! flag segment crossing threshold boundary
2654
2655	real (kind=dbl_kind), intent(out) :: intrsct_lat, intrsct_lon ! lat/lon coords of next intersect.
2656
2657	!-----------------------------------------------------------------------
2658	!
2659	! local variables
2660	!
2661	!-----------------------------------------------------------------------
2662
2663	integer (kind=int_kind) :: n, next_n, cell, gcell
2664
2665	logical (kind=log_kind) :: loutside ! flags points outside grid
2666
2667	real (kind=dbl_kind) :: pi4, rns, & ! north/south conversion
2668	x1, x2, & ! local x variables for segment
2669	y1, y2, & ! local y variables for segment
2670	begx, begy, & ! beginning x,y variables for segment
2671	endx, endy, & ! beginning x,y variables for segment
2672	begsegx, begsegy, & ! beginning x,y variables for segment
2673	grdx1, grdx2, & ! local x variables for grid cell
2674	grdy1, grdy2, & ! local y variables for grid cell
2675	vec1_y, vec1_x, & ! vectors and cross products used
2676	vec2_y, vec2_x, & ! during grid search
2677	cross_product, eps, & ! eps=small offset away from intersect
2678	s1, s2, determ, & ! variables used for linear solve to
2679	mat1, mat2, mat3, mat4, rhs1, rhs2 ! find intersection
2680
2681	real (kind=dbl_kind), dimension(srch_corners) :: srch_corner_x, & ! x of each corner of srch cells
2682	srch_corner_y ! y of each corner of srch cells
2683
2684	!***
2685	!*** save last intersection to avoid roundoff during coord
2686	!*** transformation
2687	!***
2688
2689	logical (kind=log_kind), save :: luse_last = .false.
2690
2691	real (kind=dbl_kind), save :: intrsct_x, intrsct_y ! x,y for intersection
2692
2693	!***
2694	!*** variables necessary if segment manages to hit pole
2695	!***
2696
2697	integer (kind=int_kind), save :: avoid_pole_count = 0 ! count attempts to avoid pole
2698
2699	real (kind=dbl_kind), save :: avoid_pole_offset = tiny ! endpoint offset to avoid pole
2700	!$OMP THREADPRIVATE(luse_last,intrsct_x,intrsct_y,avoid_pole_count,avoid_pole_offset)
2701
2702	!-----------------------------------------------------------------------
2703	!
2704	! initialize defaults, flags, etc.
2705	!
2706	!-----------------------------------------------------------------------
2707
2708	if (.not. lthresh) location = 0
2709	lcoinc = .false.
2710	intrsct_lat = endlat
2711	intrsct_lon = endlon
2712
2713	loutside = .false.
2714	s1 = zero
2715
2716	!-----------------------------------------------------------------------
2717	!
2718	! convert coordinates
2719	!
2720	!-----------------------------------------------------------------------
2721
2722	if (beglat > zero) then
2723	pi4 = quart*pi
2724	rns = one
2725	else
2726	pi4 = -quart*pi
2727	rns = -one
2728	endif
2729
2730	if (luse_last) then
2731	x1 = intrsct_x
2732	y1 = intrsct_y
2733	else
2734	x1 = rnstwosin(pi4 - halfbeglat)cos(beglon)
2735	y1 = twosin(pi4 - halfbeglat)*sin(beglon)
2736	luse_last = .true.
2737	endif
2738	x2 = rnstwosin(pi4 - halfendlat)cos(endlon)
2739	y2 = twosin(pi4 - halfendlat)*sin(endlon)
2740
2741	begx = x1
2742	begy = y1
2743	endx = x2
2744	endy = y2
2745	begsegx = rnstwosin(pi4 - halfbegseg(1))cos(begseg(2))
2746	begsegy = twosin(pi4 - halfbegseg(1))*sin(begseg(2))
2747	intrsct_x = endx
2748	intrsct_y = endy
2749
2750	!-----------------------------------------------------------------------
2751	!
2752	! search for location of this segment in ocean grid using cross
2753	! product method to determine whether a point is enclosed by a cell
2754	!
2755	!-----------------------------------------------------------------------
2756
2757	srch_loop: do
2758
2759	!***
2760	!*** if last segment crossed threshold, use that location
2761	!***
2762
2763	if (lthresh) then
2764	do cell=1,num_srch_cells
2765	if (srch_add(cell) == location) then
2766	eps = tiny
2767	exit srch_loop
2768	endif
2769	end do
2770	endif
2771
2772	!***
2773	!*** otherwise normal search algorithm
2774	!***
2775
2776	cell_loop: do cell=1,num_srch_cells
2777	gcell = srch_add(cell)
2778	lcoinc = .false.
2779	srch_corner_x(:) = rnstwosin(pi4 - halfgrid_corner_lat(:,gcell)) cos(grid_corner_lon(:,gcell))
2780	srch_corner_y(:) = twosin(pi4 - halfgrid_corner_lat(:,gcell)) * sin(grid_corner_lon(:,gcell))
2781	corner_loop: do n=1,srch_corners
2782	next_n = mod(n,srch_corners) + 1
2783
2784	!***
2785	!*** here we take the cross product of the vector making
2786	!*** up each cell side with the vector formed by the vertex
2787	!*** and search point. if all the cross products are
2788	!*** positive, the point is contained in the cell.
2789	!***
2790
2791	vec1_x = srch_corner_x(next_n) - srch_corner_x(n)
2792	vec1_y = srch_corner_y(next_n) - srch_corner_y(n)
2793	vec2_x = x1 - srch_corner_x(n)
2794	vec2_y = y1 - srch_corner_y(n)
2795
2796	!***
2797	!*** if endpoint coincident with vertex, offset
2798	!*** the endpoint
2799	!***
2800
2801	if (vec2_x == 0 .and. vec2_y == 0) then
2802	x1 = x1 + 1.d-10*(x2-x1)
2803	y1 = y1 + 1.d-10*(y2-y1)
2804	vec2_x = x1 - srch_corner_x(n)
2805	vec2_y = y1 - srch_corner_y(n)
2806	endif
2807
2808	cross_product = vec1_xvec2_y - vec2_xvec1_y
2809
2810	!***
2811	!*** if the cross product for a side is zero, the point
2812	!*** lies exactly on the side or the length of a side
2813	!*** is zero. if the length is zero set det > 0.
2814	!*** otherwise, perform another cross
2815	!*** product between the side and the segment itself.
2816	!*** if this cross product is also zero, the line is
2817	!*** coincident with the cell boundary - perform the
2818	!*** dot product and only choose the cell if the dot
2819	!*** product is positive (parallel vs anti-parallel).
2820	!***
2821
2822	if (cross_product == zero) then
2823	if (vec1_x /= zero .or. vec1_y /= 0) then
2824	vec2_x = x2 - x1
2825	vec2_y = y2 - y1
2826	cross_product = vec1_xvec2_y - vec2_xvec1_y
2827	else
2828	cross_product = one
2829	endif
2830
2831	if (cross_product == zero) then
2832	lcoinc = .true.
2833	cross_product = vec1_xvec2_x + vec1_yvec2_y
2834	if (lrevers) cross_product = -cross_product
2835	endif
2836	endif
2837
2838	!***
2839	!*** if cross product is less than zero, this cell
2840	!*** doesn't work
2841	!***
2842
2843	if (cross_product < zero) exit corner_loop
2844
2845	end do corner_loop
2846
2847	!***
2848	!*** if cross products all positive, we found the location
2849	!***
2850
2851	if (n > srch_corners) then
2852	location = srch_add(cell)
2853
2854	!***
2855	!*** if the beginning of this segment was outside the
2856	!*** grid, invert the segment so the intersection found
2857	!*** will be the first intersection with the grid
2858	!***
2859
2860	if (loutside) then
2861	x2 = begx
2862	y2 = begy
2863	location = 0
2864	eps = -tiny
2865	else
2866	eps = tiny
2867	endif
2868
2869	exit srch_loop
2870	endif
2871
2872	!***
2873	!*** otherwise move on to next cell
2874	!***
2875
2876	end do cell_loop
2877
2878	!***
2879	!*** if no cell found, the point lies outside the grid.
2880	!*** take some baby steps along the segment to see if any
2881	!*** part of the segment lies inside the grid.
2882	!***
2883
2884	loutside = .true.
2885	s1 = s1 + 0.001_dbl_kind
2886	x1 = begx + s1*(x2 - begx)
2887	y1 = begy + s1*(y2 - begy)
2888
2889	!***
2890	!*** reached the end of the segment and still outside the grid
2891	!*** return no intersection
2892	!***
2893
2894	if (s1 >= one) then
2895	luse_last = .false.
2896	return
2897	endif
2898
2899	end do srch_loop
2900
2901	!-----------------------------------------------------------------------
2902	!
2903	! now that a cell is found, search for the next intersection.
2904	! loop over sides of the cell to find intersection with side
2905	! must check all sides for coincidences or intersections
2906	!
2907	!-----------------------------------------------------------------------
2908
2909	gcell = srch_add(cell)
2910	srch_corner_x(:) = rnstwosin(pi4 - halfgrid_corner_lat(:,gcell)) cos(grid_corner_lon(:,gcell))
2911	srch_corner_y(:) = twosin(pi4 - halfgrid_corner_lat(:,gcell)) * sin(grid_corner_lon(:,gcell))
2912
2913	intrsct_loop: do n=1,srch_corners
2914	next_n = mod(n,srch_corners) + 1
2915
2916	grdy1 = srch_corner_y(n )
2917	grdy2 = srch_corner_y(next_n)
2918	grdx1 = srch_corner_x(n )
2919	grdx2 = srch_corner_x(next_n)
2920
2921	!***
2922	!*** set up linear system to solve for intersection
2923	!***
2924
2925	mat1 = x2 - x1
2926	mat2 = grdx1 - grdx2
2927	mat3 = y2 - y1
2928	mat4 = grdy1 - grdy2
2929	rhs1 = grdx1 - x1
2930	rhs2 = grdy1 - y1
2931
2932	determ = mat1mat4 - mat2mat3
2933
2934	!***
2935	!*** if the determinant is zero, the segments are either
2936	!*** parallel or coincident or one segment has zero length.
2937	!*** coincidences were detected above so do nothing.
2938	!*** if the determinant is non-zero, solve for the linear
2939	!*** parameters s for the intersection point on each line
2940	!*** segment.
2941	!*** if 0<s1,s2<1 then the segment intersects with this side.
2942	!*** return the point of intersection (adding a small
2943	!*** number so the intersection is off the grid line).
2944	!***
2945
2946	if (abs(determ) > 1.e-30) then
2947
2948	s1 = (rhs1mat4 - mat2rhs2)/determ
2949	s2 = (mat1rhs2 - rhs1mat3)/determ
2950
2951	!EM bug F77
2952	! if (s2 >= zero .and. s2 <= one .and. &
2953	! & s1 > zero. and. s1 <= one) then
2954	if (s2 >= zero .and. s2 <= one .and. &
2955	s1 > zero .and. s1 <= one) then
2956
2957	!***
2958	!*** recompute intersection using entire segment
2959	!*** for consistency between sweeps
2960	!***
2961
2962	if (.not. loutside) then
2963	mat1 = x2 - begsegx
2964	mat3 = y2 - begsegy
2965	rhs1 = grdx1 - begsegx
2966	rhs2 = grdy1 - begsegy
2967	else
2968	mat1 = x2 - endx
2969	mat3 = y2 - endy
2970	rhs1 = grdx1 - endx
2971	rhs2 = grdy1 - endy
2972	endif
2973
2974	determ = mat1mat4 - mat2mat3
2975
2976	!***
2977	!*** sometimes due to roundoff, the previous
2978	!*** determinant is non-zero, but the lines
2979	!*** are actually coincident. if this is the
2980	!*** case, skip the rest.
2981	!***
2982
2983	if (determ /= zero) then
2984	s1 = (rhs1mat4 - mat2rhs2)/determ
2985	s2 = (mat1rhs2 - rhs1mat3)/determ
2986
2987	if (.not. loutside) then
2988	intrsct_x = begsegx + s1*mat1
2989	intrsct_y = begsegy + s1*mat3
2990	else
2991	intrsct_x = endx + s1*mat1
2992	intrsct_y = endy + s1*mat3
2993	endif
2994
2995	!***
2996	!*** convert back to lat/lon coordinates
2997	!***
2998
2999	intrsct_lon = rns*atan2(intrsct_y,intrsct_x)
3000	if (intrsct_lon < zero) &
3001	intrsct_lon = intrsct_lon + pi2
3002
3003	if (abs(intrsct_x) > 1.d-10) then
3004	intrsct_lat = (pi4 - asin(rnshalfintrsct_x/cos(intrsct_lon)))*two
3005	else if (abs(intrsct_y) > 1.d-10) then
3006	intrsct_lat = (pi4 - asin(halfintrsct_y/sin(intrsct_lon)))two
3007	else
3008	intrsct_lat = two*pi4
3009	endif
3010
3011	!***
3012	!*** add offset in transformed space for next pass.
3013	!***
3014
3015	if (s1 - eps/determ < one) then
3016	intrsct_x = intrsct_x - mat1*(eps/determ)
3017	intrsct_y = intrsct_y - mat3*(eps/determ)
3018	else
3019	if (.not. loutside) then
3020	intrsct_x = endx
3021	intrsct_y = endy
3022	intrsct_lat = endlat
3023	intrsct_lon = endlon
3024	else
3025	intrsct_x = begsegx
3026	intrsct_y = begsegy
3027	intrsct_lat = begseg(1)
3028	intrsct_lon = begseg(2)
3029	endif
3030	endif
3031
3032	exit intrsct_loop
3033	endif
3034	endif
3035	endif
3036
3037	!***
3038	!*** no intersection this side, move on to next side
3039	!***
3040
3041	end do intrsct_loop
3042
3043	!-----------------------------------------------------------------------
3044	!
3045	! if segment manages to cross over pole, shift the beginning
3046	! endpoint in order to avoid hitting pole directly
3047	! (it is ok for endpoint to be pole point)
3048	!
3049	!-----------------------------------------------------------------------
3050
3051	if (abs(intrsct_x) < 1.e-10 .and. abs(intrsct_y) < 1.e-10 .and. &
3052	(endx /= zero .and. endy /=0)) then
3053	if (avoid_pole_count > 2) then
3054	avoid_pole_count = 0
3055	avoid_pole_offset = 10.*avoid_pole_offset
3056	endif
3057
3058	cross_product = begsegx(endy-begsegy) - begsegy(endx-begsegx)
3059	intrsct_lat = begseg(1)
3060	if (cross_product*intrsct_lat > zero) then
3061	intrsct_lon = beglon + avoid_pole_offset
3062	begseg(2) = begseg(2) + avoid_pole_offset
3063	else
3064	intrsct_lon = beglon - avoid_pole_offset
3065	begseg(2) = begseg(2) - avoid_pole_offset
3066	endif
3067
3068	avoid_pole_count = avoid_pole_count + 1
3069	luse_last = .false.
3070	else
3071	avoid_pole_count = 0
3072	avoid_pole_offset = tiny
3073	endif
3074
3075	!-----------------------------------------------------------------------
3076	!
3077	! if the segment crosses a pole threshold, reset the intersection
3078	! to be the threshold latitude and do not reuse x,y intersect
3079	! on next entry. only check if did not cross threshold last
3080	! time - sometimes the coordinate transformation can place a
3081	! segment on the other side of the threshold again
3082	!
3083	!-----------------------------------------------------------------------
3084
3085	if (lthresh) then
3086	if (intrsct_lat > north_thresh .or. intrsct_lat < south_thresh) &
3087	lthresh = .false.
3088	else if (beglat > zero .and. intrsct_lat < north_thresh) then
3089	mat4 = endlat - begseg(1)
3090	mat3 = endlon - begseg(2)
3091	if (mat3 > pi) mat3 = mat3 - pi2
3092	if (mat3 < -pi) mat3 = mat3 + pi2
3093	intrsct_lat = north_thresh - tiny
3094	s1 = (north_thresh - begseg(1))/mat4
3095	intrsct_lon = begseg(2) + s1*mat3
3096	luse_last = .false.
3097	lthresh = .true.
3098	else if (beglat < zero .and. intrsct_lat > south_thresh) then
3099	mat4 = endlat - begseg(1)
3100	mat3 = endlon - begseg(2)
3101	if (mat3 > pi) mat3 = mat3 - pi2
3102	if (mat3 < -pi) mat3 = mat3 + pi2
3103	intrsct_lat = south_thresh + tiny
3104	s1 = (south_thresh - begseg(1))/mat4
3105	intrsct_lon = begseg(2) + s1*mat3
3106	luse_last = .false.
3107	lthresh = .true.
3108	endif
3109
3110	!***
3111	!*** if reached end of segment, do not use x,y intersect
3112	!*** on next entry
3113	!***
3114
3115	if (intrsct_lat == endlat .and. intrsct_lon == endlon) then
3116	luse_last = .false.
3117	endif
3118
3119	!-----------------------------------------------------------------------
3120
3121	end subroutine pole_intersection
3122
3123	!***********************************************************************
3124
3125	subroutine line_integral(weights, num_wts, &
3126	in_phi1, in_phi2, theta1, theta2, &
3127	grid1_lon, grid2_lon)
3128
3129	!-----------------------------------------------------------------------
3130	!
3131	! this routine computes the line integral of the flux function
3132	! that results in the interpolation weights. the line is defined
3133	! by the input lat/lon of the endpoints.
3134	!
3135	!-----------------------------------------------------------------------
3136
3137	!-----------------------------------------------------------------------
3138	!
3139	! intent(in):
3140	!
3141	!-----------------------------------------------------------------------
3142
3143	integer (kind=int_kind), intent(in) :: num_wts ! number of weights to compute
3144
3145	real (kind=dbl_kind), intent(in) :: in_phi1, in_phi2, & ! longitude endpoints for the segment
3146	theta1, theta2, & ! latitude endpoints for the segment
3147	grid1_lon, & ! reference coordinates for each
3148	grid2_lon ! grid (to ensure correct 0,2pi interv.
3149
3150	!-----------------------------------------------------------------------
3151	!
3152	! intent(out):
3153	!
3154	!-----------------------------------------------------------------------
3155
3156	real (kind=dbl_kind), dimension(2*num_wts), intent(out) :: weights ! line integral contribution to weights
3157
3158	!-----------------------------------------------------------------------
3159	!
3160	! local variables
3161	!
3162	!-----------------------------------------------------------------------
3163
3164	real (kind=dbl_kind) :: dphi, sinth1, sinth2, costh1, costh2, fac, phi1, phi2, phidiff1, phidiff2
3165	real (kind=dbl_kind) :: f1, f2, fint
3166
3167	!-----------------------------------------------------------------------
3168	!
3169	! weights for the general case based on a trapezoidal approx to
3170	! the integrals.
3171	!
3172	!-----------------------------------------------------------------------
3173
3174	sinth1 = sin(theta1)
3175	sinth2 = sin(theta2)
3176	costh1 = cos(theta1)
3177	costh2 = cos(theta2)
3178
3179	dphi = in_phi1 - in_phi2
3180	if (dphi > pi) then
3181	dphi = dphi - pi2
3182	else if (dphi < -pi) then
3183	dphi = dphi + pi2
3184	endif
3185	dphi = half*dphi
3186
3187	!-----------------------------------------------------------------------
3188	!
3189	! the first weight is the area overlap integral. the second and
3190	! fourth are second-order latitude gradient weights.
3191	!
3192	!-----------------------------------------------------------------------
3193
3194	weights( 1) = dphi*(sinth1 + sinth2)
3195	weights(num_wts+1) = dphi*(sinth1 + sinth2)
3196	weights( 2) = dphi(costh1 + costh2 + (theta1sinth1 + theta2*sinth2))
3197	weights(num_wts+2) = dphi(costh1 + costh2 + (theta1sinth1 + theta2*sinth2))
3198
3199	!-----------------------------------------------------------------------
3200	!
3201	! the third and fifth weights are for the second-order phi gradient
3202	! component. must be careful of longitude range.
3203	!
3204	!-----------------------------------------------------------------------
3205
3206	f1 = half(costh1sinth1 + theta1)
3207	f2 = half(costh2sinth2 + theta2)
3208
3209	phi1 = in_phi1 - grid1_lon
3210	if (phi1 > pi) then
3211	phi1 = phi1 - pi2
3212	else if (phi1 < -pi) then
3213	phi1 = phi1 + pi2
3214	endif
3215
3216	phi2 = in_phi2 - grid1_lon
3217	if (phi2 > pi) then
3218	phi2 = phi2 - pi2
3219	else if (phi2 < -pi) then
3220	phi2 = phi2 + pi2
3221	endif
3222
3223	if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3224	weights(3) = dphi(phi1f1 + phi2*f2)
3225	else
3226	if (phi1 > zero) then
3227	fac = pi
3228	else
3229	fac = -pi
3230	endif
3231	fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
3232	weights(3) = halfphi1(phi1-fac)*f1 - &
3233	halfphi2(phi2+fac)*f2 + &
3234	halffac(phi1+phi2)*fint
3235	endif
3236
3237	phi1 = in_phi1 - grid2_lon
3238	if (phi1 > pi) then
3239	phi1 = phi1 - pi2
3240	else if (phi1 < -pi) then
3241	phi1 = phi1 + pi2
3242	endif
3243
3244	phi2 = in_phi2 - grid2_lon
3245	if (phi2 > pi) then
3246	phi2 = phi2 - pi2
3247	else if (phi2 < -pi) then
3248	phi2 = phi2 + pi2
3249	endif
3250
3251	if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3252	weights(num_wts+3) = dphi(phi1f1 + phi2*f2)
3253	else
3254	if (phi1 > zero) then
3255	fac = pi
3256	else
3257	fac = -pi
3258	endif
3259	fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
3260	weights(num_wts+3) = halfphi1(phi1-fac)*f1 - &
3261	halfphi2(phi2+fac)*f2 + &
3262	halffac(phi1+phi2)*fint
3263	endif
3264
3265	!-----------------------------------------------------------------------
3266
3267	end subroutine line_integral
3268
3269	!***********************************************************************
3270
3271	subroutine store_link_cnsrv(add1, add2, weights, id_thread)
3272
3273	!-----------------------------------------------------------------------
3274	!
3275	! this routine stores the address and weight for this link in
3276	! the appropriate address and weight arrays and resizes those
3277	! arrays if necessary.
3278	!
3279	!-----------------------------------------------------------------------
3280
3281	!-----------------------------------------------------------------------
3282	!
3283	! input variables
3284	!
3285	!-----------------------------------------------------------------------
3286
3287	integer (kind=int_kind), intent(in) :: add1, & ! address on grid1
3288	add2, & ! address on grid2
3289	id_thread ! local threaded task nb
3290
3291	real (kind=dbl_kind), dimension(:), intent(in) :: weights ! array of remapping weights for this link
3292
3293	!-----------------------------------------------------------------------
3294	!
3295	! local variables
3296	!
3297	!-----------------------------------------------------------------------
3298
3299	integer (kind=int_kind) :: nlink, min_link, max_link ! link index
3300
3301	integer (kind=int_kind), dimension(:,:), allocatable, save :: link_add1, & ! min,max link add to restrict search
3302	link_add2 ! min,max link add to restrict search
3303
3304	!-----------------------------------------------------------------------
3305	!
3306	! if all weights are zero, do not bother storing the link
3307	!
3308	!-----------------------------------------------------------------------
3309
3310	!SV if (all(weights == zero)) return
3311
3312	!-----------------------------------------------------------------------
3313	!
3314	! if the link does not yet exist, increment number of links and
3315	! check to see if remap arrays need to be increased to accomodate
3316	! the new link. then store the link.
3317	!
3318	!-----------------------------------------------------------------------
3319
3320	if (il_nbthreads .eq. 1) then
3321	num_links_map1 = num_links_map1 + 1
3322	if (num_links_map1 > max_links_map1) &
3323	call resize_remap_vars(1,resize_increment)
3324
3325	grid1_add_map1(num_links_map1) = add1
3326	grid2_add_map1(num_links_map1) = add2
3327	wts_map1 (:,num_links_map1) = weights(1:num_wts)
3328
3329	else
3330	sga_remap(id_thread)%num_links = sga_remap(id_thread)%num_links + 1
3331
3332	if (sga_remap(id_thread)%num_links > sga_remap(id_thread)%max_links) &
3333	call sga_remap(id_thread)%resize(int(0.2*real(sga_remap(id_thread)%max_links)))
3334
3335	sga_remap(id_thread)%grid1_add(sga_remap(id_thread)%num_links) = add1
3336	sga_remap(id_thread)%grid2_add(sga_remap(id_thread)%num_links) = add2
3337	sga_remap(id_thread)%wts(1:num_wts,sga_remap(id_thread)%num_links) = weights(1:num_wts)
3338
3339	endif
3340
3341
3342	!-----------------------------------------------------------------------
3343
3344	end subroutine store_link_cnsrv
3345
3346	!***********************************************************************
3347
3348	#ifdef TREAT_OVERLAY
3349	!
3350	! Modified by A. Piacentini for OASIS in order to compute only the permutation
3351	! without modifying the input arrays and to take two arrays (lon and lat) as input
3352	! The comparison is made on the lon*(lat-pi) value which is well defined because
3353	! SCRIP has forced lon in [0,2pi] and lat in [-pih,pih]
3354	!
3355	! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
3356	! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino
3357	! Copyright (C) 2001 PWSCF group
3358	!
3359	! This file is distributed under the terms of the GNU General Public
3360	! License. See the file `LICENSE' in the root directory of the
3361	! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
3362	!
3363	! Adapted from flib/hpsort_eps
3364	!---------------------------------------------------------------------
3365	subroutine hpsort_eps (n, rlon, rlat, ind, eps)
3366	!---------------------------------------------------------------------
3367	! sort an array ra(1:n) into ascending order using heapsort algorithm,
3368	! and considering two elements being equal if their values differ
3369	! for less than "eps".
3370	! n is input, ra is replaced on output by its sorted rearrangement.
3371	! create an index table (ind) by making an exchange in the index array
3372	! whenever an exchange is made on the sorted data array (ra).
3373	! in case of equal values in the data array (ra) the values in the
3374	! index array (ind) are used to order the entries.
3375	! if on input ind(1) = 0 then indices are initialized in the routine,
3376	! if on input ind(1) != 0 then indices are assumed to have been
3377	! initialized before entering the routine and these
3378	! indices are carried around during the sorting process
3379	!
3380	! no work space needed !
3381	! free us from machine-dependent sorting-routines !
3382	!
3383	! adapted from Numerical Recipes pg. 329 (new edition)
3384	!
3385	use kinds_mod
3386	implicit none
3387	!-input/output variables
3388	integer, intent(in) :: n
3389	real(kind=dbl_kind), intent(in) :: eps
3390	integer :: ind (n)
3391	real(kind=dbl_kind) :: rlon (n)
3392	real(kind=dbl_kind) :: rlat (n)
3393	!-local variables
3394	integer :: i, ir, j, l, iind
3395	real(kind=dbl_kind) :: rrlon
3396	real(kind=dbl_kind) :: rrlat
3397	!
3398	! initialize index array
3399	if (ind (1) .eq.0) then
3400	do i = 1, n
3401	ind (i) = i
3402	enddo
3403	endif
3404	! nothing to order
3405	if (n.lt.2) return
3406	! initialize indices for hiring and retirement-promotion phase
3407	l = n / 2 + 1
3408
3409	ir = n
3410
3411	sorting: do
3412
3413	! still in hiring phase
3414	if ( l .gt. 1 ) then
3415	l = l - 1
3416	iind = ind (l)
3417	rrlon = rlon (iind)
3418	rrlat = rlat (iind)
3419	! in retirement-promotion phase.
3420	else
3421	! clear a space at the end of the array
3422	iind = ind (ir)
3423	rrlon = rlon (iind)
3424	rrlat = rlat (iind)
3425	!
3426	! retire the top of the heap into it
3427	!
3428	ind (ir) = ind (1)
3429	! decrease the size of the corporation
3430	ir = ir - 1
3431	! done with the last promotion
3432	if ( ir .eq. 1 ) then
3433	! the least competent worker at all !
3434	!
3435	ind (1) = iind
3436	exit sorting
3437	endif
3438	endif
3439	! wheter in hiring or promotion phase, we
3440	i = l
3441	! set up to place rra in its proper level
3442	j = l + l
3443	!
3444	do while ( j .le. ir )
3445	if ( j .lt. ir ) then
3446	! compare to better underling
3447	if ( hslt( rlon (ind(j)), rlat (ind(j)), &
3448	& rlon (ind(j + 1)), rlat (ind(j + 1)), eps ) ) then
3449	j = j + 1
3450	endif
3451	endif
3452	! demote rra
3453	if ( hslt( rrlon, rrlat, rlon (ind(j)), rlat (ind(j)), eps ) ) then
3454	ind (i) = ind (j)
3455	i = j
3456	j = j + j
3457	! this is the right place for rra
3458	else
3459	! set j to terminate do-while loop
3460	j = ir + 1
3461	endif
3462	enddo
3463	ind (i) = iind
3464
3465	end do sorting
3466	end subroutine hpsort_eps
3467
3468	! internal function
3469	! compare two coordinates and return the result
3470
3471	logical function hslt( lo1, la1, lo2, la2, eps )
3472	real(kind=dbl_kind), intent(in) :: lo1, la1, lo2, la2
3473	real(kind=dbl_kind), intent(in) :: eps
3474
3475	if (abs(la1-la2)<eps) then
3476	if (abs(lo1-lo2)<eps) then
3477	hslt = .false.
3478	else
3479	hslt = (lo1 < lo2)
3480	end if
3481	else
3482	hslt = (la1 < la2)
3483	end if
3484	end function hslt
3485
3486	!
3487
3488	#endif
3489	end module remap_conservative
3490
3491	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

Context Navigation

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

Download in other formats: