source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/lib/scrip/src/grids.f @ 4775

Last change on this file since 4775 was 4775, checked in by aclsce, 5 years ago
  • Imported oasis3-mct from Cerfacs svn server (not suppotred anymore).

The version has been extracted from ​https://oasis3mct.cerfacs.fr/svn/branches/OASIS3-MCT_2.0_branch/oasis3-mct@1818
File size: 28.5 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 2826 2010-12-10 11:14:21Z valcke $
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      USE mod_oasis_flush
66
67      implicit none
68
69!-----------------------------------------------------------------------
70!
71!     variables that describe each grid
72!
73!-----------------------------------------------------------------------
74
75      integer (kind=int_kind), save ::
76     &             grid1_size, grid2_size, ! total points on each grid
77     &             grid1_rank, grid2_rank, ! rank of each grid
78     &             grid1_corners, grid2_corners ! number of corners
79                                                ! for each grid cell
80
81      integer (kind=int_kind), dimension(:), allocatable, save ::
82     &             grid1_dims, grid2_dims  ! size of each grid dimension
83
84      character(char_len), save :: 
85     &             grid1_name, grid2_name  ! name for each grid
86
87      character (char_len), save :: 
88     &             grid1_units, ! units for grid coords (degs/radians)
89     &             grid2_units  ! units for grid coords
90
91      real (kind=dbl_kind), parameter ::
92     &      deg2rad = pi/180.   ! conversion for deg to rads
93
94!-----------------------------------------------------------------------
95!
96!     grid coordinates and masks
97!
98!-----------------------------------------------------------------------
99
100      logical (kind=log_kind), dimension(:), allocatable, save ::
101     &             grid1_mask,        ! flag which cells participate
102     &             grid2_mask         ! flag which cells participate
103
104      real (kind=dbl_kind), dimension(:), allocatable, save ::
105     &             grid1_center_lat,  ! lat/lon coordinates for
106     &             grid1_center_lon,  ! each grid center in radians
107     &             grid2_center_lat, 
108     &             grid2_center_lon,
109     &             grid1_area,        ! tot area of each grid1 cell
110     &             grid2_area,        ! tot area of each grid2 cell
111     &             grid1_area_in,     ! area of grid1 cell from file
112     &             grid2_area_in,     ! area of grid2 cell from file
113     &             grid1_frac,        ! fractional area of grid cells
114     &             grid2_frac         ! participating in remapping
115
116      real (kind=dbl_kind), dimension(:,:), allocatable, save ::
117     &             grid1_corner_lat,  ! lat/lon coordinates for
118     &             grid1_corner_lon,  ! each grid corner in radians
119     &             grid2_corner_lat, 
120     &             grid2_corner_lon
121
122      logical (kind=log_kind), save ::
123     &             luse_grid_centers ! use centers for bounding boxes
124     &,            luse_grid1_area   ! use area from grid file
125     &,            luse_grid2_area   ! use area from grid file
126
127      real (kind=dbl_kind), dimension(:,:), allocatable, save ::
128     &             grid1_bound_box,  ! lat/lon bounding box for use
129     &             grid2_bound_box   ! in restricting grid searches
130
131!-----------------------------------------------------------------------
132!
133!     bins for restricting searches
134!
135!-----------------------------------------------------------------------
136
137      character (char_len), save ::
138     &        restrict_type  ! type of bins to use
139
140      integer (kind=int_kind), save ::
141     &        num_srch_bins,  ! num of bins for restricted srch
142     &        num_srch_red    ! num of bins for reduced case
143
144      integer (kind=int_kind), dimension(:,:), allocatable, save ::
145     &        bin_addr1, ! min,max adds for grid1 cells in this lat bin
146     &        bin_addr2  ! min,max adds for grid2 cells in this lat bin
147      integer (kind=int_kind), dimension(:,:), allocatable, save ::
148     &        bin_addr1_r  ! min,max adds for red grid1 cells
149
150      real(kind=dbl_kind), dimension(:,:), allocatable, save ::
151     &        bin_lats   ! min,max latitude for each search bin
152     &,       bin_lons   ! min,max longitude for each search bin
153      real(kind=dbl_kind), dimension(:,:), allocatable, save ::
154     &        bin_lats_r ! min,max lat for each search bin for red grid
155     &,       bin_lons_r ! min,max lon for each search bin for red grid
156
157!***********************************************************************
158
159      contains
160
161!***********************************************************************
162
163      subroutine grid_init(m_method, rst_type, n_srch_bins,
164     $                     src_size, dst_size, src_dims, dst_dims, 
165     &                     src_rank, dst_rank, ncrn_src, ncrn_dst,
166     &                     src_mask, dst_mask, src_name, dst_name,
167     &                     src_lat,  src_lon,  dst_lat,  dst_lon,
168     &                     src_corner_lat, src_corner_lon,
169     &                     dst_corner_lat, dst_corner_lon, 
170     &                     logunit)
171
172!-----------------------------------------------------------------------
173!
174!     this routine gets grid info from routine scriprmp and makes any
175!     necessary changes (e.g. for 0,2pi longitude range)
176!
177!-----------------------------------------------------------------------
178
179!-----------------------------------------------------------------------
180!
181!     input variables
182!
183!-----------------------------------------------------------------------
184
185      integer (kind=int_kind), intent(in) ::
186     &     n_srch_bins,         ! num of bins for restricted srch
187     &     src_size,            ! source grid size
188     &     dst_size,            ! target grid size
189     &     src_rank,            ! source grid rank
190     &     dst_rank,            ! target grid rank
191     &     src_dims(src_rank),  ! source grid dimensions
192     &     dst_dims(dst_rank),  ! target grid dimensions
193     &     ncrn_src,            ! number of corners of a source grid cell
194     &     ncrn_dst,            ! number of corners of a target grid cell
195     &     src_mask(src_size),  ! source grid mask (master mask)
196     &     dst_mask(dst_size)   ! target grid mask
197
198      integer(kind=int_kind), intent(in), optional ::
199     &     logunit
200
201      character*8, intent(in) ::
202     &     m_method,            ! remapping method
203     &     rst_type,            ! restriction type
204     &     src_name,            ! source grid name
205     &     dst_name             ! target grid name
206
207      real (kind=real_kind), intent (in) ::
208     &     src_lat(src_size), ! source grid latitudes
209     &     src_lon(src_size), ! sourde grid longitudes
210     &     dst_lat(dst_size),   ! target grid latitudes
211     &     dst_lon(dst_size),   ! target grid longitudes
212     &     src_corner_lat(ncrn_src,src_size), 
213     &     src_corner_lon(ncrn_src,src_size), 
214     &     dst_corner_lat(ncrn_dst,dst_size), 
215     &     dst_corner_lon(ncrn_dst,dst_size) 
216
217!-----------------------------------------------------------------------
218!
219!     local variables
220!
221!-----------------------------------------------------------------------
222
223      integer (kind=int_kind) :: 
224     &  n      ! loop counter
225     &, nele   ! element loop counter
226     &, i,j    ! logical 2d addresses
227     &, ip1,jp1
228     &, n_add, e_add, ne_add
229     &, nx, ny
230
231      real (kind=dbl_kind) :: 
232     &  dlat, dlon          ! lat/lon intervals for search bins
233
234      real (kind=dbl_kind), dimension(4) ::
235     &  tmp_lats, tmp_lons  ! temps for computing bounding boxes
236      character(len=*),parameter :: subname = 'scrip:grid_init '
237!
238!-----------------------------------------------------------------------
239!
240      if (present(logunit)) then
241         nulou = logunit
242      endif
243
244      IF (nlogprt .GE. 2) THEN
245         WRITE (UNIT = nulou,FMT = *)' '
246         WRITE (UNIT = nulou,FMT = *)'Entering routine grid_init'
247         WRITE (UNIT = nulou,FMT = *)' '
248         CALL OASIS_FLUSH_SCRIP(nulou)
249      ENDIF
250!
251!-----------------------------------------------------------------------
252!
253!     allocate grid coordinates/masks and read data
254!
255!-----------------------------------------------------------------------
256
257!      write(nulou,*) subname,trim(m_method)
258!      write(nulou,*) subname,src_size,dst_size,src_rank,dst_rank
259
260      select case(m_method)
261      case ('CONSERV')
262        luse_grid_centers = .false.
263      case ('BILINEAR')
264        luse_grid_centers = .true.
265      case ('BICUBIC')
266        luse_grid_centers = .true.
267      case ('DISTWGT')
268        luse_grid_centers = .true.
269      case ('GAUSWGT')
270        luse_grid_centers = .true.
271      case default
272        stop 'unknown mapping method'
273      end select
274
275      allocate( grid1_mask      (src_size),
276     &          grid2_mask      (dst_size),
277     &          grid1_center_lat(src_size), 
278     &          grid1_center_lon(src_size),
279     &          grid2_center_lat(dst_size), 
280     &          grid2_center_lon(dst_size),
281     &          grid1_area      (src_size),
282     &          grid2_area      (dst_size),
283     &          grid1_frac      (src_size),
284     &          grid2_frac      (dst_size),
285     &          grid1_dims      (src_rank),
286     &          grid2_dims      (dst_rank),
287     &          grid1_bound_box (4       , src_size),
288     &          grid2_bound_box (4       , dst_size))
289
290      if (.not. luse_grid_centers) then
291        allocate( grid1_corner_lat(ncrn_src, src_size),
292     &            grid1_corner_lon(ncrn_src, src_size),
293     &            grid2_corner_lat(ncrn_dst, dst_size),
294     &            grid2_corner_lon(ncrn_dst, dst_size))
295      endif
296
297!-----------------------------------------------------------------------
298!
299!     copy input data to module data
300!
301!-----------------------------------------------------------------------
302
303      restrict_type    = rst_type
304      num_srch_bins    = n_srch_bins
305      grid1_size       = src_size
306      grid2_size       = dst_size
307      grid1_dims       = src_dims
308      grid2_dims       = dst_dims 
309      grid1_rank       = src_rank
310      grid2_rank       = dst_rank
311      grid1_corners    = ncrn_src
312      grid2_corners    = ncrn_dst
313      grid1_name       = src_name
314      grid2_name       = dst_name
315      grid1_center_lat = src_lat
316      grid1_center_lon = src_lon
317      grid2_center_lat = dst_lat
318      grid2_center_lon = dst_lon
319
320      if (.not. luse_grid_centers) then
321        grid1_corner_lat = src_corner_lat
322        grid1_corner_lon = src_corner_lon
323        grid2_corner_lat = dst_corner_lat
324        grid2_corner_lon = dst_corner_lon
325      endif
326
327c      if (luse_grid1_area) then
328c        grid1_area_in
329c      endif
330c      if (luse_grid2_area) then
331c        grid2_area_in
332c      endif
333
334      grid1_area = zero
335      grid1_frac = zero
336      grid2_area = zero
337      grid2_frac = zero
338
339!-----------------------------------------------------------------------
340!
341!     initialize logical mask and convert lat/lon units if required
342!
343!-----------------------------------------------------------------------
344      where (src_mask == 1)
345        grid1_mask = .true.
346      elsewhere
347        grid1_mask = .false.
348      endwhere
349
350      where (dst_mask == 1)
351        grid2_mask = .true.
352      elsewhere
353        grid2_mask = .false.
354      endwhere
355
356C
357C* -- convert unit from degrees (OASIS unit) to radians
358C
359      grid1_center_lat = grid1_center_lat*deg2rad
360      grid1_center_lon = grid1_center_lon*deg2rad
361      grid2_center_lat = grid2_center_lat*deg2rad
362      grid2_center_lon = grid2_center_lon*deg2rad
363
364      if (.not. luse_grid_centers) then
365        grid1_corner_lat = grid1_corner_lat*deg2rad
366        grid1_corner_lon = grid1_corner_lon*deg2rad
367        grid2_corner_lat = grid2_corner_lat*deg2rad
368        grid2_corner_lon = grid2_corner_lon*deg2rad
369      endif
370
371      grid1_units='radians'
372      grid2_units='radians'
373!-----------------------------------------------------------------------
374!
375!     convert longitudes to 0,2pi interval
376!
377!-----------------------------------------------------------------------
378
379      where (grid1_center_lon .gt. pi2)  grid1_center_lon =
380     &                                   grid1_center_lon - pi2
381      where (grid1_center_lon .lt. zero) grid1_center_lon =
382     &                                   grid1_center_lon + pi2
383      where (grid2_center_lon .gt. pi2)  grid2_center_lon =
384     &                                   grid2_center_lon - pi2
385      where (grid2_center_lon .lt. zero) grid2_center_lon =
386     &                                   grid2_center_lon + pi2
387
388      if (.not. luse_grid_centers) then
389        where (grid1_corner_lon .gt. pi2)  grid1_corner_lon =
390     &                                     grid1_corner_lon - pi2
391        where (grid1_corner_lon .lt. zero) grid1_corner_lon =
392     &                                     grid1_corner_lon + pi2
393        where (grid2_corner_lon .gt. pi2)  grid2_corner_lon =
394     &                                     grid2_corner_lon - pi2
395        where (grid2_corner_lon .lt. zero) grid2_corner_lon =
396     &                                     grid2_corner_lon + pi2
397      endif
398!-----------------------------------------------------------------------
399!
400!     make sure input latitude range is within the machine values
401!     for +/- pi/2
402!
403!-----------------------------------------------------------------------
404
405      where (grid1_center_lat >  pih) grid1_center_lat =  pih
406      where (grid1_center_lat < -pih) grid1_center_lat = -pih
407      where (grid2_center_lat >  pih) grid2_center_lat =  pih
408      where (grid2_center_lat < -pih) grid2_center_lat = -pih
409
410      if (.not. luse_grid_centers) then
411        where (grid1_corner_lat >  pih) grid1_corner_lat =  pih
412        where (grid1_corner_lat < -pih) grid1_corner_lat = -pih
413        where (grid2_corner_lat >  pih) grid2_corner_lat =  pih
414        where (grid2_corner_lat < -pih) grid2_corner_lat = -pih
415      endif
416
417!-----------------------------------------------------------------------
418!
419!     compute bounding boxes for restricting future grid searches
420!
421!-----------------------------------------------------------------------
422
423      if (.not. luse_grid_centers) then
424
425        grid1_bound_box(1,:) = minval(grid1_corner_lat, DIM=1)
426        grid1_bound_box(2,:) = maxval(grid1_corner_lat, DIM=1)
427        grid1_bound_box(3,:) = minval(grid1_corner_lon, DIM=1)
428        grid1_bound_box(4,:) = maxval(grid1_corner_lon, DIM=1)
429
430        grid2_bound_box(1,:) = minval(grid2_corner_lat, DIM=1)
431        grid2_bound_box(2,:) = maxval(grid2_corner_lat, DIM=1)
432        grid2_bound_box(3,:) = minval(grid2_corner_lon, DIM=1)
433        grid2_bound_box(4,:) = maxval(grid2_corner_lon, DIM=1)
434      else
435
436        nx = grid1_dims(1)
437        ny = grid1_dims(2)
438
439        do n=1,grid1_size
440
441          !*** find N,S and NE points to this grid point
442
443          j = (n - 1)/nx +1
444          i = n - (j-1)*nx
445
446          if (i < nx) then
447            ip1 = i + 1
448          else
449            !*** assume cyclic
450            ip1 = 1
451            !*** but if it is not, correct
452            e_add = (j - 1)*nx + ip1
453            if (abs(grid1_center_lat(e_add) - 
454     &              grid1_center_lat(n   )) > pih) then
455              ip1 = i
456            endif
457          endif
458
459          if (j < ny) then
460            jp1 = j+1
461          else
462            !*** assume cyclic
463            jp1 = 1
464            !*** but if it is not, correct
465            n_add = (jp1 - 1)*nx + i
466            if (abs(grid1_center_lat(n_add) - 
467     &              grid1_center_lat(n   )) > pih) then
468              jp1 = j
469            endif
470          endif
471
472          n_add = (jp1 - 1)*nx + i
473          e_add = (j - 1)*nx + ip1
474          ne_add = (jp1 - 1)*nx + ip1
475
476          !*** find N,S and NE lat/lon coords and check bounding box
477
478          tmp_lats(1) = grid1_center_lat(n)
479          tmp_lats(2) = grid1_center_lat(e_add)
480          tmp_lats(3) = grid1_center_lat(ne_add)
481          tmp_lats(4) = grid1_center_lat(n_add)
482
483          tmp_lons(1) = grid1_center_lon(n)
484          tmp_lons(2) = grid1_center_lon(e_add)
485          tmp_lons(3) = grid1_center_lon(ne_add)
486          tmp_lons(4) = grid1_center_lon(n_add)
487
488          grid1_bound_box(1,n) = minval(tmp_lats)
489          grid1_bound_box(2,n) = maxval(tmp_lats)
490          grid1_bound_box(3,n) = minval(tmp_lons)
491          grid1_bound_box(4,n) = maxval(tmp_lons)
492        end do
493
494        nx = grid2_dims(1)
495        ny = grid2_dims(2)
496
497        do n=1,grid2_size
498
499          !*** find N,S and NE points to this grid point
500
501          j = (n - 1)/nx +1
502          i = n - (j-1)*nx
503
504          if (i < nx) then
505            ip1 = i + 1
506          else
507            !*** assume cyclic
508            ip1 = 1
509            !*** but if it is not, correct
510            e_add = (j - 1)*nx + ip1
511            if (abs(grid2_center_lat(e_add) - 
512     &              grid2_center_lat(n   )) > pih) then
513              ip1 = i
514            endif
515          endif
516
517          if (j < ny) then
518            jp1 = j+1
519          else
520            !*** assume cyclic
521            jp1 = 1
522            !*** but if it is not, correct
523            n_add = (jp1 - 1)*nx + i
524            if (abs(grid2_center_lat(n_add) - 
525     &              grid2_center_lat(n   )) > pih) then
526              jp1 = j
527            endif
528          endif
529
530          n_add = (jp1 - 1)*nx + i
531          e_add = (j - 1)*nx + ip1
532          ne_add = (jp1 - 1)*nx + ip1
533
534          !*** find N,S and NE lat/lon coords and check bounding box
535
536          tmp_lats(1) = grid2_center_lat(n)
537          tmp_lats(2) = grid2_center_lat(e_add)
538          tmp_lats(3) = grid2_center_lat(ne_add)
539          tmp_lats(4) = grid2_center_lat(n_add)
540
541          tmp_lons(1) = grid2_center_lon(n)
542          tmp_lons(2) = grid2_center_lon(e_add)
543          tmp_lons(3) = grid2_center_lon(ne_add)
544          tmp_lons(4) = grid2_center_lon(n_add)
545
546          grid2_bound_box(1,n) = minval(tmp_lats)
547          grid2_bound_box(2,n) = maxval(tmp_lats)
548          grid2_bound_box(3,n) = minval(tmp_lons)
549          grid2_bound_box(4,n) = maxval(tmp_lons)
550        end do
551
552      endif
553
554      where (abs(grid1_bound_box(4,:) - grid1_bound_box(3,:)) > pi)
555        grid1_bound_box(3,:) = zero
556        grid1_bound_box(4,:) = pi2
557      end where
558
559      where (abs(grid2_bound_box(4,:) - grid2_bound_box(3,:)) > pi)
560        grid2_bound_box(3,:) = zero
561        grid2_bound_box(4,:) = pi2
562      end where
563
564      !***
565      !*** try to check for cells that overlap poles
566      !***
567
568      where (grid1_center_lat > grid1_bound_box(2,:))
569     &  grid1_bound_box(2,:) = pih
570
571      where (grid1_center_lat < grid1_bound_box(1,:))
572     &  grid1_bound_box(1,:) = -pih
573
574      where (grid2_center_lat > grid2_bound_box(2,:))
575     &  grid2_bound_box(2,:) = pih
576
577      where (grid2_center_lat < grid2_bound_box(1,:))
578     &  grid2_bound_box(1,:) = -pih
579
580!-----------------------------------------------------------------------
581!
582!     set up and assign address ranges to search bins in order to
583!     further restrict later searches
584!
585!-----------------------------------------------------------------------
586
587      select case (restrict_type)
588
589      case ('LATITUDE')
590        write(nulou,*) 'Using latitude bins to restrict search.'
591
592        allocate(bin_addr1(2,num_srch_bins))
593        allocate(bin_addr2(2,num_srch_bins))
594        allocate(bin_lats (2,num_srch_bins))
595        allocate(bin_lons (2,num_srch_bins))
596
597        dlat = pi/num_srch_bins
598
599        do n=1,num_srch_bins
600          bin_lats(1,n) = (n-1)*dlat - pih
601          bin_lats(2,n) =     n*dlat - pih
602          bin_lons(1,n) = zero
603          bin_lons(2,n) = pi2
604          bin_addr1(1,n) = grid1_size + 1
605          bin_addr1(2,n) = 0
606          bin_addr2(1,n) = grid2_size + 1
607          bin_addr2(2,n) = 0
608        end do
609
610        do nele=1,grid1_size
611          do n=1,num_srch_bins
612            if (grid1_bound_box(1,nele) <= bin_lats(2,n) .and.
613     &          grid1_bound_box(2,nele) >= bin_lats(1,n)) then
614              bin_addr1(1,n) = min(nele,bin_addr1(1,n))
615              bin_addr1(2,n) = max(nele,bin_addr1(2,n))
616            endif
617          end do
618        end do
619
620        do nele=1,grid2_size
621          do n=1,num_srch_bins
622            if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and.
623     &          grid2_bound_box(2,nele) >= bin_lats(1,n)) then
624              bin_addr2(1,n) = min(nele,bin_addr2(1,n))
625              bin_addr2(2,n) = max(nele,bin_addr2(2,n))
626            endif
627          end do
628        end do
629
630      case ('LATLON')
631        write(nulou,*) 'Using lat/lon boxes to restrict search.'
632
633        dlat = pi /num_srch_bins
634        dlon = pi2/num_srch_bins
635
636        allocate(bin_addr1(2,num_srch_bins*num_srch_bins))
637        allocate(bin_addr2(2,num_srch_bins*num_srch_bins))
638        allocate(bin_lats (2,num_srch_bins*num_srch_bins))
639        allocate(bin_lons (2,num_srch_bins*num_srch_bins))
640
641        n = 0
642        do j=1,num_srch_bins
643          do i=1,num_srch_bins
644            n = n + 1
645           
646            bin_lats(1,n) = (j-1)*dlat - pih
647            bin_lats(2,n) =     j*dlat - pih
648            bin_lons(1,n) = (i-1)*dlon
649            bin_lons(2,n) =     i*dlon
650            bin_addr1(1,n) = grid1_size + 1
651            bin_addr1(2,n) = 0
652            bin_addr2(1,n) = grid2_size + 1
653            bin_addr2(2,n) = 0
654          end do
655        end do
656       
657        num_srch_bins = num_srch_bins**2
658
659        do nele=1,grid1_size
660          do n=1,num_srch_bins
661            if (grid1_bound_box(1,nele) <= bin_lats(2,n) .and.
662     &          grid1_bound_box(2,nele) >= bin_lats(1,n) .and.
663     &          grid1_bound_box(3,nele) <= bin_lons(2,n) .and.
664     &          grid1_bound_box(4,nele) >= bin_lons(1,n)) then
665                bin_addr1(1,n) = min(nele,bin_addr1(1,n))
666                bin_addr1(2,n) = max(nele,bin_addr1(2,n))
667            endif
668          end do
669        end do
670
671        do nele=1,grid2_size
672          do n=1,num_srch_bins
673            if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and.
674     &          grid2_bound_box(2,nele) >= bin_lats(1,n) .and.
675     &          grid2_bound_box(3,nele) <= bin_lons(2,n) .and.
676     &          grid2_bound_box(4,nele) >= bin_lons(1,n)) then
677                bin_addr2(1,n) = min(nele,bin_addr2(1,n))
678                bin_addr2(2,n) = max(nele,bin_addr2(2,n))
679            endif
680          end do
681        end do
682
683      case ('REDUCED')
684        write(nulou,*) 
685     &  '|  Using reduced bins to restrict search. Reduced grids with'
686        write(nulou,*) 
687     &  '| a maximum of 500*NBRBINS latitude circles can be treated'
688
689        allocate(bin_addr2(2,num_srch_bins))
690        allocate(bin_lats (2,num_srch_bins))
691        allocate(bin_lons (2,num_srch_bins))
692        allocate(bin_addr1_r(4,500*num_srch_bins))
693        allocate(bin_lats_r (2,500*num_srch_bins))
694        allocate(bin_lons_r (2,500*num_srch_bins))
695
696        dlat = pi/num_srch_bins
697
698        do n=1,num_srch_bins
699          bin_lats(1,n) = (n-1)*dlat - pih
700          bin_lats(2,n) =     n*dlat - pih
701          bin_lons(1,n) = zero
702          bin_lons(2,n) = pi2
703          bin_addr2(1,n) = grid2_size + 1
704          bin_addr2(2,n) = 0
705        end do
706
707        do nele=1,grid2_size
708          do n=1,num_srch_bins
709            if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and.
710     &          grid2_bound_box(2,nele) >= bin_lats(1,n)) then
711              bin_addr2(1,n) = min(nele,bin_addr2(1,n))
712              bin_addr2(2,n) = max(nele,bin_addr2(2,n))
713            endif
714          end do
715        end DO
716
717
718        bin_addr1_r(1,1) = 0
719        bin_lats_r(1,1) = grid1_center_lat(1)
720        num_srch_red = 1
721        do nele=1,grid1_size-1
722          if (grid1_center_lat(nele+1) == 
723     &        grid1_center_lat(nele)) THEN
724              bin_addr1_r(2,num_srch_red) = nele+1
725          else IF (grid1_center_lat(nele+1) /= 
726     &        grid1_center_lat(nele)) THEN
727              bin_addr1_r(1,num_srch_red+1) = nele+1
728              bin_lats_r(2,num_srch_red) = grid1_center_lat(nele+1)
729              bin_lats_r(1,num_srch_red+1) = 
730     &            bin_lats_r(2,num_srch_red)
731              num_srch_red = num_srch_red+1
732          ENDIF
733        end DO
734       
735        DO nele = 1,num_srch_red-1
736          bin_addr1_r(3,nele)=bin_addr1_r(1,nele+1)
737          bin_addr1_r(4,nele)=bin_addr1_r(2,nele+1)
738        enddo
739        bin_addr1_r(3,num_srch_red) = bin_addr1_r(1,num_srch_red-1)
740        bin_addr1_r(4,num_srch_red) = bin_addr1_r(2,num_srch_red-1)
741
742      case default
743        stop 'unknown search restriction method'
744      end select
745!
746      IF (nlogprt .GE. 2) THEN
747         WRITE (UNIT = nulou,FMT = *)' '
748         WRITE (UNIT = nulou,FMT = *)'Leaving routine grid_init'
749         WRITE (UNIT = nulou,FMT = *)' '
750         CALL OASIS_FLUSH_SCRIP(nulou)
751      ENDIF
752!
753!-----------------------------------------------------------------------
754
755      end subroutine grid_init
756
757!***********************************************************************
758
759      subroutine free_grids
760
761!-----------------------------------------------------------------------
762!
763!     subroutine to deallocate allocated arrays
764!
765!-----------------------------------------------------------------------
766!
767      IF (nlogprt .GE. 2) THEN
768         WRITE (UNIT = nulou,FMT = *)' '
769         WRITE (UNIT = nulou,FMT = *)'Entering routine free_grid'
770         WRITE (UNIT = nulou,FMT = *)' '
771         CALL OASIS_FLUSH_SCRIP(nulou)
772      ENDIF
773!
774      deallocate(grid1_mask, grid2_mask, 
775     &    grid1_center_lat, grid1_center_lon,
776     &    grid2_center_lat, grid2_center_lon,
777     &    grid1_area, grid2_area,
778     &    grid1_frac, grid2_frac,
779     &    grid1_dims, grid2_dims)
780
781      IF (restrict_TYPE == 'REDUCED') then
782          deallocate( grid1_bound_box, grid2_bound_box,
783     &        bin_addr1_r, bin_addr2,
784     &        bin_lons, bin_lats,
785     &        bin_lats_r, bin_lons_r)
786      else
787          deallocate( grid1_bound_box, grid2_bound_box,
788     &        bin_addr1, bin_addr2,
789     &        bin_lats, bin_lons)
790      endif
791
792      if (.not. luse_grid_centers) then
793        deallocate(grid1_corner_lat, grid1_corner_lon,
794     &             grid2_corner_lat, grid2_corner_lon)
795      ENDIF
796!
797      IF (nlogprt .GE. 2) THEN
798         WRITE (UNIT = nulou,FMT = *)' '
799         WRITE (UNIT = nulou,FMT = *)'Leaving routine free_grid'
800         WRITE (UNIT = nulou,FMT = *)' '
801         CALL OASIS_FLUSH_SCRIP(nulou)
802      ENDIF
803!-----------------------------------------------------------------------
804
805      end subroutine free_grids
806
807!***********************************************************************
808
809      end module grids
810
811!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
812
Note: See TracBrowser for help on using the repository browser.