Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

remap_bicubic_reduced.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: 55.2 KB

Line
1	MODULE remap_bicubic_reduced
2
3	!-----------------------------------------------------------------------
4	! BOP
5	!
6	! !MODULE: remap_bicubic_reduced
7	!
8	! !USES:
9	USE kinds_mod ! defines common data types
10	USE constants ! defines common constants
11	USE grids ! module containing grid info
12	USE timers
13	USE remap_vars ! module containing remap info
14	USE mod_oasis_flush
15	!$ use omp_lib
16
17	! !PUBLIC TYPES:
18	IMPLICIT NONE
19	! !PUBLIC MEMBER FUNCTIONS:
20	!
21	! !PUBLIC DATA MEMBERS:
22
23	real (kind=dbl_kind), dimension(:), allocatable, save :: &
24	coslat, sinlat, & ! cosine, sine of grid lats (for distance)
25	coslon, sinlon ! cosine, sine of grid lons (for distance)
26
27	integer (kind=int_kind) :: il_nbthreads = 1
28
29	integer (kind=int_kind), PARAMETER :: num_neighbors = 16
30
31	! !DESCRIPTION:
32	! This routine computes the weights for a bicubic interpolation
33	! with a reduced grid. Computes mappings from grid1 to grid2.
34	!
35	! !REVISION HISTORY:
36	! 2002.10.21 J.Ghattas created
37	!
38	! EOP
39	!-----------------------------------------------------------------------
40	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
41	! $Author: valcke $
42	!-----------------------------------------------------------------------
43
44
45
46	CONTAINS
47
48	!***********************************************************************
49
50
51	!-----------------------------------------------------------------------
52	! BOP
53	!
54	! !IROUTINE: remap_bicub_reduced
55	!
56	! !INTERFACE:
57
58	SUBROUTINE remap_bicub_reduced(ld_extrapdone, ll_nnei, &
59	mpi_comm_map, mpi_size_map, mpi_rank_map, mpi_root_map)
60
61
62	! !USES:
63
64	! !RETURN VALUE:
65
66	! !PARAMETERS:
67
68	LOGICAL, INTENT(in) :: &
69	ld_extrapdone ! logical, true if EXTRAP done on field
70
71	integer (kind=int_kind), INTENT(in) :: &
72	mpi_comm_map, mpi_rank_map, mpi_size_map, mpi_root_map
73
74	LOGICAL, INTENT(in) :: ll_nnei ! true (default) if extra search is done
75
76	INTEGER (KIND=int_kind), DIMENSION(4,4) :: &
77	ila_src_add ! address for source points non-masked
78
79	INTEGER (KIND=int_kind), DIMENSION(4) :: &
80	ila_nbr_found ! nrb of points found on each latitude
81
82	INTEGER (KIND=int_kind) :: &
83	ib_i, & ! iter index
84	ib_dst_add, & ! destination address, target point
85	il_count, & ! nbr of latitudes with found points
86	il_min, il_max, bin ! begin and end for distances calculation
87
88	REAL (KIND=dbl_kind), DIMENSION(4,4) :: &
89	rla_src_lons, & ! longitudes for the points 'ila_src_add'
90	rla_weight, & ! bicubic weights for the points 'ila_src_add'
91	rla_wght_lon ! temp. weights
92
93	REAL (KIND=dbl_kind), DIMENSION(4) :: &
94	rla_src_lats, & ! latitudes for the points 'ila_src_add'
95	rla_lats_temp, & ! temp. latitudes
96	rla_wght_lat, rla_wght_temp! temp. weights
97
98	REAL (KIND=dbl_kind) :: &
99	rl_plat, rl_plon ! latitude and longitude for destination address
100
101	REAL (KIND=dbl_kind) :: & ! variables for distances calculation
102	rl_coslat_dst, rl_sinlat_dst, &
103	rl_coslon_dst, rl_sinlon_dst, &
104	rl_distance, d_dist
105
106	REAL (KIND=dbl_kind), DIMENSION(2) :: &
107	rla_dist ! lat distances to point cible
108
109	INTEGER (KIND=int_kind), DIMENSION(4) :: &
110	ila_add_dist ! temporary addr. for distances
111
112	LOGICAL :: ll_linear ! flag
113
114	INTEGER (KIND=int_kind) :: il_splitsize
115	INTEGER (KIND=int_kind) :: ib_proc
116	INTEGER (KIND=int_kind) :: ib_thread
117	INTEGER (KIND=int_kind) :: buff_base
118	INTEGER (KIND=int_kind), DIMENSION(:), ALLOCATABLE :: ila_mpi_sz
119	INTEGER (KIND=int_kind), DIMENSION(:), ALLOCATABLE :: ila_mpi_mn
120	INTEGER (KIND=int_kind), DIMENSION(:), ALLOCATABLE :: ila_mpi_mx
121	INTEGER (KIND=int_kind), DIMENSION(:), ALLOCATABLE :: ila_thr_sz
122	INTEGER (KIND=int_kind), DIMENSION(:), ALLOCATABLE :: ila_thr_mn
123	INTEGER (KIND=int_kind), DIMENSION(:), ALLOCATABLE :: ila_thr_mx
124	INTEGER (KIND=int_kind), DIMENSION(:), ALLOCATABLE :: ila_num_links_mpi
125	INTEGER (KIND=int_kind), DIMENSION(:,:), ALLOCATABLE :: ila_req_mpi
126	INTEGER (KIND=int_kind), DIMENSION(:,:,:), ALLOCATABLE :: ila_sta_mpi
127	CHARACTER (LEN=14) :: cl_envvar
128	INTEGER (KIND=int_kind) :: il_envthreads, il_err
129
130
131	#ifdef REMAP_TIMING
132	LOGICAL :: ll_timing = .true.
133	REAL (KIND=dbl_kind), DIMENSION(:,:), ALLOCATABLE :: dla_timer
134	REAL (KIND=dbl_kind) :: dl_tstart, dl_tstop
135	INTEGER (KIND=int_kind) :: il_mythread, n
136	#endif
137
138
139	! !DESCRIPTION:
140	! This routine computes the weights for a bicubic interpolation
141	! with a reduced grid. Computes mappings from grid1 to grid2.
142	!
143	! !REVISION HISTORY:
144	! 2002.10.21 J.Ghattas created
145	!
146	! EOP
147	!-----------------------------------------------------------------------
148	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
149	! $Author: valcke $
150	!-----------------------------------------------------------------------
151	!
152	IF (nlogprt .GE. 2) THEN
153	WRITE (UNIT = nulou,FMT = *)' '
154	WRITE (UNIT = nulou,FMT = *) 'Entering routine remap_bicub_reduced'
155	CALL OASIS_FLUSH_SCRIP(nulou)
156	ENDIF
157
158	#ifdef REMAP_TIMING
159	if (ll_timing) call timer_start(2,'remap_bicubic_reduced alloc')
160	#endif
161
162	!***
163	!*** compute cos, sin of lat/lon on source grid for distance
164	!*** calculations
165	!***
166
167	allocate (coslat(grid1_size), coslon(grid1_size), sinlat(grid1_size), sinlon(grid1_size))
168
169	coslat = cos(grid1_center_lat)
170	coslon = cos(grid1_center_lon)
171	sinlat = sin(grid1_center_lat)
172	sinlon = sin(grid1_center_lon)
173
174	#ifdef REMAP_TIMING
175	if (ll_timing) call timer_stop(2)
176	#endif
177
178
179	!***
180	!*** precompute best scheduling of target grid points
181	!***
182
183	allocate(ila_mpi_mn(mpi_size_map))
184	allocate(ila_mpi_mx(mpi_size_map))
185
186	if (mpi_size_map .gt. 1) then
187
188	allocate(ila_mpi_sz(mpi_size_map))
189	il_splitsize = count(grid2_mask)
190	ila_mpi_sz(:) = floor(real(il_splitsize)/mpi_size_map)
191	ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) = ila_mpi_sz(1:il_splitsize-sum(ila_mpi_sz)) + 1
192
193	ila_mpi_mn(:) = 0
194	ib_proc = 1
195	il_splitsize = 0
196	do ib_dst_add = 1, grid2_size
197	if (grid2_mask(ib_dst_add)) then
198	if (ila_mpi_mn(ib_proc).eq.0) &
199	ila_mpi_mn(ib_proc) = ib_dst_add
200	il_splitsize = il_splitsize + 1
201	if (il_splitsize .eq. ila_mpi_sz(ib_proc)) then
202	il_splitsize = 0
203	ila_mpi_mx(ib_proc) = ib_dst_add
204	ib_proc = ib_proc + 1
205	end if
206	end if
207	end do
208	!
209	deallocate(ila_mpi_sz)
210
211	else
212
213	ila_mpi_mn(1) = 1
214	ila_mpi_mx(1) = grid2_size
215
216	endif
217
218
219	call get_environment_variable(name='OASIS_OMP_NUM_THREADS', value=cl_envvar, status=il_err)
220	if ( il_err .ne. 0) then
221	il_envthreads = 0
222	else
223	read(cl_envvar,*) il_envthreads
224	end if
225
226	!$OMP PARALLEL NUM_THREADS(il_envthreads) DEFAULT(NONE) &
227	!$OMP SHARED(ld_extrapdone) &
228	!$OMP SHARED(grid2_mask) &
229	!$OMP SHARED(il_envthreads) &
230	!$OMP SHARED(grid2_center_lat,grid2_center_lon) &
231	!$OMP SHARED(grid1_center_lat,grid1_center_lon) &
232	!$OMP SHARED(grid1_mask) &
233	!$OMP SHARED(bin_addr1_r,num_srch_red,sga_remap,num_wts) &
234	!$OMP SHARED(nulou) &
235	!$OMP SHARED(il_nbthreads) &
236	!$OMP SHARED(mpi_rank_map,mpi_root_map,ila_mpi_mn,ila_mpi_mx) &
237	!$OMP SHARED(ila_thr_sz,ila_thr_mn,ila_thr_mx) &
238	!$OMP SHARED(grid1_dims,grid1_size) &
239	!$OMP PRIVATE(nlogprt) &
240	!$OMP PRIVATE(ib_dst_add) &
241	!$OMP PRIVATE(ila_src_add) &
242	!$OMP PRIVATE(rl_coslat_dst,rl_coslon_dst,rl_sinlat_dst,rl_sinlon_dst) &
243	!$OMP PRIVATE(rl_plat,rl_plon) &
244	!$OMP PRIVATE(rla_src_lats,rla_src_lons) &
245	!$OMP PRIVATE(ila_nbr_found,ila_add_dist,bin,il_min,il_max,ib_i,ll_linear) &
246	!$OMP PRIVATE(rla_dist,rla_wght_lon,rla_wght_lat,rla_wght_temp,rla_lats_temp) &
247	!$OMP PRIVATE(d_dist) &
248	!$OMP PRIVATE(rla_weight,rl_distance,il_count) &
249	!$OMP SHARED(coslat,coslon,sinlat,sinlon) &
250	!$OMP SHARED(ll_nnei) &
251	#ifdef REMAP_TIMING
252	!$OMP PRIVATE(ib_thread,il_splitsize) &
253	!$OMP SHARED(ll_timing,dla_timer) &
254	!$OMP PRIVATE(il_mythread,dl_tstart,dl_tstop)
255
256	!$ il_mythread = OMP_GET_THREAD_NUM () + 1
257	#else
258	!$OMP PRIVATE(ib_thread,il_splitsize)
259	#endif
260
261	!$OMP SINGLE
262
263	il_nbthreads = 1
264	!$ il_nbthreads = OMP_GET_NUM_THREADS ()
265
266
267	#ifdef REMAP_TIMING
268	if (ll_timing) then
269	if (il_nbthreads.gt.1) then
270	!$ dl_tstart = OMP_GET_WTIME()
271	else
272	call timer_start(3,'remap_bicubic_reduced distr')
273	end if
274	end if
275	#endif
276
277	allocate(ila_thr_mn(il_nbthreads))
278	allocate(ila_thr_mx(il_nbthreads))
279
280	if (il_nbthreads .gt. 1) then
281
282	#ifdef REMAP_TIMING
283	if (ll_timing) then
284	allocate(dla_timer(il_nbthreads,6))
285	dla_timer(:,:) = 0.0
286	end if
287	#endif
288	nlogprt = 0
289
290	allocate(ila_thr_sz(il_nbthreads))
291	il_splitsize = COUNT(grid2_mask(ila_mpi_mn(mpi_rank_map+1):&
292	& ila_mpi_mx(mpi_rank_map+1)))
293	ila_thr_sz(:) = floor(real(il_splitsize)/il_nbthreads)
294	ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) = ila_thr_sz(1:il_splitsize-sum(ila_thr_sz)) + 1
295
296	ila_thr_mn(:) = 0
297	ib_thread = 1
298	il_splitsize = 0
299	do ib_dst_add = ila_mpi_mn(mpi_rank_map+1), ila_mpi_mx(mpi_rank_map+1)
300	if (grid2_mask(ib_dst_add)) then
301	if (ila_thr_mn(ib_thread).eq.0) &
302	ila_thr_mn(ib_thread) = ib_dst_add
303	il_splitsize = il_splitsize + 1
304	if (il_splitsize .eq. ila_thr_sz(ib_thread)) then
305	il_splitsize = 0
306	ila_thr_mx(ib_thread) = ib_dst_add
307	ib_thread = ib_thread + 1
308	end if
309	end if
310	end do
311
312	allocate(sga_remap(il_nbthreads))
313
314	do ib_thread = 1, il_nbthreads
315	il_splitsize = num_neighbors*ila_thr_sz(ib_thread)
316	sga_remap(ib_thread)%max_links = il_splitsize
317	sga_remap(ib_thread)%num_links = 0
318	allocate(sga_remap(ib_thread)%grid1_add(il_splitsize))
319	allocate(sga_remap(ib_thread)%grid2_add(il_splitsize))
320	allocate(sga_remap(ib_thread)%wts(num_wts,il_splitsize))
321	end do
322
323	deallocate(ila_thr_sz)
324
325	else
326
327	ila_thr_mn(1) = ila_mpi_mn(mpi_rank_map+1)
328	ila_thr_mx(1) = ila_mpi_mx(mpi_rank_map+1)
329
330	end if
331	#ifdef REMAP_TIMING
332	if (ll_timing) then
333	if (il_nbthreads.gt.1) then
334	!$ dl_tstop = OMP_GET_WTIME()
335	dla_timer(il_mythread,1)=dla_timer(il_mythread,1) + dl_tstop - dl_tstart
336	else
337	call timer_stop(3)
338	end if
339	end if
340	#endif
341	!$OMP END SINGLE
342
343	!***
344	!*** loop over destination grid
345	!***
346	!$OMP DO SCHEDULE(STATIC,1)
347	thread_loop: do ib_thread = 1, il_nbthreads
348
349	grid_loop1:DO ib_dst_add = ila_thr_mn(ib_thread), ila_thr_mx(ib_thread)
350
351	ll_linear=.false.
352
353	IF (.NOT. grid2_mask(ib_dst_add)) THEN
354	CYCLE grid_loop1 ! target point is masked
355	END IF
356
357	#ifdef REMAP_TIMING
358	if (ll_timing) then
359	if (il_nbthreads.gt.1) then
360	!$ dl_tstart = OMP_GET_WTIME()
361	else
362	call timer_start(4,'remap_bicubic_reduced search')
363	end if
364	end if
365	#endif
366
367	rl_plat = grid2_center_lat(ib_dst_add)
368	rl_plon = grid2_center_lon(ib_dst_add)
369
370	!
371	! Search for non-masked points among the closest 16 points
372	! on source grid (grid1)
373	!
374	CALL grid_search_16_points(ila_src_add, rla_src_lats, rla_src_lons,&
375	ila_nbr_found, bin, rl_plat, &
376	rl_plon, ld_extrapdone)
377
378	#ifdef REMAP_TIMING
379	if (ll_timing) then
380	if (il_nbthreads.gt.1) then
381	!$ dl_tstop = OMP_GET_WTIME()
382	dla_timer(ib_thread,2) = dla_timer(ib_thread,2) + dl_tstop - dl_tstart
383	else
384	call timer_stop(4)
385	end if
386	end if
387	#endif
388	!
389	! If there is no point found, search the nearest
390	! non-masked point
391	!
392	#ifdef REMAP_TIMING
393	if (ll_timing) then
394	if (il_nbthreads.gt.1) then
395	!$ dl_tstart = OMP_GET_WTIME()
396	else
397	call timer_start(5,'remap_bicubic_reduced partial neighbours')
398	end if
399	end if
400	#endif
401
402	IF (SUM(ila_nbr_found)==0) THEN
403	IF (ll_nnei .EQV. .TRUE. ) THEN
404	IF (nlogprt .GE. 2) THEN
405	WRITE(nulou,*) ' '
406	WRITE(nulou,*) &
407	'All 16 surrounding source grid points are masked'
408	WRITE(nulou,*) 'for target point ',ib_dst_add
409	WRITE(nulou,*) 'with longitude and latitude', rl_plon, rl_plat
410	WRITE(nulou,*) 'Using the nearest non-masked neighbour.'
411	WRITE(nulou,*) ' '
412	CALL OASIS_FLUSH_SCRIP(nulou)
413	ENDIF
414
415	! Search the nearest point in bin [il_min:il_max]
416	IF (bin==0 .or. bin==1) THEN
417	il_min=1
418	il_max=bin_addr1_r(2,3)
419	ELSE IF (bin==num_srch_red .or. bin==num_srch_red-1) THEN
420	il_min=bin_addr1_r(1,num_srch_red-2)
421	il_max=bin_addr1_r(2,num_srch_red)
422	ELSE
423	il_min=bin_addr1_r(1,bin-1)+1
424	il_max=bin_addr1_r(2,bin+2)
425	END IF
426
427	rl_coslat_dst = COS(rl_plat)
428	rl_sinlat_dst = SIN(rl_plat)
429	rl_coslon_dst = COS(rl_plon)
430	rl_sinlon_dst = SIN(rl_plon)
431
432	rla_weight(1,1) = bignum
433	ila_src_add(1,1) = 0
434	!cdir novector
435	DO ib_i=il_min, il_max
436	IF (grid1_mask(ib_i) .or. ld_extrapdone) THEN
437	d_dist = rl_coslat_dstcoslat(ib_i) &
438	(rl_coslon_dst*coslon(ib_i) + &
439	rl_sinlon_dst*sinlon(ib_i))+&
440	rl_sinlat_dst*sinlat(ib_i)
441	IF (d_dist < -1.0d0) THEN
442	d_dist = -1.0d0
443	ELSE IF (d_dist > 1.0d0) THEN
444	d_dist = 1.0d0
445	END IF
446	rl_distance = ACOS(d_dist)
447	IF (rl_distance < rla_weight(1,1)) THEN
448	rla_weight(1,1) = rl_distance
449	ila_src_add(1,1) = ib_i
450	END IF
451	END IF
452	END DO
453
454	IF (ila_src_add(1,1) == 0) THEN
455	rla_weight(1,1) = bignum
456	!cdir novector
457	DO ib_i=1,grid1_size
458	IF (grid1_mask(ib_i) .or. ld_extrapdone) THEN
459	d_dist = rl_coslat_dstcoslat(ib_i) &
460	(rl_coslon_dst*coslon(ib_i) + &
461	rl_sinlon_dst*sinlon(ib_i))+&
462	rl_sinlat_dst*sinlat(ib_i)
463	IF (d_dist < -1.0d0) THEN
464	d_dist = -1.0d0
465	ELSE IF (d_dist > 1.0d0) THEN
466	d_dist = 1.0d0
467	END IF
468	rl_distance = ACOS(d_dist)
469	IF (rl_distance < rla_weight(1,1)) THEN
470	rla_weight(1,1) = rl_distance
471	ila_src_add(1,1) = ib_i
472	END IF
473	END IF
474	END DO
475	ENDIF
476
477	if (ila_src_add(1,1) == 0) then
478	WRITE(nulou,*) 'Problem with neighbour identification for target grid point'
479	WRITE(nulou,*) 'with address = ',ib_dst_add
480	call abort
481	endif
482
483	rla_weight(:,:) = 0
484	rla_weight(1,1) = 1
485	CALL store_link_bicub(ib_dst_add, ila_src_add, rla_weight, ib_thread)
486	IF (nlogprt .GE. 2) THEN
487	WRITE(nulou,*) &
488	'Nearest non masked neighbour is source point ', &
489	ila_src_add(1,1)
490	WRITE(nulou,*) 'with longitude and latitude', &
491	grid1_center_lon(ila_src_add(1,1)), &
492	grid1_center_lat(ila_src_add(1,1))
493	WRITE(nulou,*) ' '
494	ENDIF
495	#ifdef REMAP_TIMING
496	if (ll_timing) then
497	if (il_nbthreads.gt.1) then
498	!$ dl_tstop = OMP_GET_WTIME()
499	dla_timer(ib_thread,3) = dla_timer(ib_thread,3) + dl_tstop - dl_tstart
500	else
501	call timer_stop(5)
502	end if
503	end if
504	#endif
505	CYCLE grid_loop1
506	ENDIF
507	END IF
508
509	rla_weight(:,:) = 0
510	! if there is only one point found, save it
511	IF (SUM(ila_nbr_found)==1) THEN
512	DO ib_i=1,4
513	IF (ila_nbr_found(ib_i)==1) THEN
514	rla_weight(ib_i,1)=1
515	EXIT
516	END IF
517	END DO
518	CALL store_link_bicub(ib_dst_add, ila_src_add, rla_weight, ib_thread)
519	#ifdef REMAP_TIMING
520	if (ll_timing) then
521	if (il_nbthreads.gt.1) then
522	!$ dl_tstop = OMP_GET_WTIME()
523	dla_timer(ib_thread,3) = dla_timer(ib_thread,3) + dl_tstop - dl_tstart
524	else
525	call timer_stop(5)
526	end if
527	end if
528	#endif
529	CYCLE grid_loop1
530	END IF
531
532	! if there are only 2 points found, distance weighted average
533	IF (SUM(ila_nbr_found)==2) THEN
534	rl_coslat_dst = COS(rl_plat)
535	rl_sinlat_dst = SIN(rl_plat)
536	rl_coslon_dst = COS(rl_plon)
537	rl_sinlon_dst = SIN(rl_plon)
538
539	rl_distance=0 ! count of total distance
540	DO ib_i=1,4
541	IF (ila_nbr_found(ib_i) > 0) THEN
542	d_dist = rl_coslat_dstCOS(rla_src_lats(ib_i)) &
543	(rl_coslon_dst*COS(rla_src_lons(ib_i,1)) + &
544	rl_sinlon_dst*SIN(rla_src_lons(ib_i,1)))+&
545	rl_sinlat_dst*SIN(rla_src_lats(ib_i))
546	IF (d_dist < -1.0d0) THEN
547	d_dist = -1.0d0
548	ELSE IF (d_dist > 1.0d0) THEN
549	d_dist = 1.0d0
550	END IF
551	rla_weight(ib_i,1) = ACOS(d_dist)
552	rl_distance = rl_distance+rla_weight(ib_i,1)
553	IF (ila_nbr_found(ib_i)==2) THEN
554	d_dist = rl_coslat_dstCOS(rla_src_lats(ib_i)) &
555	(rl_coslon_dst*COS(rla_src_lons(ib_i,2)) + &
556	rl_sinlon_dst*SIN(rla_src_lons(ib_i,2)))+&
557	rl_sinlat_dst*SIN(rla_src_lats(ib_i))
558	IF (d_dist < -1.0d0) THEN
559	d_dist = -1.0d0
560	ELSE IF (d_dist > 1.0d0) THEN
561	d_dist = 1.0d0
562	END IF
563	rla_weight(ib_i,2) = ACOS(d_dist)
564	rl_distance = rl_distance+rla_weight(ib_i,2)
565	END IF
566	END IF
567	END DO
568	rla_weight=rla_weight/rl_distance
569
570	CALL store_link_bicub(ib_dst_add, ila_src_add, rla_weight, ib_thread)
571	#ifdef REMAP_TIMING
572	if (ll_timing) then
573	if (il_nbthreads.gt.1) then
574	!$ dl_tstop = OMP_GET_WTIME()
575	dla_timer(ib_thread,3) = dla_timer(ib_thread,3) + dl_tstop - dl_tstart
576	else
577	call timer_stop(5)
578	end if
579	end if
580	#endif
581	CYCLE grid_loop1
582	END IF
583
584	#ifdef REMAP_TIMING
585	if (ll_timing) then
586	if (il_nbthreads.gt.1) then
587	!$ dl_tstop = OMP_GET_WTIME()
588	dla_timer(ib_thread,3) = dla_timer(ib_thread,3) + dl_tstop - dl_tstart
589	else
590	call timer_stop(5)
591	end if
592	end if
593	#endif
594
595	#ifdef REMAP_TIMING
596	if (ll_timing) then
597	if (il_nbthreads.gt.1) then
598	!$ dl_tstart = OMP_GET_WTIME()
599	else
600	call timer_start(6,'remap_bicubic_reduced weights')
601	end if
602	end if
603	#endif
604
605	! Some case exceptional when just one point per line found
606
607	IF (ila_nbr_found(1)==1) THEN ! elimination of point
608	ila_nbr_found(1)=0
609	ila_src_add(1,1)=0
610	END IF
611	IF (ila_nbr_found(4)==1) THEN
612	ila_nbr_found(4)=0
613	ila_src_add(4,1)=0
614	END IF
615
616	IF (ila_nbr_found(2)==1 .or. ila_nbr_found(3)==1) THEN
617	ila_add_dist(:)=4
618	rla_dist(:)=bignum
619
620	! search for the 2 nearest points or line of points
621	DO ib_i=1,4
622	IF (ila_nbr_found(ib_i) > 1) THEN
623	rl_distance=ABS(rla_src_lats(ib_i)-rl_plat)
624	ELSE IF (ila_nbr_found(ib_i)==1) THEN
625	rl_coslat_dst = COS(rl_plat)
626	rl_sinlat_dst = SIN(rl_plat)
627	rl_coslon_dst = COS(rl_plon)
628	rl_sinlon_dst = SIN(rl_plon)
629	d_dist = rl_coslat_dstCOS(rla_src_lats(ib_i)) &
630	(rl_coslon_dst*COS(rla_src_lons(ib_i,1)) + &
631	rl_sinlon_dst*SIN(rla_src_lons(ib_i,1)))+&
632	rl_sinlat_dst*SIN(rla_src_lats(ib_i))
633	IF (d_dist < -1.0d0) THEN
634	d_dist = -1.0d0
635	ELSE IF (d_dist > 1.0d0) THEN
636	d_dist = 1.0d0
637	END IF
638	rl_distance= ACOS(d_dist)
639	ELSE
640	rl_distance=bignum
641	END IF
642
643	IF (rl_distance < rla_dist(1)) THEN
644	ila_add_dist(2)=ila_add_dist(1)
645	ila_add_dist(1)=ib_i
646	rla_dist(2)=rla_dist(1)
647	rla_dist(1)=rl_distance
648	ELSE IF (rl_distance < rla_dist(2)) THEN
649	ila_add_dist(2)=ib_i
650	rla_dist(2)=rl_distance
651	END IF
652	END DO
653
654	IF (ila_nbr_found(ila_add_dist(1))>1 .and. &
655	ila_nbr_found(ila_add_dist(2))>1) THEN
656	! linearie
657	ll_linear=.true.
658	ELSE
659	! do distance weighted averege
660	rla_wght_lon(:,:)=0
661	DO ib_i=1,2
662	SELECT CASE (ila_nbr_found(ila_add_dist(ib_i)))
663	CASE (4)
664	CALL calcul_wght_irreg(rla_src_lons(ila_add_dist(ib_i),:),&
665	rl_plon, rla_wght_lon(ila_add_dist(ib_i),:))
666	rla_wght_lon(ila_add_dist(ib_i),:) = rla_wght_lon(ila_add_dist(ib_i),:)/&
667	rla_dist(ib_i)
668	CASE (3)
669	CALL calcul_wght_3(rla_src_lons(ila_add_dist(ib_i),1:3),&
670	rl_plon, rla_wght_lon(ila_add_dist(ib_i),1:3))
671	rla_wght_lon(ila_add_dist(ib_i),1:3) = rla_wght_lon(ila_add_dist(ib_i),1:3)/&
672	rla_dist(ib_i)
673	CASE (2)
674	CALL calcul_wght_2(rla_src_lons(ila_add_dist(ib_i),1:2),&
675	rl_plon, rla_wght_lon(ila_add_dist(ib_i),1:2))
676	rla_wght_lon(ila_add_dist(ib_i),1:2) = rla_wght_lon(ila_add_dist(ib_i),1:2)/&
677	rla_dist(ib_i)
678	CASE (1)
679	rla_wght_lon(ila_add_dist(ib_i),1) = 1/rla_dist(ib_i)
680	END SELECT
681	END DO
682
683	rl_distance=0
684	DO ib_i=1,4
685	rl_distance=rl_distance + sum(rla_wght_lon(ib_i,:))
686	END DO
687
688	rla_weight(:,:)=rla_wght_lon(:,:)/rl_distance
689
690	CALL store_link_bicub(ib_dst_add, ila_src_add , rla_weight, ib_thread)
691	#ifdef REMAP_TIMING
692	if (ll_timing) then
693	if (il_nbthreads.gt.1) then
694	!$ dl_tstop = OMP_GET_WTIME()
695	dla_timer(ib_thread,4) = dla_timer(ib_thread,4) + dl_tstop - dl_tstart
696	else
697	call timer_stop(6)
698	end if
699	end if
700	#endif
701	CYCLE grid_loop1
702	END IF
703	END IF
704
705	!
706	! Calculation of weights for longitudes
707	!
708
709	rla_wght_lon(:,:)=0
710	DO ib_i=1,4
711	SELECT CASE (ila_nbr_found(ib_i))
712	CASE (4)
713	CALL calcul_wght_irreg(rla_src_lons(ib_i,:), rl_plon, rla_wght_lon(ib_i,:))
714	CASE (3)
715	CALL calcul_wght_3(rla_src_lons(ib_i,1:3), rl_plon, rla_wght_lon(ib_i,1:3))
716	CASE (2)
717	CALL calcul_wght_2(rla_src_lons(ib_i,1:2), rl_plon, rla_wght_lon(ib_i,1:2))
718	END SELECT
719	END DO
720
721	IF (ll_linear) THEN
722	rla_wght_lat(:)=0
723	CALL calcul_wght_2(rla_src_lats(ila_add_dist(:)), rl_plat, rla_wght_temp(1:2))
724	rla_wght_lat(ila_add_dist(1))=rla_wght_temp(1)
725	rla_wght_lat(ila_add_dist(2))=rla_wght_temp(2)
726
727	DO ib_i=1,4
728	rla_weight(ib_i,:)=rla_wght_lat(ib_i)*rla_wght_lon(ib_i,:)
729	END DO
730
731	CALL store_link_bicub(ib_dst_add, ila_src_add , rla_weight, ib_thread)
732	#ifdef REMAP_TIMING
733	if (ll_timing) then
734	if (il_nbthreads.gt.1) then
735	!$ dl_tstop = OMP_GET_WTIME()
736	dla_timer(ib_thread,4) = dla_timer(ib_thread,4) + dl_tstop - dl_tstart
737	else
738	call timer_stop(6)
739	end if
740	end if
741	#endif
742	CYCLE grid_loop1
743	END IF
744
745	!
746	! Calculation of weights for latitudes
747	!
748	il_count=0
749	DO ib_i=1,4
750	IF (ila_nbr_found(ib_i)/=0) THEN
751	il_count=il_count+1
752	rla_lats_temp(il_count)=rla_src_lats(ib_i)
753	END IF
754	END DO
755
756	SELECT CASE (il_count)
757	CASE (4)
758	CALL calcul_wght_irreg(rla_lats_temp, rl_plat, rla_wght_temp(:))
759	CASE (3)
760	CALL calcul_wght_3(rla_lats_temp(1:3), rl_plat, rla_wght_temp(1:3))
761	CASE (2)
762	CALL calcul_wght_2(rla_lats_temp(1:2), rl_plat, rla_wght_temp(1:2))
763	CASE (1)
764	rla_wght_temp(1)=1
765	END SELECT
766
767	il_count=0
768	DO ib_i=1,4
769	IF (ila_nbr_found(ib_i)/=0) THEN
770	il_count=il_count+1
771	rla_wght_lat(ib_i)=rla_wght_temp(il_count)
772	ELSE
773	rla_wght_lat(ib_i)=0
774	END IF
775	END DO
776	!
777	! Calculation of total weight, elementwise multiplication
778	!
779
780	DO ib_i=1,4
781	rla_weight(ib_i,:)=rla_wght_lat(ib_i)*rla_wght_lon(ib_i,:)
782	END DO
783
784	CALL store_link_bicub(ib_dst_add, ila_src_add , rla_weight, ib_thread)
785	#ifdef REMAP_TIMING
786	if (ll_timing) then
787	if (il_nbthreads.gt.1) then
788	!$ dl_tstop = OMP_GET_WTIME()
789	dla_timer(ib_thread,4) = dla_timer(ib_thread,4) + dl_tstop - dl_tstart
790	else
791	call timer_stop(6)
792	end if
793	end if
794	#endif
795
796	END DO grid_loop1
797	end do thread_loop
798	!$OMP END DO
799
800	!$OMP END PARALLEL
801
802	if (il_nbthreads .gt. 1) then
803	#ifdef REMAP_TIMING
804	if (ll_timing) call timer_start(3,'remap_bicubic_reduced gather_lk')
805	#endif
806	sga_remap(1)%start_pos = 1
807	il_splitsize = sga_remap(1)%num_links
808	do ib_thread = 2, il_nbthreads
809	il_splitsize = il_splitsize + sga_remap(ib_thread)%num_links
810	sga_remap(ib_thread)%start_pos = sga_remap(ib_thread-1)%start_pos + &
811	sga_remap(ib_thread-1)%num_links
812	end do
813
814	num_links_map1 = il_splitsize
815	if (num_links_map1 > max_links_map1) &
816	call resize_remap_vars(1,num_links_map1-max_links_map1)
817
818	do ib_thread = 1, il_nbthreads
819	grid1_add_map1(sga_remap(ib_thread)%start_pos: &
820	sga_remap(ib_thread)%start_pos+ &
821	sga_remap(ib_thread)%num_links-1) = &
822	sga_remap(ib_thread)%grid1_add
823	grid2_add_map1(sga_remap(ib_thread)%start_pos: &
824	sga_remap(ib_thread)%start_pos+ &
825	sga_remap(ib_thread)%num_links-1) = &
826	sga_remap(ib_thread)%grid2_add
827	wts_map1 (:,sga_remap(ib_thread)%start_pos: &
828	sga_remap(ib_thread)%start_pos+ &
829	sga_remap(ib_thread)%num_links-1) = &
830	sga_remap(ib_thread)%wts
831
832	end do
833
834	#ifdef REMAP_TIMING
835	if (ll_timing) call timer_stop(3)
836	#endif
837	if (nlogprt.ge.2) then
838
839	do ib_thread = 1, il_nbthreads
840	if (sga_remap(ib_thread)%nb_resize.gt.0) then
841	write(nulou,*) ' Number of thread_resize_remap_vars on thread ',&
842	ib_thread, ' = ', sga_remap(ib_thread)%nb_resize
843	end if
844	end do
845
846	end if
847	#ifdef REMAP_TIMING
848	if (ll_timing.and.nlogprt.ge.2) then
849	write(nulou,*) ' On master thread '
850	call timer_print(2)
851	call timer_clear(2)
852	do ib_thread = 1,il_nbthreads
853	write(nulou,*) ' On thread ',ib_thread
854	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
855	& 'remap_bicubic_reduced distr ',dla_timer(ib_thread,1)
856	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
857	& 'remap_bicubic_reduced search ',dla_timer(ib_thread,2)
858	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
859	& 'remap_bicubic_reduced partial neighbours',dla_timer(ib_thread,3)
860	write(nulou,"(' Elapsed time for timer ',A24,':',1x,f10.4)")&
861	& 'remap_bicubic_reduced weights',dla_timer(ib_thread,4)
862	end do
863	deallocate(dla_timer)
864	write(nulou,*) ' On master thread '
865	call timer_print(3)
866	call timer_clear(3)
867	end if
868
869	else
870
871	if (ll_timing.and.nlogprt.ge.2) then
872	do n = 2, 6
873	call timer_print(n)
874	call timer_clear(n)
875	end do
876	end if
877	#endif
878	end if
879
880
881	! Gather the complete results on master proc
882
883	if (mpi_size_map .gt. 1) then
884
885	IF (mpi_rank_map == mpi_root_map) THEN
886	ALLOCATE(ila_num_links_mpi(mpi_size_map))
887	ELSE
888	ALLOCATE(ila_num_links_mpi(1))
889	END IF
890
891	CALL MPI_Gather (num_links_map1, 1,MPI_INT,&
892	& ila_num_links_mpi,1,MPI_INT,&
893	& mpi_root_map,mpi_comm_map,il_err)
894
895
896	IF (mpi_rank_map == mpi_root_map) THEN
897	num_links_map1 = SUM(ila_num_links_mpi)
898	if (num_links_map1 > max_links_map1) &
899	call resize_remap_vars(1,num_links_map1-max_links_map1)
900
901	ALLOCATE(ila_req_mpi(4,mpi_size_map-1))
902	ALLOCATE(ila_sta_mpi(MPI_STATUS_SIZE,4,mpi_size_map-1))
903
904	DO ib_i = 1, mpi_size_map-1
905
906	buff_base = SUM( ila_num_links_mpi(1:ib_i) ) + 1
907
908	CALL MPI_IRecv(grid1_add_map1(buff_base),&
909	& ila_num_links_mpi(ib_i+1),MPI_INT,ib_i,1,mpi_comm_map,&
910	& ila_req_mpi(1,ib_i),il_err)
911
912	CALL MPI_IRecv(grid2_add_map1(buff_base),&
913	& ila_num_links_mpi(ib_i+1),MPI_INT,ib_i,2,mpi_comm_map,&
914	& ila_req_mpi(2,ib_i),il_err)
915
916	CALL MPI_IRecv(wts_map1(:,buff_base),&
917	& num_wts*ila_num_links_mpi(ib_i+1),MPI_DOUBLE,ib_i,0,mpi_comm_map,&
918	& ila_req_mpi(3,ib_i),il_err)
919
920	CALL MPI_IRecv(grid2_frac(ila_mpi_mn(ib_i+1)),&
921	& ila_mpi_mx(ib_i+1)-ila_mpi_mn(ib_i+1)+1,MPI_DOUBLE,ib_i,0,mpi_comm_map,&
922	& ila_req_mpi(4,ib_i),il_err)
923
924	END DO
925
926	DO ib_i=1,4
927	CALL MPI_Waitall(mpi_size_map-1,ila_req_mpi(ib_i,:),ila_sta_mpi(1,ib_i,1),il_err)
928	END DO
929
930	DEALLOCATE(ila_req_mpi)
931	DEALLOCATE(ila_sta_mpi)
932
933	ELSE
934
935	CALL MPI_Send(grid1_add_map1,num_links_map1,MPI_INT,&
936	& mpi_root_map,1,mpi_comm_map,il_err)
937
938	CALL MPI_Send(grid2_add_map1,num_links_map1,MPI_INT,&
939	& mpi_root_map,2,mpi_comm_map,il_err)
940
941	CALL MPI_Send(wts_map1,num_wts*num_links_map1,MPI_DOUBLE,&
942	& mpi_root_map,0,mpi_comm_map,il_err)
943
944	CALL MPI_Send(grid2_frac(ila_mpi_mn(mpi_rank_map+1)),&
945	& ila_mpi_mx(mpi_rank_map+1)-ila_mpi_mn(mpi_rank_map+1)+1,MPI_DOUBLE,&
946	& mpi_root_map,0,mpi_comm_map,il_err)
947
948	END IF
949
950	deallocate(ila_num_links_mpi)
951
952	end if
953
954
955	deallocate (coslat, coslon, sinlat, sinlon)
956	deallocate(ila_mpi_mn)
957	deallocate(ila_mpi_mx)
958	deallocate(ila_thr_mn)
959	deallocate(ila_thr_mx)
960	if (il_nbthreads .gt. 1) then
961	do ib_thread = 1, il_nbthreads
962	deallocate(sga_remap(ib_thread)%grid1_add)
963	deallocate(sga_remap(ib_thread)%grid2_add)
964	deallocate(sga_remap(ib_thread)%wts)
965	end do
966	deallocate(sga_remap)
967	end if
968	!
969	IF (nlogprt .GE. 2) THEN
970	WRITE (UNIT = nulou,FMT = *)' '
971	WRITE (UNIT = nulou,FMT = *) 'Leaving routine remap_bicub_reduced'
972	CALL OASIS_FLUSH_SCRIP(nulou)
973	ENDIF
974	!
975	END SUBROUTINE remap_bicub_reduced
976
977
978	!-----------------------------------------------------------------------
979	! BOP
980	!
981	! !IROUTINE: grid_search_16_points
982	!
983	! !INTERFACE:
984	!
985	SUBROUTINE grid_search_16_points(ida_src_add, rda_src_lats, rda_src_lons,&
986	ida_nbr_found, bin, rd_plat, &
987	rd_plon, ld_extrapdone)
988	!
989	! !USES:
990	!
991	! !RETURN VALUE:
992	!
993	INTEGER (KIND=int_kind), DIMENSION(4,4), INTENT(out) :: &
994	ida_src_add ! searched addresses of the unmasked points enclosing
995	! target point
996
997	REAL (KIND=dbl_kind), DIMENSION(4,4), INTENT(out) :: &
998	rda_src_lons ! longitudes of the searched points
999
1000	REAL (KIND=dbl_kind), DIMENSION(4), INTENT(out) :: &
1001	rda_src_lats ! latitudes of the searched points
1002
1003	INTEGER (KIND=int_kind), DIMENSION(4), INTENT(out) :: &
1004	ida_nbr_found ! indicates for each line how many points found
1005
1006	INTEGER (KIND=int_kind), INTENT(out) :: &
1007	bin ! actual bin for target point
1008	!
1009	! !PARAMETERS:
1010	!
1011	REAL (KIND=dbl_kind), INTENT(in) :: &
1012	rd_plat, & ! latitude of the target point
1013	rd_plon ! longitude of the target point
1014
1015	LOGICAL, INTENT(in) :: ld_extrapdone ! true if extrapolation done
1016
1017	INTEGER (KIND=int_kind) :: &
1018	ib_k, ib_j, ib_i, & ! iteration indices
1019	il_min, il_max, il_inter ! begin and end for actual bin
1020
1021	INTEGER (KIND=int_kind), DIMENSION(4,2) :: &
1022	ila_corners ! temp addresses for bins
1023
1024	!
1025	! !DESCRIPTION:
1026	! This routine finds the location of the target point in the source
1027	! grid and returns those of the 16 nearest points that are unmasked.
1028	! The source grid is a reduced grid.
1029	!
1030	! !REVISION HISTORY:
1031	! 2002.10.21 J.Ghattas created
1032	!
1033	! EOP
1034	!-----------------------------------------------------------------------
1035	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
1036	! $Author: valcke $
1037	!-----------------------------------------------------------------------
1038
1039
1040	!
1041	! serch of actual bin
1042	!
1043
1044
1045	IF (rd_plat > bin_lats_r(1,1)) THEN ! norther of the first bin
1046	bin=0
1047	ila_corners(1:2,1:2)= 0
1048	ila_corners(3,1)= bin_addr1_r(1,1)+1
1049	ila_corners(3,2)= bin_addr1_r(2,1)
1050	ila_corners(4,1)= bin_addr1_r(1,2)
1051	ila_corners(4,2)= bin_addr1_r(2,2)
1052
1053	ELSE IF (rd_plat > bin_lats_r(1,2)) THEN ! in the first bin
1054	bin=1
1055	ila_corners(1,1:2)= 0
1056	ila_corners(2,1)= bin_addr1_r(1,1)+1
1057	ila_corners(2,2)= bin_addr1_r(2,1)
1058	ila_corners(3,1)= bin_addr1_r(1,2)
1059	ila_corners(3,2)= bin_addr1_r(2,2)
1060	ila_corners(4,1)= bin_addr1_r(1,3)
1061	ila_corners(4,2)= bin_addr1_r(2,3)
1062
1063	ELSE IF (rd_plat < bin_lats_r(1,num_srch_red)) THEN
1064	! South of the last complet bin
1065	bin=num_srch_red
1066	ila_corners(1,1) = bin_addr1_r(1,num_srch_red-1)
1067	ila_corners(1,2) = bin_addr1_r(2,num_srch_red-1)
1068	ila_corners(2,1) = bin_addr1_r(1,num_srch_red)
1069	ila_corners(2,2) = bin_addr1_r(2,num_srch_red)
1070	ila_corners(3:4,1:2) = 0
1071
1072	ELSE IF (rd_plat < bin_lats_r(1,num_srch_red-1)) THEN
1073	! in the last bin which is complet
1074	! the bin (num_srch_red-1) is the last bin which is complet
1075	bin=num_srch_red-1
1076	ila_corners(1,1) = bin_addr1_r(1,num_srch_red-2)
1077	ila_corners(1,2) = bin_addr1_r(2,num_srch_red-2)
1078	ila_corners(2,1) = bin_addr1_r(1,num_srch_red-1)
1079	ila_corners(2,2) = bin_addr1_r(2,num_srch_red-1)
1080	ila_corners(3,1) = bin_addr1_r(1,num_srch_red)
1081	ila_corners(3,2) = bin_addr1_r(2,num_srch_red)
1082	ila_corners(4,1:2) = 0
1083	ELSE
1084	il_min=2
1085	il_max=num_srch_red-1
1086	DO WHILE (il_min /= il_max-1)
1087	il_inter=(il_max-il_min)/2 + il_min
1088	IF (rd_plat <= bin_lats_r(1,il_min) .and. &
1089	rd_plat > bin_lats_r(1,il_inter)) THEN
1090	il_max=il_inter
1091	ELSE
1092	il_min=il_inter
1093	END IF
1094	END DO
1095	bin=il_min
1096	ila_corners(1,1) = bin_addr1_r(1,bin-1)
1097	ila_corners(1,2) = bin_addr1_r(2,bin-1)
1098	ila_corners(2,1) = bin_addr1_r(1,bin)
1099	ila_corners(2,2) = bin_addr1_r(2,bin)
1100	ila_corners(3,1) = bin_addr1_r(1,bin+1)
1101	ila_corners(3,2) = bin_addr1_r(2,bin+1)
1102	ila_corners(4,1) = bin_addr1_r(1,bin+2)
1103	ila_corners(4,2) = bin_addr1_r(2,bin+2)
1104
1105	IF (ila_corners(1,1)==0) THEN
1106	ila_corners(1,1)=1
1107	END IF
1108	END IF
1109
1110	DO ib_k=1,4
1111	IF (ila_corners(ib_k,1) .NE. 0) &
1112	rda_src_lats(ib_k)= grid1_center_lat(ila_corners(ib_k,1))
1113	ENDDO
1114
1115	!
1116	! now perform a more detailed search for each line
1117	!
1118
1119	ida_src_add(:,:)=0
1120	ida_nbr_found(:)=0
1121	rda_src_lons(:,:)=0
1122
1123	DO ib_k=1,4 ! for each line of found points
1124	IF (ila_corners(ib_k,1)==0) THEN
1125	CYCLE
1126	END IF
1127
1128	il_min=ila_corners(ib_k,1)
1129	il_max=ila_corners(ib_k,2)
1130
1131	IF (rd_plon < grid1_center_lon(il_min)) THEN
1132	DO ib_j=il_max-1, il_max
1133	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
1134	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
1135	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
1136	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
1137	grid1_center_lon(ib_j)-pi2
1138	END IF
1139	END DO
1140	DO ib_j=il_min, il_min+1
1141	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
1142	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
1143	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
1144	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
1145	grid1_center_lon(ib_j)
1146	END IF
1147	END DO
1148
1149	ELSE IF (rd_plon < grid1_center_lon(il_min+1)) THEN
1150	IF (grid1_mask(il_max) .or. ld_extrapdone) THEN
1151	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
1152	ida_src_add(ib_k,ida_nbr_found(ib_k)) = il_max
1153	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
1154	grid1_center_lon(il_max)-pi2
1155	END IF
1156	DO ib_j=il_min, il_min+2
1157	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
1158	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
1159	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
1160	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
1161	grid1_center_lon(ib_j)
1162	END IF
1163	END DO
1164
1165	ELSE IF (rd_plon > grid1_center_lon(il_max)) THEN
1166	DO ib_j=il_max-1, il_max
1167	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
1168	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
1169	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
1170	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
1171	grid1_center_lon(ib_j)
1172	END IF
1173	END DO
1174	DO ib_j=il_min, il_min+1
1175	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
1176	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
1177	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
1178	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
1179	grid1_center_lon(ib_j)+pi2
1180	END IF
1181	END DO
1182
1183	ELSE IF (rd_plon > grid1_center_lon(il_max-1)) THEN
1184	DO ib_j=il_max-2, il_max
1185	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
1186	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
1187	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
1188	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
1189	grid1_center_lon(ib_j)
1190	END IF
1191	END DO
1192	IF (grid1_mask(il_min) .or. ld_extrapdone) THEN
1193	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
1194	ida_src_add(ib_k,ida_nbr_found(ib_k)) = il_min
1195	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
1196	grid1_center_lon(il_min)+pi2
1197	END IF
1198
1199	ELSE
1200
1201	DO WHILE (il_min/=il_max-1)
1202	il_inter=(il_max-il_min)/2 + il_min
1203	IF (rd_plon >= grid1_center_lon(il_min) .and. &
1204	rd_plon < grid1_center_lon(il_inter)) THEN
1205	il_max=il_inter
1206	ELSE
1207	il_min=il_inter
1208	END IF
1209	END DO
1210	DO ib_i= il_min-1, il_min+2
1211	IF (grid1_mask(ib_i) .or. ld_extrapdone) THEN
1212	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
1213	ida_src_add(ib_k,ida_nbr_found(ib_k))=ib_i
1214	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
1215	grid1_center_lon(ib_i)
1216	END IF
1217	END DO
1218
1219	END IF
1220
1221	END DO
1222
1223
1224	END SUBROUTINE grid_search_16_points
1225
1226
1227
1228	!-----------------------------------------------------------------------
1229	! BOP
1230	!
1231	! !IROUTINE: calcul_wght_irreg
1232	!
1233	! !INTERFACE:
1234	!
1235	SUBROUTINE calcul_wght_irreg(rda_x, rd_pt, rda_wght)
1236	!
1237	! !USES:
1238	!
1239	! !RETURN VALUE:
1240	!
1241	REAL (KIND=dbl_kind), DIMENSION(4), INTENT(out) :: &
1242	rda_wght ! array of weights for the points x
1243	!
1244	! !PARAMETERS:
1245	!
1246	REAL (KIND=dbl_kind), DIMENSION(4), INTENT(in) :: &
1247	rda_x ! array of positions on source grid, lat or lon
1248
1249	REAL (KIND=dbl_kind),INTENT(in) :: &
1250	rd_pt ! position of target point to interpolate
1251
1252	REAL (KIND=dbl_kind) :: &
1253	rl_t1, rl_t2, rl_t3, rl_t4, rl_t5, rl_t6, rl_t7, rl_t8, rl_t9, &
1254	rl_u1, rl_u2, rl_u3, rl_u4, &
1255	rl_k1, rl_k2, rl_k3, &
1256	rl_d1, rl_d2, rl_d3, rl_d4, &
1257	rl_c1, rl_c2, rl_c3, rl_c4, &
1258	rl_b1, rl_b2, rl_b3, rl_b4, &
1259	rl_a1, rl_a2, rl_a3, rl_a4, &
1260	rl_y1, rl_y2, rl_y3, &
1261	rl_a1_y, rl_a2_y, rl_a3_y, rl_a4_y, &
1262	rl_b1_y, rl_b2_y, rl_b3_y, rl_b4_y, &
1263	rl_c1_y, rl_c2_y, rl_c3_y, rl_c4_y
1264	!
1265	! !DESCRIPTION:
1266	! Calculates a the weights of four points for a bicubic interpolation.
1267	! The distances between the points can be irregulier.
1268	!
1269	! !REVISION HISTORY:
1270	! 2002.10.21 J.Ghattas created
1271	!
1272	! EOP
1273	!-----------------------------------------------------------------------
1274	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
1275	! $Author: valcke $
1276	!-----------------------------------------------------------------------
1277
1278
1279	IF (rda_x(1)/=0.and. rda_x(2)/=0 .and. rda_x(3)/=0 .and. rda_x(4)/=0) THEN
1280
1281	rl_t1 = 1/rda_x(1) - 1/rda_x(2)
1282	rl_t2 = 1/rda_x(1)2 - 1/rda_x(2)2
1283	rl_t3 = 1/rda_x(1)3 - 1/rda_x(2)3
1284	rl_t4 = 1/rda_x(1) - 1/rda_x(3)
1285	rl_t5 = 1/rda_x(1)2 - 1/rda_x(3)2
1286	rl_t6 = 1/rda_x(1)3 - 1/rda_x(3)3
1287	rl_t7 = 1/rda_x(1) - 1/rda_x(4)
1288	rl_t8 = 1/rda_x(1)2 - 1/rda_x(4)2
1289	rl_t9 = 1/rda_x(1)3 - 1/rda_x(4)3
1290
1291	rl_u1 = rl_t2/rl_t1 - rl_t5/rl_t4
1292	rl_u2 = rl_t3/rl_t1 - rl_t6/rl_t4
1293	rl_u3 = rl_t2/rl_t1 - rl_t8/rl_t7
1294	rl_u4 = rl_t3/rl_t1 - rl_t9/rl_t7
1295
1296	rl_k1 = (1/(rl_t1rl_u1)-1/(rl_t1rl_u3)) / (rl_u2/rl_u1-rl_u4/rl_u3)
1297	rl_k2 = -1/(rl_t4*rl_u1) / (rl_u2/rl_u1-rl_u4/rl_u3)
1298	rl_k3 = 1/(rl_t7*rl_u3) / (rl_u2/rl_u1-rl_u4/rl_u3)
1299
1300
1301	rl_d1=(rl_k1+rl_k2+rl_k3)/rda_x(1)**3
1302	rl_d2 = -rl_k1/rda_x(2)**3
1303	rl_d3 = -rl_k2/rda_x(3)**3
1304	rl_d4 = -rl_k3/rda_x(4)**3
1305
1306	rl_c1 = 1/rl_u1(1/(rl_t1rda_x(1)*3)-1/(rl_t4rda_x(1)**3)- &
1307	rl_u2*rl_d1)
1308	rl_c2 = 1/rl_u1(1/(-rl_t1rda_x(2)*3)-rl_u2rl_d2)
1309	rl_c3 = 1/rl_u1(1/(rl_t4rda_x(3)*3)-rl_u2rl_d3)
1310	rl_c4 = 1/rl_u1(-rl_u2rl_d4)
1311
1312	rl_b1 = 1/rl_t1/rda_x(1)*3-rl_t2/rl_t1rl_c1-rl_t3/rl_t1*rl_d1
1313	rl_b2 = -1/rl_t1/rda_x(2)*3-rl_t2/rl_t1rl_c2-rl_t3/rl_t1*rl_d2
1314	rl_b3 = -rl_t2/rl_t1rl_c3-rl_t3/rl_t1rl_d3
1315	rl_b4 = -rl_t2/rl_t1rl_c4-rl_t3/rl_t1rl_d4
1316
1317	rl_a1 = 1/rda_x(1)*3-1/rda_x(1)rl_b1-1/rda_x(1)*2rl_c1- &
1318	1/rda_x(1)*3rl_d1
1319	rl_a2 = -1/rda_x(1)rl_b2-1/rda_x(1)2rl_c2-1/rda_x(1)*3rl_d2
1320	rl_a3 = -1/rda_x(1)rl_b3-1/rda_x(1)2rl_c3-1/rda_x(1)*3rl_d3
1321	rl_a4 = -1/rda_x(1)rl_b4-1/rda_x(1)2rl_c4-1/rda_x(1)*3rl_d4
1322
1323	! Weights
1324	rda_wght(1) = rl_a1rd_pt3 + rl_b1rd_pt*2 + rl_c1rd_pt + rl_d1
1325	rda_wght(2) = rl_a2rd_pt3 + rl_b2rd_pt*2 + rl_c2rd_pt + rl_d2
1326	rda_wght(3) = rl_a3rd_pt3 + rl_b3rd_pt*2 + rl_c3rd_pt + rl_d3
1327	rda_wght(4) = rl_a4rd_pt3 + rl_b4rd_pt*2 + rl_c4rd_pt + rl_d4
1328
1329	ELSE ! there is one point = 0
1330
1331	rl_d1=0; rl_d2=0; rl_d3=0; rl_d4=0
1332
1333	! Transformation for each case
1334	IF (rda_x(1)==0) THEN
1335	rl_y1=rda_x(2); rl_y2=rda_x(3); rl_y3=rda_x(4)
1336	rl_d1=1
1337	ELSE IF (rda_x(2)==0) THEN
1338	rl_y1=rda_x(1); rl_y2=rda_x(3); rl_y3=rda_x(4)
1339	rl_d2=1
1340	ELSE IF (rda_x(3)==0) THEN
1341	rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(4)
1342	rl_d3=1
1343	ELSE
1344	rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(3)
1345	rl_d4=1
1346	END IF
1347
1348	! Solving the system
1349	rl_t1 = 1/rl_y1-1/rl_y2
1350	rl_t2 = 1/rl_y12-1/rl_y22
1351	rl_t3 = 1/rl_y1-1/rl_y3
1352	rl_t4 = 1/rl_y12-1/rl_y32
1353
1354	rl_c1_y =(1/rl_y13/rl_t1-1/rl_y13/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3)
1355	rl_c2_y = -1/rl_y2**3/rl_t1/(rl_t2/rl_t1-rl_t4/rl_t3)
1356	rl_c3_y = 1/rl_y3**3/rl_t3/(rl_t2/rl_t1-rl_t4/rl_t3)
1357	rl_c4_y=(-1/rl_y13/rl_t1+1/rl_y23/rl_t1+ &
1358	1/rl_y13/rl_t3-1/rl_y33/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3)
1359
1360	rl_b1_y = 1/rl_y1*3/rl_t1 - rl_c1_yrl_t2/rl_t1
1361	rl_b2_y = -1/rl_y2*3/rl_t1 - rl_c2_yrl_t2/rl_t1
1362	rl_b3_y = -rl_c3_y*rl_t2/rl_t1
1363	rl_b4_y = -1/rl_y13/rl_t1 + 1/rl_y23/rl_t1 - rl_c4_y*rl_t2/rl_t1
1364
1365	rl_a1_y = 1/rl_y13 - rl_b1_y/rl_y1 - rl_c1_y/rl_y12
1366	rl_a2_y = -rl_b2_y/rl_y1 - rl_c2_y/rl_y1**2
1367	rl_a3_y = -rl_b3_y/rl_y1 - rl_c3_y/rl_y1**2
1368	rl_a4_y = -1/rl_y13 - rl_b4_y/rl_y1 - rl_c4_y/rl_y12
1369
1370	! Retransformation
1371	IF (rda_x(1)==0) THEN
1372	rl_a1=rl_a4_y; rl_a2=rl_a1_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y
1373	rl_b1=rl_b4_y; rl_b2=rl_b1_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y
1374	rl_c1=rl_c4_y; rl_c2=rl_c1_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y
1375	ELSE IF (rda_x(2)==0) THEN
1376	rl_a1=rl_a1_y; rl_a2=rl_a4_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y
1377	rl_b1=rl_b1_y; rl_b2=rl_b4_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y
1378	rl_c1=rl_c1_y; rl_c2=rl_c4_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y
1379	ELSE IF (rda_x(3)==0) THEN
1380	rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a4_y; rl_a4=rl_a3_y
1381	rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b4_y; rl_b4=rl_b3_y
1382	rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c4_y; rl_c4=rl_c3_y
1383	ELSE
1384	rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a3_y; rl_a4=rl_a4_y
1385	rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b3_y; rl_b4=rl_b4_y
1386	rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c3_y; rl_c4=rl_c4_y
1387	END IF
1388
1389	! Weights
1390	rda_wght(1) = rl_a1rd_pt3 + rl_b1rd_pt*2 + rl_c1rd_pt +rl_d1
1391	rda_wght(2) = rl_a2rd_pt3 + rl_b2rd_pt*2 + rl_c2rd_pt +rl_d2
1392	rda_wght(3) = rl_a3rd_pt3 + rl_b3rd_pt*2 + rl_c3rd_pt +rl_d3
1393	rda_wght(4) = rl_a4rd_pt3 + rl_b4rd_pt*2 + rl_c4rd_pt +rl_d4
1394
1395	END IF
1396
1397
1398	END SUBROUTINE calcul_wght_irreg
1399
1400	!-----------------------------------------------------------------------
1401	! BOP
1402	!
1403	! !IROUTINE: calcul_wght_3
1404	!
1405	! !INTERFACE:
1406
1407	SUBROUTINE calcul_wght_3(rda_x, rd_pt, rda_wght)
1408
1409	! !USES:
1410
1411	! !RETURN VALUE:
1412
1413	REAL (KIND=dbl_kind), DIMENSION(3), INTENT(out) :: &
1414	rda_wght ! array of weights for the points x
1415
1416	! !PARAMETERS:
1417
1418	REAL (KIND=dbl_kind), DIMENSION(3), INTENT(in) :: &
1419	rda_x ! array of positions on source grid, lat or lon
1420
1421	REAL (KIND=dbl_kind), INTENT(in) :: &
1422	rd_pt ! position of target point to interpolate
1423
1424	REAL (KIND=dbl_kind) :: &
1425	rl_c1, rl_c2, rl_c3, &
1426	rl_a1, rl_a2, rl_a3, &
1427	rl_b1, rl_b2, rl_b3, &
1428	rl_t1, rl_t2, rl_t3, rl_t4
1429
1430	! !DESCRIPTION:
1431	! Calculates a the weights of 3 points for a parabolic interpolation.
1432	!
1433	! !REVISION HISTORY:
1434	! 2002.10.21 J.Ghattas created
1435	!
1436	! EOP
1437	!-----------------------------------------------------------------------
1438	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
1439	! $Author: valcke $
1440	!-----------------------------------------------------------------------
1441
1442
1443	IF (rda_x(1)/=0 .and. rda_x(2)/=0 .and. rda_x(3)/=0) THEN
1444	rl_t1 = 1/rda_x(1)-1/rda_x(2)
1445	rl_t2 = 1/rda_x(1)2-1/rda_x(2)2
1446	rl_t3 = 1/rda_x(1)-1/rda_x(3)
1447	rl_t4 = 1/rda_x(1)2-1/rda_x(3)2
1448
1449	rl_c1 = (1/rda_x(1)2/rl_t1-1/rda_x(1)2/rl_t3) / &
1450	(rl_t2/rl_t1-rl_t4/rl_t3)
1451	rl_c2 = -1/rda_x(2)**2/rl_t1 / (rl_t2/rl_t1-rl_t4/rl_t3)
1452	rl_c3 = 1/rda_x(3)**2/rl_t3 / (rl_t2/rl_t1-rl_t4/rl_t3)
1453
1454	rl_b1 = 1/rda_x(1)*2/rl_t1 - rl_c1rl_t2/rl_t1
1455	rl_b2 = -1/rda_x(2)*2/rl_t1 - rl_c2rl_t2/rl_t1
1456	rl_b3 = - rl_c3*rl_t2/rl_t1
1457
1458	rl_a1 = 1/rda_x(1)2 - rl_b1/rda_x(1) - rl_c1/rda_x(1)2
1459	rl_a2 = - rl_b2/rda_x(1) - rl_c2/rda_x(1)**2
1460	rl_a3 = - rl_b3/rda_x(1) - rl_c3/rda_x(1)**2
1461
1462
1463	ELSE IF (rda_x(1)==0) THEN
1464	rl_c1 = 1; rl_c2 = 0; rl_c3 = 0
1465	rl_b1 = (-1/rda_x(2)2+1/rda_x(3)2) / (1/rda_x(2)-1/rda_x(3))
1466	rl_b2 = 1/rda_x(2)**2 / (1/rda_x(2)-1/rda_x(3))
1467	rl_b3 = -1/rda_x(3)**2 / (1/rda_x(2)-1/rda_x(3))
1468
1469	rl_a1 = -1/rda_x(2)**2 - rl_b1/rda_x(2)
1470	rl_a2 = 1/rda_x(2)**2 - rl_b2/rda_x(2)
1471	rl_a3 = - rl_b3/rda_x(2)
1472
1473	ELSE IF (rda_x(2)==0) THEN
1474
1475	rl_c1 = 0; rl_c2 = 1; rl_c3 = 0
1476	rl_b1 = 1/rda_x(1)**2 / (1/rda_x(1)-1/rda_x(3))
1477	rl_b2 = (-1/rda_x(1)2+1/rda_x(3)2) / (1/rda_x(1)-1/rda_x(3))
1478	rl_b3 = -1/rda_x(3)**2 / (1/rda_x(1)-1/rda_x(3))
1479
1480	rl_a1 = 1/rda_x(1)**2 - rl_b1/rda_x(1)
1481	rl_a2 = -1/rda_x(1)**2 - rl_b2/rda_x(1)
1482	rl_a3 = - rl_b3/rda_x(1)
1483
1484	ELSE !rda_x(3)==0
1485	rl_c1 = 0; rl_c2 = 0; rl_c3 = 1
1486	rl_b1 = 1/rda_x(1)**2 / (1/rda_x(1)-1/rda_x(2))
1487	rl_b2 = -1/rda_x(2)**2 / (1/rda_x(1)-1/rda_x(2))
1488	rl_b3 = (-1/rda_x(1)2+1/rda_x(2)2) / (1/rda_x(1)-1/rda_x(2))
1489
1490	rl_a1 = 1/rda_x(1)**2 - rl_b1/rda_x(1)
1491	rl_a2 = - rl_b2/rda_x(1)
1492	rl_a3 = -1/rda_x(1)**2 - rl_b3/rda_x(1)
1493
1494
1495	END IF
1496
1497	! Weights
1498	rda_wght(1) = rl_a1rd_pt2 + rl_b1rd_pt + rl_c1
1499	rda_wght(2) = rl_a2rd_pt2 + rl_b2rd_pt + rl_c2
1500	rda_wght(3) = rl_a3rd_pt2 + rl_b3rd_pt + rl_c3
1501
1502
1503	END SUBROUTINE calcul_wght_3
1504
1505
1506	!-----------------------------------------------------------------------
1507	! BOP
1508	!
1509	! !IROUTINE: calcul_wght_2
1510	!
1511	! !INTERFACE:
1512
1513	SUBROUTINE calcul_wght_2(rda_x, rd_pt, rda_wght)
1514
1515	! !USES:
1516
1517	! !RETURN VALUE:
1518
1519	REAL (KIND=dbl_kind), DIMENSION(2), INTENT(out) :: &
1520	rda_wght ! array of weights for the points x
1521
1522	! !PARAMETERS:
1523
1524	REAL (KIND=dbl_kind), DIMENSION(2), INTENT(in) :: &
1525	rda_x ! array of positions on source grid, lat or lon
1526
1527	REAL (KIND=dbl_kind), INTENT(in) :: &
1528	rd_pt ! position of target point to interpolate
1529
1530	REAL (KIND=dbl_kind) :: rl_b1, rl_b2, rl_a1, rl_a2
1531
1532	! !DESCRIPTION:
1533	! Calculates a the weights of 2 points for a linair interpolation.
1534	!
1535	! !REVISION HISTORY:
1536	! 2002.10.21 J.Ghattas created
1537	!
1538	! EOP
1539	!-----------------------------------------------------------------------
1540	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
1541	! $Author: valcke $
1542	!-----------------------------------------------------------------------
1543
1544
1545	IF (rda_x(1)/=0 .and. rda_x(2)/=0) THEN
1546	rl_b1 = 1/(1-rda_x(1)/rda_x(2))
1547	rl_b2 = -1/(rda_x(2)/rda_x(1)-1)
1548	rl_a1 = 1/rda_x(1) - rl_b1/rda_x(1)
1549	rl_a2 = - rl_b2/rda_x(1)
1550
1551	ELSE IF (rda_x(1)==0) THEN
1552	rl_b1=1; rl_b2=0
1553	rl_a1=-1/rda_x(2)
1554	rl_a2=1/rda_x(2)
1555	ELSE
1556	rl_b1=0; rl_b2=1
1557	rl_a1=1/rda_x(1)
1558	rl_a2=-1/rda_x(1)
1559	END IF
1560
1561	rda_wght(1) = rl_a1*rd_pt + rl_b1
1562	rda_wght(2) = rl_a2*rd_pt + rl_b2
1563
1564	END SUBROUTINE calcul_wght_2
1565
1566
1567	!-----------------------------------------------------------------------
1568	! BOP
1569	!
1570	! !IROUTINE: store_link_bicub
1571	!
1572	! !INTERFACE:
1573
1574	SUBROUTINE store_link_bicub(id_dst_add, ida_src_add, rda_wght, id_thread)
1575
1576	! !USES:
1577
1578	! !RETURN VALUE:
1579
1580	! !PARAMETERS:
1581
1582	INTEGER (KIND=int_kind), INTENT(in) :: &
1583	id_dst_add ! address on destination grid
1584
1585	INTEGER (KIND=int_kind), DIMENSION(4,4), INTENT(in) :: &
1586	ida_src_add ! addresses for links on source grid
1587
1588	REAL (KIND=dbl_kind), DIMENSION(4,4), INTENT(in) :: &
1589	rda_wght ! array of remapping weights for these links
1590
1591	INTEGER (KIND=int_kind), INTENT(in) :: &
1592	id_thread ! local threaded task
1593
1594	INTEGER (KIND=int_kind) :: ib_i, ib_j, ib_ind, &
1595	il_num_links_old ! placeholder for old link number
1596
1597	!EM INTEGER (KIND=int_kind), DIMENSION(num_neighbors) :: &
1598	!EM ila_src_add ! reshaped addresses
1599
1600	!EM REAL (KIND=dbl_kind), DIMENSION(num_neighbors) :: &
1601	!EM rla_wght ! reshaped weights
1602
1603	! !DESCRIPTION:
1604	! This routine stores the addresses and weights for 16 links associated
1605	! with one destination point in the appropriate address.
1606	!
1607	! !REVISION HISTORY:
1608	! 2002.10.21 J.Ghattas created
1609	!
1610	! EOP
1611	!-----------------------------------------------------------------------
1612	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
1613	! $Author: valcke $
1614	!-----------------------------------------------------------------------
1615
1616
1617	!
1618	! Increment number of links and check if remap arrays need
1619	! to be increased to accomodate the new link. then store the link.
1620	!
1621
1622	if (il_nbthreads .eq. 1) then
1623
1624
1625	il_num_links_old = num_links_map1
1626	num_links_map1 = il_num_links_old + num_neighbors
1627
1628	IF (num_links_map1 > max_links_map1) THEN
1629	CALL resize_remap_vars(1,MAX(resize_increment,num_neighbors))
1630	END IF
1631
1632	!EM ila_src_add=RESHAPE(ida_src_add,(/num_neighbors/))
1633	!EM rla_wght=RESHAPE(rda_wght,(/num_neighbors/))
1634
1635	ib_ind = 0
1636	do ib_j=1,4
1637	do ib_i=1,4
1638	ib_ind = ib_ind + 1
1639	grid1_add_map1(il_num_links_old+ib_ind) = ida_src_add(ib_i,ib_j)
1640	grid2_add_map1(il_num_links_old+ib_ind) = id_dst_add
1641	wts_map1(1,il_num_links_old+ib_ind) = rda_wght(ib_i,ib_j)
1642	end do
1643	end do
1644
1645	else
1646
1647	sga_remap(id_thread)%num_links = sga_remap(id_thread)%num_links + num_neighbors
1648
1649	if (sga_remap(id_thread)%num_links > sga_remap(id_thread)%max_links) &
1650	call sga_remap(id_thread)%resize(int(0.2*real(sga_remap(id_thread)%max_links)))
1651
1652	ib_ind = 0
1653	do ib_j=1,4
1654	do ib_i=1,4
1655	ib_ind = ib_ind + 1
1656	sga_remap(id_thread)%grid1_add(sga_remap(id_thread)%num_links-num_neighbors+ib_ind) = ida_src_add(ib_i,ib_j)
1657	sga_remap(id_thread)%grid2_add(sga_remap(id_thread)%num_links-num_neighbors+ib_ind) = id_dst_add
1658	sga_remap(id_thread)%wts(1,sga_remap(id_thread)%num_links-num_neighbors+ib_ind) = rda_wght(ib_i,ib_j)
1659	end do
1660	end do
1661
1662	endif
1663
1664
1665	END SUBROUTINE store_link_bicub
1666
1667
1668	END MODULE remap_bicubic_reduced
1669
1670	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1671

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_bicubic_reduced.F90 @ 6331

Download in other formats: