source: branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/SCRIP1.4/source/grids.f @ 2352

Last change on this file since 2352 was 2352, checked in by sga, 11 years ago

NEMO branch nemo_v3_3_beta
Add NOCS tools based on SCRIP package for creating weights for interpolation on the fly
These now should build with the maketools script in the TOOLS directory using the same
architecture configuration file as the model (hopefully)

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