Context Navigation

source: CPL/oasis3-mct_5.0/lib/scrip/src/remap_bi_interp.F90 @ 6328

Last change on this file since 6328 was 6328, checked in by aclsce, 17 months ago
First import of oasis3-mct_5.0 (from oasis git server, branch OASIS3-MCT_5.0)
File size: 57.3 KB

Line
1	!****
2	! *****************************
3	! * OASIS MODULE - LEVEL ? *
4	! * ------------- ------- *
5	! *****************************
6	!
7	!**** remap_biear - calculate bilinear remapping
8	!
9	! Purpose:
10	! -------
11	! Adaptation of SCRIP 1.4 remap_biear module to calculate
12	! bilinear remapping.
13	!
14	!** Interface:
15	! ---------
16	! CALL remap_bi*
17	!
18	! Input:
19	! -----
20	!
21	! Output:
22	! ------
23	!
24	! History:
25	! -------
26	! Version Programmer Date Description
27	! ------- ---------- ---- -----------
28	! 2.5 S. Valcke 2002/09 Treatment for masked point
29	!
30	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32	!
33	! this module contains necessary routines for performing an
34	! bilinear interpolation.
35	!
36	!-----------------------------------------------------------------------
37	!
38	! CVS:$Id: remap_biear.f 2826 2010-12-10 11:14:21Z valcke $
39	!
40	! Copyright (c) 1997, 1998 the Regents of the University of
41	! California.
42	!
43	! This software and ancillary information (herein called software)
44	! called SCRIP is made available under the terms described here.
45	! The software has been approved for release with associated
46	! LA-CC Number 98-45.
47	!
48	! Unless otherwise indicated, this software has been authored
49	! by an employee or employees of the University of California,
50	! operator of the Los Alamos National Laboratory under Contract
51	! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S.
52	! Government has rights to use, reproduce, and distribute this
53	! software. The public may copy and use this software without
54	! charge, provided that this Notice and any statement of authorship
55	! are reproduced on all copies. Neither the Government nor the
56	! University makes any warranty, express or implied, or assumes
57	! any liability or responsibility for the use of this software.
58	!
59	! If software is modified to produce derivative works, such modified
60	! software should be clearly marked, so as not to confuse it with
61	! the version available from Los Alamos National Laboratory.
62	!
63	!***********************************************************************
64
65	module remap_bi_interp
66
67	!-----------------------------------------------------------------------
68
69	use kinds_mod ! defines common data types
70	use constants ! defines common constants
71	use grids ! module containing grid info
72	use remap_vars ! module containing remap info
73	use timers
74	USE mod_oasis_flush
75	!$ use omp_lib
76
77	implicit NONE
78
79	!-----------------------------------------------------------------------
80
81	real (kind=dbl_kind), dimension(:), allocatable, save :: &
82	coslat, sinlat, & ! cosine, sine of grid lats (for distance)
83	coslon, sinlon ! cosine, sine of grid lons (for distance)
84
85	integer (kind=int_kind) :: il_nbthreads = 1
86
87	integer (kind=int_kind), parameter :: max_iter = 100 ! max iteration count for i,j iteration
88
89	real (kind=dbl_kind), parameter :: converge = 1.e-10_dbl_kind ! convergence criterion
90
91	!***********************************************************************
92
93	contains
94
95	!***********************************************************************
96
97	subroutine remap_bi (lextrapdone, ll_nnei, &
98	mpi_comm_map, mpi_size_map, mpi_rank_map, mpi_root_map)
99
100	!-----------------------------------------------------------------------
101	!
102	! this routine computes the weights for a bilinear interpolation.
103	!
104	!-----------------------------------------------------------------------
105	!-----------------------------------------------------------------------
106	!
107	! input variables
108	!
109	!-----------------------------------------------------------------------
110
111	integer (kind=int_kind), INTENT(IN) :: mpi_comm_map, mpi_rank_map, mpi_size_map, mpi_root_map
112
113
114	LOGICAL, INTENT(IN) :: lextrapdone ! logical, true if EXTRAP done on field
115	LOGICAL, INTENT(IN) :: ll_nnei ! true (default) if extra search is done
116
117	!-----------------------------------------------------------------------
118	!
119	! local variables
120	!
121	!-----------------------------------------------------------------------
122
123	integer (kind=int_kind) :: n,icount, i, k, &
124	dst_add, & ! destination address
125	iter, & ! iteration counter
126	nmap ! index of current map being computed
127
128	integer (kind=int_kind), dimension(4) :: src_add ! address for the four source points
129
130	real (kind=dbl_kind), dimension(4) :: src_lats, & ! latitudes of four bilinear corners
131	src_lons ! longitudes of four bilinear corners
132
133	real (kind=dbl_kind), dimension(num_wts,4) :: wgts ! bilinear/bicubic weights for four corners
134	real (kind=dbl_kind) :: plat, plon, & ! lat/lon coords of destination point
135	iguess, jguess, & ! current guess for bilinear coordinate
136	deli, delj, & ! corrections to i,j
137	dth1, dth2, dth3, & ! some latitude differences
138	dph1, dph2, dph3, & ! some longitude differences
139	dthp, dphp, & ! difference between point and sw corner
140	mat1, mat2, mat3, mat4, & ! matrix elements
141	determinant, & ! matrix determinant
142	sum_wgts ! sum of weights for normalization
143
144	real (kind=dbl_kind) :: coslat_dst, sinlat_dst, coslon_dst, sinlon_dst, &
145	dist_min, distance, & ! for computing dist-weighted avg
146	src_latsnn, d_dist
147
148	real (kind=dbl_kind) :: dl_test
149
150	integer (kind=int_kind) :: min_add, max_add, srch_add, src_addnn, num_neighbors
151
152	integer (kind=int_kind) :: il_splitsize
153	integer (kind=int_kind) :: ib_proc
154	integer (kind=int_kind) :: ib_thread
155	integer (kind=int_kind) :: buff_base
156	integer (kind=int_kind), dimension(:), allocatable :: ila_mpi_sz
157	integer (kind=int_kind), dimension(:), allocatable :: ila_mpi_mn
158	integer (kind=int_kind), dimension(:), allocatable :: ila_mpi_mx
159	integer (kind=int_kind), dimension(:), allocatable :: ila_thr_sz
160	integer (kind=int_kind), dimension(:), allocatable :: ila_thr_mn
161	integer (kind=int_kind), dimension(:), allocatable :: ila_thr_mx
162	integer (kind=int_kind), dimension(:), allocatable :: ila_num_links_mpi
163	integer (kind=int_kind), dimension(:,:), allocatable :: ila_req_mpi
164	integer (kind=int_kind), dimension(:,:,:), allocatable :: ila_sta_mpi
165	character (LEN=14) :: cl_envvar
166	integer (kind=int_kind) :: il_envthreads, il_err
167	logical :: ll_gaussreduced_grid
168	logical :: ll_bicubic
169
170	#ifdef REMAP_TIMING
171	logical :: ll_timing = .true.
172	real (kind=dbl_kind), dimension(:,:), allocatable :: dla_timer
173	real (kind=dbl_kind) :: dl_tstart, dl_tstop
174	integer (kind=int_kind) :: il_mythread
175	#endif
176
177	!-----------------------------------------------------------------------
178	!
179	! compute mappings from grid1 to grid2
180	!
181	!-----------------------------------------------------------------------
182	!
183	IF (nlogprt .GE. 2) THEN
184	WRITE (UNIT = nulou,FMT = *)' '
185	WRITE (UNIT = nulou,FMT = *)'Entering routine remap_bi'
186	CALL OASIS_FLUSH_SCRIP(nulou)
187	ENDIF
188
189	#ifdef REMAP_TIMING
190	if (ll_timing) call timer_start(2,'remap_bi alloc')
191	#endif
192
193	! num_neighbors = 4 for bilinear or bicubic interpolation
194	num_neighbors = 4
195	nmap = 1
196	ll_gaussreduced_grid = .FALSE.
197	ll_bicubic = .FALSE.
198
199
200	IF (restrict_TYPE == 'REDUCED') ll_gaussreduced_grid = .TRUE.
201	IF (map_type == map_type_bicubic ) ll_bicubic = .TRUE.
202
203	if (grid1_rank /= 2) then
204	stop 'Can not do bilinear or bicubic interpolation when grid_rank /= 2'
205	endif
206
207	!***
208	!*** compute cos, sin of lat/lon on source grid for distance
209	!*** calculations
210	!***
211
212	allocate (coslat(grid1_size), coslon(grid1_size), sinlat(grid1_size), sinlon(grid1_size))
213
214	coslat = cos(grid1_center_lat)
215	coslon = cos(grid1_center_lon)
216	sinlat = sin(grid1_center_lat)
217	sinlon = sin(grid1_center_lon)
218
219	#ifdef REMAP_TIMING
220	if (ll_timing) call timer_stop(2)
221	#endif
222
223	!***
224	!*** precompute best scheduling of target grid points
225	!***
226
227	allocate(ila_mpi_mn(mpi_size_map))
228	allocate(ila_mpi_mx(mpi_size_map))
229
230	if (mpi_size_map .gt. 1) then
231
232	allocate(ila_mpi_sz(mpi_size_map))
233	il_splitsize = count(grid2_mask)
234	ila_mpi_sz(:) = floor(real(il_splitsize)/mpi_size_map)
235	ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) = ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) + 1
236
237	ila_mpi_mn(:) = 0
238	ib_proc = 1
239	il_splitsize = 0
240	do dst_add = 1, grid2_size
241	if (grid2_mask(dst_add)) then
242	if (ila_mpi_mn(ib_proc).eq.0) &
243	ila_mpi_mn(ib_proc) = dst_add
244	il_splitsize = il_splitsize + 1
245	if (il_splitsize .eq. ila_mpi_sz(ib_proc)) then
246	il_splitsize = 0
247	ila_mpi_mx(ib_proc) = dst_add
248	ib_proc = ib_proc + 1
249	end if
250	end if
251	end do
252
253	deallocate(ila_mpi_sz)
254
255	else
256
257	ila_mpi_mn(1) = 1
258	ila_mpi_mx(1) = grid2_size
259
260	endif
261
262
263	call get_environment_variable(name='OASIS_OMP_NUM_THREADS', &
264	& value=cl_envvar, status=il_err)
265	if ( il_err .ne. 0) then
266	il_envthreads = 0
267	else
268	read(cl_envvar,*) il_envthreads
269	end if
270
271	!$OMP PARALLEL NUM_THREADS(il_envthreads) DEFAULT(NONE) &
272	!$OMP SHARED(il_envthreads) &
273	!$OMP SHARED(lextrapdone,num_neighbors) &
274	!$OMP SHARED(grid2_mask,grid2_frac) &
275	!$OMP SHARED(grid2_center_lat,grid2_center_lon) &
276	!$OMP SHARED(grid1_center_lat,grid1_center_lon) &
277	!$OMP SHARED(grid1_bound_box,grid1_bbox_per) &
278	!$OMP SHARED(grid1_mask) &
279	!$OMP SHARED(bin_addr1,nulou,sga_remap,num_wts) &
280	!$OMP SHARED(nmap,dl_test) &
281	!$OMP PRIVATE(nlogprt) &
282	!$OMP SHARED(coslat,coslon,sinlat,sinlon) &
283	!$OMP PRIVATE(wgts,sum_wgts,dst_add,icount) &
284	!$OMP PRIVATE(coslat_dst,coslon_dst,sinlat_dst,sinlon_dst) &
285	!$OMP PRIVATE(plat,plon,src_add,src_lats,src_lons,min_add,max_add) &
286	!$OMP PRIVATE(src_addnn) &
287	!$OMP PRIVATE(dth1,dth2,dth3,dph1,dph2,dph3,iguess,jguess,iter,dthp,dphp) &
288	!$OMP PRIVATE(mat1,mat2,mat3,mat4,determinant,deli,delj,d_dist) &
289	!$OMP PRIVATE(distance,dist_min,n,srch_add,src_latsnn) &
290	!$OMP SHARED(il_nbthreads) &
291	!$OMP SHARED(mpi_rank_map,mpi_root_map,ila_mpi_mn,ila_mpi_mx) &
292	!$OMP SHARED(ila_thr_sz,ila_thr_mn,ila_thr_mx) &
293	!$OMP SHARED(grid1_dims,grid1_size) &
294	!$OMP SHARED(ll_nnei,ll_gaussreduced_grid,ll_bicubic) &
295	#ifdef REMAP_TIMING
296	!$OMP PRIVATE(ib_thread,il_splitsize) &
297	!$OMP SHARED(ll_timing,dla_timer) &
298	!$OMP PRIVATE(il_mythread,dl_tstart,dl_tstop)
299
300	!$ il_mythread = OMP_GET_THREAD_NUM () + 1
301	#else
302	!$OMP PRIVATE(ib_thread,il_splitsize)
303	#endif
304
305	!$OMP SINGLE
306
307	il_nbthreads = 1
308	!$ il_nbthreads = OMP_GET_NUM_THREADS ()
309
310	#ifdef REMAP_TIMING
311	if (ll_timing) then
312	if (il_nbthreads.gt.1) then
313	!$ dl_tstart = OMP_GET_WTIME()
314	else
315	call timer_start(3,'remap_bi distr')
316	end if
317	end if
318	#endif
319	allocate(ila_thr_mn(il_nbthreads))
320	allocate(ila_thr_mx(il_nbthreads))
321
322	if (il_nbthreads .gt. 1) then
323
324	#ifdef REMAP_TIMING
325	if (ll_timing) then
326	allocate(dla_timer(il_nbthreads,8))
327	dla_timer(:,:) = 0.0
328	end if
329	#endif
330	nlogprt = 0
331
332	allocate(ila_thr_sz(il_nbthreads))
333	il_splitsize = COUNT(grid2_mask(ila_mpi_mn(mpi_rank_map+1):&
334	& ila_mpi_mx(mpi_rank_map+1)))
335	ila_thr_sz(:) = floor(real(il_splitsize)/il_nbthreads)
336	ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) = ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) + 1
337
338	ila_thr_mn(:) = 0
339	ib_thread = 1
340	il_splitsize = 0
341	DO dst_add = ila_mpi_mn(mpi_rank_map+1), ila_mpi_mx(mpi_rank_map+1)
342	if (grid2_mask(dst_add)) then
343	if (ila_thr_mn(ib_thread).eq.0) &
344	ila_thr_mn(ib_thread) = dst_add
345	il_splitsize = il_splitsize + 1
346	if (il_splitsize .eq. ila_thr_sz(ib_thread)) then
347	il_splitsize = 0
348	ila_thr_mx(ib_thread) = dst_add
349	ib_thread = ib_thread + 1
350	end if
351	end if
352	end do
353
354	allocate(sga_remap(il_nbthreads))
355
356	do ib_thread = 1, il_nbthreads
357	il_splitsize = num_neighbors*ila_thr_sz(ib_thread)
358	sga_remap(ib_thread)%max_links = il_splitsize
359	sga_remap(ib_thread)%num_links = 0
360	allocate(sga_remap(ib_thread)%grid1_add(il_splitsize))
361	allocate(sga_remap(ib_thread)%grid2_add(il_splitsize))
362	allocate(sga_remap(ib_thread)%wts(num_wts,il_splitsize))
363	end do
364
365	deallocate(ila_thr_sz)
366
367	else
368
369	ila_thr_mn(1) = ila_mpi_mn(mpi_rank_map+1)
370	ila_thr_mx(1) = ila_mpi_mx(mpi_rank_map+1)
371
372	end if
373	#ifdef REMAP_TIMING
374	if (ll_timing) then
375	if (il_nbthreads.gt.1) then
376	!$ dl_tstop = OMP_GET_WTIME()
377	dla_timer(il_mythread,1)=dla_timer(il_mythread,1) + dl_tstop - dl_tstart
378	else
379	call timer_stop(3)
380	end if
381	end if
382	#endif
383	!$OMP END SINGLE
384
385
386	!***
387	!*** loop over destination grid
388	!***
389	!$OMP DO SCHEDULE(STATIC,1)
390	thread_loop: do ib_thread = 1, il_nbthreads
391
392	grid_loop1: do dst_add = ila_thr_mn(ib_thread), ila_thr_mx(ib_thread)
393
394	if (.not. grid2_mask(dst_add)) cycle grid_loop1
395
396	#ifdef REMAP_TIMING
397	if (ll_timing) then
398	if (il_nbthreads.gt.1) then
399	!$ dl_tstart = OMP_GET_WTIME()
400	else
401	call timer_start(4,'remap_bi search')
402	end if
403	end if
404	#endif
405
406	plat = grid2_center_lat(dst_add)
407	plon = grid2_center_lon(dst_add)
408
409	!***
410	!*** find nearest square of grid points on source grid
411	!***
412
413	IF ( ll_gaussreduced_grid ) then
414
415	call grid_search_bilin_reduced(src_add, src_lats, src_lons, &
416	min_add, max_add, &
417	plat, plon, grid1_dims, &
418	grid1_center_lat, grid1_center_lon, &
419	lextrapdone)
420	ELSE
421	call grid_search_bi(src_add, src_lats, src_lons, &
422	min_add, max_add, &
423	plat, plon, grid1_dims, &
424	grid1_center_lat, grid1_center_lon, &
425	grid1_bound_box, grid1_bbox_per, &
426	bin_addr1, lextrapdone)
427	ENDIF
428
429	#ifdef REMAP_TIMING
430	if (ll_timing) then
431	if (il_nbthreads.gt.1) then
432	!$ dl_tstop = OMP_GET_WTIME()
433	dla_timer(ib_thread,2) = dla_timer(ib_thread,2) + dl_tstop - dl_tstart
434	else
435	call timer_stop(4)
436	end if
437	end if
438	#endif
439
440	if (src_add(1) > 0) THEN
441
442	!***
443	!*** if the 4 surrounding points have been found and are
444	!*** non-masked, do bilinear interpolation
445	!***
446
447	#ifdef REMAP_TIMING
448	if (ll_timing) then
449	if (il_nbthreads.gt.1) then
450	!$ dl_tstart = OMP_GET_WTIME()
451	else
452	call timer_start(5,'remap_bi weights')
453	end if
454	end if
455	#endif
456
457
458	grid2_frac(dst_add) = one
459
460
461	dth1 = src_lats(2) - src_lats(1)
462	dth2 = src_lats(4) - src_lats(1)
463	dth3 = src_lats(3) - src_lats(2) - dth2
464
465	dph1 = src_lons(2) - src_lons(1)
466	dph2 = src_lons(4) - src_lons(1)
467	dph3 = src_lons(3) - src_lons(2)
468
469	if (dph1 > three*pih) dph1 = dph1 - pi2
470	if (dph2 > three*pih) dph2 = dph2 - pi2
471	if (dph3 > three*pih) dph3 = dph3 - pi2
472	if (dph1 < -three*pih) dph1 = dph1 + pi2
473	if (dph2 < -three*pih) dph2 = dph2 + pi2
474	if (dph3 < -three*pih) dph3 = dph3 + pi2
475
476	dph3 = dph3 - dph2
477
478	iguess = half
479	jguess = half
480
481	iter_loop1: do iter=1,max_iter
482
483	dthp = plat - src_lats(1) - dth1iguess - dth2jguess - dth3iguessjguess
484	dphp = plon - src_lons(1)
485
486	if (dphp > three*pih) dphp = dphp - pi2
487	if (dphp < -three*pih) dphp = dphp + pi2
488
489	dphp = dphp - dph1iguess - dph2jguess - dph3iguessjguess
490
491	mat1 = dth1 + dth3*jguess
492	mat2 = dth2 + dth3*iguess
493	mat3 = dph1 + dph3*jguess
494	mat4 = dph2 + dph3*iguess
495
496	determinant = mat1mat4 - mat2mat3
497
498	deli = (dthpmat4 - mat2dphp)/determinant
499	delj = (mat1dphp - dthpmat3)/determinant
500
501	if (abs(deli) < converge .and. abs(delj) < converge) &
502	exit iter_loop1
503
504	iguess = iguess + deli
505	jguess = jguess + delj
506
507	end do iter_loop1
508
509	#ifdef REMAP_TIMING
510	if (ll_timing) then
511	if (il_nbthreads.gt.1) then
512	!$ dl_tstop = OMP_GET_WTIME()
513	dla_timer(ib_thread,3) = dla_timer(ib_thread,3) + dl_tstop - dl_tstart
514	else
515	call timer_stop(5)
516	end if
517	end if
518	#endif
519
520	#ifdef REMAP_TIMING
521	if (ll_timing) then
522	if (il_nbthreads.gt.1) then
523	!$ dl_tstart = OMP_GET_WTIME()
524	else
525	call timer_start(6,'remap_dist_gaus_wgt store_link')
526	end if
527	end if
528	#endif
529
530	if (iter <= max_iter) then
531
532	!***
533	!*** successfully found i,j - compute weights
534	!***
535
536	if ( ll_bicubic ) then
537	wgts(1,1) = (one - jguess*2(three-twojguess)) (one - iguess*2(three-two*iguess))
538	wgts(1,2) = (one - jguess*2(three-twojguess)) iguess*2(three-two*iguess)
539	wgts(1,3) = jguess*2(three-twojguess) iguess*2(three-two*iguess)
540	wgts(1,4) = jguess*2(three-twojguess) (one - iguess*2(three-two*iguess))
541	wgts(2,1) = (one - jguess*2(three-twojguess)) iguess(iguess-one)*2
542	wgts(2,2) = (one - jguess*2(three-twojguess)) iguess*2(iguess-one)
543	wgts(2,3) = jguess*2(three-twojguess) iguess*2(iguess-one)
544	wgts(2,4) = jguess*2(three-twojguess) iguess(iguess-one)*2
545	wgts(3,1) = jguess(jguess-one)2 (one - iguess*2(three-two*iguess))
546	wgts(3,2) = jguess(jguess-one)2 iguess*2(three-two*iguess)
547	wgts(3,3) = jguess*2(jguess-one) * iguess*2(three-two*iguess)
548	wgts(3,4) = jguess*2(jguess-one) * (one - iguess*2(three-two*iguess))
549	wgts(4,1) = iguess(iguess-one)2 jguess(jguess-one)*2
550	wgts(4,2) = iguess*2(iguess-one) * jguess(jguess-one)*2
551	wgts(4,3) = iguess*2(iguess-one) * jguess*2(jguess-one)
552	wgts(4,4) = iguess(iguess-one)2 jguess*2(jguess-one)
553	! Bilinear and reduced bilinear cases
554	else
555	wgts(1,1) = (one-iguess)*(one-jguess)
556	wgts(1,2) = iguess*(one-jguess)
557	wgts(1,3) = iguess*jguess
558	wgts(1,4) = (one-iguess)*jguess
559	end if
560
561	call store_link_bi(dst_add, src_add, wgts, nmap, ib_thread)
562
563	#ifdef REMAP_TIMING
564	if (ll_timing) then
565	if (il_nbthreads.gt.1) then
566	!$ dl_tstop = OMP_GET_WTIME()
567	dla_timer(ib_thread,4) = dla_timer(ib_thread,4) + dl_tstop - dl_tstart
568	else
569	call timer_stop(6)
570	end if
571	end if
572	#endif
573
574
575	ELSE
576	write(nulou,*)'Point coords: ',plat,plon
577	do k = 1,size(src_add)
578	if (src_add(k) > 0) then
579	write(nulou,*)'Dest grid info: ',src_add(k),src_lats(k),src_lons(k),grid1_mask(src_add(k))
580	else
581	write(nulou,*)'Dest grid info: ',src_add(k),src_lats(k),src_lons(k)
582	endif
583	enddo
584	write(nulou,*)'Current i,j : ',iguess, jguess
585	write(nulou,*)'Iteration for i,j exceed max iteration count'
586	stop
587	endif
588	!
589	else if (src_add(1) < 0) THEN
590
591	!***
592	!*** Search for neighbour search failed or at least one of the 4
593	!*** neighbours was masked. Do distance-weighted average using
594	!*** the non-masked points among the 4 closest ones.
595	!***
596
597	IF (nlogprt .eq. 2) then
598	WRITE(nulou,*) ' '
599	WRITE(nulou,*) 'WARNING: Cannot make bilinear interpolation for target point',dst_add
600	WRITE(nulou,*) 'Using non-masked points among 4 nearest neighbors.'
601	WRITE(nulou,*) ' '
602	ENDIF
603
604	#ifdef REMAP_TIMING
605	if (ll_timing) then
606	if (il_nbthreads.gt.1) then
607	!$ dl_tstart = OMP_GET_WTIME()
608	else
609	call timer_start(7,'remap_bi nneighbour')
610	end if
611	end if
612	#endif
613	!***
614	!*** Find the 4 closest points
615	!***
616
617	coslat_dst = cos(plat)
618	sinlat_dst = sin(plat)
619	coslon_dst = cos(plon)
620	sinlon_dst = sin(plon)
621	src_add = 0
622	dist_min = bignum
623	src_lats = bignum
624
625	if ( ll_gaussreduced_grid ) then
626	IF (min_add == 0) min_add = 1
627	IF (max_add > grid1_size) max_add = grid1_size
628	endif
629
630	do srch_add = min_add,max_add
631	d_dist = coslat_dstcoslat(srch_add) &
632	(coslon_dst*coslon(srch_add) + &
633	sinlon_dst*sinlon(srch_add)) + &
634	sinlat_dst*sinlat(srch_add)
635	IF (d_dist < -1.0d0) THEN
636	d_dist = -1.0d0
637	ELSE IF (d_dist > 1.0d0) THEN
638	d_dist = 1.0d0
639	END IF
640	distance=acos(d_dist)
641
642	if (distance < dist_min) then
643	sort_loop: do n=1,4
644	if (distance < src_lats(n)) then
645	do i=4,n+1,-1
646	src_add (i) = src_add (i-1)
647	src_lats(i) = src_lats(i-1)
648	end do
649	src_add (n) = srch_add
650	src_lats(n) = distance
651	dist_min = src_lats(4)
652	exit sort_loop
653	endif
654	end do sort_loop
655	endif
656	end do
657
658	src_lons = one/(src_lats + tiny)
659	distance = sum(src_lons)
660	src_lats = src_lons/distance
661
662	!***
663	!*** Among 4 closest points, keep only the non-masked ones
664	!***
665
666	icount = 0
667	do n=1,4
668	if (grid1_mask(src_add(n)) .or. &
669	(.not. grid1_mask(src_add(n)) .and. lextrapdone)) then
670	icount = icount + 1
671	else
672	src_lats(n) = zero
673	endif
674	end do
675
676	#ifdef REMAP_TIMING
677	if (ll_timing) then
678	if (il_nbthreads.gt.1) then
679	!$ dl_tstop = OMP_GET_WTIME()
680	dla_timer(ib_thread,5) = dla_timer(ib_thread,5) + dl_tstop - dl_tstart
681	else
682	call timer_stop(7)
683	end if
684	end if
685	#endif
686
687	if (icount > 0) then
688
689	#ifdef REMAP_TIMING
690	if (ll_timing) then
691	if (il_nbthreads.gt.1) then
692	!$ dl_tstart = OMP_GET_WTIME()
693	else
694	call timer_start(6,'remap_dist_gaus_wgt store_link')
695	end if
696	end if
697	#endif
698
699	!*** renormalize weights
700	sum_wgts = sum(src_lats)
701	wgts(1,1) = src_lats(1)/sum_wgts
702	wgts(1,2) = src_lats(2)/sum_wgts
703	wgts(1,3) = src_lats(3)/sum_wgts
704	wgts(1,4) = src_lats(4)/sum_wgts
705	if ( ll_bicubic ) &
706	wgts(2:4,:) = 0.
707
708	grid2_frac(dst_add) = one
709	call store_link_bi(dst_add, src_add, wgts, nmap, ib_thread)
710	#ifdef REMAP_TIMING
711	if (ll_timing) then
712	if (il_nbthreads.gt.1) then
713	!$ dl_tstop = OMP_GET_WTIME()
714	dla_timer(ib_thread,4) = dla_timer(ib_thread,4) + dl_tstop - dl_tstart
715	else
716	call timer_stop(6)
717	end if
718	end if
719	#endif
720
721	ELSE
722
723	IF (ll_nnei .eqv. .true. ) then
724
725	#ifdef REMAP_TIMING
726	if (ll_timing) then
727	if (il_nbthreads.gt.1) then
728	!$ dl_tstart = OMP_GET_WTIME()
729	else
730	call timer_start(8,'remap_dist_gaus_wgt store_lk2')
731	end if
732	end if
733	#endif
734	IF (nlogprt .ge. 2) THEN
735	WRITE(nulou,*) ' '
736	WRITE(nulou,*) 'All 4 surrounding source grid points are masked'
737	WRITE(nulou,*) 'for target point ',dst_add
738	WRITE(nulou,*) 'with longitude and latitude', plon, plat
739	WRITE(nulou,*) 'Using the nearest non-masked neighbour.'
740	CALL OASIS_FLUSH_SCRIP(nulou)
741	ENDIF
742	!
743	src_latsnn = bignum
744	src_addnn = 0
745
746	!cdir novector
747	do srch_add = min_add,max_add
748	if (grid1_mask(srch_add) .or. &
749	(.not. grid1_mask(srch_add) .and. lextrapdone)) THEN
750	d_dist = coslat_dstcoslat(srch_add) &
751	(coslon_dst*coslon(srch_add) + &
752	sinlon_dst*sinlon(srch_add))+ &
753	sinlat_dst*sinlat(srch_add)
754	IF (d_dist < -1.0d0) THEN
755	d_dist = -1.0d0
756	ELSE IF (d_dist > 1.0d0) THEN
757	d_dist = 1.0d0
758	END IF
759	distance=acos(d_dist)
760	if (distance < src_latsnn) then
761	src_addnn = srch_add
762	src_latsnn = distance
763	endif
764	endif
765	end do
766	! Special treatment for bilinear reduced
767	! consequence of different bins definition
768	! in source and target grid
769	IF (src_addnn == 0 .AND. ll_gaussreduced_grid) THEN
770	src_latsnn = bignum
771	!cdir novector
772	do srch_add = 1, grid1_size
773	if (grid1_mask(srch_add) .or. lextrapdone) THEN
774	d_dist = coslat_dstcoslat(srch_add) &
775	(coslon_dst*coslon(srch_add) + &
776	sinlon_dst*sinlon(srch_add))+ &
777	sinlat_dst*sinlat(srch_add)
778	IF (d_dist < -1.0d0) THEN
779	d_dist = -1.0d0
780	ELSE IF (d_dist > 1.0d0) THEN
781	d_dist = 1.0d0
782	END IF
783	distance=acos(d_dist)
784	if (distance < src_latsnn) then
785	src_addnn = srch_add
786	src_latsnn = distance
787	endif
788	endif
789	end DO
790	ENDIF
791	IF (src_addnn == 0) THEN
792	WRITE(nulou,*) 'Problem with target grid point'
793	WRITE(nulou,*) 'with address = ',dst_add
794	call abort
795	ENDIF
796	IF (nlogprt .ge. 2) THEN
797	WRITE(nulou,*) ' '
798	WRITE(nulou,*) 'Nearest non masked neighbour is source point ',src_addnn
799	WRITE(nulou,*) 'with longitude and latitude',grid1_center_lon(src_addnn), &
800	grid1_center_lat(src_addnn)
801	ENDIF
802
803	wgts(1,1) = 1.
804	wgts(1,2:4) = 0.
805	if ( ll_bicubic ) &
806	wgts(2:4,:) = 0.
807
808	src_add(1) = src_addnn
809	src_add(2) = 0
810	src_add(3) = 0
811	src_add(4) = 0
812
813	grid2_frac(dst_add) = one
814	call store_link_bi(dst_add, src_add, wgts, nmap, ib_thread)
815	endif
816	#ifdef REMAP_TIMING
817	if (ll_timing) then
818	if (il_nbthreads.gt.1) then
819	!$ dl_tstop = OMP_GET_WTIME()
820	dla_timer(ib_thread,6) = dla_timer(ib_thread,6) + dl_tstop - dl_tstart
821	else
822	call timer_stop(8)
823	end if
824	end if
825	#endif
826	ENDIF
827	ENDIF
828	end do grid_loop1
829
830	end do thread_loop
831	!$OMP END DO
832
833	!$OMP END PARALLEL
834
835	if (il_nbthreads .gt. 1) then
836	#ifdef REMAP_TIMING
837	if (ll_timing) call timer_start(3,'remap_dist_gaus_wgt gather_lk')
838	#endif
839	sga_remap(1)%start_pos = 1
840	il_splitsize = sga_remap(1)%num_links
841	do ib_thread = 2, il_nbthreads
842	il_splitsize = il_splitsize + sga_remap(ib_thread)%num_links
843	sga_remap(ib_thread)%start_pos = sga_remap(ib_thread-1)%start_pos + &
844	sga_remap(ib_thread-1)%num_links
845	end do
846
847	num_links_map1 = il_splitsize
848	if (num_links_map1 > max_links_map1) &
849	call resize_remap_vars(1,num_links_map1-max_links_map1)
850
851	do ib_thread = 1, il_nbthreads
852	grid1_add_map1(sga_remap(ib_thread)%start_pos: &
853	sga_remap(ib_thread)%start_pos+ &
854	sga_remap(ib_thread)%num_links-1) = &
855	sga_remap(ib_thread)%grid1_add
856	grid2_add_map1(sga_remap(ib_thread)%start_pos: &
857	sga_remap(ib_thread)%start_pos+ &
858	sga_remap(ib_thread)%num_links-1) = &
859	sga_remap(ib_thread)%grid2_add
860	wts_map1 (:,sga_remap(ib_thread)%start_pos: &
861	sga_remap(ib_thread)%start_pos+ &
862	sga_remap(ib_thread)%num_links-1) = &
863	sga_remap(ib_thread)%wts
864
865	end do
866
867	#ifdef REMAP_TIMING
868	if (ll_timing) call timer_stop(3)
869	#endif
870	if (nlogprt.ge.2) then
871
872	do ib_thread = 1, il_nbthreads
873	if (sga_remap(ib_thread)%nb_resize.gt.0) then
874	IF (nlogprt .GE. 2) &
875	write(nulou,*) ' Number of thread_resize_remap_vars on thread ',&
876	ib_thread, ' = ', sga_remap(ib_thread)%nb_resize
877	end if
878	end do
879
880	end if
881	#ifdef REMAP_TIMING
882	if (ll_timing.and.nlogprt.ge.2) then
883	write(nulou,*) ' On master thread '
884	call timer_print(2)
885	call timer_clear(2)
886	do ib_thread = 1,il_nbthreads
887	write(nulou,*) ' On thread ',ib_thread
888	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
889	& 'remap_bi distr ',dla_timer(ib_thread,1)
890	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
891	& 'remap_bi search ',dla_timer(ib_thread,2)
892	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
893	& 'remap_bi weights ',dla_timer(ib_thread,3)
894	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
895	& 'remap_bi store_link',dla_timer(ib_thread,4)
896	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
897	& 'remap_bi nneighbour',dla_timer(ib_thread,5)
898	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
899	& 'remap_bi store_lk2',dla_timer(ib_thread,6)
900	end do
901	deallocate(dla_timer)
902	write(nulou,*) ' On master thread '
903	call timer_print(3)
904	call timer_clear(3)
905	end if
906
907	else
908
909	if (ll_timing.and.nlogprt.ge.2) then
910	do n = 2, 8
911	call timer_print(n)
912	call timer_clear(n)
913	end do
914	end if
915	#endif
916	end if
917
918	! Gather the complete results on master proc
919
920	if (mpi_size_map .gt. 1) then
921
922	IF (mpi_rank_map == mpi_root_map) THEN
923	ALLOCATE(ila_num_links_mpi(mpi_size_map))
924	ELSE
925	ALLOCATE(ila_num_links_mpi(1))
926	END IF
927
928	CALL MPI_Gather (num_links_map1, 1,MPI_INT,&
929	& ila_num_links_mpi,1,MPI_INT,&
930	& mpi_root_map,mpi_comm_map,il_err)
931
932
933	IF (mpi_rank_map == mpi_root_map) THEN
934	num_links_map1 = SUM(ila_num_links_mpi)
935	if (num_links_map1 > max_links_map1) &
936	call resize_remap_vars(1,num_links_map1-max_links_map1)
937
938	ALLOCATE(ila_req_mpi(4,mpi_size_map-1))
939	ALLOCATE(ila_sta_mpi(MPI_STATUS_SIZE,4,mpi_size_map-1))
940
941
942	DO n = 1, mpi_size_map-1
943	buff_base = SUM(ila_num_links_mpi(1:n))+1
944	CALL MPI_IRecv(grid1_add_map1(buff_base),&
945	& ila_num_links_mpi(n+1),MPI_INT,n,1,mpi_comm_map,&
946	& ila_req_mpi(1,n),il_err)
947
948	CALL MPI_IRecv(grid2_add_map1(buff_base),&
949	& ila_num_links_mpi(n+1),MPI_INT,n,2,mpi_comm_map,&
950	& ila_req_mpi(2,n),il_err)
951
952	CALL MPI_IRecv(wts_map1(:,buff_base),&
953	& num_wts*ila_num_links_mpi(n+1),MPI_DOUBLE,n,0,mpi_comm_map,&
954	& ila_req_mpi(3,n),il_err)
955
956	CALL MPI_IRecv(grid2_frac(ila_mpi_mn(n+1)),&
957	& ila_mpi_mx(n+1)-ila_mpi_mn(n+1)+1,MPI_DOUBLE,n,0,mpi_comm_map,&
958	& ila_req_mpi(4,n),il_err)
959
960	END DO
961
962	DO n=1,4
963	CALL MPI_Waitall(mpi_size_map-1,ila_req_mpi(n,:),ila_sta_mpi(1,n,1),il_err)
964	END DO
965
966	DEALLOCATE(ila_req_mpi)
967	DEALLOCATE(ila_sta_mpi)
968
969	ELSE
970
971	CALL MPI_Send(grid1_add_map1,num_links_map1,MPI_INT,&
972	& mpi_root_map,1,mpi_comm_map,il_err)
973
974	CALL MPI_Send(grid2_add_map1,num_links_map1,MPI_INT,&
975	& mpi_root_map,2,mpi_comm_map,il_err)
976
977	CALL MPI_Send(wts_map1,num_wts*num_links_map1,MPI_DOUBLE,&
978	& mpi_root_map,0,mpi_comm_map,il_err)
979
980	CALL MPI_Send(grid2_frac(ila_mpi_mn(mpi_rank_map+1)),&
981	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
982	& mpi_root_map,0,mpi_comm_map,il_err)
983
984	END IF
985
986	deallocate(ila_num_links_mpi)
987
988	end if
989
990
991	deallocate (coslat, coslon, sinlat, sinlon)
992	deallocate(ila_mpi_mn)
993	deallocate(ila_mpi_mx)
994	deallocate(ila_thr_mn)
995	deallocate(ila_thr_mx)
996	if (il_nbthreads .gt. 1) then
997	do ib_thread = 1, il_nbthreads
998	deallocate(sga_remap(ib_thread)%grid1_add)
999	deallocate(sga_remap(ib_thread)%grid2_add)
1000	deallocate(sga_remap(ib_thread)%wts)
1001	end do
1002	deallocate(sga_remap)
1003	end if
1004
1005	!
1006	IF (nlogprt .GE. 2) THEN
1007	WRITE (UNIT = nulou,FMT = *)' '
1008	WRITE (UNIT = nulou,FMT = *)'Leaving routine remap_bi'
1009	CALL OASIS_FLUSH_SCRIP(nulou)
1010	ENDIF
1011	!
1012	end subroutine remap_bi
1013
1014	!***********************************************************************
1015
1016	subroutine grid_search_bi(src_add, src_lats, src_lons, &
1017	min_add, max_add, &
1018	plat, plon, src_grid_dims, &
1019	src_center_lat, src_center_lon, &
1020	src_grid_bound_box, &
1021	src_grid_bbox_per, &
1022	src_bin_add, lextrapdone)
1023
1024	!-----------------------------------------------------------------------
1025	!
1026	! this routine finds the location of the search point plat, plon
1027	! in the source grid and returns the corners needed for a bilinear
1028	! interpolation.
1029	!
1030	!-----------------------------------------------------------------------
1031
1032	!-----------------------------------------------------------------------
1033	!
1034	! output variables
1035	!
1036	!-----------------------------------------------------------------------
1037
1038	integer (kind=int_kind), dimension(4), intent(out) :: src_add ! address of each corner point enclosing P
1039
1040	real (kind=dbl_kind), dimension(4), intent(out) :: src_lats, & ! latitudes of the four corner points
1041	src_lons ! longitudes of the four corner points
1042
1043	integer (kind=int_kind) :: min_add, max_add
1044
1045	!-----------------------------------------------------------------------
1046	!
1047	! input variables
1048	!
1049	!-----------------------------------------------------------------------
1050
1051	real (kind=dbl_kind), intent(in) :: plat, & ! latitude of the search point
1052	plon ! longitude of the search point
1053
1054	integer (kind=int_kind), dimension(2), intent(in) :: src_grid_dims ! size of each src grid dimension
1055
1056	real (kind=dbl_kind), dimension(:), intent(in) :: src_center_lat, & ! latitude of each src grid center
1057	src_center_lon ! longitude of each src grid center
1058
1059	real (kind=dbl_kind), dimension(:,:), intent(in) :: src_grid_bound_box ! bound box for source grid
1060
1061	integer (kind=int_kind), dimension(:), intent(in) :: src_grid_bbox_per ! bound box periodicity for source grid
1062
1063	integer (kind=int_kind), dimension(:,:), intent(in) :: src_bin_add ! latitude bins for restricting
1064
1065	LOGICAL :: lextrapdone ! logical, true if EXTRAP done on field
1066	!-----------------------------------------------------------------------
1067	!
1068	! local variables
1069	!
1070	!-----------------------------------------------------------------------
1071
1072	integer (kind=int_kind) :: n, next_n, srch_add, ni, & ! dummy indices
1073	nx, ny, ntotmask, & ! dimensions of src grid
1074	i, j, jp1, ip1, n_add, e_add, ne_add ! addresses
1075
1076	real (kind=dbl_kind) :: & ! vectors for cross-product check
1077	vec1_lat, vec1_lon, vec2_lat, vec2_lon, cross_product, cross_product_last, &
1078	dlat, dlon
1079
1080	!-----------------------------------------------------------------------
1081	!
1082	! restrict search first using bins
1083	!
1084	!-----------------------------------------------------------------------
1085
1086	src_add = 0
1087
1088	min_add = size(src_center_lat)
1089	max_add = 1
1090
1091	do n=1,num_srch_bins
1092	if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and. &
1093	plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then
1094	min_add = min(min_add, src_bin_add(1,n))
1095	max_add = max(max_add, src_bin_add(2,n))
1096	endif
1097	end do
1098
1099	!-----------------------------------------------------------------------
1100	!
1101	! now perform a more detailed search
1102	!
1103	!-----------------------------------------------------------------------
1104
1105	nx = src_grid_dims(1)
1106	ny = src_grid_dims(2)
1107
1108	srch_loop: do srch_add = min_add,max_add
1109
1110	!*** first check bounding box
1111
1112	if (plat <= src_grid_bound_box(2,srch_add) .and. &
1113	plat >= src_grid_bound_box(1,srch_add) .and. &
1114	lf_lon_in_box (plon,src_grid_bound_box(3,srch_add),&
1115	& src_grid_bound_box(4,srch_add),src_grid_bbox_per(srch_add)) ) then
1116
1117	!***
1118	!*** we are within bounding box so get really serious
1119	!***
1120
1121	!*** determine neighbor addresses
1122
1123	j = (srch_add - 1)/nx +1
1124	i = srch_add - (j-1)*nx
1125
1126	!*** find ip1
1127	!*** Note: I do not want to take into account the number
1128	!*** of overlapping grid points, as the underlying cell
1129	!*** will be found in all cases if the grid overlaps.
1130
1131	if (i < nx) then
1132	ip1 = i + 1
1133	else
1134	!*** assume cyclic
1135	ip1 = 1
1136	!*** but if it is not, correct
1137	dlat = abs(src_center_lat(srch_add) - src_center_lat(srch_add-1))
1138	dlon = abs(src_center_lon(srch_add) - src_center_lon(srch_add-1))
1139	e_add = (j - 1)*nx + ip1
1140	if (abs(src_center_lat(e_add ) - &
1141	src_center_lat(srch_add)) > 10.0 * dlat .or. &
1142	abs(src_center_lon(e_add) - &
1143	src_center_lon(srch_add)) > 10.0 * dlon) then
1144	ip1 = i
1145	endif
1146	endif
1147
1148	if (j < ny) then
1149	jp1 = j+1
1150	else
1151	!*** assume cyclic
1152	jp1 = 1
1153	!*** but if it is not, correct
1154	dlat = abs(src_center_lat(srch_add) - src_center_lat(srch_add-nx))
1155	dlon = abs(src_center_lon(srch_add) - src_center_lon(srch_add-nx))
1156	n_add = (jp1 - 1)*nx + i
1157	if (abs(src_center_lat(n_add) - &
1158	src_center_lat(srch_add)) > 10.0 * dlat .or. &
1159	abs(src_center_lon(n_add) - &
1160	src_center_lon(srch_add)) > 10.0 * dlon) then
1161	jp1 = j
1162	endif
1163	endif
1164
1165	n_add = (jp1 - 1)*nx + i
1166	e_add = (j - 1)*nx + ip1
1167	ne_add = (jp1 - 1)*nx + ip1
1168
1169	src_lats(1) = src_center_lat(srch_add)
1170	src_lats(2) = src_center_lat(e_add)
1171	src_lats(3) = src_center_lat(ne_add)
1172	src_lats(4) = src_center_lat(n_add)
1173
1174	src_lons(1) = src_center_lon(srch_add)
1175	src_lons(2) = src_center_lon(e_add)
1176	src_lons(3) = src_center_lon(ne_add)
1177	src_lons(4) = src_center_lon(n_add)
1178
1179	!***
1180	!*** for consistency, we must make sure all lons are in
1181	!*** same 2pi interval
1182	!***
1183
1184	vec1_lon = src_lons(1) - plon
1185	if (vec1_lon > pi) then
1186	src_lons(1) = src_lons(1) - pi2
1187	else if (vec1_lon < -pi) then
1188	src_lons(1) = src_lons(1) + pi2
1189	endif
1190	do n=2,4
1191	vec1_lon = src_lons(n) - src_lons(1)
1192	if (vec1_lon > pi) then
1193	src_lons(n) = src_lons(n) - pi2
1194	else if (vec1_lon < -pi) then
1195	src_lons(n) = src_lons(n) + pi2
1196	endif
1197	end do
1198
1199	corner_loop: do n=1,4
1200	next_n = MOD(n,4) + 1
1201
1202	!***
1203	!*** here we take the cross product of the vector making
1204	!*** up each box side with the vector formed by the vertex
1205	!*** and search point. if all the cross products are
1206	!*** positive, the point is contained in the box.
1207	!***
1208
1209	vec1_lat = src_lats(next_n) - src_lats(n)
1210	vec1_lon = src_lons(next_n) - src_lons(n)
1211	vec2_lat = plat - src_lats(n)
1212	vec2_lon = plon - src_lons(n)
1213
1214	!***
1215	!*** check for 0,2pi crossings
1216	!***
1217
1218	if (vec1_lon > three*pih) then
1219	vec1_lon = vec1_lon - pi2
1220	else if (vec1_lon < -three*pih) then
1221	vec1_lon = vec1_lon + pi2
1222	endif
1223	if (vec2_lon > three*pih) then
1224	vec2_lon = vec2_lon - pi2
1225	else if (vec2_lon < -three*pih) then
1226	vec2_lon = vec2_lon + pi2
1227	endif
1228
1229	cross_product = vec1_lonvec2_lat - vec2_lonvec1_lat
1230
1231	!***
1232	!*** if cross product is less than zero, this cell
1233	!*** doesn't work
1234	!***
1235
1236	if (n == 1) cross_product_last = cross_product
1237	if (cross_product*cross_product_last < zero) &
1238	exit corner_loop
1239	if (cross_product.ne.0) cross_product_last = cross_product
1240
1241	end do corner_loop
1242
1243	!***
1244	!*** if cross products all same sign, we found the location
1245	!***
1246
1247	if (n > 4) then
1248	src_add(1) = srch_add
1249	src_add(2) = e_add
1250	src_add(3) = ne_add
1251	src_add(4) = n_add
1252
1253	! Check if one point is masked; IF so, nearest-neighbour
1254	! interpolation will be used
1255
1256	ntotmask = 0
1257	do ni=1,4
1258	if (.not. grid1_mask(src_add(ni)).and. .not. lextrapdone) &
1259	ntotmask = ntotmask + 1
1260	end DO
1261	IF (ntotmask .gt. 0) src_add(1) = -src_add(1)
1262
1263	return
1264	endif
1265
1266	!***
1267	!*** otherwise move on to next cell
1268	!***
1269
1270	endif !bounding box check
1271	end do srch_loop
1272
1273	!***
1274	!*** if no cell found, point is likely either in a box that straddles
1275	!*** either pole or is outside the grid. Put src_add = -1 so that
1276	!*** distance-weighted average of the 4 non-masked closest points
1277	!*** is done in calling routine.
1278
1279	src_add = -1
1280
1281	!-----------------------------------------------------------------------
1282
1283	contains
1284
1285	logical (kind=log_kind) function lf_lon_in_box(rd_lon,&
1286	& rd_west,rd_east,id_per)
1287
1288	real (kind=dbl_kind), intent(in) :: rd_lon
1289	real (kind=dbl_kind), intent(in) :: rd_west,rd_east
1290	integer (kind=int_kind), intent(in) :: id_per
1291
1292	select case (id_per)
1293	case (0)
1294	lf_lon_in_box = &
1295	rd_lon <= rd_east .and. rd_lon >= rd_west
1296
1297	case (1)
1298	lf_lon_in_box = &
1299	rd_lon <= rd_east .or. rd_lon >= rd_west
1300
1301	end select
1302
1303	end function lf_lon_in_box
1304
1305	!-----------------------------------------------------------------------
1306
1307	end subroutine grid_search_bi
1308
1309	!***********************************************************************
1310
1311	subroutine store_link_bi(dst_add, src_add, weights, nmap, id_thread)
1312
1313	!-----------------------------------------------------------------------
1314	!
1315	! this routine stores the address and weight for four links
1316	! associated with one destination point in the appropriate address
1317	! and weight arrays and resizes those arrays if necessary.
1318	!
1319	!-----------------------------------------------------------------------
1320
1321	!-----------------------------------------------------------------------
1322	!
1323	! input variables
1324	!
1325	!-----------------------------------------------------------------------
1326
1327	integer (kind=int_kind), intent(in) :: dst_add, & ! address on destination grid
1328	nmap, & ! identifies which direction for mapping
1329	id_thread ! local threaded task
1330
1331	integer (kind=int_kind), dimension(4), intent(in) :: src_add ! addresses on source grid
1332
1333	real (kind=dbl_kind), dimension(num_wts,4), intent(in) :: weights ! array of remapping weights for these links
1334
1335	!-----------------------------------------------------------------------
1336	!
1337	! local variables
1338	!
1339	!-----------------------------------------------------------------------
1340
1341	integer (kind=int_kind) :: n, num_links_old ! placeholder for old link number
1342
1343	!-----------------------------------------------------------------------
1344	!
1345	! increment number of links and check to see if remap arrays need
1346	! to be increased to accomodate the new link. then store the
1347	! link.
1348	!
1349	!-----------------------------------------------------------------------
1350
1351	if (il_nbthreads .eq. 1) then
1352
1353	select case (nmap)
1354	case(1)
1355
1356	num_links_old = num_links_map1
1357	num_links_map1 = num_links_old + 4
1358
1359	if (num_links_map1 > max_links_map1) &
1360	call resize_remap_vars(1,resize_increment)
1361
1362	do n=1,4
1363	grid1_add_map1(num_links_old+n) = src_add(n)
1364	grid2_add_map1(num_links_old+n) = dst_add
1365	wts_map1 (1:num_wts,num_links_old+n) = weights(1:num_wts,n)
1366	end do
1367
1368	case(2)
1369
1370	num_links_old = num_links_map2
1371	num_links_map2 = num_links_old + 4
1372
1373	if (num_links_map2 > max_links_map2) &
1374	call resize_remap_vars(2,resize_increment)
1375
1376	do n=1,4
1377	grid1_add_map2(num_links_old+n) = dst_add
1378	grid2_add_map2(num_links_old+n) = src_add(n)
1379	wts_map2 (1:num_wts,num_links_old+n) = weights(1:num_wts,n)
1380	end do
1381
1382	end select
1383
1384	else
1385
1386	sga_remap(id_thread)%num_links = sga_remap(id_thread)%num_links + 4
1387
1388	if (sga_remap(id_thread)%num_links > sga_remap(id_thread)%max_links) &
1389	call sga_remap(id_thread)%resize(int(0.2*real(sga_remap(id_thread)%max_links)))
1390
1391	do n=1,4
1392	sga_remap(id_thread)%grid1_add(sga_remap(id_thread)%num_links-4+n) = src_add(n)
1393	sga_remap(id_thread)%grid2_add(sga_remap(id_thread)%num_links-4+n) = dst_add
1394	sga_remap(id_thread)%wts(1:num_wts,sga_remap(id_thread)%num_links-4+n) = weights(1:num_wts,n)
1395	end do
1396
1397	endif
1398
1399
1400	!-----------------------------------------------------------------------
1401
1402	end subroutine store_link_bi
1403
1404	!***********************************************************************
1405
1406	subroutine grid_search_bilin_reduced(src_add, src_lats, src_lons, &
1407	min_add, max_add, &
1408	plat, plon, src_grid_dims, &
1409	src_center_lat, src_center_lon, &
1410	lextrapdone)
1411
1412	!-----------------------------------------------------------------------
1413	!
1414	! this routine finds the location of the search point plat, plon
1415	! in the source grid and returns the corners needed for a bilinear
1416	! interpolation. The target grid is a reduced grid.
1417	!
1418	!-----------------------------------------------------------------------
1419
1420	!-----------------------------------------------------------------------
1421	!
1422	! output variables
1423	!
1424	!-----------------------------------------------------------------------
1425
1426	integer (kind=int_kind), dimension(4), intent(out) :: src_add ! address of each corner point enclosing P
1427
1428	real (kind=dbl_kind), dimension(4), intent(out) :: src_lats, & ! latitudes of the four corner points
1429	src_lons ! longitudes of the four corner points
1430
1431	integer (kind=int_kind) :: min_add, max_add
1432
1433	!-----------------------------------------------------------------------
1434	!
1435	! input variables
1436	!
1437	!-----------------------------------------------------------------------
1438
1439	real (kind=dbl_kind), intent(in) :: plat, & ! latitude of the search point
1440	plon ! longitude of the search point
1441
1442	integer (kind=int_kind), dimension(2), intent(in) :: src_grid_dims ! size of each src grid dimension
1443
1444	real (kind=dbl_kind), dimension(:), intent(in) :: src_center_lat, & ! latitude of each src grid center
1445	src_center_lon ! longitude of each src grid center
1446
1447	LOGICAL :: lextrapdone ! logical, true if EXTRAP done on field
1448
1449	!-----------------------------------------------------------------------
1450	!
1451	! local variables
1452	!
1453	!-----------------------------------------------------------------------
1454
1455	integer (kind=int_kind) :: n, next_n, srch_add, ni, &
1456	nx, ny, ntotmask, & ! dimensions of src grid
1457	inter_add ! add for restricting search
1458	!
1459	integer (kind=int_kind), DIMENSION(4) :: src_bid
1460
1461	!-----------------------------------------------------------------------
1462	!
1463	! restrict search first using bins
1464	!
1465	!-----------------------------------------------------------------------
1466
1467	nx = src_grid_dims(1)
1468	inter_add = 0
1469
1470	src_add = 0
1471
1472	min_add = size(src_center_lat) + 1
1473	max_add = 1
1474	if (plat >= bin_lats_r(1,1)) then
1475	min_add = 0
1476	max_add = bin_addr1_r(4,1)
1477	inter_add = bin_addr1_r(3,1)
1478	else if (plat <= bin_lats_r(1,num_srch_red)) then
1479	max_add = nx + 1
1480	min_add = bin_addr1_r(1,num_srch_red)
1481	inter_add = bin_addr1_r(3,num_srch_red)
1482	else
1483	srch_loopred: do n=1,num_srch_red
1484	if (plat <= bin_lats_r(1,n) .and. plat >= bin_lats_r(2,n)) then
1485	min_add = bin_addr1_r(1,n)
1486	max_add = bin_addr1_r(4,n)
1487	inter_add = bin_addr1_r(3,n)
1488	exit srch_loopred
1489	endif
1490	end DO srch_loopred
1491	ENDIF
1492
1493	!-----------------------------------------------------------------------
1494	!
1495	! now perform a more detailed search
1496	!
1497	!-----------------------------------------------------------------------
1498	if (min_add .ne. 0 .and. max_add .ne. nx+1) THEN
1499	!*** The concerned bins are not the top north or south ones.
1500	!*** We should be able to find the four corners
1501	!*** for the bilinear interpolation.
1502
1503	IF ( plon <= src_center_lon(min_add) ) THEN
1504	src_add(1) = inter_add-1
1505	src_add(2) = min_add
1506	ELSE IF ( plon > src_center_lon(inter_add-1) ) then
1507	src_add(1) = inter_add-1
1508	src_add(2) = min_add
1509	else
1510	srch_loop2: do srch_add = min_add, inter_add-2
1511	if ( (plon > src_center_lon(srch_add)) .and. &
1512	(plon <= src_center_lon(srch_add+1)) )then
1513	src_add(1) = srch_add
1514	src_add(2) = srch_add+1
1515	exit srch_loop2
1516	endif
1517	end do srch_loop2
1518	ENDIF
1519	IF ( plon <= src_center_lon(inter_add) ) THEN
1520	src_add(3) = max_add
1521	src_add(4) = inter_add
1522	ELSE IF ( plon >= src_center_lon(max_add) ) then
1523	src_add(3) = max_add
1524	src_add(4) = inter_add
1525	else
1526	srch_loop3: do srch_add = inter_add, max_add
1527	if ( (plon >= src_center_lon(srch_add)) .and. &
1528	(plon <= src_center_lon(srch_add+1)) )then
1529	src_add(3) = srch_add
1530	src_add(4) = srch_add+1
1531	exit srch_loop3
1532	endif
1533	enddo srch_loop3
1534	ENDIF
1535	src_lats(1) = src_center_lat(src_add(3))
1536	src_lats(2) = src_center_lat(src_add(4))
1537	src_lats(3) = src_center_lat(src_add(2))
1538	src_lats(4) = src_center_lat(src_add(1))
1539
1540	src_lons(1) = src_center_lon(src_add(3))
1541	src_lons(2) = src_center_lon(src_add(4))
1542	src_lons(3) = src_center_lon(src_add(2))
1543	src_lons(4) = src_center_lon(src_add(1))
1544
1545	src_bid=src_add
1546
1547	src_add(1) = src_bid(3)
1548	src_add(2) = src_bid(4)
1549	src_add(3) = src_bid(2)
1550	src_add(4) = src_bid(1)
1551
1552	! Check if one point is masked; IF so, nearest-neighbour
1553	! interpolation will be used
1554
1555	ntotmask = 0
1556	do ni=1,4
1557	if (.not. grid1_mask(src_add(ni)).and. .not. lextrapdone) &
1558	ntotmask = ntotmask + 1
1559	end DO
1560	IF (ntotmask .gt. 0) src_add(1) = -src_add(1)
1561
1562	ELSE
1563
1564	!*** We are in the first or last bin. Put src_add = -1 so that
1565	!*** distance-weighted average of the 4 non-masked closest points
1566	!*** is done in calling routine.
1567
1568	src_add = -1
1569
1570	ENDIF
1571
1572	!-----------------------------------------------------------------------
1573
1574	end subroutine grid_search_bilin_reduced
1575
1576	!***********************************************************************
1577
1578	end module remap_bi_interp
1579
1580	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

Download in other formats: