!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! this module contains necessary routines for computing addresses ! and weights for a conservative interpolation between any two ! grids on a sphere. the weights are computed by performing line ! integrals around all overlap regions of the two grids. see ! Dukowicz and Kodis, SIAM J. Sci. Stat. Comput. 8, 305 (1987) and ! Jones, P.W. Monthly Weather Review (submitted). ! !----------------------------------------------------------------------- ! ! CVS:$Id$ ! ! Copyright (c) 1997, 1998 the Regents of the University of ! California. ! ! This software and ancillary information (herein called software) ! called SCRIP is made available under the terms described here. ! The software has been approved for release with associated ! LA-CC Number 98-45. ! ! Unless otherwise indicated, this software has been authored ! by an employee or employees of the University of California, ! operator of the Los Alamos National Laboratory under Contract ! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S. ! Government has rights to use, reproduce, and distribute this ! software. The public may copy and use this software without ! charge, provided that this Notice and any statement of authorship ! are reproduced on all copies. Neither the Government nor the ! University makes any warranty, express or implied, or assumes ! any liability or responsibility for the use of this software. ! ! If software is modified to produce derivative works, such modified ! software should be clearly marked, so as not to confuse it with ! the version available from Los Alamos National Laboratory. ! !*********************************************************************** module remap_conservative !----------------------------------------------------------------------- use kinds_mod ! defines common data types use constants ! defines common constants use timers ! module for timing use grids ! module containing grid information use remap_vars ! module containing remap information implicit none !----------------------------------------------------------------------- ! ! module variables ! !----------------------------------------------------------------------- integer (kind=int_kind), save :: & num_srch_cells ! num cells in restricted search arrays integer (kind=int_kind), dimension(:), allocatable, save :: & srch_add ! global address of cells in srch arrays real (kind=dbl_kind), parameter :: & north_thresh = 1.45_dbl_kind, & ! threshold for coord transf. south_thresh =-2.00_dbl_kind ! threshold for coord transf. real (kind=dbl_kind), dimension(:,:), allocatable, save :: & srch_corner_lat, & ! lat of each corner of srch cells srch_corner_lon ! lon of each corner of srch cells !*********************************************************************** contains !*********************************************************************** subroutine remap_conserv !----------------------------------------------------------------------- ! ! this routine traces the perimeters of every grid cell on each ! grid checking for intersections with the other grid and computing ! line integrals for each subsegment. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (kind=int_kind), parameter :: & max_subseg = 10000 ! max number of subsegments per segment ! to prevent infinite loop integer (kind=int_kind) :: & grid1_add, & ! current linear address for grid1 cell grid2_add, & ! current linear address for grid2 cell min_add, & ! addresses for restricting search of max_add, & ! destination grid n, nwgt, & ! generic counters corner, & ! corner of cell that segment starts from next_corn, & ! corner of cell that segment ends on num_subseg ! number of subsegments logical (kind=log_kind) :: & lcoinc, & ! flag for coincident segments lrevers, & ! flag for reversing direction of segment lbegin ! flag for first integration of a segment logical (kind=log_kind), dimension(:), allocatable :: & srch_mask ! mask for restricting searches real (kind=dbl_kind) :: & intrsct_lat, intrsct_lon, & ! lat/lon of next intersect beglat, endlat, beglon, endlon, & ! endpoints of current seg. norm_factor ! factor for normalizing wts real (kind=dbl_kind), dimension(:), allocatable :: & grid2_centroid_lat, grid2_centroid_lon, & ! centroid coords grid1_centroid_lat, grid1_centroid_lon ! on each grid real (kind=dbl_kind), dimension(2) :: begseg ! begin lat/lon for ! full segment real (kind=dbl_kind), dimension(6) :: weights ! local wgt array !----------------------------------------------------------------------- ! ! initialize centroid arrays ! !----------------------------------------------------------------------- allocate( grid1_centroid_lat(grid1_size), & grid1_centroid_lon(grid1_size), & grid2_centroid_lat(grid2_size), & grid2_centroid_lon(grid2_size)) grid1_centroid_lat = zero grid1_centroid_lon = zero grid2_centroid_lat = zero grid2_centroid_lon = zero !----------------------------------------------------------------------- ! ! integrate around each cell on grid1 ! !----------------------------------------------------------------------- allocate(srch_mask(grid2_size)) print *,'grid1 sweep ' do grid1_add = 1,grid1_size !*** !*** restrict searches first using search bins !*** call timer_start(1) min_add = grid2_size max_add = 1 do n=1,num_srch_bins if (grid1_add >= bin_addr1(1,n) .and. & grid1_add <= bin_addr1(2,n)) then min_add = min(min_add, bin_addr2(1,n)) max_add = max(max_add, bin_addr2(2,n)) endif end do !*** !*** further restrict searches using bounding boxes !*** num_srch_cells = 0 do grid2_add = min_add,max_add srch_mask(grid2_add) = (grid2_bound_box(1,grid2_add) <= & grid1_bound_box(2,grid1_add)) .and. & (grid2_bound_box(2,grid2_add) >= & grid1_bound_box(1,grid1_add)) .and. & (grid2_bound_box(3,grid2_add) <= & grid1_bound_box(4,grid1_add)) .and. & (grid2_bound_box(4,grid2_add) >= & grid1_bound_box(3,grid1_add)) if (srch_mask(grid2_add)) num_srch_cells = num_srch_cells+1 end do !*** !*** create search arrays !*** allocate(srch_add(num_srch_cells), & srch_corner_lat(grid2_corners,num_srch_cells), & srch_corner_lon(grid2_corners,num_srch_cells)) n = 0 gather1: do grid2_add = min_add,max_add if (srch_mask(grid2_add)) then n = n+1 srch_add(n) = grid2_add srch_corner_lat(:,n) = grid2_corner_lat(:,grid2_add) srch_corner_lon(:,n) = grid2_corner_lon(:,grid2_add) endif end do gather1 call timer_stop(1) !*** !*** integrate around this cell !*** do corner = 1,grid1_corners next_corn = mod(corner,grid1_corners) + 1 !*** !*** define endpoints of the current segment !*** beglat = grid1_corner_lat(corner,grid1_add) beglon = grid1_corner_lon(corner,grid1_add) endlat = grid1_corner_lat(next_corn,grid1_add) endlon = grid1_corner_lon(next_corn,grid1_add) lrevers = .false. !*** !*** to ensure exact path taken during both !*** sweeps, always integrate segments in the same !*** direction (SW to NE). !*** if ((endlat < beglat) .or. & (endlat == beglat .and. endlon < beglon)) then beglat = grid1_corner_lat(next_corn,grid1_add) beglon = grid1_corner_lon(next_corn,grid1_add) endlat = grid1_corner_lat(corner,grid1_add) endlon = grid1_corner_lon(corner,grid1_add) lrevers = .true. endif begseg(1) = beglat begseg(2) = beglon lbegin = .true. num_subseg = 0 !*** !*** if this is a constant-longitude segment, skip the rest !*** since the line integral contribution will be zero. !*** if (endlon /= beglon) then !*** !*** integrate along this segment, detecting intersections !*** and computing the line integral for each sub-segment !*** do while (beglat /= endlat .or. beglon /= endlon) !*** !*** prevent infinite loops if integration gets stuck !*** near cell or threshold boundary !*** num_subseg = num_subseg + 1 if (num_subseg > max_subseg) then stop 'integration stalled: num_subseg exceeded limit' endif !*** !*** find next intersection of this segment with a grid !*** line on grid 2. !*** call timer_start(2) call intersection(grid2_add,intrsct_lat,intrsct_lon,lcoinc, & beglat, beglon, endlat, endlon, begseg, & lbegin, lrevers) call timer_stop(2) lbegin = .false. !*** !*** compute line integral for this subsegment. !*** call timer_start(3) if (grid2_add /= 0) then call line_integral(weights, num_wts, & beglon, intrsct_lon, beglat, intrsct_lat, & grid1_center_lat(grid1_add), & grid1_center_lon(grid1_add), & grid2_center_lat(grid2_add), & grid2_center_lon(grid2_add)) else call line_integral(weights, num_wts, & beglon, intrsct_lon, beglat, intrsct_lat, & grid1_center_lat(grid1_add), & grid1_center_lon(grid1_add), & grid1_center_lat(grid1_add), & grid1_center_lon(grid1_add)) endif call timer_stop(3) !*** !*** if integrating in reverse order, change !*** sign of weights !*** if (lrevers) then weights = -weights endif !*** !*** store the appropriate addresses and weights. !*** also add contributions to cell areas and centroids. !*** !if (grid1_add == 119247) then ! print *,grid1_add,grid2_add,corner,weights(1) ! print *,grid1_corner_lat(:,grid1_add) ! print *,grid1_corner_lon(:,grid1_add) ! print *,grid2_corner_lat(:,grid2_add) ! print *,grid2_corner_lon(:,grid2_add) ! print *,beglat,beglon,intrsct_lat,intrsct_lon !endif if (grid2_add /= 0) then if (grid1_mask(grid1_add)) then call timer_start(4) call store_link_cnsrv(grid1_add, grid2_add, weights) call timer_stop(4) grid1_frac(grid1_add) = grid1_frac(grid1_add) + & weights(1) grid2_frac(grid2_add) = grid2_frac(grid2_add) + & weights(num_wts+1) endif endif grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1) grid1_centroid_lat(grid1_add) = & grid1_centroid_lat(grid1_add) + weights(2) grid1_centroid_lon(grid1_add) = & grid1_centroid_lon(grid1_add) + weights(3) !*** !*** reset beglat and beglon for next subsegment. !*** beglat = intrsct_lat beglon = intrsct_lon end do endif !*** !*** end of segment !*** end do !*** !*** finished with this cell: deallocate search array and !*** start on next cell deallocate(srch_add, srch_corner_lat, srch_corner_lon) end do deallocate(srch_mask) !----------------------------------------------------------------------- ! ! integrate around each cell on grid2 ! !----------------------------------------------------------------------- allocate(srch_mask(grid1_size)) print *,'grid2 sweep ' do grid2_add = 1,grid2_size !*** !*** restrict searches first using search bins !*** call timer_start(5) min_add = grid1_size max_add = 1 do n=1,num_srch_bins if (grid2_add >= bin_addr2(1,n) .and. & grid2_add <= bin_addr2(2,n)) then min_add = min(min_add, bin_addr1(1,n)) max_add = max(max_add, bin_addr1(2,n)) endif end do !*** !*** further restrict searches using bounding boxes !*** num_srch_cells = 0 do grid1_add = min_add, max_add srch_mask(grid1_add) = (grid1_bound_box(1,grid1_add) <= & grid2_bound_box(2,grid2_add)) .and. & (grid1_bound_box(2,grid1_add) >= & grid2_bound_box(1,grid2_add)) .and. & (grid1_bound_box(3,grid1_add) <= & grid2_bound_box(4,grid2_add)) .and. & (grid1_bound_box(4,grid1_add) >= & grid2_bound_box(3,grid2_add)) if (srch_mask(grid1_add)) num_srch_cells = num_srch_cells+1 end do allocate(srch_add(num_srch_cells), & srch_corner_lat(grid1_corners,num_srch_cells), & srch_corner_lon(grid1_corners,num_srch_cells)) n = 0 gather2: do grid1_add = min_add,max_add if (srch_mask(grid1_add)) then n = n+1 srch_add(n) = grid1_add srch_corner_lat(:,n) = grid1_corner_lat(:,grid1_add) srch_corner_lon(:,n) = grid1_corner_lon(:,grid1_add) endif end do gather2 call timer_stop(5) !*** !*** integrate around this cell !*** do corner = 1,grid2_corners next_corn = mod(corner,grid2_corners) + 1 beglat = grid2_corner_lat(corner,grid2_add) beglon = grid2_corner_lon(corner,grid2_add) endlat = grid2_corner_lat(next_corn,grid2_add) endlon = grid2_corner_lon(next_corn,grid2_add) lrevers = .false. !*** !*** to ensure exact path taken during both !*** sweeps, always integrate in the same direction !*** if ((endlat < beglat) .or. & (endlat == beglat .and. endlon < beglon)) then beglat = grid2_corner_lat(next_corn,grid2_add) beglon = grid2_corner_lon(next_corn,grid2_add) endlat = grid2_corner_lat(corner,grid2_add) endlon = grid2_corner_lon(corner,grid2_add) lrevers = .true. endif begseg(1) = beglat begseg(2) = beglon lbegin = .true. !*** !*** if this is a constant-longitude segment, skip the rest !*** since the line integral contribution will be zero. !*** if (endlon /= beglon) then num_subseg = 0 !*** !*** integrate along this segment, detecting intersections !*** and computing the line integral for each sub-segment !*** do while (beglat /= endlat .or. beglon /= endlon) !*** !*** prevent infinite loops if integration gets stuck !*** near cell or threshold boundary !*** num_subseg = num_subseg + 1 if (num_subseg > max_subseg) then stop 'integration stalled: num_subseg exceeded limit' endif !*** !*** find next intersection of this segment with a line !*** on grid 2. !*** call timer_start(6) call intersection(grid1_add,intrsct_lat,intrsct_lon,lcoinc, & beglat, beglon, endlat, endlon, begseg, & lbegin, lrevers) call timer_stop(6) lbegin = .false. !*** !*** compute line integral for this subsegment. !*** call timer_start(7) if (grid1_add /= 0) then call line_integral(weights, num_wts, & beglon, intrsct_lon, beglat, intrsct_lat, & grid1_center_lat(grid1_add), & grid1_center_lon(grid1_add), & grid2_center_lat(grid2_add), & grid2_center_lon(grid2_add)) else call line_integral(weights, num_wts, & beglon, intrsct_lon, beglat, intrsct_lat, & grid2_center_lat(grid2_add), & grid2_center_lon(grid2_add), & grid2_center_lat(grid2_add), & grid2_center_lon(grid2_add)) endif call timer_stop(7) if (lrevers) then weights = -weights endif !*** !*** store the appropriate addresses and weights. !*** also add contributions to cell areas and centroids. !*** if there is a coincidence, do not store weights !*** because they have been captured in the previous loop. !*** the grid1 mask is the master mask !*** !if (grid1_add == 119247) then ! print *,grid1_add,grid2_add,corner,weights(1) ! print *,grid1_corner_lat(:,grid1_add) ! print *,grid1_corner_lon(:,grid1_add) ! print *,grid2_corner_lat(:,grid2_add) ! print *,grid2_corner_lon(:,grid2_add) ! print *,beglat,beglon,intrsct_lat,intrsct_lon !endif if (.not. lcoinc .and. grid1_add /= 0) then if (grid1_mask(grid1_add)) then call timer_start(8) call store_link_cnsrv(grid1_add, grid2_add, weights) call timer_stop(8) grid1_frac(grid1_add) = grid1_frac(grid1_add) + & weights(1) grid2_frac(grid2_add) = grid2_frac(grid2_add) + & weights(num_wts+1) endif endif grid2_area(grid2_add) = grid2_area(grid2_add) + & weights(num_wts+1) grid2_centroid_lat(grid2_add) = & grid2_centroid_lat(grid2_add) + weights(num_wts+2) grid2_centroid_lon(grid2_add) = & grid2_centroid_lon(grid2_add) + weights(num_wts+3) !*** !*** reset beglat and beglon for next subsegment. !*** beglat = intrsct_lat beglon = intrsct_lon end do endif !*** !*** end of segment !*** end do !*** !*** finished with this cell: deallocate search array and !*** start on next cell deallocate(srch_add, srch_corner_lat, srch_corner_lon) end do deallocate(srch_mask) !----------------------------------------------------------------------- ! ! correct for situations where N/S pole not explicitly included in ! grid (i.e. as a grid corner point). if pole is missing from only ! one grid, need to correct only the area and centroid of that ! grid. if missing from both, do complete weight calculation. ! !----------------------------------------------------------------------- !*** North Pole weights(1) = pi2 weights(2) = pi*pi weights(3) = zero weights(4) = pi2 weights(5) = pi*pi weights(6) = zero grid1_add = 0 pole_loop1: do n=1,grid1_size if (grid1_area(n) < -three*pih .and. & grid1_center_lat(n) > zero) then grid1_add = n exit pole_loop1 endif end do pole_loop1 grid2_add = 0 pole_loop2: do n=1,grid2_size if (grid2_area(n) < -three*pih .and. & grid2_center_lat(n) > zero) then grid2_add = n exit pole_loop2 endif end do pole_loop2 if (grid1_add /=0) then grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1) grid1_centroid_lat(grid1_add) = & grid1_centroid_lat(grid1_add) + weights(2) grid1_centroid_lon(grid1_add) = & grid1_centroid_lon(grid1_add) + weights(3) endif if (grid2_add /=0) then grid2_area(grid2_add) = grid2_area(grid2_add) + & weights(num_wts+1) grid2_centroid_lat(grid2_add) = & grid2_centroid_lat(grid2_add) + weights(num_wts+2) grid2_centroid_lon(grid2_add) = & grid2_centroid_lon(grid2_add) + weights(num_wts+3) endif if (grid1_add /= 0 .and. grid2_add /=0) then call store_link_cnsrv(grid1_add, grid2_add, weights) grid1_frac(grid1_add) = grid1_frac(grid1_add) + & weights(1) grid2_frac(grid2_add) = grid2_frac(grid2_add) + & weights(num_wts+1) endif !*** South Pole weights(1) = pi2 weights(2) = -pi*pi weights(3) = zero weights(4) = pi2 weights(5) = -pi*pi weights(6) = zero grid1_add = 0 pole_loop3: do n=1,grid1_size if (grid1_area(n) < -three*pih .and. & grid1_center_lat(n) < zero) then grid1_add = n exit pole_loop3 endif end do pole_loop3 grid2_add = 0 pole_loop4: do n=1,grid2_size if (grid2_area(n) < -three*pih .and. & grid2_center_lat(n) < zero) then grid2_add = n exit pole_loop4 endif end do pole_loop4 if (grid1_add /=0) then grid1_area(grid1_add) = grid1_area(grid1_add) + weights(1) grid1_centroid_lat(grid1_add) = & grid1_centroid_lat(grid1_add) + weights(2) grid1_centroid_lon(grid1_add) = & grid1_centroid_lon(grid1_add) + weights(3) endif if (grid2_add /=0) then grid2_area(grid2_add) = grid2_area(grid2_add) + & weights(num_wts+1) grid2_centroid_lat(grid2_add) = & grid2_centroid_lat(grid2_add) + weights(num_wts+2) grid2_centroid_lon(grid2_add) = & grid2_centroid_lon(grid2_add) + weights(num_wts+3) endif if (grid1_add /= 0 .and. grid2_add /=0) then call store_link_cnsrv(grid1_add, grid2_add, weights) grid1_frac(grid1_add) = grid1_frac(grid1_add) + & weights(1) grid2_frac(grid2_add) = grid2_frac(grid2_add) + & weights(num_wts+1) endif !----------------------------------------------------------------------- ! ! finish centroid computation ! !----------------------------------------------------------------------- where (grid1_area /= zero) grid1_centroid_lat = grid1_centroid_lat/grid1_area grid1_centroid_lon = grid1_centroid_lon/grid1_area end where where (grid2_area /= zero) grid2_centroid_lat = grid2_centroid_lat/grid2_area grid2_centroid_lon = grid2_centroid_lon/grid2_area end where !----------------------------------------------------------------------- ! ! include centroids in weights and normalize using destination ! area if requested ! !----------------------------------------------------------------------- do n=1,num_links_map1 grid1_add = grid1_add_map1(n) grid2_add = grid2_add_map1(n) do nwgt=1,num_wts weights( nwgt) = wts_map1(nwgt,n) if (num_maps > 1) then weights(num_wts+nwgt) = wts_map2(nwgt,n) endif end do select case(norm_opt) case (norm_opt_dstarea) if (grid2_area(grid2_add) /= zero) then if (luse_grid2_area) then norm_factor = one/grid2_area_in(grid2_add) else norm_factor = one/grid2_area(grid2_add) endif else norm_factor = zero endif case (norm_opt_frcarea) if (grid2_frac(grid2_add) /= zero) then if (luse_grid2_area) then norm_factor = grid2_area(grid2_add)/ & (grid2_frac(grid2_add)* & grid2_area_in(grid2_add)) else norm_factor = one/grid2_frac(grid2_add) endif else norm_factor = zero endif case (norm_opt_none) norm_factor = one end select wts_map1(1,n) = weights(1)*norm_factor wts_map1(2,n) = (weights(2) - weights(1)* & grid1_centroid_lat(grid1_add))* & norm_factor wts_map1(3,n) = (weights(3) - weights(1)* & grid1_centroid_lon(grid1_add))* & norm_factor if (num_maps > 1) then select case(norm_opt) case (norm_opt_dstarea) if (grid1_area(grid1_add) /= zero) then if (luse_grid1_area) then norm_factor = one/grid1_area_in(grid1_add) else norm_factor = one/grid1_area(grid1_add) endif else norm_factor = zero endif case (norm_opt_frcarea) if (grid1_frac(grid1_add) /= zero) then if (luse_grid1_area) then norm_factor = grid1_area(grid1_add)/ & (grid1_frac(grid1_add)* & grid1_area_in(grid1_add)) else norm_factor = one/grid1_frac(grid1_add) endif else norm_factor = zero endif case (norm_opt_none) norm_factor = one end select wts_map2(1,n) = weights(num_wts+1)*norm_factor wts_map2(2,n) = (weights(num_wts+2) - weights(num_wts+1)* & grid2_centroid_lat(grid2_add))* & norm_factor wts_map2(3,n) = (weights(num_wts+3) - weights(num_wts+1)* & grid2_centroid_lon(grid2_add))* & norm_factor endif end do print *, 'Total number of links = ',num_links_map1 where (grid1_area /= zero) grid1_frac = grid1_frac/grid1_area where (grid2_area /= zero) grid2_frac = grid2_frac/grid2_area !----------------------------------------------------------------------- ! ! perform some error checking on final weights ! !----------------------------------------------------------------------- grid2_centroid_lat = zero grid2_centroid_lon = zero do n=1,grid1_size if (grid1_area(n) < -.01) then print *,'Grid 1 area error: ',n,grid1_area(n) endif if (grid1_centroid_lat(n) < -pih-.01 .or. & grid1_centroid_lat(n) > pih+.01) then print *,'Grid 1 centroid lat error: ',n,grid1_centroid_lat(n) endif grid1_centroid_lat(n) = zero grid1_centroid_lon(n) = zero end do do n=1,grid2_size if (grid2_area(n) < -.01) then print *,'Grid 2 area error: ',n,grid2_area(n) endif if (grid2_centroid_lat(n) < -pih-.01 .or. & grid2_centroid_lat(n) > pih+.01) then print *,'Grid 2 centroid lat error: ',n,grid2_centroid_lat(n) endif grid2_centroid_lat(n) = zero grid2_centroid_lon(n) = zero end do do n=1,num_links_map1 grid1_add = grid1_add_map1(n) grid2_add = grid2_add_map1(n) if (wts_map1(1,n) < -.01) then print *,'Map 1 weight < 0 ',grid1_add,grid2_add,wts_map1(1,n) endif if (norm_opt /= norm_opt_none .and. wts_map1(1,n) > 1.01) then print *,'Map 1 weight > 1 ',grid1_add,grid2_add,wts_map1(1,n) endif grid2_centroid_lat(grid2_add) = & grid2_centroid_lat(grid2_add) + wts_map1(1,n) if (num_maps > 1) then if (wts_map2(1,n) < -.01) then print *,'Map 2 weight < 0 ',grid1_add,grid2_add, & wts_map2(1,n) endif if (norm_opt /= norm_opt_none .and. wts_map2(1,n) > 1.01) then print *,'Map 2 weight < 0 ',grid1_add,grid2_add, & wts_map2(1,n) endif grid1_centroid_lat(grid1_add) = & grid1_centroid_lat(grid1_add) + wts_map2(1,n) endif end do do n=1,grid2_size select case(norm_opt) case (norm_opt_dstarea) norm_factor = grid2_frac(grid2_add) case (norm_opt_frcarea) norm_factor = one case (norm_opt_none) if (luse_grid2_area) then norm_factor = grid2_area_in(grid2_add) else norm_factor = grid2_area(grid2_add) endif end select if (abs(grid2_centroid_lat(grid2_add)-norm_factor) > .01) then print *,'Error: sum of wts for map1 ',grid2_add, & grid2_centroid_lat(grid2_add),norm_factor endif end do if (num_maps > 1) then do n=1,grid1_size select case(norm_opt) case (norm_opt_dstarea) norm_factor = grid1_frac(grid1_add) case (norm_opt_frcarea) norm_factor = one case (norm_opt_none) if (luse_grid1_area) then norm_factor = grid1_area_in(grid1_add) else norm_factor = grid1_area(grid1_add) endif end select if (abs(grid1_centroid_lat(grid1_add)-norm_factor) > .01) then print *,'Error: sum of wts for map2 ',grid1_add, & grid1_centroid_lat(grid1_add),norm_factor endif end do endif !----------------------------------------------------------------------- end subroutine remap_conserv !*********************************************************************** subroutine intersection(location,intrsct_lat,intrsct_lon,lcoinc, & beglat, beglon, endlat, endlon, begseg, & lbegin, lrevers) !----------------------------------------------------------------------- ! ! this routine finds the next intersection of a destination grid ! line with the line segment given by beglon, endlon, etc. ! a coincidence flag is returned if the segment is entirely ! coincident with an ocean grid line. the cells in which to search ! for an intersection must have already been restricted in the ! calling routine. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! intent(in): ! !----------------------------------------------------------------------- logical (kind=log_kind), intent(in) :: & lbegin, & ! flag for first integration along this segment lrevers ! flag whether segment integrated in reverse real (kind=dbl_kind), intent(in) :: & beglat, beglon, & ! beginning lat/lon endpoints for segment endlat, endlon ! ending lat/lon endpoints for segment real (kind=dbl_kind), dimension(2), intent(inout) :: & begseg ! begin lat/lon of full segment !----------------------------------------------------------------------- ! ! intent(out): ! !----------------------------------------------------------------------- integer (kind=int_kind), intent(out) :: & location ! address in destination array containing this ! segment logical (kind=log_kind), intent(out) :: & lcoinc ! flag segments which are entirely coincident ! with a grid line real (kind=dbl_kind), intent(out) :: & intrsct_lat, intrsct_lon ! lat/lon coords of next intersect. !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (kind=int_kind) :: n, next_n, cell, srch_corners, pole_loc integer (kind=int_kind), save :: & last_loc ! save location when crossing threshold logical (kind=log_kind) :: & loutside ! flags points outside grid logical (kind=log_kind), save :: & lthresh = .false. ! flags segments crossing threshold bndy real (kind=dbl_kind) :: & lon1, lon2, & ! local longitude variables for segment lat1, lat2, & ! local latitude variables for segment grdlon1, grdlon2, & ! local longitude variables for grid cell grdlat1, grdlat2, & ! local latitude variables for grid cell vec1_lat, vec1_lon, & ! vectors and cross products used vec2_lat, vec2_lon, & ! during grid search cross_product, & eps, offset, & ! small offset away from intersect s1, s2, determ, & ! variables used for linear solve to mat1, mat2, mat3, mat4, rhs1, rhs2 ! find intersection real (kind=dbl_kind), save :: & intrsct_lat_off, intrsct_lon_off ! lat/lon coords offset ! for next search !----------------------------------------------------------------------- ! ! initialize defaults, flags, etc. ! !----------------------------------------------------------------------- location = 0 lcoinc = .false. intrsct_lat = endlat intrsct_lon = endlon if (num_srch_cells == 0) return if (beglat > north_thresh .or. beglat < south_thresh) then if (lthresh) location = last_loc call pole_intersection(location, & intrsct_lat,intrsct_lon,lcoinc,lthresh, & beglat, beglon, endlat, endlon, begseg, lrevers) if (lthresh) then last_loc = location intrsct_lat_off = intrsct_lat intrsct_lon_off = intrsct_lon endif return endif loutside = .false. if (lbegin) then lat1 = beglat lon1 = beglon else lat1 = intrsct_lat_off lon1 = intrsct_lon_off endif lat2 = endlat lon2 = endlon if ((lon2-lon1) > three*pih) then lon2 = lon2 - pi2 else if ((lon2-lon1) < -three*pih) then lon2 = lon2 + pi2 endif s1 = zero !----------------------------------------------------------------------- ! ! search for location of this segment in ocean grid using cross ! product method to determine whether a point is enclosed by a cell ! !----------------------------------------------------------------------- call timer_start(12) srch_corners = size(srch_corner_lat,DIM=1) srch_loop: do !*** !*** if last segment crossed threshold, use that location !*** if (lthresh) then do cell=1,num_srch_cells if (srch_add(cell) == last_loc) then location = last_loc eps = tiny exit srch_loop endif end do endif !*** !*** otherwise normal search algorithm !*** cell_loop: do cell=1,num_srch_cells corner_loop: do n=1,srch_corners next_n = MOD(n,srch_corners) + 1 !*** !*** here we take the cross product of the vector making !*** up each cell side with the vector formed by the vertex !*** and search point. if all the cross products are !*** positive, the point is contained in the cell. !*** vec1_lat = srch_corner_lat(next_n,cell) - & srch_corner_lat(n ,cell) vec1_lon = srch_corner_lon(next_n,cell) - & srch_corner_lon(n ,cell) vec2_lat = lat1 - srch_corner_lat(n,cell) vec2_lon = lon1 - srch_corner_lon(n,cell) !*** !*** if endpoint coincident with vertex, offset !*** the endpoint !*** if (vec2_lat == 0 .and. vec2_lon == 0) then lat1 = lat1 + 1.d-10*(lat2-lat1) lon1 = lon1 + 1.d-10*(lon2-lon1) vec2_lat = lat1 - srch_corner_lat(n,cell) vec2_lon = lon1 - srch_corner_lon(n,cell) endif !*** !*** check for 0,2pi crossings !*** if (vec1_lon > pi) then vec1_lon = vec1_lon - pi2 else if (vec1_lon < -pi) then vec1_lon = vec1_lon + pi2 endif if (vec2_lon > pi) then vec2_lon = vec2_lon - pi2 else if (vec2_lon < -pi) then vec2_lon = vec2_lon + pi2 endif cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat !*** !*** if the cross product for a side is zero, the point !*** lies exactly on the side or the side is degenerate !*** (zero length). if degenerate, set the cross !*** product to a positive number. otherwise perform !*** another cross product between the side and the !*** segment itself. !*** if this cross product is also zero, the line is !*** coincident with the cell boundary - perform the !*** dot product and only choose the cell if the dot !*** product is positive (parallel vs anti-parallel). !*** if (cross_product == zero) then if (vec1_lat /= zero .or. vec1_lon /= zero) then vec2_lat = lat2 - lat1 vec2_lon = lon2 - lon1 if (vec2_lon > pi) then vec2_lon = vec2_lon - pi2 else if (vec2_lon < -pi) then vec2_lon = vec2_lon + pi2 endif cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat else cross_product = one endif if (cross_product == zero) then lcoinc = .true. cross_product = vec1_lon*vec2_lon + vec1_lat*vec2_lat if (lrevers) cross_product = -cross_product endif endif !*** !*** if cross product is less than zero, this cell !*** doesn't work !*** if (cross_product < zero) exit corner_loop end do corner_loop !*** !*** if cross products all positive, we found the location !*** if (n > srch_corners) then location = srch_add(cell) !*** !*** if the beginning of this segment was outside the !*** grid, invert the segment so the intersection found !*** will be the first intersection with the grid !*** if (loutside) then lat2 = beglat lon2 = beglon location = 0 eps = -tiny else eps = tiny endif exit srch_loop endif !*** !*** otherwise move on to next cell !*** end do cell_loop !*** !*** if still no cell found, the point lies outside the grid. !*** take some baby steps along the segment to see if any !*** part of the segment lies inside the grid. !*** loutside = .true. s1 = s1 + 0.001_dbl_kind lat1 = beglat + s1*(endlat - beglat) lon1 = beglon + s1*(lon2 - beglon) !*** !*** reached the end of the segment and still outside the grid !*** return no intersection !*** if (s1 >= one) return end do srch_loop call timer_stop(12) !----------------------------------------------------------------------- ! ! now that a cell is found, search for the next intersection. ! loop over sides of the cell to find intersection with side ! must check all sides for coincidences or intersections ! !----------------------------------------------------------------------- call timer_start(13) intrsct_loop: do n=1,srch_corners next_n = mod(n,srch_corners) + 1 grdlon1 = srch_corner_lon(n ,cell) grdlon2 = srch_corner_lon(next_n,cell) grdlat1 = srch_corner_lat(n ,cell) grdlat2 = srch_corner_lat(next_n,cell) !*** !*** set up linear system to solve for intersection !*** mat1 = lat2 - lat1 mat2 = grdlat1 - grdlat2 mat3 = lon2 - lon1 mat4 = grdlon1 - grdlon2 rhs1 = grdlat1 - lat1 rhs2 = grdlon1 - lon1 if (mat3 > pi) then mat3 = mat3 - pi2 else if (mat3 < -pi) then mat3 = mat3 + pi2 endif if (mat4 > pi) then mat4 = mat4 - pi2 else if (mat4 < -pi) then mat4 = mat4 + pi2 endif if (rhs2 > pi) then rhs2 = rhs2 - pi2 else if (rhs2 < -pi) then rhs2 = rhs2 + pi2 endif determ = mat1*mat4 - mat2*mat3 !*** !*** if the determinant is zero, the segments are either !*** parallel or coincident. coincidences were detected !*** above so do nothing. !*** if the determinant is non-zero, solve for the linear !*** parameters s for the intersection point on each line !*** segment. !*** if 0 1.e-30) then s1 = (rhs1*mat4 - mat2*rhs2)/determ s2 = (mat1*rhs2 - rhs1*mat3)/determ if (s2 >= zero .and. s2 <= one .and. & s1 > zero .and. s1 <= one) then !*** !*** recompute intersection based on full segment !*** so intersections are consistent for both sweeps !*** if (.not. loutside) then mat1 = lat2 - begseg(1) mat3 = lon2 - begseg(2) rhs1 = grdlat1 - begseg(1) rhs2 = grdlon1 - begseg(2) else mat1 = begseg(1) - endlat mat3 = begseg(2) - endlon rhs1 = grdlat1 - endlat rhs2 = grdlon1 - endlon endif if (mat3 > pi) then mat3 = mat3 - pi2 else if (mat3 < -pi) then mat3 = mat3 + pi2 endif if (rhs2 > pi) then rhs2 = rhs2 - pi2 else if (rhs2 < -pi) then rhs2 = rhs2 + pi2 endif determ = mat1*mat4 - mat2*mat3 !*** !*** sometimes due to roundoff, the previous !*** determinant is non-zero, but the lines !*** are actually coincident. if this is the !*** case, skip the rest. !*** if (determ /= zero) then s1 = (rhs1*mat4 - mat2*rhs2)/determ s2 = (mat1*rhs2 - rhs1*mat3)/determ offset = s1 + eps/determ if (offset > one) offset = one if (.not. loutside) then intrsct_lat = begseg(1) + mat1*s1 intrsct_lon = begseg(2) + mat3*s1 intrsct_lat_off = begseg(1) + mat1*offset intrsct_lon_off = begseg(2) + mat3*offset else intrsct_lat = endlat + mat1*s1 intrsct_lon = endlon + mat3*s1 intrsct_lat_off = endlat + mat1*offset intrsct_lon_off = endlon + mat3*offset endif exit intrsct_loop endif endif endif !*** !*** no intersection this side, move on to next side !*** end do intrsct_loop call timer_stop(13) !----------------------------------------------------------------------- ! ! if the segment crosses a pole threshold, reset the intersection ! to be the threshold latitude. only check if this was not a ! threshold segment since sometimes coordinate transform can end ! up on other side of threshold again. ! !----------------------------------------------------------------------- if (lthresh) then if (intrsct_lat < north_thresh .or. intrsct_lat > south_thresh) & lthresh = .false. else if (lat1 > zero .and. intrsct_lat > north_thresh) then intrsct_lat = north_thresh + tiny intrsct_lat_off = north_thresh + eps*mat1 s1 = (intrsct_lat - begseg(1))/mat1 intrsct_lon = begseg(2) + s1*mat3 intrsct_lon_off = begseg(2) + (s1+eps)*mat3 last_loc = location lthresh = .true. else if (lat1 < zero .and. intrsct_lat < south_thresh) then intrsct_lat = south_thresh - tiny intrsct_lat_off = south_thresh + eps*mat1 s1 = (intrsct_lat - begseg(1))/mat1 intrsct_lon = begseg(2) + s1*mat3 intrsct_lon_off = begseg(2) + (s1+eps)*mat3 last_loc = location lthresh = .true. endif !----------------------------------------------------------------------- end subroutine intersection !*********************************************************************** subroutine pole_intersection(location, & intrsct_lat,intrsct_lon,lcoinc,lthresh, & beglat, beglon, endlat, endlon, begseg, lrevers) !----------------------------------------------------------------------- ! ! this routine is identical to the intersection routine except ! that a coordinate transformation (using a Lambert azimuthal ! equivalent projection) is performed to treat polar cells more ! accurately. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! intent(in): ! !----------------------------------------------------------------------- real (kind=dbl_kind), intent(in) :: & beglat, beglon, & ! beginning lat/lon endpoints for segment endlat, endlon ! ending lat/lon endpoints for segment real (kind=dbl_kind), dimension(2), intent(inout) :: & begseg ! begin lat/lon of full segment logical (kind=log_kind), intent(in) :: & lrevers ! flag true if segment integrated in reverse !----------------------------------------------------------------------- ! ! intent(out): ! !----------------------------------------------------------------------- integer (kind=int_kind), intent(inout) :: & location ! address in destination array containing this ! segment -- also may contain last location on ! entry logical (kind=log_kind), intent(out) :: & lcoinc ! flag segment coincident with grid line logical (kind=log_kind), intent(inout) :: & lthresh ! flag segment crossing threshold boundary real (kind=dbl_kind), intent(out) :: & intrsct_lat, intrsct_lon ! lat/lon coords of next intersect. !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (kind=int_kind) :: n, next_n, cell, srch_corners, pole_loc logical (kind=log_kind) :: loutside ! flags points outside grid real (kind=dbl_kind) :: pi4, rns, & ! north/south conversion x1, x2, & ! local x variables for segment y1, y2, & ! local y variables for segment begx, begy, & ! beginning x,y variables for segment endx, endy, & ! beginning x,y variables for segment begsegx, begsegy, & ! beginning x,y variables for segment grdx1, grdx2, & ! local x variables for grid cell grdy1, grdy2, & ! local y variables for grid cell vec1_y, vec1_x, & ! vectors and cross products used vec2_y, vec2_x, & ! during grid search cross_product, eps, & ! eps=small offset away from intersect s1, s2, determ, & ! variables used for linear solve to mat1, mat2, mat3, mat4, rhs1, rhs2 ! find intersection real (kind=dbl_kind), dimension(:,:), allocatable :: & srch_corner_x, & ! x of each corner of srch cells srch_corner_y ! y of each corner of srch cells !*** !*** save last intersection to avoid roundoff during coord !*** transformation !*** logical (kind=log_kind), save :: luse_last = .false. real (kind=dbl_kind), save :: & intrsct_x, intrsct_y ! x,y for intersection !*** !*** variables necessary if segment manages to hit pole !*** integer (kind=int_kind), save :: & avoid_pole_count = 0 ! count attempts to avoid pole real (kind=dbl_kind), save :: & avoid_pole_offset = tiny ! endpoint offset to avoid pole !----------------------------------------------------------------------- ! ! initialize defaults, flags, etc. ! !----------------------------------------------------------------------- if (.not. lthresh) location = 0 lcoinc = .false. intrsct_lat = endlat intrsct_lon = endlon loutside = .false. s1 = zero !----------------------------------------------------------------------- ! ! convert coordinates ! !----------------------------------------------------------------------- allocate(srch_corner_x(size(srch_corner_lat,DIM=1), & size(srch_corner_lat,DIM=2)), & srch_corner_y(size(srch_corner_lat,DIM=1), & size(srch_corner_lat,DIM=2))) if (beglat > zero) then pi4 = quart*pi rns = one else pi4 = -quart*pi rns = -one endif if (luse_last) then x1 = intrsct_x y1 = intrsct_y else x1 = rns*two*sin(pi4 - half*beglat)*cos(beglon) y1 = two*sin(pi4 - half*beglat)*sin(beglon) luse_last = .true. endif x2 = rns*two*sin(pi4 - half*endlat)*cos(endlon) y2 = two*sin(pi4 - half*endlat)*sin(endlon) srch_corner_x = rns*two*sin(pi4 - half*srch_corner_lat)* & cos(srch_corner_lon) srch_corner_y = two*sin(pi4 - half*srch_corner_lat)* & sin(srch_corner_lon) begx = x1 begy = y1 endx = x2 endy = y2 begsegx = rns*two*sin(pi4 - half*begseg(1))*cos(begseg(2)) begsegy = two*sin(pi4 - half*begseg(1))*sin(begseg(2)) intrsct_x = endx intrsct_y = endy !----------------------------------------------------------------------- ! ! search for location of this segment in ocean grid using cross ! product method to determine whether a point is enclosed by a cell ! !----------------------------------------------------------------------- call timer_start(12) srch_corners = size(srch_corner_lat,DIM=1) srch_loop: do !*** !*** if last segment crossed threshold, use that location !*** if (lthresh) then do cell=1,num_srch_cells if (srch_add(cell) == location) then eps = tiny exit srch_loop endif end do endif !*** !*** otherwise normal search algorithm !*** cell_loop: do cell=1,num_srch_cells corner_loop: do n=1,srch_corners next_n = MOD(n,srch_corners) + 1 !*** !*** here we take the cross product of the vector making !*** up each cell side with the vector formed by the vertex !*** and search point. if all the cross products are !*** positive, the point is contained in the cell. !*** vec1_x = srch_corner_x(next_n,cell) - & srch_corner_x(n ,cell) vec1_y = srch_corner_y(next_n,cell) - & srch_corner_y(n ,cell) vec2_x = x1 - srch_corner_x(n,cell) vec2_y = y1 - srch_corner_y(n,cell) !*** !*** if endpoint coincident with vertex, offset !*** the endpoint !*** if (vec2_x == 0 .and. vec2_y == 0) then x1 = x1 + 1.d-10*(x2-x1) y1 = y1 + 1.d-10*(y2-y1) vec2_x = x1 - srch_corner_x(n,cell) vec2_y = y1 - srch_corner_y(n,cell) endif cross_product = vec1_x*vec2_y - vec2_x*vec1_y !*** !*** if the cross product for a side is zero, the point !*** lies exactly on the side or the length of a side !*** is zero. if the length is zero set det > 0. !*** otherwise, perform another cross !*** product between the side and the segment itself. !*** if this cross product is also zero, the line is !*** coincident with the cell boundary - perform the !*** dot product and only choose the cell if the dot !*** product is positive (parallel vs anti-parallel). !*** if (cross_product == zero) then if (vec1_x /= zero .or. vec1_y /= 0) then vec2_x = x2 - x1 vec2_y = y2 - y1 cross_product = vec1_x*vec2_y - vec2_x*vec1_y else cross_product = one endif if (cross_product == zero) then lcoinc = .true. cross_product = vec1_x*vec2_x + vec1_y*vec2_y if (lrevers) cross_product = -cross_product endif endif !*** !*** if cross product is less than zero, this cell !*** doesn't work !*** if (cross_product < zero) exit corner_loop end do corner_loop !*** !*** if cross products all positive, we found the location !*** if (n > srch_corners) then location = srch_add(cell) !*** !*** if the beginning of this segment was outside the !*** grid, invert the segment so the intersection found !*** will be the first intersection with the grid !*** if (loutside) then x2 = begx y2 = begy location = 0 eps = -tiny else eps = tiny endif exit srch_loop endif !*** !*** otherwise move on to next cell !*** end do cell_loop !*** !*** if no cell found, the point lies outside the grid. !*** take some baby steps along the segment to see if any !*** part of the segment lies inside the grid. !*** loutside = .true. s1 = s1 + 0.001_dbl_kind x1 = begx + s1*(x2 - begx) y1 = begy + s1*(y2 - begy) !*** !*** reached the end of the segment and still outside the grid !*** return no intersection !*** if (s1 >= one) then deallocate(srch_corner_x, srch_corner_y) luse_last = .false. return endif end do srch_loop call timer_stop(12) !----------------------------------------------------------------------- ! ! now that a cell is found, search for the next intersection. ! loop over sides of the cell to find intersection with side ! must check all sides for coincidences or intersections ! !----------------------------------------------------------------------- call timer_start(13) intrsct_loop: do n=1,srch_corners next_n = mod(n,srch_corners) + 1 grdy1 = srch_corner_y(n ,cell) grdy2 = srch_corner_y(next_n,cell) grdx1 = srch_corner_x(n ,cell) grdx2 = srch_corner_x(next_n,cell) !*** !*** set up linear system to solve for intersection !*** mat1 = x2 - x1 mat2 = grdx1 - grdx2 mat3 = y2 - y1 mat4 = grdy1 - grdy2 rhs1 = grdx1 - x1 rhs2 = grdy1 - y1 determ = mat1*mat4 - mat2*mat3 !*** !*** if the determinant is zero, the segments are either !*** parallel or coincident or one segment has zero length. !*** coincidences were detected above so do nothing. !*** if the determinant is non-zero, solve for the linear !*** parameters s for the intersection point on each line !*** segment. !*** if 0 1.e-30) then s1 = (rhs1*mat4 - mat2*rhs2)/determ s2 = (mat1*rhs2 - rhs1*mat3)/determ if (s2 >= zero .and. s2 <= one .and. & s1 > zero .and. s1 <= one) then !*** !*** recompute intersection using entire segment !*** for consistency between sweeps !*** if (.not. loutside) then mat1 = x2 - begsegx mat3 = y2 - begsegy rhs1 = grdx1 - begsegx rhs2 = grdy1 - begsegy else mat1 = x2 - endx mat3 = y2 - endy rhs1 = grdx1 - endx rhs2 = grdy1 - endy endif determ = mat1*mat4 - mat2*mat3 !*** !*** sometimes due to roundoff, the previous !*** determinant is non-zero, but the lines !*** are actually coincident. if this is the !*** case, skip the rest. !*** if (determ /= zero) then s1 = (rhs1*mat4 - mat2*rhs2)/determ s2 = (mat1*rhs2 - rhs1*mat3)/determ if (.not. loutside) then intrsct_x = begsegx + s1*mat1 intrsct_y = begsegy + s1*mat3 else intrsct_x = endx + s1*mat1 intrsct_y = endy + s1*mat3 endif !*** !*** convert back to lat/lon coordinates !*** intrsct_lon = rns*atan2(intrsct_y,intrsct_x) if (intrsct_lon < zero) & intrsct_lon = intrsct_lon + pi2 if (abs(intrsct_x) > 1.d-10) then intrsct_lat = (pi4 - & asin(rns*half*intrsct_x/cos(intrsct_lon)))*two else if (abs(intrsct_y) > 1.d-10) then intrsct_lat = (pi4 - & asin(half*intrsct_y/sin(intrsct_lon)))*two else intrsct_lat = two*pi4 endif !*** !*** add offset in transformed space for next pass. !*** if (s1 - eps/determ < one) then intrsct_x = intrsct_x - mat1*(eps/determ) intrsct_y = intrsct_y - mat3*(eps/determ) else if (.not. loutside) then intrsct_x = endx intrsct_y = endy intrsct_lat = endlat intrsct_lon = endlon else intrsct_x = begsegx intrsct_y = begsegy intrsct_lat = begseg(1) intrsct_lon = begseg(2) endif endif exit intrsct_loop endif endif endif !*** !*** no intersection this side, move on to next side !*** end do intrsct_loop call timer_stop(13) deallocate(srch_corner_x, srch_corner_y) !----------------------------------------------------------------------- ! ! if segment manages to cross over pole, shift the beginning ! endpoint in order to avoid hitting pole directly ! (it is ok for endpoint to be pole point) ! !----------------------------------------------------------------------- if (abs(intrsct_x) < 1.e-10 .and. abs(intrsct_y) < 1.e-10 .and. & (endx /= zero .and. endy /=0)) then if (avoid_pole_count > 2) then avoid_pole_count = 0 avoid_pole_offset = 10.*avoid_pole_offset endif cross_product = begsegx*(endy-begsegy) - begsegy*(endx-begsegx) intrsct_lat = begseg(1) if (cross_product*intrsct_lat > zero) then intrsct_lon = beglon + avoid_pole_offset begseg(2) = begseg(2) + avoid_pole_offset else intrsct_lon = beglon - avoid_pole_offset begseg(2) = begseg(2) - avoid_pole_offset endif avoid_pole_count = avoid_pole_count + 1 luse_last = .false. else avoid_pole_count = 0 avoid_pole_offset = tiny endif !----------------------------------------------------------------------- ! ! if the segment crosses a pole threshold, reset the intersection ! to be the threshold latitude and do not reuse x,y intersect ! on next entry. only check if did not cross threshold last ! time - sometimes the coordinate transformation can place a ! segment on the other side of the threshold again ! !----------------------------------------------------------------------- if (lthresh) then if (intrsct_lat > north_thresh .or. intrsct_lat < south_thresh) & lthresh = .false. else if (beglat > zero .and. intrsct_lat < north_thresh) then mat4 = endlat - begseg(1) mat3 = endlon - begseg(2) if (mat3 > pi) mat3 = mat3 - pi2 if (mat3 < -pi) mat3 = mat3 + pi2 intrsct_lat = north_thresh - tiny s1 = (north_thresh - begseg(1))/mat4 intrsct_lon = begseg(2) + s1*mat3 luse_last = .false. lthresh = .true. else if (beglat < zero .and. intrsct_lat > south_thresh) then mat4 = endlat - begseg(1) mat3 = endlon - begseg(2) if (mat3 > pi) mat3 = mat3 - pi2 if (mat3 < -pi) mat3 = mat3 + pi2 intrsct_lat = south_thresh + tiny s1 = (south_thresh - begseg(1))/mat4 intrsct_lon = begseg(2) + s1*mat3 luse_last = .false. lthresh = .true. endif !*** !*** if reached end of segment, do not use x,y intersect !*** on next entry !*** if (intrsct_lat == endlat .and. intrsct_lon == endlon) then luse_last = .false. endif !----------------------------------------------------------------------- end subroutine pole_intersection !*********************************************************************** subroutine line_integral(weights, num_wts, & in_phi1, in_phi2, theta1, theta2, & grid1_lat, grid1_lon, grid2_lat, grid2_lon) !----------------------------------------------------------------------- ! ! this routine computes the line integral of the flux function ! that results in the interpolation weights. the line is defined ! by the input lat/lon of the endpoints. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! intent(in): ! !----------------------------------------------------------------------- integer (kind=int_kind), intent(in) :: & num_wts ! number of weights to compute real (kind=dbl_kind), intent(in) :: & in_phi1, in_phi2, & ! longitude endpoints for the segment theta1, theta2, & ! latitude endpoints for the segment grid1_lat, grid1_lon, & ! reference coordinates for each grid2_lat, grid2_lon ! grid (to ensure correct 0,2pi interv. !----------------------------------------------------------------------- ! ! intent(out): ! !----------------------------------------------------------------------- real (kind=dbl_kind), dimension(2*num_wts), intent(out) :: & weights ! line integral contribution to weights !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- real (kind=dbl_kind) :: dphi, sinth1, sinth2, costh1, costh2, fac, & phi1, phi2, phidiff1, phidiff2, sinint real (kind=dbl_kind) :: f1, f2, fint !----------------------------------------------------------------------- ! ! weights for the general case based on a trapezoidal approx to ! the integrals. ! !----------------------------------------------------------------------- sinth1 = SIN(theta1) sinth2 = SIN(theta2) costh1 = COS(theta1) costh2 = COS(theta2) dphi = in_phi1 - in_phi2 if (dphi > pi) then dphi = dphi - pi2 else if (dphi < -pi) then dphi = dphi + pi2 endif dphi = half*dphi !----------------------------------------------------------------------- ! ! the first weight is the area overlap integral. the second and ! fourth are second-order latitude gradient weights. ! !----------------------------------------------------------------------- weights( 1) = dphi*(sinth1 + sinth2) weights(num_wts+1) = dphi*(sinth1 + sinth2) weights( 2) = dphi*(costh1 + costh2 + (theta1*sinth1 + & theta2*sinth2)) weights(num_wts+2) = dphi*(costh1 + costh2 + (theta1*sinth1 + & theta2*sinth2)) !----------------------------------------------------------------------- ! ! the third and fifth weights are for the second-order phi gradient ! component. must be careful of longitude range. ! !----------------------------------------------------------------------- f1 = half*(costh1*sinth1 + theta1) f2 = half*(costh2*sinth2 + theta2) phi1 = in_phi1 - grid1_lon if (phi1 > pi) then phi1 = phi1 - pi2 else if (phi1 < -pi) then phi1 = phi1 + pi2 endif phi2 = in_phi2 - grid1_lon if (phi2 > pi) then phi2 = phi2 - pi2 else if (phi2 < -pi) then phi2 = phi2 + pi2 endif if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then weights(3) = dphi*(phi1*f1 + phi2*f2) else if (phi1 > zero) then fac = pi else fac = -pi endif fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi) weights(3) = half*phi1*(phi1-fac)*f1 - & half*phi2*(phi2+fac)*f2 + & half*fac*(phi1+phi2)*fint endif phi1 = in_phi1 - grid2_lon if (phi1 > pi) then phi1 = phi1 - pi2 else if (phi1 < -pi) then phi1 = phi1 + pi2 endif phi2 = in_phi2 - grid2_lon if (phi2 > pi) then phi2 = phi2 - pi2 else if (phi2 < -pi) then phi2 = phi2 + pi2 endif if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then weights(num_wts+3) = dphi*(phi1*f1 + phi2*f2) else if (phi1 > zero) then fac = pi else fac = -pi endif fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi) weights(num_wts+3) = half*phi1*(phi1-fac)*f1 - & half*phi2*(phi2+fac)*f2 + & half*fac*(phi1+phi2)*fint endif !----------------------------------------------------------------------- end subroutine line_integral !*********************************************************************** subroutine store_link_cnsrv(add1, add2, weights) !----------------------------------------------------------------------- ! ! this routine stores the address and weight for this link in ! the appropriate address and weight arrays and resizes those ! arrays if necessary. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! input variables ! !----------------------------------------------------------------------- integer (kind=int_kind), intent(in) :: & add1, & ! address on grid1 add2 ! address on grid2 real (kind=dbl_kind), dimension(:), intent(in) :: & weights ! array of remapping weights for this link !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (kind=int_kind) :: nlink, min_link, max_link ! link index integer (kind=int_kind), dimension(:,:), allocatable, save :: & link_add1, & ! min,max link add to restrict search link_add2 ! min,max link add to restrict search logical (kind=log_kind), save :: first_call = .true. !----------------------------------------------------------------------- ! ! if all weights are zero, do not bother storing the link ! !----------------------------------------------------------------------- if (all(weights == zero)) return !----------------------------------------------------------------------- ! ! restrict the range of links to search for existing links ! !----------------------------------------------------------------------- if (first_call) then allocate(link_add1(2,grid1_size), link_add2(2,grid2_size)) link_add1 = 0 link_add2 = 0 first_call = .false. min_link = 1 max_link = 0 else min_link = min(link_add1(1,add1),link_add2(1,add2)) max_link = max(link_add1(2,add1),link_add2(2,add2)) if (min_link == 0) then min_link = 1 max_link = 0 endif endif !----------------------------------------------------------------------- ! ! if the link already exists, add the weight to the current weight ! arrays ! !----------------------------------------------------------------------- do nlink=min_link,max_link if (add1 == grid1_add_map1(nlink)) then if (add2 == grid2_add_map1(nlink)) then wts_map1(:,nlink) = wts_map1(:,nlink) + weights(1:num_wts) if (num_maps == 2) then wts_map2(:,nlink) = wts_map2(:,nlink) + & weights(num_wts+1:2*num_wts) endif return endif endif end do !----------------------------------------------------------------------- ! ! if the link does not yet exist, increment number of links and ! check to see if remap arrays need to be increased to accomodate ! the new link. then store the link. ! !----------------------------------------------------------------------- num_links_map1 = num_links_map1 + 1 if (num_links_map1 > max_links_map1) & call resize_remap_vars(1,resize_increment) grid1_add_map1(num_links_map1) = add1 grid2_add_map1(num_links_map1) = add2 wts_map1 (:,num_links_map1) = weights(1:num_wts) if (num_maps > 1) then num_links_map2 = num_links_map2 + 1 if (num_links_map2 > max_links_map2) & call resize_remap_vars(2,resize_increment) grid1_add_map2(num_links_map2) = add1 grid2_add_map2(num_links_map2) = add2 wts_map2 (:,num_links_map2) = weights(num_wts+1:2*num_wts) endif if (link_add1(1,add1) == 0) link_add1(1,add1) = num_links_map1 if (link_add2(1,add2) == 0) link_add2(1,add2) = num_links_map1 link_add1(2,add1) = num_links_map1 link_add2(2,add2) = num_links_map1 !----------------------------------------------------------------------- end subroutine store_link_cnsrv !*********************************************************************** end module remap_conservative !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!