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

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

Imported oasis3 (tag ipslcm5a) from cvs server to svn server (igcmg project).

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