!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! this module contains necessary routines for performing an ! bilinear interpolation. ! !----------------------------------------------------------------------- ! ! 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_bilinear !----------------------------------------------------------------------- use kinds_mod ! defines common data types use constants ! defines common constants use grids ! module containing grid info use remap_vars ! module containing remap info implicit none !----------------------------------------------------------------------- integer (kind=int_kind), parameter :: & max_iter = 100 ! max iteration count for i,j iteration real (kind=dbl_kind), parameter :: & converge = 1.e-10_dbl_kind ! convergence criterion !*********************************************************************** contains !*********************************************************************** subroutine remap_bilin !----------------------------------------------------------------------- ! ! this routine computes the weights for a bilinear interpolation. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (kind=int_kind) :: n,icount, & dst_add, & ! destination address iter, & ! iteration counter nmap ! index of current map being computed integer (kind=int_kind), dimension(4) :: & src_add ! address for the four source points real (kind=dbl_kind), dimension(4) :: & src_lats, & ! latitudes of four bilinear corners src_lons, & ! longitudes of four bilinear corners wgts ! bilinear weights for four corners real (kind=dbl_kind) :: & plat, plon, & ! lat/lon coords of destination point iguess, jguess, & ! current guess for bilinear coordinate thguess, phguess, & ! current guess for lat/lon coordinate deli, delj, & ! corrections to i,j dth1, dth2, dth3, & ! some latitude differences dph1, dph2, dph3, & ! some longitude differences dthp, dphp, & ! difference between point and sw corner mat1, mat2, mat3, mat4, & ! matrix elements determinant, & ! matrix determinant sum_wgts ! sum of weights for normalization !----------------------------------------------------------------------- ! ! compute mappings from grid1 to grid2 ! !----------------------------------------------------------------------- nmap = 1 if (grid1_rank /= 2) then stop 'Can not do bilinear interpolation when grid_rank /= 2' endif !*** !*** loop over destination grid !*** grid_loop1: do dst_add = 1, grid2_size if (.not. grid2_mask(dst_add)) cycle grid_loop1 plat = grid2_center_lat(dst_add) plon = grid2_center_lon(dst_add) !*** !*** find nearest square of grid points on source grid !*** call grid_search_bilin(src_add, src_lats, src_lons, & plat, plon, grid1_dims, & grid1_center_lat, grid1_center_lon, & grid1_bound_box, bin_addr1, bin_addr2) !*** !*** check to see if points are land points !*** if (src_add(1) > 0) then do n=1,4 if (.not. grid1_mask(src_add(n))) src_add(1) = 0 end do endif !*** !*** if point found, find local i,j coordinates for weights !*** if (src_add(1) > 0) then grid2_frac(dst_add) = one !*** !*** iterate to find i,j for bilinear approximation !*** dth1 = src_lats(2) - src_lats(1) dth2 = src_lats(4) - src_lats(1) dth3 = src_lats(3) - src_lats(2) - dth2 dph1 = src_lons(2) - src_lons(1) dph2 = src_lons(4) - src_lons(1) dph3 = src_lons(3) - src_lons(2) if (dph1 > three*pih) dph1 = dph1 - pi2 if (dph2 > three*pih) dph2 = dph2 - pi2 if (dph3 > three*pih) dph3 = dph3 - pi2 if (dph1 < -three*pih) dph1 = dph1 + pi2 if (dph2 < -three*pih) dph2 = dph2 + pi2 if (dph3 < -three*pih) dph3 = dph3 + pi2 dph3 = dph3 - dph2 iguess = half jguess = half iter_loop1: do iter=1,max_iter dthp = plat - src_lats(1) - dth1*iguess - & dth2*jguess - dth3*iguess*jguess dphp = plon - src_lons(1) if (dphp > three*pih) dphp = dphp - pi2 if (dphp < -three*pih) dphp = dphp + pi2 dphp = dphp - dph1*iguess - dph2*jguess - & dph3*iguess*jguess mat1 = dth1 + dth3*jguess mat2 = dth2 + dth3*iguess mat3 = dph1 + dph3*jguess mat4 = dph2 + dph3*iguess determinant = mat1*mat4 - mat2*mat3 deli = (dthp*mat4 - mat2*dphp)/determinant delj = (mat1*dphp - dthp*mat3)/determinant if (abs(deli) < converge .and. & abs(delj) < converge) exit iter_loop1 iguess = iguess + deli jguess = jguess + delj end do iter_loop1 if (iter <= max_iter) then !*** !*** successfully found i,j - compute weights !*** wgts(1) = (one-iguess)*(one-jguess) wgts(2) = iguess*(one-jguess) wgts(3) = iguess*jguess wgts(4) = (one-iguess)*jguess call store_link_bilin(dst_add, src_add, wgts, nmap) else print *,'Point coords: ',plat,plon print *,'Dest grid lats: ',src_lats print *,'Dest grid lons: ',src_lons print *,'Dest grid addresses: ',src_add print *,'Current i,j : ',iguess, jguess stop 'Iteration for i,j exceed max iteration count' endif !*** !*** search for bilinear failed - use a distance-weighted !*** average instead (this is typically near the pole) !*** else if (src_add(1) < 0) then src_add = abs(src_add) icount = 0 do n=1,4 if (grid1_mask(src_add(n))) then icount = icount + 1 else src_lats(n) = zero endif end do if (icount > 0) then !*** renormalize weights sum_wgts = sum(src_lats) wgts(1) = src_lats(1)/sum_wgts wgts(2) = src_lats(2)/sum_wgts wgts(3) = src_lats(3)/sum_wgts wgts(4) = src_lats(4)/sum_wgts grid2_frac(dst_add) = one call store_link_bilin(dst_add, src_add, wgts, nmap) endif endif end do grid_loop1 !----------------------------------------------------------------------- ! ! compute mappings from grid2 to grid1 if necessary ! !----------------------------------------------------------------------- if (num_maps > 1) then nmap = 2 if (grid2_rank /= 2) then stop 'Can not do bilinear interpolation when grid_rank /= 2' endif !*** !*** loop over destination grid !*** grid_loop2: do dst_add = 1, grid1_size if (.not. grid1_mask(dst_add)) cycle grid_loop2 plat = grid1_center_lat(dst_add) plon = grid1_center_lon(dst_add) !*** !*** find nearest square of grid points on source grid !*** call grid_search_bilin(src_add, src_lats, src_lons, & plat, plon, grid2_dims, & grid2_center_lat, grid2_center_lon, & grid2_bound_box, bin_addr2, bin_addr1) !*** !*** check to see if points are land points !*** if (src_add(1) > 0) then do n=1,4 if (.not. grid2_mask(src_add(n))) src_add(1) = 0 end do endif !*** !*** if point found, find i,j coordinates for weights !*** if (src_add(1) > 0) then grid1_frac(dst_add) = one !*** !*** iterate to find i,j for bilinear approximation !*** dth1 = src_lats(2) - src_lats(1) dth2 = src_lats(4) - src_lats(1) dth3 = src_lats(3) - src_lats(2) - dth2 dph1 = src_lons(2) - src_lons(1) dph2 = src_lons(4) - src_lons(1) dph3 = src_lons(3) - src_lons(2) if (dph1 > pi) dph1 = dph1 - pi2 if (dph2 > pi) dph2 = dph2 - pi2 if (dph3 > pi) dph3 = dph3 - pi2 if (dph1 < -pi) dph1 = dph1 + pi2 if (dph2 < -pi) dph2 = dph2 + pi2 if (dph3 < -pi) dph3 = dph3 + pi2 dph3 = dph3 - dph2 iguess = zero jguess = zero iter_loop2: do iter=1,max_iter dthp = plat - src_lats(1) - dth1*iguess - & dth2*jguess - dth3*iguess*jguess dphp = plon - src_lons(1) if (dphp > pi) dphp = dphp - pi2 if (dphp < -pi) dphp = dphp + pi2 dphp = dphp - dph1*iguess - dph2*jguess - & dph3*iguess*jguess mat1 = dth1 + dth3*jguess mat2 = dth2 + dth3*iguess mat3 = dph1 + dph3*jguess mat4 = dph2 + dph3*iguess determinant = mat1*mat4 - mat2*mat3 deli = (dthp*mat4 - mat2*dphp)/determinant delj = (mat1*dphp - dthp*mat3)/determinant if (abs(deli) < converge .and. & abs(delj) < converge) exit iter_loop2 iguess = iguess + deli jguess = jguess + delj end do iter_loop2 if (iter <= max_iter) then !*** !*** successfully found i,j - compute weights !*** wgts(1) = (one-iguess)*(one-jguess) wgts(2) = iguess*(one-jguess) wgts(3) = iguess*jguess wgts(4) = (one-iguess)*jguess call store_link_bilin(dst_add, src_add, wgts, nmap) else print *,'Point coords: ',plat,plon print *,'Dest grid lats: ',src_lats print *,'Dest grid lons: ',src_lons print *,'Dest grid addresses: ',src_add print *,'Current i,j : ',iguess, jguess stop 'Iteration for i,j exceed max iteration count' endif !*** !*** search for bilinear failed - us a distance-weighted !*** average instead !*** else if (src_add(1) < 0) then src_add = abs(src_add) icount = 0 do n=1,4 if (grid2_mask(src_add(n))) then icount = icount + 1 else src_lats(n) = zero endif end do if (icount > 0) then !*** renormalize weights sum_wgts = sum(src_lats) wgts(1) = src_lats(1)/sum_wgts wgts(2) = src_lats(2)/sum_wgts wgts(3) = src_lats(3)/sum_wgts wgts(4) = src_lats(4)/sum_wgts grid1_frac(dst_add) = one call store_link_bilin(dst_add, src_add, wgts, nmap) endif endif end do grid_loop2 endif ! nmap=2 !----------------------------------------------------------------------- end subroutine remap_bilin !*********************************************************************** subroutine grid_search_bilin(src_add, src_lats, src_lons, & plat, plon, src_grid_dims, & src_center_lat, src_center_lon, & src_grid_bound_box, & src_bin_add, dst_bin_add) !----------------------------------------------------------------------- ! ! this routine finds the location of the search point plat, plon ! in the source grid and returns the corners needed for a bilinear ! interpolation. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! output variables ! !----------------------------------------------------------------------- integer (kind=int_kind), dimension(4), intent(out) :: & src_add ! address of each corner point enclosing P real (kind=dbl_kind), dimension(4), intent(out) :: & src_lats, & ! latitudes of the four corner points src_lons ! longitudes of the four corner points !----------------------------------------------------------------------- ! ! input variables ! !----------------------------------------------------------------------- real (kind=dbl_kind), intent(in) :: & plat, & ! latitude of the search point plon ! longitude of the search point integer (kind=int_kind), dimension(2), intent(in) :: & src_grid_dims ! size of each src grid dimension real (kind=dbl_kind), dimension(:), intent(in) :: & src_center_lat, & ! latitude of each src grid center src_center_lon ! longitude of each src grid center real (kind=dbl_kind), dimension(:,:), intent(in) :: & src_grid_bound_box ! bound box for source grid integer (kind=int_kind), dimension(:,:), intent(in) :: & src_bin_add, & ! latitude bins for restricting dst_bin_add ! searches !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (kind=int_kind) :: n, next_n, srch_add, & ! dummy indices nx, ny, & ! dimensions of src grid min_add, max_add, & ! addresses for restricting search i, j, jp1, ip1, n_add, e_add, ne_add ! addresses real (kind=dbl_kind) :: & ! vectors for cross-product check vec1_lat, vec1_lon, & vec2_lat, vec2_lon, cross_product, cross_product_last, & coslat_dst, sinlat_dst, coslon_dst, sinlon_dst, & dist_min, distance ! for computing dist-weighted avg !----------------------------------------------------------------------- ! ! restrict search first using bins ! !----------------------------------------------------------------------- src_add = 0 min_add = size(src_center_lat) max_add = 1 do n=1,num_srch_bins if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and. & plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then min_add = min(min_add, src_bin_add(1,n)) max_add = max(max_add, src_bin_add(2,n)) endif end do !----------------------------------------------------------------------- ! ! now perform a more detailed search ! !----------------------------------------------------------------------- nx = src_grid_dims(1) ny = src_grid_dims(2) srch_loop: do srch_add = min_add,max_add !*** first check bounding box if (plat <= src_grid_bound_box(2,srch_add) .and. & plat >= src_grid_bound_box(1,srch_add) .and. & plon <= src_grid_bound_box(4,srch_add) .and. & plon >= src_grid_bound_box(3,srch_add)) then !*** !*** we are within bounding box so get really serious !*** !*** determine neighbor addresses j = (srch_add - 1)/nx +1 i = srch_add - (j-1)*nx if (i < nx) then ip1 = i + 1 else ip1 = 1 endif if (j < ny) then jp1 = j+1 else jp1 = 1 endif n_add = (jp1 - 1)*nx + i e_add = (j - 1)*nx + ip1 ne_add = (jp1 - 1)*nx + ip1 src_lats(1) = src_center_lat(srch_add) src_lats(2) = src_center_lat(e_add) src_lats(3) = src_center_lat(ne_add) src_lats(4) = src_center_lat(n_add) src_lons(1) = src_center_lon(srch_add) src_lons(2) = src_center_lon(e_add) src_lons(3) = src_center_lon(ne_add) src_lons(4) = src_center_lon(n_add) !*** !*** for consistency, we must make sure all lons are in !*** same 2pi interval !*** vec1_lon = src_lons(1) - plon if (vec1_lon > pi) then src_lons(1) = src_lons(1) - pi2 else if (vec1_lon < -pi) then src_lons(1) = src_lons(1) + pi2 endif do n=2,4 vec1_lon = src_lons(n) - src_lons(1) if (vec1_lon > pi) then src_lons(n) = src_lons(n) - pi2 else if (vec1_lon < -pi) then src_lons(n) = src_lons(n) + pi2 endif end do corner_loop: do n=1,4 next_n = MOD(n,4) + 1 !*** !*** here we take the cross product of the vector making !*** up each box side with the vector formed by the vertex !*** and search point. if all the cross products are !*** positive, the point is contained in the box. !*** vec1_lat = src_lats(next_n) - src_lats(n) vec1_lon = src_lons(next_n) - src_lons(n) vec2_lat = plat - src_lats(n) vec2_lon = plon - src_lons(n) !*** !*** check for 0,2pi crossings !*** if (vec1_lon > three*pih) then vec1_lon = vec1_lon - pi2 else if (vec1_lon < -three*pih) then vec1_lon = vec1_lon + pi2 endif if (vec2_lon > three*pih) then vec2_lon = vec2_lon - pi2 else if (vec2_lon < -three*pih) then vec2_lon = vec2_lon + pi2 endif cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat !*** !*** if cross product is less than zero, this cell !*** doesn't work !*** if (n == 1) cross_product_last = cross_product if (cross_product*cross_product_last < zero) & exit corner_loop cross_product_last = cross_product end do corner_loop !*** !*** if cross products all same sign, we found the location !*** if (n > 4) then src_add(1) = srch_add src_add(2) = e_add src_add(3) = ne_add src_add(4) = n_add return endif !*** !*** otherwise move on to next cell !*** endif !bounding box check end do srch_loop !*** !*** if no cell found, point is likely either in a box that !*** straddles either pole or is outside the grid. fall back !*** to a distance-weighted average of the four closest !*** points. go ahead and compute weights here, but store !*** in src_lats and return -add to prevent the parent !*** routine from computing bilinear weights !*** !print *,'Could not find location for ',plat,plon !print *,'Using nearest-neighbor average for this point' coslat_dst = cos(plat) sinlat_dst = sin(plat) coslon_dst = cos(plon) sinlon_dst = sin(plon) dist_min = bignum src_lats = bignum do srch_add = min_add,max_add distance = acos(coslat_dst*cos(src_center_lat(srch_add))* & (coslon_dst*cos(src_center_lon(srch_add)) + & sinlon_dst*sin(src_center_lon(srch_add)))+ & sinlat_dst*sin(src_center_lat(srch_add))) if (distance < dist_min) then sort_loop: do n=1,4 if (distance < src_lats(n)) then do i=4,n+1,-1 src_add (i) = src_add (i-1) src_lats(i) = src_lats(i-1) end do src_add (n) = -srch_add src_lats(n) = distance dist_min = src_lats(4) exit sort_loop endif end do sort_loop endif end do src_lons = one/(src_lats + tiny) distance = sum(src_lons) src_lats = src_lons/distance !----------------------------------------------------------------------- end subroutine grid_search_bilin !*********************************************************************** subroutine store_link_bilin(dst_add, src_add, weights, nmap) !----------------------------------------------------------------------- ! ! this routine stores the address and weight for four links ! associated with one destination point in the appropriate address ! and weight arrays and resizes those arrays if necessary. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! input variables ! !----------------------------------------------------------------------- integer (kind=int_kind), intent(in) :: & dst_add, & ! address on destination grid nmap ! identifies which direction for mapping integer (kind=int_kind), dimension(4), intent(in) :: & src_add ! addresses on source grid real (kind=dbl_kind), dimension(4), intent(in) :: & weights ! array of remapping weights for these links !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (kind=int_kind) :: n, & ! dummy index num_links_old ! placeholder for old link number !----------------------------------------------------------------------- ! ! increment number of links and check to see if remap arrays need ! to be increased to accomodate the new link. then store the ! link. ! !----------------------------------------------------------------------- select case (nmap) case(1) num_links_old = num_links_map1 num_links_map1 = num_links_old + 4 if (num_links_map1 > max_links_map1) & call resize_remap_vars(1,resize_increment) do n=1,4 grid1_add_map1(num_links_old+n) = src_add(n) grid2_add_map1(num_links_old+n) = dst_add wts_map1 (1,num_links_old+n) = weights(n) end do case(2) num_links_old = num_links_map2 num_links_map2 = num_links_old + 4 if (num_links_map2 > max_links_map2) & call resize_remap_vars(2,resize_increment) do n=1,4 grid1_add_map2(num_links_old+n) = dst_add grid2_add_map2(num_links_old+n) = src_add(n) wts_map2 (1,num_links_old+n) = weights(n) end do end select !----------------------------------------------------------------------- end subroutine store_link_bilin !*********************************************************************** end module remap_bilinear !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!