source: utils/tools/NESTING/src/bicubic_interp.f90 @ 12253

Last change on this file since 12253 was 2143, checked in by rblod, 10 years ago

Improvement of FCM branch

  • Property svn:keywords set to Id
File size: 34.5 KB
Line 
1!
2MODULE bicubic_interp
3  !
4  USE agrif_modutil
5  !
6  !************************************************************************
7  !                           *
8  ! MODULE  BICUBIC INTERP                *
9  !                           *
10  ! bicubic interpolation routines from SCRIP package       *     
11  !                           *
12  !http://climate.lanl.gov/Software/SCRIP/            *
13  !                           *
14  !Bicubic remapping                   *
15  !                           *
16  !************************************************************************
17  !                 
18  IMPLICIT NONE
19  !
20  INTEGER :: grid1_size,grid2_size,grid1_rank, grid2_rank     
21  INTEGER, DIMENSION(:), POINTER :: grid1_dims, grid2_dims 
22  LOGICAL, DIMENSION(:), POINTER :: grid1_mask,grid2_mask       
23  REAL*8,DIMENSION(:),POINTER :: &
24       grid1_center_lat,  &
25       grid1_center_lon,  & 
26       grid2_center_lat,  &
27       grid2_center_lon,  &
28       grid1_frac,        & 
29       grid2_frac           
30  REAL*8,DIMENSION(:,:), POINTER :: grid1_bound_box,grid2_bound_box   
31  INTEGER, PARAMETER :: num_srch_bins = 90 
32  INTEGER,DIMENSION(:,:),POINTER :: bin_addr1,bin_addr2 
33  REAL*8, DIMENSION(:,:),POINTER :: bin_lats,bin_lons       
34  REAL*8, PARAMETER :: zero   = 0.0,  &
35       one    = 1.0,  &
36       two    = 2.0,  &
37       three  = 3.0,  &
38       four   = 4.0,  &
39       five   = 5.0,  & 
40       half   = 0.5,  &
41       quart  = 0.25, &
42       bignum = 1.e+20, &
43       tiny   = 1.e-14, &
44       pi     = 3.14159265359, &
45       pi2    = two*pi, &
46       pih    = half*pi       
47
48  REAL*8, PARAMETER :: deg2rad = pi/180.   
49  INTEGER , PARAMETER :: max_iter = 100   
50  REAL*8, PARAMETER :: converge = 1.e-10 
51  INTEGER, PARAMETER :: norm_opt_none    = 1 &
52       ,norm_opt_dstarea = 2 &
53       ,norm_opt_frcarea = 3
54  !
55  INTEGER, PARAMETER :: map_type_conserv  = 1 &
56       ,map_type_bilinear = 2 &
57       ,map_type_bicubic  = 3 &
58       ,map_type_distwgt  = 4
59  !
60  INTEGER :: max_links_map1  &  ! current size of link arrays
61       ,num_links_map1  &  ! actual number of links for remapping
62       ,max_links_map2  &  ! current size of link arrays
63       ,num_links_map2  &  ! actual number of links for remapping
64       ,num_maps        &  ! num of remappings for this grid pair
65       ,num_wts         &  ! num of weights used in remapping
66       ,map_type        &  ! identifier for remapping method
67       ,norm_opt        &  ! option for normalization (conserv only)
68       ,resize_increment ! default amount to increase array size
69
70  INTEGER , DIMENSION(:), POINTER :: &
71       bicubic_grid1_add_map1, &  ! grid1 address for each link in mapping 1
72       bicubic_grid2_add_map1, &  ! grid2 address for each link in mapping 1
73       grid1_add_map2, &  ! grid1 address for each link in mapping 2
74       grid2_add_map2   ! grid2 address for each link in mapping 2
75
76  REAL*8, DIMENSION(:,:), POINTER ::   &
77       bicubic_wts_map1, &   ! map weights for each link (num_wts,max_links)
78       wts_map2 => NULL()    ! map weights for each link (num_wts,max_links)
79
80
81
82
83
84
85CONTAINS
86
87
88
89
90
91
92
93  !***********************************************************************
94  SUBROUTINE get_remap_bicub(grid1_lat,grid2_lat,grid1_lon,grid2_lon,mask, &
95       remap_matrix,source_add,destination_add)
96
97    !-----------------------------------------------------------------------
98    !
99    !     this routine computes the weights for a bicubic interpolation.
100    !
101    !-----------------------------------------------------------------------
102    REAL*8,DIMENSION(:,:),POINTER :: grid1_lat,grid2_lat,grid1_lon,grid2_lon
103    LOGICAL,DIMENSION(:,:) :: mask
104    !     
105    INTEGER,DIMENSION(:),POINTER :: source_add,destination_add
106    REAL*8,DIMENSION(:,:),POINTER :: remap_matrix 
107    !-----------------------------------------------------------------------
108    !
109    !     local variables
110    !
111    !-----------------------------------------------------------------------
112
113    INTEGER :: n,icount,dst_add,iter,nx,ny,i,j,ip1,e_add,jp1,n_add,ne_add,nele
114    REAL*8, DIMENSION(4) :: tmp_lats, tmp_lons
115    INTEGER, DIMENSION(4) :: src_add
116
117    REAL*8, DIMENSION(4)  :: src_lats,src_lons
118    INTEGER lastsrc_add
119    REAL*8, DIMENSION(4,4)  :: wgts
120    REAL*8 :: dlat,dlon
121    REAL*8 :: plat, plon,iguess, jguess,thguess, phguess,deli, delj, &
122         dth1, dth2, dth3,dph1, dph2, dph3,dthp, dphp,          &
123         mat1, mat2, mat3, mat4,determinant,sum_wgts,           &
124         w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,w11,w12,w13,w14,w15,w16
125
126
127    !-----------------------------------------------------------------------
128    !
129    !     compute mappings from grid1 to grid2
130    !
131    !-----------------------------------------------------------------------
132    !
133    ALLOCATE(grid1_dims(2),grid2_dims(2))
134    grid1_dims(1) = SIZE(grid1_lat,2)
135    grid1_dims(2) = SIZE(grid1_lat,1)
136    grid2_dims(1) = SIZE(grid2_lat,2)
137    grid2_dims(2) = SIZE(grid2_lat,1)
138    grid1_size = SIZE(grid1_lat,2) * SIZE(grid1_lat,1)
139    grid2_size = SIZE(grid2_lat,2) * SIZE(grid2_lat,1) 
140    !     
141    ALLOCATE( grid2_mask(grid2_size),         &
142         grid1_bound_box (4,grid1_size), &
143         grid2_bound_box (4,grid2_size), &
144         grid1_frac      (grid1_size),   &
145         grid2_frac      (grid2_size))
146
147    grid1_frac = zero
148    grid2_frac = zero
149
150    !
151    ! 2D array -> 1D array
152    !
153    ALLOCATE(grid1_center_lat(SIZE(grid1_lat,1)*SIZE(grid1_lat,2)))
154    CALL tab2Dto1D(grid1_lat,grid1_center_lat)
155
156    ALLOCATE(grid1_center_lon(SIZE(grid1_lon,1)*SIZE(grid1_lon,2)))
157    CALL tab2Dto1D(grid1_lon,grid1_center_lon)
158
159    ALLOCATE(grid2_center_lat(SIZE(grid2_lat,1)*SIZE(grid2_lat,2)))
160    CALL tab2Dto1D(grid2_lat,grid2_center_lat)
161
162    ALLOCATE(grid2_center_lon(SIZE(grid2_lon,1)*SIZE(grid2_lon,2)))     
163    CALL tab2Dto1D(grid2_lon,grid2_center_lon) 
164
165    ALLOCATE(grid1_mask(SIZE(grid1_lat,1)*SIZE(grid1_lat,2)))
166    CALL logtab2Dto1D(mask,grid1_mask)
167    !     
168    !      Write(*,*) ,'grid1_mask = ',grid1_mask
169    !             
170    !
171    ! degrees to radian
172    !
173    grid1_center_lat = grid1_center_lat*deg2rad
174    grid1_center_lon = grid1_center_lon*deg2rad
175    grid2_center_lat = grid2_center_lat*deg2rad
176    grid2_center_lon = grid2_center_lon*deg2rad
177
178    !-----------------------------------------------------------------------
179    !     convert longitudes to 0,2pi interval
180    !-----------------------------------------------------------------------
181
182    WHERE (grid1_center_lon .GT. pi2)  grid1_center_lon =       &
183         grid1_center_lon - pi2
184    WHERE (grid1_center_lon .LT. zero) grid1_center_lon =       &
185         grid1_center_lon + pi2
186    WHERE (grid2_center_lon .GT. pi2)  grid2_center_lon =       &
187         grid2_center_lon - pi2
188    WHERE (grid2_center_lon .LT. zero) grid2_center_lon =       &
189         grid2_center_lon + pi2
190
191    !-----------------------------------------------------------------------
192    !
193    !     make sure input latitude range is within the machine values
194    !     for +/- pi/2
195    !
196    !-----------------------------------------------------------------------
197
198    WHERE (grid1_center_lat >  pih) grid1_center_lat =  pih
199    WHERE (grid1_center_lat < -pih) grid1_center_lat = -pih
200    WHERE (grid2_center_lat >  pih) grid2_center_lat =  pih
201    WHERE (grid2_center_lat < -pih) grid2_center_lat = -pih
202    !
203
204    !
205    !-----------------------------------------------------------------------
206    !
207    !     compute bounding boxes for restricting future grid searches
208    !
209    !-----------------------------------------------------------------------
210    !
211    nx = grid1_dims(1)
212    ny = grid1_dims(2)
213
214    DO n=1,grid1_size
215
216       !*** find N,S and NE points to this grid point
217
218       j = (n - 1)/nx +1
219       i = n - (j-1)*nx
220
221       IF (i < nx) THEN
222          ip1 = i + 1
223       ELSE
224          !*** assume cyclic
225          ip1 = 1
226          !*** but if it is not, correct
227          e_add = (j - 1)*nx + ip1
228          IF (ABS(grid1_center_lat(e_add) -     &
229               grid1_center_lat(n   )) > pih) THEN
230             ip1 = i
231          ENDIF
232          ip1=nx
233       ENDIF
234
235       IF (j < ny) THEN
236          jp1 = j+1
237       ELSE
238          !*** assume cyclic
239          jp1 = 1
240          !*** but if it is not, correct
241          n_add = (jp1 - 1)*nx + i
242          IF (ABS(grid1_center_lat(n_add) -             &
243               grid1_center_lat(n   )) > pih) THEN
244             jp1 = j
245          ENDIF
246          jp1=ny
247       ENDIF
248
249       n_add = (jp1 - 1)*nx + i
250       e_add = (j - 1)*nx + ip1
251       ne_add = (jp1 - 1)*nx + ip1
252
253       !*** find N,S and NE lat/lon coords and check bounding box
254
255       tmp_lats(1) = grid1_center_lat(n)
256       tmp_lats(2) = grid1_center_lat(e_add)
257       tmp_lats(3) = grid1_center_lat(ne_add)
258       tmp_lats(4) = grid1_center_lat(n_add)
259
260       tmp_lons(1) = grid1_center_lon(n)
261       tmp_lons(2) = grid1_center_lon(e_add)
262       tmp_lons(3) = grid1_center_lon(ne_add)
263       tmp_lons(4) = grid1_center_lon(n_add)
264
265       grid1_bound_box(1,n) = MINVAL(tmp_lats)
266       grid1_bound_box(2,n) = MAXVAL(tmp_lats)
267
268       grid1_bound_box(3,n) = MINVAL(tmp_lons)
269       grid1_bound_box(4,n) = MAXVAL(tmp_lons)
270    END DO
271
272    nx = grid2_dims(1)
273    ny = grid2_dims(2)
274
275    DO n=1,grid2_size
276
277       !*** find N,S and NE points to this grid point
278
279       j = (n - 1)/nx +1
280       i = n - (j-1)*nx
281
282       IF (i < nx) THEN
283          ip1 = i + 1
284       ELSE
285          !*** assume cyclic
286          ip1 = 1
287          !*** but if it is not, correct
288          e_add = (j - 1)*nx + ip1
289          IF (ABS(grid2_center_lat(e_add) -  &
290               grid2_center_lat(n   )) > pih) THEN
291             ip1 = i
292          ENDIF
293       ENDIF
294
295       IF (j < ny) THEN
296          jp1 = j+1
297       ELSE
298          !*** assume cyclic
299          jp1 = 1
300          !*** but if it is not, correct
301          n_add = (jp1 - 1)*nx + i
302          IF (ABS(grid2_center_lat(n_add) -  &
303               grid2_center_lat(n   )) > pih) THEN
304             jp1 = j
305          ENDIF
306       ENDIF
307
308       n_add = (jp1 - 1)*nx + i
309       e_add = (j - 1)*nx + ip1
310       ne_add = (jp1 - 1)*nx + ip1
311
312       !*** find N,S and NE lat/lon coords and check bounding box
313
314       tmp_lats(1) = grid2_center_lat(n)
315       tmp_lats(2) = grid2_center_lat(e_add)
316       tmp_lats(3) = grid2_center_lat(ne_add)
317       tmp_lats(4) = grid2_center_lat(n_add)
318
319       tmp_lons(1) = grid2_center_lon(n)
320       tmp_lons(2) = grid2_center_lon(e_add)
321       tmp_lons(3) = grid2_center_lon(ne_add)
322       tmp_lons(4) = grid2_center_lon(n_add)
323
324       grid2_bound_box(1,n) = MINVAL(tmp_lats)
325       grid2_bound_box(2,n) = MAXVAL(tmp_lats)
326       grid2_bound_box(3,n) = MINVAL(tmp_lons)
327       grid2_bound_box(4,n) = MAXVAL(tmp_lons)
328    END DO
329    !
330    !
331    !
332    WHERE (ABS(grid1_bound_box(4,:) - grid1_bound_box(3,:)) > pi)
333       grid1_bound_box(3,:) = zero
334       grid1_bound_box(4,:) = pi2
335    END WHERE
336
337    WHERE (ABS(grid2_bound_box(4,:) - grid2_bound_box(3,:)) > pi)
338       grid2_bound_box(3,:) = zero
339       grid2_bound_box(4,:) = pi2
340    END WHERE
341
342    !***
343    !*** try to check for cells that overlap poles
344    !***
345
346    WHERE (grid1_center_lat > grid1_bound_box(2,:)) &
347         grid1_bound_box(2,:) = pih
348
349    WHERE (grid1_center_lat < grid1_bound_box(1,:)) &
350         grid1_bound_box(1,:) = -pih
351
352    WHERE (grid2_center_lat > grid2_bound_box(2,:)) &
353         grid2_bound_box(2,:) = pih
354
355    WHERE (grid2_center_lat < grid2_bound_box(1,:)) &
356         grid2_bound_box(1,:) = -pih
357
358    !-----------------------------------------------------------------------
359    !     set up and assign address ranges to search bins in order to
360    !     further restrict later searches
361    !-----------------------------------------------------------------------
362
363    ALLOCATE(bin_addr1(2,num_srch_bins))
364    ALLOCATE(bin_addr2(2,num_srch_bins))
365    ALLOCATE(bin_lats (2,num_srch_bins))
366    ALLOCATE(bin_lons (2,num_srch_bins))
367
368    dlat = pi/num_srch_bins
369
370    DO n=1,num_srch_bins
371       bin_lats(1,n) = (n-1)*dlat - pih
372       bin_lats(2,n) =     n*dlat - pih
373       bin_lons(1,n) = zero
374       bin_lons(2,n) = pi2
375       bin_addr1(1,n) = grid1_size + 1
376       bin_addr1(2,n) = 0
377       bin_addr2(1,n) = grid2_size + 1
378       bin_addr2(2,n) = 0
379    END DO
380
381    DO nele=1,grid1_size
382       DO n=1,num_srch_bins
383          IF (grid1_bound_box(1,nele) <= bin_lats(2,n) .AND.   &
384               grid1_bound_box(2,nele) >= bin_lats(1,n)) THEN
385             bin_addr1(1,n) = MIN(nele,bin_addr1(1,n))
386             bin_addr1(2,n) = MAX(nele,bin_addr1(2,n))
387          ENDIF
388       END DO
389    END DO
390
391    DO nele=1,grid2_size
392       DO n=1,num_srch_bins
393          IF (grid2_bound_box(1,nele) <= bin_lats(2,n) .AND.    &
394               grid2_bound_box(2,nele) >= bin_lats(1,n)) THEN
395             bin_addr2(1,n) = MIN(nele,bin_addr2(1,n))
396             bin_addr2(2,n) = MAX(nele,bin_addr2(2,n))
397          ENDIF
398       END DO
399    END DO
400
401    CALL init_remap_vars   
402    !
403    !***
404    !*** loop over destination grid
405    !***
406    !     
407    grid2_mask = .TRUE.
408    !
409    lastsrc_add=1
410    !         
411    WRITE(*,*) 'Bicubic remapping weights computation'
412    !     
413    DO dst_add = 1, grid2_size
414
415       plat = grid2_center_lat(dst_add)
416       plon = grid2_center_lon(dst_add)         
417       !           
418       !           
419       CALL grid_search_bicub(src_add, src_lats, src_lons,     &
420            plat, plon, grid1_dims,              &
421            grid1_center_lat, grid1_center_lon,  &
422            grid1_bound_box, bin_addr1, bin_addr2,lastsrc_add)
423       !       
424       !           
425       IF (src_add(1) > 0) THEN
426          DO n=1,4
427             IF (.NOT. grid1_mask(src_add(n))) src_add(1) = 0
428          END DO
429       ENDIF
430
431       !-----------------------------------------------------------------------
432       !
433       !       if point found, find local i,j coordinates for weights
434       !
435       !-----------------------------------------------------------------------
436
437       IF (src_add(1) > 0) THEN
438
439          grid2_frac(dst_add) = one
440          !
441          dth1 = src_lats(2) - src_lats(1)
442          dth2 = src_lats(4) - src_lats(1)
443          dth3 = src_lats(3) - src_lats(2) - dth2
444
445          dph1 = src_lons(2) - src_lons(1)
446          dph2 = src_lons(4) - src_lons(1)
447          dph3 = src_lons(3) - src_lons(2)
448
449          IF (dph1 >  three*pih) dph1 = dph1 - pi2
450          IF (dph2 >  three*pih) dph2 = dph2 - pi2
451          IF (dph3 >  three*pih) dph3 = dph3 - pi2
452          IF (dph1 < -three*pih) dph1 = dph1 + pi2
453          IF (dph2 < -three*pih) dph2 = dph2 + pi2
454          IF (dph3 < -three*pih) dph3 = dph3 + pi2
455
456          dph3 = dph3 - dph2
457          iguess = half
458          jguess = half
459
460          DO iter=1,max_iter
461
462             dthp = plat-src_lats(1)-dth1*iguess- &
463                  dth2*jguess-dth3*iguess*jguess
464             dphp = plon - src_lons(1)
465
466             IF (dphp >  three*pih) dphp = dphp - pi2
467             IF (dphp < -three*pih) dphp = dphp + pi2
468
469             dphp = dphp - dph1*iguess - dph2*jguess - & 
470                  dph3*iguess*jguess
471             !
472             mat1 = dth1 + dth3*jguess
473             mat2 = dth2 + dth3*iguess
474             mat3 = dph1 + dph3*jguess
475             mat4 = dph2 + dph3*iguess
476             !
477             determinant = mat1*mat4 - mat2*mat3
478             !
479             deli = (dthp*mat4 - mat2*dphp)/determinant
480             delj = (mat1*dphp - dthp*mat3)/determinant
481             !
482             IF (ABS(deli) < converge .AND.            & 
483                  ABS(delj) < converge) EXIT
484             !
485             iguess = iguess + deli
486             jguess = jguess + delj
487
488          END DO
489
490          IF (iter <= max_iter) THEN
491             !
492             wgts(1,1) = (one - jguess**2*(three-two*jguess))*   &
493                  (one - iguess**2*(three-two*iguess))
494
495             wgts(1,2) = (one - jguess**2*(three-two*jguess))*   &
496                  iguess**2*(three-two*iguess)
497
498             wgts(1,3) = jguess**2*(three-two*jguess)*           &
499                  iguess**2*(three-two*iguess)
500
501             wgts(1,4) = jguess**2*(three-two*jguess)*           &
502                  (one - iguess**2*(three-two*iguess))
503
504             wgts(2,1) = (one - jguess**2*(three-two*jguess))*   &
505                  iguess*(iguess-one)**2
506
507             wgts(2,2) = (one - jguess**2*(three-two*jguess))*   &
508                  iguess**2*(iguess-one)
509
510             wgts(2,3) = jguess**2*(three-two*jguess)*           &
511                  iguess**2*(iguess-one)
512
513             wgts(2,4) =  jguess**2*(three-two*jguess)*          &
514                  iguess*(iguess-one)**2
515
516             wgts(3,1) =  jguess*(jguess-one)**2*                &
517                  (one - iguess**2*(three-two*iguess))
518
519             wgts(3,2) =  jguess*(jguess-one)**2*                &
520                  iguess**2*(three-two*iguess)                           
521
522             wgts(3,3) =  jguess**2*(jguess-one)*                &
523                  iguess**2*(three-two*iguess)
524
525             wgts(3,4) =  jguess**2*(jguess-one)*                &
526                  (one - iguess**2*(three-two*iguess))
527
528             wgts(4,1) =  iguess*(iguess-one)**2*                &
529                  jguess*(jguess-one)**2
530
531             wgts(4,2) =  iguess**2*(iguess-one)*                &
532                  jguess*(jguess-one)**2
533
534             wgts(4,3) =  iguess**2*(iguess-one)*                &
535                  jguess**2*(jguess-one)
536
537             wgts(4,4) =  iguess*(iguess-one)**2*                &
538                  jguess**2*(jguess-one)
539
540             CALL store_link_bicub(dst_add, src_add, wgts)
541
542          ELSE
543             STOP 'Iteration for i,j exceed max iteration count'
544          ENDIF
545          !
546       ELSE IF (src_add(1) < 0) THEN
547
548          src_add = ABS(src_add)
549          !
550          icount = 0
551          !         
552          DO n=1,4
553             IF (grid1_mask(src_add(n))) THEN
554                icount = icount + 1
555             ELSE
556                src_lats(n) = zero
557             ENDIF
558          END DO
559
560          IF (icount > 0) THEN
561
562             sum_wgts = SUM(src_lats)
563             wgts(1,1) = src_lats(1)/sum_wgts
564             wgts(1,2) = src_lats(2)/sum_wgts
565             wgts(1,3) = src_lats(3)/sum_wgts
566             wgts(1,4) = src_lats(4)/sum_wgts
567             wgts(2:4,:) = zero
568
569             CALL store_link_bicub(dst_add, src_add, wgts)
570
571          ENDIF
572
573       ENDIF
574       !
575    END DO
576    !             
577    ALLOCATE(remap_matrix(SIZE(bicubic_wts_map1,1),SIZE(bicubic_wts_map1,2)), &
578         source_add(num_links_map1),       &
579         destination_add(num_links_map1))
580
581    DO j = 1,SIZE(bicubic_wts_map1,2)
582       DO i = 1,SIZE(bicubic_wts_map1,1)
583
584          remap_matrix(i,j) = bicubic_wts_map1(i,j)
585
586       END DO
587    END DO
588    !       
589    source_add(:) = bicubic_grid1_add_map1(1:num_links_map1)
590    destination_add(:) = bicubic_grid2_add_map1(1:num_links_map1)
591    !               
592    !
593    WHERE(destination_add == 0)
594       destination_add = 1
595    END WHERE
596
597    WHERE(source_add == 0)
598       source_add = 1
599    END WHERE
600    !               
601    DEALLOCATE(grid1_bound_box,grid2_bound_box,grid1_center_lat,grid1_center_lon)
602    DEALLOCATE(grid2_center_lat,grid2_center_lon,bicubic_grid1_add_map1,bicubic_wts_map1)
603    DEALLOCATE(grid1_frac,grid2_frac,grid1_dims,grid2_dims,grid2_mask)
604    DEALLOCATE(bin_addr1,bin_addr2,bin_lats,bin_lons)
605    DEALLOCATE(grid1_mask,bicubic_grid2_add_map1) 
606    !
607    !-----------------------------------------------------------------------
608
609  END SUBROUTINE get_remap_bicub
610
611  !
612  !***********************************************************************
613
614  !***********************************************************************
615!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
616  !************************************************************************
617  !   SUBROUTINE GRID_SEARCH_BILIN
618  !************************************************************************
619!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
620  !
621  SUBROUTINE grid_search_bicub(src_add, src_lats, src_lons,   &
622       plat, plon, src_grid_dims,      &
623       src_center_lat, src_center_lon, & 
624       src_grid_bound_box,             &
625       src_bin_add, dst_bin_add,lastsrc_add)
626
627    !-----------------------------------------------------------------------
628    !
629    !     this routine finds the location of the search point plat, plon
630    !     in the source grid and returns the corners needed for a bilinear
631    !     interpolation.
632    !
633    !-----------------------------------------------------------------------
634
635    !-----------------------------------------------------------------------
636    !     output variables
637    !-----------------------------------------------------------------------
638    !
639    ! address of each corner point enclosing P
640    !
641    INTEGER,DIMENSION(4) :: src_add 
642    REAL*8,DIMENSION(4) :: src_lats,src_lons 
643    !     
644    !-----------------------------------------------------------------------
645    !     input variables
646    !-----------------------------------------------------------------------
647    !
648    ! latitude, longitude of the search point
649    !
650    REAL*8, INTENT(in) :: plat,plon   
651    !
652    ! size of each src grid dimension
653    !
654    INTEGER, DIMENSION(2), INTENT(in) :: src_grid_dims 
655    !
656    ! latitude, longitude of each src grid center
657    !
658    REAL*8, DIMENSION(:), INTENT(in) :: src_center_lat,src_center_lon 
659    !
660    ! bound box for source grid
661    !
662    REAL*8, DIMENSION(:,:), INTENT(in) :: src_grid_bound_box 
663    !
664    ! latitude bins for restricting searches
665    !
666    INTEGER, DIMENSION(:,:), INTENT(in) ::src_bin_add,dst_bin_add 
667
668    INTEGER,OPTIONAL :: lastsrc_add
669    INTEGER :: loopsrc,l1,l2
670    !     
671    !-----------------------------------------------------------------------
672    !     local variables
673    !-----------------------------------------------------------------------
674    !
675    INTEGER :: n,next_n,srch_add,nx, ny,min_add, max_add,  &
676         i, j, jp1, ip1, n_add, e_add, ne_add
677
678
679    REAL*8 ::  vec1_lat, vec1_lon,vec2_lat, vec2_lon, cross_product,  &
680         cross_product_last,coslat_dst, sinlat_dst, coslon_dst, &
681         sinlon_dst,dist_min, distance 
682
683    !-----------------------------------------------------------------------
684    !     restrict search first using bins
685    !-----------------------------------------------------------------------
686
687    src_add = 0
688    min_add = SIZE(src_center_lat)
689    max_add = 1
690    DO n=1,num_srch_bins
691       IF (plat >= bin_lats(1,n) .AND. plat <= bin_lats(2,n) .AND. &
692            plon >= bin_lons(1,n) .AND. plon <= bin_lons(2,n)) THEN
693          min_add = MIN(min_add, src_bin_add(1,n))
694          max_add = MAX(max_add, src_bin_add(2,n))
695       ENDIF
696    END DO
697
698    !-----------------------------------------------------------------------
699    !     now perform a more detailed search
700    !-----------------------------------------------------------------------
701
702    nx = src_grid_dims(1)
703    ny = src_grid_dims(2)
704
705    loopsrc=0
706    DO WHILE (loopsrc <= max_add)
707
708
709       l1=MAX(min_add,lastsrc_add-loopsrc)
710       l2=MIN(max_add,lastsrc_add+loopsrc)     
711
712       loopsrc = loopsrc+1
713
714       srch_loop: DO srch_add = l1,l2,MAX(l2-l1,1)
715
716          !*** first check bounding box
717
718          IF (plat <= src_grid_bound_box(2,srch_add) .AND. & 
719               plat >= src_grid_bound_box(1,srch_add) .AND.  &
720               plon <= src_grid_bound_box(4,srch_add) .AND.  &
721               plon >= src_grid_bound_box(3,srch_add)) THEN
722
723
724             !***
725             !*** we are within bounding box so get really serious
726             !***
727             !*** determine neighbor addresses
728             !
729             j = (srch_add - 1)/nx +1
730             i = srch_add - (j-1)*nx
731             !
732             IF (i < nx) THEN
733                ip1 = i + 1
734             ELSE
735                ip1 = 1
736             ENDIF
737             !
738             IF (j < ny) THEN
739                jp1 = j+1
740             ELSE
741                jp1 = 1
742             ENDIF
743             !
744             n_add = (jp1 - 1)*nx + i
745             e_add = (j - 1)*nx + ip1
746             ne_add = (jp1 - 1)*nx + ip1
747             !
748             src_lats(1) = src_center_lat(srch_add)
749             src_lats(2) = src_center_lat(e_add)
750             src_lats(3) = src_center_lat(ne_add)
751             src_lats(4) = src_center_lat(n_add)
752             !
753             src_lons(1) = src_center_lon(srch_add)
754             src_lons(2) = src_center_lon(e_add)
755             src_lons(3) = src_center_lon(ne_add)
756             src_lons(4) = src_center_lon(n_add)
757             !
758             !***
759             !*** for consistency, we must make sure all lons are in
760             !*** same 2pi interval
761             !***
762             !
763             vec1_lon = src_lons(1) - plon
764             IF (vec1_lon >  pi) THEN
765                src_lons(1) = src_lons(1) - pi2
766             ELSE IF (vec1_lon < -pi) THEN
767                src_lons(1) = src_lons(1) + pi2
768             ENDIF
769             DO n=2,4
770                vec1_lon = src_lons(n) - src_lons(1)
771                IF (vec1_lon >  pi) THEN
772                   src_lons(n) = src_lons(n) - pi2
773                ELSE IF (vec1_lon < -pi) THEN
774                   src_lons(n) = src_lons(n) + pi2
775                ENDIF
776             END DO
777             !
778             corner_loop: DO n=1,4
779                next_n = MOD(n,4) + 1
780                !***
781                !*** here we take the cross product of the vector making
782                !*** up each box side with the vector formed by the vertex
783                !*** and search point.  if all the cross products are
784                !*** positive, the point is contained in the box.
785                !***
786                vec1_lat = src_lats(next_n) - src_lats(n)
787                vec1_lon = src_lons(next_n) - src_lons(n)
788                vec2_lat = plat - src_lats(n)
789                vec2_lon = plon - src_lons(n)
790                !***
791                !*** check for 0,2pi crossings
792                !***
793                IF (vec1_lon >  three*pih) THEN
794                   vec1_lon = vec1_lon - pi2
795                ELSE IF (vec1_lon < -three*pih) THEN
796                   vec1_lon = vec1_lon + pi2
797                ENDIF
798                IF (vec2_lon >  three*pih) THEN
799                   vec2_lon = vec2_lon - pi2
800                ELSE IF (vec2_lon < -three*pih) THEN
801                   vec2_lon = vec2_lon + pi2
802                ENDIF
803                !
804                cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
805                !
806                !***
807                !*** if cross product is less than zero, this cell
808                !*** doesn't work
809                !***
810                IF (n == 1) cross_product_last = cross_product
811                IF (cross_product*cross_product_last < zero) &
812                     EXIT corner_loop
813                cross_product_last = cross_product
814                !
815             END DO corner_loop
816             !***
817             !*** if cross products all same sign, we found the location
818             !***
819             IF (n > 4) THEN
820                src_add(1) = srch_add
821                src_add(2) = e_add
822                src_add(3) = ne_add
823                src_add(4) = n_add
824                !
825                lastsrc_add = srch_add
826                RETURN
827             ENDIF
828             !***
829             !*** otherwise move on to next cell
830             !***
831          ENDIF !bounding box check
832       END DO srch_loop
833
834
835    ENDDO
836
837    !***
838    !*** if no cell found, point is likely either in a box that
839    !*** straddles either pole or is outside the grid.  fall back
840    !*** to a distance-weighted average of the four closest
841    !*** points.  go ahead and compute weights here, but store
842    !*** in src_lats and return -add to prevent the parent
843    !*** routine from computing bilinear weights
844    !***
845    !print *,'Could not find location for ',plat,plon
846    !print *,'Using nearest-neighbor average for this point'
847    !
848    coslat_dst = COS(plat)
849    sinlat_dst = SIN(plat)
850    coslon_dst = COS(plon)
851    sinlon_dst = SIN(plon)
852    !
853    dist_min = bignum
854    src_lats = bignum
855    DO srch_add = min_add,max_add
856       distance = ACOS(coslat_dst*COS(src_center_lat(srch_add))*   &
857            (coslon_dst*COS(src_center_lon(srch_add)) +   &
858            sinlon_dst*SIN(src_center_lon(srch_add)))+   &
859            sinlat_dst*SIN(src_center_lat(srch_add)))
860
861       IF (distance < dist_min) THEN
862          sort_loop: DO n=1,4
863             IF (distance < src_lats(n)) THEN
864                DO i=4,n+1,-1
865                   src_add (i) = src_add (i-1)
866                   src_lats(i) = src_lats(i-1)
867                END DO
868                src_add (n) = -srch_add
869                src_lats(n) = distance
870                dist_min = src_lats(4)
871                EXIT sort_loop
872             ENDIF
873          END DO sort_loop
874       ENDIF
875    END DO
876    !
877    src_lons = one/(src_lats + tiny)
878    distance = SUM(src_lons)
879    src_lats = src_lons/distance
880
881    !-----------------------------------------------------------------------
882
883  END SUBROUTINE grid_search_bicub
884
885
886
887
888
889
890
891
892
893
894
895
896
897  !-----------------------------------------------------------------------
898
899  SUBROUTINE store_link_bicub(dst_add, src_add, weights)
900
901    !-----------------------------------------------------------------------
902    !
903    !     this routine stores the address and weight for four links
904    !     associated with one destination point in the appropriate address
905    !     and weight arrays and resizes those arrays if necessary.
906    !
907    !-----------------------------------------------------------------------
908
909    !-----------------------------------------------------------------------
910    !
911    !     input variables
912    !
913    !-----------------------------------------------------------------------
914
915    INTEGER :: dst_add
916    INTEGER, DIMENSION(4) :: src_add
917    REAL*8, DIMENSION(4,4) :: weights 
918
919    !-----------------------------------------------------------------------
920    !
921    !     local variables
922    !
923    !-----------------------------------------------------------------------
924
925    INTEGER :: n,num_links_old
926
927    !-----------------------------------------------------------------------
928    !
929    !     increment number of links and check to see if remap arrays need
930    !     to be increased to accomodate the new link.  then store the
931    !     link.
932    !
933    !-----------------------------------------------------------------------
934
935
936    num_links_old  = num_links_map1
937    num_links_map1 = num_links_old + 4
938
939    IF (num_links_map1 > max_links_map1)    & 
940         CALL resize_remap_vars(resize_increment)
941
942    DO n=1,4
943       bicubic_grid1_add_map1(num_links_old+n) = src_add(n)
944       bicubic_grid2_add_map1(num_links_old+n) = dst_add
945       bicubic_wts_map1    (:,num_links_old+n) = weights(:,n)
946    END DO
947
948
949
950
951
952    !-----------------------------------------------------------------------
953
954  END SUBROUTINE store_link_bicub
955
956  !***********************************************************************
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971  !
972  !************************************************************************
973  !   SUBROUTINE INIT_REMAP_VAR
974  !************************************************************************
975  !
976  SUBROUTINE init_remap_vars
977    !
978    num_wts = 4
979    !
980    num_links_map1 = 0
981    max_links_map1 = 4*grid2_size
982    IF (num_maps > 1) THEN
983       num_links_map2 = 0
984       max_links_map1 = MAX(4*grid1_size,4*grid2_size)
985       max_links_map2 = max_links_map1
986    ENDIF
987
988    resize_increment = 0.1*MAX(grid1_size,grid2_size)
989    !
990    !-----------------------------------------------------------------------
991    !
992    !     allocate address and weight arrays for mapping 1
993    !   
994    !-----------------------------------------------------------------------
995    !
996    ALLOCATE (bicubic_grid1_add_map1(max_links_map1),  &
997         bicubic_grid2_add_map1(max_links_map1),  &
998         bicubic_wts_map1(num_wts, max_links_map1))
999
1000    !-----------------------------------------------------------------------
1001    !
1002    !     allocate address and weight arrays for mapping 2 if necessary
1003    !   
1004    !-----------------------------------------------------------------------
1005
1006    IF (num_maps > 1) THEN
1007       ALLOCATE (grid1_add_map2(max_links_map2),    &
1008            grid2_add_map2(max_links_map2),    &
1009            wts_map2(num_wts, max_links_map2))
1010    ENDIF
1011
1012    !-----------------------------------------------------------------------
1013
1014  END SUBROUTINE init_remap_vars
1015
1016  !***********************************************************************
1017  !************************************************************************
1018  !   SUBROUTINE RESIZE_REMAP_VAR
1019  !************************************************************************
1020
1021  SUBROUTINE resize_remap_vars(increment)
1022
1023    !-----------------------------------------------------------------------
1024    !     this routine resizes remapping arrays by increasing(decreasing)
1025    !     the max_links by increment
1026    !-----------------------------------------------------------------------
1027
1028    !-----------------------------------------------------------------------
1029    !     input variables
1030    !-----------------------------------------------------------------------
1031
1032    INTEGER :: increment
1033
1034    !-----------------------------------------------------------------------
1035    !     local variables
1036    !-----------------------------------------------------------------------
1037
1038    INTEGER ::    &
1039         ierr,    &  ! error flag
1040         mxlinks   ! size of link arrays
1041
1042    INTEGER, DIMENSION(:), POINTER ::    &
1043         add1_tmp,   & ! temp array for resizing address arrays
1044         add2_tmp  ! temp array for resizing address arrays
1045    !
1046    ! temp array for resizing weight arrays
1047    !
1048    REAL*8, DIMENSION(:,:), POINTER :: wts_tmp   
1049    !
1050    !-----------------------------------------------------------------------
1051    !***
1052    !*** allocate temporaries to hold original values
1053    !***
1054    mxlinks = SIZE(bicubic_grid1_add_map1)
1055    ALLOCATE (add1_tmp(mxlinks), add2_tmp(mxlinks), &
1056         wts_tmp(num_wts,mxlinks))
1057
1058    add1_tmp = bicubic_grid1_add_map1
1059    add2_tmp = bicubic_grid2_add_map1
1060    wts_tmp  = bicubic_wts_map1
1061
1062    !***
1063    !*** deallocate originals and increment max_links then
1064    !*** reallocate arrays at new size
1065    !***
1066
1067    DEALLOCATE (bicubic_grid1_add_map1, bicubic_grid2_add_map1, bicubic_wts_map1)
1068    max_links_map1 = mxlinks + increment
1069    ALLOCATE (bicubic_grid1_add_map1(max_links_map1),    &
1070         bicubic_grid2_add_map1(max_links_map1),    &
1071         bicubic_wts_map1(num_wts,max_links_map1))
1072    !***
1073    !*** restore original values from temp arrays and
1074    !*** deallocate temps
1075    !***
1076    mxlinks = MIN(mxlinks, max_links_map1)
1077    bicubic_grid1_add_map1(1:mxlinks) = add1_tmp (1:mxlinks)
1078    bicubic_grid2_add_map1(1:mxlinks) = add2_tmp (1:mxlinks)
1079    bicubic_wts_map1    (:,1:mxlinks) = wts_tmp(:,1:mxlinks)
1080    DEALLOCATE(add1_tmp, add2_tmp, wts_tmp)
1081
1082    !-----------------------------------------------------------------------
1083    !
1084  END SUBROUTINE resize_remap_vars
1085  !
1086  !************************************************************************
1087  !
1088
1089END MODULE bicubic_interp
1090
1091!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.