WAVEWATCH III  beta 0.0.1
scrip_remap_conservative.f
Go to the documentation of this file.
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: remap_conserv.f,v 1.10 2001/08/21 21:05:13 pwjones Exp $
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 ! This code has been modified from the version available from
38 ! Los Alamos National Laboratory, for the purpose of running it
39 ! within WW3. Primary modifications:
40 ! - renamed many variables to be unique across the code
41 ! - "save" variables moved from subroutine to module so that
42 ! we can "clear" them later.
43 ! - print statements added.
44 ! - phi_or_theta = 2 instead of phi_or_theta = 1 (important!)
45 !
46 !***********************************************************************
47 ! Modifications introduced by M. Dutour (MD) for
48 ! running with WAVEWATCH III ... see below
49 !
50 !
51 ! BE CAREFUL ABOUT EXPLICIT INITIALIZATION OF VARIABLES IN
52 ! MULTI-THREADED VERSION OF THE CODE - INLINE INITIALIZATION OF
53 ! A VARIABLE IN FORTRAN 90/95 MAKES THE VARIABLE IMPLICITLY STATIC.
54 ! OPENMP FORCES _ALL_ FORTRAN IMPLEMENTATIONS TO MAKE THE VARIABLE
55 ! STATIC (OR OF THE TYPE SAVE) IF IT IS INITIALIZED IN THE
56 ! DECLARATION LINE
57 !
58 !
59 
60  module scrip_remap_conservative
61 
62 !-----------------------------------------------------------------------
63 
64  use scrip_kindsmod ! defines common data types
65  use scrip_constants ! defines common constants
66  use scrip_timers ! module for timing
67  use scrip_grids ! module containing grid information
68  use scrip_remap_vars ! module containing remap information
69  use omp_lib
70 
71  implicit none
72 
73  integer (SCRIP_i4) :: nthreads=2 ! Number of parallel threads
74 
75 !............variables that needed to be moved from "local level" to
76 !............ "module level" in order that we can clear them later.
77 !............These are all local variables that had the "save" attribute
78 !............in the standard version of SCRIP
79 
80  integer (SCRIP_i4), save ::
81  & avoid_pole_count = 0 ! count attempts to avoid pole
82 
83  real (scrip_r8), save ::
84  & avoid_pole_offset = tiny ! endpoint offset to avoid pole
85 
86  integer (SCRIP_i4), dimension(:,:), allocatable, save ::
87  & link_add1, ! min,max link add to restrict search
88  & link_add2 ! min,max link add to restrict search
89 
90  logical (SCRIP_logical), save ::
91  & first_call_store_link_cnsrv = .true.
92 
93  logical (SCRIP_logical), save ::
94  & first_call_locate_segstart= .true.
95 
96  integer (SCRIP_i4), save ::
97  & last_cell_locate_segstart=0, ! save the search parameters
98  & last_cell_grid_num_locate_segstart=0, ! if unchanged, reuse
99  ! search lists
100  & last_srch_grid_num_locate_segstart=0
101 
102  integer (SCRIP_i4), save ::
103  & num_srch_cells_locate_segstart=0,
104  & srch_corners_locate_segstart ! number of corners for
105  ! each cell
106 
107  integer (SCRIP_i4), dimension(:), allocatable, save ::
108  & srch_add_locate_segstart ! global address of cells
109  ! in srch arrays
110 
111  real (scrip_r8), dimension(:,:), allocatable, save ::
112  & srch_corner_lat_locate_segstart, ! lat of each corner of
113  ! srch cells
114  & srch_corner_lon_locate_segstart ! lon of each corner of
115  ! srch cells
116 
117  real(scrip_r8), dimension(:), allocatable, save ::
118  & srch_center_lat_locate_segstart,! lat of center of srch cells
119  & srch_center_lon_locate_segstart ! lon of center of srch cells
120 
121  logical (SCRIP_logical), save ::
122  & first_call_locate_point= .true.
123 
124  integer (SCRIP_i4), save ::
125  & last_cell_locate_point=0, ! save the search parameters
126  & last_cell_grid_num_locate_point=0, ! if unchanged, reuse
127  ! search lists
128  & last_srch_grid_num_locate_point=0
129 
130  integer (SCRIP_i4), save ::
131  & num_srch_cell_locate_points=0,
132  & srch_corners_locate_point ! number of corners for each cell
133 
134  integer (SCRIP_i4), dimension(:), allocatable, save ::
135  & srch_add_locate_point ! global address of cells in
136  ! srch arrays
137 
138  real (scrip_r8), dimension(:,:), allocatable, save ::
139  & srch_corner_lat_locate_point, ! lat of each corner of srch
140  ! cells
141  & srch_corner_lon_locate_point ! lon of each corner of srch
142  ! cells
143 
144  real (scrip_r8), dimension(:), allocatable, save ::
145  & srch_center_lat_locate_point, ! lat of center of srch cells
146  & srch_center_lon_locate_point ! lon of center of srch cells
147 
148  integer (SCRIP_i4), save ::
149  & num_srch_cells_loc_get_srch_cells, ! Number of srch cells
150  ! found
151  & srch_corners_loc_get_srch_cells ! Number of corners for
152  ! search cells
153 
154  integer (SCRIP_i4), dimension(:), allocatable, save ::
155  & srch_add_loc_get_srch_cells ! Global addresses of
156  ! search cells
157 
158  real (scrip_r8), dimension(:,:), allocatable, save ::
159  & srch_corner_lat_loc_get_srch_cells,
160  & srch_corner_lon_loc_get_srch_cells
161 
162  real (scrip_r8), dimension(:), allocatable, save ::
163  & srch_center_lat_loc_get_srch_cells,
164  & srch_center_lon_loc_get_srch_cells
165 
166  integer (SCRIP_i4), save ::
167  & last_cell_add_get_srch_cells,
168  & last_cell_grid_num_get_srch_cells,
169  & last_srch_grid_num_get_srch_cells
170 
171  logical (SCRIP_logical), save ::
172  & first_call_get_srch_cells=.true.
173 
174  logical (SCRIP_logical), save ::
175  & first_call_find_adj_cell=.true.
176 
177  logical (SCRIP_logical), private :: is_master
178  ! module's equivalent of "l_master"
179 
180  integer (SCRIP_i4), save ::
181  & last_cell_find_adj_cell,
182  & last_cell_grid_num_find_adj_cell,
183  & num_srch_cells_find_adj_cell,
184  & srch_corners_find_adj_cell
185 
186  integer (SCRIP_i4), dimension(:), allocatable, save ::
187  & srch_add_find_adj_cell
188  real (scrip_r8), dimension(:,:), allocatable, save ::
189  & srch_corner_lat_find_adj_cell, srch_corner_lon_find_adj_cell
190 
191  real (scrip_r8), dimension(:), allocatable, save ::
192  & srch_center_lat_find_adj_cell, srch_center_lon_find_adj_cell
193 
194 C$OMP THREADPRIVATE(last_cell_grid_num_get_srch_cells,
195 C$OMP& last_srch_grid_num_get_srch_cells,
196 C$OMP& first_call_get_srch_cells,
197 C$OMP& last_cell_add_get_srch_cells,
198 C$OMP& num_srch_cells_loc_get_srch_cells,
199 C$OMP& srch_corners_loc_get_srch_cells,
200 C$OMP& srch_add_loc_get_srch_cells,
201 C$OMP& srch_corner_lat_loc_get_srch_cells,
202 C$OMP& srch_corner_lon_loc_get_srch_cells,
203 C$OMP& srch_center_lat_loc_get_srch_cells,
204 C$OMP& srch_center_lon_loc_get_srch_cells)
205 
206 C$OMP THREADPRIVATE(first_call_locate_segstart,
207 C$OMP& last_cell_locate_segstart,
208 C$OMP& last_cell_grid_num_locate_segstart,
209 C$OMP& last_srch_grid_num_locate_segstart,
210 C$OMP& num_srch_cells_locate_segstart,
211 C$OMP& srch_corners_locate_segstart,
212 C$OMP& srch_add_locate_segstart,
213 C$OMP& srch_corner_lat_locate_segstart,
214 C$OMP& srch_corner_lon_locate_segstart,
215 C$OMP& srch_center_lat_locate_segstart,
216 C$OMP& srch_center_lon_locate_segstart)
217 
218 C$OMP THREADPRIVATE(first_call_locate_point,
219 C$OMP& last_cell_locate_point,
220 C$OMP& last_cell_grid_num_locate_point,
221 C$OMP& last_srch_grid_num_locate_point,
222 C$OMP& num_srch_cell_locate_points,
223 C$OMP& srch_add_locate_point,srch_corner_lat_locate_point,
224 C$OMP& srch_corner_lon_locate_point,
225 C$OMP& srch_center_lat_locate_point,
226 C$OMP& srch_center_lon_locate_point)
227 
228 C$OMP THREADPRIVATE(first_call_find_adj_cell,
229 C$OMP& last_cell_find_adj_cell,
230 C$OMP& last_cell_grid_num_find_adj_cell,
231 C$OMP& num_srch_cells_find_adj_cell,
232 C$OMP& srch_corners_find_adj_cell,
233 C$OMP& srch_add_find_adj_cell,
234 C$OMP& srch_corner_lat_find_adj_cell,
235 C$OMP& srch_corner_lon_find_adj_cell,
236 C$OMP& srch_center_lat_find_adj_cell,
237 C$OMP& srch_center_lon_find_adj_cell)
238 
239 !***********************************************************************
240 
241  contains
242 
243 !***********************************************************************
244 
245  subroutine remap_conserv(l_master, l_test)
246 
247 !-----------------------------------------------------------------------
248 !
249 ! this routine traces the perimeters of every grid cell on each
250 ! grid checking for intersections with the other grid and computing
251 ! line integrals for each subsegment.
252 !
253 !-----------------------------------------------------------------------
254 
255  logical(SCRIP_Logical), intent(in) :: l_master ! Am I the master
256  ! processor (do I/O)?
257  logical(SCRIP_Logical), intent(in) :: l_test ! Whether to
258  !include test output
259 
260 !-----------------------------------------------------------------------
261 !
262 ! local variables
263 !
264 !-----------------------------------------------------------------------
265 
266  integer (SCRIP_i4), parameter ::
267  & phi_or_theta = 2 ! integrate w.r.t. phi (1) or theta (2)
268 
269 
270  integer (SCRIP_i4) ::
271  & i, inext, !
272  & n, nwgt,
273  & grid1_add, ! Current linear address for grid1 cell
274  & grid2_add, ! Current linear address for grid2 cell
275  & grid_num, ! Index (1,2) of grid that we are
276  ! processing
277  & opp_grid_num, ! Index of opposite grid (2,1)
278  & maxrd_cell, ! cell with the max. relative difference
279  ! in area
280  & progint, ! Intervals at which progress is to be
281  ! printed
282  & icount ! for counting
283 
284  real (SCRIP_r8) ::
285  & norm_factor ! factor for normalizing wts
286 
287  real (SCRIP_r8), dimension(6) ::
288  & weights ! Weights array
289 
290  real (SCRIP_r8) ::
291  & beglat, beglon,
292  & endlat, endlon,
293  & ave_reldiff, ! Average rel. diff. in areas
294  & max_reldiff, ! Maximum rel. diff in areas
295  & maxrd_area, ! Computed area for cell with max rel
296  ! diff
297  & maxrd_true ! True area for cell with max rel diff
298 
299  real (SCRIP_r8), dimension(:), allocatable ::
300  & reldiff, ! Relative difference in computed
301  ! and true area
302  & ref_area ! Area of cell as computed by direct
303  ! integration around its boundaries
304 
305 ! call OMP_SET_DYNAMIC(.FALSE.)
306 
307 !-----------------------------------------------------------------------
308 !
309 ! integrate around each cell on grid1
310 !
311 !-----------------------------------------------------------------------
312 
313  is_master=l_master ! set module variable using subroutine input
314  ! argument variable.
315  ! Use the former subsequently.
316 
317  if(is_master)print *,'grid1 sweep'
318 
319 !NRL Progress is slow when the other grid (grid 2) is large, so we use
320 !NRL that. Really, it would be a better to do this with a timer...
321  if (grid2_size > 500000) then
322  progint = 1000
323  elseif (grid2_size > 250000) then
324  progint = 2000
325  elseif (grid2_size > 100000) then
326  progint = 5000
327  else
328  progint = 10000
329  endif
330 
331  grid_num = 1
332  opp_grid_num = 2
333 
334  call timer_start(1)
335 
336 C$OMP PARALLEL DEFAULT(SHARED) PRIVATE(grid1_add) NUM_THREADS(nthreads)
337 
338 C$OMP DO SCHEDULE(DYNAMIC)
339 
340  do grid1_add = 1,grid1_size
341 
342  if (mod(grid1_add,progint) .eq. 0 .and. is_master) then
343  print *, grid1_add,' of ',grid1_size,' cells processed ...'
344  endif
345 
346  call cell_integrate(grid1_add, grid_num, phi_or_theta)
347 
348  end do ! do grid1_add=...
349 
350 C$OMP END DO
351 
352 C$OMP END PARALLEL
353 
354 !-----------------------------------------------------------------------
355 !
356 ! integrate around each cell on grid2
357 !
358 !-----------------------------------------------------------------------
359 
360  if(is_master)print *,'grid2 sweep '
361 
362 !NRL Progress is slow when the other grid (grid 1) is large, so we use
363 !NRL that.
364  if (grid1_size > 500000) then
365  progint = 1000
366  elseif (grid1_size > 250000) then
367  progint = 2000
368  elseif (grid1_size > 100000) then
369  progint = 5000
370  else
371  progint = 10000
372  endif
373 
374  grid_num = 2
375  opp_grid_num = 1
376 
377  call timer_start(2)
378 
379 C$OMP PARALLEL DEFAULT(SHARED) PRIVATE(grid2_add) NUM_THREADS(nthreads)
380 
381 C$OMP DO SCHEDULE(DYNAMIC)
382 
383  do grid2_add = 1,grid2_size
384 
385  if (mod(grid2_add,progint) .eq. 0 .and. is_master) then
386  print *, grid2_add,' of ',grid2_size,' cells processed ...'
387  endif
388 
389  call cell_integrate(grid2_add, grid_num, phi_or_theta)
390 
391  end do ! do grid2_add=...
392 
393 C$OMP END DO
394 
395 C$OMP END PARALLEL
396 
397  call timer_stop(2)
398 
399 !-----------------------------------------------------------------------
400 !
401 ! correct for situations where N/S pole not explicitly included in
402 ! grid (i.e. as a grid corner point). if pole is missing from only
403 ! one grid, need to correct only the area and centroid of that
404 ! grid. if missing from both, do complete weight calculation.
405 ! This is necessary only when integrating w.r.t. phi (longitude)
406 !
407 !-----------------------------------------------------------------------
408 
409  if (phi_or_theta .eq. 1) then
410 
411  !*** North Pole
412  weights(1) = pi2
413  weights(2) = pi*pi
414  weights(3) = zero
415  weights(4) = pi2
416  weights(5) = pi*pi
417  weights(6) = zero
418 
419  if (grid1_npole_cell /=0) then
420  grid1_area(grid1_npole_cell) = grid1_area(grid1_npole_cell)
421  & + weights(1)
422  grid1_centroid_lat(grid1_npole_cell) =
423  & grid1_centroid_lat(grid1_npole_cell) + weights(2)
424  grid1_centroid_lon(grid1_npole_cell) =
425  & grid1_centroid_lon(grid1_npole_cell) + weights(3)
426  endif
427 
428  if (grid2_npole_cell /=0) then
429  grid2_area(grid2_npole_cell) = grid2_area(grid2_npole_cell)
430  & + weights(num_wts+1)
431  grid2_centroid_lat(grid2_npole_cell) =
432  & grid2_centroid_lat(grid2_npole_cell) +
433  & weights(num_wts+2)
434  grid2_centroid_lon(grid2_npole_cell) =
435  & grid2_centroid_lon(grid2_npole_cell) +
436  & weights(num_wts+3)
437  endif
438 
439  if (grid1_npole_cell /= 0 .and. grid2_npole_cell /=0) then
440  call store_link_cnsrv(grid1_npole_cell,
441  & grid2_npole_cell, weights)
442 
443  grid1_frac(grid1_npole_cell) = grid1_frac(grid1_npole_cell)
444  & + weights(1)
445  grid2_frac(grid2_npole_cell) = grid2_frac(grid2_npole_cell)
446  & + weights(num_wts+1)
447  endif
448 
449 
450  !*** South Pole
451  weights(1) = pi2
452  weights(2) = -pi*pi
453  weights(3) = zero
454  weights(4) = pi2
455  weights(5) = -pi*pi
456  weights(6) = zero
457 
458  if (grid1_spole_cell /=0) then
459  grid1_area(grid1_spole_cell) = grid1_area(grid1_spole_cell)
460  & + weights(1)
461  grid1_centroid_lat(grid1_spole_cell) =
462  & grid1_centroid_lat(grid1_spole_cell) + weights(2)
463  grid1_centroid_lon(grid1_spole_cell) =
464  & grid1_centroid_lon(grid1_spole_cell) + weights(3)
465  endif
466 
467  if (grid2_spole_cell /=0) then
468  grid2_area(grid2_spole_cell) = grid2_area(grid2_spole_cell)
469  & + weights(num_wts+1)
470  grid2_centroid_lat(grid2_spole_cell) =
471  & grid2_centroid_lat(grid2_spole_cell) +
472  & weights(num_wts+2)
473  grid2_centroid_lon(grid2_spole_cell) =
474  & grid2_centroid_lon(grid2_spole_cell) +
475  & weights(num_wts+3)
476  endif
477 
478  if (grid1_spole_cell /= 0 .and. grid2_spole_cell /=0) then
479  call store_link_cnsrv(grid1_spole_cell,
480  & grid2_spole_cell, weights)
481 
482  grid1_frac(grid1_spole_cell) = grid1_frac(grid1_spole_cell)
483  & + weights(1)
484  grid2_frac(grid2_spole_cell) = grid2_frac(grid2_spole_cell)
485  & + weights(num_wts+1)
486  endif
487  endif
488 
489 
490 
491  if(is_master)print *, 'Grid sweeps completed'
492 
493 
494 !-----------------------------------------------------------------------
495 !
496 ! finish centroid computation
497 !
498 !-----------------------------------------------------------------------
499 
500  call timer_start(3)
501 
502 C$OMP PARALLEL
503 C$OMP WORKSHARE
504  where (grid1_area /= zero)
505  grid1_centroid_lat = grid1_centroid_lat/grid1_area
506  grid1_centroid_lon = grid1_centroid_lon/grid1_area
507  end where
508 C$OMP END WORKSHARE
509 
510 C$OMP WORKSHARE
511  where (grid2_area /= zero)
512  grid2_centroid_lat = grid2_centroid_lat/grid2_area
513  grid2_centroid_lon = grid2_centroid_lon/grid2_area
514  end where
515 C$OMP END WORKSHARE
516 C$OMP END PARALLEL
517 
518 
519 !-----------------------------------------------------------------------
520 !
521 ! include centroids in weights and normalize using destination
522 ! area if requested
523 !
524 !-----------------------------------------------------------------------
525 
526 C$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads)
527 C$OMP& PRIVATE(n,grid1_add,grid2_add,nwgt,weights,norm_factor)
528 
529 C$OMP DO SCHEDULE(DYNAMIC)
530 
531  do n=1,num_links_map1
532  grid1_add = grid1_add_map1(n)
533  grid2_add = grid2_add_map1(n)
534  do nwgt=1,num_wts
535  weights( nwgt) = wts_map1(nwgt,n)
536  if (num_maps > 1) then
537  weights(num_wts+nwgt) = wts_map2(nwgt,n)
538  endif
539  end do
540 
541  select case(norm_opt)
542  case (norm_opt_dstarea)
543  if (grid2_area(grid2_add) /= zero) then
544  if (luse_grid2_area) then
545  norm_factor = one/grid2_area_in(grid2_add)
546  else
547  norm_factor = one/grid2_area(grid2_add)
548  endif
549  else
550  norm_factor = zero
551  endif
552  case (norm_opt_frcarea)
553  if (grid2_frac(grid2_add) /= zero) then
554  if (luse_grid2_area) then
555  norm_factor = grid2_area(grid2_add)/
556  & (grid2_frac(grid2_add)*
557  & grid2_area_in(grid2_add))
558  else
559  norm_factor = one/grid2_frac(grid2_add)
560  endif
561  else
562  norm_factor = zero
563  endif
564  case (norm_opt_none)
565  norm_factor = one
566  end select
567 
568  wts_map1(1,n) = weights(1)*norm_factor
569  wts_map1(2,n) = (weights(2) - weights(1)*
570  & grid1_centroid_lat(grid1_add))*
571  & norm_factor
572  wts_map1(3,n) = (weights(3) - weights(1)*
573  & grid1_centroid_lon(grid1_add))*
574  & norm_factor
575 
576  if (num_maps > 1) then
577  select case(norm_opt)
578  case (norm_opt_dstarea)
579  if (grid1_area(grid1_add) /= zero) then
580  if (luse_grid1_area) then
581  norm_factor = one/grid1_area_in(grid1_add)
582  else
583  norm_factor = one/grid1_area(grid1_add)
584  endif
585  else
586  norm_factor = zero
587  endif
588  case (norm_opt_frcarea)
589  if (grid1_frac(grid1_add) /= zero) then
590  if (luse_grid1_area) then
591  norm_factor = grid1_area(grid1_add)/
592  & (grid1_frac(grid1_add)*
593  & grid1_area_in(grid1_add))
594  else
595  norm_factor = one/grid1_frac(grid1_add)
596  endif
597  else
598  norm_factor = zero
599  endif
600  case (norm_opt_none)
601  norm_factor = one
602  end select
603 
604  wts_map2(1,n) = weights(num_wts+1)*norm_factor
605  wts_map2(2,n) = (weights(num_wts+2) - weights(num_wts+1)*
606  & grid2_centroid_lat(grid2_add))*
607  & norm_factor
608  wts_map2(3,n) = (weights(num_wts+3) - weights(num_wts+1)*
609  & grid2_centroid_lon(grid2_add))*
610  & norm_factor
611  endif
612 
613  end do
614 
615 C$OMP END DO
616 
617 C$OMP END PARALLEL
618 
619  if(is_master)print *, 'Total number of links = ',num_links_map1
620 
621 C$OMP PARALLEL
622 C$OMP WORKSHARE
623  where (grid1_area /= zero) grid1_frac = grid1_frac/grid1_area
624 C$OMP END WORKSHARE
625 C$OMP WORKSHARE
626  where (grid2_area /= zero) grid2_frac = grid2_frac/grid2_area
627 C$OMP END WORKSHARE
628 C$OMP END PARALLEL
629 
630  call timer_stop(3)
631 
632 !-----------------------------------------------------------------------
633 !
634 ! perform some error checking on final weights
635 !
636 !-----------------------------------------------------------------------
637 
638  allocate(ref_area(grid1_size))
639  allocate(reldiff(grid1_size))
640 
641 C$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads)
642 C$OMP& PRIVATE(n, i, inext, beglat, beglon, endlat, endlon, weights)
643 C$OMP DO SCHEDULE(DYNAMIC)
644 
645  do n=1,grid1_size
646  if (grid1_area(n) < -.01 .and. is_master) then
647  print *,'Grid 1 area error: ',n,grid1_area(n)
648  endif
649  if ((grid1_centroid_lat(n) < -pih-.01 .or.
650  & grid1_centroid_lat(n) > pih+.01) .and. is_master) then
651  print *,'Grid 1 centroid lat error: ',n,grid1_centroid_lat(n)
652  endif
653 
654  ref_area(n) = 0.0
655  do i = 1, grid1_corners
656  inext = 1 + mod(i,grid1_corners)
657 
658  beglat = grid1_corner_lat(i,n)
659  beglon = grid1_corner_lon(i,n)
660  endlat = grid1_corner_lat(inext,n)
661  endlon = grid1_corner_lon(inext,n)
662 
663  if ((phi_or_theta .eq. 1 .and. beglon .eq. endlon) .or.
664  & (phi_or_theta .eq. 2 .and. beglat .eq. endlat)) cycle
665 
666  call line_integral(phi_or_theta, weights, num_wts, beglon,
667  & endlon, beglat, endlat, grid1_center_lat(n),
668  & grid1_center_lon(n), grid1_center_lat(n),
669  & grid1_center_lon(n))
670 
671  ref_area(n) = ref_area(n) + weights(1)
672  enddo
673  enddo
674 C$OMP END DO
675 C$OMP END PARALLEL
676 
677 
678 ! Correct for polar cells
679 
680  if (phi_or_theta .eq. 1) then
681 
682  !*** North Pole
683  weights(1) = pi2
684 
685  if (grid1_npole_cell /=0) then
686  ref_area(grid1_npole_cell) = ref_area(grid1_npole_cell)
687  & + weights(1)
688  endif
689 
690  !*** South Pole
691  weights(1) = pi2
692 
693  if (grid1_spole_cell /=0) then
694  ref_area(grid1_spole_cell) = ref_area(grid1_spole_cell)
695  & + weights(1)
696  endif
697 
698  endif
699 
700 
701  ave_reldiff = 0.0
702  max_reldiff = -1.0
703 
704  do n = 1, grid1_size
705  if(ref_area(n).gt.0.0)then ! added May 21 2013
706  reldiff(n) = abs(ref_area(n)-grid1_area(n))/abs(ref_area(n))
707  endif
708  ave_reldiff = ave_reldiff + reldiff(n)
709  if (reldiff(n) > max_reldiff) then
710  max_reldiff = reldiff(n)
711  maxrd_cell = n
712  maxrd_area = grid1_area(n)
713  maxrd_true = ref_area(n)
714  endif
715  end do
716 
717  ave_reldiff = ave_reldiff/grid1_size
718 
719  if(is_master.and.l_test)then
720  print *
721  print *
722  print *,'Grid 1: Ave. rel. diff. in areas: ',
723  & ave_reldiff
724  print *,' rel. diff. = abs(area-refarea)/refarea'
725  print *
726  print *,'Grid 1: Max. rel. diff. in areas: ',
727  & max_reldiff
728  print *, 'Max rel. diff. is in cell ',maxrd_cell
729  print *, 'Computed Area: ', maxrd_area
730  print *, 'Reference Area: ',maxrd_true
731  print *
732  endif
733 
734  deallocate(ref_area, reldiff)
735 
736 
737 
738  allocate(ref_area(grid2_size))
739  allocate(reldiff(grid2_size))
740 
741 C$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads)
742 C$OMP& PRIVATE(n, i, inext, beglat, beglon, endlat, endlon, weights)
743 C$OMP DO SCHEDULE(DYNAMIC)
744 
745  do n=1,grid2_size
746  if (grid2_area(n) < -.01 .and. is_master) then
747  print *,'Grid 2 area error: ',n,grid2_area(n)
748  endif
749  if ((grid2_centroid_lat(n) < -pih-.01 .or.
750  & grid2_centroid_lat(n) > pih+.01) .and. is_master) then
751  print *,'Grid 2 centroid lat error: ',n,grid2_centroid_lat(n)
752  endif
753 
754  ref_area(n) = 0.0
755  do i = 1, grid2_corners
756  inext = 1 + mod(i,grid2_corners)
757 
758  beglat = grid2_corner_lat(i,n)
759  beglon = grid2_corner_lon(i,n)
760  endlat = grid2_corner_lat(inext,n)
761  endlon = grid2_corner_lon(inext,n)
762 
763  if ((phi_or_theta .eq. 1 .and. beglon .eq. endlon) .or.
764  & (phi_or_theta .eq. 2 .and. beglat .eq. endlat)) cycle
765 
766  call line_integral(phi_or_theta, weights, num_wts, beglon,
767  & endlon, beglat, endlat, grid2_center_lat(n),
768  & grid2_center_lon(n), grid2_center_lat(n),
769  & grid2_center_lon(n))
770 
771  ref_area(n) = ref_area(n) + weights(1)
772  enddo
773  enddo
774 C$OMP END DO
775 C$OMP END PARALLEL
776 
777 
778 ! Correct for polar cells
779 
780  if (phi_or_theta .eq. 1) then
781 
782  !*** North Pole
783  weights(1) = pi2
784 
785  if (grid2_npole_cell /=0) then
786  ref_area(grid2_npole_cell) = ref_area(grid2_npole_cell)
787  & + weights(1)
788  endif
789 
790  !*** South Pole
791  weights(1) = pi2
792 
793  if (grid2_spole_cell /=0) then
794  ref_area(grid2_spole_cell) = ref_area(grid2_spole_cell)
795  & + weights(1)
796  endif
797 
798  endif
799 
800 
801  ave_reldiff = 0.0
802  max_reldiff = -1.0
803 
804  do n = 1, grid2_size
805  reldiff(n) = abs(ref_area(n)-grid2_area(n))/abs(ref_area(n))
806  ave_reldiff = ave_reldiff + reldiff(n)
807  if (reldiff(n) > max_reldiff) then
808  max_reldiff = reldiff(n)
809  maxrd_cell = n
810  maxrd_area = grid2_area(n)
811  maxrd_true = ref_area(n)
812  endif
813  end do
814 
815  ave_reldiff = ave_reldiff/grid2_size
816 
817  if(is_master.and.l_test)then
818  print *
819  print *,'Grid 2: Ave. rel. diff. in areas: ',
820  & ave_reldiff
821  print *,' rel. diff. = abs(area-refarea)/refarea'
822  print *
823  print *,'Grid 2: Max. rel. diff. in areas: ',
824  & max_reldiff
825  print *, 'Max rel. diff. is in cell ',maxrd_cell
826  print *, 'Computed Area: ', maxrd_area
827  print *, 'Reference Area: ',maxrd_true
828  print *
829  endif
830 
831  deallocate(ref_area,reldiff)
832 
833  if(is_master.and.l_test)then
834  print *, 'Computed area = Area of cell computed by adding areas'
835  print *, ' of intersection with other cells'
836  print *, 'Reference area = Area of cell by direct integration'
837  print *
838  endif
839 
840  !***
841  !*** In the following code, gridN_centroid_lat is being used to
842  !*** store running tallies of the cell areas - so it is a
843  !*** misnomer used to avoid allocation of a new variable
844  !***
845 
846  grid1_centroid_lat = zero
847  grid2_centroid_lat = zero
848  icount=0
849 
850 C$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads)
851 C$OMP& PRIVATE(n,grid1_add,grid2_add,nwgt,weights)
852 C$OMP DO SCHEDULE(DYNAMIC)
853 
854  do n=1,num_links_map1
855  grid1_add = grid1_add_map1(n)
856  grid2_add = grid2_add_map1(n)
857 
858  do nwgt=1,num_wts
859  weights( nwgt) = wts_map1(nwgt,n)
860  if (num_maps > 1) then
861  weights(num_wts+nwgt) = wts_map2(nwgt,n)
862  endif
863  end do
864 
865 ! count warnings about weights that will be excluded
866  if (grid2_frac(grid2_add).gt.frac_lowest .and.
867  & grid2_frac(grid2_add).lt.frac_highest .and. is_master) then
868  if ( (wts_map1(1,n) < wt_lowest) )then
869  icount=icount+1
870 ! print statements that were here have been moved to another routine...
871  endif
872  if (norm_opt /= norm_opt_none .and. wts_map1(1,n) >
873  & wt_highest)then
874  icount=icount+1
875 ! print statements that were here have been moved to another routine...
876  endif
877  endif
878 C$OMP CRITICAL
879  grid2_centroid_lat(grid2_add) =
880  & grid2_centroid_lat(grid2_add) + wts_map1(1,n)
881 C$OMP END CRITICAL
882 
883  if (num_maps > 1) then
884  if (wts_map2(1,n) < -.01 .and. is_master) then
885  print *,'Map 2 weight < 0 ',grid1_add,grid2_add,
886  & wts_map2(1,n)
887  endif
888  if (norm_opt /= norm_opt_none .and. wts_map2(1,n) > 1.01
889  & .and. is_master) then
890  print *,'Map 2 weight > 1 ',grid1_add,grid2_add,
891  & wts_map2(1,n)
892  endif
893 C$OMP CRITICAL
894  grid1_centroid_lat(grid1_add) =
895  & grid1_centroid_lat(grid1_add) + wts_map2(1,n)
896 C$OMP END CRITICAL
897  endif
898  end do
899 
900 C$OMP END DO
901 C$OMP END PARALLEL
902 
903  if(icount.gt.0.and.is_master)then
904  print *,'We had problems in ',icount,' points.'
905  endif
906 ! stop condition was here...has been moved to another routine...
907 
908  !***
909  !*** If grid1 has masks, links between some cells of grid1 and
910  !*** grid2 do not exist even though they overlap. In such a case,
911  !*** the following code will generate errors even though nothing
912  !*** is wrong (grid1_centroid_lat or grid2_centroid_lat are never
913  !*** updated in the above loop)
914  !***
915 
916  do n=1,grid2_size
917  select case(norm_opt)
918  case (norm_opt_dstarea)
919  norm_factor = grid2_frac(n)
920  case (norm_opt_frcarea)
921  norm_factor = one
922  case (norm_opt_none)
923  if (luse_grid2_area) then
924  norm_factor = grid2_area_in(n)
925  else
926  norm_factor = grid2_area(n)
927  endif
928  end select
929 ! if (abs(grid2_centroid_lat(n)-norm_factor) > .01
930 ! & .and. is_master) then
931 ! print *,'Warning: sum of wts for map1 ',n,
932 ! & grid2_centroid_lat(n),norm_factor
933 ! endif
934 ! write(501,*)n,grid2_centroid_lat(n)
935  end do
936 
937 
938  if (num_maps > 1) then
939  do n=1,grid1_size
940  select case(norm_opt)
941  case (norm_opt_dstarea)
942  norm_factor = grid1_frac(n)
943  case (norm_opt_frcarea)
944  norm_factor = one
945  case (norm_opt_none)
946  if (luse_grid1_area) then
947  norm_factor = grid1_area_in(n)
948  else
949  norm_factor = grid1_area(n)
950  endif
951  end select
952  if (abs(grid1_centroid_lat(n)-norm_factor) > .01
953  & .and. is_master) then
954  print *,'Error: sum of wts for map2 ',n,
955  & grid1_centroid_lat(n),norm_factor
956  endif
957  end do
958  endif
959 !-----------------------------------------------------------------------
960 
961  call timer_stop(4)
962 
963  if(is_master)print *, 'Finished Conservative Remapping'
964 
965  if(l_test)then
966  call timer_print(1)
967  call timer_print(2)
968  call timer_print(3)
969  call timer_print(4)
970  endif
971 
972  end subroutine remap_conserv
973 
974 !***********************************************************************
975 
976 
977 
978 !***********************************************************************
979 
980  subroutine cellblock_integrate(ibegin, iend, grid_num,
981  & phi_or_theta)
982 
983  integer (SCRIP_i4) :: ibegin, iend, grid_num, phi_or_theta
984 
985  integer (SCRIP_i4) :: cell_add
986 
987 
988  do cell_add = ibegin, iend
989 
990  call cell_integrate(cell_add, grid_num, phi_or_theta)
991 
992  enddo
993 
994 
995 
996  end subroutine cellblock_integrate
997 
998 !***********************************************************************
999 
1000 
1001 
1002 !***********************************************************************
1003 
1004  subroutine cell_integrate(cell_add, grid_num, phi_or_theta)
1005 
1006 !-----------------------------------------------------------------------
1007 !
1008 ! Integrate around cell while finding intersecting with opposite
1009 ! grid cells and finding segments of cell boundary lying in cells
1010 ! of opposite grid
1011 !
1012 !-----------------------------------------------------------------------
1013 
1014 !-----------------------------------------------------------------------
1015 !
1016 ! Input variables
1017 !
1018 !-----------------------------------------------------------------------
1019 
1020  integer (SCRIP_i4) ::
1021  & cell_add, ! cell to be processed
1022  & grid_num, ! grid that the cell belongs to
1023  & phi_or_theta ! Integration var :
1024  ! phi (lon) or theta (lat)
1025 
1026 
1027 !-----------------------------------------------------------------------
1028 !
1029 ! local variables
1030 !
1031 !-----------------------------------------------------------------------
1032 
1033  integer (SCRIP_i4), parameter ::
1034  & max_subseg = 500 ! max number of subsegments per segment
1035  ! to prevent infinite loop
1036 
1037 
1038  integer (SCRIP_i4) ::
1039  & i, inext, !
1040  & j, jnext, ! generic counters
1041  & ic, k, ns, !
1042  & n, next_n, !
1043  & nwgt, it, !
1044  & oppcell_add, ! Cell from opposite grid we are
1045  ! intersecting
1046  & opp_grid_num, ! Index of opposite grid (2,1)
1047  & min_add, ! addresses for restricting search of
1048  & max_add, ! destination grid
1049  & corner, ! corner of cell that segment starts
1050  ! from
1051  & next_corn, ! corner of cell that segment ends on
1052  & nseg, ! number of segments to use to represent
1053  ! edges near the pole
1054  & num_subseg, ! number of subsegments
1055  & bedgeid1, !
1056  & bedgeid2, ! ID of edge that a point is on
1057  & bedgeid3, !
1058  & intedge, ! ID of intersected edge
1059  & last_add, ! Address of last cell we were in
1060  & next_add, ! Address of next cell we will go into
1061  & adj_add ! Address of cell adjacent to current
1062  ! one
1063 
1064  logical (SCRIP_logical) ::
1065  & lcoinc, ! Are segments coincident?
1066  & lrevers, ! Are we integrating segment in reverse?
1067  & lboundary1,
1068  & lboundary2, ! Is point is on cell boundary?
1069  & lboundary3,
1070  & last_lboundary, ! Is last point is on cell bdry?
1071  & loutside, ! Is point outside the grid?
1072  & lthresh, ! Has segment crossed threshold?
1073  & srch_success, ! Was search for segment start
1074  ! successful?
1075  & intrsct_success, ! Was intersection of segment with
1076  ! opposite grid successful?
1077  & inpoly, ! Is point is in polygon
1078  & last_endpt_inpoly, ! Was end point of last segment in cell
1079  & last_endpt_onedge, ! Was end point of last segment on edge
1080  ! of cell
1081  & lstuck, ! Is the walk stuck inside a cell
1082  & seg_outside, ! Is segment completely outside the grid
1083  & bndedge, ! Is segment on the boundary of the grid
1084  & search, ! Do we have to search to locate point
1085  ! in grid
1086  & inpolar, ! Are we in the polar region?
1087  & special_cell, ! Is this a special cell
1088  ! (only 1 vtx at pole)
1089  & l_exit_do ! Do we need to escape from infinite
1090  ! loop? (NRL)
1091 
1092  real (SCRIP_r8) ::
1093  & intrsct_lat, ! lat of next intersection point
1094  & intrsct_lon, ! lon of next intersection point
1095  & beglat, beglon, ! start point of current sub seg
1096  & endlat, endlon, ! endpoint of current seg
1097  ! (chg to endseg?)
1098  & endlat1, endlon1, ! endpoint of current subseg
1099  & norm_factor ! factor for normalizing wts
1100 
1101  real (SCRIP_r8), dimension(2) ::
1102  & begseg ! begin lat/lon for full segment
1103 
1104  real (SCRIP_r8), dimension(6) ::
1105  & weights, ! local wgt array
1106  & rev_weights ! Weights for grid1 and grid2 flipped
1107 
1108  real (SCRIP_r8) ::
1109  & vec1_lat, vec1_lon, ! vectors, products
1110  & vec2_lat, vec2_lon, ! used in grid search
1111  & vec1_len, dp,
1112  & midlat, midlon, ! Midpoint of segment
1113  & tmplat, tmplon,
1114  & srchpt_lat, ! Search point (offset from seg. start)
1115  & srchpt_lon,
1116  & offset, delta, ! Offset and offset increase for search
1117  & sinang2, ! Square of sine of angle b/w two
1118  ! segments
1119  & dist2, ! Square of distance b/w two points
1120  & fullseg_len2, ! Square of full segment length
1121  & partseg_len2, ! Square of length of segment integrated
1122  ! so far
1123  & fullseg_dlat, ! Lat diff of full segment endpoints
1124  & fullseg_dlon, ! Lon diff of full segment endpoints
1125  & prevlon,
1126  & nextlon,
1127  & pole_lat,
1128  & cell_center_lat,
1129  & cell_center_lon,
1130  & oppcell_center_lat,
1131  & oppcell_center_lon
1132 
1133  real (SCRIP_r8), dimension(:), allocatable ::
1134  & cell_corner_lat, ! Local copies of cell coordinates
1135  & cell_corner_lon, ! May be augmented for computational
1136  ! reasons
1137  & oppcell_corner_lat,
1138  & oppcell_corner_lon
1139 
1140  integer (SCRIP_i4) ::
1141  & ncorners, ! Number of corners in local copy of cell
1142  & ncorners_opp, ! Number of corners in local copy of oppcell
1143  & nalloc, ! Allocation for the cell_corner_* array
1144  & nalloc_opp ! Allocation for the oppcell_corner_* array
1145 
1146  real (SCRIP_r8) ::
1147  & tmpwt1, tmpwt2
1148 
1149  integer (SCRIP_i4) ::
1150  & ncorners_at_pole,
1151  & previdx,
1152  & nextidx
1153 
1154 
1155  if (grid_num .eq. 1) then
1156 
1157  !***
1158  !*** Set up a local copy of the cell with room to add
1159  !*** degenerate edges
1160  !***
1161 
1162  ncorners = grid1_corners
1163  nalloc = min(ncorners + 2,
1164  & size(grid1_corner_lat(:,1)))
1165  allocate (cell_corner_lat(nalloc),
1166  & cell_corner_lon(nalloc))
1167 
1168  do corner = 1, ncorners
1169  cell_corner_lat(corner) = grid1_corner_lat(corner,cell_add)
1170  cell_corner_lon(corner) = grid1_corner_lon(corner,cell_add)
1171  enddo
1172 
1173  cell_center_lat = grid1_center_lat(cell_add)
1174  cell_center_lon = grid1_center_lon(cell_add)
1175 
1176  special_cell = special_polar_cell1(cell_add)
1177 
1178 
1179  !***
1180  !*** Also, allocate storage for the cell from the opposite grid
1181  !***
1182 
1183  opp_grid_num = 2
1184  ncorners_opp = grid2_corners
1185  nalloc_opp = ncorners_opp+2
1186  allocate (oppcell_corner_lat(nalloc_opp),
1187  & oppcell_corner_lon(nalloc_opp))
1188 
1189  else
1190 
1191  !***
1192  !*** Set up the cell info with room to add degenerate edges
1193  !***
1194 
1195  ncorners = grid2_corners
1196  nalloc = min(ncorners + 2,
1197  & size(grid2_corner_lat(:,1)))
1198  allocate (cell_corner_lat(nalloc),
1199  & cell_corner_lon(nalloc))
1200 
1201  do corner = 1, ncorners
1202  cell_corner_lat(corner) = grid2_corner_lat(corner,cell_add)
1203  cell_corner_lon(corner) = grid2_corner_lon(corner,cell_add)
1204  enddo
1205 
1206  cell_center_lat = grid2_center_lat(cell_add)
1207  cell_center_lon = grid2_center_lon(cell_add)
1208 
1209  special_cell = special_polar_cell2(cell_add)
1210 
1211  !***
1212  !*** Also, allocate storage for the cell from the opposite grid
1213  !***
1214 
1215  opp_grid_num = 1
1216  ncorners_opp = grid1_corners
1217  nalloc_opp = ncorners_opp + 2
1218  allocate (oppcell_corner_lat(nalloc_opp),
1219  & oppcell_corner_lon(nalloc_opp))
1220 
1221  endif
1222 
1223  if (special_cell) then
1224 
1225  !***
1226  !*** Special cell with only one corner at the pole Such cells
1227  !*** can have an artificially extreme distortion of the edges
1228  !*** when mapped to the Lambert projection because of the span
1229  !*** of longitudes on the edges So we will augment such cells
1230  !*** with degenerate edges at the pole so that the edges coming
1231  !*** off the pole will actually have the same longitude values
1232  !*** at both ends
1233  !***
1234  !*** lon_p lon_p+ lon_p lon_p-
1235  !*** pi/2 pi/2 pi/2 pi/2
1236  !*** * *--------*------*
1237  !*** / \ | |
1238  !*** / \ | |
1239  !*** / \ | |
1240  !*** * * * *
1241  !*** lon_p+ lon_p- lon_p+ lon_p-
1242  !*** lat_p+ lat_p- lat_p+ lat_p-
1243  !***
1244 
1245  call modify_polar_cell(ncorners,nalloc,cell_corner_lat,
1246  & cell_corner_lon)
1247 
1248  endif
1249 
1250  !***
1251  !*** Cell info set up - Now process the cell
1252  !***
1253 
1254  do corner = 1, ncorners
1255  next_corn = mod(corner,ncorners) + 1
1256 
1257  !***
1258  !*** define endpoints of the current segment
1259  !***
1260 
1261  beglat = cell_corner_lat(corner)
1262  beglon = cell_corner_lon(corner)
1263  endlat = cell_corner_lat(next_corn)
1264  endlon = cell_corner_lon(next_corn)
1265  lrevers = .false.
1266 
1267  !***
1268  !*** if this is a constant-longitude segment, skip the rest
1269  !*** since the line integral contribution will be zero.
1270  !***
1271 
1272  if ((phi_or_theta == 1 .and. endlon == beglon) .or.
1273  & (phi_or_theta == 2 .and. endlat == beglat)) cycle
1274 
1275  !***
1276  !*** to ensure exact path taken during both
1277  !*** sweeps, always integrate segments in the same
1278  !*** direction (SW to NE).
1279  !***
1280 
1281  if ((endlat < beglat) .or.
1282  & (endlat == beglat .and. endlon < beglon)) then
1283  tmplat = beglat
1284  beglat = endlat
1285  endlat = tmplat
1286  tmplon = beglon
1287  beglon = endlon
1288  endlon = tmplon
1289  lrevers = .not. lrevers
1290  endif
1291 
1292  !*** But if one of the segment ends is in the polar region,
1293  !*** we want to start from that (makes some logic easier)
1294 
1295  if ((beglat < north_thresh .and. endlat > north_thresh) .or.
1296  & (beglat > south_thresh .and. endlat < south_thresh))
1297  & then
1298  tmplat = beglat
1299  beglat = endlat
1300  endlat = tmplat
1301  tmplon = beglon
1302  beglon = endlon
1303  endlon = tmplon
1304  lrevers = .not. lrevers
1305  endif
1306 
1307  begseg(1) = beglat
1308  begseg(2) = beglon
1309 
1310  fullseg_dlat = endlat-beglat
1311  fullseg_dlon = endlon-beglon
1312  if (fullseg_dlon > pi) fullseg_dlon = fullseg_dlon - pi2
1313  if (fullseg_dlon < -pi) fullseg_dlon = fullseg_dlon + pi2
1314  fullseg_len2 = fullseg_dlat*fullseg_dlat +
1315  & fullseg_dlon*fullseg_dlon
1316 
1317  partseg_len2 = 0.0
1318 
1319  !***
1320  !*** Is this an edge on the boundary of the grid or
1321  !*** on the boundary of the active cells
1322  !***
1323 
1324 ! Commented out by MD
1325 ! call find_adj_cell(cell_add, corner, grid_num, adj_add)
1326 ! if (grid_num .eq. 1) then
1327 ! if (adj_add .eq. 0 .or. .not. grid1_mask(adj_add)) then
1328 ! bndedge = .true.
1329 ! else
1330 ! bndedge = .false.
1331 ! endif
1332 ! else
1333 ! if (adj_add .eq. 0 .or. .not. grid2_mask(adj_add)) then
1334 ! bndedge = .true.
1335 ! else
1336 ! bndedge = .false.
1337 ! endif
1338 ! endif
1339 
1340  call find_adj_cell(cell_add, corner, grid_num, adj_add)
1341  bndedge = .false.
1342  if (grid_num .eq. 1) then
1343  if (adj_add .eq. 0) then
1344  bndedge = .true.
1345  else
1346  if (.not. grid1_mask(adj_add)) then
1347  bndedge = .true.
1348  endif
1349  endif
1350  else
1351  if (adj_add .eq. 0) then
1352  bndedge = .true.
1353  else
1354  if (.not. grid2_mask(adj_add)) then
1355  bndedge = .true.
1356  endif
1357  endif
1358  endif
1359 
1360  !***
1361  !*** integrate along this segment, detecting intersections
1362  !*** and computing the line integral for each sub-segment
1363  !***
1364 
1365  if (beglat .gt. north_thresh .or. beglat .lt. south_thresh)
1366  & then
1367  nseg = npseg ! Use multiple subsegments near the pole
1368  inpolar = .true.
1369  else
1370  nseg = 1
1371  inpolar = .false.
1372  endif
1373 
1374 
1375  last_add = 0
1376  last_lboundary = .false.
1377  last_endpt_inpoly = .false.
1378  last_endpt_onedge = .false.
1379  next_add = 0
1380  search = .true.
1381  ns = 1
1382 
1383 ! outer "do while"
1384 
1385  do while (beglat /= endlat .or. beglon /= endlon)
1386 
1387  l_exit_do=.false. !NRL
1388 
1389  if ((ns .eq. nseg) .or. (inpolar .eqv. .false.)) then
1390  !
1391  ! Last subseg or out of the polar region
1392  ! Go directly to end of segment
1393  !
1394  endlat1 = endlat
1395  endlon1 = endlon
1396  else
1397  endlat1 = begseg(1) + ns*(fullseg_dlat)/nseg
1398  endlon1 = begseg(2) + ns*(fullseg_dlon)/nseg
1399  endif
1400 
1401  num_subseg = 0
1402 
1403 ! inner "do while"
1404 
1405  do while (beglat /= endlat1 .or. beglon /= endlon1)
1406 
1407  !***
1408  !*** If we integrated to the end or just past it (due to
1409  !*** numerical errors), we are done with this segment
1410  !***
1411 
1412 !NRL see notes below re: infinite "do while" loop
1413  l_exit_do=.false. !NRL
1414  if (partseg_len2 .ge. fullseg_len2) then
1415  write(*,*).ge.'partseg_len2 fullseg_len2'
1416  write(*,*)'beglat,beglon = ',beglat,beglon
1417  write(*,*)'endlat,endlon = ',endlat,endlon
1418  write(*,*)'endlat1,endlon1 = ',endlat1,endlon1
1419  write(*,*)'exiting inner do while loop'
1420  l_exit_do=.true. !NRL
1421  exit
1422  end if
1423 
1424  !******************************************************
1425  !*** Try to find which cell of the opposite grid this
1426  !*** segment is starting in and where it is exiting this
1427  !*** cell
1428  !******************************************************
1429 
1430  vec1_lat = endlat1-beglat
1431  vec1_lon = endlon1-beglon
1432  if (vec1_lon > pi) vec1_lon = vec1_lon - pi2
1433  if (vec1_lon < -pi) vec1_lon = vec1_lon + pi2
1434  vec1_len = sqrt(vec1_lat*vec1_lat+vec1_lon*vec1_lon)
1435  vec1_lat = vec1_lat/vec1_len
1436  vec1_lon = vec1_lon/vec1_len
1437 
1438  offset = 100.0*tiny
1439  oppcell_add = 0
1440  delta = 10*tiny
1441  intrsct_success = .false.
1442  loutside = .false.
1443  lstuck = .false.
1444  lboundary1 = .false.
1445  lboundary2 = .false.
1446  lboundary3 = .false.
1447 
1448  do while (.not. intrsct_success)
1449 
1450  !*************************************************
1451  !*** Find out which cell the segment starts in
1452  !*************************************************
1453 
1454  srch_success = .false.
1455  if (search) then
1456 
1457  !***
1458  !*** Offset the start point in ever increasing
1459  !*** amounts until we are able to reliably locate
1460  !*** the point in a cell of grid2. Inability to locate
1461  !*** the point causes the offset amount to increase
1462 
1463  it = 0
1464  do while (.not. srch_success)
1465 
1466  srchpt_lat = beglat + offset*vec1_lat
1467  srchpt_lon = beglon + offset*vec1_lon
1468 
1469  call locate_point(srchpt_lat, srchpt_lon,
1470  & cell_add, grid_num, opp_grid_num,
1471  & oppcell_add, lboundary1, bedgeid1)
1472 
1473  if (oppcell_add .eq. 0) then
1474  loutside = .true.
1475 ! lcoinc added by MD
1476  lcoinc = .false.
1477  exit ! exit the search loop
1478  else
1479  if (oppcell_add .ne. last_add .or. lthresh)
1480  & then
1481  srch_success = .true.
1482  else
1483  offset = offset + delta
1484  if (offset .ge. vec1_len) then
1485  exit
1486  endif
1487  if (it .gt. 3) then
1488  delta = 2.0*delta
1489  it = 0
1490  endif
1491  endif
1492  endif
1493 
1494  it = it + 1
1495  enddo ! do while (.not. srch_success)
1496 
1497  else
1498  if (last_endpt_inpoly) then
1499 
1500  !*** We know the grid cell the end of the last
1501  !*** segment (which is the beginning of this
1502  !*** segment)
1503 
1504  oppcell_add = last_add
1505  lboundary1 = last_lboundary
1506 
1507  else if (next_add .ne. 0) then
1508 
1509  !*** We know the edge of the grid2 cell that the
1510  !*** last segment intersected, so we move into
1511  !*** the adjacent cell
1512 
1513  oppcell_add = next_add
1514  lboundary1 = .true.
1515 
1516  endif
1517 
1518  srch_success = .true.
1519 
1520  endif
1521 
1522  !*****************************************************
1523  !*** Find where the segment exits this cell, if at all
1524  !*****************************************************
1525 
1526  if (srch_success) then
1527 
1528  !***
1529  !*** First setup local copy of oppcell with room for
1530  !*** adding degenerate edges
1531  !***
1532 
1533  if (grid_num .eq. 1) then
1534  ncorners_opp = grid2_corners
1535  do i = 1, ncorners_opp
1536  oppcell_corner_lat(i) =
1537  & grid2_corner_lat(i,oppcell_add)
1538  oppcell_corner_lon(i) =
1539  & grid2_corner_lon(i,oppcell_add)
1540  enddo
1541  oppcell_center_lat =
1542  & grid2_center_lat(oppcell_add)
1543  oppcell_center_lon =
1544  & grid2_center_lon(oppcell_add)
1545 
1546  special_cell = special_polar_cell2(oppcell_add)
1547  else
1548  ncorners_opp = grid1_corners
1549  do i = 1, ncorners_opp
1550  oppcell_corner_lat(i) =
1551  & grid1_corner_lat(i,oppcell_add)
1552  oppcell_corner_lon(i) =
1553  & grid1_corner_lon(i,oppcell_add)
1554  enddo
1555  oppcell_center_lat =
1556  & grid1_center_lat(oppcell_add)
1557  oppcell_center_lon =
1558  & grid1_center_lon(oppcell_add)
1559 
1560  special_cell = special_polar_cell1(oppcell_add)
1561  endif
1562 
1563  if (special_cell) then
1564  call modify_polar_cell(ncorners_opp, nalloc_opp,
1565  & oppcell_corner_lat, oppcell_corner_lon)
1566  endif
1567 
1568  !***
1569  !*** First see if the segment end is
1570  !*** in the same cell
1571  !***
1572 
1573  call ptincell(endlat1,endlon1, oppcell_add,
1574  & ncorners_opp,
1575  & oppcell_corner_lat,oppcell_corner_lon,
1576  & oppcell_center_lat,oppcell_center_lon,
1577  & opp_grid_num,inpoly,
1578  & lboundary2,bedgeid2)
1579 
1580  if (inpoly) then
1581  intrsct_lat = endlat1
1582  intrsct_lon = endlon1
1583  intrsct_success = .true.
1584  search = .false.
1585  next_add = 0
1586  last_add = oppcell_add ! for next subseg
1587  last_lboundary = lboundary2
1588  last_endpt_inpoly = .true.
1589 
1590  if (lboundary1 .and. lboundary2) then
1591 
1592  !*** This is a edge on the boundary of the
1593  !*** active mesh and both of its endpoints
1594  !*** are on the boundary of the containing
1595  !*** cell. Check if the the segment is also
1596  !*** on the boundary
1597 
1598  midlat = (beglat+endlat1)/2.0
1599  if (abs(beglon-endlon1) .ge. pi) then
1600  midlon = (beglon+endlon1)/2.0 - pi
1601  else
1602  midlon = (beglon+endlon1)/2.0
1603  endif
1604 
1605  call ptincell(midlat,midlon, oppcell_add,
1606  & ncorners_opp,
1607  & oppcell_corner_lat, oppcell_corner_lon,
1608  & oppcell_center_lat, oppcell_center_lon,
1609  & opp_grid_num, inpoly, lboundary3,
1610  & bedgeid3)
1611 
1612  if (inpoly .and. lboundary3) then
1613  lcoinc = .true.
1614  intedge = bedgeid3
1615  endif
1616 
1617  else
1618  lcoinc = .false.
1619  endif
1620 
1621  else
1622 
1623  !***
1624  !*** Do an intersection to find out where the
1625  !*** segment exits the cell
1626  !***
1627 
1628  call intersection(cell_add,grid_num,
1629  & beglat, beglon, endlat1, endlon1,
1630  & begseg,
1631  & bedgeid1,
1632  & oppcell_add, ncorners_opp,
1633  & oppcell_corner_lat, oppcell_corner_lon,
1634  & opp_grid_num,
1635  & intrsct_lat, intrsct_lon, intedge,
1636  & sinang2, lcoinc, lthresh)
1637 
1638  if (intedge /= 0) then
1639  intrsct_success = .true.
1640  last_add = oppcell_add ! for next subseg
1641  last_endpt_onedge = .true.
1642  last_endpt_inpoly = .false.
1643  last_lboundary = .true.
1644 
1645  if (.not. lthresh) then
1646  call find_adj_cell(oppcell_add,intedge,
1647  & opp_grid_num,next_add)
1648  if (next_add .ne. 0) then
1649  search = .false.
1650  else
1651  search = .true.
1652  endif
1653  else
1654  search = .true.
1655  endif
1656  endif
1657 
1658  endif
1659 
1660  if (.not. intrsct_success) then
1661 
1662  !*** Offset point and try again
1663 
1664  search = .true.
1665  delta = 2.0*delta
1666  offset = offset + delta
1667  if (offset .gt. vec1_len) then
1668 
1669  ! Punt - exit the intersection loop
1670 
1671  intrsct_lat = endlat1
1672  intrsct_lon = endlon1
1673  last_add = 0
1674  last_lboundary = .false.
1675  exit
1676 
1677  endif
1678  endif
1679 
1680 !NRL if (lcoinc .and. .not. bndedge) then
1681 
1682  if (lcoinc .and. .not. bndedge !NRL
1683  & .and. intedge /= 0) then !NRL
1684 
1685  !***
1686  !*** Segment is coincident with edge of other grid
1687  !*** which means it could belong to one of 2 cells
1688  !*** Choose the cell such that edge that is
1689  !*** coincident with the segment is in the same
1690  !*** dir as the segment
1691 
1692  i = intedge
1693  inext = mod(i,ncorners_opp)+1
1694  vec2_lat = oppcell_corner_lat(inext) -
1695  & oppcell_corner_lat(i)
1696  vec2_lon = oppcell_corner_lon(inext) -
1697  & oppcell_corner_lon(i)
1698 
1699  if (vec2_lon > pi) vec2_lon = vec2_lon - pi2
1700  if (vec2_lon < -pi) vec2_lon = vec2_lon + pi2
1701 
1702  dp = vec1_lat*vec2_lat + vec1_lon*vec2_lon
1703 
1704  if ((.not. lrevers .and. dp .lt. 0) .or.
1705  & (lrevers .and. dp .gt. 0)) then
1706 
1707  !*** Integrals from this segment must be
1708  !*** assigned to the adjacent cell of
1709  !*** opcell_add but only if such an adjacent
1710  !*** cell exists
1711 
1712  call find_adj_cell(oppcell_add, intedge,
1713  & opp_grid_num, adj_add)
1714 
1715  if (adj_add .gt. 0) then
1716  oppcell_add = adj_add
1717 
1718  if (grid_num .eq. 1) then
1719  ncorners_opp = grid2_corners
1720  do i = 1, ncorners_opp
1721  oppcell_corner_lat(i) =
1722  & grid2_corner_lat(i,oppcell_add)
1723  oppcell_corner_lon(i) =
1724  & grid2_corner_lon(i,oppcell_add)
1725  enddo
1726  oppcell_center_lat =
1727  & grid2_center_lat(oppcell_add)
1728  oppcell_center_lon =
1729  & grid2_center_lon(oppcell_add)
1730 
1731  special_cell =
1732  & special_polar_cell2(oppcell_add)
1733  else
1734  ncorners_opp = grid1_corners
1735  do i = 1, ncorners_opp
1736  oppcell_corner_lat(i) =
1737  & grid1_corner_lat(i,oppcell_add)
1738  oppcell_corner_lon(i) =
1739  & grid1_corner_lon(i,oppcell_add)
1740  enddo
1741  oppcell_center_lat =
1742  & grid1_center_lat(oppcell_add)
1743  oppcell_center_lon =
1744  & grid1_center_lon(oppcell_add)
1745 
1746  special_cell =
1747  & special_polar_cell1(oppcell_add)
1748  endif
1749 
1750  if (special_cell) then
1751  call modify_polar_cell(ncorners_opp,
1752  & nalloc_opp, oppcell_corner_lat,
1753  & oppcell_corner_lon)
1754  endif
1755 
1756  endif
1757 
1758  endif
1759 
1760  endif
1761 
1762  else
1763 
1764  !***
1765  !*** Could not locate a viable cell for the segment
1766  !*** start
1767  !***
1768 
1769  if (oppcell_add .eq. 0) then
1770  loutside = .true.
1771 ! lcoinc added by MD
1772  lcoinc = .false.
1773 
1774  !***
1775  !*** Take baby steps to see if any part of the
1776  !*** segment is inside a cell of the other grid
1777  !***
1778 
1779  seg_outside = .false.
1780  delta = vec1_len/100.00
1781  offset = delta
1782  do while (.not. srch_success)
1783 
1784  srchpt_lat = beglat + offset*vec1_lat
1785  srchpt_lon = beglon + offset*vec1_lon
1786 
1787  call locate_point(srchpt_lat, srchpt_lon,
1788  & cell_add, grid_num, opp_grid_num,
1789  & oppcell_add, lboundary1, bedgeid1)
1790 
1791  if (oppcell_add /= 0) then
1792  srch_success = .true.
1793 
1794  !***
1795  !*** Found a point of the segment in the
1796  !*** cell. Do a bisection method to find
1797  !*** the starting point of the segment
1798  !*** in the cell
1799  !***
1800 
1801  call converge_to_bdry(oppcell_add,
1802  & opp_grid_num, ncorners_opp,
1803  & oppcell_corner_lat,
1804  & oppcell_corner_lon,
1805  & oppcell_center_lat,
1806  & oppcell_center_lon,
1807  & srchpt_lat, srchpt_lon,
1808  & beglat, beglon,
1809  & intrsct_lat, intrsct_lon,
1810  & bedgeid1)
1811 
1812  search = .false.
1813  last_endpt_onedge = .true.
1814  next_add = oppcell_add
1815  last_lboundary = .true.
1816 
1817  oppcell_add = 0
1818 
1819  else
1820 
1821  offset = offset + delta
1822 
1823  if (offset .ge. vec1_len) then
1824 ! print *,
1825 ! & 'Segment fully outside grid2'
1826 ! print *, 'Segment of grid1_add',
1827 ! & grid1_add
1828 ! print *, beglat,beglon
1829 ! print *, endlat1,endlon1
1830 
1831  seg_outside = .true.
1832 
1833  intrsct_lat = endlat1
1834  intrsct_lon = endlon1
1835 
1836  search = .true.
1837  last_add = 0
1838  last_lboundary = .false.
1839 
1840  exit ! leave search loop
1841  endif
1842  endif
1843 
1844  enddo
1845 
1846  ! int. loop
1847  if (srch_success .or. seg_outside) exit
1848 
1849  else
1850 
1851  if(is_master)then
1852  print *, 'Unable to move out of last cell'
1853  print *, 'Segment of edge ',corner,
1854  & ' of grid cell ',cell_add
1855  print *, 'Stuck in opposite grid cell ',
1856  & oppcell_add
1857  dist2 =
1858  & (endlat1-begseg(1))*(endlat1-begseg(1)) +
1859  & (endlon1-begseg(2))*(endlon1-begseg(2))
1860  print *, 'Fraction of segment left ',
1861  & vec1_len/sqrt(dist2)
1862  endif
1863  lstuck = .true.
1864 
1865  !***
1866  !*** Punt - just assign the rest of the segment
1867  !*** to the current cell it is stuck in by
1868  !*** tagging the segment endpoint as the
1869  !*** intersection point
1870  !***
1871 
1872  intrsct_lat = endlat1
1873  intrsct_lon = endlon1
1874 
1875  search = .true.
1876  last_add = 0
1877  last_lboundary = .false.
1878 
1879  endif
1880 
1881  exit ! exit the intersection loop
1882 
1883  endif ! if (srch_success) then ... else ....
1884 
1885  end do ! do while (.not. intrsct_success)
1886 
1887  !********************************************************
1888  !*** Compute the line integrals for this subsegment
1889  !********************************************************
1890 
1891  if (oppcell_add /= 0) then
1892  call line_integral(phi_or_theta, weights, num_wts,
1893  & beglon, intrsct_lon, beglat, intrsct_lat,
1894  & cell_center_lat, cell_center_lon,
1895  & oppcell_center_lat, oppcell_center_lon)
1896  else
1897  call line_integral(phi_or_theta, weights, num_wts,
1898  & beglon, intrsct_lon, beglat, intrsct_lat,
1899  & cell_center_lat, cell_center_lon,
1900  & cell_center_lat, cell_center_lon)
1901  endif
1902 
1903  !***
1904  !*** if integrating in reverse order, change
1905  !*** sign of weights
1906  !***
1907 
1908  if (lrevers) then
1909  weights = -weights
1910  endif
1911 
1912  !***
1913  !*** store the appropriate addresses and weights.
1914  !*** also add contributions to cell areas and centroids.
1915  !***
1916 
1917  if (grid_num .eq. 1) then
1918 
1919  if (oppcell_add /= 0) then
1920  if (grid1_mask(cell_add)) then
1921  call store_link_cnsrv(cell_add, oppcell_add,
1922  & weights)
1923 
1924 C$OMP CRITICAL(block1)
1925 !
1926 ! Could have another thread that found an intersection between that
1927 ! cell address and oppcell_add in which case it will try to write
1928 ! into this address - we have to block that until we are finished
1929 !
1930  grid1_frac(cell_add) =
1931  & grid1_frac(cell_add) + weights(1)
1932 
1933  grid2_frac(oppcell_add) =
1934  & grid2_frac(oppcell_add) + weights(num_wts+1)
1935 C$OMP END CRITICAL(block1)
1936  endif
1937 
1938  endif
1939 
1940 C$OMP CRITICAL(block2)
1941  grid1_area(cell_add) = grid1_area(cell_add) +
1942  & weights(1)
1943  grid1_centroid_lat(cell_add) =
1944  & grid1_centroid_lat(cell_add) + weights(2)
1945  grid1_centroid_lon(cell_add) =
1946  & grid1_centroid_lon(cell_add) + weights(3)
1947 C$OMP END CRITICAL(block2)
1948 
1949  else
1950 
1951  !*** swap weights because in store_link_cnsrv
1952  !*** we are always sending in grid1 weights first
1953  !*** and then grid2 weights
1954 
1955  do i = 1, num_wts
1956  rev_weights(num_wts+i) = weights(i)
1957  rev_weights(i) = weights(num_wts+i)
1958  enddo
1959 
1960  if (.not. lcoinc .and. oppcell_add /= 0) then
1961  if (grid1_mask(oppcell_add)) then
1962  call store_link_cnsrv(oppcell_add, cell_add,
1963  & rev_weights)
1964 
1965 C$OMP CRITICAL(block3)
1966 !
1967 ! Could have another thread that found an intersection between that
1968 ! cell address and oppcell_add in which case it will try to write
1969 ! into this address - we have to block that until we are finished
1970 !
1971  grid2_frac(cell_add) =
1972  & grid2_frac(cell_add) + weights(1)
1973 
1974  grid1_frac(oppcell_add) =
1975  & grid1_frac(oppcell_add) + weights(num_wts+1)
1976 C$OMP END CRITICAL(block3)
1977 
1978  endif
1979 
1980  endif
1981 
1982 C$OMP CRITICAL(block4)
1983  grid2_area(cell_add) = grid2_area(cell_add) +
1984  & weights(1)
1985  grid2_centroid_lat(cell_add) =
1986  & grid2_centroid_lat(cell_add) + weights(2)
1987  grid2_centroid_lon(cell_add) =
1988  & grid2_centroid_lon(cell_add) + weights(3)
1989 C$OMP END CRITICAL(block4)
1990  endif
1991 
1992  !***
1993  !*** reset beglat and beglon for next subsegment.
1994  !***
1995 
1996  beglat = intrsct_lat
1997  beglon = intrsct_lon
1998 
1999  !***
2000  !*** How far have we come from the start of the segment
2001  !***
2002 
2003  vec2_lat = intrsct_lat-begseg(1)
2004  vec2_lon = intrsct_lon-begseg(2)
2005  if (vec2_lon > pi) vec2_lon = vec2_lon - pi2
2006  if (vec2_lon < -pi) vec2_lon = vec2_lon + pi2
2007 
2008  partseg_len2 = vec2_lat*vec2_lat + vec2_lon*vec2_lon
2009 
2010  !***
2011  !*** prevent infinite loops if integration gets stuck
2012  !*** near cell or threshold boundary
2013  !***
2014 
2015  num_subseg = num_subseg + 1
2016  if (num_subseg > max_subseg) then
2017  print *,
2018  & 'integration stalled: num_subseg exceeded limit'
2019  print *, 'Cell ',cell_add
2020  print *, 'Edge ',corner
2021  print *, 'Grid ',1
2022  dist2 = (endlat1-begseg(1))*(endlat1-begseg(1)) +
2023  & (endlon1-begseg(2))*(endlon1-begseg(2))
2024  print *, 'Fraction of segment left ',
2025  & vec1_len/sqrt(dist2)
2026 ! exit ! Give up and exit
2027  stop ! Give up and stop
2028  endif
2029 
2030 ! inner "do while"
2031  end do ! do while (beglat /= endlat1 ...
2032 
2033 !NRL We add an exit to outer do similar to exit of inner do:
2034 !NRL This was an apparent bug: exit statement would escape
2035 !NRL inner do but then computation could not get out of
2036 !NRL outer do since beglat, beglon controlling outer do
2037 !NRL never changed b/c it never gets to the part of the
2038 !NRL code that changes beglat, beglon, b/c it keeps
2039 !NRL exiting inner do.
2040 
2041 !NRL This should happen very rarely, so we have a print
2042 !NRL statement to notify user.
2043 
2044  if (l_exit_do)then ! NRL
2045  write(*,*)'partseg_len2,fullseg_len2 = ', ! NRL
2046  & partseg_len2,fullseg_len2 ! NRL
2047  write(*,*)'exiting outer do while loop' ! NRL
2048  exit ! NRL
2049  endif ! NRL
2050 
2051  ns = ns + 1
2052  if ((beglat > 0 .and. beglat < north_thresh) .or.
2053  & (beglat < 0 .and. beglat > south_thresh))
2054  & then
2055  inpolar = .false.
2056  endif
2057 
2058 ! outer "do while"
2059  end do ! do while (beglat /= endlat ....
2060 
2061  call line_integral(phi_or_theta, weights, num_wts,
2062  & begseg(2), endlon, begseg(1), endlat,
2063  & cell_center_lat,
2064  & cell_center_lon,
2065  & cell_center_lat,
2066  & cell_center_lon)
2067 
2068  !***
2069  !*** end of segment
2070  !***
2071 
2072  end do ! do corner=....
2073 
2074  end subroutine cell_integrate
2075 !***********************************************************************
2076 
2077 
2078 !***********************************************************************
2079 
2080  subroutine modify_polar_cell(ncorners, nalloc, cell_corner_lat,
2081  & cell_corner_lon)
2082 
2083  !*** Input variables
2084 
2085  integer (SCRIP_i4), intent(in) ::
2086  & nalloc
2087 
2088  !*** In/Out Variables
2089 
2090  integer (SCRIP_i4), intent(inout) ::
2091  & ncorners
2092  real (SCRIP_r8), dimension(:), intent(inout) ::
2093  & cell_corner_lat(nalloc),
2094  & cell_corner_lon(nalloc)
2095 
2096  !*** Local variables
2097 
2098  integer (SCRIP_i4) ::
2099  & npcorners, ! Number of polar corners
2100  & pcorner, ! Index of the polar corner
2101  ! (if only 1 is found)
2102  & corner, ! Corner iterator variable
2103  & previdx, ! Index of previous corner to polar corner
2104  & nextidx ! Index of next corner to polar corner
2105 
2106  real (SCRIP_r8) ::
2107  & pole_lat, ! Latitude considered to be pole
2108  & prevlon, ! Latitude of previous corner to polar corner
2109  & nextlon ! Latitude of next corner to polar corner
2110 
2111 
2112  !***
2113  !*** Modify special cell with only one corner at the pole. Such
2114  !*** cells can have an artificially extreme distortion of the
2115  !*** edges when mapped to the Lambert projection because of the
2116  !*** span of longitudes on the edges So we will augment such
2117  !*** cells with degenerate edges at the pole so that the edges
2118  !*** coming off the pole will actually have the same longitude
2119  !*** values at both ends
2120  !***
2121  !*** lon_p lon_p+ lon_p lon_p-
2122  !*** pi/2 pi/2 pi/2 pi/2
2123  !*** * *--------*------*
2124  !*** / \ | |
2125  !*** / \ | |
2126  !*** / \ | |
2127  !*** * * * *
2128  !*** lon_p+ lon_p- lon_p+ lon_p-
2129  !*** lat_p+ lat_p- lat_p+ lat_p-
2130  !***
2131 
2132 
2133  !***
2134  !*** MAJOR ASSUMPTION HERE IS THAT CELL_CORNER_LAT AND
2135  !*** CELL_CORNER_LON HAVE ROOM TO GROW
2136  !***
2137  if (ncorners .ge. nalloc) return ! ** * No room to grow
2138 
2139  pcorner = 0
2140  npcorners = 0
2141  do corner = 1, ncorners
2142  if (abs(abs(cell_corner_lat(corner))-pih) .le. 1.0e-05) then
2143  pcorner = corner
2144  pole_lat = cell_corner_lat(corner)
2145  npcorners = npcorners + 1
2146  endif
2147  enddo
2148 
2149 
2150  if (npcorners .ne. 1) return !*** Not the kind of cell we want
2151 
2152  previdx = mod((pcorner-1)-1+ncorners,ncorners) + 1
2153  prevlon = cell_corner_lon(previdx)
2154 
2155  nextidx = mod(pcorner,ncorners) + 1
2156  nextlon = cell_corner_lon(nextidx)
2157 
2158  !*** Move entries from pcorner+1 on back by one
2159 
2160  do corner = ncorners, pcorner+1, -1
2161  cell_corner_lat(corner+1) = cell_corner_lat(corner)
2162  cell_corner_lon(corner+1) = cell_corner_lon(corner)
2163  enddo
2164 
2165  !*** Add a corner after pcorner
2166 
2167  cell_corner_lat(pcorner+1) = pole_lat
2168  cell_corner_lon(pcorner+1) = nextlon
2169 
2170  ncorners = ncorners+1
2171 
2172  !*** Move entries from pcorner on back by one
2173 
2174  do corner = ncorners, pcorner, -1
2175  cell_corner_lat(corner+1) = cell_corner_lat(corner)
2176  cell_corner_lon(corner+1) = cell_corner_lon(corner)
2177  enddo
2178 
2179  !*** Add a corner before pcorner
2180 
2181  cell_corner_lat(pcorner) = pole_lat
2182  cell_corner_lon(pcorner) = prevlon
2183 
2184  ncorners = ncorners+1
2185 
2186  end subroutine modify_polar_cell
2187 
2188 
2189 !***********************************************************************
2190 
2191  subroutine intersection(seg_cell_id, seg_grid_id,
2192  & beglat, beglon, endlat, endlon, begseg, begedge,
2193  & cell_id, ncorners, cell_corner_lat,
2194  & cell_corner_lon, cell_grid_id, intrsct_lat, intrsct_lon,
2195  & intedge, sinang2, lcoinc, lthresh)
2196 
2197 !-----------------------------------------------------------------------
2198 !
2199 ! this routine finds the intersection of a line segment given by
2200 ! beglon, endlon, etc. with a cell from another grid
2201 ! A coincidence flag is returned if the segment is entirely
2202 ! coincident with an edge of the opposite.
2203 !
2204 !-----------------------------------------------------------------------
2205 
2206 !-----------------------------------------------------------------------
2207 !
2208 ! intent(in):
2209 !
2210 !-----------------------------------------------------------------------
2211 
2212  integer (SCRIP_i4), intent(in) ::
2213  & seg_cell_id ! ID of cell that intersecting segment is from
2214 
2215  integer (SCRIP_i4), intent(in) ::
2216  & seg_grid_id ! ID of grid that intersecting segment is from
2217 
2218  real (SCRIP_r8), intent(in) ::
2219  & beglat, beglon,! beginning lat/lon endpoints for segment
2220  & endlat, endlon ! ending lat/lon endpoints for segment
2221 
2222  real (SCRIP_r8), dimension(2), intent(inout) ::
2223  & begseg ! begin lat/lon of full segment
2224 
2225  integer (SCRIP_i4), intent(in) ::
2226  & begedge ! edge that beginning point is on (can be 0)
2227 
2228  integer (SCRIP_i4), intent(in) ::
2229  & cell_id ! cell to intersect with
2230 
2231  integer (SCRIP_i4), intent(in) ::
2232  & ncorners ! number of corners of cell
2233 
2234  real (SCRIP_r8), dimension(ncorners), intent(in) ::
2235  & cell_corner_lat, ! coordinates of cell corners
2236  & cell_corner_lon
2237 
2238  integer (SCRIP_i4), intent(in) ::
2239  & cell_grid_id ! which grid is the cell from?
2240 
2241 
2242 !-----------------------------------------------------------------------
2243 !
2244 ! intent(out):
2245 !
2246 !-----------------------------------------------------------------------
2247 
2248  real (SCRIP_r8), intent(out) ::
2249  & intrsct_lat,
2250  & intrsct_lon ! lat/lon coords of intersection
2251 
2252  real (SCRIP_r8), intent(out) ::
2253  & sinang2 ! square of sine of angle between
2254  ! intersecting lines
2255 
2256  integer (SCRIP_i4), intent(out) ::
2257  & intedge ! edge that is intersected
2258 
2259  logical (SCRIP_logical), intent(out) ::
2260  & lcoinc ! True if segment is coincident with
2261  ! a cell edge
2262 
2263  logical (SCRIP_logical), intent(out) ::
2264  & lthresh ! True if segment crosses threshold
2265 
2266 !-----------------------------------------------------------------------
2267 !
2268 ! local variables
2269 !
2270 !-----------------------------------------------------------------------
2271 
2272  integer (SCRIP_i4) ::
2273  & n, next_n
2274 
2275  logical (SCRIP_logical) ::
2276  & found, first
2277 
2278  real (SCRIP_r8) ::
2279  & lon1, lon2, ! local longitude variables for segment
2280  & lat1, lat2, ! local latitude variables for segment
2281  & grdlon1, grdlon2, ! local longitude variables for grid cell
2282  & grdlat1, grdlat2, ! local latitude variables for grid cell
2283  & vec1_lat, vec1_lon,
2284  & vec2_lat, vec2_lon, !
2285  & vec3_lat, vec3_lon, ! vectors and vector products used
2286  & cross_product, ! during grid search
2287  & dot_product, !
2288  & lensqr1, lensqr2, !
2289  & lensqr3, !
2290  & s1, s2, determ,
2291  & mat1, mat2, ! variables used for linear solve to
2292  & mat3, mat4, ! find intersection
2293  & rhs1, rhs2, !
2294  & denom,
2295  & begsegloc(2), ! local copy of full segment start
2296  & dist2, ! distance from start pt to intersection
2297  ! pt
2298  & maxdist2, ! max dist from start pt to any
2299  ! intersection pt
2300  & max_intrsct_lat, ! latitude of farthest intersection point
2301  & max_intrsct_lon, ! longitude of farthest intersection
2302  ! point
2303  & minlat, maxlat, ! min and max latitudes of segment
2304  & minlon, maxlon, ! min and max longitudes of segment
2305  & tmplat, tmplon
2306 
2307 
2308 !-----------------------------------------------------------------------
2309 !
2310 ! initialize defaults, flags, etc.
2311 !
2312 !-----------------------------------------------------------------------
2313 
2314  lcoinc = .false.
2315  lthresh = .false.
2316  intedge = 0
2317  first = .true.
2318 
2319  lat1 = beglat
2320  lon1 = beglon
2321  lat2 = endlat
2322  lon2 = endlon
2323 
2324  ! No edge is allowed to span more than pi radians
2325  ! Accordingly transform one or the other end point
2326 
2327  if ((lon2-lon1) > pi) then
2328  lon2 = lon2 - pi2
2329  else if ((lon2-lon1) < -pi) then
2330  lon1 = lon1 - pi2
2331  endif
2332  s1 = zero
2333 
2334 !-----------------------------------------------------------------------
2335 !
2336 ! loop over sides of the cell to find intersection with side
2337 ! must check all sides for coincidences or intersections
2338 !
2339 !-----------------------------------------------------------------------
2340 
2341  if (beglat > north_thresh .or. beglat < south_thresh) then
2342 
2343  !*** Special intersection routine for cells near the pole
2344  !*** Intersection is done in a transformed space using
2345  !*** multi-segmented representation of the cell
2346 
2347  call pole_intersection(cell_id,ncorners,
2348  & cell_corner_lat,cell_corner_lon,cell_grid_id,
2349  & beglat, beglon, endlat,
2350  & endlon, begseg, begedge,
2351  & intedge,intrsct_lat,intrsct_lon,
2352  & sinang2,lcoinc,lthresh)
2353 
2354  return
2355 
2356  endif
2357 
2358 
2359  maxdist2 = -9999999.0
2360 
2361  begsegloc(1) = begseg(1)
2362  begsegloc(2) = begseg(2)
2363 
2364  lthresh = .false.
2365  intrsct_loop: do n=1,ncorners
2366  next_n = mod(n,ncorners) + 1
2367 
2368  grdlat1 = cell_corner_lat(n)
2369  grdlon1 = cell_corner_lon(n)
2370  grdlat2 = cell_corner_lat(next_n)
2371  grdlon2 = cell_corner_lon(next_n)
2372 
2373  lensqr2 = (grdlat1-grdlat2)*(grdlat1-grdlat2) +
2374  & (grdlon1-grdlon2)*(grdlon1-grdlon2)
2375 
2376  if (lensqr2 .le. tiny*tiny) cycle ! degenerate edge
2377 
2378  ! No edge can span more than pi radians
2379 
2380  if (grdlon2-grdlon1 > pi) then
2381  grdlon2 = grdlon2 - pi2
2382  else if (grdlon2-grdlon1 < -pi) then
2383  grdlon1 = grdlon1 - pi2
2384  endif
2385 
2386  ! Also the two intersecting segments together
2387  ! cannot span more than 2*pi radians
2388 
2389  minlon = min(lon1,lon2)
2390  maxlon = max(grdlon1,grdlon2)
2391  if (maxlon-minlon > pi2) then
2392  grdlon1 = grdlon1 - pi2
2393  grdlon2 = grdlon2 - pi2
2394  else
2395  minlon = min(grdlon1,grdlon2)
2396  maxlon = max(lon1,lon2)
2397  if (maxlon-minlon > pi2) then
2398  grdlon1 = grdlon1 + pi2
2399  grdlon2 = grdlon2 + pi2
2400  endif
2401  endif
2402 
2403 
2404  !***
2405  !*** set up linear system to solve for intersection
2406  !***
2407 
2408  mat1 = lat2 - lat1
2409  mat2 = grdlat1 - grdlat2
2410  mat3 = lon2 - lon1
2411  mat4 = grdlon1 - grdlon2
2412  rhs1 = grdlat1 - lat1
2413  rhs2 = grdlon1 - lon1
2414 
2415  determ = mat1*mat4 - mat2*mat3
2416 
2417  !***
2418  !*** if the determinant is zero, the segments are either
2419  !*** parallel or coincident. coincidences were detected
2420  !*** above so do nothing.
2421 
2422  if (abs(determ) > tiny*tiny) then
2423 
2424  !*** if the determinant is non-zero, solve for the linear
2425  !*** parameters s for the intersection point on each line
2426  !*** segment.
2427  !*** if 0<s1,s2<1 then the segment intersects with this side.
2428  !*** return the point of intersection (adding a small
2429  !*** number so the intersection is off the grid line).
2430  !***
2431 
2432  s1 = (rhs1*mat4 - mat2*rhs2)/determ
2433  s2 = (mat1*rhs2 - rhs1*mat3)/determ
2434 
2435  if (s2 >= zero .and. s2 <= one .and.
2436  & s1 > zero .and. s1 <= one) then
2437 
2438  !***
2439  !*** recompute intersection based on full segment
2440  !*** so intersections are consistent for both sweeps
2441  !***
2442 
2443  if (lon2-begsegloc(2) > pi) then
2444  lon2 = lon2 - pi2
2445  else if (lon2-begsegloc(2) < -pi) then
2446  begsegloc(2) = begsegloc(2) - pi2
2447  endif
2448 
2449 
2450  ! Also the two intersecting segments together
2451  ! cannot span more than 2*pi radians
2452 
2453  minlon = min(begsegloc(2),lon2)
2454  maxlon = max(grdlon1,grdlon2)
2455  if (maxlon-minlon > pi2) then
2456  grdlon1 = grdlon1 - pi2
2457  grdlon2 = grdlon2 - pi2
2458  else
2459  minlon = min(grdlon1,grdlon2)
2460  maxlon = max(begsegloc(2),lon2)
2461  if (maxlon-minlon > pi2) then
2462  grdlon1 = grdlon1 + pi2
2463  grdlon2 = grdlon2 + pi2
2464  endif
2465  endif
2466 
2467 
2468  mat1 = lat2 - begsegloc(1)
2469  mat3 = lon2 - begsegloc(2)
2470  rhs1 = grdlat1 - begsegloc(1)
2471  rhs2 = grdlon1 - begsegloc(2)
2472 
2473  determ = mat1*mat4 - mat2*mat3
2474 
2475  !***
2476  !*** sometimes due to roundoff, the previous
2477  !*** determinant is non-zero, but the lines
2478  !*** are actually coincident. if this is the
2479  !*** case, skip the rest.
2480  !***
2481 
2482  if (determ /= zero) then
2483  s1 = (rhs1*mat4 - mat2*rhs2)/determ
2484  s2 = (mat1*rhs2 - rhs1*mat3)/determ
2485 
2486  intrsct_lat = begsegloc(1) + mat1*s1
2487  intrsct_lon = begsegloc(2) + mat3*s1
2488 
2489  if (intrsct_lon < 0.0) then
2490  intrsct_lon = intrsct_lon + pi2
2491  else if (intrsct_lon > pi2) then
2492  intrsct_lon = intrsct_lon - pi2
2493  endif
2494 
2495  !***
2496  !*** Make sure the intersection point is not within
2497  !*** tolerance of the starting point
2498  !***
2499 
2500  if (first) then
2501  max_intrsct_lat = intrsct_lat
2502  max_intrsct_lon = intrsct_lon
2503 
2504  vec1_lat = intrsct_lat-beglat
2505  vec1_lon = intrsct_lon-beglon
2506  if (vec1_lon > pi) then
2507  vec1_lon = vec1_lon - pi2
2508  else if (vec1_lon < -pi) then
2509  vec1_lon = vec1_lon + pi2
2510  endif
2511 
2512  maxdist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
2513  dist2 = maxdist2
2514 
2515  denom = (mat1*mat1+mat2*mat2)*(mat3*mat3+mat4*mat4)
2516  sinang2 = determ*determ/denom
2517  intedge = n
2518  first = .false.
2519  else
2520  vec1_lat = intrsct_lat-beglat
2521  vec1_lon = intrsct_lon-beglon
2522  if (vec1_lon > pi) then
2523  vec1_lon = vec1_lon - pi2
2524  else if (vec1_lon < -pi) then
2525  vec1_lon = vec1_lon + pi2
2526  endif
2527 
2528  dist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
2529 
2530  if (dist2 > maxdist2) then
2531  if (begedge .eq. 0 .or. begedge .ne. n) then
2532  max_intrsct_lat = intrsct_lat
2533  max_intrsct_lon = intrsct_lon
2534  maxdist2 = dist2
2535 
2536  denom =
2537  & (mat1*mat1+mat2*mat2)*(mat3*mat3+mat4*mat4)
2538  sinang2 = determ*determ/denom
2539  intedge = n
2540  endif
2541  endif
2542  endif
2543 
2544  else
2545  print *, 'DEBUG: zero determ'
2546  stop
2547  endif
2548 
2549  endif
2550 
2551  else
2552 
2553  !***
2554  !*** Coincident lines or parallel lines
2555  !***
2556 
2557  cross_product = mat2*rhs2 - mat4*rhs1
2558 
2559  !***
2560  !*** If area of triangle formed by endlat,endlon and
2561  !*** the gridline is negligible then the lines are coincident
2562  !***
2563 
2564 
2565  if (abs(cross_product) < tiny) then
2566 
2567  dot_product = mat1*(-mat2) + mat3*(-mat4)
2568 
2569  lensqr1 = mat1*mat1 + mat3*mat3 ! length sqrd of input
2570  ! segment
2571 
2572  if (dot_product < zero) then
2573 
2574  !***
2575  !*** Segments oriented in the same direction
2576  !***
2577 
2578 
2579  tmplat = grdlat2
2580  tmplon = grdlon2
2581  grdlat2 = grdlat1
2582  grdlon2 = grdlon1
2583  grdlat1 = tmplat
2584  grdlon1 = tmplon
2585 
2586  endif
2587 
2588 
2589  vec2_lat = grdlat1 - lat1
2590  vec2_lon = grdlon1 - lon1
2591  if (vec2_lon > pi) vec2_lon = vec2_lon - pi2
2592  if (vec2_lon < -pi) vec2_lon = vec2_lon + pi2
2593 
2594  lensqr2 = vec2_lat*vec2_lat + vec2_lon*vec2_lon
2595 
2596  if (vec2_lat*mat1 + vec2_lon*mat3 < 0) then
2597  lensqr2 = -lensqr2
2598  endif
2599 
2600  vec3_lat = grdlat2 - lat1
2601  vec3_lon = grdlon2 - lon1
2602  if (vec3_lon > pi) vec3_lon = vec3_lon - pi2
2603  if (vec3_lon < -pi) vec3_lon = vec3_lon + pi2
2604 
2605  lensqr3 = (vec3_lat*vec3_lat+vec3_lon*vec3_lon)
2606 
2607  if (vec3_lat*mat1 + vec3_lon*mat3 < 0) then
2608  lensqr3 = -lensqr3
2609  endif
2610 
2611  found = .false.
2612 
2613  if (lensqr2 > 0) then
2614  if (lensqr2 <= lensqr1) then
2615  intrsct_lat = grdlat1
2616  intrsct_lon = grdlon1
2617  found = .true.
2618  endif
2619  else
2620  if (lensqr3 > 0) then
2621  if (lensqr3 > lensqr1) then
2622  intrsct_lat = lat2
2623  intrsct_lon = lon2
2624  found = .true.
2625  else
2626  intrsct_lat = grdlat2
2627  intrsct_lon = grdlon2
2628  found = .true.
2629  endif
2630  endif
2631  endif
2632 
2633  if (found) then
2634 
2635  dist2 = (intrsct_lat-beglat)*(intrsct_lat-beglat)+
2636  & (intrsct_lon-beglon)*(intrsct_lon-beglon)
2637 
2638  !*** Coincidence intersection always wins
2639 
2640  max_intrsct_lat = intrsct_lat
2641  max_intrsct_lon = intrsct_lon
2642  maxdist2 = dist2
2643  sinang2 = 0
2644  intedge = n
2645  lcoinc = .true.
2646 
2647  exit intrsct_loop
2648  endif
2649 
2650  endif
2651 
2652  endif
2653 
2654  !*** restore lon1 and lon2 in case it got modified
2655 
2656  lon1 = beglon
2657  lon2 = endlon
2658  begsegloc(2) = begseg(2)
2659  if ((lon2-lon1) > pi) then
2660  lon2 = lon2 - pi2
2661  else if ((lon2-lon1) < -pi) then
2662  lon1 = lon1 - pi2
2663  endif
2664 
2665  end do intrsct_loop
2666 
2667  if (intedge .eq. 0) then
2668  return
2669  else
2670  if (maxdist2 < 1e6*tiny*tiny) then
2671  intedge = 0
2672  return
2673  else
2674  intrsct_lat = max_intrsct_lat
2675  intrsct_lon = max_intrsct_lon
2676  endif
2677  endif
2678 
2679 !-----------------------------------------------------------------------
2680 !
2681 ! if the segment crosses a pole threshold, reset the intersection
2682 ! to be the threshold latitude. only check if this was not a
2683 ! threshold segment since sometimes coordinate transform can end
2684 ! up on other side of threshold again.
2685 !
2686 !-----------------------------------------------------------------------
2687 
2688  if (lthresh) then
2689  if (intrsct_lat < north_thresh .or. intrsct_lat > south_thresh)
2690  & lthresh = .false.
2691  else if (lat1 > zero .and. intrsct_lat > north_thresh) then
2692 ! intrsct_lat = north_thresh + tiny
2693  intrsct_lat = north_thresh
2694  mat1 = lat2 - begsegloc(1)
2695  mat3 = lon2 - begsegloc(2)
2696  s1 = (intrsct_lat - begsegloc(1))/mat1
2697  intrsct_lon = begsegloc(2) + s1*mat3
2698  lthresh = .true.
2699  else if (lat1 < zero .and. intrsct_lat < south_thresh) then
2700 ! intrsct_lat = south_thresh - tiny
2701  intrsct_lat = south_thresh
2702  mat1 = lat2 - begsegloc(1)
2703  mat3 = lon2 - begsegloc(2)
2704  s1 = (intrsct_lat - begsegloc(1))/mat1
2705  intrsct_lon = begsegloc(2) + s1*mat3
2706  lthresh = .true.
2707  endif
2708 
2709  if (intrsct_lon < 0.0) then
2710  intrsct_lon = intrsct_lon + pi2
2711  else if (intrsct_lon > pi2) then
2712  intrsct_lon = intrsct_lon - pi2
2713  endif
2714 
2715 
2716 !-----------------------------------------------------------------------
2717 
2718  end subroutine intersection
2719 
2720 !***********************************************************************
2721 
2722 
2723  subroutine pole_intersection(location,ncorners,
2724  & cell_corners_lat,cell_corners_lon,cell_grid_id,
2725  & beglat, beglon, endlat, endlon, begseg, begedge,
2726  & intedge, intrsct_lat, intrsct_lon,
2727  & sinang2, lcoinc, lthresh)
2728 
2729 !-----------------------------------------------------------------------
2730 !
2731 ! Special intersection routine for line segment in cell close to
2732 ! poles
2733 ! A coordinate transformation (using a Lambert azimuthal
2734 ! equivalent projection) is performed to perform the intersection
2735 ! Also, since a straight line in lat-lon space is a curve in this
2736 ! transformed space, we represent each edge of the cell as having
2737 ! 'npseg' segments whose endpoints are mapped using the Lambert
2738 ! projection
2739 !
2740 !-----------------------------------------------------------------------
2741 
2742 !-----------------------------------------------------------------------
2743 !
2744 ! intent(in):
2745 !
2746 !-----------------------------------------------------------------------
2747 
2748  integer (SCRIP_i4), intent(in) ::
2749  & location ! cell to intersect segment with
2750 
2751  integer (SCRIP_i4), intent(in) ::
2752  & ncorners ! Number of cell corners
2753 
2754  real (SCRIP_r8), dimension(ncorners), intent(in) ::
2755  & cell_corners_lat, ! Cell corner coordinates
2756  & cell_corners_lon
2757 
2758  integer (SCRIP_i4), intent(in) ::
2759  & cell_grid_id ! which grid is the cell from?
2760 
2761  real (SCRIP_r8), intent(in) ::
2762  & beglat, beglon, ! beginning lat/lon coords for segment
2763  & endlat, endlon ! ending lat/lon coords for segment
2764 
2765  real (SCRIP_r8), dimension(2), intent(inout) ::
2766  & begseg ! begin lat/lon of full segment
2767 
2768  integer (SCRIP_i4), intent(in) ::
2769  & begedge ! edge on which segment start is on
2770  ! (can be 0)
2771 
2772 !-----------------------------------------------------------------------
2773 !
2774 ! intent(out):
2775 !
2776 !-----------------------------------------------------------------------
2777 
2778  integer (SCRIP_i4), intent(out) ::
2779  & intedge ! Edge that segment intersects
2780 
2781  real (SCRIP_r8), intent(out) ::
2782  & intrsct_lat, ! lat/lon coords of intersection
2783  & intrsct_lon
2784 
2785  real (SCRIP_r8), intent(out) ::
2786  & sinang2 ! square of sine of angle between
2787  ! intersecting line segments
2788 
2789  logical (SCRIP_logical), intent(out) ::
2790  & lcoinc ! True if segment is coincident with
2791  ! a cell edge
2792 
2793  logical (SCRIP_logical), intent(inout) ::
2794  & lthresh ! True if segment crosses threshold
2795 
2796 
2797 !-----------------------------------------------------------------------
2798 !
2799 ! local variables
2800 !
2801 !-----------------------------------------------------------------------
2802 
2803  integer (SCRIP_i4) ::
2804  & n, n1, next_n, prev_n,
2805  & it, i, j,
2806  & ncorners2,
2807  & intedge1
2808 
2809  logical (SCRIP_logical) ::
2810  & first,
2811  & found
2812 
2813  real (SCRIP_r8) ::
2814  & pi4, rns, ! north/south conversion
2815  & x1, x2, ! local x variables for segment
2816  & y1, y2, ! local y variables for segment
2817  & grdx1, grdx2, ! local x variables for grid cell
2818  & grdy1, grdy2, ! local y variables for grid cell
2819  & grdlat1, grdlat2, ! latitude vars for grid cell
2820  & grdlon1, grdlon2, ! longitude vars for grid cell
2821  & vec1_y, vec1_x, !
2822  & vec2_y, vec2_x, ! vectors and cross products used
2823  & vec3_y, vec3_x, !
2824  & vec1_lat, vec1_lon, !
2825  & vec2_lat, vec2_lon, !
2826  & vec3_lon, !
2827  & cross_product, !
2828  & dot_product, !
2829  & s1, s2, determ, ! variables used for linear solve to
2830  & mat1, mat2, !
2831  & mat3, mat4, ! find intersection
2832  & rhs1, rhs2, !
2833  & denom, !
2834  & intrsct_x, intrsct_y, ! intersection coordinates in
2835  ! transformed space
2836  & max_intrsct_lat, ! intersection point at max distance
2837  & max_intrsct_lon, ! from the start point
2838  & dist2, ! dist of intersection point from start
2839  ! point
2840  & maxdist2, ! max dist of intersection point from
2841  ! start pnt
2842  & lensqr1, lensqr2, ! various segment lengths
2843  & lensqr3,
2844  & tmpx, tmpy,
2845  & tmplat, tmplon,
2846  & ldummy
2847 
2848  !***
2849  !*** variables necessary if segment manages to hit pole
2850  !***
2851 
2852  real (SCRIP_r8), dimension(npseg*ncorners) ::
2853  & cell_corners_lat_loc,! Lat/Lon coordinates of multi-segmented
2854  & cell_corners_lon_loc ! version of cell
2855 
2856 
2857 
2858 !-----------------------------------------------------------------------
2859 !
2860 ! initialize defaults, flags, etc.
2861 !
2862 !-----------------------------------------------------------------------
2863 
2864  max_intrsct_lat = pi ! intersection point at max distance
2865  max_intrsct_lon = 4*pi ! from the start point
2866 
2867  intedge = 0
2868  first = .true.
2869  maxdist2 = -999999.00
2870 
2871  s1 = zero
2872 
2873 !-----------------------------------------------------------------------
2874 !
2875 ! convert coordinates
2876 !
2877 !-----------------------------------------------------------------------
2878 
2879  if (beglat > zero) then
2880  pi4 = quart*pi
2881  rns = one
2882  else
2883  pi4 = -quart*pi
2884  rns = -one
2885  endif
2886 
2887  x1 = rns*two*sin(pi4 - half*beglat)*cos(beglon)
2888  y1 = two*sin(pi4 - half*beglat)*sin(beglon)
2889  x2 = rns*two*sin(pi4 - half*endlat)*cos(endlon)
2890  y2 = two*sin(pi4 - half*endlat)*sin(endlon)
2891 
2892  intrsct_x = x2
2893  intrsct_y = y2
2894 
2895 
2896 !-----------------------------------------------------------------------
2897 !
2898 ! now that a cell is found, search for the next intersection.
2899 ! loop over sides of the cell to find intersection with side
2900 ! must check all sides for coincidences or intersections
2901 !
2902 !-----------------------------------------------------------------------
2903 
2904 
2905  if (abs(x1) .le. tiny .and. abs(y1) .le. tiny .and.
2906  & abs(x2) .le. tiny .and. abs(y2) .le. tiny) then
2907 
2908  !***
2909  !*** The segment is a polar segment which is degenerate
2910  !*** in the transformed Lambert space. Find out which
2911  !*** cell edge it is coincident with and find the
2912  !*** point where the segment exits this cell (if at all)
2913  !*** NOTE 1: THIS MUST BE DONE IN LAT-LON SPACE
2914  !*** NOTE 2: CODE RELEVANT ONLY FOR INTEGRATION W.R.T. phi
2915  !***
2916 
2917  intrsct_loop1: do n = 1, ncorners
2918  next_n = mod(n,ncorners) + 1
2919 
2920  grdlat1 = cell_corners_lat(n)
2921  grdlon1 = cell_corners_lon(n)
2922  grdlat2 = cell_corners_lat(next_n)
2923  grdlon2 = cell_corners_lon(next_n)
2924  grdx1 = rns*two*sin(pi4 - half*grdlat1)*cos(grdlon1)
2925  grdy1 = two*sin(pi4 - half*grdlat1)*sin(grdlon1)
2926  grdx2 = rns*two*sin(pi4 - half*grdlat2)*cos(grdlon2)
2927  grdy2 = two*sin(pi4 - half*grdlat2)*sin(grdlon2)
2928 
2929  if (abs(grdx1) .le. tiny .and. abs(grdy1) .le. tiny .and.
2930  & abs(grdx2) .le. tiny .and. abs(grdy2) .le. tiny) then
2931 
2932  !***
2933  !*** Found polar segment in cell
2934  !***
2935 
2936  vec1_lon = endlon-beglon
2937  if (vec1_lon .gt. pi) vec1_lon = vec1_lon - pi2
2938  if (vec1_lon .lt. -pi) vec1_lon = vec1_lon + pi2
2939 
2940  vec2_lon = grdlon2-grdlon1
2941  if (vec2_lon .gt. pi) vec2_lon = vec2_lon - pi2
2942  if (vec2_lon .lt. -pi) vec2_lon = vec2_lon + pi2
2943 
2944  if (vec1_lon*vec2_lon .lt. 0) then
2945 
2946  !*** switch coordinates to simplify logic below
2947 
2948  tmplat = grdlat2
2949  tmplon = grdlon2
2950  grdlat2 = grdlat1
2951  grdlon2 = grdlon1
2952  grdlat1 = tmplat
2953  grdlon1 = tmplon
2954  endif
2955 
2956  vec2_lon = grdlon1 - beglon
2957  if (vec2_lon .gt. pi) vec2_lon = vec2_lon - pi2
2958  if (vec2_lon .lt. -pi) vec2_lon = vec2_lon + pi2
2959 
2960  vec3_lon = grdlon2 - beglon
2961  if (vec3_lon .gt. pi) vec3_lon = vec3_lon - pi2
2962  if (vec3_lon .lt. -pi) vec3_lon = vec3_lon + pi2
2963 
2964  found = .false.
2965 
2966  if (vec2_lon*vec1_lon > 0) then
2967  if (abs(vec3_lon) < abs(vec1_lon)) then
2968  intrsct_lon = grdlon2
2969  found = .true.
2970  else if (abs(vec2_lon) < abs(vec1_lon)) then
2971  intrsct_lon = grdlon1 ! Shouldn't be here
2972  found = .true.
2973  endif
2974  else
2975  if (vec3_lon*vec1_lon > 0) then
2976  if (abs(vec3_lon) < abs(vec1_lon)) then
2977  intrsct_lon = grdlon2
2978  found = .true.
2979  endif
2980  endif
2981 
2982  endif
2983 
2984  if (found) then
2985  intrsct_lat = endlat
2986  lcoinc = .true.
2987  sinang2 = 0.0
2988  intedge = n
2989  exit intrsct_loop1
2990  endif
2991 
2992  endif
2993 
2994  end do intrsct_loop1
2995 
2996  return
2997  endif
2998 
2999 
3000 
3001 
3002  !****
3003  !**** General intersection
3004  !****
3005 
3006 
3007 
3008  !***
3009  !*** Construct multi-segmented version of the cell
3010  !***
3011 
3012  i = 0
3013  do n = ncorners, 1, -1
3014  i = i+1
3015  n1 = mod(n,ncorners)+1
3016  cell_corners_lat_loc(i) = cell_corners_lat(n1)
3017  cell_corners_lon_loc(i) = cell_corners_lon(n1)
3018 
3019  prev_n = n1-1
3020  if (prev_n .eq. 0) prev_n = ncorners ! how do we do (j-1+n)%n
3021  ! in F90 ?
3022 
3023  vec1_lat = cell_corners_lat(prev_n)-cell_corners_lat(n1)
3024  vec1_lon = cell_corners_lon(prev_n)-cell_corners_lon(n1)
3025  if (vec1_lon > pi) then
3026  vec1_lon = vec1_lon - pi2
3027  else if (vec1_lon < -pi) then
3028  vec1_lon = vec1_lon + pi2
3029  endif
3030 
3031  do j = 1, npseg-1
3032  i = i+1
3033  cell_corners_lat_loc(i) =
3034  & cell_corners_lat(n1) + j*vec1_lat/npseg
3035  cell_corners_lon_loc(i) =
3036  & cell_corners_lon(n1) + j*vec1_lon/npseg
3037  enddo
3038  enddo
3039 
3040  ncorners2 = npseg*ncorners
3041 
3042 
3043 
3044  !***
3045  !*** Now intersect segment with multi-segmented version of cell
3046  !***
3047 
3048 
3049  intrsct_loop2: do n= 1, ncorners2
3050 
3051  next_n = mod(n,ncorners2) + 1
3052 
3053  grdlat1 = cell_corners_lat_loc(n)
3054  grdlon1 = cell_corners_lon_loc(n)
3055  grdlat2 = cell_corners_lat_loc(next_n)
3056  grdlon2 = cell_corners_lon_loc(next_n)
3057  grdx1 = rns*two*sin(pi4 - half*grdlat1)*cos(grdlon1)
3058  grdy1 = two*sin(pi4 - half*grdlat1)*sin(grdlon1)
3059  grdx2 = rns*two*sin(pi4 - half*grdlat2)*cos(grdlon2)
3060  grdy2 = two*sin(pi4 - half*grdlat2)*sin(grdlon2)
3061 
3062  if ((grdx1-grdx2)*(grdx1-grdx2)+(grdy1-grdy2)*(grdy1-grdy2) .le.
3063  & tiny*tiny) cycle
3064 
3065 
3066  !***
3067  !*** set up linear system to solve for intersection
3068  !***
3069 
3070  mat1 = x2 - x1
3071  mat2 = grdx1 - grdx2
3072  mat3 = y2 - y1
3073  mat4 = grdy1 - grdy2
3074  rhs1 = grdx1 - x1
3075  rhs2 = grdy1 - y1
3076 
3077  determ = mat1*mat4 - mat2*mat3
3078 
3079  !***
3080  !*** if the determinant is zero, the segments are either
3081  !*** parallel or coincident or one segment has zero length.
3082 
3083  !*** if the determinant is non-zero, solve for the linear
3084  !*** parameters s for the intersection point on each line
3085  !*** segment.
3086  !*** if 0<s1,s2<1 then the segment intersects with this side.
3087  !*** return the point of intersection (adding a small
3088  !*** number so the intersection is off the grid line).
3089  !***
3090 
3091  if (abs(determ) > 1.e-30) then
3092 
3093  s1 = (rhs1*mat4 - mat2*rhs2)/determ
3094  s2 = (mat1*rhs2 - rhs1*mat3)/determ
3095 
3096  if (s2 >= zero .and. s2 <= one .and.
3097  & s1 > tiny .and. s1 <= one) then
3098 
3099  intrsct_x = x1 + s1*mat1
3100  intrsct_y = y1 + s1*mat3
3101 
3102  !***
3103  !*** convert back to lat/lon coordinates
3104  !***
3105 
3106  if (abs(intrsct_x) .gt. tiny .or.
3107  & abs(intrsct_y) .gt. tiny) then
3108 
3109  intrsct_lon = rns*atan2(intrsct_y,intrsct_x)
3110 
3111  else
3112 
3113  !*** Degenerate case - we don't have a good way of
3114  !*** finding out what the longitude corresponding
3115  !*** to a (0,0) intersection is. So we take the
3116  !*** the intersection as one of the two endpoints of
3117  !*** the grid segment
3118 
3119  if (abs(abs(grdlat1)-pih) .lt. 1e-5 .and.
3120  & abs(abs(grdlat2)-pih) .lt. 1e-5) then
3121 
3122  !*** Both endpoints of the grid segment are at the pole
3123  !*** but at different longitudes
3124 
3125  vec1_lat = grdlat1-beglat
3126  vec1_lon = grdlon1-beglon
3127  if (vec1_lon > pi) then
3128  vec1_lon = vec1_lon - pi2
3129  else if (vec1_lon < -pi) then
3130  vec1_lon = vec1_lon + pi2
3131  endif
3132  dist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
3133 
3134  vec2_lat = grdlat2-beglat
3135  vec2_lon = grdlon2-beglon
3136  if (vec2_lon > pi) then
3137  vec2_lon = vec2_lon - pi2
3138  else if (vec2_lon < -pi) then
3139  vec2_lon = vec2_lon + pi2
3140  endif
3141 
3142  !*** pick the endpoint of the grid segment that is
3143  !*** farthest from the beg point of the segment
3144 
3145  if ((vec1_lat*vec1_lat + vec1_lon*vec1_lon) .ge.
3146  & (vec2_lat*vec2_lat + vec2_lon*vec2_lon)) then
3147  intrsct_lon = grdlon1
3148  else
3149  intrsct_lon = grdlon2
3150  endif
3151 
3152  else if (abs(abs(grdlat1)-pih) .lt. 1e-5) then
3153  intrsct_lon = grdlon1
3154  else if (abs(abs(grdlat2)-pih) .lt. 1e-5) then
3155  intrsct_lon = grdlon2
3156  endif
3157 
3158  !*** Make sure this longitude is not outside the
3159  !*** beglon,endlon range
3160 
3161  vec1_lon = endlon-intrsct_lon
3162  if (vec1_lon > pi) then
3163  vec1_lon = vec1_lon - pi2
3164  else if (vec1_lon < -pi) then
3165  vec1_lon = vec1_lon + pi2
3166  endif
3167 
3168  vec2_lon = beglon-intrsct_lon
3169  if (vec2_lon > pi) then
3170  vec2_lon = vec2_lon - pi2
3171  else if (vec2_lon < -pi) then
3172  vec2_lon = vec2_lon + pi2
3173  endif
3174 
3175  !*** if vec1_lon and vec2_lon are of the same sign
3176  !*** then intrsct_lon is outside the beglon,endlon
3177  !*** range
3178 
3179  if (vec1_lon*vec2_lon > 0) cycle
3180 
3181  endif
3182 
3183  if (intrsct_lon < zero)
3184  & intrsct_lon = intrsct_lon + pi2
3185 
3186  if (abs(intrsct_x) > 1.d-10) then
3187  intrsct_lat = (pi4 -
3188  & asin(rns*half*intrsct_x/cos(intrsct_lon)))*two
3189  ldummy = two*(pi4 -
3190  & asin(sqrt(intrsct_x*intrsct_x+intrsct_y*intrsct_y)/2.))
3191  else if (abs(intrsct_y) > 1.d-10) then
3192  intrsct_lat = (pi4 -
3193  & asin(half*intrsct_y/sin(intrsct_lon)))*two
3194  ldummy = two*(pi4 -
3195  & asin(sqrt(intrsct_x*intrsct_x+intrsct_y*intrsct_y)/2.))
3196  else
3197  intrsct_lat = two*pi4
3198  endif
3199 
3200 
3201  !***
3202  !*** If there are multiple intersection points, accept the
3203  !*** one that is not on the edge we started from but is
3204  !*** closest to the start point - need this for
3205  !*** intersection to work for non-convex edges
3206  !***
3207 
3208  if (first) then
3209 
3210  intedge1 = (n-1)/npseg + 1
3211  intedge1 = ncorners - intedge1 + 1 ! dir of edges was
3212  ! reversed
3213  if (intedge1 .ne. begedge) then
3214 
3215  max_intrsct_lat = intrsct_lat
3216  max_intrsct_lon = intrsct_lon
3217 
3218  vec1_lat = intrsct_lat-beglat
3219  vec1_lon = intrsct_lon-beglon
3220  if (vec1_lon > pi) then
3221  vec1_lon = vec1_lon - pi2
3222  else if (vec1_lon < -pi) then
3223  vec1_lon = vec1_lon + pi2
3224  endif
3225  maxdist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
3226  dist2 = maxdist2
3227 
3228  denom = (mat1*mat1+mat2*mat2)*(mat3*mat3+mat4*mat4)
3229  sinang2 = determ*determ/denom
3230  intedge = intedge1
3231 
3232  first = .false.
3233  endif
3234 
3235  else
3236  vec1_lat = intrsct_lat-beglat
3237  vec1_lon = intrsct_lon-beglon
3238  if (vec1_lon > pi) then
3239  vec1_lon = vec1_lon - pi2
3240  else if (vec1_lon < -pi) then
3241  vec1_lon = vec1_lon + pi2
3242  endif
3243  dist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
3244 
3245  !*** if the first intersection was on the same edge
3246  !*** as the starting edge or
3247  !*** the current intersection point is not on the
3248  !*** starting edge and the distance to the beginning
3249  !*** point is less than that of the previous
3250  !*** intersection accept this intersection
3251 
3252  intedge1 = (n-1)/npseg + 1
3253  intedge1 = ncorners - intedge1 + 1 ! dir of edges was
3254  ! reversed
3255  if (dist2 > maxdist2) then
3256  if (begedge == 0 .or. intedge1 .ne. begedge) then
3257  max_intrsct_lat = intrsct_lat
3258  max_intrsct_lon = intrsct_lon
3259  maxdist2 = dist2
3260 
3261  denom =
3262  & (mat1*mat1+mat2*mat2)*(mat3*mat3+mat4*mat4)
3263  sinang2 = determ*determ/denom
3264  intedge = intedge1
3265  endif
3266  endif
3267  endif
3268  endif
3269 
3270  else
3271 
3272  !***
3273  !*** Coincident lines or parallel lines
3274  !***
3275 
3276  cross_product = mat2*rhs2 - mat4*rhs1
3277 
3278  if (abs(cross_product) < tiny) then
3279 
3280  dot_product = mat1*(-mat2) + mat3*(-mat4)
3281 
3282  !***
3283  !*** If area of triangle formed by x2,y2 and the gridline
3284  !*** is negligible then the lines are coincident
3285  !***
3286 
3287  lensqr1 = mat1*mat1 + mat3*mat3 ! length sqrd of input
3288  ! segment
3289 
3290  if (dot_product < zero) then
3291  tmpx = grdx2
3292  tmpy = grdy2
3293  tmplat = grdlat2
3294  tmplon = grdlon2
3295  grdx2 = grdx1
3296  grdy2 = grdy1
3297  grdlat2 = grdlat1
3298  grdlon2 = grdlon1
3299  grdx1 = tmpx
3300  grdy1 = tmpy
3301  grdlat1 = tmplat
3302  grdlon1 = tmplon
3303  endif
3304 
3305 
3306  vec2_x = grdx1 - x1
3307  vec2_y = grdy1 - y1
3308  lensqr2 = vec2_x*vec2_x + vec2_y*vec2_y
3309  if (vec2_x*mat1+vec2_y*mat3 < 0) then
3310  lensqr2 = -lensqr2
3311  endif
3312 
3313 
3314  vec3_x = grdx2 - x1
3315  vec3_y = grdy2 - y1
3316  lensqr3 = (vec3_x*vec3_x+vec3_y*vec3_y)
3317  if (vec3_x*mat1+vec3_y*mat3 < 0) then
3318  lensqr3 = -lensqr3
3319  endif
3320 
3321  found = .false.
3322 
3323  if (lensqr2 > 0) then
3324  if (lensqr2 <= lensqr1) then
3325  intrsct_x = grdx1
3326  intrsct_y = grdy1
3327  intrsct_lat = grdlat1
3328  intrsct_lon = grdlon1
3329  found = .true.
3330  endif
3331  else
3332  if (lensqr3 > 0) then
3333  if (lensqr3 > lensqr1) then
3334  intrsct_x = x2
3335  intrsct_y = y2
3336  intrsct_lat = endlat
3337  intrsct_lon = endlon
3338  found = .true.
3339  else
3340  intrsct_x = grdx2
3341  intrsct_y = grdy2
3342  intrsct_lat = grdlat2
3343  intrsct_lon = grdlon2
3344  found = .true.
3345  endif
3346  endif
3347  endif
3348 
3349  if (found) then
3350  dist2 = (intrsct_lat-beglat)*(intrsct_lat-beglat)+
3351  & (intrsct_lon-beglon)*(intrsct_lon-beglon)
3352 
3353  if (dist2 > tiny*tiny) then
3354 
3355  !*** Coincidence intersection always wins
3356 
3357  max_intrsct_lat = intrsct_lat
3358  max_intrsct_lon = intrsct_lon
3359  maxdist2 = dist2
3360  sinang2 = 0
3361  intedge = (n-1)/npseg + 1
3362  intedge = ncorners - intedge + 1
3363  lcoinc = .true.
3364 
3365  exit intrsct_loop2
3366  endif
3367  endif
3368 
3369 
3370  endif ! if (abs(cross_product) < tiny)
3371 
3372  endif ! if (abs(determ) > 1.e-30) .. else .. endif
3373 
3374  end do intrsct_loop2
3375 
3376  if (maxdist2 < 1e6*tiny*tiny) then
3377  intedge = 0
3378  return
3379  else
3380  intrsct_lat = max_intrsct_lat
3381  intrsct_lon = max_intrsct_lon
3382  endif
3383 
3384 !-----------------------------------------------------------------------
3385 !
3386 ! if segment manages to cross over pole, shift the beginning
3387 ! endpoint in order to avoid hitting pole directly
3388 ! (it is ok for endpoint to be pole point)
3389 !
3390 !-----------------------------------------------------------------------
3391 
3392  if (abs(intrsct_x) < 1.e-10 .and. abs(intrsct_y) < 1.e-10 .and.
3393  & (x2 /= zero .and. y2 /=0)) then
3394  if (avoid_pole_count > 2) then
3395  avoid_pole_count = 0
3396  avoid_pole_offset = 10.*avoid_pole_offset
3397  endif
3398 
3399  cross_product = x1*(y2-y1) - y1*(x2-x1)
3400  intrsct_lat = beglat
3401  if (cross_product*intrsct_lat > zero) then
3402  intrsct_lon = beglon + avoid_pole_offset
3403  else
3404  intrsct_lon = beglon - avoid_pole_offset
3405  endif
3406 
3407  avoid_pole_count = avoid_pole_count + 1
3408  else
3409  avoid_pole_count = 0
3410  avoid_pole_offset = tiny
3411  endif
3412 
3413 !-----------------------------------------------------------------------
3414 !
3415 ! if the segment crosses a pole threshold, reset the intersection
3416 ! to be the threshold latitude and do not reuse x,y intersect
3417 ! on next entry. only check if did not cross threshold last
3418 ! time - sometimes the coordinate transformation can place a
3419 ! segment on the other side of the threshold again
3420 !
3421 !-----------------------------------------------------------------------
3422 
3423  if (lthresh) then
3424  if (intrsct_lat > north_thresh .or. intrsct_lat < south_thresh)
3425  & lthresh = .false.
3426  else if (beglat > zero .and. intrsct_lat < north_thresh) then
3427  mat4 = endlat - begseg(1)
3428  mat3 = endlon - begseg(2)
3429  if (mat3 > pi) mat3 = mat3 - pi2
3430  if (mat3 < -pi) mat3 = mat3 + pi2
3431 ! intrsct_lat = north_thresh - tiny
3432  intrsct_lat = north_thresh
3433  s1 = (north_thresh - begseg(1))/mat4
3434  intrsct_lon = begseg(2) + s1*mat3
3435  lthresh = .true.
3436  else if (beglat < zero .and. intrsct_lat > south_thresh) then
3437  mat4 = endlat - begseg(1)
3438  mat3 = endlon - begseg(2)
3439  if (mat3 > pi) mat3 = mat3 - pi2
3440  if (mat3 < -pi) mat3 = mat3 + pi2
3441 ! intrsct_lat = south_thresh + tiny
3442  intrsct_lat = south_thresh
3443  s1 = (south_thresh - begseg(1))/mat4
3444  intrsct_lon = begseg(2) + s1*mat3
3445  lthresh = .true.
3446  endif
3447 
3448 
3449 !-----------------------------------------------------------------------
3450 
3451  end subroutine pole_intersection
3452 
3453 !***********************************************************************
3454 
3455 
3456 
3457  subroutine line_integral(phi_or_theta, weights, num_wts,
3458  & in_phi1, in_phi2, theta1, theta2,
3459  & grid1_lat, grid1_lon, grid2_lat, grid2_lon)
3460 
3461 !-----------------------------------------------------------------------
3462 !
3463 ! this routine computes the line integral of the flux function
3464 ! that results in the interpolation weights. the line is defined
3465 ! by the input lat/lon of the endpoints.
3466 !
3467 !-----------------------------------------------------------------------
3468 
3469 !-----------------------------------------------------------------------
3470 !
3471 ! intent(in):
3472 !
3473 !-----------------------------------------------------------------------
3474 
3475  integer (SCRIP_i4), intent(in) ::
3476  & phi_or_theta ! Integration variable (lat or lon)
3477 
3478  integer (SCRIP_i4), intent(in) ::
3479  & num_wts ! number of weights to compute
3480 
3481  real (SCRIP_r8), intent(in) ::
3482  & in_phi1, in_phi2, ! longitude endpoints for the segment
3483  & theta1, theta2, ! latitude endpoints for the segment
3484  & grid1_lat, grid1_lon, ! reference coordinates for each
3485  & grid2_lat, grid2_lon ! grid (to ensure correct 0,2pi interv.
3486 
3487 !-----------------------------------------------------------------------
3488 !
3489 ! intent(out):
3490 !
3491 !-----------------------------------------------------------------------
3492 
3493  real (SCRIP_r8), dimension(2*num_wts), intent(out) ::
3494  & weights ! line integral contribution to weights
3495 
3496 
3497 ! write(*,*)'subroutine line_integral'
3498  if (phi_or_theta .eq. 1) then
3499  call line_integral_phi(weights, num_wts, in_phi1, in_phi2,
3500  & theta1, theta2, grid1_lat, grid1_lon,
3501  & grid2_lat, grid2_lon)
3502  else
3503  call line_integral_theta(weights, num_wts,in_phi1,in_phi2,
3504  & theta1, theta2, grid1_lat, grid1_lon,
3505  & grid2_lat, grid2_lon)
3506  endif
3507 
3508 
3509  return
3510 
3511 !-----------------------------------------------------------------------
3512 
3513  end subroutine line_integral
3514 
3515 !***********************************************************************
3516 
3517 
3518 
3519  subroutine line_integral_phi(weights, num_wts,
3520  & in_phi1, in_phi2, theta1, theta2,
3521  & grid1_lat, grid1_lon, grid2_lat, grid2_lon)
3522 
3523 !-----------------------------------------------------------------------
3524 !
3525 ! this routine computes the line integral of the flux function
3526 ! that results in the interpolation weights. the line is defined
3527 ! by the input lat/lon of the endpoints. Integration is w.r.t. lon
3528 !
3529 !-----------------------------------------------------------------------
3530 
3531 !-----------------------------------------------------------------------
3532 !
3533 ! intent(in):
3534 !
3535 !-----------------------------------------------------------------------
3536 
3537  integer (SCRIP_i4), intent(in) ::
3538  & num_wts ! number of weights to compute
3539 
3540  real (SCRIP_r8), intent(in) ::
3541  & in_phi1, in_phi2, ! longitude endpoints for the segment
3542  & theta1, theta2, ! latitude endpoints for the segment
3543  & grid1_lat, grid1_lon, ! reference coordinates for each
3544  & grid2_lat, grid2_lon ! grid (to ensure correct 0,2pi interv.
3545 
3546 !-----------------------------------------------------------------------
3547 !
3548 ! intent(out):
3549 !
3550 !-----------------------------------------------------------------------
3551 
3552  real (SCRIP_r8), dimension(2*num_wts), intent(out) ::
3553  & weights ! line integral contribution to weights
3554 
3555 !-----------------------------------------------------------------------
3556 !
3557 ! local variables
3558 !
3559 !-----------------------------------------------------------------------
3560 
3561  real (SCRIP_r8) :: dphi, sinth1, sinth2, costh1, costh2, fac,
3562  & phi1, phi2
3563  real (SCRIP_r8) :: f1, f2, fint
3564 
3565 !-----------------------------------------------------------------------
3566 !
3567 ! weights for the general case based on a trapezoidal approx to
3568 ! the integrals.
3569 !
3570 !-----------------------------------------------------------------------
3571 
3572 
3573 ! write(*,*)'subroutine line_integral_phi'
3574 
3575  sinth1 = sin(theta1)
3576  sinth2 = sin(theta2)
3577  costh1 = cos(theta1)
3578  costh2 = cos(theta2)
3579 
3580  dphi = in_phi1 - in_phi2
3581  if (dphi > pi) then
3582  dphi = dphi - pi2
3583  else if (dphi < -pi) then
3584  dphi = dphi + pi2
3585  endif
3586  dphi = half*dphi
3587 
3588 !-----------------------------------------------------------------------
3589 !
3590 ! the first weight is the area overlap integral. the second and
3591 ! fourth are second-order latitude gradient weights.
3592 !
3593 !-----------------------------------------------------------------------
3594 
3595  weights( 1) = dphi*(sinth1 + sinth2)
3596  write(401,*)weights(1),' % A'
3597  weights(num_wts+1) = dphi*(sinth1 + sinth2)
3598  weights( 2) = dphi*(costh1 + costh2 + (theta1*sinth1 +
3599  & theta2*sinth2))
3600  weights(num_wts+2) = dphi*(costh1 + costh2 + (theta1*sinth1 +
3601  & theta2*sinth2))
3602 
3603 !-----------------------------------------------------------------------
3604 !
3605 ! the third and fifth weights are for the second-order phi gradient
3606 ! component. must be careful of longitude range.
3607 !
3608 !-----------------------------------------------------------------------
3609 
3610  f1 = half*(costh1*sinth1 + theta1)
3611  f2 = half*(costh2*sinth2 + theta2)
3612 
3613  phi1 = in_phi1 - grid1_lon
3614  if (phi1 > pi) then
3615  phi1 = phi1 - pi2
3616  else if (phi1 < -pi) then
3617  phi1 = phi1 + pi2
3618  endif
3619 
3620  phi2 = in_phi2 - grid1_lon
3621  if (phi2 > pi) then
3622  phi2 = phi2 - pi2
3623  else if (phi2 < -pi) then
3624  phi2 = phi2 + pi2
3625  endif
3626 
3627  if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3628  weights(3) = dphi*(phi1*f1 + phi2*f2)
3629  else
3630  if (phi1 > zero) then
3631  fac = pi
3632  else
3633  fac = -pi
3634  endif
3635  fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
3636  weights(3) = half*phi1*(phi1-fac)*f1 -
3637  & half*phi2*(phi2+fac)*f2 +
3638  & half*fac*(phi1+phi2)*fint
3639  endif
3640 
3641  phi1 = in_phi1 - grid2_lon
3642  if (phi1 > pi) then
3643  phi1 = phi1 - pi2
3644  else if (phi1 < -pi) then
3645  phi1 = phi1 + pi2
3646  endif
3647 
3648  phi2 = in_phi2 - grid2_lon
3649  if (phi2 > pi) then
3650  phi2 = phi2 - pi2
3651  else if (phi2 < -pi) then
3652  phi2 = phi2 + pi2
3653  endif
3654 
3655  if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3656  weights(num_wts+3) = dphi*(phi1*f1 + phi2*f2)
3657  else
3658  if (phi1 > zero) then
3659  fac = pi
3660  else
3661  fac = -pi
3662  endif
3663  fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
3664  weights(num_wts+3) = half*phi1*(phi1-fac)*f1 -
3665  & half*phi2*(phi2+fac)*f2 +
3666  & half*fac*(phi1+phi2)*fint
3667  endif
3668 
3669 !-----------------------------------------------------------------------
3670 
3671  end subroutine line_integral_phi
3672 
3673 !***********************************************************************
3674 
3675 
3676 
3677 !***********************************************************************
3678 
3679  subroutine line_integral_theta(weights, num_wts,
3680  & in_phi1, in_phi2, theta1, theta2,
3681  & grid1_lat, grid1_lon, grid2_lat, grid2_lon)
3682 
3683 !-----------------------------------------------------------------------
3684 !
3685 ! this routine computes the line integral of the flux function
3686 ! that results in the interpolation weights. the line is defined
3687 ! by the input lat/lon of the endpoints. Integration is w.r.t. lat
3688 !
3689 ! Needed to use Simpson rule for this integration to get lower errors
3690 !-----------------------------------------------------------------------
3691 
3692 !-----------------------------------------------------------------------
3693 !
3694 ! intent(in):
3695 !
3696 !-----------------------------------------------------------------------
3697 
3698  integer (SCRIP_i4), intent(in) ::
3699  & num_wts ! number of weights to compute
3700 
3701  real (SCRIP_r8), intent(in) ::
3702  & in_phi1, in_phi2, ! longitude endpoints for the segment
3703  & theta1, theta2, ! latitude endpoints for the segment
3704  & grid1_lat, grid1_lon, ! reference coordinates for each
3705  & grid2_lat, grid2_lon ! grid (to ensure correct 0,2pi interv.
3706 
3707 !-----------------------------------------------------------------------
3708 !
3709 ! intent(out):
3710 !
3711 !-----------------------------------------------------------------------
3712 
3713  real (SCRIP_r8), dimension(2*num_wts), intent(out) ::
3714  & weights ! line integral contribution to weights
3715 
3716 !-----------------------------------------------------------------------
3717 !
3718 ! local variables
3719 !
3720 !-----------------------------------------------------------------------
3721 
3722  real (SCRIP_r8) :: dtheta, dtheta2, costh1, costh2, costhpi,
3723  & phi1, phi2, theta_pi, f1, f2, fpi,
3724  & fm, costhm, part1, part2
3725 
3726 !-----------------------------------------------------------------------
3727 !
3728 ! weights for the general case based on a trapezoidal approx to
3729 ! the integrals.
3730 !
3731 !-----------------------------------------------------------------------
3732 
3733  costh1 = cos(theta1)
3734  costh2 = cos(theta2)
3735  costhm = cos(half*(theta1+theta2))
3736 
3737  dtheta = theta2 - theta1
3738  dtheta2 = half*dtheta
3739 
3740 ! write(*,*)' subroutine line_integral_theta'
3741 
3742 
3743 !-----------------------------------------------------------------------
3744 !
3745 ! Need to account for double value of longitude in calculations of
3746 ! all the weights. First we transform all the phis to be relative
3747 ! to the grid center This takes care of a good number of cases where
3748 ! the the phis span the periodic boundary in the longitudinal
3749 ! direction. If we still have a line that spans the periodic
3750 ! boundary then we have to integrate along the line in two parts -
3751 ! from point 1 to the periodic boundary and from the periodic
3752 ! boundary to the second point
3753 !
3754 ! Example: Consider a line which has points at phi1 = -100 and phi2
3755 ! = 100 degrees and say the grid center is at phi_c = 0
3756 ! degrees. Then phi1-phi_c > -180 and phi2-phi_c < 180. But
3757 ! phi2-phi1 > 180.
3758 !
3759 ! *********************************************!!!!!!!!!!!
3760 ! If we are doing the second step anyway, why are we normalizing the
3761 ! coordinates with respect to the grid centers?
3762 !
3763 ! We need it particularly in this integration because phi figures
3764 ! explicitly in the expressions - so if a cell straddles the 0,2pi
3765 ! boundary, we integrate some edges with phi values close to zero
3766 ! and others with phi values close to 2pi leading to errors
3767 ! *********************************************!!!!!!!!!!!
3768 !
3769 !-----------------------------------------------------------------------
3770 
3771  phi1 = in_phi1 - grid1_lon
3772  if (phi1 > pi) then
3773  phi1 = phi1 - pi2
3774  else if (phi1 < -pi) then
3775  phi1 = phi1 + pi2
3776  endif
3777 
3778  phi2 = in_phi2 - grid1_lon
3779  if (phi2 > pi) then
3780  phi2 = phi2 - pi2
3781  else if (phi2 < -pi) then
3782  phi2 = phi2 + pi2
3783  endif
3784 
3785  f1 = phi1*costh1
3786  f2 = phi2*costh2
3787 
3788  if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3789 
3790  fm = half*(phi1+phi2)*costhm
3791 
3792  weights(1) = dtheta*(f1 + 4*fm + f2)/6.0
3793 ! write(401,*)weights(1),' % A'
3794 
3795  weights(2) = dtheta2*(theta1*f1 + theta2*f2)
3796 
3797  weights(3) = half*dtheta2*(f1*f1 + f2*f2)
3798 
3799  else
3800  if (phi1 > zero) then ! Means phi2-phi1 < -pi
3801 
3802 ! theta at phi = pi
3803  theta_pi = theta1 + (pi - phi1)*dtheta/(phi2 + pi2 - phi1)
3804 ! print *, ''
3805 ! print *, 'phi1',phi1,' phi2',phi2
3806 ! print *, 'theta1',theta1,' theta2',theta2
3807 ! print *, 'theta_pi',theta_pi
3808 
3809  costhpi = cos(theta_pi)
3810  fpi = pi*costhpi
3811 
3812  fm = half*(phi1+pi)*cos(half*(theta1+theta_pi))
3813  part1 = (theta_pi - theta1)*(f1 + 4*fm + fpi)/6.0
3814 
3815  fm = half*(phi2-pi)*cos(half*(theta1+theta_pi))
3816  part2 = 0.5*(theta2 - theta_pi)*(-fpi + 4*fm + f2)/6.0
3817 
3818  weights(1) = part1 + part2
3819 ! write(401,*)weights(1),' % B'
3820 
3821  part1 = 0.5*(theta_pi - theta1)*(theta1*f1 + theta_pi*fpi)
3822  part2 = 0.5*(theta2 - theta_pi)*(-theta_pi*fpi + theta2*f2)
3823  weights(2) = part1 + part2
3824 
3825 
3826  else ! Means phi2-phi1 > pi
3827 
3828 ! theta at phi = -pi
3829  theta_pi = theta1 + (-pi - phi1)*dtheta/(phi2 - pi2 - phi1)
3830 ! print *, ''
3831 ! print *, 'phi1',phi1,' phi2',phi2
3832 ! print *, 'theta1',theta1,' theta2',theta2
3833 ! print *, 'theta_pi',theta_pi
3834 
3835  costhpi = cos(theta_pi)
3836  fpi = pi*costhpi
3837 
3838  fm = half*(phi1-pi)*cos(half*(theta1+theta_pi))
3839  part1 = 0.5*(theta_pi - theta1)*(f1 + 4*fm - fpi)/6.0
3840 
3841  fm = half*(pi+phi2)*cos(half*(theta2+theta_pi))
3842  part2 = 0.5*(theta2 - theta_pi)*(fpi + 4*fm + f2)/6.0
3843  weights(1) = part1 + part2
3844 ! write(401,*)weights(1),' % C'
3845 
3846  part1 = 0.5*(theta_pi - theta1)*(theta1*f1 - theta_pi*fpi)
3847  part2 = 0.5*(theta2 - theta_pi)*(theta_pi*fpi + theta2*f2)
3848  weights(2) = part1 + part2
3849 
3850 
3851  endif
3852 
3853  part1 = 0.25*(theta_pi - theta1)*(f1*f1 + fpi*fpi)
3854  part2 = 0.25*(theta2 - theta_pi)*(fpi*fpi + f2*f2)
3855  weights(3) = part1 + part2
3856 
3857  endif
3858 
3859 
3860  phi1 = in_phi1 - grid2_lon
3861  if (phi1 > pi) then
3862  phi1 = phi1 - pi2
3863  else if (phi1 < -pi) then
3864  phi1 = phi1 + pi2
3865  endif
3866 
3867  phi2 = in_phi2 - grid2_lon
3868  if (phi2 > pi) then
3869  phi2 = phi2 - pi2
3870  else if (phi2 < -pi) then
3871  phi2 = phi2 + pi2
3872  endif
3873 
3874 
3875  f1 = phi1*costh1
3876  f2 = phi2*costh2
3877 
3878  if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3879 
3880  fm = half*(phi1+phi2)*costhm
3881 
3882  weights(num_wts+1) = dtheta2*(f1 + f2)
3883 
3884  weights(num_wts+2) = dtheta2*(theta1*f1 + theta2*f2)
3885 
3886  weights(num_wts+3) = half*dtheta2*(f1*f1 + f2*f2)
3887 
3888  else
3889  if (phi1 > zero) then
3890 
3891  theta_pi = theta1 + (pi - phi1)*dtheta/(phi2 + pi2 - phi1)
3892 ! print *, ''
3893 ! print *, 'phi1',phi1,' phi2',phi2
3894 ! print *, 'theta1',theta1,' theta2',theta2
3895 ! print *, 'theta_pi',theta_pi
3896 
3897  costhpi = cos(theta_pi)
3898  fpi = pi*costhpi
3899 
3900  fm = half*(phi1+pi)*cos(half*(theta1+theta_pi))
3901  part1 = (theta_pi - theta1)*(f1 + 4*fm + fpi)/6.0
3902 
3903  fm = half*(-pi+phi2)*cos(half*(theta2+theta_pi))
3904  part2 = (theta2 - theta_pi)*(-fpi + 4*fm + f2)/6.0
3905  weights(num_wts+1) = part1 + part2
3906 
3907  part1 = 0.5*(theta_pi - theta1)*(theta1*f1 + theta_pi*fpi)
3908  part2 = 0.5*(theta2 - theta_pi)*(-theta_pi*fpi + theta2*f2)
3909  weights(num_wts+2) = part1 + part2
3910 
3911 
3912  else
3913 
3914  theta_pi = theta1 + (-pi - phi1)*dtheta/(phi2 - pi2 - phi1)
3915 ! print *, ''
3916 ! print *, 'phi1',phi1,' phi2',phi2
3917 ! print *, 'theta1',theta1,' theta2',theta2
3918 ! print *, 'theta_pi',theta_pi
3919 
3920  costhpi = cos(theta_pi)
3921  fpi = pi*costhpi
3922 
3923  fm = half*(phi1-pi)*cos(half*(theta1+theta_pi))
3924  part1 = (theta_pi - theta1)*(f1 +4*fm - fpi)/6.0
3925 
3926  fm = half*(phi2+pi)*cos(half*(theta2+theta_pi))
3927  part2 = 0.5*(theta2 - theta_pi)*(fpi + 4*fm + f2)/6.0
3928  weights(num_wts+1) = part1 + part2
3929 
3930  part1 = 0.5*(theta_pi - theta1)*(theta1*f1 - theta_pi*fpi)
3931  part2 = 0.5*(theta2 - theta_pi)*(theta_pi*fpi + theta2*f2)
3932  weights(num_wts+2) = part1 + part2
3933 
3934  endif
3935 
3936  part1 = 0.25*(theta_pi - theta1)*(f1*f1 + fpi*fpi)
3937  part2 = 0.25*(theta2 - theta_pi)*(fpi*fpi + f2*f2)
3938  weights(num_wts+3) = part1 + part2
3939 
3940  endif
3941 
3942 !-----------------------------------------------------------------------
3943 
3944  end subroutine line_integral_theta
3945 
3946 !***********************************************************************
3947 
3948 
3949 
3950  subroutine store_link_cnsrv(add1, add2, weights)
3951 
3952 !-----------------------------------------------------------------------
3953 !
3954 ! this routine stores the address and weight for this link in
3955 ! the appropriate address and weight arrays and resizes those
3956 ! arrays if necessary.
3957 !
3958 !-----------------------------------------------------------------------
3959 
3960 !-----------------------------------------------------------------------
3961 !
3962 ! input variables
3963 !
3964 !-----------------------------------------------------------------------
3965 
3966  integer (SCRIP_i4), intent(in) ::
3967  & add1, ! address on grid1
3968  & add2 ! address on grid2
3969 
3970  real (SCRIP_r8), dimension(:), intent(in) ::
3971  & weights ! array of remapping weights for this link
3972 
3973 !-----------------------------------------------------------------------
3974 !
3975 ! local variables
3976 !
3977 !-----------------------------------------------------------------------
3978 
3979  integer (SCRIP_i4) :: nlink, min_link, max_link ! link index
3980 
3981  logical (SCRIP_logical) :: found
3982 
3983 
3984 !-----------------------------------------------------------------------
3985 !
3986 ! if all weights are zero, do not bother storing the link
3987 !
3988 !-----------------------------------------------------------------------
3989 
3990  if (all(weights == zero)) return
3991 
3992 !-----------------------------------------------------------------------
3993 !
3994 ! restrict the range of links to search for existing links
3995 !
3996 !-----------------------------------------------------------------------
3997 
3998 C$OMP CRITICAL(block5)
3999 ! first_call should be within critical block or else multiple
4000 ! threads will see it as true the first time around
4001 
4002  if (first_call_store_link_cnsrv) then
4003  allocate(link_add1(2,grid1_size), link_add2(2,grid2_size))
4004  link_add1 = 0
4005  link_add2 = 0
4006  first_call_store_link_cnsrv = .false.
4007  min_link = 1
4008  max_link = 0
4009  else
4010  min_link = min(link_add1(1,add1),link_add2(1,add2))
4011  max_link = max(link_add1(2,add1),link_add2(2,add2))
4012  if (min_link == 0) then
4013  min_link = 1
4014  max_link = 0
4015  endif
4016  endif
4017 C$OMP END CRITICAL(block5)
4018 
4019 !-----------------------------------------------------------------------
4020 !
4021 ! if the link already exists, add the weight to the current weight
4022 ! arrays
4023 !
4024 !-----------------------------------------------------------------------
4025 
4026  found = .false.
4027 
4028  do nlink=min_link,max_link
4029  if (add1 == grid1_add_map1(nlink)) then
4030  if (add2 == grid2_add_map1(nlink)) then
4031 
4032 C$OMP CRITICAL(block3a)
4033  wts_map1(:,nlink) = wts_map1(:,nlink) + weights(1:num_wts)
4034  if (num_maps == 2) then
4035  wts_map2(:,nlink) = wts_map2(:,nlink) +
4036  & weights(num_wts+1:2*num_wts)
4037  endif
4038 C$OMP END CRITICAL(block3a)
4039  found = .true.
4040  exit
4041 
4042  endif
4043  endif
4044  end do
4045 
4046 
4047  if (found) return
4048 
4049 !-----------------------------------------------------------------------
4050 !
4051 ! if the link does not yet exist, increment number of links and
4052 ! check to see if remap arrays need to be increased to accomodate
4053 ! the new link. then store the link.
4054 !
4055 !-----------------------------------------------------------------------
4056 
4057 C$OMP CRITICAL(block6)
4058 
4059  num_links_map1 = num_links_map1 + 1
4060  if (num_links_map1 > max_links_map1)
4061  & call resize_remap_vars(1,resize_increment)
4062 
4063  grid1_add_map1(num_links_map1) = add1
4064  grid2_add_map1(num_links_map1) = add2
4065  wts_map1(:,num_links_map1) = weights(1:num_wts)
4066 
4067  if (num_maps > 1) then
4068  num_links_map2 = num_links_map2 + 1
4069  if (num_links_map2 > max_links_map2)
4070  & call resize_remap_vars(2,resize_increment)
4071 
4072  grid1_add_map2(num_links_map2) = add1
4073  grid2_add_map2(num_links_map2) = add2
4074  wts_map2(:,num_links_map2) = weights(num_wts+1:2*num_wts)
4075  endif
4076 
4077  if (link_add1(1,add1) == 0) link_add1(1,add1) = num_links_map1
4078  if (link_add2(1,add2) == 0) link_add2(1,add2) = num_links_map1
4079  link_add1(2,add1) = num_links_map1
4080  link_add2(2,add2) = num_links_map1
4081 
4082 C$OMP END CRITICAL(block6)
4083 
4084 !-----------------------------------------------------------------------
4085 
4086  end subroutine store_link_cnsrv
4087 
4088 !***********************************************************************
4089 
4090 
4091 
4092 
4093  subroutine locate_segstart(cell_grid_num, cell,
4094  & beglat, beglon, endlat, endlon, offset,
4095  & srch_grid_num, cont_cell, lboundary, edgeid)
4096 
4097 !-----------------------------------------------------------------------
4098 !
4099 ! Find the cell containing the given point
4100 !
4101 !-----------------------------------------------------------------------
4102 
4103  implicit none
4104 
4105 !-----------------------------------------------------------------------
4106 !
4107 ! intent(in):
4108 !
4109 !-----------------------------------------------------------------------
4110 
4111  real (SCRIP_r8), intent(in) ::
4112  & beglat, beglon, ! beginning and end points of segment
4113  & endlat, endlon ! on which the point to be located lies
4114 
4115  real (SCRIP_r8), intent(in) ::
4116  & offset ! Offset to calculate the search point
4117 
4118  integer (SCRIP_i4), intent(in) ::
4119  & cell, ! Cell from which point originates
4120  ! Point will be on boundary of orig_cell
4121  & cell_grid_num ! Index of grid to which cell belongs
4122 
4123  integer (SCRIP_i4), intent(in) ::
4124  & srch_grid_num ! num indicating if we are locating a
4125  ! grid1 point in a cell of grid2 (num=2)
4126  ! or a grid2 point in a cell of grid1
4127  ! (num=1)
4128 
4129 !-----------------------------------------------------------------------
4130 !
4131 ! intent(out):
4132 !
4133 !-----------------------------------------------------------------------
4134 
4135  integer (SCRIP_i4), intent(out) ::
4136  & cont_cell ! grid cell containing this point
4137 
4138  logical (SCRIP_logical), intent(out) ::
4139  & lboundary ! flag points that lie on the boundary
4140  ! of the cell
4141 
4142  integer (SCRIP_i4), intent(out) ::
4143  & edgeid ! if point is on boundary, which local
4144  ! edge is it on? (0 otherwise)
4145 
4146 !-----------------------------------------------------------------------
4147 !
4148 ! local variables
4149 !
4150 !-----------------------------------------------------------------------
4151 
4152  integer (SCRIP_i4) :: i, j, k, n, ic
4153  integer (SCRIP_i4) :: whichpole, srch_cell_add,
4154  & grid1_add, grid2_add, min_add, max_add
4155 
4156  real (SCRIP_r8), dimension(:), allocatable ::
4157  & cell_corner_x, cell_corner_y
4158 
4159  logical (SCRIP_logical) :: inpoly, latlon
4160 
4161  real (SCRIP_r8) ::
4162  & vec1_x, vec1_y, vec1_lenx, vec1_lat, vec1_lon, vec1_len,
4163  & begx, begy, endx, endy, ptx, pty, rns, pi4, ptlat, ptlon,
4164  & lat, lon, cell_center_x, cell_center_y
4165 
4166 !-----------------------------------------------------------------------
4167 !
4168 ! initialize defaults, flags, etc.
4169 !
4170 !-----------------------------------------------------------------------
4171 
4172  lboundary = .false.
4173  edgeid = 0
4174  cont_cell = 0
4175 
4176 
4177  if (cell /= last_cell_locate_segstart .or.
4178  & cell_grid_num /= last_cell_grid_num_locate_segstart
4179  & .or. srch_grid_num /= last_srch_grid_num_locate_segstart)
4180  & then
4181 
4182  last_cell_locate_segstart = cell
4183  last_cell_grid_num_locate_segstart = cell_grid_num
4184  last_srch_grid_num_locate_segstart = srch_grid_num
4185 
4186  if (first_call_locate_segstart) then
4187  first_call_locate_segstart = .false.
4188  last_cell_locate_segstart = 0
4189  last_cell_grid_num_locate_segstart = 0
4190  last_srch_grid_num_locate_segstart = 0
4191  num_srch_cells_locate_segstart = 0
4192  else
4193  if (num_srch_cells_locate_segstart .gt. 0) then
4194  deallocate(srch_add_locate_segstart,
4195  & srch_corner_lat_locate_segstart,
4196  & srch_corner_lon_locate_segstart,
4197  & srch_center_lat_locate_segstart,
4198  & srch_center_lon_locate_segstart)
4199  endif
4200  endif
4201 
4202  call get_srch_cells(cell, cell_grid_num, srch_grid_num,
4203  & num_srch_cells_locate_segstart, srch_add_locate_segstart,
4204  & srch_corners_locate_segstart,
4205  & srch_corner_lat_locate_segstart,
4206  & srch_corner_lon_locate_segstart,
4207  & srch_center_lat_locate_segstart,
4208  & srch_center_lon_locate_segstart)
4209 
4210  endif
4211 
4212  if (num_srch_cells_locate_segstart == 0) return
4213 
4214 
4215  do ic=1,num_srch_cells_locate_segstart
4216 
4217  srch_cell_add = srch_add_locate_segstart(ic)
4218 
4219 
4220 
4221  !**** CAN WE ACCOMPLISH THE FOLLOWING THROUGH A SUBROUTINE
4222  !**** CALLED SEGSTART_INCELL ??
4223 
4224 
4225  !*** IF POINT IS IN POLAR REGION, CHECK IN A TRANSFORMED SPACE
4226  !*** HOWEVER, POINTS THAT ARE PRACTICALLY AT THE POLE CANNOT
4227  !*** BE CORRECTLY LOCATED THIS WAY BECAUSE THE POLE IS A SINGULARITY
4228  !*** AND CONTAINMENT IN ANY CELL INCIDENT ON THE POLE WILL GIVE US A
4229  !*** POSITIVE ANSWER. FOR THESE POINTS REVERT TO THE LATLON SPACE
4230  !***
4231 
4232 
4233 
4234  vec1_lat = endlat-beglat
4235  vec1_lon = endlon-beglon
4236  if (vec1_lon > pi) then
4237  vec1_lon = vec1_lon - pi2
4238  else if (vec1_lon < -pi) then
4239  vec1_lon = vec1_lon + pi2
4240  endif
4241  vec1_len = sqrt(vec1_lat*vec1_lat+vec1_lon*vec1_lon)
4242  vec1_lat = vec1_lat/vec1_len
4243  vec1_lon = vec1_lon/vec1_len
4244 
4245  ptlat = beglat + offset*vec1_lat
4246  ptlon = beglon + offset*vec1_lon
4247 
4248 
4249  if ((ptlat .gt. north_thresh .and. abs(ptlat-pih) .ge. 0.001)
4250  & .or.
4251  & (ptlat .lt. south_thresh .and. abs(ptlat+pih) .ge. 0.001))
4252  & then
4253 
4254  if (ptlat > zero) then
4255  pi4 = quart*pi
4256  rns = one
4257  else
4258  pi4 = -quart*pi
4259  rns = -one
4260  endif
4261 
4262 
4263  begx = rns*two*sin(pi4 - half*beglat)*cos(beglon)
4264  begy = two*sin(pi4 - half*beglat)*sin(beglon)
4265  endx = rns*two*sin(pi4 - half*endlat)*cos(endlon)
4266  endy = two*sin(pi4 - half*endlat)*sin(endlon)
4267 
4268  vec1_x = endx-begx
4269  vec1_y = endy-begy
4270 
4271  vec1_lenx = sqrt(vec1_x*vec1_x + vec1_y*vec1_y)
4272  vec1_x = vec1_x/vec1_lenx
4273  vec1_y = vec1_y/vec1_lenx
4274 
4275 
4276  !*** Must calculate ptx and pty as an offset on straight
4277  !*** line in polar space rather than calculating it on a
4278  !*** straight line in latlon space an offset point in latlon
4279  !*** space will be off the straight line in polar space
4280 
4281  ptx = begx + offset*vec1_x
4282  pty = begy + offset*vec1_y
4283 
4284  latlon = .false.
4285 
4286  ! Since we want greater fidelity for locating the points
4287  ! we send in the mid-points of the polygon edges too
4288  ! BUT THAT MAKES THE POLYGON NON-CONVEX SOMETIMES AND
4289  ! THE CROSS-PRODUCT CHECK FAILS. SO USE CODE TO CHECK GENERAL
4290  ! POLYGONS
4291 
4292 
4293  allocate(cell_corner_x(npseg*srch_corners_locate_segstart),
4294  & cell_corner_y(npseg*srch_corners_locate_segstart))
4295 
4296 
4297  k = 0
4298  do i = srch_corners_locate_segstart, 1, -1
4299  k = k+1
4300  lat = srch_corner_lat_locate_segstart(i,ic)
4301  lon = srch_corner_lon_locate_segstart(i,ic)
4302  cell_corner_x(k) = rns*two*sin(pi4-half*lat)*cos(lon)
4303  cell_corner_y(k) = two*sin(pi4-half*lat)*sin(lon)
4304 
4305  j = i-1
4306  if (j .eq. 0) j = srch_corners_locate_segstart ! how do
4307  ! we do (j-1+n)%n in F90?
4308 
4309  vec1_lat = srch_corner_lat_locate_segstart(j,ic)
4310  & -srch_corner_lat_locate_segstart(i,ic)
4311  vec1_lon = srch_corner_lon_locate_segstart(j,ic)
4312  & -srch_corner_lon_locate_segstart(i,ic)
4313  if (vec1_lon > pi) then
4314  vec1_lon = vec1_lon - pi2
4315  else if (vec1_lon < -pi) then
4316  vec1_lon = vec1_lon + pi2
4317  endif
4318 
4319  do j = 1, npseg-1
4320  k = k+1
4321  lat = srch_corner_lat_locate_segstart(i,ic)
4322  & + j*vec1_lat/npseg
4323  lon = srch_corner_lon_locate_segstart(i,ic)
4324  & + j*vec1_lon/npseg
4325  cell_corner_x(k) = rns*two*sin(pi4-half*lat)*cos(lon)
4326  cell_corner_y(k) = two*sin(pi4-half*lat)*sin(lon)
4327  enddo
4328  enddo
4329 
4330 
4331  call ptinpolygen2(ptx, pty, k, cell_corner_x,
4332  & cell_corner_y, latlon, inpoly, lboundary, edgeid)
4333 
4334  if (lboundary) then
4335  edgeid = (edgeid-1)/npseg + 1 ! convert from index in
4336  ! multi-segmented to regular cell
4337  endif
4338 
4339  deallocate(cell_corner_x, cell_corner_y)
4340 
4341  else
4342 
4343  latlon = .true.
4344 
4345  whichpole = 0
4346  if (srch_grid_num .eq. 1 .and.
4347  & srch_cell_add .eq. grid1_spole_cell) then
4348 
4349  whichpole = -1 ! S pole
4350  call ptinpolarpoly(ptlat, ptlon,
4351  & srch_corners_locate_segstart,
4352  & srch_corner_lat_locate_segstart(:,ic),
4353  & srch_corner_lon_locate_segstart(:,ic),
4354  & latlon, whichpole, inpoly, lboundary, edgeid)
4355 
4356  else if (srch_grid_num .eq. 1 .and.
4357  & srch_cell_add .eq. grid1_npole_cell) then
4358 
4359  whichpole = 1 ! N pole
4360  call ptinpolarpoly(ptlat, ptlon,
4361  & srch_corners_locate_segstart,
4362  & srch_corner_lat_locate_segstart(:,ic),
4363  & srch_corner_lon_locate_segstart(:,ic),
4364  & latlon, whichpole, inpoly, lboundary, edgeid)
4365 
4366  else if (srch_grid_num .eq. 2 .and.
4367  & srch_cell_add .eq. grid2_spole_cell) then
4368 
4369  whichpole = -1 ! S pole
4370  call ptinpolarpoly(ptlat, ptlon,
4371  & srch_corners_locate_segstart,
4372  & srch_corner_lat_locate_segstart(:,ic),
4373  & srch_corner_lon_locate_segstart(:,ic),
4374  & latlon, whichpole, inpoly, lboundary, edgeid)
4375 
4376  else if (srch_grid_num .eq. 2 .and.
4377  & srch_cell_add .eq. grid2_npole_cell) then
4378 
4379  whichpole = 1 ! N pole
4380  call ptinpolarpoly(ptlat, ptlon,
4381  & srch_corners_locate_segstart,
4382  & srch_corner_lat_locate_segstart(:,ic),
4383  & srch_corner_lon_locate_segstart(:,ic),
4384  & latlon, whichpole, inpoly, lboundary, edgeid)
4385 
4386  else
4387 
4388  !***
4389  !*** General cell
4390  !***
4391 
4392  call ptinpoly(ptlat, ptlon, srch_corners_locate_segstart,
4393  & srch_corner_lat_locate_segstart(:,ic),
4394  & srch_corner_lon_locate_segstart(:,ic),
4395  & latlon, inpoly, lboundary, edgeid)
4396 
4397  endif
4398 
4399  endif
4400 
4401  if (inpoly) then
4402  cont_cell = srch_cell_add
4403  exit
4404  endif
4405 
4406  end do
4407 
4408  return
4409 
4410 !----------------------------------------------------------------------
4411 
4412  end subroutine locate_segstart
4413 
4414 !**********************************************************************
4415 
4416 
4417 
4418 
4419 !**********************************************************************
4420 
4421  subroutine locate_point(ptlat, ptlon, cell, cell_grid_num,
4422  & srch_grid_num, cont_cell, lboundary, edgeid)
4423 
4424 !-----------------------------------------------------------------------
4425 !
4426 ! Find the cell containing the given point
4427 !
4428 !-----------------------------------------------------------------------
4429 
4430  implicit none
4431 
4432 !-----------------------------------------------------------------------
4433 !
4434 ! intent(in):
4435 !
4436 !-----------------------------------------------------------------------
4437 
4438  real (SCRIP_r8), intent(in) ::
4439  & ptlat, ptlon ! Point to locate
4440 
4441  integer (SCRIP_i4), intent(in) ::
4442  & cell, ! Cell from which point originates
4443  ! Point will be on boundary of orig_cell
4444  & cell_grid_num ! Index of grid to which cell belongs
4445 
4446  integer (SCRIP_i4), intent(in) ::
4447  & srch_grid_num ! num indicating if we are locating a
4448  ! grid1 point in a cell of grid2 (num=2)
4449  ! or a grid2 point in a cell of grid1
4450  ! (num=1)
4451 
4452 !-----------------------------------------------------------------------
4453 !
4454 ! intent(out):
4455 !
4456 !-----------------------------------------------------------------------
4457 
4458  integer (SCRIP_i4), intent(out) ::
4459  & cont_cell ! grid cell containing this point
4460 
4461  logical (SCRIP_logical), intent(out) ::
4462  & lboundary ! flag points that lie on the boundary
4463  ! of the cell
4464 
4465  integer (SCRIP_i4), intent(out) ::
4466  & edgeid ! if point is on boundary, which local
4467  ! edge is it on? (0 otherwise)
4468 
4469 !-----------------------------------------------------------------------
4470 !
4471 ! local variables
4472 !
4473 !-----------------------------------------------------------------------
4474 
4475  integer (SCRIP_i4) :: i, j, n, ic
4476  integer (SCRIP_i4) :: whichpole, srch_cell_add,
4477  & grid1_add, grid2_add, min_add, max_add,
4478  & previdx, nextidx, pcorner, corner,
4479  & ncorners, nalloc
4480 
4481  real (SCRIP_r8), dimension(:), allocatable ::
4482  & cell_corner_lat,
4483  & cell_corner_lon
4484 
4485  real (SCRIP_r8) ::
4486  & prevlon,
4487  & nextlon,
4488  & polelat,
4489  & cell_center_lat,
4490  & cell_center_lon
4491 
4492 
4493  logical (SCRIP_logical) :: inpoly, latlon
4494  logical (SCRIP_logical) :: test
4495 
4496 !-----------------------------------------------------------------------
4497 !
4498 ! initialize defaults, flags, etc.
4499 !
4500 !-----------------------------------------------------------------------
4501 
4502  lboundary = .false.
4503  edgeid = 0
4504  cont_cell = 0
4505 
4506  if (cell /= last_cell_locate_point .or. cell_grid_num /=
4507  & last_cell_grid_num_locate_point
4508  & .or. srch_grid_num /= last_srch_grid_num_locate_point) then
4509 
4510  last_cell_locate_point = cell
4511  last_cell_grid_num_locate_point = cell_grid_num
4512  last_srch_grid_num_locate_point = srch_grid_num
4513 
4514  if (first_call_locate_point) then
4515  first_call_locate_point = .false.
4516  last_cell_locate_point = 0
4517  last_cell_grid_num_locate_point = 0
4518  last_srch_grid_num_locate_point = 0
4519  num_srch_cell_locate_points = 0
4520  else
4521  if (num_srch_cell_locate_points .gt. 0) then
4522  deallocate(srch_add_locate_point,
4523  & srch_corner_lat_locate_point,
4524  & srch_corner_lon_locate_point)
4525  endif
4526  endif
4527 
4528  call get_srch_cells(cell, cell_grid_num, srch_grid_num,
4529  & num_srch_cell_locate_points, srch_add_locate_point,
4530  & srch_corners_locate_point,
4531  & srch_corner_lat_locate_point,srch_corner_lon_locate_point,
4532  & srch_center_lat_locate_point,srch_center_lon_locate_point)
4533 
4534  endif
4535 
4536  if (num_srch_cell_locate_points == 0) return
4537 
4538 
4539  ncorners = srch_corners_locate_point
4540  nalloc = ncorners+2
4541  allocate(cell_corner_lat(nalloc),
4542  & cell_corner_lon(nalloc))
4543 
4544 
4545  do ic=1,num_srch_cell_locate_points
4546 
4547  srch_cell_add = srch_add_locate_point(ic)
4548 
4549  do i = 1, ncorners
4550  cell_corner_lat(i) = srch_corner_lat_locate_point(i,ic)
4551  cell_corner_lon(i) = srch_corner_lon_locate_point(i,ic)
4552  enddo
4553 
4554  cell_center_lat = srch_center_lat_locate_point(ic)
4555  cell_center_lon = srch_center_lon_locate_point(ic)
4556 
4557 ! if ((srch_grid_num .eq. 1 .and.
4558 ! & (special_polar_cell1(srch_cell_add))) .or.
4559 ! & (srch_grid_num .eq. 2 .and.
4560 ! & (special_polar_cell2(srch_cell_add)))) then
4561 !
4562 ! Modified by MD
4563  test=.false.
4564  if (srch_grid_num .eq. 1) then
4565  if (special_polar_cell1(srch_cell_add)) then
4566  test=.true.
4567  endif
4568  else
4569  if (special_polar_cell2(srch_cell_add)) then
4570  test=.true.
4571  endif
4572  endif
4573  if (test) then
4574  call modify_polar_cell(ncorners, nalloc, cell_corner_lat,
4575  & cell_corner_lon)
4576 
4577  endif
4578 
4579  call ptincell(ptlat, ptlon, srch_cell_add, ncorners,
4580  & cell_corner_lat, cell_corner_lon,
4581  & cell_center_lat, cell_center_lon,
4582  & srch_grid_num, inpoly, lboundary, edgeid)
4583 
4584 
4585  if (inpoly) then
4586  cont_cell = srch_cell_add
4587  exit
4588  endif
4589 
4590  ncorners = srch_corners_locate_point ! reset it for other srch
4591  !cells
4592  end do
4593 
4594 !----------------------------------------------------------------------
4595 
4596  end subroutine locate_point
4597 
4598 !**********************************************************************
4599 
4600 
4601 
4602 !**********************************************************************
4603 
4604  subroutine ptincell(ptlat, ptlon, cell_add, ncorners,
4605  & cell_corner_lat, cell_corner_lon,
4606  & cell_center_lat, cell_center_lon,
4607  & cell_grid_id, inpoly, lboundary, edgeid)
4608 
4609 !----------------------------------------------------------------------
4610 
4611  implicit none
4612 
4613 !-----------------------------------------------------------------------
4614 !
4615 ! intent(in):
4616 !
4617 !-----------------------------------------------------------------------
4618 
4619  real (SCRIP_r8), intent(in) ::
4620  & ptlat, ptlon ! Point to locate
4621 
4622  integer (SCRIP_i4), intent(in) ::
4623  & cell_add ! ID of cell
4624 
4625  integer (SCRIP_i4), intent(in) ::
4626  & ncorners
4627 
4628  real (SCRIP_r8), dimension(ncorners), intent(in) ::
4629  & cell_corner_lat, cell_corner_lon
4630 
4631  real (SCRIP_r8), intent(in) ::
4632  & cell_center_lat,
4633  & cell_center_lon
4634 
4635  integer (SCRIP_i4), intent(in) ::
4636  & cell_grid_id ! num indicating if we are locating a grid1
4637  ! point in a cell of grid2 (num = 2) or
4638  ! a grid2 point in a cell of grid1 (num = 1)
4639 
4640 
4641 !-----------------------------------------------------------------------
4642 !
4643 ! intent(out):
4644 !
4645 !-----------------------------------------------------------------------
4646 
4647  logical (SCRIP_logical), intent(out) ::
4648  & inpoly ! is point in polygon?
4649 
4650  logical (SCRIP_logical), intent(out) ::
4651  & lboundary ! flag points that lie on the boundary
4652  ! of the cell
4653 
4654  integer (SCRIP_i4), intent(out) ::
4655  & edgeid ! if point is on boundary, which local
4656  ! edge is it on? (0 otherwise)
4657 !-----------------------------------------------------------------------
4658 !
4659 ! local variables
4660 !
4661 !-----------------------------------------------------------------------
4662 
4663  integer (SCRIP_i4) :: i, j, k, ic
4664  integer (SCRIP_i4) :: whichpole
4665 
4666  real (SCRIP_r8) :: rns, pi4, ptx, pty, lat, lon,
4667  & cell_center_x, cell_center_y, vec1_lat, vec1_lon
4668 
4669  logical (SCRIP_logical) ::
4670  & latlon
4671 
4672  real (kind=scrip_r8), dimension(npseg*ncorners) ::
4673  & cell_corner_x, ! x of each corner of cell
4674  & cell_corner_y ! y of each corner of cell
4675 
4676 !----------------------------------------------------------------------
4677 
4678  edgeid = 0
4679 
4680 
4681  !*** IF POINTS ARE ABOVE THE THRESHOLD, CHECK THEM IN A TRANSFORMED
4682  !*** SPACE
4683  !*** HOWEVER, POINTS THAT ARE PRACTICALLY AT THE POLE CANNOT
4684  !*** BE CORRECTLY LOCATED THIS WAY BECAUSE THE POLE IS A SINGULARITY
4685  !*** AND CONTAINMENT IN ANY CELL INCIDENT ON THE POLE WILL GIVE US A
4686  !*** POSITIVE ANSWER. FOR THESE POINTS REVERT TO THE LATLON SPACE
4687  !***
4688 
4689  if ((ptlat .gt. north_thresh .and. abs(ptlat-pih) .ge. 0.001) .or.
4690  & (ptlat .lt. south_thresh .and. abs(ptlat+pih) .ge. 0.001))
4691  & then
4692 
4693  if (ptlat > zero) then
4694  pi4 = quart*pi
4695  rns = one
4696  else
4697  pi4 = -quart*pi
4698  rns = -one
4699  endif
4700 
4701  ptx = rns*two*sin(pi4 - half*ptlat)*cos(ptlon)
4702  pty = two*sin(pi4 - half*ptlat)*sin(ptlon)
4703 
4704  latlon = .false.
4705 
4706  ! Since we want greater fidelity for locating the points
4707  ! we send in the mid-points of the polygon edges too
4708  ! BUT THAT MAKES THE POLYGON NON-CONVEX SOMETIMES AND
4709  ! THE CROSS-PRODUCT CHECK FAILS. SO USE CODE TO CHECK GENERAL
4710  ! POLYGONS
4711 
4712 
4713  k = 0
4714  do i = ncorners, 1, -1
4715  k = k+1
4716  lat = cell_corner_lat(i)
4717  lon = cell_corner_lon(i)
4718  cell_corner_x(k) = rns*two*sin(pi4-half*lat)*cos(lon)
4719  cell_corner_y(k) = two*sin(pi4-half*lat)*sin(lon)
4720 
4721  j = i-1
4722  if (j .eq. 0) j = ncorners ! how do we do (j-1+n)%n in F90?
4723 
4724  vec1_lat = cell_corner_lat(j)-cell_corner_lat(i)
4725  vec1_lon = cell_corner_lon(j)-cell_corner_lon(i)
4726  if (vec1_lon > pi) then
4727  vec1_lon = vec1_lon - pi2
4728  else if (vec1_lon < -pi) then
4729  vec1_lon = vec1_lon + pi2
4730  endif
4731 
4732  do j = 1, npseg-1
4733  k = k+1
4734  lat = cell_corner_lat(i) + j*vec1_lat/npseg
4735  lon = cell_corner_lon(i) + j*vec1_lon/npseg
4736  cell_corner_x(k) = rns*two*sin(pi4-half*lat)*cos(lon)
4737  cell_corner_y(k) = two*sin(pi4-half*lat)*sin(lon)
4738  enddo
4739  enddo
4740 
4741  !*** cell is so non-convex that no feasible center exists
4742  !*** we have to fall back on a different algorithm
4743 
4744  call ptinpolygen2(ptx, pty, k, cell_corner_x,
4745  & cell_corner_y, latlon, inpoly, lboundary, edgeid)
4746 
4747  if (lboundary) then
4748  edgeid = (edgeid-1)/npseg + 1 ! convert from index in
4749  ! multi-segmented cell to
4750  ! regular cell
4751  endif
4752  else
4753 
4754  latlon = .true.
4755 
4756  whichpole = 0
4757  if (cell_grid_id .eq. 1 .and.
4758  & cell_add .eq. grid1_spole_cell) then
4759 
4760  whichpole = -1 ! S pole
4761  call ptinpolarpoly(ptlat, ptlon, ncorners,
4762  & cell_corner_lat, cell_corner_lon,
4763  & latlon, whichpole, inpoly, lboundary, edgeid)
4764 
4765  else if (cell_grid_id .eq. 1 .and.
4766  & cell_add .eq. grid1_npole_cell) then
4767 
4768  whichpole = 1 ! N pole
4769  call ptinpolarpoly(ptlat, ptlon, ncorners,
4770  & cell_corner_lat, cell_corner_lon,
4771  & latlon, whichpole, inpoly, lboundary, edgeid)
4772 
4773  else if (cell_grid_id .eq. 2 .and.
4774  & cell_add .eq. grid2_spole_cell) then
4775 
4776  whichpole = -1 ! S pole
4777  call ptinpolarpoly(ptlat, ptlon, ncorners,
4778  & cell_corner_lat, cell_corner_lon,
4779  & latlon, whichpole, inpoly, lboundary, edgeid)
4780 
4781  else if (cell_grid_id .eq. 2 .and.
4782  & cell_add .eq. grid2_npole_cell) then
4783 
4784  whichpole = 1 ! N pole
4785  call ptinpolarpoly(ptlat, ptlon, ncorners,
4786  & cell_corner_lat, cell_corner_lon,
4787  & latlon, whichpole, inpoly, lboundary, edgeid)
4788 
4789  else
4790 
4791  !***
4792  !*** General cell
4793  !***
4794 
4795  call ptinpoly(ptlat, ptlon, ncorners,
4796  & cell_corner_lat, cell_corner_lon,
4797  & latlon, inpoly, lboundary, edgeid)
4798 
4799  endif
4800 
4801  endif
4802 
4803  return
4804 
4805 !----------------------------------------------------------------------
4806 
4807  end subroutine ptincell
4808 
4809 !**********************************************************************
4810 
4811 !**********************************************************************
4812 
4813  subroutine ptinpoly(ptx, pty, ncorners, cell_corner_x,
4814  & cell_corner_y, latlon, inpoly, lboundary, edgeid)
4815 
4816 !----------------------------------------------------------------------
4817 !
4818 ! Check if point is in (convex) polygonal cell
4819 !
4820 !----------------------------------------------------------------------
4821 
4822  implicit none
4823 
4824 !----------------------------------------------------------------------
4825 !
4826 ! Input arguments
4827 !
4828 !----------------------------------------------------------------------
4829 
4830  real (SCRIP_r8), intent(in) ::
4831  & ptx, pty ! Point to check
4832 
4833  integer (SCRIP_i4), intent(in) ::
4834  & ncorners ! Number of polygon corners
4835 
4836  real (SCRIP_r8), dimension(ncorners), intent(in) ::
4837  & cell_corner_x, ! Coordinates of cell corners
4838  & cell_corner_y ! Could be x-y or lat-lon or ...
4839 
4840  logical (SCRIP_logical), intent(in) ::
4841  & latlon ! Are coordinates in latlon space?
4842 
4843 !----------------------------------------------------------------------
4844 !
4845 ! Output arguments
4846 !
4847 !----------------------------------------------------------------------
4848 
4849  logical (SCRIP_logical), intent(out) ::
4850  & inpoly ! Is point in the polygon?
4851 
4852  logical (SCRIP_logical), intent(out) ::
4853  & lboundary ! Is point on the boundary of the polygon?
4854 
4855  integer (SCRIP_i4), intent(out) ::
4856  & edgeid ! if point is on boundary, which local
4857  ! edge is it on? (0 otherwise)
4858 
4859 !----------------------------------------------------------------------
4860 !
4861 ! Local variables
4862 !
4863 !----------------------------------------------------------------------
4864 
4865  integer (SCRIP_i4) :: n, next_n
4866 
4867  real (SCRIP_r8) :: x1, y1, x2, y2, vec1_x, vec1_y, vec2_x, vec2_y,
4868  & cross_product, minlon, maxlon, ptx_loc, pty_loc
4869 
4870  real (SCRIP_r8), dimension(ncorners) ::
4871  & cell_corner_lat_loc, cell_corner_lon_loc
4872 
4873 
4874  !***********************************************************
4875  !*** We should just remove the latlon argument since that is
4876  !*** the only coordinate system we are using it for
4877  !***********************************************************
4878 
4879 
4880  !***
4881  !*** here we take the cross product of the vector making
4882  !*** up each cell side with the vector formed by the vertex
4883  !*** and search point. if all the cross products are
4884  !*** positive, the point is contained in the cell.
4885  !***
4886 
4887  inpoly = .false.
4888  lboundary = .false.
4889  edgeid = 0
4890 
4891  if (.not. latlon) then
4892 
4893  do n = 1, ncorners
4894  next_n = mod(n,ncorners) + 1
4895 
4896  x1 = cell_corner_x(n)
4897  y1 = cell_corner_y(n)
4898  x2 = cell_corner_x(next_n)
4899  y2 = cell_corner_y(next_n)
4900 
4901  vec1_x = x2 - x1
4902  vec1_y = y2 - y1
4903  vec2_x = ptx - x1
4904  vec2_y = pty - y1
4905 
4906  cross_product = vec1_y*vec2_x - vec2_y*vec1_x
4907 
4908  !***
4909  !*** if the cross product for a side is zero, the point
4910  !*** lies exactly on the side or the side is degenerate
4911  !*** (zero length). if degenerate, set the cross
4912  !*** product to a positive number.
4913  !***
4914 
4915  if (abs(cross_product) < tiny) then
4916  if (vec1_x*vec1_x + vec1_y*vec1_y .le. tiny*tiny) then
4917  cross_product = one
4918  else
4919  lboundary = .true.
4920  edgeid = n
4921  endif
4922  else
4923 
4924  !***
4925  !*** if cross product is less than zero, this cell
4926  !*** doesn't work
4927  !***
4928  !*** Should we say "if (cp < zero .and. abs(cp) > tiny)" ?
4929 
4930  if (cross_product < zero) then
4931  inpoly = .false.
4932  lboundary = .false.
4933  return
4934  endif
4935  endif
4936 
4937  end do
4938 
4939  else
4940 
4941  !*** Checking in latlon space
4942  !*** If the grid cell coordinates spans more than pi radians
4943  !*** transform the coordinates so that they don't
4944 
4945  cell_corner_lat_loc = cell_corner_x
4946  cell_corner_lon_loc = cell_corner_y
4947 
4948  minlon = 9999.0
4949  maxlon = -9999.0
4950  do n = 1, ncorners
4951  if (cell_corner_lon_loc(n) < minlon) then
4952  minlon = cell_corner_lon_loc(n)
4953  endif
4954  if (cell_corner_lon_loc(n) > maxlon) then
4955  maxlon = cell_corner_lon_loc(n)
4956  endif
4957  enddo
4958 
4959  if (maxlon-minlon > pi) then
4960 
4961  do n = 1, ncorners
4962  if (cell_corner_lon_loc(n)-minlon > pi) then
4963  cell_corner_lon_loc(n) = cell_corner_lon_loc(n)-pi2
4964  endif
4965  enddo
4966 
4967  endif
4968 
4969  ptx_loc = ptx
4970  pty_loc = pty
4971  if (pty_loc - minlon > pi) then
4972  pty_loc = pty_loc - pi2
4973  else if (pty_loc - minlon < -pi) then
4974  pty_loc = pty_loc + pi2
4975  endif
4976 
4977 
4978  do n = 1, ncorners
4979  next_n = mod(n,ncorners) + 1
4980 
4981  x1 = cell_corner_lat_loc(n)
4982  y1 = cell_corner_lon_loc(n)
4983  x2 = cell_corner_lat_loc(next_n)
4984  y2 = cell_corner_lon_loc(next_n)
4985 
4986  vec1_x = x2 - x1
4987  vec1_y = y2 - y1
4988  vec2_x = ptx_loc - x1
4989  vec2_y = pty_loc - y1
4990 
4991  cross_product = vec1_y*vec2_x - vec2_y*vec1_x
4992 
4993  !***
4994  !*** if the cross product for a side is zero, the point
4995  !*** lies exactly on the side or the side is degenerate
4996  !*** (zero length). if degenerate, set the cross
4997  !*** product to a positive number.
4998  !***
4999 
5000  if (abs(cross_product) < tiny) then
5001  if (vec1_x*vec1_x + vec1_y*vec1_y .le. tiny*tiny) then
5002  cross_product = one
5003  else
5004  lboundary = .true.
5005  edgeid = n
5006  endif
5007  else
5008 
5009  !***
5010  !*** if cross product is less than zero, this cell
5011  !*** doesn't work
5012  !***
5013  !*** Should we say "if (cp < zero .and. abs(cp) > tiny)" ?
5014 
5015  if (cross_product < zero) then
5016  inpoly = .false.
5017  lboundary = .false.
5018  return
5019  endif
5020  endif
5021 
5022  end do
5023 
5024  endif
5025  !***
5026  !*** if cross products all positive, we found the location
5027  !***
5028 
5029  inpoly = .true.
5030  return
5031 
5032 !----------------------------------------------------------------------
5033 
5034  end subroutine ptinpoly
5035 
5036 !**********************************************************************
5037 
5038 
5039 
5040  subroutine ptinpolarpoly(ptx, pty, ncorners, cell_corner_x,
5041  & cell_corner_y, latlon, whichpole, inpoly, lboundary, edgeid)
5042 
5043 !----------------------------------------------------------------------
5044 !
5045 ! Check if point is in polygonal cell overlapping the pole
5046 ! Cannot check the containment as is in latlon space - We have
5047 ! to check by connecting each edge of the polygon to the pole
5048 ! and check containment in the resulting quadrilateral in latlon
5049 ! space
5050 ! The cell can be non-convex as long as the pole is 'visible' to
5051 ! all the edges of the polygon, i.e., we can connect the pole to
5052 ! each edge of the polygon and form a triangle with positive area
5053 !
5054 !----------------------------------------------------------------------
5055 
5056  implicit none
5057 
5058 !----------------------------------------------------------------------
5059 !
5060 ! Input arguments
5061 !
5062 !----------------------------------------------------------------------
5063 
5064  real (SCRIP_r8), intent(in) ::
5065  & ptx, pty ! Point to check
5066 
5067  integer (SCRIP_i4), intent(in) ::
5068  & ncorners ! Number of polygon corners
5069 
5070  real (SCRIP_r8), dimension(ncorners), intent(in) ::
5071  & cell_corner_x, ! Coordinates of cell corners
5072  & cell_corner_y ! Could be x-y or lat-lon or ...
5073 
5074  logical (SCRIP_logical), intent(in) ::
5075  & latlon ! Are coordinates in latlon space?
5076 
5077  integer (SCRIP_i4), intent(in) ::
5078  & whichpole ! South or North pole
5079 
5080 !----------------------------------------------------------------------
5081 !
5082 ! Output arguments
5083 !
5084 !----------------------------------------------------------------------
5085 
5086  logical (SCRIP_logical), intent(out) ::
5087  & inpoly ! Is point in the polygon?
5088 
5089  logical (SCRIP_logical), intent(out) ::
5090  & lboundary ! Is point on the boundary of the polygon?
5091 
5092  integer (SCRIP_i4), intent(out) ::
5093  & edgeid ! if point is on boundary, which local
5094  ! edge is it on? (0 otherwise)
5095 
5096 !----------------------------------------------------------------------
5097 !
5098 ! Local variables
5099 !
5100 !----------------------------------------------------------------------
5101 
5102  integer (SCRIP_i4) :: n, next_n, ledgeid
5103 
5104  real (SCRIP_r8), dimension(4) ::
5105  & pquad_corner_x, ! Coordinates of polar quad
5106  & pquad_corner_y
5107 
5108  real (SCRIP_r8) :: x1, y1, x2, y2, vec1_x, vec1_y, vec2_x, vec2_y,
5109  & cross_product, pole_lat
5110 
5111  pole_lat = whichpole*pih
5112 
5113  !***
5114  !*** This is a polygon that overlaps the pole
5115  !*** A normal point in polygon check could fail
5116  !*** So, with each edge of the polygon form a quadrilateral
5117  !*** in latlon space using the polar latitude and the longitude
5118  !*** values of the endpoints of the edge. Then check containment
5119  !*** of the point in this quadrilateral
5120  !***
5121 
5122  inpoly = .false.
5123  lboundary = .false.
5124 
5125  do n = 1, ncorners
5126  next_n = mod(n,ncorners) + 1
5127 
5128  pquad_corner_x(1) = cell_corner_x(n)
5129  pquad_corner_y(1) = cell_corner_y(n)
5130  pquad_corner_x(2) = cell_corner_x(next_n)
5131  pquad_corner_y(2) = cell_corner_y(next_n)
5132  pquad_corner_x(3) = pole_lat
5133  pquad_corner_y(3) = cell_corner_y(next_n)
5134  pquad_corner_x(4) = pole_lat
5135  pquad_corner_y(4) = cell_corner_y(n)
5136 
5137 
5138  call ptinpoly(ptx,pty,4,pquad_corner_x,pquad_corner_y,
5139  & latlon,inpoly,lboundary, ledgeid)
5140 
5141  if (inpoly) then
5142 
5143  if (lboundary) then
5144 
5145  !***
5146  !*** Check to see if the lboundary flag is being
5147  !*** triggered by the outer edge of the polygon or
5148  !*** by one of the artificial internal edges
5149  !***
5150 
5151  vec1_x = pquad_corner_x(2) - pquad_corner_x(1)
5152  vec1_y = pquad_corner_y(2) - pquad_corner_y(1)
5153  vec2_x = ptx - pquad_corner_x(1)
5154  vec2_y = pty - pquad_corner_y(1)
5155 
5156 
5157  if (latlon) then
5158 
5159  !***
5160  !*** check for 0,2pi crossings
5161  !***
5162 
5163  if (vec1_y > pi) vec1_y = vec1_y - pi2
5164  if (vec1_y < -pi) vec1_y = vec1_y + pi2
5165  if (vec2_y > pi) vec2_y = vec2_y - pi2
5166  if (vec2_y < -pi) vec2_y = vec2_y + pi2
5167 
5168  endif
5169 
5170  cross_product = vec1_y*vec2_x - vec2_y*vec1_x
5171 
5172  !***
5173  !*** if the cross product for a side is zero, the point
5174  !*** lies exactly on the side or the side is degenerate
5175  !*** (zero length). if degenerate, set the cross
5176  !*** product to a positive number.
5177  !***
5178 
5179  if (abs(cross_product) < tiny) then
5180  if (vec1_x .eq. zero .and. vec1_y .eq. zero) then
5181  cross_product = one
5182  lboundary = .false.
5183  else
5184  edgeid = n
5185  lboundary = .true.
5186  endif
5187  else
5188  lboundary = .false.
5189  endif
5190  endif ! if (lboundary)
5191 
5192  return ! pt in polygon
5193 
5194  endif ! if (inpoly)
5195 
5196  end do
5197 
5198  return ! pt outside polygon
5199 
5200 !----------------------------------------------------------------------
5201 
5202  end subroutine ptinpolarpoly
5203 
5204 !**********************************************************************
5205 
5206 
5207 
5208  subroutine ptinpolygen(ptx, pty, ncorners, cell_corner_x,
5209  & cell_corner_y, cell_center_x, cell_center_y,
5210  & latlon, inpoly, lboundary, edgeid)
5211 
5212 !----------------------------------------------------------------------
5213 !
5214 ! Check if point is in general (convex or mildly non-convex)
5215 ! polygonal cell by connecting each edge of the polygon to a
5216 ! a central point (average of vertices) and check containment in
5217 ! the resulting triangle
5218 !
5219 ! The cell can be non-convex as long as the 'center' is 'visible' to
5220 ! all the edges of the polygon, i.e., we can connect the 'center' to
5221 ! each edge of the polygon and form a triangle with positive area
5222 !
5223 !----------------------------------------------------------------------
5224 
5225  implicit none
5226 
5227 !----------------------------------------------------------------------
5228 !
5229 ! Input arguments
5230 !
5231 !----------------------------------------------------------------------
5232 
5233  real (SCRIP_r8), intent(in) ::
5234  & ptx, pty ! Point to check
5235 
5236  integer (SCRIP_i4), intent(in) ::
5237  & ncorners ! Number of polygon corners
5238 
5239  real (SCRIP_r8), dimension(ncorners), intent(in) ::
5240  & cell_corner_x, ! Coordinates of cell corners
5241  & cell_corner_y ! Could be x-y or lat-lon or ...
5242 
5243  real (SCRIP_r8), intent(in) ::
5244  & cell_center_x,
5245  & cell_center_y
5246 
5247  logical (SCRIP_logical), intent(in) ::
5248  & latlon ! Are coordinates in latlon space?
5249 
5250 !----------------------------------------------------------------------
5251 !
5252 ! Output arguments
5253 !
5254 !----------------------------------------------------------------------
5255 
5256  logical (SCRIP_logical), intent(out) ::
5257  & inpoly ! Is point in the polygon?
5258 
5259  logical (SCRIP_logical), intent(out) ::
5260  & lboundary ! Is point on the boundary of the polygon?
5261 
5262  integer (SCRIP_i4), intent(out) ::
5263  & edgeid ! if point is on boundary, which local
5264  ! edge is it on? (0 otherwise)
5265 
5266 !----------------------------------------------------------------------
5267 !
5268 ! Local variables
5269 !
5270 !----------------------------------------------------------------------
5271 
5272  integer (SCRIP_i4) :: n, next_n, ledgeid
5273 
5274  real (SCRIP_r8), dimension(3) ::
5275  & tri_corner_x, ! Coordinates of triangle
5276  & tri_corner_y
5277 
5278  real (SCRIP_r8) :: x1, y1, x2, y2, vec1_x, vec1_y, vec2_x, vec2_y,
5279  & cross_product
5280 
5281 
5282  !***
5283  !*** So, with each edge of the polygon form a triangle
5284  !*** by connecting a 'central' point to the endpoints of
5285  !*** the edge. Then check containment of the point in this tri
5286  !***
5287 
5288  inpoly = .false.
5289  lboundary = .false.
5290 
5291  do n = 1, ncorners
5292  next_n = mod(n,ncorners) + 1
5293 
5294  tri_corner_x(1) = cell_corner_x(n)
5295  tri_corner_y(1) = cell_corner_y(n)
5296  tri_corner_x(2) = cell_corner_x(next_n)
5297  tri_corner_y(2) = cell_corner_y(next_n)
5298  tri_corner_x(3) = cell_center_x
5299  tri_corner_y(3) = cell_center_y
5300 
5301  vec1_x = tri_corner_x(2) - tri_corner_x(1)
5302  vec1_y = tri_corner_y(2) - tri_corner_y(1)
5303 
5304  !*** Skip triangles arising from degenerate edges
5305 
5306  if (vec1_x*vec1_x+vec1_y*vec1_y .le. tiny*tiny) cycle
5307 
5308  call ptinpoly(ptx,pty,3,tri_corner_x,tri_corner_y,
5309  & latlon,inpoly,lboundary, ledgeid)
5310 
5311  if (inpoly) then
5312 
5313  if (lboundary) then
5314 
5315  !***
5316  !*** Check to see if the lboundary flag is being
5317  !*** triggered by the outer edge of the polygon or
5318  !*** by one of the artificial internal edges
5319  !***
5320 
5321  vec2_x = ptx - tri_corner_x(1)
5322  vec2_y = pty - tri_corner_y(1)
5323 
5324 
5325  if (latlon) then
5326 
5327  !***
5328  !*** check for 0,2pi crossings
5329  !***
5330 
5331  if (vec1_y > pi) vec1_y = vec1_y - pi2
5332  if (vec1_y < -pi) vec1_y = vec1_y + pi2
5333  if (vec2_y > pi) vec2_y = vec2_y - pi2
5334  if (vec2_y < -pi) vec2_y = vec2_y + pi2
5335 
5336  endif
5337 
5338  cross_product = vec1_y*vec2_x - vec2_y*vec1_x
5339 
5340  !***
5341  !*** if the cross product for a side is zero, the point
5342  !*** lies exactly on the side or the side is degenerate
5343  !*** (zero length). if degenerate, set the cross
5344  !*** product to a positive number.
5345  !***
5346 
5347  if (abs(cross_product) < tiny) then
5348  if (vec1_x*vec1_x+vec1_y*vec1_y .le. tiny*tiny) then
5349  cross_product = one
5350  lboundary = .false.
5351  else
5352  edgeid = n
5353  lboundary = .true.
5354  endif
5355  else
5356  lboundary = .false.
5357  endif
5358  endif ! if (lboundary)
5359 
5360  return ! pt in polygon
5361 
5362  endif ! if (inpoly)
5363 
5364  end do
5365 
5366  return ! pt outside polygon
5367 
5368 !----------------------------------------------------------------------
5369 
5370  end subroutine ptinpolygen
5371 
5372 !**********************************************************************
5373 
5374 
5375 
5376  subroutine ptinpolygen2(ptx, pty, ncorners, cell_corner_x,
5377  & cell_corner_y, latlon, inpoly, lboundary, edgeid)
5378 
5379 !----------------------------------------------------------------------
5380 !
5381 ! Check if point is in general (convex or mildly non-convex)
5382 ! polygonal cell by connecting each edge of the polygon to a
5383 ! a central point (average of vertices) and check containment in
5384 ! the resulting triangle
5385 !
5386 ! The cell can be non-convex as long as the 'center' is 'visible' to
5387 ! all the edges of the polygon, i.e., we can connect the 'center' to
5388 ! each edge of the polygon and form a triangle with positive area
5389 !
5390 !----------------------------------------------------------------------
5391 
5392  implicit none
5393 
5394 !----------------------------------------------------------------------
5395 !
5396 ! Input arguments
5397 !
5398 !----------------------------------------------------------------------
5399 
5400  real (SCRIP_r8), intent(in) ::
5401  & ptx, pty ! Point to check
5402 
5403  integer (SCRIP_i4), intent(in) ::
5404  & ncorners ! Number of polygon corners
5405 
5406  real (SCRIP_r8), dimension(ncorners), intent(in) ::
5407  & cell_corner_x, ! Coordinates of cell corners
5408  & cell_corner_y ! Could be x-y or lat-lon or ...
5409 
5410  logical (SCRIP_logical), intent(in) ::
5411  & latlon ! Are coordinates in latlon space?
5412 
5413 !----------------------------------------------------------------------
5414 !
5415 ! Output arguments
5416 !
5417 !----------------------------------------------------------------------
5418 
5419  logical (SCRIP_logical), intent(out) ::
5420  & inpoly ! Is point in the polygon?
5421 
5422  logical (SCRIP_logical), intent(out) ::
5423  & lboundary ! Is point on the boundary of the polygon?
5424 
5425  integer (SCRIP_i4), intent(out) ::
5426  & edgeid ! if point is on boundary, which local
5427  ! edge is it on? (0 otherwise)
5428 
5429 !----------------------------------------------------------------------
5430 !
5431 ! Local variables
5432 !
5433 !----------------------------------------------------------------------
5434 
5435  integer (SCRIP_i4) :: c, n, next_n
5436 
5437  real (SCRIP_r8) :: x1, y1, x2, y2, vec1_x, vec1_y, vec2_x, vec2_y,
5438  & vec3_x, vec3_y, vec1_len, vec2_len, vec3_len,
5439  & cross_product, dot_product
5440 
5441 
5442  !***
5443  !*** So, with each edge of the polygon form a triangle
5444  !*** by connecting a 'central' point to the endpoints of
5445  !*** the edge. Then check containment of the point in this tri
5446  !***
5447 
5448  inpoly = .false.
5449  lboundary = .false.
5450 
5451  c = 0
5452  do n = 1, ncorners
5453  next_n = mod(n,ncorners) + 1
5454 
5455  x1 = cell_corner_x(n)
5456  y1 = cell_corner_y(n)
5457  x2 = cell_corner_x(next_n)
5458  y2 = cell_corner_y(next_n)
5459 
5460  if (((y1 > pty .and. y2 <= pty) .or.
5461  & (y2 > pty .and. y1 <= pty)) .and.
5462  & (ptx <= (x1 + (pty-y1)*(x2-x1)/(y2-y1)))) then
5463 
5464  c = 1 - c
5465 
5466  endif
5467  enddo
5468 
5469  if (c .eq. 1) inpoly = .true.
5470 
5471 
5472  !*** Check if the point is on the boundary of the polygon
5473 
5474  do n = 1, ncorners
5475 
5476  next_n = mod(n,ncorners) + 1
5477 
5478  x1 = cell_corner_x(n)
5479  y1 = cell_corner_y(n)
5480  x2 = cell_corner_x(next_n)
5481  y2 = cell_corner_y(next_n)
5482 
5483  vec1_x = x2 - x1
5484  vec1_y = y2 - y1
5485  vec1_len = sqrt(vec1_x*vec1_x + vec1_y*vec1_y)
5486  vec1_x = vec1_x/vec1_len
5487  vec1_y = vec1_y/vec1_len
5488 
5489  vec2_x = ptx - x1
5490  vec2_y = pty - y1
5491  vec2_len = sqrt(vec2_x*vec2_x + vec2_y*vec2_y)
5492 
5493  cross_product = vec1_x*vec2_y - vec2_x*vec1_y
5494  if (abs(cross_product) > tiny .and. vec2_len > tiny) then
5495  cross_product = cross_product/vec2_len
5496  endif
5497 
5498  if (abs(cross_product) < 1e5*tiny .and.
5499  & abs(cross_product) > 10*tiny) then
5500 
5501  !*** Sometimes when the point is too close to a vertex
5502  !*** then the cross product computation has errors due
5503  !*** to subtraction of two small numbers - So check w.r.t.
5504  !*** other vertex of the segment as well
5505 
5506  vec3_x = ptx - x2
5507  vec3_y = pty - y2
5508  vec3_len = sqrt(vec3_x*vec3_x + vec3_y*vec3_y)
5509 
5510  cross_product = -vec1_x*vec3_y + vec1_y*vec3_x
5511  if (abs(cross_product) > tiny .and. vec3_len > tiny) then
5512  !***
5513  !*** Normalize only if we won't be dividing two small
5514  !*** numbers
5515  cross_product = cross_product/vec3_len
5516  endif
5517  endif
5518 
5519  if (abs(cross_product) < 10*tiny) then
5520 
5521  if (vec1_x*vec1_x+vec1_y*vec1_y .le. tiny*tiny) then
5522  cross_product = one
5523  else
5524  dot_product = vec1_x*vec2_x + vec1_y*vec2_y
5525 
5526  if (dot_product >= 0 .and. dot_product <= vec1_len) then
5527  inpoly = .true.
5528  lboundary = .true.
5529  edgeid = n
5530  exit
5531  endif
5532  endif
5533  endif
5534 
5535  enddo
5536 
5537  return
5538 
5539 !----------------------------------------------------------------------
5540 
5541  end subroutine ptinpolygen2
5542 
5543 !**********************************************************************
5544 
5545 
5546 
5547 
5548 
5549 !----------------------------------------------------------------------
5550 !
5551 !----------------------------------------------------------------------
5552 
5553  subroutine get_srch_cells(cell_add, cell_grid_num, srch_grid_num,
5554  & num_srch_cells, srch_add, srch_corners,
5555  & srch_corner_lat, srch_corner_lon,
5556  & srch_center_lat, srch_center_lon)
5557 
5558 !----------------------------------------------------------------------
5559 !
5560 ! Input arguments
5561 !
5562 !----------------------------------------------------------------------
5563 
5564  integer (SCRIP_i4), intent(in) ::
5565  & cell_add, ! cell in whose nbrhood we must find other
5566  & cell_grid_num, ! cells grid number from which 'cell_add'
5567  & srch_grid_num ! is grid number in which we must find
5568  ! search cells
5569 
5570 !----------------------------------------------------------------------
5571 !
5572 ! Output arguments
5573 !
5574 !----------------------------------------------------------------------
5575 
5576  integer (SCRIP_i4), intent(out) ::
5577  & num_srch_cells,
5578  & srch_corners ! Number of corners for search cells
5579 
5580  integer (SCRIP_i4), dimension(:), allocatable, intent(out) ::
5581  & srch_add ! Global addresses of search cells
5582 
5583  real (SCRIP_r8), dimension(:,:), allocatable, intent(out) ::
5584  & srch_corner_lat, srch_corner_lon
5585 
5586  real (SCRIP_r8), dimension(:), allocatable, intent(out) ::
5587  & srch_center_lat, srch_center_lon
5588 
5589 
5590 !-----------------------------------------------------------------------
5591 !
5592 ! Local arguments
5593 !
5594 !-----------------------------------------------------------------------
5595 
5596  logical (SCRIP_logical), dimension(:), allocatable ::
5597  & srch_mask
5598 
5599  integer (SCRIP_i4) :: grid1_add, grid2_add, max_add, min_add,
5600  & n
5601 
5602 !-----------------------------------------------------------------------
5603 
5604  num_srch_cells = 0
5605 
5606  !***
5607  !*** restrict searches first using search bins
5608  !***
5609 
5610  if (last_cell_add_get_srch_cells /= cell_add .or.
5611  & last_cell_grid_num_get_srch_cells /= cell_grid_num .or.
5612  & last_srch_grid_num_get_srch_cells /= srch_grid_num) then
5613 
5614  if (first_call_get_srch_cells) then
5615  first_call_get_srch_cells = .false.
5616  num_srch_cells_loc_get_srch_cells = 0
5617  srch_corners_loc_get_srch_cells = 0
5618  last_cell_add_get_srch_cells = 0
5619  last_cell_grid_num_get_srch_cells = 0
5620  last_srch_grid_num_get_srch_cells = 0
5621  else
5622  if (num_srch_cells_loc_get_srch_cells .gt. 0) then
5623  deallocate(srch_add_loc_get_srch_cells,
5624  & srch_corner_lat_loc_get_srch_cells,
5625  & srch_corner_lon_loc_get_srch_cells,
5626  & srch_center_lat_loc_get_srch_cells,
5627  & srch_center_lon_loc_get_srch_cells)
5628  endif
5629 
5630  endif
5631 
5632 
5633  last_cell_add_get_srch_cells = cell_add
5634  last_cell_grid_num_get_srch_cells = cell_grid_num
5635  last_srch_grid_num_get_srch_cells = srch_grid_num
5636 
5637 
5638  if (cell_grid_num == 1) then
5639 
5640  if (srch_grid_num == 1) then
5641 
5642  !*** Grid 1 neighbors of grid 1 cell
5643 
5644  allocate(srch_mask(grid1_size))
5645 
5646  min_add = grid1_size
5647  max_add = 1
5648  do n=1,num_srch_bins
5649  if (cell_add >= bin_addr1(1,n) .and.
5650  & cell_add <= bin_addr1(2,n)) then
5651  min_add = min(min_add, bin_addr1(1,n))
5652  max_add = max(max_add, bin_addr1(2,n))
5653  endif
5654  end do
5655 
5656  !***
5657  !*** further restrict searches using bounding boxes
5658  !***
5659 
5660  num_srch_cells_loc_get_srch_cells = 0
5661  do grid1_add = min_add,max_add
5662  srch_mask(grid1_add) =
5663  & (grid1_bound_box(1,grid1_add) <=
5664  & grid1_bound_box(2,cell_add)) .and.
5665  & (grid1_bound_box(2,grid1_add) >=
5666  & grid1_bound_box(1,cell_add)) .and.
5667  & (grid1_bound_box(3,grid1_add) <=
5668  & grid1_bound_box(4,cell_add)) .and.
5669  & (grid1_bound_box(4,grid1_add) >=
5670  & grid1_bound_box(3,cell_add))
5671 
5672  if (srch_mask(grid1_add))
5673  & num_srch_cells_loc_get_srch_cells =
5674  & num_srch_cells_loc_get_srch_cells+1
5675  end do
5676 
5677  if (num_srch_cells_loc_get_srch_cells /= 0) then
5678 
5679  !***
5680  !*** create search arrays
5681  !***
5682 
5683  allocate(srch_add_loc_get_srch_cells
5684  & (num_srch_cells_loc_get_srch_cells),
5685  & srch_corner_lat_loc_get_srch_cells
5686  & (grid1_corners,num_srch_cells_loc_get_srch_cells),
5687  & srch_corner_lon_loc_get_srch_cells
5688  & (grid1_corners,num_srch_cells_loc_get_srch_cells),
5689  & srch_center_lat_loc_get_srch_cells
5690  & (num_srch_cells_loc_get_srch_cells),
5691  & srch_center_lon_loc_get_srch_cells
5692  & (num_srch_cells_loc_get_srch_cells))
5693 
5694  n = 0
5695  do grid1_add = min_add,max_add
5696  if (srch_mask(grid1_add)) then
5697  n = n+1
5698  srch_add_loc_get_srch_cells(n) = grid1_add
5699  srch_corner_lat_loc_get_srch_cells(:,n) =
5700  & grid1_corner_lat(:,grid1_add)
5701  srch_corner_lon_loc_get_srch_cells(:,n) =
5702  & grid1_corner_lon(:,grid1_add)
5703  srch_center_lat_loc_get_srch_cells(n) =
5704  & grid1_center_lat(grid1_add)
5705  srch_center_lon_loc_get_srch_cells(n) =
5706  & grid1_center_lon(grid1_add)
5707  endif
5708  end do
5709 
5710  srch_corners_loc_get_srch_cells = grid1_corners
5711  endif
5712 
5713  deallocate(srch_mask)
5714 
5715  else
5716 
5717  !*** Grid 2 neighbors of grid 1 cell
5718 
5719  allocate(srch_mask(grid2_size))
5720 
5721  min_add = grid2_size
5722  max_add = 1
5723  do n=1,num_srch_bins
5724  if (cell_add >= bin_addr1(1,n) .and.
5725  & cell_add <= bin_addr1(2,n)) then
5726  min_add = min(min_add, bin_addr2(1,n))
5727  max_add = max(max_add, bin_addr2(2,n))
5728  endif
5729  end do
5730 
5731  !***
5732  !*** further restrict searches using bounding boxes
5733  !***
5734 
5735  num_srch_cells_loc_get_srch_cells = 0
5736  do grid2_add = min_add,max_add
5737  srch_mask(grid2_add) =
5738  & (grid2_bound_box(1,grid2_add) <=
5739  & grid1_bound_box(2,cell_add)) .and.
5740  & (grid2_bound_box(2,grid2_add) >=
5741  & grid1_bound_box(1,cell_add)) .and.
5742  & (grid2_bound_box(3,grid2_add) <=
5743  & grid1_bound_box(4,cell_add)) .and.
5744  & (grid2_bound_box(4,grid2_add) >=
5745  & grid1_bound_box(3,cell_add))
5746 
5747 
5748 
5749  if (srch_mask(grid2_add))
5750  & num_srch_cells_loc_get_srch_cells =
5751  & num_srch_cells_loc_get_srch_cells+1
5752  end do
5753 
5754 
5755  if (num_srch_cells_loc_get_srch_cells /= 0) then
5756 
5757  !***
5758  !*** create search arrays
5759  !***
5760 
5761  allocate(srch_add_loc_get_srch_cells(
5762  & num_srch_cells_loc_get_srch_cells),
5763  & srch_corner_lat_loc_get_srch_cells(
5764  & grid2_corners,num_srch_cells_loc_get_srch_cells),
5765  & srch_corner_lon_loc_get_srch_cells(
5766  & grid2_corners,num_srch_cells_loc_get_srch_cells),
5767  & srch_center_lat_loc_get_srch_cells(
5768  & num_srch_cells_loc_get_srch_cells),
5769  & srch_center_lon_loc_get_srch_cells(
5770  & num_srch_cells_loc_get_srch_cells))
5771 
5772  n = 0
5773  do grid2_add = min_add,max_add
5774  if (srch_mask(grid2_add)) then
5775  n = n+1
5776  srch_add_loc_get_srch_cells(n) = grid2_add
5777  srch_corner_lat_loc_get_srch_cells(:,n) =
5778  & grid2_corner_lat(:,grid2_add)
5779  srch_corner_lon_loc_get_srch_cells(:,n) =
5780  & grid2_corner_lon(:,grid2_add)
5781  srch_center_lat_loc_get_srch_cells(n) =
5782  & grid2_center_lat(grid2_add)
5783  srch_center_lon_loc_get_srch_cells(n) =
5784  & grid2_center_lon(grid2_add)
5785  endif
5786  end do
5787 
5788  srch_corners_loc_get_srch_cells = grid2_corners
5789  endif
5790 
5791  deallocate(srch_mask)
5792  endif
5793 
5794  else
5795 
5796  if (srch_grid_num == 1) then
5797 
5798  !*** Grid 1 neighbors of grid 2 cell
5799 
5800  allocate(srch_mask(grid1_size))
5801 
5802  min_add = grid1_size
5803  max_add = 1
5804  do n=1,num_srch_bins
5805  if (cell_add >= bin_addr2(1,n) .and.
5806  & cell_add <= bin_addr2(2,n)) then
5807  min_add = min(min_add, bin_addr1(1,n))
5808  max_add = max(max_add, bin_addr1(2,n))
5809  endif
5810  end do
5811 
5812  !***
5813  !*** further restrict searches using bounding boxes
5814  !***
5815 
5816  num_srch_cells_loc_get_srch_cells = 0
5817  do grid1_add = min_add,max_add
5818  srch_mask(grid1_add) =
5819  & (grid1_bound_box(1,grid1_add) <=
5820  & grid2_bound_box(2,cell_add)) .and.
5821  & (grid1_bound_box(2,grid1_add) >=
5822  & grid2_bound_box(1,cell_add)) .and.
5823  & (grid1_bound_box(3,grid1_add) <=
5824  & grid2_bound_box(4,cell_add)) .and.
5825  & (grid1_bound_box(4,grid1_add) >=
5826  & grid2_bound_box(3,cell_add))
5827 
5828  if (srch_mask(grid1_add))
5829  & num_srch_cells_loc_get_srch_cells =
5830  & num_srch_cells_loc_get_srch_cells+1
5831  end do
5832 
5833 
5834  if (num_srch_cells_loc_get_srch_cells /= 0) then
5835 
5836  !***
5837  !*** create search arrays
5838  !***
5839 
5840  allocate(srch_add_loc_get_srch_cells(
5841  & num_srch_cells_loc_get_srch_cells),
5842  & srch_corner_lat_loc_get_srch_cells(
5843  & grid1_corners,num_srch_cells_loc_get_srch_cells),
5844  & srch_corner_lon_loc_get_srch_cells(
5845  & grid1_corners,num_srch_cells_loc_get_srch_cells),
5846  & srch_center_lat_loc_get_srch_cells(
5847  & num_srch_cells_loc_get_srch_cells),
5848  & srch_center_lon_loc_get_srch_cells(
5849  & num_srch_cells_loc_get_srch_cells))
5850 
5851  n = 0
5852  do grid1_add = min_add,max_add
5853  if (srch_mask(grid1_add)) then
5854  n = n+1
5855  srch_add_loc_get_srch_cells(n) = grid1_add
5856  srch_corner_lat_loc_get_srch_cells(:,n) =
5857  & grid1_corner_lat(:,grid1_add)
5858  srch_corner_lon_loc_get_srch_cells(:,n) =
5859  & grid1_corner_lon(:,grid1_add)
5860  srch_center_lat_loc_get_srch_cells(n) =
5861  & grid1_center_lat(grid1_add)
5862  srch_center_lon_loc_get_srch_cells(n) =
5863  & grid1_center_lon(grid1_add)
5864  endif
5865  end do
5866 
5867  srch_corners_loc_get_srch_cells = grid1_corners
5868  endif
5869 
5870  deallocate(srch_mask)
5871 
5872  else
5873 
5874  !*** Grid 2 neighbors of grid 2 cell
5875 
5876  allocate(srch_mask(grid2_size))
5877 
5878  min_add = grid2_size
5879  max_add = 1
5880  do n=1,num_srch_bins
5881  if (cell_add >= bin_addr2(1,n) .and.
5882  & cell_add <= bin_addr2(2,n)) then
5883  min_add = min(min_add, bin_addr2(1,n))
5884  max_add = max(max_add, bin_addr2(2,n))
5885  endif
5886  end do
5887 
5888  !***
5889  !*** further restrict searches using bounding boxes
5890  !***
5891 
5892  num_srch_cells_loc_get_srch_cells = 0
5893  do grid2_add = min_add,max_add
5894  srch_mask(grid2_add) =
5895  & (grid2_bound_box(1,grid2_add) <=
5896  & grid2_bound_box(2,cell_add)) .and.
5897  & (grid2_bound_box(2,grid2_add) >=
5898  & grid2_bound_box(1,cell_add)) .and.
5899  & (grid2_bound_box(3,grid2_add) <=
5900  & grid2_bound_box(4,cell_add)) .and.
5901  & (grid2_bound_box(4,grid2_add) >=
5902  & grid2_bound_box(3,cell_add))
5903 
5904  if (srch_mask(grid2_add))
5905  & num_srch_cells_loc_get_srch_cells =
5906  & num_srch_cells_loc_get_srch_cells+1
5907  end do
5908 
5909 
5910  if (num_srch_cells_loc_get_srch_cells /= 0) then
5911 
5912  !***
5913  !*** create search arrays
5914  !***
5915 
5916  allocate(srch_add_loc_get_srch_cells(
5917  & num_srch_cells_loc_get_srch_cells),
5918  & srch_corner_lat_loc_get_srch_cells(
5919  & grid2_corners,num_srch_cells_loc_get_srch_cells),
5920  & srch_corner_lon_loc_get_srch_cells(
5921  & grid2_corners,num_srch_cells_loc_get_srch_cells),
5922  & srch_center_lat_loc_get_srch_cells(
5923  & num_srch_cells_loc_get_srch_cells),
5924  & srch_center_lon_loc_get_srch_cells(
5925  & num_srch_cells_loc_get_srch_cells))
5926 
5927  n = 0
5928  do grid2_add = min_add,max_add
5929  if (srch_mask(grid2_add)) then
5930  n = n+1
5931  srch_add_loc_get_srch_cells(n) = grid2_add
5932  srch_corner_lat_loc_get_srch_cells(:,n) =
5933  & grid2_corner_lat(:,grid2_add)
5934  srch_corner_lon_loc_get_srch_cells(:,n) =
5935  & grid2_corner_lon(:,grid2_add)
5936  srch_center_lat_loc_get_srch_cells(n) =
5937  & grid2_center_lat(grid2_add)
5938  srch_center_lon_loc_get_srch_cells(n) =
5939  & grid2_center_lon(grid2_add)
5940  endif
5941  end do
5942 
5943  srch_corners_loc_get_srch_cells = grid2_corners
5944  endif
5945 
5946  deallocate(srch_mask)
5947 
5948  endif
5949 
5950  endif
5951 
5952  endif
5953 
5954 
5955  num_srch_cells = num_srch_cells_loc_get_srch_cells
5956 
5957  if (num_srch_cells .eq. 0) then
5958  return
5959  endif
5960 
5961  srch_corners = srch_corners_loc_get_srch_cells
5962  allocate(srch_add(num_srch_cells),
5963  & srch_corner_lat(srch_corners,num_srch_cells),
5964  & srch_corner_lon(srch_corners,num_srch_cells),
5965  & srch_center_lat(num_srch_cells),
5966  & srch_center_lon(num_srch_cells))
5967  srch_add = srch_add_loc_get_srch_cells
5968  srch_corner_lat = srch_corner_lat_loc_get_srch_cells
5969  srch_corner_lon = srch_corner_lon_loc_get_srch_cells
5970  srch_center_lat = srch_center_lat_loc_get_srch_cells
5971  srch_center_lon = srch_center_lon_loc_get_srch_cells
5972 
5973  end subroutine get_srch_cells
5974 
5975 
5976 !**********************************************************************
5977 
5978 
5979 !----------------------------------------------------------------------
5980 !
5981 ! Find cell adjacent to edge (edge_id) of given cell (cell_add)
5982 !
5983 !----------------------------------------------------------------------
5984 
5985  subroutine find_adj_cell(cell_add, edge_id, cell_grid_num,
5986  & adj_add)
5987 
5988 !----------------------------------------------------------------------
5989 !
5990 ! Input variables
5991 !
5992 !----------------------------------------------------------------------
5993 
5994  integer (SCRIP_i4), intent(in) ::
5995  & cell_add, ! cell whose edge we are checking
5996  & edge_id, ! index of edge that we are check
5997  & cell_grid_num ! grid to which cell belongs
5998 
5999 !----------------------------------------------------------------------
6000 !
6001 ! Output variables
6002 !
6003 !----------------------------------------------------------------------
6004 
6005  integer (SCRIP_i4), intent(out) :: adj_add
6006 
6007 !----------------------------------------------------------------------
6008 !
6009 ! Local variables
6010 !
6011 !----------------------------------------------------------------------
6012 
6013  integer (SCRIP_i4) :: i, inx, n, global_add
6014  logical (SCRIP_logical) :: found
6015  real (SCRIP_r8) :: lat1, lon1, lat2, lon2
6016 
6017  adj_add = 0
6018 
6019  if (cell_grid_num .eq. 1) then
6020 
6021  i = edge_id
6022  inx = 1 + mod(edge_id,grid1_corners)
6023 
6024  lat1 = grid1_corner_lat(i,cell_add)
6025  lon1 = grid1_corner_lon(i,cell_add)
6026  lat2 = grid1_corner_lat(inx,cell_add)
6027  lon2 = grid1_corner_lon(inx,cell_add)
6028 
6029  !***
6030  !*** Often the cell with the next or previous index is
6031  !*** the adjacent cell. Check that first
6032  !***
6033 
6034  if (cell_add .lt. grid1_size) then
6035 
6036  global_add = cell_add + 1
6037 
6038  do i = 1, grid1_corners
6039  inx = mod(i,grid1_corners)+1
6040  if (abs(grid1_corner_lat(inx,global_add)-lat1) .le. tiny
6041  & .and.
6042  & abs(grid1_corner_lat(i,global_add)-lat2) .le. tiny
6043  & .and.
6044  & abs(grid1_corner_lon(inx,global_add)-lon1) .le. tiny
6045  & .and.
6046  & abs(grid1_corner_lon(i,global_add)-lon2) .le. tiny)
6047  & then
6048 
6049  adj_add = global_add
6050  return
6051  endif
6052  enddo
6053 
6054  endif
6055 
6056  if (cell_add .gt. 1) then
6057 
6058  global_add = cell_add - 1
6059 
6060  do i = 1, grid1_corners
6061  inx = mod(i,grid1_corners)+1
6062  if (abs(grid1_corner_lat(inx,global_add)-lat1) .le. tiny
6063  & .and.
6064  & abs(grid1_corner_lat(i,global_add)-lat2) .le. tiny
6065  & .and.
6066  & abs(grid1_corner_lon(inx,global_add)-lon1) .le. tiny
6067  & .and.
6068  & abs(grid1_corner_lon(i,global_add)-lon2) .le. tiny)
6069  & then
6070 
6071  adj_add = global_add
6072  return
6073  endif
6074  enddo
6075 
6076  endif
6077 
6078 
6079 
6080  else
6081 
6082  i = edge_id
6083  inx = 1 + mod(edge_id,grid2_corners)
6084 
6085  lat1 = grid2_corner_lat(i,cell_add)
6086  lon1 = grid2_corner_lon(i,cell_add)
6087  lat2 = grid2_corner_lat(inx,cell_add)
6088  lon2 = grid2_corner_lon(inx,cell_add)
6089 
6090 
6091  !***
6092  !*** Often the cell with the next or previous index is
6093  !*** the adjacent cell. Check that first
6094  !***
6095 
6096  if (cell_add .lt. grid2_size) then
6097 
6098  global_add = cell_add + 1
6099 
6100  do i = 1, grid2_corners
6101  inx = mod(i,grid2_corners)+1
6102  if (abs(grid2_corner_lat(inx,global_add)-lat1) .le. tiny
6103  & .and.
6104  & abs(grid2_corner_lat(i,global_add)-lat2) .le. tiny
6105  & .and.
6106  & abs(grid2_corner_lon(inx,global_add)-lon1) .le. tiny
6107  & .and.
6108  & abs(grid2_corner_lon(i,global_add)-lon2) .le. tiny)
6109  & then
6110 
6111  adj_add = global_add
6112  return
6113  endif
6114  enddo
6115 
6116  endif
6117 
6118  if (cell_add .gt. 1) then
6119 
6120  global_add = cell_add - 1
6121 
6122  do i = 1, grid2_corners
6123  inx = mod(i,grid2_corners)+1
6124  if (abs(grid2_corner_lat(inx,global_add)-lat1) .le. tiny
6125  & .and.
6126  & abs(grid2_corner_lat(i,global_add)-lat2) .le. tiny
6127  & .and.
6128  & abs(grid2_corner_lon(inx,global_add)-lon1) .le. tiny
6129  & .and.
6130  & abs(grid2_corner_lon(i,global_add)-lon2) .le. tiny)
6131  & then
6132 
6133  adj_add = global_add
6134  return
6135  endif
6136  enddo
6137 
6138  endif
6139 
6140 
6141  endif
6142 
6143 
6144  if (cell_add /= last_cell_find_adj_cell .or.
6145  & cell_grid_num /= last_cell_grid_num_find_adj_cell) then
6146 
6147  last_cell_find_adj_cell = cell_add
6148  last_cell_grid_num_find_adj_cell = cell_grid_num
6149 
6150  if (first_call_find_adj_cell) then
6151  first_call_find_adj_cell = .false.
6152  last_cell_find_adj_cell = 0
6153  last_cell_grid_num_find_adj_cell = 0
6154  else
6155  if (num_srch_cells_find_adj_cell .gt. 0) then
6156  deallocate(srch_add_find_adj_cell,
6157  & srch_corner_lat_find_adj_cell,
6158  & srch_corner_lon_find_adj_cell,
6159  & srch_center_lat_find_adj_cell,
6160  & srch_center_lon_find_adj_cell)
6161  endif
6162  endif
6163 
6164  call get_srch_cells(cell_add, cell_grid_num, cell_grid_num,
6165  & num_srch_cells_find_adj_cell, srch_add_find_adj_cell,
6166  & srch_corners_find_adj_cell, srch_corner_lat_find_adj_cell,
6167  & srch_corner_lon_find_adj_cell,
6168  & srch_center_lat_find_adj_cell,
6169  & srch_center_lon_find_adj_cell)
6170 
6171  endif
6172 
6173 
6174  found = .false.
6175  do n = 1, num_srch_cells_find_adj_cell
6176 
6177  global_add = srch_add_find_adj_cell(n)
6178 
6179  do i = 1, srch_corners_find_adj_cell
6180  inx = mod(i,srch_corners_find_adj_cell)+1
6181  if (abs(srch_corner_lat_find_adj_cell(inx,n)-lat1) .le. tiny
6182  & .and.
6183  & abs(srch_corner_lat_find_adj_cell(i,n)-lat2) .le. tiny
6184  & .and.
6185  & abs(srch_corner_lon_find_adj_cell(inx,n)-lon1) .le.tiny
6186  & .and.
6187  & abs(srch_corner_lon_find_adj_cell(i,n)-lon2) .le. tiny)
6188  & then
6189 
6190  adj_add = global_add
6191  found = .true.
6192 
6193  exit
6194  endif
6195  enddo
6196 
6197  if (found) exit
6198 
6199  enddo
6200 
6201  return
6202  end subroutine find_adj_cell
6203 
6204 
6205 !----------------------------------------------------------------------
6206 !
6207 ! Given points inside and outside a cell, converge to the boundary
6208 !
6209 !----------------------------------------------------------------------
6210 
6211 
6212  subroutine converge_to_bdry(cell_add, cell_grid_num,
6213  & ncorners, cell_corner_lat,
6214  & cell_corner_lon, cell_center_lat, cell_center_lon,
6215  & inpt_lat, inpt_lon, outpt_lat, outpt_lon,
6216  & bpt_lat, bpt_lon, bedgeid)
6217 
6218 !----------------------------------------------------------------------
6219 !
6220 ! Input variables
6221 !
6222 !----------------------------------------------------------------------
6223 
6224  integer (SCRIP_i4), intent(in) ::
6225  & cell_add, ! Cell in which we are operating
6226  & cell_grid_num, ! Grid to which cell belongs
6227  & ncorners ! Number of corners in cell
6228 
6229  real (SCRIP_r8), dimension(ncorners), intent(in) ::
6230  & cell_corner_lat, ! Latitude values of cell corners
6231  & cell_corner_lon ! Longitude values of cell corners
6232 
6233  real (SCRIP_r8), intent(in) ::
6234  & cell_center_lat, ! Latitude of cell center
6235  & cell_center_lon, ! Longitude of cell center,
6236  & inpt_lat, ! Latitude of inside point
6237  & inpt_lon, ! Longitude of inside point
6238  & outpt_lat, ! Latitude of outside point
6239  & outpt_lon ! Longitude of outside point
6240 
6241 
6242 !----------------------------------------------------------------------
6243 !
6244 ! Output variables
6245 !
6246 !----------------------------------------------------------------------
6247 
6248  real (SCRIP_r8), intent(out) ::
6249  & bpt_lat, ! Latitude of boundary point
6250  & bpt_lon ! Longitude of boundary point
6251 
6252  integer (SCRIP_i4), intent(out) ::
6253  & bedgeid ! ID of edge that point converged to
6254 
6255 !----------------------------------------------------------------------
6256 !
6257 ! Local variables
6258 !
6259 !----------------------------------------------------------------------
6260 
6261  logical (SCRIP_logical) ::
6262  & converged,
6263  & lboundary,
6264  & inpoly
6265 
6266  integer (SCRIP_i4) ::
6267  & it
6268 
6269  real (SCRIP_r8) ::
6270  & lat1, lon1,
6271  & lat2, lon2,
6272  & midlat, midlon
6273 
6274  bedgeid = 0
6275 
6276  lat1 = inpt_lat
6277  lon1 = inpt_lon
6278  lat2 = outpt_lat
6279  lon2 = outpt_lon
6280 
6281 
6282  converged = .false.
6283  it = 0
6284  do while (.not. converged)
6285 
6286  midlat = (lat1+lat2)/2.0
6287  if (abs(lon1-lon2) < pi) then
6288  midlon = (lon1+lon2)/2.0
6289  else
6290  midlon = (lon1+lon2)/2.0 - pi2
6291  endif
6292 
6293 
6294  call ptincell(midlat, midlon,
6295  & cell_add, ncorners,
6296  & cell_corner_lat, cell_corner_lon,
6297  & cell_center_lat, cell_center_lon,
6298  & cell_grid_num,
6299  & inpoly, lboundary, bedgeid)
6300 
6301  if (inpoly) then
6302  lat1 = midlat
6303  lon1 = midlon
6304  else
6305  lat2 = midlat
6306  lon2 = midlon
6307  endif
6308 
6309  if (abs(lat1-lat2) < tiny .and.
6310  & abs(lon1-lon2) < tiny .and. lboundary) then
6311  converged = .true.
6312  endif
6313 
6314  if (it > 100) then
6315  exit
6316  endif
6317 
6318  it = it + 1
6319  enddo ! do while (not converged)
6320 
6321 
6322  bpt_lat = midlat
6323  bpt_lon = midlon
6324 
6325  end subroutine converge_to_bdry
6326 
6327 
6328 
6329 
6330  end module scrip_remap_conservative
6331 
6332 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
scrip_remap_conservative::srch_corners_locate_point
integer(scrip_i4), save srch_corners_locate_point
Definition: scrip_remap_conservative.f:130
scrip_remap_conservative::avoid_pole_offset
real(scrip_r8), save avoid_pole_offset
Definition: scrip_remap_conservative.f:83
scrip_grids::grid2_centroid_lat
real(scrip_r8), dimension(:), allocatable, target, save grid2_centroid_lat
Definition: scrip_grids.f:103
scrip_remap_vars::norm_opt_none
integer(scrip_i4), parameter norm_opt_none
Definition: scrip_remap_vars.f:54
scrip_grids::grid1_corners
integer(scrip_i4), save grid1_corners
Definition: scrip_grids.f:68
scrip_remap_vars::grid1_add_map1
integer(scrip_i4), dimension(:), allocatable, save grid1_add_map1
Definition: scrip_remap_vars.f:77
scrip_constants::half
real(kind=scrip_r8), parameter half
Definition: scrip_constants.f:50
scrip_grids::grid1_mask
logical(scrip_logical), dimension(:), allocatable, target, save grid1_mask
Definition: scrip_grids.f:93
scrip_remap_conservative::srch_corners_locate_segstart
integer(scrip_i4), save srch_corners_locate_segstart
Definition: scrip_remap_conservative.f:102
scrip_remap_conservative::srch_center_lat_locate_point
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_locate_point
Definition: scrip_remap_conservative.f:144
scrip_remap_conservative::srch_corner_lon_locate_segstart
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_locate_segstart
Definition: scrip_remap_conservative.f:111
scrip_remap_conservative::ptinpolygen
subroutine ptinpolygen(ptx, pty, ncorners, cell_corner_x, cell_corner_y, cell_center_x, cell_center_y, latlon, inpoly, lboundary, edgeid)
Definition: scrip_remap_conservative.f:5211
scrip_remap_conservative::srch_center_lon_locate_segstart
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_locate_segstart
Definition: scrip_remap_conservative.f:117
scrip_remap_conservative::last_cell_grid_num_get_srch_cells
integer(scrip_i4), save last_cell_grid_num_get_srch_cells
Definition: scrip_remap_conservative.f:166
scrip_grids::grid1_centroid_lat
real(scrip_r8), dimension(:), allocatable, target, save grid1_centroid_lat
Definition: scrip_grids.f:103
scrip_remap_conservative::find_adj_cell
subroutine find_adj_cell(cell_add, edge_id, cell_grid_num, adj_add)
Definition: scrip_remap_conservative.f:5987
scrip_grids::grid1_spole_cell
integer(scrip_i4), save grid1_spole_cell
Definition: scrip_grids.f:135
scrip_remap_vars::wt_lowest
real(kind=scrip_r8) wt_lowest
Definition: scrip_remap_vars.f:88
scrip_remap_vars::num_maps
integer(scrip_i4), save num_maps
Definition: scrip_remap_vars.f:66
scrip_remap_conservative::last_cell_locate_segstart
integer(scrip_i4), save last_cell_locate_segstart
Definition: scrip_remap_conservative.f:96
scrip_remap_vars::wt_highest
real(kind=scrip_r8) wt_highest
Definition: scrip_remap_vars.f:88
scrip_remap_conservative::converge_to_bdry
subroutine converge_to_bdry(cell_add, cell_grid_num, ncorners, cell_corner_lat, cell_corner_lon, cell_center_lat, cell_center_lon, inpt_lat, inpt_lon, outpt_lat, outpt_lon, bpt_lat, bpt_lon, bedgeid)
Definition: scrip_remap_conservative.f:6217
scrip_grids::luse_grid2_area
logical(scrip_logical), save luse_grid2_area
Definition: scrip_grids.f:126
scrip_remap_conservative::first_call_locate_point
logical(scrip_logical), save first_call_locate_point
Definition: scrip_remap_conservative.f:121
scrip_remap_conservative::link_add2
integer(scrip_i4), dimension(:,:), allocatable, save link_add2
Definition: scrip_remap_conservative.f:86
scrip_timers::timer_start
subroutine timer_start(timer)
Definition: scrip_timers.f:216
scrip_grids::grid2_centroid_lon
real(scrip_r8), dimension(:), allocatable, target, save grid2_centroid_lon
Definition: scrip_grids.f:103
scrip_remap_conservative::srch_corner_lon_loc_get_srch_cells
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_loc_get_srch_cells
Definition: scrip_remap_conservative.f:158
scrip_grids::num_srch_bins
integer(scrip_i4), save num_srch_bins
Definition: scrip_grids.f:152
scrip_remap_conservative::ptincell
subroutine ptincell(ptlat, ptlon, cell_add, ncorners, cell_corner_lat, cell_corner_lon, cell_center_lat, cell_center_lon, cell_grid_id, inpoly, lboundary, edgeid)
Definition: scrip_remap_conservative.f:4608
scrip_remap_conservative::last_srch_grid_num_locate_segstart
integer(scrip_i4), save last_srch_grid_num_locate_segstart
Definition: scrip_remap_conservative.f:96
scrip_remap_conservative::store_link_cnsrv
subroutine store_link_cnsrv(add1, add2, weights)
Definition: scrip_remap_conservative.f:3951
scrip_remap_conservative::last_cell_grid_num_locate_point
integer(scrip_i4), save last_cell_grid_num_locate_point
Definition: scrip_remap_conservative.f:124
scrip_remap_vars::num_links_map1
integer(scrip_i4), save num_links_map1
Definition: scrip_remap_vars.f:66
scrip_timers::timer_print
subroutine timer_print(timer)
Definition: scrip_timers.f:176
scrip_remap_conservative::remap_conserv
subroutine remap_conserv(l_master, l_test)
Definition: scrip_remap_conservative.f:246
scrip_grids::grid2_mask
logical(scrip_logical), dimension(:), allocatable, target, save grid2_mask
Definition: scrip_grids.f:93
scrip_remap_vars::norm_opt
integer(scrip_i4), save norm_opt
Definition: scrip_remap_vars.f:66
scrip_grids::grid2_center_lat
real(scrip_r8), dimension(:), allocatable, target, save grid2_center_lat
Definition: scrip_grids.f:103
scrip_remap_conservative::cell_integrate
subroutine cell_integrate(cell_add, grid_num, phi_or_theta)
Definition: scrip_remap_conservative.f:1005
scrip_constants::pi2
real(kind=scrip_r8), parameter pi2
Definition: scrip_constants.f:50
scrip_remap_conservative::srch_add_find_adj_cell
integer(scrip_i4), dimension(:), allocatable, save srch_add_find_adj_cell
Definition: scrip_remap_conservative.f:186
scrip_remap_conservative::cellblock_integrate
subroutine cellblock_integrate(ibegin, iend, grid_num, phi_or_theta)
Definition: scrip_remap_conservative.f:982
scrip_remap_vars::max_links_map2
integer(scrip_i4), save max_links_map2
Definition: scrip_remap_vars.f:66
scrip_remap_conservative::last_cell_locate_point
integer(scrip_i4), save last_cell_locate_point
Definition: scrip_remap_conservative.f:124
scrip_remap_conservative::first_call_locate_segstart
logical(scrip_logical), save first_call_locate_segstart
Definition: scrip_remap_conservative.f:93
scrip_remap_conservative::srch_center_lon_locate_point
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_locate_point
Definition: scrip_remap_conservative.f:144
scrip_constants::one
real(kind=scrip_r8), parameter one
Definition: scrip_constants.f:50
scrip_remap_vars::num_links_map2
integer(scrip_i4), save num_links_map2
Definition: scrip_remap_vars.f:66
scrip_grids::grid1_area
real(scrip_r8), dimension(:), allocatable, target, save grid1_area
Definition: scrip_grids.f:103
scrip_grids::north_thresh
real(scrip_r8), save north_thresh
Definition: scrip_grids.f:169
scrip_grids::grid1_npole_cell
integer(scrip_i4), save grid1_npole_cell
Definition: scrip_grids.f:135
scrip_remap_conservative::srch_corner_lat_locate_point
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_locate_point
Definition: scrip_remap_conservative.f:138
scrip_remap_conservative::get_srch_cells
subroutine get_srch_cells(cell_add, cell_grid_num, srch_grid_num, num_srch_cells, srch_add, srch_corners, srch_corner_lat, srch_corner_lon, srch_center_lat, srch_center_lon)
Definition: scrip_remap_conservative.f:5557
scrip_grids
Definition: scrip_grids.f:49
scrip_constants::tiny
real(kind=scrip_r8), parameter tiny
Definition: scrip_constants.f:50
scrip_remap_conservative::srch_corner_lon_find_adj_cell
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_find_adj_cell
Definition: scrip_remap_conservative.f:188
scrip_grids::special_polar_cell2
logical(scrip_logical), dimension(:), allocatable, target, save special_polar_cell2
Definition: scrip_grids.f:93
scrip_remap_conservative::last_srch_grid_num_get_srch_cells
integer(scrip_i4), save last_srch_grid_num_get_srch_cells
Definition: scrip_remap_conservative.f:166
scrip_grids::grid2_corners
integer(scrip_i4), save grid2_corners
Definition: scrip_grids.f:68
scrip_grids::grid1_area_in
real(scrip_r8), dimension(:), allocatable, target, save grid1_area_in
Definition: scrip_grids.f:103
scrip_remap_conservative::num_srch_cell_locate_points
integer(scrip_i4), save num_srch_cell_locate_points
Definition: scrip_remap_conservative.f:130
scrip_constants::zero
real(kind=scrip_r8), parameter zero
Definition: scrip_constants.f:50
scrip_kindsmod::scrip_r8
integer, parameter, public scrip_r8
Definition: scrip_kindsmod.f90:38
scrip_remap_conservative::locate_segstart
subroutine locate_segstart(cell_grid_num, cell, beglat, beglon, endlat, endlon, offset, srch_grid_num, cont_cell, lboundary, edgeid)
Definition: scrip_remap_conservative.f:4096
scrip_remap_conservative::last_cell_find_adj_cell
integer(scrip_i4), save last_cell_find_adj_cell
Definition: scrip_remap_conservative.f:180
scrip_constants::quart
real(kind=scrip_r8), parameter quart
Definition: scrip_constants.f:50
scrip_grids::grid1_centroid_lon
real(scrip_r8), dimension(:), allocatable, target, save grid1_centroid_lon
Definition: scrip_grids.f:103
scrip_grids::grid1_bound_box
real(scrip_r8), dimension(:,:), allocatable, target, save grid1_bound_box
Definition: scrip_grids.f:131
scrip_remap_vars::grid2_add_map2
integer(scrip_i4), dimension(:), allocatable, save grid2_add_map2
Definition: scrip_remap_vars.f:77
scrip_grids::grid2_corner_lat
real(scrip_r8), dimension(:,:), allocatable, target, save grid2_corner_lat
Definition: scrip_grids.f:120
scrip_remap_vars::norm_opt_frcarea
integer(scrip_i4), parameter norm_opt_frcarea
Definition: scrip_remap_vars.f:54
scrip_remap_vars::max_links_map1
integer(scrip_i4), save max_links_map1
Definition: scrip_remap_vars.f:66
scrip_remap_conservative::num_srch_cells_locate_segstart
integer(scrip_i4), save num_srch_cells_locate_segstart
Definition: scrip_remap_conservative.f:102
scrip_remap_conservative::srch_center_lat_find_adj_cell
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_find_adj_cell
Definition: scrip_remap_conservative.f:191
scrip_grids::luse_grid1_area
logical(scrip_logical), save luse_grid1_area
Definition: scrip_grids.f:126
scrip_constants::pih
real(kind=scrip_r8), parameter pih
Definition: scrip_constants.f:50
scrip_grids::grid1_center_lon
real(scrip_r8), dimension(:), allocatable, target, save grid1_center_lon
Definition: scrip_grids.f:103
scrip_constants::two
real(kind=scrip_r8), parameter two
Definition: scrip_constants.f:50
scrip_grids::grid2_npole_cell
integer(scrip_i4), save grid2_npole_cell
Definition: scrip_grids.f:135
scrip_grids::grid2_area
real(scrip_r8), dimension(:), allocatable, target, save grid2_area
Definition: scrip_grids.f:103
scrip_remap_conservative::intersection
subroutine intersection(seg_cell_id, seg_grid_id, beglat, beglon, endlat, endlon, begseg, begedge, cell_id, ncorners, cell_corner_lat, cell_corner_lon, cell_grid_id, intrsct_lat, intrsct_lon, intedge, sinang2, lcoinc, lthresh)
Definition: scrip_remap_conservative.f:2196
scrip_grids::grid1_corner_lat
real(scrip_r8), dimension(:,:), allocatable, target, save grid1_corner_lat
Definition: scrip_grids.f:120
scrip_remap_conservative::srch_center_lon_loc_get_srch_cells
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_loc_get_srch_cells
Definition: scrip_remap_conservative.f:162
scrip_remap_vars::frac_highest
real(kind=scrip_r8) frac_highest
Definition: scrip_remap_vars.f:87
scrip_remap_conservative::srch_corner_lat_find_adj_cell
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_find_adj_cell
Definition: scrip_remap_conservative.f:188
scrip_grids::grid2_size
integer(scrip_i4), save grid2_size
Definition: scrip_grids.f:68
scrip_grids::bin_addr2
integer(scrip_i4), dimension(:,:), allocatable, save bin_addr2
Definition: scrip_grids.f:155
scrip_remap_vars::frac_lowest
real(kind=scrip_r8) frac_lowest
Definition: scrip_remap_vars.f:87
scrip_remap_conservative::last_cell_grid_num_locate_segstart
integer(scrip_i4), save last_cell_grid_num_locate_segstart
Definition: scrip_remap_conservative.f:96
scrip_remap_conservative::ptinpoly
subroutine ptinpoly(ptx, pty, ncorners, cell_corner_x, cell_corner_y, latlon, inpoly, lboundary, edgeid)
Definition: scrip_remap_conservative.f:4815
scrip_remap_conservative::nthreads
integer(scrip_i4) nthreads
Definition: scrip_remap_conservative.f:73
scrip_remap_conservative::last_cell_add_get_srch_cells
integer(scrip_i4), save last_cell_add_get_srch_cells
Definition: scrip_remap_conservative.f:166
scrip_remap_vars::wts_map1
real(scrip_r8), dimension(:,:), allocatable, save wts_map1
Definition: scrip_remap_vars.f:83
scrip_remap_conservative::link_add1
integer(scrip_i4), dimension(:,:), allocatable, save link_add1
Definition: scrip_remap_conservative.f:86
scrip_remap_conservative::srch_add_locate_point
integer(scrip_i4), dimension(:), allocatable, save srch_add_locate_point
Definition: scrip_remap_conservative.f:134
scrip_remap_vars::resize_remap_vars
subroutine resize_remap_vars(nmap, increment)
Definition: scrip_remap_vars.f:173
scrip_remap_conservative::line_integral
subroutine line_integral(phi_or_theta, weights, num_wts, in_phi1, in_phi2, theta1, theta2, grid1_lat, grid1_lon, grid2_lat, grid2_lon)
Definition: scrip_remap_conservative.f:3460
scrip_constants
Definition: scrip_constants.f:38
scrip_remap_conservative::last_cell_grid_num_find_adj_cell
integer(scrip_i4), save last_cell_grid_num_find_adj_cell
Definition: scrip_remap_conservative.f:180
scrip_grids::grid2_area_in
real(scrip_r8), dimension(:), allocatable, target, save grid2_area_in
Definition: scrip_grids.f:103
scrip_remap_conservative
Definition: scrip_remap_conservative.f:60
scrip_remap_conservative::num_srch_cells_loc_get_srch_cells
integer(scrip_i4), save num_srch_cells_loc_get_srch_cells
Definition: scrip_remap_conservative.f:148
scrip_grids::grid1_center_lat
real(scrip_r8), dimension(:), allocatable, target, save grid1_center_lat
Definition: scrip_grids.f:103
scrip_remap_conservative::srch_corner_lat_loc_get_srch_cells
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_loc_get_srch_cells
Definition: scrip_remap_conservative.f:158
scrip_kindsmod
Definition: scrip_kindsmod.f90:3
scrip_remap_conservative::num_srch_cells_find_adj_cell
integer(scrip_i4), save num_srch_cells_find_adj_cell
Definition: scrip_remap_conservative.f:180
scrip_remap_conservative::modify_polar_cell
subroutine modify_polar_cell(ncorners, nalloc, cell_corner_lat, cell_corner_lon)
Definition: scrip_remap_conservative.f:2082
scrip_remap_conservative::ptinpolarpoly
subroutine ptinpolarpoly(ptx, pty, ncorners, cell_corner_x, cell_corner_y, latlon, whichpole, inpoly, lboundary, edgeid)
Definition: scrip_remap_conservative.f:5042
scrip_grids::grid1_corner_lon
real(scrip_r8), dimension(:,:), allocatable, target, save grid1_corner_lon
Definition: scrip_grids.f:120
scrip_remap_conservative::srch_corner_lon_locate_point
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_locate_point
Definition: scrip_remap_conservative.f:138
scrip_grids::npseg
integer(scrip_i4), save npseg
Definition: scrip_grids.f:178
scrip_remap_vars::grid1_add_map2
integer(scrip_i4), dimension(:), allocatable, save grid1_add_map2
Definition: scrip_remap_vars.f:77
scrip_remap_vars::grid2_add_map1
integer(scrip_i4), dimension(:), allocatable, save grid2_add_map1
Definition: scrip_remap_vars.f:77
scrip_remap_conservative::pole_intersection
subroutine pole_intersection(location, ncorners, cell_corners_lat, cell_corners_lon, cell_grid_id, beglat, beglon, endlat, endlon, begseg, begedge, intedge, intrsct_lat, intrsct_lon, sinang2, lcoinc, lthresh)
Definition: scrip_remap_conservative.f:2728
scrip_grids::south_thresh
real(scrip_r8), save south_thresh
Definition: scrip_grids.f:169
scrip_remap_conservative::srch_corner_lat_locate_segstart
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_locate_segstart
Definition: scrip_remap_conservative.f:111
scrip_grids::grid1_size
integer(scrip_i4), save grid1_size
Definition: scrip_grids.f:68
scrip_remap_conservative::first_call_find_adj_cell
logical(scrip_logical), save first_call_find_adj_cell
Definition: scrip_remap_conservative.f:174
scrip_remap_conservative::srch_center_lon_find_adj_cell
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_find_adj_cell
Definition: scrip_remap_conservative.f:191
scrip_remap_vars::wts_map2
real(scrip_r8), dimension(:,:), allocatable, save wts_map2
Definition: scrip_remap_vars.f:83
scrip_grids::grid2_bound_box
real(scrip_r8), dimension(:,:), allocatable, target, save grid2_bound_box
Definition: scrip_grids.f:131
scrip_timers
Definition: scrip_timers.f:39
scrip_grids::grid2_corner_lon
real(scrip_r8), dimension(:,:), allocatable, target, save grid2_corner_lon
Definition: scrip_grids.f:120
scrip_grids::grid1_frac
real(scrip_r8), dimension(:), allocatable, target, save grid1_frac
Definition: scrip_grids.f:103
scrip_remap_conservative::srch_corners_find_adj_cell
integer(scrip_i4), save srch_corners_find_adj_cell
Definition: scrip_remap_conservative.f:180
scrip_remap_vars
Definition: scrip_remap_vars.f:40
scrip_remap_conservative::srch_center_lat_loc_get_srch_cells
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_loc_get_srch_cells
Definition: scrip_remap_conservative.f:162
scrip_remap_conservative::srch_corners_loc_get_srch_cells
integer(scrip_i4), save srch_corners_loc_get_srch_cells
Definition: scrip_remap_conservative.f:148
scrip_remap_conservative::srch_add_loc_get_srch_cells
integer(scrip_i4), dimension(:), allocatable, save srch_add_loc_get_srch_cells
Definition: scrip_remap_conservative.f:154
scrip_remap_conservative::line_integral_theta
subroutine line_integral_theta(weights, num_wts, in_phi1, in_phi2, theta1, theta2, grid1_lat, grid1_lon, grid2_lat, grid2_lon)
Definition: scrip_remap_conservative.f:3682
scrip_remap_conservative::line_integral_phi
subroutine line_integral_phi(weights, num_wts, in_phi1, in_phi2, theta1, theta2, grid1_lat, grid1_lon, grid2_lat, grid2_lon)
Definition: scrip_remap_conservative.f:3522
scrip_remap_vars::norm_opt_dstarea
integer(scrip_i4), parameter norm_opt_dstarea
Definition: scrip_remap_vars.f:54
scrip_remap_conservative::ptinpolygen2
subroutine ptinpolygen2(ptx, pty, ncorners, cell_corner_x, cell_corner_y, latlon, inpoly, lboundary, edgeid)
Definition: scrip_remap_conservative.f:5378
scrip_remap_conservative::srch_add_locate_segstart
integer(scrip_i4), dimension(:), allocatable, save srch_add_locate_segstart
Definition: scrip_remap_conservative.f:107
scrip_grids::grid2_spole_cell
integer(scrip_i4), save grid2_spole_cell
Definition: scrip_grids.f:135
scrip_remap_conservative::avoid_pole_count
integer(scrip_i4), save avoid_pole_count
Definition: scrip_remap_conservative.f:80
scrip_remap_conservative::last_srch_grid_num_locate_point
integer(scrip_i4), save last_srch_grid_num_locate_point
Definition: scrip_remap_conservative.f:124
scrip_constants::pi
real(kind=scrip_r8), parameter pi
Definition: scrip_constants.f:50
scrip_remap_conservative::first_call_store_link_cnsrv
logical(scrip_logical), save first_call_store_link_cnsrv
Definition: scrip_remap_conservative.f:90
scrip_grids::special_polar_cell1
logical(scrip_logical), dimension(:), allocatable, target, save special_polar_cell1
Definition: scrip_grids.f:93
scrip_remap_vars::num_wts
integer(scrip_i4), save num_wts
Definition: scrip_remap_vars.f:66
scrip_remap_conservative::first_call_get_srch_cells
logical(scrip_logical), save first_call_get_srch_cells
Definition: scrip_remap_conservative.f:171
scrip_timers::timer_stop
subroutine timer_stop(timer)
Definition: scrip_timers.f:250
scrip_grids::grid2_frac
real(scrip_r8), dimension(:), allocatable, target, save grid2_frac
Definition: scrip_grids.f:103
scrip_remap_conservative::locate_point
subroutine locate_point(ptlat, ptlon, cell, cell_grid_num, srch_grid_num, cont_cell, lboundary, edgeid)
Definition: scrip_remap_conservative.f:4423
scrip_grids::bin_addr1
integer(scrip_i4), dimension(:,:), allocatable, save bin_addr1
Definition: scrip_grids.f:155
scrip_grids::grid2_center_lon
real(scrip_r8), dimension(:), allocatable, target, save grid2_center_lon
Definition: scrip_grids.f:103
scrip_remap_conservative::srch_center_lat_locate_segstart
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_locate_segstart
Definition: scrip_remap_conservative.f:117
scrip_remap_vars::resize_increment
integer(scrip_i4), save resize_increment
Definition: scrip_remap_vars.f:66