source: CPL/oasis3/trunk/src/lib/scrip/src/grids.f @ 1677

Last change on this file since 1677 was 1677, checked in by aclsce, 12 years ago

Imported oasis3 (tag ipslcm5a) from cvs server to svn server (igcmg project).
File size: 27.2 KB
Line 
1C****
2C                    ************************
3C                    *     OASIS MODULE     *
4C                    *     ------------     *
5C                    ************************
6C****
7C***********************************************************************
8C     This module belongs to the SCRIP library. It is modified to run
9C     within OASIS. 
10C     Modifications:
11C       - routine does not read SCRIP grid description files, but gets 
12C         arrays from the calling routine
13C       - unit conversion : only from degrees (OASIS unit) to radians
14C       - some allocated array will be freed in the end to allow multiple
15C         calls of SCRIP
16C       - map-methods and restriction-types are written in capital letters
17C       - added bin definition for reduced grid
18C
19C     Modified by            V. Gayler,  M&D                  20.09.2001
20C     Modified by            D. Declat,  CERFACS              27.06.2002
21C***********************************************************************
22!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23!
24!     This module reads in and initializes two grids for remapping.
25!     NOTE: grid1 must be the master grid -- the grid that determines
26!           which cells participate (e.g. land mask) and the fractional
27!           area of grid2 cells that participate in the remapping.
28!
29!-----------------------------------------------------------------------
30!
31!     CVS:$Id: grids.f,v 1.1.1.1 2005/03/23 16:01:11 adm Exp $
32!
33!     Copyright (c) 1997, 1998 the Regents of the University of
34!       California.
35!
36!     This software and ancillary information (herein called software)
37!     called SCRIP is made available under the terms described here. 
38!     The software has been approved for release with associated
39!     LA-CC Number 98-45.
40!
41!     Unless otherwise indicated, this software has been authored
42!     by an employee or employees of the University of California,
43!     operator of the Los Alamos National Laboratory under Contract
44!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
45!     Government has rights to use, reproduce, and distribute this
46!     software.  The public may copy and use this software without
47!     charge, provided that this Notice and any statement of authorship
48!     are reproduced on all copies.  Neither the Government nor the
49!     University makes any warranty, express or implied, or assumes
50!     any liability or responsibility for the use of this software.
51!
52!     If software is modified to produce derivative works, such modified
53!     software should be clearly marked, so as not to confuse it with
54!     the version available from Los Alamos National Laboratory.
55!
56!***********************************************************************
57
58      module grids
59
60!-----------------------------------------------------------------------
61
62      use kinds_mod    ! defines data types
63      use constants    ! common constants
64      use iounits      ! I/O unit manager
65
66      implicit none
67
68!-----------------------------------------------------------------------
69!
70!     variables that describe each grid
71!
72!-----------------------------------------------------------------------
73
74      integer (kind=int_kind), save ::
75     &             grid1_size, grid2_size, ! total points on each grid
76     &             grid1_rank, grid2_rank, ! rank of each grid
77     &             grid1_corners, grid2_corners ! number of corners
78                                                ! for each grid cell
79
80      integer (kind=int_kind), dimension(:), allocatable, save ::
81     &             grid1_dims, grid2_dims  ! size of each grid dimension
82
83      character(char_len), save :: 
84     &             grid1_name, grid2_name  ! name for each grid
85
86      character (char_len), save :: 
87     &             grid1_units, ! units for grid coords (degs/radians)
88     &             grid2_units  ! units for grid coords
89
90      real (kind=dbl_kind), parameter ::
91     &      deg2rad = pi/180.   ! conversion for deg to rads
92
93!-----------------------------------------------------------------------
94!
95!     grid coordinates and masks
96!
97!-----------------------------------------------------------------------
98
99      logical (kind=log_kind), dimension(:), allocatable, save ::
100     &             grid1_mask,        ! flag which cells participate
101     &             grid2_mask         ! flag which cells participate
102
103      real (kind=dbl_kind), dimension(:), allocatable, save ::
104     &             grid1_center_lat,  ! lat/lon coordinates for
105     &             grid1_center_lon,  ! each grid center in radians
106     &             grid2_center_lat, 
107     &             grid2_center_lon,
108     &             grid1_area,        ! tot area of each grid1 cell
109     &             grid2_area,        ! tot area of each grid2 cell
110     &             grid1_area_in,     ! area of grid1 cell from file
111     &             grid2_area_in,     ! area of grid2 cell from file
112     &             grid1_frac,        ! fractional area of grid cells
113     &             grid2_frac         ! participating in remapping
114
115      real (kind=dbl_kind), dimension(:,:), allocatable, save ::
116     &             grid1_corner_lat,  ! lat/lon coordinates for
117     &             grid1_corner_lon,  ! each grid corner in radians
118     &             grid2_corner_lat, 
119     &             grid2_corner_lon
120
121      logical (kind=log_kind), save ::
122     &             luse_grid_centers ! use centers for bounding boxes
123     &,            luse_grid1_area   ! use area from grid file
124     &,            luse_grid2_area   ! use area from grid file
125
126      real (kind=dbl_kind), dimension(:,:), allocatable, save ::
127     &             grid1_bound_box,  ! lat/lon bounding box for use
128     &             grid2_bound_box   ! in restricting grid searches
129
130!-----------------------------------------------------------------------
131!
132!     bins for restricting searches
133!
134!-----------------------------------------------------------------------
135
136      character (char_len), save ::
137     &        restrict_type  ! type of bins to use
138
139      integer (kind=int_kind), save ::
140     &        num_srch_bins,  ! num of bins for restricted srch
141     &        num_srch_red    ! num of bins for reduced case
142
143      integer (kind=int_kind), dimension(:,:), allocatable, save ::
144     &        bin_addr1, ! min,max adds for grid1 cells in this lat bin
145     &        bin_addr2  ! min,max adds for grid2 cells in this lat bin
146      integer (kind=int_kind), dimension(:,:), allocatable, save ::
147     &        bin_addr1_r  ! min,max adds for red grid1 cells
148
149      real(kind=dbl_kind), dimension(:,:), allocatable, save ::
150     &        bin_lats   ! min,max latitude for each search bin
151     &,       bin_lons   ! min,max longitude for each search bin
152      real(kind=dbl_kind), dimension(:,:), allocatable, save ::
153     &        bin_lats_r ! min,max lat for each search bin for red grid
154     &,       bin_lons_r ! min,max lon for each search bin for red grid
155
156!***********************************************************************
157
158      contains
159
160!***********************************************************************
161
162      subroutine grid_init(m_method, rst_type, n_srch_bins,
163     $                     src_size, dst_size, src_dims, dst_dims, 
164     &                     src_rank, dst_rank, ncrn_src, ncrn_dst,
165     &                     src_mask, dst_mask, src_name, dst_name,
166     &                     src_lat,  src_lon,  dst_lat,  dst_lon,
167     &                     src_corner_lat, src_corner_lon,
168     &                     dst_corner_lat, dst_corner_lon)
169
170!-----------------------------------------------------------------------
171!
172!     this routine gets grid info from routine scriprmp and makes any
173!     necessary changes (e.g. for 0,2pi longitude range)
174!
175!-----------------------------------------------------------------------
176
177!-----------------------------------------------------------------------
178!
179!     input variables
180!
181!-----------------------------------------------------------------------
182
183      integer (kind=int_kind), intent(in) ::
184     &     n_srch_bins,         ! num of bins for restricted srch
185     &     src_size,            ! source grid size
186     &     dst_size,            ! target grid size
187     &     src_rank,            ! source grid rank
188     &     dst_rank,            ! target grid rank
189     &     src_dims(src_rank),  ! source grid dimensions
190     &     dst_dims(dst_rank),  ! target grid dimensions
191     &     ncrn_src,            ! number of corners of a source grid cell
192     &     ncrn_dst,            ! number of corners of a target grid cell
193     &     src_mask(src_size),  ! source grid mask (master mask)
194     &     dst_mask(dst_size)   ! target grid mask
195
196      character*8, intent(in) ::
197     &     m_method,            ! remapping method
198     &     rst_type,            ! restriction type
199     &     src_name,            ! source grid name
200     &     dst_name             ! target grid name
201
202      real (kind=real_kind), intent (in) ::
203     &     src_lat(src_size), ! source grid latitudes
204     &     src_lon(src_size), ! sourde grid longitudes
205     &     dst_lat(dst_size),   ! target grid latitudes
206     &     dst_lon(dst_size),   ! target grid longitudes
207     &     src_corner_lat(ncrn_src,src_size), 
208     &     src_corner_lon(ncrn_src,src_size), 
209     &     dst_corner_lat(ncrn_dst,dst_size), 
210     &     dst_corner_lon(ncrn_dst,dst_size) 
211
212!-----------------------------------------------------------------------
213!
214!     local variables
215!
216!-----------------------------------------------------------------------
217
218      integer (kind=int_kind) :: 
219     &  n      ! loop counter
220     &, nele   ! element loop counter
221     &, i,j    ! logical 2d addresses
222     &, ip1,jp1
223     &, n_add, e_add, ne_add
224     &, nx, ny
225
226      real (kind=dbl_kind) :: 
227     &  dlat, dlon          ! lat/lon intervals for search bins
228
229      real (kind=dbl_kind), dimension(4) ::
230     &  tmp_lats, tmp_lons  ! temps for computing bounding boxes
231
232!-----------------------------------------------------------------------
233!
234!     allocate grid coordinates/masks and read data
235!
236!-----------------------------------------------------------------------
237     
238      select case(m_method)
239      case ('CONSERV')
240        luse_grid_centers = .false.
241      case ('BILINEAR')
242        luse_grid_centers = .true.
243      case ('BICUBIC')
244        luse_grid_centers = .true.
245      case ('DISTWGT')
246        luse_grid_centers = .true.
247      case ('GAUSWGT')
248        luse_grid_centers = .true.
249      case default
250        stop 'unknown mapping method'
251      end select
252
253      allocate( grid1_mask      (src_size),
254     &          grid2_mask      (dst_size),
255     &          grid1_center_lat(src_size), 
256     &          grid1_center_lon(src_size),
257     &          grid2_center_lat(dst_size), 
258     &          grid2_center_lon(dst_size),
259     &          grid1_area      (src_size),
260     &          grid2_area      (dst_size),
261     &          grid1_frac      (src_size),
262     &          grid2_frac      (dst_size),
263     &          grid1_dims      (src_rank),
264     &          grid2_dims      (dst_rank),
265     &          grid1_bound_box (4       , src_size),
266     &          grid2_bound_box (4       , dst_size))
267
268      if (.not. luse_grid_centers) then
269        allocate( grid1_corner_lat(ncrn_src, src_size),
270     &            grid1_corner_lon(ncrn_src, src_size),
271     &            grid2_corner_lat(ncrn_dst, dst_size),
272     &            grid2_corner_lon(ncrn_dst, dst_size))
273      endif
274
275!-----------------------------------------------------------------------
276!
277!     copy input data to module data
278!
279!-----------------------------------------------------------------------
280     
281      restrict_type    = rst_type
282      num_srch_bins    = n_srch_bins
283      grid1_size       = src_size
284      grid2_size       = dst_size
285      grid1_dims       = src_dims
286      grid2_dims       = dst_dims 
287      grid1_rank       = src_rank
288      grid2_rank       = dst_rank
289      grid1_corners    = ncrn_src
290      grid2_corners    = ncrn_dst
291      grid1_name       = src_name
292      grid2_name       = dst_name
293      grid1_center_lat = src_lat
294      grid1_center_lon = src_lon
295      grid2_center_lat = dst_lat
296      grid2_center_lon = dst_lon
297
298      if (.not. luse_grid_centers) then
299        grid1_corner_lat = src_corner_lat
300        grid1_corner_lon = src_corner_lon
301        grid2_corner_lat = dst_corner_lat
302        grid2_corner_lon = dst_corner_lon
303      endif
304
305c      if (luse_grid1_area) then
306c        grid1_area_in
307c      endif
308c      if (luse_grid2_area) then
309c        grid2_area_in
310c      endif
311
312      grid1_area = zero
313      grid1_frac = zero
314      grid2_area = zero
315      grid2_frac = zero
316
317
318!-----------------------------------------------------------------------
319!
320!     initialize logical mask and convert lat/lon units if required
321!
322!-----------------------------------------------------------------------
323      where (src_mask == 1)
324        grid1_mask = .true.
325      elsewhere
326        grid1_mask = .false.
327      endwhere
328
329      where (dst_mask == 1)
330        grid2_mask = .true.
331      elsewhere
332        grid2_mask = .false.
333      endwhere
334
335C
336C* -- convert unit from degrees (OASIS unit) to radians
337C
338      grid1_center_lat = grid1_center_lat*deg2rad
339      grid1_center_lon = grid1_center_lon*deg2rad
340      grid2_center_lat = grid2_center_lat*deg2rad
341      grid2_center_lon = grid2_center_lon*deg2rad
342
343      if (.not. luse_grid_centers) then
344        grid1_corner_lat = grid1_corner_lat*deg2rad
345        grid1_corner_lon = grid1_corner_lon*deg2rad
346        grid2_corner_lat = grid2_corner_lat*deg2rad
347        grid2_corner_lon = grid2_corner_lon*deg2rad
348      endif
349
350      grid1_units='radians'
351      grid2_units='radians'
352!-----------------------------------------------------------------------
353!
354!     convert longitudes to 0,2pi interval
355!
356!-----------------------------------------------------------------------
357
358      where (grid1_center_lon .gt. pi2)  grid1_center_lon =
359     &                                   grid1_center_lon - pi2
360      where (grid1_center_lon .lt. zero) grid1_center_lon =
361     &                                   grid1_center_lon + pi2
362      where (grid2_center_lon .gt. pi2)  grid2_center_lon =
363     &                                   grid2_center_lon - pi2
364      where (grid2_center_lon .lt. zero) grid2_center_lon =
365     &                                   grid2_center_lon + pi2
366
367      if (.not. luse_grid_centers) then
368        where (grid1_corner_lon .gt. pi2)  grid1_corner_lon =
369     &                                     grid1_corner_lon - pi2
370        where (grid1_corner_lon .lt. zero) grid1_corner_lon =
371     &                                     grid1_corner_lon + pi2
372        where (grid2_corner_lon .gt. pi2)  grid2_corner_lon =
373     &                                     grid2_corner_lon - pi2
374        where (grid2_corner_lon .lt. zero) grid2_corner_lon =
375     &                                     grid2_corner_lon + pi2
376      endif
377!-----------------------------------------------------------------------
378!
379!     make sure input latitude range is within the machine values
380!     for +/- pi/2
381!
382!-----------------------------------------------------------------------
383
384      where (grid1_center_lat >  pih) grid1_center_lat =  pih
385      where (grid1_center_lat < -pih) grid1_center_lat = -pih
386      where (grid2_center_lat >  pih) grid2_center_lat =  pih
387      where (grid2_center_lat < -pih) grid2_center_lat = -pih
388
389      if (.not. luse_grid_centers) then
390        where (grid1_corner_lat >  pih) grid1_corner_lat =  pih
391        where (grid1_corner_lat < -pih) grid1_corner_lat = -pih
392        where (grid2_corner_lat >  pih) grid2_corner_lat =  pih
393        where (grid2_corner_lat < -pih) grid2_corner_lat = -pih
394      endif
395
396!-----------------------------------------------------------------------
397!
398!     compute bounding boxes for restricting future grid searches
399!
400!-----------------------------------------------------------------------
401
402      if (.not. luse_grid_centers) then
403
404        grid1_bound_box(1,:) = minval(grid1_corner_lat, DIM=1)
405        grid1_bound_box(2,:) = maxval(grid1_corner_lat, DIM=1)
406        grid1_bound_box(3,:) = minval(grid1_corner_lon, DIM=1)
407        grid1_bound_box(4,:) = maxval(grid1_corner_lon, DIM=1)
408
409        grid2_bound_box(1,:) = minval(grid2_corner_lat, DIM=1)
410        grid2_bound_box(2,:) = maxval(grid2_corner_lat, DIM=1)
411        grid2_bound_box(3,:) = minval(grid2_corner_lon, DIM=1)
412        grid2_bound_box(4,:) = maxval(grid2_corner_lon, DIM=1)
413      else
414
415        nx = grid1_dims(1)
416        ny = grid1_dims(2)
417
418        do n=1,grid1_size
419
420          !*** find N,S and NE points to this grid point
421
422          j = (n - 1)/nx +1
423          i = n - (j-1)*nx
424
425          if (i < nx) then
426            ip1 = i + 1
427          else
428            !*** assume cyclic
429            ip1 = 1
430            !*** but if it is not, correct
431            e_add = (j - 1)*nx + ip1
432            if (abs(grid1_center_lat(e_add) - 
433     &              grid1_center_lat(n   )) > pih) then
434              ip1 = i
435            endif
436          endif
437
438          if (j < ny) then
439            jp1 = j+1
440          else
441            !*** assume cyclic
442            jp1 = 1
443            !*** but if it is not, correct
444            n_add = (jp1 - 1)*nx + i
445            if (abs(grid1_center_lat(n_add) - 
446     &              grid1_center_lat(n   )) > pih) then
447              jp1 = j
448            endif
449          endif
450
451          n_add = (jp1 - 1)*nx + i
452          e_add = (j - 1)*nx + ip1
453          ne_add = (jp1 - 1)*nx + ip1
454
455          !*** find N,S and NE lat/lon coords and check bounding box
456
457          tmp_lats(1) = grid1_center_lat(n)
458          tmp_lats(2) = grid1_center_lat(e_add)
459          tmp_lats(3) = grid1_center_lat(ne_add)
460          tmp_lats(4) = grid1_center_lat(n_add)
461
462          tmp_lons(1) = grid1_center_lon(n)
463          tmp_lons(2) = grid1_center_lon(e_add)
464          tmp_lons(3) = grid1_center_lon(ne_add)
465          tmp_lons(4) = grid1_center_lon(n_add)
466
467          grid1_bound_box(1,n) = minval(tmp_lats)
468          grid1_bound_box(2,n) = maxval(tmp_lats)
469          grid1_bound_box(3,n) = minval(tmp_lons)
470          grid1_bound_box(4,n) = maxval(tmp_lons)
471        end do
472
473        nx = grid2_dims(1)
474        ny = grid2_dims(2)
475
476        do n=1,grid2_size
477
478          !*** find N,S and NE points to this grid point
479
480          j = (n - 1)/nx +1
481          i = n - (j-1)*nx
482
483          if (i < nx) then
484            ip1 = i + 1
485          else
486            !*** assume cyclic
487            ip1 = 1
488            !*** but if it is not, correct
489            e_add = (j - 1)*nx + ip1
490            if (abs(grid2_center_lat(e_add) - 
491     &              grid2_center_lat(n   )) > pih) then
492              ip1 = i
493            endif
494          endif
495
496          if (j < ny) then
497            jp1 = j+1
498          else
499            !*** assume cyclic
500            jp1 = 1
501            !*** but if it is not, correct
502            n_add = (jp1 - 1)*nx + i
503            if (abs(grid2_center_lat(n_add) - 
504     &              grid2_center_lat(n   )) > pih) then
505              jp1 = j
506            endif
507          endif
508
509          n_add = (jp1 - 1)*nx + i
510          e_add = (j - 1)*nx + ip1
511          ne_add = (jp1 - 1)*nx + ip1
512
513          !*** find N,S and NE lat/lon coords and check bounding box
514
515          tmp_lats(1) = grid2_center_lat(n)
516          tmp_lats(2) = grid2_center_lat(e_add)
517          tmp_lats(3) = grid2_center_lat(ne_add)
518          tmp_lats(4) = grid2_center_lat(n_add)
519
520          tmp_lons(1) = grid2_center_lon(n)
521          tmp_lons(2) = grid2_center_lon(e_add)
522          tmp_lons(3) = grid2_center_lon(ne_add)
523          tmp_lons(4) = grid2_center_lon(n_add)
524
525          grid2_bound_box(1,n) = minval(tmp_lats)
526          grid2_bound_box(2,n) = maxval(tmp_lats)
527          grid2_bound_box(3,n) = minval(tmp_lons)
528          grid2_bound_box(4,n) = maxval(tmp_lons)
529        end do
530
531      endif
532
533      where (abs(grid1_bound_box(4,:) - grid1_bound_box(3,:)) > pi)
534        grid1_bound_box(3,:) = zero
535        grid1_bound_box(4,:) = pi2
536      end where
537
538      where (abs(grid2_bound_box(4,:) - grid2_bound_box(3,:)) > pi)
539        grid2_bound_box(3,:) = zero
540        grid2_bound_box(4,:) = pi2
541      end where
542
543      !***
544      !*** try to check for cells that overlap poles
545      !***
546
547      where (grid1_center_lat > grid1_bound_box(2,:))
548     &  grid1_bound_box(2,:) = pih
549
550      where (grid1_center_lat < grid1_bound_box(1,:))
551     &  grid1_bound_box(1,:) = -pih
552
553      where (grid2_center_lat > grid2_bound_box(2,:))
554     &  grid2_bound_box(2,:) = pih
555
556      where (grid2_center_lat < grid2_bound_box(1,:))
557     &  grid2_bound_box(1,:) = -pih
558
559!-----------------------------------------------------------------------
560!
561!     set up and assign address ranges to search bins in order to
562!     further restrict later searches
563!
564!-----------------------------------------------------------------------
565
566      select case (restrict_type)
567
568      case ('LATITUDE')
569        write(stdout,*) 'Using latitude bins to restrict search.'
570
571        allocate(bin_addr1(2,num_srch_bins))
572        allocate(bin_addr2(2,num_srch_bins))
573        allocate(bin_lats (2,num_srch_bins))
574        allocate(bin_lons (2,num_srch_bins))
575
576        dlat = pi/num_srch_bins
577
578        do n=1,num_srch_bins
579          bin_lats(1,n) = (n-1)*dlat - pih
580          bin_lats(2,n) =     n*dlat - pih
581          bin_lons(1,n) = zero
582          bin_lons(2,n) = pi2
583          bin_addr1(1,n) = grid1_size + 1
584          bin_addr1(2,n) = 0
585          bin_addr2(1,n) = grid2_size + 1
586          bin_addr2(2,n) = 0
587        end do
588
589        do nele=1,grid1_size
590          do n=1,num_srch_bins
591            if (grid1_bound_box(1,nele) <= bin_lats(2,n) .and.
592     &          grid1_bound_box(2,nele) >= bin_lats(1,n)) then
593              bin_addr1(1,n) = min(nele,bin_addr1(1,n))
594              bin_addr1(2,n) = max(nele,bin_addr1(2,n))
595            endif
596          end do
597        end do
598
599        do nele=1,grid2_size
600          do n=1,num_srch_bins
601            if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and.
602     &          grid2_bound_box(2,nele) >= bin_lats(1,n)) then
603              bin_addr2(1,n) = min(nele,bin_addr2(1,n))
604              bin_addr2(2,n) = max(nele,bin_addr2(2,n))
605            endif
606          end do
607        end do
608
609      case ('LATLON')
610        write(stdout,*) 'Using lat/lon boxes to restrict search.'
611
612        dlat = pi /num_srch_bins
613        dlon = pi2/num_srch_bins
614
615        allocate(bin_addr1(2,num_srch_bins*num_srch_bins))
616        allocate(bin_addr2(2,num_srch_bins*num_srch_bins))
617        allocate(bin_lats (2,num_srch_bins*num_srch_bins))
618        allocate(bin_lons (2,num_srch_bins*num_srch_bins))
619
620        n = 0
621        do j=1,num_srch_bins
622          do i=1,num_srch_bins
623            n = n + 1
624           
625            bin_lats(1,n) = (j-1)*dlat - pih
626            bin_lats(2,n) =     j*dlat - pih
627            bin_lons(1,n) = (i-1)*dlon
628            bin_lons(2,n) =     i*dlon
629            bin_addr1(1,n) = grid1_size + 1
630            bin_addr1(2,n) = 0
631            bin_addr2(1,n) = grid2_size + 1
632            bin_addr2(2,n) = 0
633          end do
634        end do
635       
636        num_srch_bins = num_srch_bins**2
637
638        do nele=1,grid1_size
639          do n=1,num_srch_bins
640            if (grid1_bound_box(1,nele) <= bin_lats(2,n) .and.
641     &          grid1_bound_box(2,nele) >= bin_lats(1,n) .and.
642     &          grid1_bound_box(3,nele) <= bin_lons(2,n) .and.
643     &          grid1_bound_box(4,nele) >= bin_lons(1,n)) then
644                bin_addr1(1,n) = min(nele,bin_addr1(1,n))
645                bin_addr1(2,n) = max(nele,bin_addr1(2,n))
646            endif
647          end do
648        end do
649
650        do nele=1,grid2_size
651          do n=1,num_srch_bins
652            if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and.
653     &          grid2_bound_box(2,nele) >= bin_lats(1,n) .and.
654     &          grid2_bound_box(3,nele) <= bin_lons(2,n) .and.
655     &          grid2_bound_box(4,nele) >= bin_lons(1,n)) then
656                bin_addr2(1,n) = min(nele,bin_addr2(1,n))
657                bin_addr2(2,n) = max(nele,bin_addr2(2,n))
658            endif
659          end do
660        end do
661
662      case ('REDUCED')
663        write(stdout,*) 
664     &  '|  Using reduced bins to restrict search. Reduced grids with'
665        write(stdout,*) 
666     &  '| a maximum of 500*NBRBINS latitude circles can be treated'
667
668        allocate(bin_addr2(2,num_srch_bins))
669        allocate(bin_lats (2,num_srch_bins))
670        allocate(bin_lons (2,num_srch_bins))
671        allocate(bin_addr1_r(4,500*num_srch_bins))
672        allocate(bin_lats_r (2,500*num_srch_bins))
673        allocate(bin_lons_r (2,500*num_srch_bins))
674
675        dlat = pi/num_srch_bins
676
677        do n=1,num_srch_bins
678          bin_lats(1,n) = (n-1)*dlat - pih
679          bin_lats(2,n) =     n*dlat - pih
680          bin_lons(1,n) = zero
681          bin_lons(2,n) = pi2
682          bin_addr2(1,n) = grid2_size + 1
683          bin_addr2(2,n) = 0
684        end do
685
686        do nele=1,grid2_size
687          do n=1,num_srch_bins
688            if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and.
689     &          grid2_bound_box(2,nele) >= bin_lats(1,n)) then
690              bin_addr2(1,n) = min(nele,bin_addr2(1,n))
691              bin_addr2(2,n) = max(nele,bin_addr2(2,n))
692            endif
693          end do
694        end DO
695
696
697        bin_addr1_r(1,1) = 0
698        bin_lats_r(1,1) = grid1_center_lat(1)
699        num_srch_red = 1
700        do nele=1,grid1_size-1
701          if (grid1_center_lat(nele+1) == 
702     &        grid1_center_lat(nele)) THEN
703              bin_addr1_r(2,num_srch_red) = nele+1
704          else IF (grid1_center_lat(nele+1) /= 
705     &        grid1_center_lat(nele)) THEN
706              bin_addr1_r(1,num_srch_red+1) = nele+1
707              bin_lats_r(2,num_srch_red) = grid1_center_lat(nele+1)
708              bin_lats_r(1,num_srch_red+1) = 
709     &            bin_lats_r(2,num_srch_red)
710              num_srch_red = num_srch_red+1
711          ENDIF
712        end DO
713       
714        DO nele = 1,num_srch_red-1
715          bin_addr1_r(3,nele)=bin_addr1_r(1,nele+1)
716          bin_addr1_r(4,nele)=bin_addr1_r(2,nele+1)
717        enddo
718        bin_addr1_r(3,num_srch_red) = bin_addr1_r(1,num_srch_red-1)
719        bin_addr1_r(4,num_srch_red) = bin_addr1_r(2,num_srch_red-1)
720
721      case default
722        stop 'unknown search restriction method'
723      end select
724 
725!-----------------------------------------------------------------------
726
727      end subroutine grid_init
728
729!***********************************************************************
730
731      subroutine free_grids
732
733!-----------------------------------------------------------------------
734!
735!     subroutine to deallocate allocated arrays
736!
737!-----------------------------------------------------------------------
738
739      deallocate(grid1_mask, grid2_mask, 
740     &    grid1_center_lat, grid1_center_lon,
741     &    grid2_center_lat, grid2_center_lon,
742     &    grid1_area, grid2_area,
743     &    grid1_frac, grid2_frac,
744     &    grid1_dims, grid2_dims)
745
746      IF (restrict_TYPE == 'REDUCED') then
747          deallocate( grid1_bound_box, grid2_bound_box,
748     &        bin_addr1_r, bin_addr2,
749     &        bin_lons, bin_lats,
750     &        bin_lats_r, bin_lons_r)
751      else
752          deallocate( grid1_bound_box, grid2_bound_box,
753     &        bin_addr1, bin_addr2,
754     &        bin_lats, bin_lons)
755      endif
756
757      if (.not. luse_grid_centers) then
758        deallocate(grid1_corner_lat, grid1_corner_lon,
759     &             grid2_corner_lat, grid2_corner_lon)
760      endif
761
762!-----------------------------------------------------------------------
763
764      end subroutine free_grids
765
766!***********************************************************************
767
768      end module grids
769
770!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
771
Note: See TracBrowser for help on using the repository browser.