Context Navigation

source: CPL/oasis3-mct_5.0/lib/scrip/src/remap_dist_gaus_wgt.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: 44.5 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	! Modifications:
11	! - restrict types are written in capital letters
12	! - bug line 386: bin_lons(3,n) instead of bin_lons(2,n)
13	!
14	! Modified by V. Gayler, M&D 20.09.2001
15	! Modified by D. Declat, CERFACS 27.06.2002
16	! Modified by E. Maisonnave, CERFACS 21.12.2017
17	! Modified by A. Piacentini, CERFACS 21.12.2017
18	!***********************************************************************
19	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20	!
21	! this module contains necessary routines for performing an
22	! interpolation using a distance-weighted average of n nearest
23	! neighbors.
24	!
25	!-----------------------------------------------------------------------
26	!
27	! CVS:$Id: remap_dist_gaus_wgt.F 2826 2010-12-10 11:14:21Z valcke $
28	!
29	! Copyright (c) 1997, 1998 the Regents of the University of
30	! California.
31	!
32	! This software and ancillary information (herein called software)
33	! called SCRIP is made available under the terms described here.
34	! The software has been approved for release with associated
35	! LA-CC Number 98-45.
36	!
37	! Unless otherwise indicated, this software has been authored
38	! by an employee or employees of the University of California,
39	! operator of the Los Alamos National Laboratory under Contract
40	! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S.
41	! Government has rights to use, reproduce, and distribute this
42	! software. The public may copy and use this software without
43	! charge, provided that this Notice and any statement of authorship
44	! are reproduced on all copies. Neither the Government nor the
45	! University makes any warranty, express or implied, or assumes
46	! any liability or responsibility for the use of this software.
47	!
48	! If software is modified to produce derivative works, such modified
49	! software should be clearly marked, so as not to confuse it with
50	! the version available from Los Alamos National Laboratory.
51	!
52	!***********************************************************************
53
54	module remap_distance_gaussian_weight
55
56	!-----------------------------------------------------------------------
57
58	use kinds_mod ! defines common data types
59	use constants ! defines common constants
60	use grids ! module containing grid info
61	use remap_vars ! module containing remap info
62	use timers
63	use mod_oasis_flush
64	!$ use omp_lib
65
66	implicit none
67
68	!-----------------------------------------------------------------------
69	!
70	! module variables
71	!
72	!-----------------------------------------------------------------------
73
74	real (kind=dbl_kind), dimension(:), allocatable, save :: &
75	coslat, sinlat, & ! cosine, sine of grid lats (for distance)
76	coslon, sinlon, & ! cosine, sine of grid lons (for distance)
77	wgtstmp ! an array to hold the link weight
78
79	logical :: ll_nnei_dgw=.true.
80
81	integer (kind=int_kind) :: il_nbthreads = 1
82
83	!***********************************************************************
84
85	contains
86
87	!***********************************************************************
88
89	subroutine remap_dist_gaus_wgt (lextrapdone, ll_nnei_in, num_neighbors, &
90	mpi_comm_map, mpi_size_map, mpi_rank_map, mpi_root_map, &
91	r_varmul )
92
93	!-----------------------------------------------------------------------
94	!
95	! this routine computes the inverse-distance weights for a
96	! nearest-neighbor interpolation.
97	!
98	!-----------------------------------------------------------------------
99	!-----------------------------------------------------------------------
100	!
101	! input variables
102	!
103	!-----------------------------------------------------------------------
104
105	logical, INTENT(in) :: lextrapdone ! logical, true if EXTRAP done on field
106
107	logical, INTENT(in) :: ll_nnei_in ! logical if nearest neighbor fill to be done
108
109	real (kind=dbl_kind), INTENT(in), optional :: r_varmul ! Gaussian variance
110
111	integer (kind=int_kind), INTENT(in) :: num_neighbors ! number of neighbours
112
113	integer (kind=int_kind), INTENT(in) :: mpi_comm_map, mpi_rank_map, mpi_size_map, mpi_root_map
114
115	!-----------------------------------------------------------------------
116	!
117	! local variables
118	!
119	!-----------------------------------------------------------------------
120
121	logical (kind=log_kind), dimension(num_neighbors) :: nbr_mask ! mask at nearest neighbors
122
123	logical (kind=log_kind) :: ll_debug ! for debug outputs
124
125	integer (kind=int_kind) :: i_n, &
126	dst_add, & ! destination address
127	nmap, & ! index of current map being computed
128	icount ! number of non masked nearest neighbour
129
130	integer (kind=int_kind), dimension(num_neighbors) :: nbr_add ! source address at nearest neighbors
131
132	real (kind=dbl_kind), dimension(num_neighbors) :: nbr_dist ! angular distance four nearest neighbors
133
134	real (kind=dbl_kind) :: coslat_dst, & ! cos(lat) of destination grid point
135	coslon_dst, & ! cos(lon) of destination grid point
136	sinlat_dst, & ! sin(lat) of destination grid point
137	sinlon_dst, & ! sin(lon) of destination grid point
138	dist_tot, & ! sum of neighbor distances (for normalizing)
139	dist_average ! average distance between the source points
140
141	real (kind=dbl_kind) :: dl_test
142
143	real (kind=dbl_kind) :: distance, plat, plon, src_latsnn, dl_dist ! angular distance
144	real (kind=dbl_kind), dimension (1) :: wgts_new
145	real (kind=dbl_kind) :: rl_varmul ! local Gaussian variance
146	real (kind=dbl_kind) :: rl_local_variance ! local mean variance
147	real (kind=dbl_kind) :: rl_nb_distance ! local neighbour distance
148
149
150	integer (kind=int_kind) :: src_addnn, il_debug_add
151
152	integer (kind=int_kind) :: il_splitsize
153	integer (kind=int_kind) :: ib_proc
154	integer (kind=int_kind) :: ib_thread, ib_nb
155	integer (kind=int_kind) :: il_nb_pairs
156	integer (kind=int_kind) :: buff_base
157	integer (kind=int_kind), dimension(:), allocatable :: ila_mpi_sz
158	integer (kind=int_kind), dimension(:), allocatable :: ila_mpi_mn
159	integer (kind=int_kind), dimension(:), allocatable :: ila_mpi_mx
160	integer (kind=int_kind), dimension(:), allocatable :: ila_thr_sz
161	integer (kind=int_kind), dimension(:), allocatable :: ila_thr_mn
162	integer (kind=int_kind), dimension(:), allocatable :: ila_thr_mx
163	integer (kind=int_kind), dimension(:), allocatable :: ila_num_links_mpi
164	integer (kind=int_kind), dimension(:,:), allocatable :: ila_req_mpi
165	integer (kind=int_kind), dimension(:,:,:), allocatable :: ila_sta_mpi
166	character (LEN=14) :: cl_envvar
167	integer (kind=int_kind) :: il_envthreads, il_err, il_n1_add, il_n2_add
168	logical :: ll_gaussian_weights
169	#ifdef REMAP_TIMING
170	logical :: ll_timing = .true.
171	real (kind=8), dimension(:,:), allocatable :: dla_timer
172	real (kind=8) :: dl_tstart, dl_tstop
173	integer (kind=int_kind) :: il_mythread
174	#endif
175
176	!-----------------------------------------------------------------------
177	!
178	if (nlogprt .ge. 2) then
179	write (UNIT = nulou,FMT = *)' '
180	write (UNIT = nulou,FMT = *)'Entering routine remap_dist_gaus_wgt'
181	call OASIS_FLUSH_SCRIP(nulou)
182	endif
183
184	! copy to module data
185	ll_nnei_dgw = ll_nnei_in
186
187	!
188	!-----------------------------------------------------------------------
189	!
190	! compute mappings from grid1 to grid2
191	!
192	!-----------------------------------------------------------------------
193	#ifdef REMAP_TIMING
194	if (ll_timing) call timer_start(2,'remap_dist_gaus_wgt alloc-res_ave')
195	#endif
196
197	if (PRESENT(r_varmul )) then
198
199	! Check that variance is not zero
200	IF ( r_varmul < epsilon(1.) ) then
201	write (UNIT = nulou,FMT = ) '** ERROR ***'
202	write (UNIT = nulou,FMT = *) 'Namcouple GAUSWGT variance $VAR cannot be zero'
203	call OASIS_FLUSH_SCRIP(nulou)
204	stop
205	ENDIF
206
207	! if num_neighbors = 1, GAUSWGT not applicable
208	IF ( num_neighbors < 2 ) then
209	write (UNIT = nulou,FMT = ) '** WARNING ***'
210	write (UNIT = nulou,FMT = *) 'For GAUSWGT, $NV has to be bigger than 1'
211	write (UNIT = nulou,FMT = *) 'DISTWGT performed instead of GAUSWGT'
212	call OASIS_FLUSH_SCRIP(nulou)
213
214	ll_gaussian_weights = .false.
215	rl_varmul = 0.
216	ELSE
217	ll_gaussian_weights = .true.
218	rl_varmul = r_varmul
219	ENDIF
220
221	else
222
223	ll_gaussian_weights = .false.
224	rl_varmul = 0.
225
226	endif
227
228	dl_test = 0.0
229	nmap = 1
230
231	!***
232	!*** allocate wgtstmp to be consistent with store_link interface
233	!***
234
235	allocate (wgtstmp(num_wts))
236
237	!***
238	!*** compute cos, sin of lat/lon on source grid for distance
239	!*** calculations
240	!***
241
242	allocate (coslat(grid1_size), coslon(grid1_size), sinlat(grid1_size), sinlon(grid1_size))
243	coslat = cos(grid1_center_lat)
244	coslon = cos(grid1_center_lon)
245	sinlat = sin(grid1_center_lat)
246	sinlon = sin(grid1_center_lon)
247
248	#ifdef REMAP_TIMING
249	if (ll_timing) call timer_stop(2)
250	#endif
251
252	!***
253	!*** precompute best scheduling of target grid points
254	!***
255
256	allocate(ila_mpi_mn(mpi_size_map))
257	allocate(ila_mpi_mx(mpi_size_map))
258
259	if (mpi_size_map .gt. 1) then
260
261	allocate(ila_mpi_sz(mpi_size_map))
262	il_splitsize = count(grid2_mask)
263	ila_mpi_sz(:) = floor(real(il_splitsize)/mpi_size_map)
264	ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) = ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) + 1
265
266	ila_mpi_mn(:) = 0
267	ib_proc = 1
268	il_splitsize = 0
269	do dst_add = 1, grid2_size
270	if (grid2_mask(dst_add)) then
271	if (ila_mpi_mn(ib_proc).eq.0) &
272	ila_mpi_mn(ib_proc) = dst_add
273	il_splitsize = il_splitsize + 1
274	if (il_splitsize .eq. ila_mpi_sz(ib_proc)) then
275	il_splitsize = 0
276	ila_mpi_mx(ib_proc) = dst_add
277	ib_proc = ib_proc + 1
278	end if
279	end if
280	end do
281
282	deallocate(ila_mpi_sz)
283
284	else
285
286	ila_mpi_mn(1) = 1
287	ila_mpi_mx(1) = grid2_size
288
289	endif
290
291
292	call get_environment_variable(name='OASIS_OMP_NUM_THREADS', &
293	& value=cl_envvar, status=il_err)
294	if ( il_err .ne. 0) then
295	il_envthreads = 0
296	else
297	read(cl_envvar,*) il_envthreads
298	end if
299
300	!$OMP PARALLEL NUM_THREADS(il_envthreads) DEFAULT(NONE) &
301	!$OMP SHARED(il_envthreads) &
302	!$OMP SHARED(lextrapdone,ll_nnei_dgw,num_neighbors) &
303	!$OMP SHARED(grid2_mask,grid2_frac) &
304	!$OMP SHARED(grid2_center_lat,grid2_center_lon) &
305	!$OMP SHARED(grid1_mask) &
306	!$OMP SHARED(bin_addr1,nulou,sga_remap,num_wts) &
307	!$OMP SHARED(nmap,dl_test) &
308	!$OMP PRIVATE(nlogprt) &
309	!$OMP SHARED(coslat,coslon,sinlat,sinlon) &
310	!$OMP PRIVATE(nbr_mask,wgtstmp,dst_add,icount,nbr_add,nbr_dist) &
311	!$OMP PRIVATE(coslat_dst,coslon_dst,sinlat_dst,sinlon_dst,dist_tot) &
312	!$OMP PRIVATE(ll_debug,src_addnn, il_debug_add) &
313	!$OMP SHARED(il_nbthreads) &
314	!$OMP SHARED(mpi_rank_map,mpi_root_map,ila_mpi_mn,ila_mpi_mx) &
315	!$OMP SHARED(ila_thr_sz,ila_thr_mn,ila_thr_mx) &
316	!$OMP SHARED(ll_gaussian_weights,rl_varmul) &
317	!$OMP PRIVATE(rl_local_variance,ib_nb,i_n,rl_nb_distance,il_nb_pairs) &
318	!$OMP PRIVATE(il_n1_add, il_n2_add) &
319	#ifdef REMAP_TIMING
320	!$OMP PRIVATE(ib_thread,il_splitsize) &
321	!$OMP SHARED(ll_timing,dla_timer) &
322	!$OMP PRIVATE(il_mythread,dl_tstart,dl_tstop)
323
324	!$ il_mythread = OMP_GET_THREAD_NUM () + 1
325	#else
326	!$OMP PRIVATE(ib_thread,il_splitsize)
327	#endif
328
329	!$OMP SINGLE
330
331	il_nbthreads = 1
332	!$ il_nbthreads = OMP_GET_NUM_THREADS ()
333
334	#ifdef REMAP_TIMING
335	if (ll_timing) then
336	if (il_nbthreads.gt.1) then
337	!$ dl_tstart = OMP_GET_WTIME()
338	else
339	call timer_start(3,'remap_dist_gaus_wgt distr')
340	end if
341	end if
342	#endif
343	allocate(ila_thr_mn(il_nbthreads))
344	allocate(ila_thr_mx(il_nbthreads))
345
346	if (il_nbthreads .gt. 1) then
347
348	#ifdef REMAP_TIMING
349	if (ll_timing) then
350	allocate(dla_timer(il_nbthreads,4))
351	dla_timer(:,:) = 0.0
352	end if
353	#endif
354	nlogprt = 0
355
356	allocate(ila_thr_sz(il_nbthreads))
357	il_splitsize = COUNT(grid2_mask(ila_mpi_mn(mpi_rank_map+1):&
358	& ila_mpi_mx(mpi_rank_map+1)))
359	ila_thr_sz(:) = floor(real(il_splitsize)/il_nbthreads)
360	ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) = ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) + 1
361
362	ila_thr_mn(:) = 0
363	ib_thread = 1
364	il_splitsize = 0
365	DO dst_add = ila_mpi_mn(mpi_rank_map+1), ila_mpi_mx(mpi_rank_map+1)
366	if (grid2_mask(dst_add)) then
367	if (ila_thr_mn(ib_thread).eq.0) &
368	ila_thr_mn(ib_thread) = dst_add
369	il_splitsize = il_splitsize + 1
370	if (il_splitsize .eq. ila_thr_sz(ib_thread)) then
371	il_splitsize = 0
372	ila_thr_mx(ib_thread) = dst_add
373	ib_thread = ib_thread + 1
374	end if
375	end if
376	end do
377
378	allocate(sga_remap(il_nbthreads))
379
380	do ib_thread = 1, il_nbthreads
381	il_splitsize = num_neighbors*ila_thr_sz(ib_thread)
382	sga_remap(ib_thread)%max_links = il_splitsize
383	sga_remap(ib_thread)%num_links = 0
384	allocate(sga_remap(ib_thread)%grid1_add(il_splitsize))
385	allocate(sga_remap(ib_thread)%grid2_add(il_splitsize))
386	allocate(sga_remap(ib_thread)%wts(num_wts,il_splitsize))
387	end do
388
389	deallocate(ila_thr_sz)
390
391	else
392
393	ila_thr_mn(1) = ila_mpi_mn(mpi_rank_map+1)
394	ila_thr_mx(1) = ila_mpi_mx(mpi_rank_map+1)
395
396	end if
397	#ifdef REMAP_TIMING
398	if (ll_timing) then
399	if (il_nbthreads.gt.1) then
400	!$ dl_tstop = OMP_GET_WTIME()
401	dla_timer(il_mythread,1)=dla_timer(il_mythread,1) + dl_tstop - dl_tstart
402	else
403	call timer_stop(3)
404	end if
405	end if
406	#endif
407	!$OMP END SINGLE
408
409
410	!***
411	!*** loop over destination grid
412	!***
413	!$OMP DO SCHEDULE(STATIC,1)
414	thread_loop: do ib_thread = 1, il_nbthreads
415
416	grid_loop1: do dst_add =ila_thr_mn(ib_thread), ila_thr_mx(ib_thread)
417
418	if (.not. grid2_mask(dst_add)) cycle grid_loop1
419
420	#ifdef REMAP_TIMING
421	if (ll_timing) then
422	if (il_nbthreads.gt.1) then
423	!$ dl_tstart = OMP_GET_WTIME()
424	else
425	call timer_start(4,'remap_dist_gaus_wgt search')
426	end if
427	end if
428	#endif
429	coslat_dst = cos(grid2_center_lat(dst_add))
430	coslon_dst = cos(grid2_center_lon(dst_add))
431	sinlat_dst = sin(grid2_center_lat(dst_add))
432	sinlon_dst = sin(grid2_center_lon(dst_add))
433
434	!***
435	!*** find nearest grid points on source grid
436	!*** or non masked nearest neighbour
437	!*** and distances to each point
438	!***
439
440	call grid_search_nbr(src_addnn, nbr_add, nbr_dist, &
441	grid2_center_lat(dst_add), &
442	grid2_center_lon(dst_add), &
443	coslat_dst, coslon_dst, sinlat_dst, sinlon_dst, &
444	bin_addr1, num_neighbors, lextrapdone, dst_add)
445
446	ll_debug = .false.
447	if (ll_debug) then
448	il_debug_add = 15248
449	if (dst_add == il_debug_add) then
450	write(nulou,*) 'nbr_add = ', nbr_add(:)
451	write(nulou,*) 'nbr_dist = ', nbr_dist(:)
452	call OASIS_FLUSH_SCRIP(nulou)
453	endif
454	endif
455
456	#ifdef REMAP_TIMING
457	if (ll_timing) then
458	if (il_nbthreads.gt.1) then
459	!$ dl_tstop = OMP_GET_WTIME()
460	dla_timer(ib_thread,2) = dla_timer(ib_thread,2) + dl_tstop - dl_tstart
461	else
462	call timer_stop(4)
463	end if
464	end if
465	#endif
466	!***
467	!*** compute weights based on inverse distance
468	!*** if mask is false, eliminate those points
469	!***
470
471	#ifdef REMAP_TIMING
472	if (ll_timing) then
473	if (il_nbthreads.gt.1) then
474	!$ dl_tstart = OMP_GET_WTIME()
475	else
476	call timer_start(5,'remap_dist_gaus_wgt weights')
477	end if
478	end if
479	#endif
480	! now , if gaussian weights, compute local variance in parallel
481	if ( ll_gaussian_weights ) THEN
482	rl_local_variance = 0.
483	il_nb_pairs = 0
484	do i_n = 1, num_neighbors-1
485	do ib_nb = i_n+1, num_neighbors
486	il_n1_add = nbr_add(i_n)
487	il_n2_add = nbr_add(ib_nb)
488	if ( grid1_mask(il_n1_add) .and. grid1_mask(il_n2_add) ) then
489	rl_nb_distance = sinlat(il_n2_add)*sinlat(il_n1_add) + &
490	coslat(il_n2_add)coslat(il_n1_add) &
491	( coslon(il_n2_add)*coslon(il_n1_add) + &
492	sinlon(il_n2_add)*sinlon(il_n1_add) )
493	rl_nb_distance = MAX ( MIN ( rl_nb_distance, 1.d0), -1.d0)
494	rl_nb_distance = acos ( rl_nb_distance )
495	rl_local_variance = rl_local_variance + rl_nb_distance
496	il_nb_pairs = il_nb_pairs + 1
497	end if
498	enddo
499	enddo
500	! if less than two neighbours, DISTWGT with one neighbour will be performed
501	! if not , average distance and calculate variance
502	if ( il_nb_pairs > 0 ) &
503	rl_local_variance = rl_varmul * ( rl_local_variance / REAL(il_nb_pairs) ) ** 2
504	end if
505
506	dist_tot = zero
507	do i_n=1,num_neighbors
508	if (grid1_mask(nbr_add(i_n)) .or. lextrapdone) then
509	!
510	! Change distance if gausswgt option and more than one neighbour
511	if ( ll_gaussian_weights .and. il_nb_pairs > 0 ) &
512	nbr_dist(i_n) = exp(.5nbr_dist(i_n)nbr_dist(i_n)/rl_local_variance)
513	nbr_dist(i_n) = one/(nbr_dist(i_n) + epsilon(dl_test))
514	dist_tot = dist_tot + nbr_dist(i_n)
515	nbr_mask(i_n) = .true.
516	else
517	nbr_mask(i_n) = .false.
518	endif
519	end do
520
521	if (ll_debug) then
522	if (dst_add == il_debug_add) then
523	write(nulou,*) 'nbr_mask = ', nbr_mask(:)
524	write(nulou,*) 'nbr_dist = ', nbr_dist(:)
525	write(nulou,*) 'dist_tot = ', dist_tot
526	endif
527	endif
528
529	#ifdef REMAP_TIMING
530	if (ll_timing) then
531	if (il_nbthreads.gt.1) then
532	!$ dl_tstop = OMP_GET_WTIME()
533	dla_timer(ib_thread,3) = dla_timer(ib_thread,3) + dl_tstop - dl_tstart
534	else
535	call timer_stop(5)
536	end if
537	end if
538	#endif
539	!***
540	!*** normalize weights and store the link
541	!***
542
543	#ifdef REMAP_TIMING
544	if (ll_timing) then
545	if (il_nbthreads.gt.1) then
546	!$ dl_tstart = OMP_GET_WTIME()
547	else
548	call timer_start(6,'remap_dist_gaus_wgt store_link')
549	end if
550	end if
551	#endif
552	icount = 0
553	do i_n=1,num_neighbors
554	if (nbr_mask(i_n)) then
555	wgtstmp(1) = nbr_dist(i_n)/dist_tot
556	if (ll_debug) then
557	if (dst_add == il_debug_add) then
558	write(nulou,*) 'wgtstmp = ', wgtstmp(1)
559	write(nulou,*) 'nbr_dist = ', nbr_dist(:)
560	endif
561	endif
562	call store_link_nbr(nbr_add(i_n), dst_add, wgtstmp, nmap, ib_thread)
563	grid2_frac(dst_add) = one
564	icount = icount + 1
565	endif
566	end do
567	if (ll_nnei_dgw) then
568	if (icount == 0) then
569	if (nlogprt .ge. 2) then
570	write(nulou,*) ' '
571	write(nulou,*) 'Using the nearest non-masked neighbour because all 4 &
572	& surrounding source points are masked for target point ',dst_add
573	call OASIS_FLUSH_SCRIP(nulou)
574	endif
575	wgtstmp(1) = 1.
576	grid2_frac(dst_add) = one
577	call store_link_nbr(src_addnn, dst_add, wgtstmp, nmap, ib_thread)
578	endif
579	endif
580	#ifdef REMAP_TIMING
581	if (ll_timing) then
582	if (il_nbthreads.gt.1) then
583	!$ dl_tstop = OMP_GET_WTIME()
584	dla_timer(ib_thread,4) = dla_timer(ib_thread,4) + dl_tstop - dl_tstart
585	else
586	call timer_stop(6)
587	end if
588	end if
589	#endif
590	end do grid_loop1
591
592	end do thread_loop
593	!$OMP END DO
594
595	!$OMP END PARALLEL
596
597	if (il_nbthreads .gt. 1) then
598	#ifdef REMAP_TIMING
599	if (ll_timing) call timer_start(3,'remap_dist_gaus_wgt gather_lk')
600	#endif
601	sga_remap(1)%start_pos = 1
602	il_splitsize = sga_remap(1)%num_links
603	do ib_thread = 2, il_nbthreads
604	il_splitsize = il_splitsize + sga_remap(ib_thread)%num_links
605	sga_remap(ib_thread)%start_pos = sga_remap(ib_thread-1)%start_pos + &
606	sga_remap(ib_thread-1)%num_links
607	end do
608
609	num_links_map1 = il_splitsize
610	if (num_links_map1 > max_links_map1) &
611	call resize_remap_vars(1,num_links_map1-max_links_map1)
612
613	do ib_thread = 1, il_nbthreads
614	grid1_add_map1(sga_remap(ib_thread)%start_pos: &
615	sga_remap(ib_thread)%start_pos+ &
616	sga_remap(ib_thread)%num_links-1) = &
617	sga_remap(ib_thread)%grid1_add
618	grid2_add_map1(sga_remap(ib_thread)%start_pos: &
619	sga_remap(ib_thread)%start_pos+ &
620	sga_remap(ib_thread)%num_links-1) = &
621	sga_remap(ib_thread)%grid2_add
622	wts_map1 (:,sga_remap(ib_thread)%start_pos: &
623	sga_remap(ib_thread)%start_pos+ &
624	sga_remap(ib_thread)%num_links-1) = &
625	sga_remap(ib_thread)%wts
626
627	end do
628
629	#ifdef REMAP_TIMING
630	if (ll_timing) call timer_stop(3)
631	#endif
632	if (nlogprt.ge.2) then
633
634	do ib_thread = 1, il_nbthreads
635	if (sga_remap(ib_thread)%nb_resize.gt.0) then
636	write(nulou,*) ' Number of thread_resize_remap_vars on thread ',&
637	ib_thread, ' = ', sga_remap(ib_thread)%nb_resize
638	end if
639	end do
640
641	end if
642
643	#ifdef REMAP_TIMING
644	if (ll_timing.and.nlogprt.ge.2) then
645	write(nulou,*) ' On master thread '
646	call timer_print(2)
647	call timer_clear(2)
648	do ib_thread = 1,il_nbthreads
649	write(nulou,*) ' On thread ',ib_thread
650	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
651	& 'remap_dist_gaus_wgt distr ',dla_timer(ib_thread,1)
652	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
653	& 'remap_dist_gaus_wgt search ',dla_timer(ib_thread,2)
654	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
655	& 'remap_dist_gaus_wgt weights ',dla_timer(ib_thread,3)
656	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
657	& 'remap_dist_gaus_wgt store_link',dla_timer(ib_thread,4)
658	end do
659	deallocate(dla_timer)
660	write(nulou,*) ' On master thread '
661	call timer_print(3)
662	call timer_clear(3)
663	end if
664
665	else
666
667	if (ll_timing.and.nlogprt.ge.2) then
668	do i_n = 2, 6
669	call timer_print(i_n)
670	call timer_clear(i_n)
671	end do
672	end if
673	#endif
674	end if
675
676	! Gather the complete results on master proc
677
678	if (mpi_size_map .gt. 1) then
679
680	IF (mpi_rank_map == mpi_root_map) THEN
681	ALLOCATE(ila_num_links_mpi(mpi_size_map))
682	ELSE
683	ALLOCATE(ila_num_links_mpi(1))
684	END IF
685
686	CALL MPI_Gather (num_links_map1, 1,MPI_INT,&
687	& ila_num_links_mpi,1,MPI_INT,&
688	& mpi_root_map,mpi_comm_map,il_err)
689
690
691	IF (mpi_rank_map == mpi_root_map) THEN
692	num_links_map1 = SUM(ila_num_links_mpi)
693	if (num_links_map1 > max_links_map1) &
694	call resize_remap_vars(1,num_links_map1-max_links_map1)
695
696	ALLOCATE(ila_req_mpi(4,mpi_size_map-1))
697	ALLOCATE(ila_sta_mpi(MPI_STATUS_SIZE,4,mpi_size_map-1))
698
699	DO i_n = 1, mpi_size_map-1
700	buff_base = SUM(ila_num_links_mpi(1:i_n))+1
701	CALL MPI_IRecv(grid1_add_map1(buff_base),&
702	& ila_num_links_mpi(i_n+1),MPI_INT,i_n,1,mpi_comm_map,&
703	& ila_req_mpi(1,i_n),il_err)
704
705	CALL MPI_IRecv(grid2_add_map1(buff_base),&
706	& ila_num_links_mpi(i_n+1),MPI_INT,i_n,2,mpi_comm_map,&
707	& ila_req_mpi(2,i_n),il_err)
708
709	CALL MPI_IRecv(wts_map1(1,buff_base),&
710	& num_wts*ila_num_links_mpi(i_n+1),MPI_DOUBLE,i_n,0,mpi_comm_map,&
711	& ila_req_mpi(3,i_n),il_err)
712
713	CALL MPI_IRecv(grid2_frac(ila_mpi_mn(i_n+1)),&
714	& ila_mpi_mx(i_n+1)-ila_mpi_mn(i_n+1)+1,MPI_DOUBLE,i_n,0,mpi_comm_map,&
715	& ila_req_mpi(4,i_n),il_err)
716
717	END DO
718
719	DO i_n=1,4
720	CALL MPI_Waitall(mpi_size_map-1,ila_req_mpi(i_n,:),ila_sta_mpi(1,i_n,1),il_err)
721	END DO
722
723	DEALLOCATE(ila_req_mpi)
724	DEALLOCATE(ila_sta_mpi)
725
726	ELSE
727
728	CALL MPI_Send(grid1_add_map1,num_links_map1,MPI_INT,&
729	& mpi_root_map,1,mpi_comm_map,il_err)
730
731	CALL MPI_Send(grid2_add_map1,num_links_map1,MPI_INT,&
732	& mpi_root_map,2,mpi_comm_map,il_err)
733
734	CALL MPI_Send(wts_map1,num_wts*num_links_map1,MPI_DOUBLE,&
735	& mpi_root_map,0,mpi_comm_map,il_err)
736
737	CALL MPI_Send(grid2_frac(ila_mpi_mn(mpi_rank_map+1)),&
738	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
739	& mpi_root_map,0,mpi_comm_map,il_err)
740
741
742	END IF
743
744
745
746
747	! Reduction (SUM) of num_links_map1 on master proc
748	! Most probably it will be better to gather the single num_links_map1
749	! Because they are useful for storing the gatherv or Irecv results
750	!
751	! if (ll_master_remap) then
752	! allocate(ila_num_links_mpi(mpi_size_map))
753	! !! MPI receive the num_links_map1 of every proc in right pos
754	! num_links_map1 = SUM(ila_num_links_mpi)
755	! if (num_links_map1 > max_links_map1) &
756	! call resize_remap_vars(1,num_links_map1-max_links_map1)
757	! else
758	! !! MPI send of num_links_map1
759	! end IF
760	!
761	!
762	! Gather (ou similar) with positions derived from partial sums
763	! of ila_num_links mpi :
764	!grid1_add_map1
765	!grid2_add_map1
766	!wts_map1
767	!
768	! Gather (ou reduction sum) avec attention au decalage d'index de
769	!grid2_frac
770	!
771	! if (ll_master_remap) deallocate(ila_num_links_mpi)
772
773	deallocate(ila_num_links_mpi)
774
775	end if
776
777
778	deallocate (coslat, coslon, sinlat, sinlon)
779	deallocate(wgtstmp)
780	deallocate(ila_mpi_mn)
781	deallocate(ila_mpi_mx)
782	deallocate(ila_thr_mn)
783	deallocate(ila_thr_mx)
784	if (il_nbthreads .gt. 1) then
785	do ib_thread = 1, il_nbthreads
786	deallocate(sga_remap(ib_thread)%grid1_add)
787	deallocate(sga_remap(ib_thread)%grid2_add)
788	deallocate(sga_remap(ib_thread)%wts)
789	end do
790	deallocate(sga_remap)
791	end if
792	!
793	if (nlogprt .ge. 2) then
794	write (UNIT = nulou,FMT = *)' '
795	write (UNIT = nulou,FMT = *)'Leaving routine remap_dist_gaus_wgt'
796	call OASIS_FLUSH_SCRIP(nulou)
797	endif
798	!
799	!-----------------------------------------------------------------------
800
801	end subroutine remap_dist_gaus_wgt
802
803	!***********************************************************************
804
805	subroutine grid_search_nbr(src_addnn, nbr_add, nbr_dist, plat, plon, &
806	coslat_dst, coslon_dst, sinlat_dst, sinlon_dst, &
807	src_bin_add, num_neighbors, lextrapdone, dst_add )
808
809	!-----------------------------------------------------------------------
810	!
811	! this routine finds the closest num_neighbor points to a search
812	! point and computes a distance to each of the neighbors.
813	! it excludes redundant points as specified by the tolerance
814	! defined by dist_chk.
815	!
816	!-----------------------------------------------------------------------
817
818	!-----------------------------------------------------------------------
819	!
820	! input variables
821	!
822	!-----------------------------------------------------------------------
823
824	integer (kind=int_kind), dimension(:,:), intent(in) :: src_bin_add ! search bins for restricting search
825
826	integer (kind=int_kind), intent(in) :: dst_add
827
828	real (kind=dbl_kind), intent(in) :: plat, & ! latitude of the search point
829	plon, & ! longitude of the search point
830	coslat_dst, & ! cos(lat) of the search point
831	coslon_dst, & ! cos(lon) of the search point
832	sinlat_dst, & ! sin(lat) of the search point
833	sinlon_dst ! sin(lon) of the search point
834
835	logical, intent(in) :: lextrapdone ! logical, true if EXTRAP done on field
836
837	integer (kind=int_kind) :: num_neighbors ! number of neighbours
838
839	logical :: ll_debug
840
841
842	!-----------------------------------------------------------------------
843	!
844	! output variables
845	!
846	!-----------------------------------------------------------------------
847
848	integer (kind=int_kind), dimension(num_neighbors), intent(out) :: nbr_add ! address of each of the closest points
849
850	real (kind=dbl_kind), dimension(num_neighbors), intent(out) :: nbr_dist ! distance to each of the closest points
851
852	integer (kind=int_kind), intent(out) :: src_addnn
853
854	!-----------------------------------------------------------------------
855	!
856	! local variables
857	!
858	!-----------------------------------------------------------------------
859
860	integer (kind=int_kind) :: n, nmax, nadd, nchk, & ! dummy indices
861	nm1, np1, i, j, ip1, im1, jp1, jm1, &
862	min_add, max_add, &
863	il_debug_add, nbr_count
864
865	real (kind=dbl_kind) :: distance, rl_dist, src_latsnn, & ! angular distance
866	dist1, dist2 ! temporary dist calcs
867
868	real (kind=dbl_kind), dimension(num_neighbors) :: &
869	nbr_coslat, nbr_sinlat, nbr_coslon, nbr_sinlon ! stored nbr lon/lat
870
871	real (kind=dbl_kind), parameter :: dist_chk = 1.0e-7 ! delta distance limit radians
872
873	logical (kind=log_kind) :: nchkflag ! flag for nchk
874
875	!-----------------------------------------------------------------------
876	!
877	! loop over source grid and find nearest neighbors
878	!
879	!-----------------------------------------------------------------------
880
881	!***
882	!*** restrict the search using search bins
883	!*** expand the bins to catch neighbors
884	!***
885
886	select case (restrict_type)
887	case('LATITUDE')
888
889	do n=1,num_srch_bins
890	if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n)) then
891	min_add = src_bin_add(1,n)
892	max_add = src_bin_add(2,n)
893
894	nm1 = max(n-1,1)
895	np1 = min(n+1,num_srch_bins)
896
897	min_add = min(min_add,src_bin_add(1,nm1))
898	max_add = max(max_add,src_bin_add(2,nm1))
899	min_add = min(min_add,src_bin_add(1,np1))
900	max_add = max(max_add,src_bin_add(2,np1))
901	endif
902	end do
903
904	case('LATLON')
905
906	n = 0
907	nmax = nint(sqrt(real(num_srch_bins)))
908	do j=1,nmax
909	jp1 = min(j+1,nmax)
910	jm1 = max(j-1,1)
911	do i=1,nmax
912	ip1 = min(i+1,nmax)
913	im1 = max(i-1,1)
914
915	n = n+1
916	if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and. &
917	plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then
918	min_add = src_bin_add(1,n)
919	max_add = src_bin_add(2,n)
920
921	nm1 = (jm1-1)*nmax + im1
922	np1 = (jp1-1)*nmax + ip1
923	nm1 = max(nm1,1)
924	np1 = min(np1,num_srch_bins)
925
926	min_add = min(min_add,src_bin_add(1,nm1))
927	max_add = max(max_add,src_bin_add(2,nm1))
928	min_add = min(min_add,src_bin_add(1,np1))
929	max_add = max(max_add,src_bin_add(2,np1))
930	endif
931	end do
932	end do
933
934	end select
935
936	!***
937	!*** initialize distance and address arrays
938	!***
939
940	nbr_count = 0
941	nbr_add = 0
942	nbr_dist = bignum
943	nbr_coslat = bignum
944	nbr_sinlat = bignum
945	nbr_coslon = bignum
946	nbr_sinlon = bignum
947	src_latsnn = bignum
948	src_addnn = 0
949
950	do nadd=min_add,max_add
951
952	!***
953	!*** find distance to this point
954	!***
955	rl_dist = sinlat_dst*sinlat(nadd) + &
956	coslat_dstcoslat(nadd) &
957	(coslon_dst*coslon(nadd) + &
958	sinlon_dst*sinlon(nadd))
959	if (rl_dist < -1.0d0) then
960	rl_dist = -1.0d0
961	else if (rl_dist > 1.0d0) then
962	rl_dist = 1.0d0
963	end if
964	distance = acos(rl_dist)
965	ll_debug = .false.
966	if (ll_debug) then
967	il_debug_add = 15248
968	if (dst_add == il_debug_add) then
969	write(nulou,1009) nadd, distance
970	endif
971	endif
972
973	!*** store the address and distance if this is one of the
974	!*** smallest four so far
975	!***
976
977	if (distance .lt. nbr_dist(num_neighbors)) then
978	! compute nchk, the first neighbor with ge distance than nadd
979	nchk = num_neighbors
980	nchkflag = .true.
981	do while (nchkflag)
982	if (distance .lt. nbr_dist(nchk-1)) then
983	nchk = nchk - 1
984	else
985	nchkflag = .false.
986	endif
987	if (nchk == 1) nchkflag = .false.
988	enddo
989
990	! check that points are not the same, need to compare nchk-1 and nchk
991	! only compare against neighbors that have been initialized (nchk <= nbr_count)
992	! nchk-1 should always be an initialized point if nchk > 1
993	! nchk may be initialized or not
994	if (nchk > 1) then
995	dist1 = nbr_sinlat(nchk-1)*sinlat(nadd) + &
996	nbr_coslat(nchk-1)coslat(nadd) &
997	(nbr_coslon(nchk-1)*coslon(nadd) + &
998	nbr_sinlon(nchk-1)*sinlon(nadd))
999	dist1 = MAX ( MIN ( dist1, 1.d0), -1.d0)
1000	!dist1 = acos(dist1)
1001	! avoid cost of acos above since we're comparing to a small number
1002	! use cos(x) = 1-x^2/2, x = sqrt(2 * (1 - cos(x)))
1003	dist1 = sqrt(2.0d0*(1.0d0-abs(dist1)))
1004	else
1005	dist1 = bignum
1006	endif
1007
1008	if (dist1 > dist_chk) then
1009	if (nchk <= nbr_count) then
1010	dist2 = nbr_sinlat(nchk)*sinlat(nadd) + &
1011	nbr_coslat(nchk)coslat(nadd) &
1012	(nbr_coslon(nchk)*coslon(nadd) + &
1013	nbr_sinlon(nchk)*sinlon(nadd))
1014	dist2 = MAX ( MIN ( dist2, 1.d0), -1.d0)
1015	!dist2 = acos(dist2)
1016	! avoid cost of acos above since we're comparing to a small number
1017	! use cos(x) = 1-x^2/2, x = sqrt(2 * (1 - cos(x)))
1018	dist2 = sqrt(2.0d0*(1.0d0-abs(dist2)))
1019	else
1020	dist2 = bignum
1021	endif
1022
1023	if (dist2 > dist_chk) then
1024	nbr_count = min(nbr_count + 1, num_neighbors)
1025	do n=num_neighbors,nchk+1,-1
1026	nbr_add(n) = nbr_add(n-1)
1027	nbr_dist(n) = nbr_dist(n-1)
1028	nbr_coslat(n) = nbr_coslat(n-1)
1029	nbr_sinlat(n) = nbr_sinlat(n-1)
1030	nbr_coslon(n) = nbr_coslon(n-1)
1031	nbr_sinlon(n) = nbr_sinlon(n-1)
1032	end do
1033	nbr_add(nchk) = nadd
1034	nbr_dist(nchk) = distance
1035	nbr_coslat(nchk) = coslat(nadd)
1036	nbr_sinlat(nchk) = sinlat(nadd)
1037	nbr_coslon(nchk) = coslon(nadd)
1038	nbr_sinlon(nchk) = sinlon(nadd)
1039	else ! dist2
1040	! if (nlogprt .ge. 2) then
1041	! write(nulou,*) 'remap_dist_gaus_wgt nbr_search skip point2: ',dist2
1042	! write(nulou,*) ' ',nadd,nbr_add(nchk)
1043	! write(nulou,*) ' ',nbr_coslat(nchk),coslat(nadd)
1044	! write(nulou,*) ' ',nbr_sinlat(nchk),sinlat(nadd)
1045	! write(nulou,*) ' ',nbr_coslon(nchk),coslon(nadd)
1046	! write(nulou,*) ' ',nbr_sinlon(nchk),sinlon(nadd)
1047	! endif
1048	endif ! dist2
1049	else ! dist1
1050	! if (nlogprt .ge. 2) then
1051	! write(nulou,*) 'remap_dist_gaus_wgt nbr_search skip point1: ',dist1
1052	! write(nulou,*) ' ',nadd,nbr_add(nchk-1)
1053	! write(nulou,*) ' ',nbr_coslat(nchk-1),coslat(nadd)
1054	! write(nulou,*) ' ',nbr_sinlat(nchk-1),sinlat(nadd)
1055	! write(nulou,*) ' ',nbr_coslon(nchk-1),coslon(nadd)
1056	! write(nulou,*) ' ',nbr_sinlon(nchk-1),sinlon(nadd)
1057	! endif
1058	endif ! dist1
1059	endif ! distance
1060
1061	if (ll_debug) then
1062	if (dst_add == il_debug_add) then
1063	write(nulou,1010) nadd, distance
1064	endif
1065	endif
1066
1067	if (ll_nnei_dgw) then
1068	!***
1069	!*** store the non-masked closest neighbour
1070	!***
1071	if (distance < src_latsnn) then
1072	if (grid1_mask(nadd) .or. lextrapdone) then
1073	src_addnn = nadd
1074	src_latsnn = distance
1075	endif
1076	if (ll_debug) then
1077	if (dst_add == il_debug_add) then
1078	write(nulou,1010) nadd, distance
1079	endif
1080	endif
1081	endif
1082	endif
1083	end do
1084
1085	1009 format (1X, 'nadd0 =', 1X, I6, 2X, 'distance0 =', 1X, F18.16)
1086	1010 format (1X, 'nadd1 =', 1X, I6, 2X, 'distance0 =', 1X, F18.16)
1087	1011 format (1X, 'nadd2 =', 1X, I6, 2X, 'distance2 =', 1X, F18.16)
1088
1089	!-----------------------------------------------------------------------
1090
1091	end subroutine grid_search_nbr
1092
1093	!***********************************************************************
1094
1095	subroutine store_link_nbr(add1, add2, weights, nmap, id_thread)
1096
1097	!-----------------------------------------------------------------------
1098	!
1099	! this routine stores the address and weight for this link in
1100	! the appropriate address and weight arrays and resizes those
1101	! arrays if necessary.
1102	!
1103	!-----------------------------------------------------------------------
1104
1105	!-----------------------------------------------------------------------
1106	!
1107	! input variables
1108	!
1109	!-----------------------------------------------------------------------
1110
1111	integer (kind=int_kind), intent(in) :: add1, & ! address on grid1
1112	add2, & ! address on grid2
1113	nmap, & ! identifies which direction for mapping
1114	id_thread ! local threaded task
1115
1116	real (kind=dbl_kind), dimension(:), intent(in) :: weights ! array of remapping weights for this link
1117
1118	!-----------------------------------------------------------------------
1119	!
1120	! increment number of links and check to see if remap arrays need
1121	! to be increased to accomodate the new link. then store the
1122	! link.
1123	!
1124	!-----------------------------------------------------------------------
1125
1126	if (il_nbthreads .eq. 1) then
1127
1128
1129	select case (nmap)
1130	case(1)
1131
1132	num_links_map1 = num_links_map1 + 1
1133
1134	if (num_links_map1 > max_links_map1) &
1135	call resize_remap_vars(1,resize_increment)
1136
1137	grid1_add_map1(num_links_map1) = add1
1138	grid2_add_map1(num_links_map1) = add2
1139	wts_map1 (:,num_links_map1) = weights
1140
1141	case(2)
1142
1143	num_links_map2 = num_links_map2 + 1
1144
1145	if (num_links_map2 > max_links_map2) &
1146	call resize_remap_vars(2,resize_increment)
1147
1148	grid1_add_map2(num_links_map2) = add1
1149	grid2_add_map2(num_links_map2) = add2
1150	wts_map2 (:,num_links_map2) = weights
1151
1152	end select
1153
1154	else
1155
1156	sga_remap(id_thread)%num_links = sga_remap(id_thread)%num_links + 1
1157
1158	if (sga_remap(id_thread)%num_links > sga_remap(id_thread)%max_links) &
1159	call sga_remap(id_thread)%resize(int(0.2*real(sga_remap(id_thread)%max_links)))
1160
1161	sga_remap(id_thread)%grid1_add(sga_remap(id_thread)%num_links) = add1
1162	sga_remap(id_thread)%grid2_add(sga_remap(id_thread)%num_links) = add2
1163	sga_remap(id_thread)%wts(:,sga_remap(id_thread)%num_links) = weights
1164
1165	endif
1166
1167	!-----------------------------------------------------------------------
1168
1169	end subroutine store_link_nbr
1170
1171	!***********************************************************************
1172
1173	end module remap_distance_gaussian_weight
1174
1175	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1176
1177
1178

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

Download in other formats: