From d1ce50ad463a0de402a2078c18d46d5046f956ba Mon Sep 17 00:00:00 2001 From: Ben Johnson Date: Wed, 6 Nov 2024 13:38:58 -0700 Subject: [PATCH] Separate subroutines for barycentric interpolation into barycentric_utilities_mod.f90 --- models/utilities/barycentric_utils_mod.f90 | 201 +++++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 models/utilities/barycentric_utils_mod.f90 diff --git a/models/utilities/barycentric_utils_mod.f90 b/models/utilities/barycentric_utils_mod.f90 new file mode 100644 index 0000000000..495da9aff5 --- /dev/null +++ b/models/utilities/barycentric_utils_mod.f90 @@ -0,0 +1,201 @@ +!> Barycentric interpolation routines for triangular structures within grids. These were originally +!> Developed for the mpas_atm model_mod. They're being separated into their own utilities module +!> so they can be used by multiple model_mods. + +module barycentric_utils_mod + + use types_mod, only : r8, DEG2RAD, radius => earth_radius + + implicit none + private + + public :: inside_triangle, & + get_3d_weights, & + get_barycentric_weights, & + latlon_to_xyz, & + barycentric_average + + real(r8), parameter :: roundoff = 1.0e-12_r8 + + contains + + subroutine inside_triangle(t1, t2, t3, r, lat, lon, inside, weights) + + ! given 3 corners of a triangle and an xyz point, compute whether + ! the point is inside the triangle. this assumes r is coplanar + ! with the triangle - the caller must have done the lat/lon to + ! xyz conversion with a constant radius and then this will be + ! true (enough). sets inside to true/false, and returns the + ! weights if true. weights are set to 0 if false. + + real(r8), intent(in) :: t1(3), t2(3), t3(3) + real(r8), intent(in) :: r(3), lat, lon + logical, intent(out) :: inside + real(r8), intent(out) :: weights(3) + + ! check for degenerate cases first - is the test point located + ! directly on one of the vertices? (this case may be common + ! if we're computing on grid point locations.) + if (all(abs(r - t1) < roundoff)) then + inside = .true. + weights = (/ 1.0_r8, 0.0_r8, 0.0_r8 /) + return + else if (all(abs(r - t2) < roundoff)) then + inside = .true. + weights = (/ 0.0_r8, 1.0_r8, 0.0_r8 /) + return + else if (all(abs(r - t3) < roundoff)) then + inside = .true. + weights = (/ 0.0_r8, 0.0_r8, 1.0_r8 /) + return + endif + + ! not a vertex. compute the weights. if any are + ! negative, the point is outside. since these are + ! real valued computations define a lower bound for + ! numerical roundoff error and be sure that the + ! weights are not just *slightly* negative. + call get_3d_weights(r, t1, t2, t3, lat, lon, weights) + + if (any(weights < -roundoff)) then + inside = .false. + weights = 0.0_r8 + return + endif + + ! truncate barely negative values to 0 + inside = .true. + where (weights < 0.0_r8) weights = 0.0_r8 + return + + end subroutine inside_triangle + + subroutine get_3d_weights(p, v1, v2, v3, lat, lon, weights) + + ! Given a point p (x,y,z) inside a triangle, and the (x,y,z) + ! coordinates of the triangle corner points (v1, v2, v3), + ! find the weights for a barycentric interpolation. this + ! computation only needs two of the three coordinates, so figure + ! out which quadrant of the sphere the triangle is in and pick + ! the 2 axes which are the least planar: + ! (x,y) near the poles, + ! (y,z) near 0 and 180 longitudes near the equator, + ! (x,z) near 90 and 270 longitude near the equator. + ! (lat/lon are the coords of p. we could compute them here + ! but since in all cases we already have them, pass them + ! down for efficiency) + + real(r8), intent(in) :: p(3) + real(r8), intent(in) :: v1(3), v2(3), v3(3) + real(r8), intent(in) :: lat, lon + real(r8), intent(out) :: weights(3) + + real(r8) :: cxs(3), cys(3) + + ! above or below 45 in latitude, where -90 < lat < 90: + if (lat >= 45.0_r8 .or. lat <= -45.0_r8) then + cxs(1) = v1(1) + cxs(2) = v2(1) + cxs(3) = v3(1) + cys(1) = v1(2) + cys(2) = v2(2) + cys(3) = v3(2) + call get_barycentric_weights(p(1), p(2), cxs, cys, weights) + return + endif + + ! nearest 0 or 180 in longitude, where 0 < lon < 360: + if ( lon <= 45.0_r8 .or. lon >= 315.0_r8 .or. & + (lon >= 135.0_r8 .and. lon <= 225.0_r8)) then + cxs(1) = v1(2) + cxs(2) = v2(2) + cxs(3) = v3(2) + cys(1) = v1(3) + cys(2) = v2(3) + cys(3) = v3(3) + call get_barycentric_weights(p(2), p(3), cxs, cys, weights) + return + endif + + ! last option, nearest 90 or 270 in lon: + cxs(1) = v1(1) + cxs(2) = v2(1) + cxs(3) = v3(1) + cys(1) = v1(3) + cys(2) = v2(3) + cys(3) = v3(3) + call get_barycentric_weights(p(1), p(3), cxs, cys, weights) + + end subroutine get_3d_weights + + subroutine get_barycentric_weights(x, y, cxs, cys, weights) + + ! Computes the barycentric weights for a 2d interpolation point + ! (x,y) in a 2d triangle with the given (cxs,cys) corners. + + real(r8), intent(in) :: x, y, cxs(3), cys(3) + real(r8), intent(out) :: weights(3) + + real(r8) :: denom + + ! Get denominator + denom = (cys(2) - cys(3)) * (cxs(1) - cxs(3)) + & + (cxs(3) - cxs(2)) * (cys(1) - cys(3)) + + weights(1) = ((cys(2) - cys(3)) * (x - cxs(3)) + & + (cxs(3) - cxs(2)) * (y - cys(3))) / denom + + weights(2) = ((cys(3) - cys(1)) * (x - cxs(3)) + & + (cxs(1) - cxs(3)) * (y - cys(3))) / denom + + weights(3) = 1.0_r8 - weights(1) - weights(2) + + if (any(abs(weights) < roundoff)) then + where (abs(weights) < roundoff) weights = 0.0_r8 + where (abs(1.0_r8 - abs(weights)) < roundoff) weights = 1.0_r8 + endif + + end subroutine get_barycentric_weights + + subroutine latlon_to_xyz(lat, lon, x, y, z) + + ! Given a lat, lon in degrees, return the cartesian x,y,z coordinate + ! on the surface of a specified radius relative to the origin + ! at the center of the earth. + + real(r8), intent(in) :: lat, lon + real(r8), intent(out) :: x, y, z + + real(r8) :: rlat, rlon + + rlat = lat * DEG2RAD + rlon = lon * DEG2RAD + + x = radius * cos(rlon) * cos(rlat) + y = radius * sin(rlon) * cos(rlat) + z = radius * sin(rlat) + + end subroutine latlon_to_xyz + + function barycentric_average(nitems, weights, vertex_temp_values) result (averaged_values) + + integer, intent(in) :: nitems + real(r8), dimension(3), intent(in) :: weights + real(r8), dimension(3, nitems), intent(in) :: vertex_temp_values + + real(r8), dimension(nitems) :: averaged_values + + integer :: iweight, iitem + + averaged_values(:) = 0 + + do iitem = 1, nitems + do iweight = 1, 3 + averaged_values(iitem) = averaged_values(iitem) + weights(iweight)*vertex_temp_values(iweight, iitem) + end do + end do + + end function barycentric_average + +end module barycentric_utils_mod + \ No newline at end of file