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.
agrif_gridsearch.f90 in branches/UKMO/dev_r5518_clean_shutdown/NEMOGCM/TOOLS/NESTING/src – NEMO

source: branches/UKMO/dev_r5518_clean_shutdown/NEMOGCM/TOOLS/NESTING/src/agrif_gridsearch.f90 @ 5674

Last change on this file since 5674 was 5674, checked in by dancopsey, 9 years ago

Stripped out SVN keywords.

File size: 25.3 KB
Line 
1!
2MODULE agrif_gridsearch
3  !
4  USE agrif_modutil
5  !
6  !************************************************************************
7  !                           *
8  ! MODULE  AGRIF_GRIDSEARCH                 *
9  !                           *
10  !************************************************************************
11  !     
12  !-----------------------------------------------------------------------
13  IMPLICIT NONE
14
15  !-----------------------------------------------------------------------
16  !     variables that describe each grid
17  !-----------------------------------------------------------------------
18  !
19  !      integer :: grid1_size,grid2_size,grid1_rank, grid2_rank     
20  !      integer, dimension(:), pointer :: grid1_dims, grid2_dims 
21  !      logical, dimension(:), pointer :: grid1_mask,grid2_mask       
22  !      real*8,dimension(:),pointer :: grid1_center_lat,grid1_center_lon,grid2_center_lat,grid2_center_lon
23  !
24  !      real*8,dimension(:,:), pointer :: grid1_bound_box,grid2_bound_box   
25  !      integer, parameter :: num_srch_bins = 90 
26  !      integer,dimension(:,:),pointer :: bin_addr1,bin_addr2
27  !      real*8, dimension(:,:),pointer :: bin_lats,bin_lons
28  REAL*8, PARAMETER :: zero   = 0.0,  &
29       one    = 1.0,  &
30       two    = 2.0,  &
31       three  = 3.0,  &
32       four   = 4.0,  &
33       five   = 5.0,  & 
34       half   = 0.5,  &
35       quart  = 0.25, &
36       bignum = 1.e+20, &
37       tiny   = 1.e-14, &
38       pi     = 3.14159265359, &
39       pi2    = two*pi, &
40       pih    = half*pi       
41  !     
42  REAL*8, PARAMETER :: deg2rad = pi/180.
43  !
44  !
45CONTAINS 
46  !     
47  !
48  SUBROUTINE get_detected_pts(grid1_lat,grid2_lat,grid1_lon,grid2_lon,maskC,maskF,detected_pts)
49    !
50    !-----------------------------------------------------------------------
51    !this routine makes any necessary changes (e.g. for 0,2pi longitude range)
52    !-----------------------------------------------------------------------
53    !
54    REAL*8,DIMENSION(:,:),POINTER :: grid1_lat,grid2_lat,grid1_lon,grid2_lon
55    LOGICAL, POINTER, DIMENSION(:,:) :: masksrc,maskdst 
56    LOGICAL, DIMENSION(:,:) :: detected_pts 
57    REAL*8,DIMENSION(:,:) :: maskF,maskC
58    LOGICAL,POINTER,DIMENSION(:) :: detected_pts_1D     
59    REAL*8 :: plat,plon
60    INTEGER :: lastsrc_add
61    INTEGER, DIMENSION(4) :: src_add 
62    REAL*8,DIMENSION(4) :: src_lats,src_lons     
63    INTEGER :: grid1_size,grid2_size,grid1_rank, grid2_rank     
64    INTEGER, DIMENSION(:), POINTER :: grid1_dims, grid2_dims 
65    LOGICAL, DIMENSION(:), POINTER :: grid1_mask,grid2_mask       
66    REAL*8,DIMENSION(:),POINTER :: grid1_center_lat,grid1_center_lon,grid2_center_lat,grid2_center_lon
67
68    REAL*8,DIMENSION(:,:), POINTER :: grid1_bound_box,grid2_bound_box   
69    !      integer, parameter :: num_srch_bins = 90 
70    INTEGER,DIMENSION(:,:),POINTER :: bin_addr1,bin_addr2 
71    REAL*8, DIMENSION(:,:),POINTER :: bin_lats,bin_lons 
72
73    REAL*8, PARAMETER :: zero   = 0.0,  &
74         one    = 1.0,  &
75         two    = 2.0,  &
76         three  = 3.0,  &
77         four   = 4.0,  &
78         five   = 5.0,  & 
79         half   = 0.5,  &
80         quart  = 0.25, &
81         bignum = 1.e+20, &
82         tiny   = 1.e-14, &
83         pi     = 3.14159265359, &
84         pi2    = two*pi, &
85         pih    = half*pi       
86
87    REAL*8, PARAMETER :: deg2rad = pi/180.
88    !
89    !-----------------------------------------------------------------------
90    ! local variables
91    !-----------------------------------------------------------------------
92    !
93    INTEGER :: n,nele,i,j,ip1,jp1,n_add,e_add,ne_add,nx,ny
94    INTEGER :: xpos,ypos,dst_add
95    !         
96    ! integer mask
97    !
98    INTEGER, DIMENSION(:), POINTER :: imask 
99    !
100    ! lat/lon intervals for search bins
101    !
102    REAL*8 :: dlat,dlon           
103    !     
104    ! temps for computing bounding boxes
105    !
106    REAL*8, DIMENSION(4) :: tmp_lats, tmp_lons 
107    !
108    !      write(*,*)'proceed to Bilinear interpolation ...'
109    !
110    ALLOCATE(grid1_dims(2),grid2_dims(2))
111    !     
112    grid1_dims(1) = SIZE(grid1_lat,2)
113    grid1_dims(2) = SIZE(grid1_lat,1)
114    grid2_dims(1) = SIZE(grid2_lat,2)
115    grid2_dims(2) = SIZE(grid2_lat,1)
116    grid1_size = SIZE(grid1_lat,2) * SIZE(grid1_lat,1)
117    grid2_size = SIZE(grid2_lat,2) * SIZE(grid2_lat,1) 
118    !     
119    !-----------------------------------------------------------------------
120    !     allocate grid coordinates/masks and read data
121    !-----------------------------------------------------------------------
122    !     
123    ALLOCATE( grid1_bound_box (4,grid1_size),grid2_bound_box (4,grid2_size))
124
125    ALLOCATE(detected_pts_1D(grid1_size))
126    ALLOCATE(masksrc(SIZE(maskC,1),SIZE(maskC,2)))
127    ALLOCATE(maskdst(SIZE(maskF,1),SIZE(maskF,2)))
128    !     
129    WHERE(maskC == 1.)
130       masksrc= .TRUE.
131    ELSEWHERE
132       masksrc= .FALSE.
133    END WHERE
134    !     
135    WHERE(maskF == 1.)
136       maskdst= .TRUE.
137    ELSEWHERE
138       maskdst= .FALSE.
139    END WHERE
140    !                 
141    !     
142    !
143    ! 2D array -> 1D array
144    !
145    ALLOCATE(grid1_center_lat(SIZE(grid1_lat,1)*SIZE(grid1_lat,2)))
146    CALL tab2Dto1D(grid1_lat,grid1_center_lat)
147
148    ALLOCATE(grid1_center_lon(SIZE(grid1_lon,1)*SIZE(grid1_lon,2)))
149    CALL tab2Dto1D(grid1_lon,grid1_center_lon)
150
151    ALLOCATE(grid2_center_lat(SIZE(grid2_lat,1)*SIZE(grid2_lat,2)))
152    CALL tab2Dto1D(grid2_lat,grid2_center_lat)
153
154    ALLOCATE(grid2_center_lon(SIZE(grid2_lon,1)*SIZE(grid2_lon,2)))     
155    CALL tab2Dto1D(grid2_lon,grid2_center_lon) 
156
157    ALLOCATE(grid1_mask(SIZE(grid1_lat,1)*SIZE(grid1_lat,2)))
158    CALL logtab2Dto1D(masksrc,grid1_mask)
159
160    ALLOCATE(grid2_mask(SIZE(grid2_lat,1)*SIZE(grid2_lat,2)))     
161    CALL logtab2Dto1D(maskdst,grid2_mask)
162    !                     
163    !
164    ! degrees to radian
165    !
166    grid1_center_lat = grid1_center_lat*deg2rad
167    grid1_center_lon = grid1_center_lon*deg2rad
168    grid2_center_lat = grid2_center_lat*deg2rad
169    grid2_center_lon = grid2_center_lon*deg2rad
170
171    !-----------------------------------------------------------------------
172    !     convert longitudes to 0,2pi interval
173    !-----------------------------------------------------------------------
174
175    WHERE (grid1_center_lon .GT. pi2)  grid1_center_lon =       &
176         grid1_center_lon - pi2
177    WHERE (grid1_center_lon .LT. zero) grid1_center_lon =       &
178         grid1_center_lon + pi2
179    WHERE (grid2_center_lon .GT. pi2)  grid2_center_lon =       &
180         grid2_center_lon - pi2
181    WHERE (grid2_center_lon .LT. zero) grid2_center_lon =       &
182         grid2_center_lon + pi2
183
184    !-----------------------------------------------------------------------
185    !
186    !     make sure input latitude range is within the machine values
187    !     for +/- pi/2
188    !
189    !-----------------------------------------------------------------------
190
191    WHERE (grid1_center_lat >  pih) grid1_center_lat =  pih
192    WHERE (grid1_center_lat < -pih) grid1_center_lat = -pih
193    WHERE (grid2_center_lat >  pih) grid2_center_lat =  pih
194    WHERE (grid2_center_lat < -pih) grid2_center_lat = -pih
195
196    !-----------------------------------------------------------------------
197    !
198    !     compute bounding boxes for restricting future grid searches
199    !
200    !-----------------------------------------------------------------------
201    !
202    nx = grid1_dims(1)
203    ny = grid1_dims(2)
204
205    DO n=1,grid1_size
206
207       !*** find N,S and NE points to this grid point
208
209       j = (n - 1)/nx +1
210       i = n - (j-1)*nx
211
212       IF (i < nx) THEN
213          ip1 = i + 1
214       ELSE
215          !*** assume cyclic
216          ip1 = 1
217          !*** but if it is not, correct
218          e_add = (j - 1)*nx + ip1
219          IF (ABS(grid1_center_lat(e_add) -     &
220               grid1_center_lat(n   )) > pih) THEN
221             ip1 = i
222          ENDIF
223          ip1=nx
224       ENDIF
225
226       IF (j < ny) THEN
227          jp1 = j+1
228       ELSE
229          !*** assume cyclic
230          jp1 = 1
231          !*** but if it is not, correct
232          n_add = (jp1 - 1)*nx + i
233          IF (ABS(grid1_center_lat(n_add) -             &
234               grid1_center_lat(n   )) > pih) THEN
235             jp1 = j
236          ENDIF
237          jp1=ny
238       ENDIF
239
240       n_add = (jp1 - 1)*nx + i
241       e_add = (j - 1)*nx + ip1
242       ne_add = (jp1 - 1)*nx + ip1
243
244       !*** find N,S and NE lat/lon coords and check bounding box
245
246       tmp_lats(1) = grid1_center_lat(n)
247       tmp_lats(2) = grid1_center_lat(e_add)
248       tmp_lats(3) = grid1_center_lat(ne_add)
249       tmp_lats(4) = grid1_center_lat(n_add)
250
251       tmp_lons(1) = grid1_center_lon(n)
252       tmp_lons(2) = grid1_center_lon(e_add)
253       tmp_lons(3) = grid1_center_lon(ne_add)
254       tmp_lons(4) = grid1_center_lon(n_add)
255
256       grid1_bound_box(1,n) = MINVAL(tmp_lats)
257       grid1_bound_box(2,n) = MAXVAL(tmp_lats)
258
259       grid1_bound_box(3,n) = MINVAL(tmp_lons)
260       grid1_bound_box(4,n) = MAXVAL(tmp_lons)
261    END DO
262
263    nx = grid2_dims(1)
264    ny = grid2_dims(2)
265
266    DO n=1,grid2_size
267
268       !*** find N,S and NE points to this grid point
269
270       j = (n - 1)/nx +1
271       i = n - (j-1)*nx
272
273       IF (i < nx) THEN
274          ip1 = i + 1
275       ELSE
276          !*** assume cyclic
277          ip1 = 1
278          !*** but if it is not, correct
279          e_add = (j - 1)*nx + ip1
280          IF (ABS(grid2_center_lat(e_add) -  &
281               grid2_center_lat(n   )) > pih) THEN
282             ip1 = i
283          ENDIF
284       ENDIF
285
286       IF (j < ny) THEN
287          jp1 = j+1
288       ELSE
289          !*** assume cyclic
290          jp1 = 1
291          !*** but if it is not, correct
292          n_add = (jp1 - 1)*nx + i
293          IF (ABS(grid2_center_lat(n_add) -  &
294               grid2_center_lat(n   )) > pih) THEN
295             jp1 = j
296          ENDIF
297       ENDIF
298
299       n_add = (jp1 - 1)*nx + i
300       e_add = (j - 1)*nx + ip1
301       ne_add = (jp1 - 1)*nx + ip1
302
303       !*** find N,S and NE lat/lon coords and check bounding box
304
305       tmp_lats(1) = grid2_center_lat(n)
306       tmp_lats(2) = grid2_center_lat(e_add)
307       tmp_lats(3) = grid2_center_lat(ne_add)
308       tmp_lats(4) = grid2_center_lat(n_add)
309
310       tmp_lons(1) = grid2_center_lon(n)
311       tmp_lons(2) = grid2_center_lon(e_add)
312       tmp_lons(3) = grid2_center_lon(ne_add)
313       tmp_lons(4) = grid2_center_lon(n_add)
314
315       grid2_bound_box(1,n) = MINVAL(tmp_lats)
316       grid2_bound_box(2,n) = MAXVAL(tmp_lats)
317       grid2_bound_box(3,n) = MINVAL(tmp_lons)
318       grid2_bound_box(4,n) = MAXVAL(tmp_lons)
319    END DO
320    !
321    !
322    !
323    WHERE (ABS(grid1_bound_box(4,:) - grid1_bound_box(3,:)) > pi)
324       grid1_bound_box(3,:) = zero
325       grid1_bound_box(4,:) = pi2
326    END WHERE
327
328    WHERE (ABS(grid2_bound_box(4,:) - grid2_bound_box(3,:)) > pi)
329       grid2_bound_box(3,:) = zero
330       grid2_bound_box(4,:) = pi2
331    END WHERE
332
333    !***
334    !*** try to check for cells that overlap poles
335    !***
336
337    WHERE (grid1_center_lat > grid1_bound_box(2,:)) &
338         grid1_bound_box(2,:) = pih
339
340    WHERE (grid1_center_lat < grid1_bound_box(1,:)) &
341         grid1_bound_box(1,:) = -pih
342
343    WHERE (grid2_center_lat > grid2_bound_box(2,:)) &
344         grid2_bound_box(2,:) = pih
345
346    WHERE (grid2_center_lat < grid2_bound_box(1,:)) &
347         grid2_bound_box(1,:) = -pih
348
349    !-----------------------------------------------------------------------
350    !     set up and assign address ranges to search bins in order to
351    !     further restrict later searches
352    !-----------------------------------------------------------------------
353
354    ALLOCATE(bin_addr1(2,90))
355    ALLOCATE(bin_addr2(2,90))
356    ALLOCATE(bin_lats (2,90))
357    ALLOCATE(bin_lons (2,90))
358
359    dlat = pi/90
360
361    DO n=1,90
362       bin_lats(1,n) = (n-1)*dlat - pih
363       bin_lats(2,n) =     n*dlat - pih
364       bin_lons(1,n) = zero
365       bin_lons(2,n) = pi2
366       bin_addr1(1,n) = grid1_size + 1
367       bin_addr1(2,n) = 0
368       bin_addr2(1,n) = grid2_size + 1
369       bin_addr2(2,n) = 0
370    END DO
371
372    DO nele=1,grid1_size
373       DO n=1,90
374          IF (grid1_bound_box(1,nele) <= bin_lats(2,n) .AND.   &
375               grid1_bound_box(2,nele) >= bin_lats(1,n)) THEN
376             bin_addr1(1,n) = MIN(nele,bin_addr1(1,n))
377             bin_addr1(2,n) = MAX(nele,bin_addr1(2,n))
378          ENDIF
379       END DO
380    END DO
381
382    DO nele=1,grid2_size
383       DO n=1,90
384          IF (grid2_bound_box(1,nele) <= bin_lats(2,n) .AND.    &
385               grid2_bound_box(2,nele) >= bin_lats(1,n)) THEN
386             bin_addr2(1,n) = MIN(nele,bin_addr2(1,n))
387             bin_addr2(2,n) = MAX(nele,bin_addr2(2,n))
388          ENDIF
389       END DO
390    END DO
391    !
392    ! Call init_remap_vars
393    !
394
395    lastsrc_add=1
396    detected_pts_1D = .FALSE.
397    !
398    DO dst_add = 1, grid2_size
399       !       
400       plat = grid2_center_lat(dst_add)
401       plon = grid2_center_lon(dst_add)
402       !***
403       !*** find nearest square of grid points on source grid
404       !***
405       CALL grid_search_bilin(bin_lons,bin_lats,src_add, src_lats, src_lons,          &
406            plat, plon, grid1_dims,               &
407            grid1_center_lat, grid1_center_lon,   & 
408            grid1_bound_box, bin_addr1, bin_addr2,lastsrc_add)
409       !
410       IF (src_add(1) > 0) THEN
411          !
412          IF(grid2_mask(dst_add)) THEN                  !mask true on destination grid     
413             DO n=1,4
414                IF(.NOT. grid1_mask(src_add(n))) THEN
415                   detected_pts_1D(src_add(n)) = .TRUE.
416                ENDIF
417             END DO
418          ENDIF
419       ENDIF
420    END DO
421    !
422    !
423    CALL logtab1Dto2D(detected_pts_1D,detected_pts,SIZE(detected_pts,2),SIZE(detected_pts,1))   
424    !
425    DEALLOCATE(detected_pts_1D,grid1_bound_box,grid2_bound_box)
426    DEALLOCATE(grid1_center_lon,grid1_center_lat,grid2_center_lon,grid2_center_lat)
427    DEALLOCATE(grid1_mask,grid2_mask,masksrc,maskdst)
428    DEALLOCATE(bin_addr1,bin_addr2,bin_lats,bin_lons)
429    DEALLOCATE(grid1_dims,grid2_dims)
430    !
431    !-----------------------------------------------------------------------
432
433  END SUBROUTINE get_detected_pts
434
435  !**********************************************************************
436!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
437  !************************************************************************
438  !   SUBROUTINE GRID_SEARCH_BILIN
439  !************************************************************************
440!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
441  !
442  SUBROUTINE grid_search_bilin(bin_lons,bin_lats,src_add, src_lats, src_lons,   &
443       plat, plon, src_grid_dims,      &
444       src_center_lat, src_center_lon, & 
445       src_grid_bound_box,             &
446       src_bin_add, dst_bin_add,lastsrc_add)
447
448    !-----------------------------------------------------------------------
449    !
450    !     this routine finds the location of the search point plat, plon
451    !     in the source grid and returns the corners needed for a bilinear
452    !     interpolation.
453    !
454    !-----------------------------------------------------------------------
455
456    !-----------------------------------------------------------------------
457    !     output variables
458    !-----------------------------------------------------------------------
459    !
460    ! address of each corner point enclosing P
461    !
462    INTEGER,DIMENSION(4) :: src_add 
463    REAL*8,DIMENSION(4) :: src_lats,src_lons 
464    !     
465    !-----------------------------------------------------------------------
466    !     input variables
467    !-----------------------------------------------------------------------
468    !
469    ! latitude, longitude of the search point
470    !
471    REAL*8, DIMENSION(:,:) :: bin_lats,bin_lons 
472    REAL*8, INTENT(in) :: plat,plon   
473    !
474    ! size of each src grid dimension
475    !
476    INTEGER, DIMENSION(2), INTENT(in) :: src_grid_dims 
477    !
478    ! latitude, longitude of each src grid center
479    !
480    REAL*8, DIMENSION(:), INTENT(in) :: src_center_lat,src_center_lon 
481    !
482    ! bound box for source grid
483    !
484    REAL*8, DIMENSION(:,:), INTENT(in) :: src_grid_bound_box 
485    !
486    ! latitude bins for restricting searches
487    !
488    INTEGER, DIMENSION(:,:), INTENT(in) ::src_bin_add,dst_bin_add 
489
490    INTEGER,OPTIONAL :: lastsrc_add
491    INTEGER :: loopsrc,l1,l2
492    !     
493    !-----------------------------------------------------------------------
494    !     local variables
495    !-----------------------------------------------------------------------
496    !
497    INTEGER :: n,next_n,srch_add,nx, ny,min_add, max_add,  &
498         i, j, jp1, ip1, n_add, e_add, ne_add
499
500
501    REAL*8 ::  vec1_lat, vec1_lon,vec2_lat, vec2_lon, cross_product,  &
502         cross_product_last,coslat_dst, sinlat_dst, coslon_dst, &
503         sinlon_dst,dist_min, distance 
504
505    !-----------------------------------------------------------------------
506    !     restrict search first using bins
507    !-----------------------------------------------------------------------
508
509    src_add = 0
510
511    min_add = SIZE(src_center_lat)
512    max_add = 1
513    DO n=1,90
514       IF (plat >= bin_lats(1,n) .AND. plat <= bin_lats(2,n) .AND. &
515            plon >= bin_lons(1,n) .AND. plon <= bin_lons(2,n)) THEN
516          min_add = MIN(min_add, src_bin_add(1,n))
517          max_add = MAX(max_add, src_bin_add(2,n))
518       ENDIF
519    END DO
520
521    !-----------------------------------------------------------------------
522    !     now perform a more detailed search
523    !-----------------------------------------------------------------------
524
525    nx = src_grid_dims(1)
526    ny = src_grid_dims(2)
527
528    loopsrc=0
529    DO WHILE (loopsrc <= max_add)
530
531
532       l1=MAX(min_add,lastsrc_add-loopsrc)
533       l2=MIN(max_add,lastsrc_add+loopsrc)     
534
535       loopsrc = loopsrc+1
536
537       srch_loop: DO srch_add = l1,l2,MAX(l2-l1,1)
538
539          !*** first check bounding box
540
541          IF (plat <= src_grid_bound_box(2,srch_add) .AND. & 
542               plat >= src_grid_bound_box(1,srch_add) .AND.  &
543               plon <= src_grid_bound_box(4,srch_add) .AND.  &
544               plon >= src_grid_bound_box(3,srch_add)) THEN
545             !***
546             !*** we are within bounding box so get really serious
547             !***
548             !*** determine neighbor addresses
549             !
550             j = (srch_add - 1)/nx +1
551             i = srch_add - (j-1)*nx
552             !
553             IF (i < nx) THEN
554                ip1 = i + 1
555             ELSE
556                ip1 = 1
557             ENDIF
558             !
559             IF (j < ny) THEN
560                jp1 = j+1
561             ELSE
562                jp1 = 1
563             ENDIF
564             !
565             n_add = (jp1 - 1)*nx + i
566             e_add = (j - 1)*nx + ip1
567             ne_add = (jp1 - 1)*nx + ip1
568             !
569             src_lats(1) = src_center_lat(srch_add)
570             src_lats(2) = src_center_lat(e_add)
571             src_lats(3) = src_center_lat(ne_add)
572             src_lats(4) = src_center_lat(n_add)
573             !
574             src_lons(1) = src_center_lon(srch_add)
575             src_lons(2) = src_center_lon(e_add)
576             src_lons(3) = src_center_lon(ne_add)
577             src_lons(4) = src_center_lon(n_add)
578             !
579             !***
580             !*** for consistency, we must make sure all lons are in
581             !*** same 2pi interval
582             !***
583             !
584             vec1_lon = src_lons(1) - plon
585             IF (vec1_lon >  pi) THEN
586                src_lons(1) = src_lons(1) - pi2
587             ELSE IF (vec1_lon < -pi) THEN
588                src_lons(1) = src_lons(1) + pi2
589             ENDIF
590             DO n=2,4
591                vec1_lon = src_lons(n) - src_lons(1)
592                IF (vec1_lon >  pi) THEN
593                   src_lons(n) = src_lons(n) - pi2
594                ELSE IF (vec1_lon < -pi) THEN
595                   src_lons(n) = src_lons(n) + pi2
596                ENDIF
597             END DO
598             !
599             corner_loop: DO n=1,4
600                next_n = MOD(n,4) + 1
601                !***
602                !*** here we take the cross product of the vector making
603                !*** up each box side with the vector formed by the vertex
604                !*** and search point.  if all the cross products are
605                !*** positive, the point is contained in the box.
606                !***
607                vec1_lat = src_lats(next_n) - src_lats(n)
608                vec1_lon = src_lons(next_n) - src_lons(n)
609                vec2_lat = plat - src_lats(n)
610                vec2_lon = plon - src_lons(n)
611                !***
612                !*** check for 0,2pi crossings
613                !***
614                IF (vec1_lon >  three*pih) THEN
615                   vec1_lon = vec1_lon - pi2
616                ELSE IF (vec1_lon < -three*pih) THEN
617                   vec1_lon = vec1_lon + pi2
618                ENDIF
619                IF (vec2_lon >  three*pih) THEN
620                   vec2_lon = vec2_lon - pi2
621                ELSE IF (vec2_lon < -three*pih) THEN
622                   vec2_lon = vec2_lon + pi2
623                ENDIF
624                !
625                cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
626                !
627                !***
628                !*** if cross product is less than zero, this cell
629                !*** doesn't work
630                !***
631                IF (n == 1) cross_product_last = cross_product
632                IF (cross_product*cross_product_last < zero) &
633                     EXIT corner_loop
634                cross_product_last = cross_product
635                !
636             END DO corner_loop
637             !***
638             !*** if cross products all same sign, we found the location
639             !***
640             IF (n > 4) THEN
641                src_add(1) = srch_add
642                src_add(2) = e_add
643                src_add(3) = ne_add
644                src_add(4) = n_add
645                !
646                lastsrc_add = srch_add
647                RETURN
648             ENDIF
649             !***
650             !*** otherwise move on to next cell
651             !***
652          ENDIF !bounding box check
653       END DO srch_loop
654
655
656    ENDDO
657
658
659    !***
660    !*** if no cell found, point is likely either in a box that
661    !*** straddles either pole or is outside the grid.  fall back
662    !*** to a distance-weighted average of the four closest
663    !*** points.  go ahead and compute weights here, but store
664    !*** in src_lats and return -add to prevent the parent
665    !*** routine from computing bilinear weights
666    !***
667    !print *,'Could not find location for ',plat,plon
668    !print *,'Using nearest-neighbor average for this point'
669    !
670    coslat_dst = COS(plat)
671    sinlat_dst = SIN(plat)
672    coslon_dst = COS(plon)
673    sinlon_dst = SIN(plon)
674    !
675    dist_min = bignum
676    src_lats = bignum
677    DO srch_add = min_add,max_add
678       distance = ACOS(coslat_dst*COS(src_center_lat(srch_add))*   &
679            (coslon_dst*COS(src_center_lon(srch_add)) +   &
680            sinlon_dst*SIN(src_center_lon(srch_add)))+   &
681            sinlat_dst*SIN(src_center_lat(srch_add)))
682
683       IF (distance < dist_min) THEN
684          sort_loop: DO n=1,4
685             IF (distance < src_lats(n)) THEN
686                DO i=4,n+1,-1
687                   src_add (i) = src_add (i-1)
688                   src_lats(i) = src_lats(i-1)
689                END DO
690                src_add (n) = -srch_add
691                src_lats(n) = distance
692                dist_min = src_lats(4)
693                EXIT sort_loop
694             ENDIF
695          END DO sort_loop
696       ENDIF
697    END DO
698    !
699    src_lons = one/(src_lats + tiny)
700    distance = SUM(src_lons)
701    src_lats = src_lons/distance
702
703    !-----------------------------------------------------------------------
704
705  END SUBROUTINE grid_search_bilin
706
707!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
708  !************************************************************************
709  !   SUBROUTINE INIT_REMAP_VARS
710  !************************************************************************
711!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
712  !
713  !      subroutine init_remap_vars
714  !
715  !-----------------------------------------------------------------------
716  !
717  !     this routine initializes some variables and provides an initial
718  !     allocation of arrays (fairly large so frequent resizing
719  !     unnecessary).
720  !
721  !-----------------------------------------------------------------------
722  !
723  !-----------------------------------------------------------------------
724  !     determine the number of weights
725  !-----------------------------------------------------------------------
726  !        num_wts = 1     ! bilinear interpolation
727  !-----------------------------------------------------------------------
728  !     initialize num_links and set max_links to four times the largest
729  !     of the destination grid sizes initially (can be changed later).
730  !     set a default resize increment to increase the size of link
731  !     arrays if the number of links exceeds the initial size 
732  !!-----------------------------------------------------------------------
733  !     
734  !      num_links_map1 = 0
735  !      max_links_map1 = 4*grid2_size
736  !      if (num_maps > 1) then
737  !        num_links_map2 = 0
738  !        max_links_map1 = max(4*grid1_size,4*grid2_size)
739  !        max_links_map2 = max_links_map1
740  !      endif
741  !
742  !      resize_increment = 0.1*max(grid1_size,grid2_size)
743  !
744  !-----------------------------------------------------------------------
745  !     allocate address and weight arrays for mapping 1 
746  !-----------------------------------------------------------------------
747  !
748  !      allocate (grid1_add_map1(max_links_map1),    &
749  !                grid2_add_map1(max_links_map1),    &
750  !                wts_map1(num_wts, max_links_map1))
751  !     
752  !   grid1_add_map1 = 0. 
753  !   grid2_add_map1 = 0.
754  !   wts_map1 = 0.!
755  !
756  !!-----------------------------------------------------------------------
757  !
758  !      end subroutine init_remap_vars
759
760END MODULE agrif_gridsearch
Note: See TracBrowser for help on using the repository browser.