Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

remap_bicubic_reduced.F90 @ 4775

Last change on this file since 4775 was 4775, checked in by aclsce, 5 years ago

Imported oasis3-mct from Cerfacs svn server (not suppotred anymore).

The version has been extracted from https://oasis3mct.cerfacs.fr/svn/branches/OASIS3-MCT_2.0_branch/oasis3-mct@1818

File size: 33.5 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 remap_vars ! module containing remap info
13	USE mod_oasis_flush
14
15	! !PUBLIC TYPES:
16	IMPLICIT NONE
17	! !PUBLIC MEMBER FUNCTIONS:
18	!
19	! !PUBLIC DATA MEMBERS:
20
21	! !DESCRIPTION:
22	! This routine computes the weights for a bicubic interpolation
23	! with a reduced grid. Computes mappings from grid1 to grid2.
24	!
25	! !REVISION HISTORY:
26	! 2002.10.21 J.Ghattas created
27	!
28	! EOP
29	!-----------------------------------------------------------------------
30	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
31	! $Author: valcke $
32	!-----------------------------------------------------------------------
33
34
35
36	CONTAINS
37
38	!***********************************************************************
39
40
41	!-----------------------------------------------------------------------
42	! BOP
43	!
44	! !IROUTINE: remap_bicub_reduced
45	!
46	! !INTERFACE:
47
48	SUBROUTINE remap_bicub_reduced(ld_extrapdone)
49
50	! !USES:
51
52	! !RETURN VALUE:
53
54	! !PARAMETERS:
55
56	LOGICAL, INTENT(in) :: &
57	ld_extrapdone ! logical, true if EXTRAP done on field
58	LOGICAL :: ll_nnei ! true (default) if extra search is done
59
60	INTEGER (KIND=int_kind), DIMENSION(4,4) :: &
61	ila_src_add ! address for source points non-masked
62
63	INTEGER (KIND=int_kind), DIMENSION(4) :: &
64	ila_nbr_found ! nrb of points found on each latitude
65
66	INTEGER (KIND=int_kind) :: &
67	ib_i, & ! iter index
68	ib_dst_add, & ! destination address, target point
69	il_count, & ! nbr of latitudes with found points
70	il_min, il_max, bin ! begin and end for distances calculation
71
72	REAL (KIND=dbl_kind), DIMENSION(4,4) :: &
73	rla_src_lons, & ! longitudes for the points 'ila_src_add'
74	rla_weight, & ! bicubic weights for the points 'ila_src_add'
75	rla_wght_lon ! temp. weights
76
77	REAL (KIND=dbl_kind), DIMENSION(4) :: &
78	rla_src_lats, & ! latitudes for the points 'ila_src_add'
79	rla_lats_temp, & ! temp. latitudes
80	rla_wght_lat, rla_wght_temp! temp. weights
81
82	REAL (KIND=dbl_kind) :: &
83	rl_plat, rl_plon ! latitude and longitude for destination address
84
85	REAL (KIND=dbl_kind) :: & ! variables for distances calculation
86	rl_coslat_dst, rl_sinlat_dst, &
87	rl_coslon_dst, rl_sinlon_dst, &
88	rl_distance, arg
89
90	REAL (KIND=dbl_kind), DIMENSION(2) :: &
91	rla_dist ! lat distances to point cible
92
93	INTEGER (KIND=int_kind), DIMENSION(4) :: &
94	ila_add_dist ! temporary addr. for distances
95
96	LOGICAL :: ll_linear ! flag
97
98
99	! !DESCRIPTION:
100	! This routine computes the weights for a bicubic interpolation
101	! with a reduced grid. Computes mappings from grid1 to grid2.
102	!
103	! !REVISION HISTORY:
104	! 2002.10.21 J.Ghattas created
105	!
106	! EOP
107	!-----------------------------------------------------------------------
108	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
109	! $Author: valcke $
110	!-----------------------------------------------------------------------
111	!
112	IF (nlogprt .GE. 2) THEN
113	WRITE (UNIT = nulou,FMT = *)' '
114	WRITE (UNIT = nulou,FMT = *) 'Entering routine remap_bicub_reduced'
115	CALL OASIS_FLUSH_SCRIP(nulou)
116	ENDIF
117	!
118	ll_nnei = .true.
119	!
120	! Loop over destination grid
121	!
122
123	DO ib_dst_add = 1, grid2_size ! for each target point
124	ll_linear=.false.
125	IF (.NOT. grid2_mask(ib_dst_add)) THEN
126	CYCLE ! target point is masked
127	END IF
128
129	rl_plat = grid2_center_lat(ib_dst_add)
130	rl_plon = grid2_center_lon(ib_dst_add)
131
132
133	!
134	! Search for non-masked points among the closest 16 points
135	! on source grid (grid1)
136	!
137
138	CALL grid_search_16_points(ila_src_add, rla_src_lats, rla_src_lons,&
139	ila_nbr_found, bin, rl_plat, &
140	rl_plon, ld_extrapdone)
141
142
143	!
144	! If there is no point found, search the neaerst
145	! non-masked point
146	!
147
148	IF (SUM(ila_nbr_found)==0) THEN
149	IF (ll_nnei .EQV. .TRUE. ) THEN
150	IF (nlogprt .GE. 2) THEN
151	WRITE(nulou,*) ' '
152	WRITE(nulou,*) &
153	'All 16 surrounding source grid points are masked'
154	WRITE(nulou,*) 'for target point ',ib_dst_add
155	WRITE(nulou,*) 'with longitude and latitude', rl_plon, rl_plat
156	WRITE(nulou,*) 'Using the nearest non-masked neighbour.'
157	WRITE(nulou,*) ' '
158	CALL OASIS_FLUSH_SCRIP(nulou)
159	ENDIF
160
161	! Search the nearest point in bin [il_min:il_max]
162	IF (bin==0 .or. bin==1) THEN
163	il_min=1
164	il_max=bin_addr1_r(2,3)
165	ELSE IF (bin==num_srch_red .or. bin==num_srch_red-1) THEN
166	il_min=bin_addr1_r(1,num_srch_red-2)
167	il_max=bin_addr1_r(2,num_srch_red)
168	ELSE
169	il_min=bin_addr1_r(1,bin-1)+1
170	il_max=bin_addr1_r(2,bin+2)
171	END IF
172
173	rl_coslat_dst = COS(rl_plat)
174	rl_sinlat_dst = SIN(rl_plat)
175	rl_coslon_dst = COS(rl_plon)
176	rl_sinlon_dst = SIN(rl_plon)
177
178	rla_weight(1,1) = bignum
179	ila_src_add(1,1) = 0
180	!cdir novector
181	DO ib_i=il_min, il_max
182	IF (grid1_mask(ib_i) .or. ld_extrapdone) THEN
183	arg = rl_coslat_dstCOS(grid1_center_lat(ib_i)) &
184	(rl_coslon_dst*COS(grid1_center_lon(ib_i)) + &
185	rl_sinlon_dst*SIN(grid1_center_lon(ib_i)))+&
186	rl_sinlat_dst*SIN(grid1_center_lat(ib_i))
187	IF (arg < -1.0d0) THEN
188	arg = -1.0d0
189	ELSE IF (arg > 1.0d0) THEN
190	arg = 1.0d0
191	END IF
192	rl_distance = ACOS(arg)
193	IF (rl_distance < rla_weight(1,1)) THEN
194	rla_weight(1,1) = rl_distance
195	ila_src_add(1,1) = ib_i
196	END IF
197	END IF
198	END DO
199	rla_weight(:,:) = 0
200	rla_weight(1,1) = 1
201
202	CALL store_link_bicub(ib_dst_add, ila_src_add, rla_weight)
203	IF (nlogprt .GE. 2) THEN
204	WRITE(nulou,*) &
205	'Nearest non masked neighbour is source point ', &
206	ila_src_add(1,1)
207	WRITE(nulou,*) 'with longitude and latitude', &
208	grid1_center_lon(ila_src_add(1,1)), &
209	grid1_center_lat(ila_src_add(1,1))
210	WRITE(nulou,*) ' '
211	ENDIF
212	CYCLE
213	ENDIF
214	END IF
215
216	rla_weight(:,:) = 0
217	! if there is only one point found, save it
218	IF (SUM(ila_nbr_found)==1) THEN
219	DO ib_i=1,4
220	IF (ila_nbr_found(ib_i)==1) THEN
221	rla_weight(ib_i,1)=1
222	EXIT
223	END IF
224	END DO
225	CALL store_link_bicub(ib_dst_add, ila_src_add, rla_weight)
226	CYCLE
227	END IF
228
229	! if there are only 2 points found, distance weighted average
230	IF (SUM(ila_nbr_found)==2) THEN
231	rl_coslat_dst = COS(rl_plat)
232	rl_sinlat_dst = SIN(rl_plat)
233	rl_coslon_dst = COS(rl_plon)
234	rl_sinlon_dst = SIN(rl_plon)
235
236	rl_distance=0 ! count of total distance
237	DO ib_i=1,4
238	IF (ila_nbr_found(ib_i) > 0) THEN
239	arg = rl_coslat_dstCOS(rla_src_lats(ib_i)) &
240	(rl_coslon_dst*COS(rla_src_lons(ib_i,1)) + &
241	rl_sinlon_dst*SIN(rla_src_lons(ib_i,1)))+&
242	rl_sinlat_dst*SIN(rla_src_lats(ib_i))
243	IF (arg < -1.0d0) THEN
244	arg = -1.0d0
245	ELSE IF (arg > 1.0d0) THEN
246	arg = 1.0d0
247	END IF
248	rla_weight(ib_i,1) = ACOS(arg)
249	rl_distance = rl_distance+rla_weight(ib_i,1)
250	IF (ila_nbr_found(ib_i)==2) THEN
251	arg = rl_coslat_dstCOS(rla_src_lats(ib_i)) &
252	(rl_coslon_dst*COS(rla_src_lons(ib_i,2)) + &
253	rl_sinlon_dst*SIN(rla_src_lons(ib_i,2)))+&
254	rl_sinlat_dst*SIN(rla_src_lats(ib_i))
255	IF (arg < -1.0d0) THEN
256	arg = -1.0d0
257	ELSE IF (arg > 1.0d0) THEN
258	arg = 1.0d0
259	END IF
260	rla_weight(ib_i,2) = ACOS(arg)
261	rl_distance = rl_distance+rla_weight(ib_i,2)
262	END IF
263	END IF
264	END DO
265	rla_weight=rla_weight/rl_distance
266
267	CALL store_link_bicub(ib_dst_add, ila_src_add, rla_weight)
268	CYCLE
269	END IF
270
271	! Some case exceptional when just one point per line found
272
273	IF (ila_nbr_found(1)==1) THEN ! elimination of point
274	ila_nbr_found(1)=0
275	ila_src_add(1,1)=0
276	END IF
277	IF (ila_nbr_found(4)==1) THEN
278	ila_nbr_found(4)=0
279	ila_src_add(4,1)=0
280	END IF
281
282
283
284	IF (ila_nbr_found(2)==1 .or. ila_nbr_found(3)==1) THEN
285	ila_add_dist(:)=4
286	rla_dist(:)=bignum
287
288	! search for the 2 nearest points or line of points
289	DO ib_i=1,4
290	IF (ila_nbr_found(ib_i) > 1) THEN
291	rl_distance=ABS(rla_src_lats(ib_i)-rl_plat)
292	ELSE IF (ila_nbr_found(ib_i)==1) THEN
293	rl_coslat_dst = COS(rl_plat)
294	rl_sinlat_dst = SIN(rl_plat)
295	rl_coslon_dst = COS(rl_plon)
296	rl_sinlon_dst = SIN(rl_plon)
297	arg = rl_coslat_dstCOS(rla_src_lats(ib_i)) &
298	(rl_coslon_dst*COS(rla_src_lons(ib_i,1)) + &
299	rl_sinlon_dst*SIN(rla_src_lons(ib_i,1)))+&
300	rl_sinlat_dst*SIN(rla_src_lats(ib_i))
301	IF (arg < -1.0d0) THEN
302	arg = -1.0d0
303	ELSE IF (arg > 1.0d0) THEN
304	arg = 1.0d0
305	END IF
306	rl_distance= ACOS(arg)
307	ELSE
308	rl_distance=bignum
309	END IF
310
311	IF (rl_distance < rla_dist(1)) THEN
312	ila_add_dist(2)=ila_add_dist(1)
313	ila_add_dist(1)=ib_i
314	rla_dist(2)=rla_dist(1)
315	rla_dist(1)=rl_distance
316	ELSE IF (rl_distance < rla_dist(2)) THEN
317	ila_add_dist(2)=ib_i
318	rla_dist(2)=rl_distance
319	END IF
320	END DO
321
322	IF (ila_nbr_found(ila_add_dist(1))>1 .and. &
323	ila_nbr_found(ila_add_dist(2))>1) THEN
324	! linearie
325	ll_linear=.true.
326	ELSE
327	! do distance weighted averege
328	rla_wght_lon(:,:)=0
329	DO ib_i=1,2
330	SELECT CASE (ila_nbr_found(ila_add_dist(ib_i)))
331	CASE (4)
332	CALL calcul_wght_irreg(rla_src_lons(ila_add_dist(ib_i),:),&
333	rl_plon, rla_wght_lon(ila_add_dist(ib_i),:))
334	rla_wght_lon(ila_add_dist(ib_i),:)=&
335	rla_wght_lon(ila_add_dist(ib_i),:)/&
336	rla_dist(ib_i)
337	CASE (3)
338	CALL calcul_wght_3(rla_src_lons(ila_add_dist(ib_i),1:3),&
339	rl_plon, rla_wght_lon(ila_add_dist(ib_i),1:3))
340	rla_wght_lon(ila_add_dist(ib_i),1:3)=&
341	rla_wght_lon(ila_add_dist(ib_i),1:3)/&
342	rla_dist(ib_i)
343	CASE (2)
344	CALL calcul_wght_2(rla_src_lons(ila_add_dist(ib_i),1:2),&
345	rl_plon, rla_wght_lon(ila_add_dist(ib_i),1:2))
346	rla_wght_lon(ila_add_dist(ib_i),1:2)=&
347	rla_wght_lon(ila_add_dist(ib_i),1:2)/&
348	rla_dist(ib_i)
349	CASE (1)
350	rla_wght_lon(ila_add_dist(ib_i),1)=1/rla_dist(ib_i)
351	END SELECT
352	END DO
353	rl_distance=0
354	DO ib_i=1,4
355	rl_distance=rl_distance + sum(rla_wght_lon(ib_i,:))
356	END DO
357	rla_weight(:,:)=rla_wght_lon(:,:)/rl_distance
358
359	CALL store_link_bicub(ib_dst_add, ila_src_add , rla_weight)
360	CYCLE
361	END IF
362	END IF
363
364	!
365	! Calculation of weights for longitudes
366	!
367
368	rla_wght_lon(:,:)=0
369	DO ib_i=1,4
370	SELECT CASE (ila_nbr_found(ib_i))
371	CASE (4)
372	CALL calcul_wght_irreg(rla_src_lons(ib_i,:), rl_plon, &
373	rla_wght_lon(ib_i,:))
374	CASE (3)
375	CALL calcul_wght_3(rla_src_lons(ib_i,1:3), rl_plon, &
376	rla_wght_lon(ib_i,1:3))
377	CASE (2)
378	CALL calcul_wght_2(rla_src_lons(ib_i,1:2), rl_plon, &
379	rla_wght_lon(ib_i,1:2))
380	END SELECT
381	END DO
382
383
384	IF (ll_linear) THEN
385	rla_wght_lat(:)=0
386	CALL calcul_wght_2(rla_src_lats(ila_add_dist(:)), rl_plat, &
387	rla_wght_temp(1:2))
388	rla_wght_lat(ila_add_dist(1))=rla_wght_temp(1)
389	rla_wght_lat(ila_add_dist(2))=rla_wght_temp(2)
390	DO ib_i=1,4
391	rla_weight(ib_i,:)=rla_wght_lat(ib_i)*rla_wght_lon(ib_i,:)
392	END DO
393
394	CALL store_link_bicub(ib_dst_add, ila_src_add , rla_weight)
395	CYCLE
396	END IF
397
398
399	!
400	! Calculation of weights for latitudes
401	!
402
403	il_count=0
404	DO ib_i=1,4
405	IF (ila_nbr_found(ib_i)/=0) THEN
406	il_count=il_count+1
407	rla_lats_temp(il_count)=rla_src_lats(ib_i)
408	END IF
409	END DO
410
411	SELECT CASE (il_count)
412	CASE (4)
413	CALL calcul_wght_irreg(rla_lats_temp, rl_plat, rla_wght_temp(:))
414	CASE (3)
415	CALL calcul_wght_3(rla_lats_temp(1:3), rl_plat, rla_wght_temp(1:3))
416	CASE (2)
417	CALL calcul_wght_2(rla_lats_temp(1:2), rl_plat, rla_wght_temp(1:2))
418	CASE (1)
419	rla_wght_temp(1)=1
420	END SELECT
421
422	il_count=0
423	DO ib_i=1,4
424	IF (ila_nbr_found(ib_i)/=0) THEN
425	il_count=il_count+1
426	rla_wght_lat(ib_i)=rla_wght_temp(il_count)
427	ELSE
428	rla_wght_lat(ib_i)=0
429	END IF
430	END DO
431
432	!
433	! Calculation of total weight, elementwise multiplication
434	!
435
436	DO ib_i=1,4
437	rla_weight(ib_i,:)=rla_wght_lat(ib_i)*rla_wght_lon(ib_i,:)
438	END DO
439
440	CALL store_link_bicub(ib_dst_add, ila_src_add , rla_weight)
441
442	END DO
443	!
444	IF (nlogprt .GE. 2) THEN
445	WRITE (UNIT = nulou,FMT = *)' '
446	WRITE (UNIT = nulou,FMT = *) 'Leaving routine remap_bicub_reduced'
447	CALL OASIS_FLUSH_SCRIP(nulou)
448	ENDIF
449	!
450	END SUBROUTINE remap_bicub_reduced
451
452
453	!-----------------------------------------------------------------------
454	! BOP
455	!
456	! !IROUTINE: grid_search_16_points
457	!
458	! !INTERFACE:
459	!
460	SUBROUTINE grid_search_16_points(ida_src_add, rda_src_lats, rda_src_lons,&
461	ida_nbr_found, bin, rd_plat, &
462	rd_plon, ld_extrapdone)
463	!
464	! !USES:
465	!
466	! !RETURN VALUE:
467	!
468	INTEGER (KIND=int_kind), DIMENSION(4,4), INTENT(out) :: &
469	ida_src_add ! searched addresses of the unmasked points enclosing
470	! target point
471
472	REAL (KIND=dbl_kind), DIMENSION(4,4), INTENT(out) :: &
473	rda_src_lons ! longitudes of the searched points
474
475	REAL (KIND=dbl_kind), DIMENSION(4), INTENT(out) :: &
476	rda_src_lats ! latitudes of the searched points
477
478	INTEGER (KIND=int_kind), DIMENSION(4), INTENT(out) :: &
479	ida_nbr_found ! indicates for each line how many points found
480
481	INTEGER (KIND=int_kind), INTENT(out) :: &
482	bin ! actual bin for target point
483	!
484	! !PARAMETERS:
485	!
486	REAL (KIND=dbl_kind), INTENT(in) :: &
487	rd_plat, & ! latitude of the target point
488	rd_plon ! longitude of the target point
489
490	LOGICAL, INTENT(in) :: ld_extrapdone ! true if extrapolation done
491
492	INTEGER (KIND=int_kind) :: &
493	ib_k, ib_j, ib_i, & ! iteration indices
494	il_min, il_max, il_inter ! begin and end for actual bin
495
496	INTEGER (KIND=int_kind), DIMENSION(4,2) :: &
497	ila_corners ! temp addresses for bins
498
499	!
500	! !DESCRIPTION:
501	! This routine finds the location of the target point in the source
502	! grid and returns those of the 16 nearest points that are unmasked.
503	! The source grid is a reduced grid.
504	!
505	! !REVISION HISTORY:
506	! 2002.10.21 J.Ghattas created
507	!
508	! EOP
509	!-----------------------------------------------------------------------
510	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
511	! $Author: valcke $
512	!-----------------------------------------------------------------------
513
514
515	!
516	! serch of actual bin
517	!
518
519
520	IF (rd_plat > bin_lats_r(1,1)) THEN ! norther of the first bin
521	bin=0
522	ila_corners(1:2,1:2)= 0
523	ila_corners(3,1)= bin_addr1_r(1,1)+1
524	ila_corners(3,2)= bin_addr1_r(2,1)
525	ila_corners(4,1)= bin_addr1_r(1,2)
526	ila_corners(4,2)= bin_addr1_r(2,2)
527
528	ELSE IF (rd_plat > bin_lats_r(1,2)) THEN ! in the first bin
529	bin=1
530	ila_corners(1,1:2)= 0
531	ila_corners(2,1)= bin_addr1_r(1,1)+1
532	ila_corners(2,2)= bin_addr1_r(2,1)
533	ila_corners(3,1)= bin_addr1_r(1,2)
534	ila_corners(3,2)= bin_addr1_r(2,2)
535	ila_corners(4,1)= bin_addr1_r(1,3)
536	ila_corners(4,2)= bin_addr1_r(2,3)
537
538	ELSE IF (rd_plat < bin_lats_r(1,num_srch_red)) THEN
539	! South of the last complet bin
540	bin=num_srch_red
541	ila_corners(1,1) = bin_addr1_r(1,num_srch_red-1)
542	ila_corners(1,2) = bin_addr1_r(2,num_srch_red-1)
543	ila_corners(2,1) = bin_addr1_r(1,num_srch_red)
544	ila_corners(2,2) = bin_addr1_r(2,num_srch_red)
545	ila_corners(3:4,1:2) = 0
546
547	ELSE IF (rd_plat < bin_lats_r(1,num_srch_red-1)) THEN
548	! in the last bin which is complet
549	! the bin (num_srch_red-1) is the last bin which is complet
550	bin=num_srch_red-1
551	ila_corners(1,1) = bin_addr1_r(1,num_srch_red-2)
552	ila_corners(1,2) = bin_addr1_r(2,num_srch_red-2)
553	ila_corners(2,1) = bin_addr1_r(1,num_srch_red-1)
554	ila_corners(2,2) = bin_addr1_r(2,num_srch_red-1)
555	ila_corners(3,1) = bin_addr1_r(1,num_srch_red)
556	ila_corners(3,2) = bin_addr1_r(2,num_srch_red)
557	ila_corners(4,1:2) = 0
558	ELSE
559	il_min=2
560	il_max=num_srch_red-1
561	DO WHILE (il_min /= il_max-1)
562	il_inter=(il_max-il_min)/2 + il_min
563	IF (rd_plat <= bin_lats_r(1,il_min) .and. &
564	rd_plat > bin_lats_r(1,il_inter)) THEN
565	il_max=il_inter
566	ELSE
567	il_min=il_inter
568	END IF
569	END DO
570	bin=il_min
571	ila_corners(1,1) = bin_addr1_r(1,bin-1)
572	ila_corners(1,2) = bin_addr1_r(2,bin-1)
573	ila_corners(2,1) = bin_addr1_r(1,bin)
574	ila_corners(2,2) = bin_addr1_r(2,bin)
575	ila_corners(3,1) = bin_addr1_r(1,bin+1)
576	ila_corners(3,2) = bin_addr1_r(2,bin+1)
577	ila_corners(4,1) = bin_addr1_r(1,bin+2)
578	ila_corners(4,2) = bin_addr1_r(2,bin+2)
579
580	IF (ila_corners(1,1)==0) THEN
581	ila_corners(1,1)=1
582	END IF
583	END IF
584
585	DO ib_k=1,4
586	IF (ila_corners(ib_k,1) .NE. 0) &
587	rda_src_lats(ib_k)= grid1_center_lat(ila_corners(ib_k,1))
588	ENDDO
589
590	!
591	! now perform a more detailed search for each line
592	!
593
594	ida_src_add(:,:)=0
595	ida_nbr_found(:)=0
596	rda_src_lons(:,:)=0
597
598	DO ib_k=1,4 ! for each line of found points
599	IF (ila_corners(ib_k,1)==0) THEN
600	CYCLE
601	END IF
602
603	il_min=ila_corners(ib_k,1)
604	il_max=ila_corners(ib_k,2)
605
606	IF (rd_plon < grid1_center_lon(il_min)) THEN
607	DO ib_j=il_max-1, il_max
608	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
609	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
610	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
611	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
612	grid1_center_lon(ib_j)-pi2
613	END IF
614	END DO
615	DO ib_j=il_min, il_min+1
616	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
617	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
618	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
619	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
620	grid1_center_lon(ib_j)
621	END IF
622	END DO
623
624	ELSE IF (rd_plon < grid1_center_lon(il_min+1)) THEN
625	IF (grid1_mask(il_max) .or. ld_extrapdone) THEN
626	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
627	ida_src_add(ib_k,ida_nbr_found(ib_k)) = il_max
628	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
629	grid1_center_lon(il_max)-pi2
630	END IF
631	DO ib_j=il_min, il_min+2
632	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
633	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
634	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
635	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
636	grid1_center_lon(ib_j)
637	END IF
638	END DO
639
640	ELSE IF (rd_plon > grid1_center_lon(il_max)) THEN
641	DO ib_j=il_max-1, il_max
642	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
643	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
644	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
645	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
646	grid1_center_lon(ib_j)
647	END IF
648	END DO
649	DO ib_j=il_min, il_min+1
650	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
651	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
652	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
653	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
654	grid1_center_lon(ib_j)+pi2
655	END IF
656	END DO
657
658	ELSE IF (rd_plon > grid1_center_lon(il_max-1)) THEN
659	DO ib_j=il_max-2, il_max
660	IF (grid1_mask(ib_j) .or. ld_extrapdone) THEN
661	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
662	ida_src_add(ib_k,ida_nbr_found(ib_k)) = ib_j
663	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
664	grid1_center_lon(ib_j)
665	END IF
666	END DO
667	IF (grid1_mask(il_min) .or. ld_extrapdone) THEN
668	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
669	ida_src_add(ib_k,ida_nbr_found(ib_k)) = il_min
670	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
671	grid1_center_lon(il_min)+pi2
672	END IF
673
674	ELSE
675
676	DO WHILE (il_min/=il_max-1)
677	il_inter=(il_max-il_min)/2 + il_min
678	IF (rd_plon >= grid1_center_lon(il_min) .and. &
679	rd_plon < grid1_center_lon(il_inter)) THEN
680	il_max=il_inter
681	ELSE
682	il_min=il_inter
683	END IF
684	END DO
685	DO ib_i= il_min-1, il_min+2
686	IF (grid1_mask(ib_i) .or. ld_extrapdone) THEN
687	ida_nbr_found(ib_k)=ida_nbr_found(ib_k)+1
688	ida_src_add(ib_k,ida_nbr_found(ib_k))=ib_i
689	rda_src_lons(ib_k,ida_nbr_found(ib_k))= &
690	grid1_center_lon(ib_i)
691	END IF
692	END DO
693
694	END IF
695
696	END DO
697
698
699	END SUBROUTINE grid_search_16_points
700
701
702
703	!-----------------------------------------------------------------------
704	! BOP
705	!
706	! !IROUTINE: calcul_wght_irreg
707	!
708	! !INTERFACE:
709	!
710	SUBROUTINE calcul_wght_irreg(rda_x, rd_pt, rda_wght)
711	!
712	! !USES:
713	!
714	! !RETURN VALUE:
715	!
716	REAL (KIND=dbl_kind), DIMENSION(4), INTENT(out) :: &
717	rda_wght ! array of weights for the points x
718	!
719	! !PARAMETERS:
720	!
721	REAL (KIND=dbl_kind), DIMENSION(4), INTENT(in) :: &
722	rda_x ! array of positions on source grid, lat or lon
723
724	REAL (KIND=dbl_kind),INTENT(in) :: &
725	rd_pt ! position of target point to interpolate
726
727	REAL (KIND=dbl_kind) :: &
728	rl_t1, rl_t2, rl_t3, rl_t4, rl_t5, rl_t6, rl_t7, rl_t8, rl_t9, &
729	rl_u1, rl_u2, rl_u3, rl_u4, &
730	rl_k1, rl_k2, rl_k3, &
731	rl_d1, rl_d2, rl_d3, rl_d4, &
732	rl_c1, rl_c2, rl_c3, rl_c4, &
733	rl_b1, rl_b2, rl_b3, rl_b4, &
734	rl_a1, rl_a2, rl_a3, rl_a4, &
735	rl_y1, rl_y2, rl_y3, &
736	rl_a1_y, rl_a2_y, rl_a3_y, rl_a4_y, &
737	rl_b1_y, rl_b2_y, rl_b3_y, rl_b4_y, &
738	rl_c1_y, rl_c2_y, rl_c3_y, rl_c4_y
739	!
740	! !DESCRIPTION:
741	! Calculates a the weights of four points for a bicubic interpolation.
742	! The distances between the points can be irregulier.
743	!
744	! !REVISION HISTORY:
745	! 2002.10.21 J.Ghattas created
746	!
747	! EOP
748	!-----------------------------------------------------------------------
749	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
750	! $Author: valcke $
751	!-----------------------------------------------------------------------
752
753
754	IF (rda_x(1)/=0.and. rda_x(2)/=0 .and. rda_x(3)/=0 .and. rda_x(4)/=0) THEN
755
756	rl_t1 = 1/rda_x(1) - 1/rda_x(2)
757	rl_t2 = 1/rda_x(1)2 - 1/rda_x(2)2
758	rl_t3 = 1/rda_x(1)3 - 1/rda_x(2)3
759	rl_t4 = 1/rda_x(1) - 1/rda_x(3)
760	rl_t5 = 1/rda_x(1)2 - 1/rda_x(3)2
761	rl_t6 = 1/rda_x(1)3 - 1/rda_x(3)3
762	rl_t7 = 1/rda_x(1) - 1/rda_x(4)
763	rl_t8 = 1/rda_x(1)2 - 1/rda_x(4)2
764	rl_t9 = 1/rda_x(1)3 - 1/rda_x(4)3
765
766	rl_u1 = rl_t2/rl_t1 - rl_t5/rl_t4
767	rl_u2 = rl_t3/rl_t1 - rl_t6/rl_t4
768	rl_u3 = rl_t2/rl_t1 - rl_t8/rl_t7
769	rl_u4 = rl_t3/rl_t1 - rl_t9/rl_t7
770
771	rl_k1 = (1/(rl_t1rl_u1)-1/(rl_t1rl_u3)) / (rl_u2/rl_u1-rl_u4/rl_u3)
772	rl_k2 = -1/(rl_t4*rl_u1) / (rl_u2/rl_u1-rl_u4/rl_u3)
773	rl_k3 = 1/(rl_t7*rl_u3) / (rl_u2/rl_u1-rl_u4/rl_u3)
774
775
776	rl_d1=(rl_k1+rl_k2+rl_k3)/rda_x(1)**3
777	rl_d2 = -rl_k1/rda_x(2)**3
778	rl_d3 = -rl_k2/rda_x(3)**3
779	rl_d4 = -rl_k3/rda_x(4)**3
780
781	rl_c1 = 1/rl_u1(1/(rl_t1rda_x(1)*3)-1/(rl_t4rda_x(1)**3)- &
782	rl_u2*rl_d1)
783	rl_c2 = 1/rl_u1(1/(-rl_t1rda_x(2)*3)-rl_u2rl_d2)
784	rl_c3 = 1/rl_u1(1/(rl_t4rda_x(3)*3)-rl_u2rl_d3)
785	rl_c4 = 1/rl_u1(-rl_u2rl_d4)
786
787	rl_b1 = 1/rl_t1/rda_x(1)*3-rl_t2/rl_t1rl_c1-rl_t3/rl_t1*rl_d1
788	rl_b2 = -1/rl_t1/rda_x(2)*3-rl_t2/rl_t1rl_c2-rl_t3/rl_t1*rl_d2
789	rl_b3 = -rl_t2/rl_t1rl_c3-rl_t3/rl_t1rl_d3
790	rl_b4 = -rl_t2/rl_t1rl_c4-rl_t3/rl_t1rl_d4
791
792	rl_a1 = 1/rda_x(1)*3-1/rda_x(1)rl_b1-1/rda_x(1)*2rl_c1- &
793	1/rda_x(1)*3rl_d1
794	rl_a2 = -1/rda_x(1)rl_b2-1/rda_x(1)2rl_c2-1/rda_x(1)*3rl_d2
795	rl_a3 = -1/rda_x(1)rl_b3-1/rda_x(1)2rl_c3-1/rda_x(1)*3rl_d3
796	rl_a4 = -1/rda_x(1)rl_b4-1/rda_x(1)2rl_c4-1/rda_x(1)*3rl_d4
797
798	! Weights
799	rda_wght(1) = rl_a1rd_pt3 + rl_b1rd_pt*2 + rl_c1rd_pt + rl_d1
800	rda_wght(2) = rl_a2rd_pt3 + rl_b2rd_pt*2 + rl_c2rd_pt + rl_d2
801	rda_wght(3) = rl_a3rd_pt3 + rl_b3rd_pt*2 + rl_c3rd_pt + rl_d3
802	rda_wght(4) = rl_a4rd_pt3 + rl_b4rd_pt*2 + rl_c4rd_pt + rl_d4
803
804	ELSE ! there is one point = 0
805
806	rl_d1=0; rl_d2=0; rl_d3=0; rl_d4=0
807
808	! Transformation for each case
809	IF (rda_x(1)==0) THEN
810	rl_y1=rda_x(2); rl_y2=rda_x(3); rl_y3=rda_x(4)
811	rl_d1=1
812	ELSE IF (rda_x(2)==0) THEN
813	rl_y1=rda_x(1); rl_y2=rda_x(3); rl_y3=rda_x(4)
814	rl_d2=1
815	ELSE IF (rda_x(3)==0) THEN
816	rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(4)
817	rl_d3=1
818	ELSE
819	rl_y1=rda_x(1); rl_y2=rda_x(2); rl_y3=rda_x(3)
820	rl_d4=1
821	END IF
822
823	! Solving the system
824	rl_t1 = 1/rl_y1-1/rl_y2
825	rl_t2 = 1/rl_y12-1/rl_y22
826	rl_t3 = 1/rl_y1-1/rl_y3
827	rl_t4 = 1/rl_y12-1/rl_y32
828
829	rl_c1_y =(1/rl_y13/rl_t1-1/rl_y13/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3)
830	rl_c2_y = -1/rl_y2**3/rl_t1/(rl_t2/rl_t1-rl_t4/rl_t3)
831	rl_c3_y = 1/rl_y3**3/rl_t3/(rl_t2/rl_t1-rl_t4/rl_t3)
832	rl_c4_y=(-1/rl_y13/rl_t1+1/rl_y23/rl_t1+ &
833	1/rl_y13/rl_t3-1/rl_y33/rl_t3)/(rl_t2/rl_t1-rl_t4/rl_t3)
834
835	rl_b1_y = 1/rl_y1*3/rl_t1 - rl_c1_yrl_t2/rl_t1
836	rl_b2_y = -1/rl_y2*3/rl_t1 - rl_c2_yrl_t2/rl_t1
837	rl_b3_y = -rl_c3_y*rl_t2/rl_t1
838	rl_b4_y = -1/rl_y13/rl_t1 + 1/rl_y23/rl_t1 - rl_c4_y*rl_t2/rl_t1
839
840	rl_a1_y = 1/rl_y13 - rl_b1_y/rl_y1 - rl_c1_y/rl_y12
841	rl_a2_y = -rl_b2_y/rl_y1 - rl_c2_y/rl_y1**2
842	rl_a3_y = -rl_b3_y/rl_y1 - rl_c3_y/rl_y1**2
843	rl_a4_y = -1/rl_y13 - rl_b4_y/rl_y1 - rl_c4_y/rl_y12
844
845	! Retransformation
846	IF (rda_x(1)==0) THEN
847	rl_a1=rl_a4_y; rl_a2=rl_a1_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y
848	rl_b1=rl_b4_y; rl_b2=rl_b1_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y
849	rl_c1=rl_c4_y; rl_c2=rl_c1_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y
850	ELSE IF (rda_x(2)==0) THEN
851	rl_a1=rl_a1_y; rl_a2=rl_a4_y; rl_a3=rl_a2_y; rl_a4=rl_a3_y
852	rl_b1=rl_b1_y; rl_b2=rl_b4_y; rl_b3=rl_b2_y; rl_b4=rl_b3_y
853	rl_c1=rl_c1_y; rl_c2=rl_c4_y; rl_c3=rl_c2_y; rl_c4=rl_c3_y
854	ELSE IF (rda_x(3)==0) THEN
855	rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a4_y; rl_a4=rl_a3_y
856	rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b4_y; rl_b4=rl_b3_y
857	rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c4_y; rl_c4=rl_c3_y
858	ELSE
859	rl_a1=rl_a1_y; rl_a2=rl_a2_y; rl_a3=rl_a3_y; rl_a4=rl_a4_y
860	rl_b1=rl_b1_y; rl_b2=rl_b2_y; rl_b3=rl_b3_y; rl_b4=rl_b4_y
861	rl_c1=rl_c1_y; rl_c2=rl_c2_y; rl_c3=rl_c3_y; rl_c4=rl_c4_y
862	END IF
863
864	! Weights
865	rda_wght(1) = rl_a1rd_pt3 + rl_b1rd_pt*2 + rl_c1rd_pt +rl_d1
866	rda_wght(2) = rl_a2rd_pt3 + rl_b2rd_pt*2 + rl_c2rd_pt +rl_d2
867	rda_wght(3) = rl_a3rd_pt3 + rl_b3rd_pt*2 + rl_c3rd_pt +rl_d3
868	rda_wght(4) = rl_a4rd_pt3 + rl_b4rd_pt*2 + rl_c4rd_pt +rl_d4
869
870	END IF
871
872
873	END SUBROUTINE calcul_wght_irreg
874
875	!-----------------------------------------------------------------------
876	! BOP
877	!
878	! !IROUTINE: calcul_wght_3
879	!
880	! !INTERFACE:
881
882	SUBROUTINE calcul_wght_3(rda_x, rd_pt, rda_wght)
883
884	! !USES:
885
886	! !RETURN VALUE:
887
888	REAL (KIND=dbl_kind), DIMENSION(3), INTENT(out) :: &
889	rda_wght ! array of weights for the points x
890
891	! !PARAMETERS:
892
893	REAL (KIND=dbl_kind), DIMENSION(3), INTENT(in) :: &
894	rda_x ! array of positions on source grid, lat or lon
895
896	REAL (KIND=dbl_kind), INTENT(in) :: &
897	rd_pt ! position of target point to interpolate
898
899	REAL (KIND=dbl_kind) :: &
900	rl_c1, rl_c2, rl_c3, &
901	rl_a1, rl_a2, rl_a3, &
902	rl_b1, rl_b2, rl_b3, &
903	rl_t1, rl_t2, rl_t3, rl_t4
904
905	! !DESCRIPTION:
906	! Calculates a the weights of 3 points for a parabolic interpolation.
907	!
908	! !REVISION HISTORY:
909	! 2002.10.21 J.Ghattas created
910	!
911	! EOP
912	!-----------------------------------------------------------------------
913	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
914	! $Author: valcke $
915	!-----------------------------------------------------------------------
916
917
918	IF (rda_x(1)/=0 .and. rda_x(2)/=0 .and. rda_x(3)/=0) THEN
919	rl_t1 = 1/rda_x(1)-1/rda_x(2)
920	rl_t2 = 1/rda_x(1)2-1/rda_x(2)2
921	rl_t3 = 1/rda_x(1)-1/rda_x(3)
922	rl_t4 = 1/rda_x(1)2-1/rda_x(3)2
923
924	rl_c1 = (1/rda_x(1)2/rl_t1-1/rda_x(1)2/rl_t3) / &
925	(rl_t2/rl_t1-rl_t4/rl_t3)
926	rl_c2 = -1/rda_x(2)**2/rl_t1 / (rl_t2/rl_t1-rl_t4/rl_t3)
927	rl_c3 = 1/rda_x(3)**2/rl_t3 / (rl_t2/rl_t1-rl_t4/rl_t3)
928
929	rl_b1 = 1/rda_x(1)*2/rl_t1 - rl_c1rl_t2/rl_t1
930	rl_b2 = -1/rda_x(2)*2/rl_t1 - rl_c2rl_t2/rl_t1
931	rl_b3 = - rl_c3*rl_t2/rl_t1
932
933	rl_a1 = 1/rda_x(1)2 - rl_b1/rda_x(1) - rl_c1/rda_x(1)2
934	rl_a2 = - rl_b2/rda_x(1) - rl_c2/rda_x(1)**2
935	rl_a3 = - rl_b3/rda_x(1) - rl_c3/rda_x(1)**2
936
937
938	ELSE IF (rda_x(1)==0) THEN
939	rl_c1 = 1; rl_c2 = 0; rl_c3 = 0
940	rl_b1 = (-1/rda_x(2)2+1/rda_x(3)2) / (1/rda_x(2)-1/rda_x(3))
941	rl_b2 = 1/rda_x(2)**2 / (1/rda_x(2)-1/rda_x(3))
942	rl_b3 = -1/rda_x(3)**2 / (1/rda_x(2)-1/rda_x(3))
943
944	rl_a1 = -1/rda_x(2)**2 - rl_b1/rda_x(2)
945	rl_a2 = 1/rda_x(2)**2 - rl_b2/rda_x(2)
946	rl_a3 = - rl_b3/rda_x(2)
947
948	ELSE IF (rda_x(2)==0) THEN
949
950	rl_c1 = 0; rl_c2 = 1; rl_c3 = 0
951	rl_b1 = 1/rda_x(1)**2 / (1/rda_x(1)-1/rda_x(3))
952	rl_b2 = (-1/rda_x(1)2+1/rda_x(3)2) / (1/rda_x(1)-1/rda_x(3))
953	rl_b3 = -1/rda_x(3)**2 / (1/rda_x(1)-1/rda_x(3))
954
955	rl_a1 = 1/rda_x(1)**2 - rl_b1/rda_x(1)
956	rl_a2 = -1/rda_x(1)**2 - rl_b2/rda_x(1)
957	rl_a3 = - rl_b3/rda_x(1)
958
959	ELSE !rda_x(3)==0
960	rl_c1 = 0; rl_c2 = 0; rl_c3 = 1
961	rl_b1 = 1/rda_x(1)**2 / (1/rda_x(1)-1/rda_x(2))
962	rl_b2 = -1/rda_x(2)**2 / (1/rda_x(1)-1/rda_x(2))
963	rl_b3 = (-1/rda_x(1)2+1/rda_x(2)2) / (1/rda_x(1)-1/rda_x(2))
964
965	rl_a1 = 1/rda_x(1)**2 - rl_b1/rda_x(1)
966	rl_a2 = - rl_b2/rda_x(1)
967	rl_a3 = -1/rda_x(1)**2 - rl_b3/rda_x(1)
968
969
970	END IF
971
972	! Weights
973	rda_wght(1) = rl_a1rd_pt2 + rl_b1rd_pt + rl_c1
974	rda_wght(2) = rl_a2rd_pt2 + rl_b2rd_pt + rl_c2
975	rda_wght(3) = rl_a3rd_pt2 + rl_b3rd_pt + rl_c3
976
977
978	END SUBROUTINE calcul_wght_3
979
980
981	!-----------------------------------------------------------------------
982	! BOP
983	!
984	! !IROUTINE: calcul_wght_2
985	!
986	! !INTERFACE:
987
988	SUBROUTINE calcul_wght_2(rda_x, rd_pt, rda_wght)
989
990	! !USES:
991
992	! !RETURN VALUE:
993
994	REAL (KIND=dbl_kind), DIMENSION(2), INTENT(out) :: &
995	rda_wght ! array of weights for the points x
996
997	! !PARAMETERS:
998
999	REAL (KIND=dbl_kind), DIMENSION(2), INTENT(in) :: &
1000	rda_x ! array of positions on source grid, lat or lon
1001
1002	REAL (KIND=dbl_kind), INTENT(in) :: &
1003	rd_pt ! position of target point to interpolate
1004
1005	REAL (KIND=dbl_kind) :: rl_b1, rl_b2, rl_a1, rl_a2
1006
1007	! !DESCRIPTION:
1008	! Calculates a the weights of 2 points for a linair interpolation.
1009	!
1010	! !REVISION HISTORY:
1011	! 2002.10.21 J.Ghattas created
1012	!
1013	! EOP
1014	!-----------------------------------------------------------------------
1015	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
1016	! $Author: valcke $
1017	!-----------------------------------------------------------------------
1018
1019
1020	IF (rda_x(1)/=0 .and. rda_x(2)/=0) THEN
1021	rl_b1 = 1/(1-rda_x(1)/rda_x(2))
1022	rl_b2 = -1/(rda_x(2)/rda_x(1)-1)
1023	rl_a1 = 1/rda_x(1) - rl_b1/rda_x(1)
1024	rl_a2 = - rl_b2/rda_x(1)
1025
1026	ELSE IF (rda_x(1)==0) THEN
1027	rl_b1=1; rl_b2=0
1028	rl_a1=-1/rda_x(2)
1029	rl_a2=1/rda_x(2)
1030	ELSE
1031	rl_b1=0; rl_b2=1
1032	rl_a1=1/rda_x(1)
1033	rl_a2=-1/rda_x(1)
1034	END IF
1035
1036	rda_wght(1) = rl_a1*rd_pt + rl_b1
1037	rda_wght(2) = rl_a2*rd_pt + rl_b2
1038
1039	END SUBROUTINE calcul_wght_2
1040
1041
1042	!-----------------------------------------------------------------------
1043	! BOP
1044	!
1045	! !IROUTINE: store_link_bicub
1046	!
1047	! !INTERFACE:
1048
1049	SUBROUTINE store_link_bicub(id_dst_add, ida_src_add, rda_wght)
1050
1051	! !USES:
1052
1053	! !RETURN VALUE:
1054
1055	! !PARAMETERS:
1056
1057	INTEGER (KIND=int_kind), INTENT(in) :: &
1058	id_dst_add ! address on destination grid
1059
1060	INTEGER (KIND=int_kind), DIMENSION(4,4), INTENT(in) :: &
1061	ida_src_add ! addresses for links on source grid
1062
1063	REAL (KIND=dbl_kind), DIMENSION(4,4), INTENT(in) :: &
1064	rda_wght ! array of remapping weights for these links
1065
1066	INTEGER (KIND=int_kind) :: ib_i, &
1067	il_num_links_old ! placeholder for old link number
1068
1069	INTEGER (KIND=int_kind), DIMENSION(16) :: &
1070	ila_src_add ! reshaped addresses
1071
1072	REAL (KIND=dbl_kind), DIMENSION(16) :: &
1073	rla_wght ! reshaped weights
1074
1075	! !DESCRIPTION:
1076	! This routine stores the addresses and weights for 16 links associated
1077	! with one destination point in the appropriate address.
1078	!
1079	! !REVISION HISTORY:
1080	! 2002.10.21 J.Ghattas created
1081	!
1082	! EOP
1083	!-----------------------------------------------------------------------
1084	! $Id: remap_bicubic_reduced.F90 2826 2010-12-10 11:14:21Z valcke $
1085	! $Author: valcke $
1086	!-----------------------------------------------------------------------
1087
1088
1089	!
1090	! Increment number of links and check if remap arrays need
1091	! to be increased to accomodate the new link. then store the link.
1092	!
1093
1094	il_num_links_old = num_links_map1
1095	num_links_map1 = il_num_links_old + 16
1096
1097	IF (num_links_map1 > max_links_map1) THEN
1098	CALL resize_remap_vars(1,MAX(resize_increment,16))
1099	END IF
1100
1101	ila_src_add=RESHAPE(ida_src_add,(/16/))
1102	rla_wght=RESHAPE(rda_wght,(/16/))
1103
1104	DO ib_i=1,16
1105	grid1_add_map1(il_num_links_old+ib_i) = ila_src_add(ib_i)
1106	grid2_add_map1(il_num_links_old+ib_i) = id_dst_add
1107	wts_map1(1,il_num_links_old+ib_i) = rla_wght(ib_i)
1108	END DO
1109
1110	END SUBROUTINE store_link_bicub
1111
1112
1113	END MODULE remap_bicubic_reduced
1114
1115	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1116

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

Context Navigation

source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/lib/scrip/src/remap_bicubic_reduced.F90 @ 4775

Download in other formats: