source: branches/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/WEIGHTS/src/grids.f90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 29.7 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$
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.