New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
bilinear_interp.f90 in branches/UKMO/dev_r5107_mld_zint/NEMOGCM/TOOLS/NESTING/src – NEMO

source: branches/UKMO/dev_r5107_mld_zint/NEMOGCM/TOOLS/NESTING/src/bilinear_interp.f90 @ 5450

Last change on this file since 5450 was 5450, checked in by davestorkey, 9 years ago

Clear SVn keywords from UKMO/dev_r5107_mld_zint branch.

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