New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
remap_conserv.f90 in branches/UKMO/r8395_mix-lyr_diag/NEMOGCM/TOOLS/WEIGHTS/src – NEMO

source: branches/UKMO/r8395_mix-lyr_diag/NEMOGCM/TOOLS/WEIGHTS/src/remap_conserv.f90 @ 11290

Last change on this file since 11290 was 11290, checked in by jcastill, 5 years ago

Remove svn keywords

File size: 73.3 KB
Line 
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!
3!     this module contains necessary routines for computing addresses
4!     and weights for a conservative interpolation  between any two
5!     grids on a sphere.  the weights are computed by performing line
6!     integrals around all overlap regions of the two grids.  see
7!     Dukowicz and Kodis, SIAM J. Sci. Stat. Comput. 8, 305 (1987) and
8!     Jones, P.W. Monthly Weather Review (submitted).
9!
10!-----------------------------------------------------------------------
11!
12!     CVS:$Id$
13!
14!     Copyright (c) 1997, 1998 the Regents of the University of
15!       California.
16!
17!     This software and ancillary information (herein called software)
18!     called SCRIP is made available under the terms described here. 
19!     The software has been approved for release with associated
20!     LA-CC Number 98-45.
21!
22!     Unless otherwise indicated, this software has been authored
23!     by an employee or employees of the University of California,
24!     operator of the Los Alamos National Laboratory under Contract
25!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
26!     Government has rights to use, reproduce, and distribute this
27!     software.  The public may copy and use this software without
28!     charge, provided that this Notice and any statement of authorship
29!     are reproduced on all copies.  Neither the Government nor the
30!     University makes any warranty, express or implied, or assumes
31!     any liability or responsibility for the use of this software.
32!
33!     If software is modified to produce derivative works, such modified
34!     software should be clearly marked, so as not to confuse it with
35!     the version available from Los Alamos National Laboratory.
36!
37!***********************************************************************
38
39      module remap_conservative
40
41!-----------------------------------------------------------------------
42
43      use kinds_mod    ! defines common data types
44      use constants    ! defines common constants
45      use timers       ! module for timing
46      use grids        ! module containing grid information
47      use remap_vars   ! module containing remap information
48
49      implicit none
50
51!-----------------------------------------------------------------------
52!
53!     module variables
54!
55!-----------------------------------------------------------------------
56
57      integer (kind=int_kind), save ::  &
58              num_srch_cells ! num cells in restricted search arrays
59
60      integer (kind=int_kind), dimension(:), allocatable, save ::  &
61              srch_add       ! global address of cells in srch arrays
62
63      real (kind=dbl_kind), parameter ::  &
64           north_thresh = 1.45_dbl_kind,  & ! threshold for coord transf.
65           south_thresh =-2.00_dbl_kind  ! threshold for coord transf.
66
67      real (kind=dbl_kind), dimension(:,:), allocatable, save :: &
68           srch_corner_lat,   & ! lat of each corner of srch cells
69           srch_corner_lon   ! lon of each corner of srch cells
70
71!***********************************************************************
72
73      contains
74
75!***********************************************************************
76
77      subroutine remap_conserv
78
79!-----------------------------------------------------------------------
80!
81!     this routine traces the perimeters of every grid cell on each
82!     grid checking for intersections with the other grid and computing
83!     line integrals for each subsegment.
84!
85!-----------------------------------------------------------------------
86
87!-----------------------------------------------------------------------
88!
89!     local variables
90!
91!-----------------------------------------------------------------------
92
93      integer (kind=int_kind), parameter ::  &
94              max_subseg = 10000 ! max number of subsegments per segment
95                                 ! to prevent infinite loop
96
97      integer (kind=int_kind) ::  &
98              grid1_add,   & ! current linear address for grid1 cell
99              grid2_add,   & ! current linear address for grid2 cell
100              min_add,     & ! addresses for restricting search of
101              max_add,     & !   destination grid
102              n, nwgt,     & ! generic counters
103              corner,      & ! corner of cell that segment starts from
104              next_corn,   & ! corner of cell that segment ends on
105              num_subseg  ! number of subsegments
106
107      logical (kind=log_kind) ::  &
108              lcoinc,   & ! flag for coincident segments
109              lrevers,  & ! flag for reversing direction of segment
110              lbegin   ! flag for first integration of a segment
111
112      logical (kind=log_kind), dimension(:), allocatable :: &
113              srch_mask   ! mask for restricting searches
114
115      real (kind=dbl_kind) :: &
116           intrsct_lat, intrsct_lon,        & ! lat/lon of next intersect
117           beglat, endlat, beglon, endlon,  & ! endpoints of current seg.
118           norm_factor                     ! factor for normalizing wts
119
120      real (kind=dbl_kind), dimension(:), allocatable :: &
121             grid2_centroid_lat, grid2_centroid_lon,  & ! centroid coords
122             grid1_centroid_lat, grid1_centroid_lon  ! on each grid
123
124      real (kind=dbl_kind), dimension(2) :: begseg ! begin lat/lon for
125                                                   ! full segment
126
127      real (kind=dbl_kind), dimension(6) :: weights ! local wgt array
128
129!-----------------------------------------------------------------------
130!
131!     initialize centroid arrays
132!
133!-----------------------------------------------------------------------
134
135      allocate( grid1_centroid_lat(grid1_size), &
136                grid1_centroid_lon(grid1_size), &
137                grid2_centroid_lat(grid2_size), &
138                grid2_centroid_lon(grid2_size))
139
140      grid1_centroid_lat = zero
141      grid1_centroid_lon = zero
142      grid2_centroid_lat = zero
143      grid2_centroid_lon = zero
144
145!-----------------------------------------------------------------------
146!
147!     integrate around each cell on grid1
148!
149!-----------------------------------------------------------------------
150
151      allocate(srch_mask(grid2_size))
152
153      print *,'grid1 sweep '
154      do grid1_add = 1,grid1_size
155
156        !***
157        !*** restrict searches first using search bins
158        !***
159
160        call timer_start(1)
161        min_add = grid2_size
162        max_add = 1
163        do n=1,num_srch_bins
164          if (grid1_add >= bin_addr1(1,n) .and. &
165              grid1_add <= bin_addr1(2,n)) then
166            min_add = min(min_add, bin_addr2(1,n))
167            max_add = max(max_add, bin_addr2(2,n))
168          endif
169        end do
170
171        !***
172        !*** further restrict searches using bounding boxes
173        !***
174
175        num_srch_cells = 0
176        do grid2_add = min_add,max_add
177          srch_mask(grid2_add) = (grid2_bound_box(1,grid2_add) <=  &
178                                  grid1_bound_box(2,grid1_add)) .and. &
179                                 (grid2_bound_box(2,grid2_add) >=  &
180                                  grid1_bound_box(1,grid1_add)) .and. &
181                                 (grid2_bound_box(3,grid2_add) <=  &
182                                  grid1_bound_box(4,grid1_add)) .and. &
183                                 (grid2_bound_box(4,grid2_add) >=  &
184                                  grid1_bound_box(3,grid1_add))
185
186          if (srch_mask(grid2_add)) num_srch_cells = num_srch_cells+1
187        end do
188
189        !***
190        !*** create search arrays
191        !***
192
193        allocate(srch_add(num_srch_cells), &
194                 srch_corner_lat(grid2_corners,num_srch_cells), &
195                 srch_corner_lon(grid2_corners,num_srch_cells))
196
197        n = 0
198        gather1: do grid2_add = min_add,max_add
199          if (srch_mask(grid2_add)) then
200            n = n+1
201            srch_add(n) = grid2_add
202            srch_corner_lat(:,n) = grid2_corner_lat(:,grid2_add)
203            srch_corner_lon(:,n) = grid2_corner_lon(:,grid2_add)
204          endif
205        end do gather1
206        call timer_stop(1)
207
208        !***
209        !*** integrate around this cell
210        !***
211
212        do corner = 1,grid1_corners
213          next_corn = mod(corner,grid1_corners) + 1
214
215          !***
216          !*** define endpoints of the current segment
217          !***
218
219          beglat = grid1_corner_lat(corner,grid1_add)
220          beglon = grid1_corner_lon(corner,grid1_add)
221          endlat = grid1_corner_lat(next_corn,grid1_add)
222          endlon = grid1_corner_lon(next_corn,grid1_add)
223          lrevers = .false.
224
225          !***
226          !*** to ensure exact path taken during both
227          !*** sweeps, always integrate segments in the same
228          !*** direction (SW to NE).
229          !***
230
231          if ((endlat < beglat) .or. &
232              (endlat == beglat .and. endlon < beglon)) then
233            beglat = grid1_corner_lat(next_corn,grid1_add)
234            beglon = grid1_corner_lon(next_corn,grid1_add)
235            endlat = grid1_corner_lat(corner,grid1_add)
236            endlon = grid1_corner_lon(corner,grid1_add)
237            lrevers = .true.
238          endif
239
240          begseg(1) = beglat
241          begseg(2) = beglon
242          lbegin = .true.
243          num_subseg = 0
244
245          !***
246          !*** if this is a constant-longitude segment, skip the rest
247          !*** since the line integral contribution will be zero.
248          !***
249
250          if (endlon /= beglon) then
251
252          !***
253          !*** integrate along this segment, detecting intersections
254          !*** and computing the line integral for each sub-segment
255          !***
256
257          do while (beglat /= endlat .or. beglon /= endlon)
258
259            !***
260            !*** prevent infinite loops if integration gets stuck
261            !*** near cell or threshold boundary
262            !***
263
264            num_subseg = num_subseg + 1
265            if (num_subseg > max_subseg) then
266              stop 'integration stalled: num_subseg exceeded limit'
267            endif
268
269            !***
270            !*** find next intersection of this segment with a grid
271            !*** line on grid 2.
272            !***
273
274            call timer_start(2)
275            call intersection(grid2_add,intrsct_lat,intrsct_lon,lcoinc, &
276                              beglat, beglon, endlat, endlon, begseg,  &
277                              lbegin, lrevers)
278            call timer_stop(2)
279            lbegin = .false.
280
281            !***
282            !*** compute line integral for this subsegment.
283            !***
284
285            call timer_start(3)
286            if (grid2_add /= 0) then
287              call line_integral(weights, num_wts, &
288                               beglon, intrsct_lon, beglat, intrsct_lat, &
289                               grid1_center_lat(grid1_add),  &
290                               grid1_center_lon(grid1_add), &
291                               grid2_center_lat(grid2_add),  &
292                               grid2_center_lon(grid2_add))
293            else
294              call line_integral(weights, num_wts, &
295                               beglon, intrsct_lon, beglat, intrsct_lat, &
296                               grid1_center_lat(grid1_add),  &
297                               grid1_center_lon(grid1_add), &
298                               grid1_center_lat(grid1_add),  &
299                               grid1_center_lon(grid1_add))
300            endif
301            call timer_stop(3)
302
303            !***
304            !*** if integrating in reverse order, change
305            !*** sign of weights
306            !***
307
308            if (lrevers) then
309              weights = -weights
310            endif
311
312            !***
313            !*** store the appropriate addresses and weights.
314            !*** also add contributions to cell areas and centroids.
315            !***
316
317            !if (grid1_add == 119247) then
318            !  print *,grid1_add,grid2_add,corner,weights(1)
319            !  print *,grid1_corner_lat(:,grid1_add)
320            !  print *,grid1_corner_lon(:,grid1_add)
321            !  print *,grid2_corner_lat(:,grid2_add)
322            !  print *,grid2_corner_lon(:,grid2_add)
323            !  print *,beglat,beglon,intrsct_lat,intrsct_lon
324            !endif
325
326            if (grid2_add /= 0) then
327              if (grid1_mask(grid1_add)) then
328                call timer_start(4)
329                call store_link_cnsrv(grid1_add, grid2_add, weights)
330                call timer_stop(4)
331                grid1_frac(grid1_add) = grid1_frac(grid1_add) +  &
332                                        weights(1)
333                grid2_frac(grid2_add) = grid2_frac(grid2_add) +  &
334                                        weights(num_wts+1)
335              endif
336
337            endif
338
339            grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1)
340            grid1_centroid_lat(grid1_add) =  &
341            grid1_centroid_lat(grid1_add) + weights(2)
342            grid1_centroid_lon(grid1_add) =  &
343            grid1_centroid_lon(grid1_add) + weights(3)
344
345            !***
346            !*** reset beglat and beglon for next subsegment.
347            !***
348
349            beglat = intrsct_lat
350            beglon = intrsct_lon
351          end do
352
353          endif
354
355          !***
356          !*** end of segment
357          !***
358
359        end do
360
361        !***
362        !*** finished with this cell: deallocate search array and
363        !*** start on next cell
364
365        deallocate(srch_add, srch_corner_lat, srch_corner_lon)
366
367      end do
368
369      deallocate(srch_mask)
370
371!-----------------------------------------------------------------------
372!
373!     integrate around each cell on grid2
374!
375!-----------------------------------------------------------------------
376
377      allocate(srch_mask(grid1_size))
378
379      print *,'grid2 sweep '
380      do grid2_add = 1,grid2_size
381
382        !***
383        !*** restrict searches first using search bins
384        !***
385
386        call timer_start(5)
387        min_add = grid1_size
388        max_add = 1
389        do n=1,num_srch_bins
390          if (grid2_add >= bin_addr2(1,n) .and. &
391              grid2_add <= bin_addr2(2,n)) then
392            min_add = min(min_add, bin_addr1(1,n))
393            max_add = max(max_add, bin_addr1(2,n))
394          endif
395        end do
396
397        !***
398        !*** further restrict searches using bounding boxes
399        !***
400
401        num_srch_cells = 0
402        do grid1_add = min_add, max_add
403          srch_mask(grid1_add) = (grid1_bound_box(1,grid1_add) <=  &
404                                  grid2_bound_box(2,grid2_add)) .and. &
405                                 (grid1_bound_box(2,grid1_add) >=  &
406                                  grid2_bound_box(1,grid2_add)) .and. &
407                                 (grid1_bound_box(3,grid1_add) <=  &
408                                  grid2_bound_box(4,grid2_add)) .and. &
409                                 (grid1_bound_box(4,grid1_add) >=  &
410                                  grid2_bound_box(3,grid2_add))
411
412          if (srch_mask(grid1_add)) num_srch_cells = num_srch_cells+1
413        end do
414
415        allocate(srch_add(num_srch_cells), &
416                 srch_corner_lat(grid1_corners,num_srch_cells), &
417                 srch_corner_lon(grid1_corners,num_srch_cells))
418
419        n = 0
420        gather2: do grid1_add = min_add,max_add
421          if (srch_mask(grid1_add)) then
422            n = n+1
423            srch_add(n) = grid1_add
424            srch_corner_lat(:,n) = grid1_corner_lat(:,grid1_add)
425            srch_corner_lon(:,n) = grid1_corner_lon(:,grid1_add)
426          endif
427        end do gather2
428        call timer_stop(5)
429
430        !***
431        !*** integrate around this cell
432        !***
433
434        do corner = 1,grid2_corners
435          next_corn = mod(corner,grid2_corners) + 1
436
437          beglat = grid2_corner_lat(corner,grid2_add)
438          beglon = grid2_corner_lon(corner,grid2_add)
439          endlat = grid2_corner_lat(next_corn,grid2_add)
440          endlon = grid2_corner_lon(next_corn,grid2_add)
441          lrevers = .false.
442
443          !***
444          !*** to ensure exact path taken during both
445          !*** sweeps, always integrate in the same direction
446          !***
447
448          if ((endlat < beglat) .or. &
449              (endlat == beglat .and. endlon < beglon)) then
450            beglat = grid2_corner_lat(next_corn,grid2_add)
451            beglon = grid2_corner_lon(next_corn,grid2_add)
452            endlat = grid2_corner_lat(corner,grid2_add)
453            endlon = grid2_corner_lon(corner,grid2_add)
454            lrevers = .true.
455          endif
456
457          begseg(1) = beglat
458          begseg(2) = beglon
459          lbegin = .true.
460
461          !***
462          !*** if this is a constant-longitude segment, skip the rest
463          !*** since the line integral contribution will be zero.
464          !***
465
466          if (endlon /= beglon) then
467          num_subseg = 0
468
469          !***
470          !*** integrate along this segment, detecting intersections
471          !*** and computing the line integral for each sub-segment
472          !***
473
474          do while (beglat /= endlat .or. beglon /= endlon)
475
476            !***
477            !*** prevent infinite loops if integration gets stuck
478            !*** near cell or threshold boundary
479            !***
480
481            num_subseg = num_subseg + 1
482            if (num_subseg > max_subseg) then
483              stop 'integration stalled: num_subseg exceeded limit'
484            endif
485
486            !***
487            !*** find next intersection of this segment with a line
488            !*** on grid 2.
489            !***
490
491            call timer_start(6)
492            call intersection(grid1_add,intrsct_lat,intrsct_lon,lcoinc, &
493                              beglat, beglon, endlat, endlon, begseg, &
494                              lbegin, lrevers)
495            call timer_stop(6)
496            lbegin = .false.
497
498            !***
499            !*** compute line integral for this subsegment.
500            !***
501
502            call timer_start(7)
503            if (grid1_add /= 0) then
504              call line_integral(weights, num_wts, &
505                               beglon, intrsct_lon, beglat, intrsct_lat, &
506                               grid1_center_lat(grid1_add),  &
507                               grid1_center_lon(grid1_add), &
508                               grid2_center_lat(grid2_add),  &
509                               grid2_center_lon(grid2_add))
510            else
511              call line_integral(weights, num_wts, &
512                               beglon, intrsct_lon, beglat, intrsct_lat, &
513                               grid2_center_lat(grid2_add),  &
514                               grid2_center_lon(grid2_add), &
515                               grid2_center_lat(grid2_add),  &
516                               grid2_center_lon(grid2_add))
517            endif
518            call timer_stop(7)
519
520            if (lrevers) then
521              weights = -weights
522            endif
523
524            !***
525            !*** store the appropriate addresses and weights.
526            !*** also add contributions to cell areas and centroids.
527            !*** if there is a coincidence, do not store weights
528            !*** because they have been captured in the previous loop.
529            !*** the grid1 mask is the master mask
530            !***
531
532            !if (grid1_add == 119247) then
533            !  print *,grid1_add,grid2_add,corner,weights(1)
534            !  print *,grid1_corner_lat(:,grid1_add)
535            !  print *,grid1_corner_lon(:,grid1_add)
536            !  print *,grid2_corner_lat(:,grid2_add)
537            !  print *,grid2_corner_lon(:,grid2_add)
538            !  print *,beglat,beglon,intrsct_lat,intrsct_lon
539            !endif
540
541            if (.not. lcoinc .and. grid1_add /= 0) then
542              if (grid1_mask(grid1_add)) then
543                call timer_start(8)
544                call store_link_cnsrv(grid1_add, grid2_add, weights)
545                call timer_stop(8)
546                grid1_frac(grid1_add) = grid1_frac(grid1_add) +  &
547                                        weights(1)
548                grid2_frac(grid2_add) = grid2_frac(grid2_add) +  &
549                                        weights(num_wts+1)
550              endif
551
552            endif
553
554            grid2_area(grid2_add) = grid2_area(grid2_add) +  &
555                                            weights(num_wts+1)
556            grid2_centroid_lat(grid2_add) =  &
557            grid2_centroid_lat(grid2_add) + weights(num_wts+2)
558            grid2_centroid_lon(grid2_add) =  &
559            grid2_centroid_lon(grid2_add) + weights(num_wts+3)
560
561            !***
562            !*** reset beglat and beglon for next subsegment.
563            !***
564
565            beglat = intrsct_lat
566            beglon = intrsct_lon
567          end do
568
569          endif
570
571          !***
572          !*** end of segment
573          !***
574
575        end do
576
577        !***
578        !*** finished with this cell: deallocate search array and
579        !*** start on next cell
580
581        deallocate(srch_add, srch_corner_lat, srch_corner_lon)
582
583      end do
584
585      deallocate(srch_mask)
586
587!-----------------------------------------------------------------------
588!
589!     correct for situations where N/S pole not explicitly included in
590!     grid (i.e. as a grid corner point). if pole is missing from only
591!     one grid, need to correct only the area and centroid of that
592!     grid.  if missing from both, do complete weight calculation.
593!
594!-----------------------------------------------------------------------
595
596      !*** North Pole
597      weights(1) =  pi2
598      weights(2) =  pi*pi
599      weights(3) =  zero
600      weights(4) =  pi2
601      weights(5) =  pi*pi
602      weights(6) =  zero
603
604      grid1_add = 0
605      pole_loop1: do n=1,grid1_size
606        if (grid1_area(n) < -three*pih .and. &
607            grid1_center_lat(n) > zero) then
608          grid1_add = n
609          exit pole_loop1
610        endif
611      end do pole_loop1
612
613      grid2_add = 0
614      pole_loop2: do n=1,grid2_size
615        if (grid2_area(n) < -three*pih .and. &
616            grid2_center_lat(n) > zero) then
617          grid2_add = n
618          exit pole_loop2
619        endif
620      end do pole_loop2
621
622      if (grid1_add /=0) then
623        grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1)
624        grid1_centroid_lat(grid1_add) =  &
625        grid1_centroid_lat(grid1_add) + weights(2)
626        grid1_centroid_lon(grid1_add) = &
627        grid1_centroid_lon(grid1_add) + weights(3)
628      endif
629
630      if (grid2_add /=0) then
631        grid2_area(grid2_add) = grid2_area(grid2_add) +  &
632                                        weights(num_wts+1)
633        grid2_centroid_lat(grid2_add) =  &
634        grid2_centroid_lat(grid2_add) + weights(num_wts+2)
635        grid2_centroid_lon(grid2_add) = &
636        grid2_centroid_lon(grid2_add) + weights(num_wts+3)
637      endif
638
639      if (grid1_add /= 0 .and. grid2_add /=0) then
640        call store_link_cnsrv(grid1_add, grid2_add, weights)
641
642        grid1_frac(grid1_add) = grid1_frac(grid1_add) +  &
643                                weights(1)
644        grid2_frac(grid2_add) = grid2_frac(grid2_add) +  &
645                                weights(num_wts+1)
646      endif
647
648      !*** South Pole
649      weights(1) =  pi2
650      weights(2) = -pi*pi
651      weights(3) =  zero
652      weights(4) =  pi2
653      weights(5) = -pi*pi
654      weights(6) =  zero
655
656      grid1_add = 0
657      pole_loop3: do n=1,grid1_size
658        if (grid1_area(n) < -three*pih .and. &
659            grid1_center_lat(n) < zero) then
660          grid1_add = n
661          exit pole_loop3
662        endif
663      end do pole_loop3
664
665      grid2_add = 0
666      pole_loop4: do n=1,grid2_size
667        if (grid2_area(n) < -three*pih .and. &
668            grid2_center_lat(n) < zero) then
669          grid2_add = n
670          exit pole_loop4
671        endif
672      end do pole_loop4
673
674      if (grid1_add /=0) then
675        grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1)
676        grid1_centroid_lat(grid1_add) =  &
677        grid1_centroid_lat(grid1_add) + weights(2)
678        grid1_centroid_lon(grid1_add) = &
679        grid1_centroid_lon(grid1_add) + weights(3)
680      endif
681
682      if (grid2_add /=0) then
683        grid2_area(grid2_add) = grid2_area(grid2_add) +  &
684                                        weights(num_wts+1)
685        grid2_centroid_lat(grid2_add) =  &
686        grid2_centroid_lat(grid2_add) + weights(num_wts+2)
687        grid2_centroid_lon(grid2_add) = &
688        grid2_centroid_lon(grid2_add) + weights(num_wts+3)
689      endif
690
691      if (grid1_add /= 0 .and. grid2_add /=0) then
692        call store_link_cnsrv(grid1_add, grid2_add, weights)
693
694        grid1_frac(grid1_add) = grid1_frac(grid1_add) +  &
695                                weights(1)
696        grid2_frac(grid2_add) = grid2_frac(grid2_add) +  &
697                                weights(num_wts+1)
698      endif
699
700!-----------------------------------------------------------------------
701!
702!     finish centroid computation
703!
704!-----------------------------------------------------------------------
705
706      where (grid1_area /= zero)
707        grid1_centroid_lat = grid1_centroid_lat/grid1_area
708        grid1_centroid_lon = grid1_centroid_lon/grid1_area
709      end where
710
711      where (grid2_area /= zero)
712        grid2_centroid_lat = grid2_centroid_lat/grid2_area
713        grid2_centroid_lon = grid2_centroid_lon/grid2_area
714      end where
715
716!-----------------------------------------------------------------------
717!
718!     include centroids in weights and normalize using destination
719!     area if requested
720!
721!-----------------------------------------------------------------------
722
723      do n=1,num_links_map1
724        grid1_add = grid1_add_map1(n)
725        grid2_add = grid2_add_map1(n)
726        do nwgt=1,num_wts
727          weights(        nwgt) = wts_map1(nwgt,n)
728          if (num_maps > 1) then
729            weights(num_wts+nwgt) = wts_map2(nwgt,n)
730          endif
731        end do
732
733        select case(norm_opt)
734        case (norm_opt_dstarea)
735          if (grid2_area(grid2_add) /= zero) then
736            if (luse_grid2_area) then
737              norm_factor = one/grid2_area_in(grid2_add)
738            else
739              norm_factor = one/grid2_area(grid2_add)
740            endif
741          else
742            norm_factor = zero
743          endif
744        case (norm_opt_frcarea)
745          if (grid2_frac(grid2_add) /= zero) then
746            if (luse_grid2_area) then
747              norm_factor = grid2_area(grid2_add)/ &
748                           (grid2_frac(grid2_add)* &
749                            grid2_area_in(grid2_add))
750            else
751              norm_factor = one/grid2_frac(grid2_add)
752            endif
753          else
754            norm_factor = zero
755          endif
756        case (norm_opt_none)
757          norm_factor = one
758        end select
759
760        wts_map1(1,n) =  weights(1)*norm_factor
761        wts_map1(2,n) = (weights(2) - weights(1)* &
762                                    grid1_centroid_lat(grid1_add))* &
763                                    norm_factor
764        wts_map1(3,n) = (weights(3) - weights(1)* &
765                                    grid1_centroid_lon(grid1_add))* &
766                                    norm_factor
767
768        if (num_maps > 1) then
769          select case(norm_opt)
770          case (norm_opt_dstarea)
771            if (grid1_area(grid1_add) /= zero) then
772              if (luse_grid1_area) then
773                norm_factor = one/grid1_area_in(grid1_add)
774              else
775                norm_factor = one/grid1_area(grid1_add)
776              endif
777            else
778              norm_factor = zero
779            endif
780          case (norm_opt_frcarea)
781            if (grid1_frac(grid1_add) /= zero) then
782              if (luse_grid1_area) then
783                norm_factor = grid1_area(grid1_add)/ &
784                             (grid1_frac(grid1_add)* &
785                              grid1_area_in(grid1_add))
786              else
787                norm_factor = one/grid1_frac(grid1_add)
788              endif
789            else
790              norm_factor = zero
791            endif
792          case (norm_opt_none)
793            norm_factor = one
794          end select
795
796          wts_map2(1,n) =  weights(num_wts+1)*norm_factor
797          wts_map2(2,n) = (weights(num_wts+2) - weights(num_wts+1)* &
798                                      grid2_centroid_lat(grid2_add))* &
799                                      norm_factor
800          wts_map2(3,n) = (weights(num_wts+3) - weights(num_wts+1)* &
801                                      grid2_centroid_lon(grid2_add))* &
802                                      norm_factor
803        endif
804
805      end do
806
807      print *, 'Total number of links = ',num_links_map1
808
809      where (grid1_area /= zero) grid1_frac = grid1_frac/grid1_area
810      where (grid2_area /= zero) grid2_frac = grid2_frac/grid2_area
811
812!-----------------------------------------------------------------------
813!
814!     perform some error checking on final weights
815!
816!-----------------------------------------------------------------------
817
818      grid2_centroid_lat = zero
819      grid2_centroid_lon = zero
820
821      do n=1,grid1_size
822        if (grid1_area(n) < -.01) then
823          print *,'Grid 1 area error: ',n,grid1_area(n)
824        endif
825        if (grid1_centroid_lat(n) < -pih-.01 .or. &
826            grid1_centroid_lat(n) >  pih+.01) then
827          print *,'Grid 1 centroid lat error: ',n,grid1_centroid_lat(n)
828        endif
829        grid1_centroid_lat(n) = zero
830        grid1_centroid_lon(n) = zero
831      end do
832
833      do n=1,grid2_size
834        if (grid2_area(n) < -.01) then
835          print *,'Grid 2 area error: ',n,grid2_area(n)
836        endif
837        if (grid2_centroid_lat(n) < -pih-.01 .or. &
838            grid2_centroid_lat(n) >  pih+.01) then
839          print *,'Grid 2 centroid lat error: ',n,grid2_centroid_lat(n)
840        endif
841        grid2_centroid_lat(n) = zero
842        grid2_centroid_lon(n) = zero
843      end do
844
845      do n=1,num_links_map1
846        grid1_add = grid1_add_map1(n)
847        grid2_add = grid2_add_map1(n)
848       
849        if (wts_map1(1,n) < -.01) then
850          print *,'Map 1 weight < 0 ',grid1_add,grid2_add,wts_map1(1,n)
851        endif
852        if (norm_opt /= norm_opt_none .and. wts_map1(1,n) > 1.01) then
853          print *,'Map 1 weight > 1 ',grid1_add,grid2_add,wts_map1(1,n)
854        endif
855        grid2_centroid_lat(grid2_add) =  &
856        grid2_centroid_lat(grid2_add) + wts_map1(1,n)
857
858        if (num_maps > 1) then
859          if (wts_map2(1,n) < -.01) then
860            print *,'Map 2 weight < 0 ',grid1_add,grid2_add, &
861                                        wts_map2(1,n)
862          endif
863          if (norm_opt /= norm_opt_none .and. wts_map2(1,n) > 1.01) then
864            print *,'Map 2 weight < 0 ',grid1_add,grid2_add, &
865                                        wts_map2(1,n)
866          endif
867          grid1_centroid_lat(grid1_add) =  &
868          grid1_centroid_lat(grid1_add) + wts_map2(1,n)
869        endif
870      end do
871
872      do n=1,grid2_size
873        select case(norm_opt)
874        case (norm_opt_dstarea)
875          norm_factor = grid2_frac(grid2_add)
876        case (norm_opt_frcarea)
877          norm_factor = one
878        case (norm_opt_none)
879          if (luse_grid2_area) then
880            norm_factor = grid2_area_in(grid2_add)
881          else
882            norm_factor = grid2_area(grid2_add)
883          endif
884        end select
885        if (abs(grid2_centroid_lat(grid2_add)-norm_factor) > .01) then
886          print *,'Error: sum of wts for map1 ',grid2_add, &
887                  grid2_centroid_lat(grid2_add),norm_factor
888        endif
889      end do
890
891      if (num_maps > 1) then
892        do n=1,grid1_size
893          select case(norm_opt)
894          case (norm_opt_dstarea)
895            norm_factor = grid1_frac(grid1_add)
896          case (norm_opt_frcarea)
897            norm_factor = one
898          case (norm_opt_none)
899            if (luse_grid1_area) then
900              norm_factor = grid1_area_in(grid1_add)
901            else
902              norm_factor = grid1_area(grid1_add)
903            endif
904          end select
905          if (abs(grid1_centroid_lat(grid1_add)-norm_factor) > .01) then
906            print *,'Error: sum of wts for map2 ',grid1_add, &
907                    grid1_centroid_lat(grid1_add),norm_factor
908          endif
909        end do
910      endif
911!-----------------------------------------------------------------------
912
913      end subroutine remap_conserv
914
915!***********************************************************************
916
917      subroutine intersection(location,intrsct_lat,intrsct_lon,lcoinc, &
918                              beglat, beglon, endlat, endlon, begseg, &
919                              lbegin, lrevers)
920
921!-----------------------------------------------------------------------
922!
923!     this routine finds the next intersection of a destination grid
924!     line with the line segment given by beglon, endlon, etc.
925!     a coincidence flag is returned if the segment is entirely
926!     coincident with an ocean grid line.  the cells in which to search
927!     for an intersection must have already been restricted in the
928!     calling routine.
929!
930!-----------------------------------------------------------------------
931
932!-----------------------------------------------------------------------
933!
934!     intent(in):
935!
936!-----------------------------------------------------------------------
937
938      logical (kind=log_kind), intent(in) :: &
939           lbegin,  & ! flag for first integration along this segment
940           lrevers ! flag whether segment integrated in reverse
941
942      real (kind=dbl_kind), intent(in) ::  &
943           beglat, beglon,   & ! beginning lat/lon endpoints for segment
944           endlat, endlon   ! ending    lat/lon endpoints for segment
945
946      real (kind=dbl_kind), dimension(2), intent(inout) ::  &
947           begseg ! begin lat/lon of full segment
948
949!-----------------------------------------------------------------------
950!
951!     intent(out):
952!
953!-----------------------------------------------------------------------
954
955      integer (kind=int_kind), intent(out) :: &
956              location  ! address in destination array containing this
957                        ! segment
958
959      logical (kind=log_kind), intent(out) :: &
960              lcoinc    ! flag segments which are entirely coincident
961                        ! with a grid line
962
963      real (kind=dbl_kind), intent(out) :: &
964           intrsct_lat, intrsct_lon ! lat/lon coords of next intersect.
965
966!-----------------------------------------------------------------------
967!
968!     local variables
969!
970!-----------------------------------------------------------------------
971
972      integer (kind=int_kind) :: n, next_n, cell, srch_corners, pole_loc
973
974      integer (kind=int_kind), save ::  &
975           last_loc  ! save location when crossing threshold
976
977      logical (kind=log_kind) ::  &
978           loutside  ! flags points outside grid
979
980      logical (kind=log_kind), save ::  &
981           lthresh = .false.  ! flags segments crossing threshold bndy
982
983      real (kind=dbl_kind) :: &
984           lon1, lon2,        & ! local longitude variables for segment
985           lat1, lat2,        & ! local latitude  variables for segment
986           grdlon1, grdlon2,  & ! local longitude variables for grid cell
987           grdlat1, grdlat2,  & ! local latitude  variables for grid cell
988           vec1_lat, vec1_lon,  & ! vectors and cross products used
989           vec2_lat, vec2_lon,  & ! during grid search
990           cross_product,  &
991           eps, offset,         & ! small offset away from intersect
992           s1, s2, determ,      & ! variables used for linear solve to
993           mat1, mat2, mat3, mat4, rhs1, rhs2  ! find intersection
994
995      real (kind=dbl_kind), save :: &
996           intrsct_lat_off, intrsct_lon_off ! lat/lon coords offset
997                                            ! for next search
998
999!-----------------------------------------------------------------------
1000!
1001!     initialize defaults, flags, etc.
1002!
1003!-----------------------------------------------------------------------
1004
1005      location = 0
1006      lcoinc = .false.
1007      intrsct_lat = endlat
1008      intrsct_lon = endlon
1009
1010      if (num_srch_cells == 0) return
1011
1012      if (beglat > north_thresh .or. beglat < south_thresh) then
1013
1014        if (lthresh) location = last_loc
1015        call pole_intersection(location, &
1016                     intrsct_lat,intrsct_lon,lcoinc,lthresh, &
1017                     beglat, beglon, endlat, endlon, begseg, lrevers)
1018        if (lthresh) then
1019          last_loc = location
1020          intrsct_lat_off = intrsct_lat
1021          intrsct_lon_off = intrsct_lon
1022        endif
1023        return
1024
1025      endif
1026
1027      loutside = .false.
1028      if (lbegin) then
1029        lat1 = beglat
1030        lon1 = beglon
1031      else
1032        lat1 = intrsct_lat_off
1033        lon1 = intrsct_lon_off
1034      endif
1035      lat2 = endlat
1036      lon2 = endlon
1037      if ((lon2-lon1) > three*pih) then
1038        lon2 = lon2 - pi2
1039      else if ((lon2-lon1) < -three*pih) then
1040        lon2 = lon2 + pi2
1041      endif
1042      s1 = zero
1043
1044!-----------------------------------------------------------------------
1045!
1046!     search for location of this segment in ocean grid using cross
1047!     product method to determine whether a point is enclosed by a cell
1048!
1049!-----------------------------------------------------------------------
1050
1051      call timer_start(12)
1052      srch_corners = size(srch_corner_lat,DIM=1)
1053      srch_loop: do
1054
1055        !***
1056        !*** if last segment crossed threshold, use that location
1057        !***
1058
1059        if (lthresh) then
1060          do cell=1,num_srch_cells
1061            if (srch_add(cell) == last_loc) then
1062              location = last_loc
1063              eps = tiny
1064              exit srch_loop
1065            endif
1066          end do
1067        endif
1068
1069        !***
1070        !*** otherwise normal search algorithm
1071        !***
1072
1073        cell_loop: do cell=1,num_srch_cells
1074          corner_loop: do n=1,srch_corners
1075            next_n = MOD(n,srch_corners) + 1
1076
1077            !***
1078            !*** here we take the cross product of the vector making
1079            !*** up each cell side with the vector formed by the vertex
1080            !*** and search point.  if all the cross products are
1081            !*** positive, the point is contained in the cell.
1082            !***
1083
1084            vec1_lat = srch_corner_lat(next_n,cell) -  &
1085                       srch_corner_lat(n     ,cell)
1086            vec1_lon = srch_corner_lon(next_n,cell) -  &
1087                       srch_corner_lon(n     ,cell)
1088            vec2_lat = lat1 - srch_corner_lat(n,cell)
1089            vec2_lon = lon1 - srch_corner_lon(n,cell)
1090
1091            !***
1092            !*** if endpoint coincident with vertex, offset
1093            !*** the endpoint
1094            !***
1095
1096            if (vec2_lat == 0 .and. vec2_lon == 0) then
1097              lat1 = lat1 + 1.d-10*(lat2-lat1)
1098              lon1 = lon1 + 1.d-10*(lon2-lon1)
1099              vec2_lat = lat1 - srch_corner_lat(n,cell)
1100              vec2_lon = lon1 - srch_corner_lon(n,cell)
1101            endif
1102
1103            !***
1104            !*** check for 0,2pi crossings
1105            !***
1106
1107            if (vec1_lon >  pi) then
1108              vec1_lon = vec1_lon - pi2
1109            else if (vec1_lon < -pi) then
1110              vec1_lon = vec1_lon + pi2
1111            endif
1112            if (vec2_lon >  pi) then
1113              vec2_lon = vec2_lon - pi2
1114            else if (vec2_lon < -pi) then
1115              vec2_lon = vec2_lon + pi2
1116            endif
1117
1118            cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
1119
1120            !***
1121            !*** if the cross product for a side is zero, the point
1122            !***   lies exactly on the side or the side is degenerate
1123            !***   (zero length).  if degenerate, set the cross
1124            !***   product to a positive number.  otherwise perform
1125            !***   another cross product between the side and the
1126            !***   segment itself.
1127            !*** if this cross product is also zero, the line is
1128            !***   coincident with the cell boundary - perform the
1129            !***   dot product and only choose the cell if the dot
1130            !***   product is positive (parallel vs anti-parallel).
1131            !***
1132
1133            if (cross_product == zero) then
1134              if (vec1_lat /= zero .or. vec1_lon /= zero) then
1135                vec2_lat = lat2 - lat1
1136                vec2_lon = lon2 - lon1
1137
1138                if (vec2_lon >  pi) then
1139                  vec2_lon = vec2_lon - pi2
1140                else if (vec2_lon < -pi) then
1141                  vec2_lon = vec2_lon + pi2
1142                endif
1143
1144                cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
1145              else
1146                cross_product = one
1147              endif
1148
1149              if (cross_product == zero) then
1150                lcoinc = .true.
1151                cross_product = vec1_lon*vec2_lon + vec1_lat*vec2_lat
1152                if (lrevers) cross_product = -cross_product
1153              endif
1154            endif
1155
1156            !***
1157            !*** if cross product is less than zero, this cell
1158            !*** doesn't work
1159            !***
1160
1161            if (cross_product < zero) exit corner_loop
1162
1163          end do corner_loop
1164
1165          !***
1166          !*** if cross products all positive, we found the location
1167          !***
1168
1169          if (n > srch_corners) then
1170            location = srch_add(cell)
1171
1172            !***
1173            !*** if the beginning of this segment was outside the
1174            !*** grid, invert the segment so the intersection found
1175            !*** will be the first intersection with the grid
1176            !***
1177
1178            if (loutside) then
1179              lat2 = beglat
1180              lon2 = beglon
1181              location = 0
1182              eps  = -tiny
1183            else
1184              eps  = tiny
1185            endif
1186
1187            exit srch_loop
1188          endif
1189
1190          !***
1191          !*** otherwise move on to next cell
1192          !***
1193
1194        end do cell_loop
1195
1196        !***
1197        !*** if still no cell found, the point lies outside the grid.
1198        !***   take some baby steps along the segment to see if any
1199        !***   part of the segment lies inside the grid. 
1200        !***
1201
1202        loutside = .true.
1203        s1 = s1 + 0.001_dbl_kind
1204        lat1 = beglat + s1*(endlat - beglat)
1205        lon1 = beglon + s1*(lon2   - beglon)
1206
1207        !***
1208        !*** reached the end of the segment and still outside the grid
1209        !*** return no intersection
1210        !***
1211
1212        if (s1 >= one) return
1213
1214      end do srch_loop
1215      call timer_stop(12)
1216
1217!-----------------------------------------------------------------------
1218!
1219!     now that a cell is found, search for the next intersection.
1220!     loop over sides of the cell to find intersection with side
1221!     must check all sides for coincidences or intersections
1222!
1223!-----------------------------------------------------------------------
1224
1225      call timer_start(13)
1226      intrsct_loop: do n=1,srch_corners
1227        next_n = mod(n,srch_corners) + 1
1228
1229        grdlon1 = srch_corner_lon(n     ,cell)
1230        grdlon2 = srch_corner_lon(next_n,cell)
1231        grdlat1 = srch_corner_lat(n     ,cell)
1232        grdlat2 = srch_corner_lat(next_n,cell)
1233
1234        !***
1235        !*** set up linear system to solve for intersection
1236        !***
1237
1238        mat1 = lat2 - lat1
1239        mat2 = grdlat1 - grdlat2
1240        mat3 = lon2 - lon1
1241        mat4 = grdlon1 - grdlon2
1242        rhs1 = grdlat1 - lat1
1243        rhs2 = grdlon1 - lon1
1244
1245        if (mat3 >  pi) then
1246          mat3 = mat3 - pi2
1247        else if (mat3 < -pi) then
1248          mat3 = mat3 + pi2
1249        endif
1250        if (mat4 >  pi) then
1251          mat4 = mat4 - pi2
1252        else if (mat4 < -pi) then
1253          mat4 = mat4 + pi2
1254        endif
1255        if (rhs2 >  pi) then
1256          rhs2 = rhs2 - pi2
1257        else if (rhs2 < -pi) then
1258          rhs2 = rhs2 + pi2
1259        endif
1260
1261        determ = mat1*mat4 - mat2*mat3
1262
1263        !***
1264        !*** if the determinant is zero, the segments are either
1265        !***   parallel or coincident.  coincidences were detected
1266        !***   above so do nothing.
1267        !*** if the determinant is non-zero, solve for the linear
1268        !***   parameters s for the intersection point on each line
1269        !***   segment.
1270        !*** if 0<s1,s2<1 then the segment intersects with this side.
1271        !***   return the point of intersection (adding a small
1272        !***   number so the intersection is off the grid line).
1273        !***
1274
1275        if (abs(determ) > 1.e-30) then
1276
1277          s1 = (rhs1*mat4 - mat2*rhs2)/determ
1278          s2 = (mat1*rhs2 - rhs1*mat3)/determ
1279
1280          if (s2 >= zero .and. s2 <= one .and. &
1281              s1 >  zero .and. s1 <= one) then
1282
1283            !***
1284            !*** recompute intersection based on full segment
1285            !*** so intersections are consistent for both sweeps
1286            !***
1287
1288            if (.not. loutside) then
1289              mat1 = lat2 - begseg(1)
1290              mat3 = lon2 - begseg(2)
1291              rhs1 = grdlat1 - begseg(1)
1292              rhs2 = grdlon1 - begseg(2)
1293            else
1294              mat1 = begseg(1) - endlat
1295              mat3 = begseg(2) - endlon
1296              rhs1 = grdlat1 - endlat
1297              rhs2 = grdlon1 - endlon
1298            endif
1299
1300            if (mat3 >  pi) then
1301              mat3 = mat3 - pi2
1302            else if (mat3 < -pi) then
1303              mat3 = mat3 + pi2
1304            endif
1305            if (rhs2 >  pi) then
1306              rhs2 = rhs2 - pi2
1307            else if (rhs2 < -pi) then
1308              rhs2 = rhs2 + pi2
1309            endif
1310
1311            determ = mat1*mat4 - mat2*mat3
1312
1313            !***
1314            !*** sometimes due to roundoff, the previous
1315            !*** determinant is non-zero, but the lines
1316            !*** are actually coincident.  if this is the
1317            !*** case, skip the rest.
1318            !***
1319
1320            if (determ /= zero) then
1321              s1 = (rhs1*mat4 - mat2*rhs2)/determ
1322              s2 = (mat1*rhs2 - rhs1*mat3)/determ
1323
1324              offset = s1 + eps/determ
1325              if (offset > one) offset = one
1326
1327              if (.not. loutside) then
1328                intrsct_lat = begseg(1) + mat1*s1
1329                intrsct_lon = begseg(2) + mat3*s1
1330                intrsct_lat_off = begseg(1) + mat1*offset
1331                intrsct_lon_off = begseg(2) + mat3*offset
1332              else
1333                intrsct_lat = endlat + mat1*s1
1334                intrsct_lon = endlon + mat3*s1
1335                intrsct_lat_off = endlat + mat1*offset
1336                intrsct_lon_off = endlon + mat3*offset
1337              endif
1338              exit intrsct_loop
1339            endif
1340
1341          endif
1342        endif
1343
1344        !***
1345        !*** no intersection this side, move on to next side
1346        !***
1347
1348      end do intrsct_loop
1349      call timer_stop(13)
1350
1351!-----------------------------------------------------------------------
1352!
1353!     if the segment crosses a pole threshold, reset the intersection
1354!     to be the threshold latitude.  only check if this was not a
1355!     threshold segment since sometimes coordinate transform can end
1356!     up on other side of threshold again.
1357!
1358!-----------------------------------------------------------------------
1359
1360      if (lthresh) then
1361        if (intrsct_lat < north_thresh .or. intrsct_lat > south_thresh) &
1362            lthresh = .false.
1363      else if (lat1 > zero .and. intrsct_lat > north_thresh) then
1364        intrsct_lat = north_thresh + tiny
1365        intrsct_lat_off = north_thresh + eps*mat1
1366        s1 = (intrsct_lat - begseg(1))/mat1
1367        intrsct_lon     = begseg(2) + s1*mat3
1368        intrsct_lon_off = begseg(2) + (s1+eps)*mat3
1369        last_loc = location
1370        lthresh = .true.
1371      else if (lat1 < zero .and. intrsct_lat < south_thresh) then
1372        intrsct_lat = south_thresh - tiny
1373        intrsct_lat_off = south_thresh + eps*mat1
1374        s1 = (intrsct_lat - begseg(1))/mat1
1375        intrsct_lon     = begseg(2) + s1*mat3
1376        intrsct_lon_off = begseg(2) + (s1+eps)*mat3
1377        last_loc = location
1378        lthresh = .true.
1379      endif
1380
1381!-----------------------------------------------------------------------
1382
1383      end subroutine intersection
1384
1385!***********************************************************************
1386
1387      subroutine pole_intersection(location, &
1388                       intrsct_lat,intrsct_lon,lcoinc,lthresh, &
1389                       beglat, beglon, endlat, endlon, begseg, lrevers)
1390
1391!-----------------------------------------------------------------------
1392!
1393!     this routine is identical to the intersection routine except
1394!     that a coordinate transformation (using a Lambert azimuthal
1395!     equivalent projection) is performed to treat polar cells more
1396!     accurately.
1397!
1398!-----------------------------------------------------------------------
1399
1400!-----------------------------------------------------------------------
1401!
1402!     intent(in):
1403!
1404!-----------------------------------------------------------------------
1405
1406      real (kind=dbl_kind), intent(in) ::  &
1407           beglat, beglon,   & ! beginning lat/lon endpoints for segment
1408           endlat, endlon   ! ending    lat/lon endpoints for segment
1409
1410      real (kind=dbl_kind), dimension(2), intent(inout) ::  &
1411           begseg ! begin lat/lon of full segment
1412
1413      logical (kind=log_kind), intent(in) :: &
1414              lrevers   ! flag true if segment integrated in reverse
1415
1416!-----------------------------------------------------------------------
1417!
1418!     intent(out):
1419!
1420!-----------------------------------------------------------------------
1421
1422      integer (kind=int_kind), intent(inout) :: &
1423              location  ! address in destination array containing this
1424                        ! segment -- also may contain last location on
1425                        ! entry
1426
1427      logical (kind=log_kind), intent(out) :: &
1428              lcoinc    ! flag segment coincident with grid line
1429
1430      logical (kind=log_kind), intent(inout) :: &
1431              lthresh   ! flag segment crossing threshold boundary
1432
1433      real (kind=dbl_kind), intent(out) :: &
1434           intrsct_lat, intrsct_lon ! lat/lon coords of next intersect.
1435
1436!-----------------------------------------------------------------------
1437!
1438!     local variables
1439!
1440!-----------------------------------------------------------------------
1441
1442      integer (kind=int_kind) :: n, next_n, cell, srch_corners, pole_loc
1443
1444      logical (kind=log_kind) :: loutside ! flags points outside grid
1445
1446      real (kind=dbl_kind) :: pi4, rns,  & ! north/south conversion
1447           x1, x2,        & ! local x variables for segment
1448           y1, y2,        & ! local y variables for segment
1449           begx, begy,    & ! beginning x,y variables for segment
1450           endx, endy,    & ! beginning x,y variables for segment
1451           begsegx, begsegy,    & ! beginning x,y variables for segment
1452           grdx1, grdx2,  & ! local x variables for grid cell
1453           grdy1, grdy2,  & ! local y variables for grid cell
1454           vec1_y, vec1_x,  & ! vectors and cross products used
1455           vec2_y, vec2_x,  & ! during grid search
1456           cross_product, eps,  & ! eps=small offset away from intersect
1457           s1, s2, determ,      & ! variables used for linear solve to
1458           mat1, mat2, mat3, mat4, rhs1, rhs2  ! find intersection
1459
1460      real (kind=dbl_kind), dimension(:,:), allocatable :: &
1461           srch_corner_x,   & ! x of each corner of srch cells
1462           srch_corner_y   ! y of each corner of srch cells
1463
1464      !***
1465      !*** save last intersection to avoid roundoff during coord
1466      !*** transformation
1467      !***
1468
1469      logical (kind=log_kind), save :: luse_last = .false.
1470
1471      real (kind=dbl_kind), save ::  &
1472           intrsct_x, intrsct_y  ! x,y for intersection
1473
1474      !***
1475      !*** variables necessary if segment manages to hit pole
1476      !***
1477
1478      integer (kind=int_kind), save ::  &
1479           avoid_pole_count = 0  ! count attempts to avoid pole
1480
1481      real (kind=dbl_kind), save ::  &
1482           avoid_pole_offset = tiny  ! endpoint offset to avoid pole
1483
1484!-----------------------------------------------------------------------
1485!
1486!     initialize defaults, flags, etc.
1487!
1488!-----------------------------------------------------------------------
1489
1490      if (.not. lthresh) location = 0
1491      lcoinc = .false.
1492      intrsct_lat = endlat
1493      intrsct_lon = endlon
1494
1495      loutside = .false.
1496      s1 = zero
1497
1498!-----------------------------------------------------------------------
1499!
1500!     convert coordinates
1501!
1502!-----------------------------------------------------------------------
1503
1504      allocate(srch_corner_x(size(srch_corner_lat,DIM=1), &
1505                             size(srch_corner_lat,DIM=2)), &
1506               srch_corner_y(size(srch_corner_lat,DIM=1), &
1507                             size(srch_corner_lat,DIM=2)))
1508
1509      if (beglat > zero) then
1510        pi4 = quart*pi
1511        rns = one
1512      else
1513        pi4 = -quart*pi
1514        rns = -one
1515      endif
1516
1517      if (luse_last) then
1518        x1 = intrsct_x
1519        y1 = intrsct_y
1520      else
1521        x1 = rns*two*sin(pi4 - half*beglat)*cos(beglon)
1522        y1 =     two*sin(pi4 - half*beglat)*sin(beglon)
1523        luse_last = .true.
1524      endif
1525      x2 = rns*two*sin(pi4 - half*endlat)*cos(endlon)
1526      y2 =     two*sin(pi4 - half*endlat)*sin(endlon)
1527      srch_corner_x = rns*two*sin(pi4 - half*srch_corner_lat)* &
1528                              cos(srch_corner_lon)
1529      srch_corner_y =     two*sin(pi4 - half*srch_corner_lat)* &
1530                              sin(srch_corner_lon)
1531
1532      begx = x1
1533      begy = y1
1534      endx = x2
1535      endy = y2
1536      begsegx = rns*two*sin(pi4 - half*begseg(1))*cos(begseg(2))
1537      begsegy =     two*sin(pi4 - half*begseg(1))*sin(begseg(2))
1538      intrsct_x = endx
1539      intrsct_y = endy
1540
1541!-----------------------------------------------------------------------
1542!
1543!     search for location of this segment in ocean grid using cross
1544!     product method to determine whether a point is enclosed by a cell
1545!
1546!-----------------------------------------------------------------------
1547
1548      call timer_start(12)
1549      srch_corners = size(srch_corner_lat,DIM=1)
1550      srch_loop: do
1551
1552        !***
1553        !*** if last segment crossed threshold, use that location
1554        !***
1555
1556        if (lthresh) then
1557          do cell=1,num_srch_cells
1558            if (srch_add(cell) == location) then
1559              eps = tiny
1560              exit srch_loop
1561            endif
1562          end do
1563        endif
1564
1565        !***
1566        !*** otherwise normal search algorithm
1567        !***
1568
1569        cell_loop: do cell=1,num_srch_cells
1570          corner_loop: do n=1,srch_corners
1571            next_n = MOD(n,srch_corners) + 1
1572
1573            !***
1574            !*** here we take the cross product of the vector making
1575            !*** up each cell side with the vector formed by the vertex
1576            !*** and search point.  if all the cross products are
1577            !*** positive, the point is contained in the cell.
1578            !***
1579
1580            vec1_x = srch_corner_x(next_n,cell) -  &
1581                     srch_corner_x(n     ,cell)
1582            vec1_y = srch_corner_y(next_n,cell) -  &
1583                     srch_corner_y(n     ,cell)
1584            vec2_x = x1 - srch_corner_x(n,cell)
1585            vec2_y = y1 - srch_corner_y(n,cell)
1586
1587            !***
1588            !*** if endpoint coincident with vertex, offset
1589            !*** the endpoint
1590            !***
1591
1592            if (vec2_x == 0 .and. vec2_y == 0) then
1593              x1 = x1 + 1.d-10*(x2-x1)
1594              y1 = y1 + 1.d-10*(y2-y1)
1595              vec2_x = x1 - srch_corner_x(n,cell)
1596              vec2_y = y1 - srch_corner_y(n,cell)
1597            endif
1598
1599            cross_product = vec1_x*vec2_y - vec2_x*vec1_y
1600
1601            !***
1602            !*** if the cross product for a side is zero, the point
1603            !***   lies exactly on the side or the length of a side
1604            !***   is zero.  if the length is zero set det > 0.
1605            !***   otherwise, perform another cross
1606            !***   product between the side and the segment itself.
1607            !*** if this cross product is also zero, the line is
1608            !***   coincident with the cell boundary - perform the
1609            !***   dot product and only choose the cell if the dot
1610            !***   product is positive (parallel vs anti-parallel).
1611            !***
1612
1613            if (cross_product == zero) then
1614              if (vec1_x /= zero .or. vec1_y /= 0) then
1615                vec2_x = x2 - x1
1616                vec2_y = y2 - y1
1617                cross_product = vec1_x*vec2_y - vec2_x*vec1_y
1618              else
1619                cross_product = one
1620              endif
1621
1622              if (cross_product == zero) then
1623                lcoinc = .true.
1624                cross_product = vec1_x*vec2_x + vec1_y*vec2_y
1625                if (lrevers) cross_product = -cross_product
1626              endif
1627            endif
1628
1629            !***
1630            !*** if cross product is less than zero, this cell
1631            !*** doesn't work
1632            !***
1633
1634            if (cross_product < zero) exit corner_loop
1635
1636          end do corner_loop
1637
1638          !***
1639          !*** if cross products all positive, we found the location
1640          !***
1641
1642          if (n > srch_corners) then
1643            location = srch_add(cell)
1644
1645            !***
1646            !*** if the beginning of this segment was outside the
1647            !*** grid, invert the segment so the intersection found
1648            !*** will be the first intersection with the grid
1649            !***
1650
1651            if (loutside) then
1652              x2 = begx
1653              y2 = begy
1654              location = 0
1655              eps  = -tiny
1656            else
1657              eps  = tiny
1658            endif
1659
1660            exit srch_loop
1661          endif
1662
1663          !***
1664          !*** otherwise move on to next cell
1665          !***
1666
1667        end do cell_loop
1668
1669        !***
1670        !*** if no cell found, the point lies outside the grid.
1671        !***   take some baby steps along the segment to see if any
1672        !***   part of the segment lies inside the grid. 
1673        !***
1674
1675        loutside = .true.
1676        s1 = s1 + 0.001_dbl_kind
1677        x1 = begx + s1*(x2 - begx)
1678        y1 = begy + s1*(y2 - begy)
1679
1680        !***
1681        !*** reached the end of the segment and still outside the grid
1682        !*** return no intersection
1683        !***
1684
1685        if (s1 >= one) then
1686          deallocate(srch_corner_x, srch_corner_y)
1687          luse_last = .false.
1688          return
1689        endif
1690
1691      end do srch_loop
1692      call timer_stop(12)
1693
1694!-----------------------------------------------------------------------
1695!
1696!     now that a cell is found, search for the next intersection.
1697!     loop over sides of the cell to find intersection with side
1698!     must check all sides for coincidences or intersections
1699!
1700!-----------------------------------------------------------------------
1701
1702      call timer_start(13)
1703      intrsct_loop: do n=1,srch_corners
1704        next_n = mod(n,srch_corners) + 1
1705
1706        grdy1 = srch_corner_y(n     ,cell)
1707        grdy2 = srch_corner_y(next_n,cell)
1708        grdx1 = srch_corner_x(n     ,cell)
1709        grdx2 = srch_corner_x(next_n,cell)
1710
1711        !***
1712        !*** set up linear system to solve for intersection
1713        !***
1714
1715        mat1 = x2 - x1
1716        mat2 = grdx1 - grdx2
1717        mat3 = y2 - y1
1718        mat4 = grdy1 - grdy2
1719        rhs1 = grdx1 - x1
1720        rhs2 = grdy1 - y1
1721
1722        determ = mat1*mat4 - mat2*mat3
1723
1724        !***
1725        !*** if the determinant is zero, the segments are either
1726        !***   parallel or coincident or one segment has zero length. 
1727        !***   coincidences were detected above so do nothing.
1728        !*** if the determinant is non-zero, solve for the linear
1729        !***   parameters s for the intersection point on each line
1730        !***   segment.
1731        !*** if 0<s1,s2<1 then the segment intersects with this side.
1732        !***   return the point of intersection (adding a small
1733        !***   number so the intersection is off the grid line).
1734        !***
1735
1736        if (abs(determ) > 1.e-30) then
1737
1738          s1 = (rhs1*mat4 - mat2*rhs2)/determ
1739          s2 = (mat1*rhs2 - rhs1*mat3)/determ
1740
1741          if (s2 >= zero .and. s2 <= one .and. &
1742              s1 >  zero .and. s1 <= one) then
1743
1744            !***
1745            !*** recompute intersection using entire segment
1746            !*** for consistency between sweeps
1747            !***
1748
1749            if (.not. loutside) then
1750              mat1 = x2 - begsegx
1751              mat3 = y2 - begsegy
1752              rhs1 = grdx1 - begsegx
1753              rhs2 = grdy1 - begsegy
1754            else
1755              mat1 = x2 - endx
1756              mat3 = y2 - endy
1757              rhs1 = grdx1 - endx
1758              rhs2 = grdy1 - endy
1759            endif
1760
1761            determ = mat1*mat4 - mat2*mat3
1762
1763            !***
1764            !*** sometimes due to roundoff, the previous
1765            !*** determinant is non-zero, but the lines
1766            !*** are actually coincident.  if this is the
1767            !*** case, skip the rest.
1768            !***
1769
1770            if (determ /= zero) then
1771              s1 = (rhs1*mat4 - mat2*rhs2)/determ
1772              s2 = (mat1*rhs2 - rhs1*mat3)/determ
1773
1774              if (.not. loutside) then
1775                intrsct_x = begsegx + s1*mat1
1776                intrsct_y = begsegy + s1*mat3
1777              else
1778                intrsct_x = endx + s1*mat1
1779                intrsct_y = endy + s1*mat3
1780              endif
1781
1782              !***
1783              !*** convert back to lat/lon coordinates
1784              !***
1785
1786              intrsct_lon = rns*atan2(intrsct_y,intrsct_x)
1787              if (intrsct_lon < zero)  &
1788                intrsct_lon = intrsct_lon + pi2
1789
1790              if (abs(intrsct_x) > 1.d-10) then
1791                intrsct_lat = (pi4 -  &
1792                  asin(rns*half*intrsct_x/cos(intrsct_lon)))*two
1793              else if (abs(intrsct_y) > 1.d-10) then
1794                intrsct_lat = (pi4 -  &
1795                  asin(half*intrsct_y/sin(intrsct_lon)))*two
1796              else
1797                intrsct_lat = two*pi4
1798              endif
1799
1800              !***
1801              !*** add offset in transformed space for next pass.
1802              !***
1803
1804              if (s1 - eps/determ < one) then
1805                intrsct_x = intrsct_x - mat1*(eps/determ)
1806                intrsct_y = intrsct_y - mat3*(eps/determ)
1807              else
1808                if (.not. loutside) then
1809                  intrsct_x = endx
1810                  intrsct_y = endy
1811                  intrsct_lat = endlat
1812                  intrsct_lon = endlon
1813                else
1814                  intrsct_x = begsegx
1815                  intrsct_y = begsegy
1816                  intrsct_lat = begseg(1)
1817                  intrsct_lon = begseg(2)
1818                endif
1819              endif
1820
1821              exit intrsct_loop
1822            endif
1823          endif
1824        endif
1825
1826        !***
1827        !*** no intersection this side, move on to next side
1828        !***
1829
1830      end do intrsct_loop
1831      call timer_stop(13)
1832
1833      deallocate(srch_corner_x, srch_corner_y)
1834
1835!-----------------------------------------------------------------------
1836!
1837!     if segment manages to cross over pole, shift the beginning
1838!     endpoint in order to avoid hitting pole directly
1839!     (it is ok for endpoint to be pole point)
1840!
1841!-----------------------------------------------------------------------
1842
1843      if (abs(intrsct_x) < 1.e-10 .and. abs(intrsct_y) < 1.e-10 .and. &
1844          (endx /= zero .and. endy /=0)) then
1845        if (avoid_pole_count > 2) then
1846           avoid_pole_count = 0
1847           avoid_pole_offset = 10.*avoid_pole_offset
1848        endif
1849
1850        cross_product = begsegx*(endy-begsegy) - begsegy*(endx-begsegx)
1851        intrsct_lat = begseg(1)
1852        if (cross_product*intrsct_lat > zero) then
1853          intrsct_lon = beglon    + avoid_pole_offset
1854          begseg(2)   = begseg(2) + avoid_pole_offset
1855        else
1856          intrsct_lon = beglon    - avoid_pole_offset
1857          begseg(2)   = begseg(2) - avoid_pole_offset
1858        endif
1859
1860        avoid_pole_count = avoid_pole_count + 1
1861        luse_last = .false.
1862      else
1863        avoid_pole_count = 0
1864        avoid_pole_offset = tiny
1865      endif
1866
1867!-----------------------------------------------------------------------
1868!
1869!     if the segment crosses a pole threshold, reset the intersection
1870!     to be the threshold latitude and do not reuse x,y intersect
1871!     on next entry.  only check if did not cross threshold last
1872!     time - sometimes the coordinate transformation can place a
1873!     segment on the other side of the threshold again
1874!
1875!-----------------------------------------------------------------------
1876
1877      if (lthresh) then
1878        if (intrsct_lat > north_thresh .or. intrsct_lat < south_thresh) &
1879          lthresh = .false.
1880      else if (beglat > zero .and. intrsct_lat < north_thresh) then
1881        mat4 = endlat - begseg(1)
1882        mat3 = endlon - begseg(2)
1883        if (mat3 >  pi) mat3 = mat3 - pi2
1884        if (mat3 < -pi) mat3 = mat3 + pi2
1885        intrsct_lat = north_thresh - tiny
1886        s1 = (north_thresh - begseg(1))/mat4
1887        intrsct_lon = begseg(2) + s1*mat3
1888        luse_last = .false.
1889        lthresh = .true.
1890      else if (beglat < zero .and. intrsct_lat > south_thresh) then
1891        mat4 = endlat - begseg(1)
1892        mat3 = endlon - begseg(2)
1893        if (mat3 >  pi) mat3 = mat3 - pi2
1894        if (mat3 < -pi) mat3 = mat3 + pi2
1895        intrsct_lat = south_thresh + tiny
1896        s1 = (south_thresh - begseg(1))/mat4
1897        intrsct_lon = begseg(2) + s1*mat3
1898        luse_last = .false.
1899        lthresh = .true.
1900      endif
1901
1902      !***
1903      !*** if reached end of segment, do not use x,y intersect
1904      !*** on next entry
1905      !***
1906
1907      if (intrsct_lat == endlat .and. intrsct_lon == endlon) then
1908        luse_last = .false.
1909      endif
1910
1911!-----------------------------------------------------------------------
1912
1913      end subroutine pole_intersection
1914
1915!***********************************************************************
1916
1917      subroutine line_integral(weights, num_wts,  &
1918                             in_phi1, in_phi2, theta1, theta2, &
1919                             grid1_lat, grid1_lon, grid2_lat, grid2_lon)
1920
1921!-----------------------------------------------------------------------
1922!
1923!     this routine computes the line integral of the flux function
1924!     that results in the interpolation weights.  the line is defined
1925!     by the input lat/lon of the endpoints.
1926!
1927!-----------------------------------------------------------------------
1928
1929!-----------------------------------------------------------------------
1930!
1931!     intent(in):
1932!
1933!-----------------------------------------------------------------------
1934
1935      integer (kind=int_kind), intent(in) :: &
1936              num_wts  ! number of weights to compute
1937
1938      real (kind=dbl_kind), intent(in) ::  &
1939           in_phi1, in_phi2,      & ! longitude endpoints for the segment
1940           theta1, theta2,        & ! latitude  endpoints for the segment
1941           grid1_lat, grid1_lon,  & ! reference coordinates for each
1942           grid2_lat, grid2_lon  ! grid (to ensure correct 0,2pi interv.
1943
1944!-----------------------------------------------------------------------
1945!
1946!     intent(out):
1947!
1948!-----------------------------------------------------------------------
1949
1950      real (kind=dbl_kind), dimension(2*num_wts), intent(out) :: &
1951           weights   ! line integral contribution to weights
1952
1953!-----------------------------------------------------------------------
1954!
1955!     local variables
1956!
1957!-----------------------------------------------------------------------
1958
1959      real (kind=dbl_kind) :: dphi, sinth1, sinth2, costh1, costh2, fac, &
1960                              phi1, phi2, phidiff1, phidiff2, sinint
1961      real (kind=dbl_kind) :: f1, f2, fint
1962
1963!-----------------------------------------------------------------------
1964!
1965!     weights for the general case based on a trapezoidal approx to
1966!     the integrals.
1967!
1968!-----------------------------------------------------------------------
1969
1970      sinth1 = SIN(theta1)
1971      sinth2 = SIN(theta2)
1972      costh1 = COS(theta1)
1973      costh2 = COS(theta2)
1974
1975      dphi = in_phi1 - in_phi2
1976      if (dphi >  pi) then
1977        dphi = dphi - pi2
1978      else if (dphi < -pi) then
1979        dphi = dphi + pi2
1980      endif
1981      dphi = half*dphi
1982
1983!-----------------------------------------------------------------------
1984!
1985!     the first weight is the area overlap integral. the second and
1986!     fourth are second-order latitude gradient weights.
1987!
1988!-----------------------------------------------------------------------
1989
1990      weights(        1) = dphi*(sinth1 + sinth2)
1991      weights(num_wts+1) = dphi*(sinth1 + sinth2)
1992      weights(        2) = dphi*(costh1 + costh2 + (theta1*sinth1 + &
1993                                                    theta2*sinth2))
1994      weights(num_wts+2) = dphi*(costh1 + costh2 + (theta1*sinth1 + &
1995                                                    theta2*sinth2))
1996
1997!-----------------------------------------------------------------------
1998!
1999!     the third and fifth weights are for the second-order phi gradient
2000!     component.  must be careful of longitude range.
2001!
2002!-----------------------------------------------------------------------
2003
2004      f1 = half*(costh1*sinth1 + theta1)
2005      f2 = half*(costh2*sinth2 + theta2)
2006
2007      phi1 = in_phi1 - grid1_lon
2008      if (phi1 >  pi) then
2009        phi1 = phi1 - pi2
2010      else if (phi1 < -pi) then
2011        phi1 = phi1 + pi2
2012      endif
2013
2014      phi2 = in_phi2 - grid1_lon
2015      if (phi2 >  pi) then
2016        phi2 = phi2 - pi2
2017      else if (phi2 < -pi) then
2018        phi2 = phi2 + pi2
2019      endif
2020
2021      if ((phi2-phi1) <  pi .and. (phi2-phi1) > -pi) then
2022        weights(3) = dphi*(phi1*f1 + phi2*f2)
2023      else
2024        if (phi1 > zero) then
2025          fac = pi
2026        else
2027          fac = -pi
2028        endif
2029        fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
2030        weights(3) = half*phi1*(phi1-fac)*f1 - &
2031                     half*phi2*(phi2+fac)*f2 + &
2032                     half*fac*(phi1+phi2)*fint
2033      endif
2034
2035      phi1 = in_phi1 - grid2_lon
2036      if (phi1 >  pi) then
2037        phi1 = phi1 - pi2
2038      else if (phi1 < -pi) then
2039        phi1 = phi1 + pi2
2040      endif
2041
2042      phi2 = in_phi2 - grid2_lon
2043      if (phi2 >  pi) then
2044        phi2 = phi2 - pi2
2045      else if (phi2 < -pi) then
2046        phi2 = phi2 + pi2
2047      endif
2048
2049      if ((phi2-phi1) <  pi .and. (phi2-phi1) > -pi) then
2050        weights(num_wts+3) = dphi*(phi1*f1 + phi2*f2)
2051      else
2052        if (phi1 > zero) then
2053          fac = pi
2054        else
2055          fac = -pi
2056        endif
2057        fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
2058        weights(num_wts+3) = half*phi1*(phi1-fac)*f1 - &
2059                             half*phi2*(phi2+fac)*f2 + &
2060                             half*fac*(phi1+phi2)*fint
2061      endif
2062
2063!-----------------------------------------------------------------------
2064
2065      end subroutine line_integral
2066
2067!***********************************************************************
2068
2069      subroutine store_link_cnsrv(add1, add2, weights)
2070
2071!-----------------------------------------------------------------------
2072!
2073!     this routine stores the address and weight for this link in
2074!     the appropriate address and weight arrays and resizes those
2075!     arrays if necessary.
2076!
2077!-----------------------------------------------------------------------
2078
2079!-----------------------------------------------------------------------
2080!
2081!     input variables
2082!
2083!-----------------------------------------------------------------------
2084
2085      integer (kind=int_kind), intent(in) :: &
2086              add1,   & ! address on grid1
2087              add2   ! address on grid2
2088
2089      real (kind=dbl_kind), dimension(:), intent(in) :: &
2090              weights ! array of remapping weights for this link
2091
2092!-----------------------------------------------------------------------
2093!
2094!     local variables
2095!
2096!-----------------------------------------------------------------------
2097
2098      integer (kind=int_kind) :: nlink, min_link, max_link ! link index
2099
2100      integer (kind=int_kind), dimension(:,:), allocatable, save :: &
2101              link_add1,   & ! min,max link add to restrict search
2102              link_add2   ! min,max link add to restrict search
2103
2104      logical (kind=log_kind), save :: first_call = .true.
2105
2106!-----------------------------------------------------------------------
2107!
2108!     if all weights are zero, do not bother storing the link
2109!
2110!-----------------------------------------------------------------------
2111
2112      if (all(weights == zero)) return
2113
2114!-----------------------------------------------------------------------
2115!
2116!     restrict the range of links to search for existing links
2117!
2118!-----------------------------------------------------------------------
2119
2120      if (first_call) then
2121        allocate(link_add1(2,grid1_size), link_add2(2,grid2_size))
2122        link_add1 = 0
2123        link_add2 = 0
2124        first_call = .false.
2125        min_link = 1
2126        max_link = 0
2127      else
2128        min_link = min(link_add1(1,add1),link_add2(1,add2))
2129        max_link = max(link_add1(2,add1),link_add2(2,add2))
2130        if (min_link == 0) then
2131          min_link = 1
2132          max_link = 0
2133        endif
2134      endif
2135
2136!-----------------------------------------------------------------------
2137!
2138!     if the link already exists, add the weight to the current weight
2139!     arrays
2140!
2141!-----------------------------------------------------------------------
2142
2143      do nlink=min_link,max_link
2144        if (add1 == grid1_add_map1(nlink)) then
2145        if (add2 == grid2_add_map1(nlink)) then
2146
2147          wts_map1(:,nlink) = wts_map1(:,nlink) + weights(1:num_wts)
2148          if (num_maps == 2) then
2149            wts_map2(:,nlink) = wts_map2(:,nlink) +  &
2150                                        weights(num_wts+1:2*num_wts)
2151          endif
2152          return
2153
2154        endif
2155        endif
2156      end do
2157
2158!-----------------------------------------------------------------------
2159!
2160!     if the link does not yet exist, increment number of links and
2161!     check to see if remap arrays need to be increased to accomodate
2162!     the new link.  then store the link.
2163!
2164!-----------------------------------------------------------------------
2165
2166      num_links_map1  = num_links_map1 + 1
2167      if (num_links_map1 > max_links_map1)  &
2168         call resize_remap_vars(1,resize_increment)
2169
2170      grid1_add_map1(num_links_map1) = add1
2171      grid2_add_map1(num_links_map1) = add2
2172      wts_map1    (:,num_links_map1) = weights(1:num_wts)
2173
2174      if (num_maps > 1) then
2175        num_links_map2  = num_links_map2 + 1
2176        if (num_links_map2 > max_links_map2)  &
2177           call resize_remap_vars(2,resize_increment)
2178
2179        grid1_add_map2(num_links_map2) = add1
2180        grid2_add_map2(num_links_map2) = add2
2181        wts_map2    (:,num_links_map2) = weights(num_wts+1:2*num_wts)
2182      endif
2183
2184      if (link_add1(1,add1) == 0) link_add1(1,add1) = num_links_map1
2185      if (link_add2(1,add2) == 0) link_add2(1,add2) = num_links_map1
2186      link_add1(2,add1) = num_links_map1
2187      link_add2(2,add2) = num_links_map1
2188
2189!-----------------------------------------------------------------------
2190
2191      end subroutine store_link_cnsrv
2192
2193!***********************************************************************
2194
2195      end module remap_conservative
2196
2197!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.