From 3d05ddd05a35e2f65adb7f9f49440b8d82bb709e Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Thu, 18 Nov 2021 16:28:00 -0500 Subject: [PATCH 1/9] Fixes #1199 --- CHANGELOG.md | 5 +- base/HorizontalFluxRegridder.F90 | 66 +++- base/MAPL_SphericalGeometry.F90 | 607 ++++++++++++++++--------------- 3 files changed, 372 insertions(+), 306 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5fc7fb8043e7..15ad87898131 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,10 +12,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed ### Fixed +- Corrected bug in HorizontalFluxRegridder. Fluxes need to be + multiplied by edge length for correct treatment. -### Removed - -### Deprecated ## [2.37.0] - 2023-04-03 diff --git a/base/HorizontalFluxRegridder.F90 b/base/HorizontalFluxRegridder.F90 index 0cb8fadc79f2..5d355ed03795 100644 --- a/base/HorizontalFluxRegridder.F90 +++ b/base/HorizontalFluxRegridder.F90 @@ -10,6 +10,7 @@ module mapl_HorizontalFluxRegridder use mapl_KeywordEnforcerMod use mapl_ErrorHandlingMod use mapl_BaseMod + use mapl_SphericalGeometry implicit none private @@ -20,6 +21,8 @@ module mapl_HorizontalFluxRegridder integer :: resolution_ratio = -1 integer :: im_in, jm_in integer :: im_out, jm_out + real, allocatable :: dx_in(:,:), dy_in(:,:) + real, allocatable :: dx_out(:,:), dy_out(:,:) contains procedure, nopass :: supports procedure :: initialize_subclass @@ -70,6 +73,8 @@ subroutine initialize_subclass(this, unusable, rc) integer :: counts(5) integer :: status + integer :: units ! unused + real, pointer :: lons(:,:), lats(:,:) _UNUSED_DUMMY(unusable) spec = this%get_spec() @@ -91,8 +96,36 @@ subroutine initialize_subclass(this, unusable, rc) _ASSERT((IM_in / IM_out) == (JM_in / JM_out), 'inconsistent aspect ratio') this%resolution_ratio = (IM_in / IM_out) + + call ESMF_GridGetCoord(grid_in, coordDim=1, farrayPtr=lons, & + localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__) + call ESMF_GridGetCoord(grid_in, coordDim=2, farrayPtr=lats, & + localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__) + + this%dx_in = distance( & + lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & + lons(2:IM_in+1,1:JM_in), lats(2:IM_in+1,1:JM_in)) + + this%dy_in = distance( & + lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & + lons(1:IM_in,2:JM_in+1), lats(1:IM_in,2:JM_in+1)) + + call ESMF_GridGetCoord(grid_out, coordDim=1, farrayPtr=lons, & + localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__) + call ESMF_GridGetCoord(grid_out, coordDim=2, farrayPtr=lats, & + localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__) + + this%dx_out = distance( & + lons(1:IM_out,1:JM_out), lats(1:IM_out,1:JM_out), & + lons(2:IM_out+1,1:JM_out), lats(2:IM_out+1,1:JM_out)) + + this%dy_out = distance( & + lons(1:IM_out,1:JM_out), lats(1:IM_out,1:JM_out), & + lons(1:IM_out,2:JM_out+1), lats(1:IM_out,2:JM_out+1)) + end associate end associate + _RETURN(_SUCCESS) end subroutine initialize_subclass @@ -129,9 +162,14 @@ subroutine regrid_vector_2d_real32(this, u_in, v_in, u_out, v_out, rotate, rc) do i = 1, IM m_y = 0 do ii = 1 + (i-1)*N, i*N - m_y = m_y + v_in(ii,jj) + associate (d_in => this%dx_in(ii,jj)) + m_y = m_y + v_in(ii,jj) * d_in + end associate end do - v_out(i,j) = m_y + + associate (d_out => this%dx_out(i,j)) + v_out(i,j) = m_y / d_out + end associate end do end do @@ -141,9 +179,13 @@ subroutine regrid_vector_2d_real32(this, u_in, v_in, u_out, v_out, rotate, rc) do j = 1, JM m_x = 0 do jj = 1 + (j-1)*N, j*N - m_x = m_x + u_in(ii,jj) + associate (d_in => this%dy_in(ii,jj)) + m_x = m_x + u_in(ii,jj) * d_in + end associate end do - u_out(i,j) = m_x + associate (d_out => this%dy_out(i,j)) + u_out(i,j) = m_x / d_out + end associate end do end do @@ -186,9 +228,13 @@ subroutine regrid_vector_2d_real64(this, u_in, v_in, u_out, v_out, rotate, rc) do i = 1, IM m_y = 0 do ii = 1 + (i-1)*N, i*N - m_y = m_y + v_in(ii,jj) + associate (d_in => this%dx_in(ii,jj)) + m_y = m_y + v_in(ii,jj) * d_in + end associate end do - v_out(i,j) = m_y + associate (d_out => this%dx_out(i,j)) + v_out(i,j) = m_y / d_out + end associate end do end do @@ -198,9 +244,13 @@ subroutine regrid_vector_2d_real64(this, u_in, v_in, u_out, v_out, rotate, rc) do j = 1, JM m_x = 0 do jj = 1 + (j-1)*N, j*N - m_x = m_x + u_in(ii,jj) + associate (d_in => this%dy_in(ii,jj)) + m_x = m_x + u_in(ii,jj) * d_in + end associate end do - u_out(i,j) = m_x + associate (d_out => this%dy_out(i,j)) + u_out(i,j) = m_x / d_out + end associate end do end do diff --git a/base/MAPL_SphericalGeometry.F90 b/base/MAPL_SphericalGeometry.F90 index 7e61564bb072..f901cf1383df 100644 --- a/base/MAPL_SphericalGeometry.F90 +++ b/base/MAPL_SphericalGeometry.F90 @@ -1,307 +1,324 @@ #include "MAPL_ErrLog.h" -module MAPL_SphericalGeometry - use MAPL_KeywordEnforcerMod +module mapl_SphericalGeometry + use mapl_KeywordEnforcerMod use mapl_ErrorHandlingMod use ESMF - use MAPL_Constants + use mapl_Constants use, intrinsic :: iso_fortran_env, only: REAL64,REAL32 + implicit none + private + + public get_points_in_spherical_domain + public get_area_spherical_polygon + public :: distance -implicit none -private -public get_points_in_spherical_domain -public get_area_spherical_polygon contains - ! get area of spherical rectangle given the four corners - ! p4 ------ p3 - ! | | - ! | | - ! | | - ! p1 ------ p2 - function get_area_spherical_polygon(p1,p4,p2,p3) result(area) - real(real64) :: area - real(real64), intent(in) :: p1(2),p2(2),p3(2),p4(2) - - real(real64) :: e1(3),e2(3),e3(3) - real(real64) :: ang1,ang2,ang3,ang4 - - e1 = convert_to_cart(p1) - e2 = convert_to_cart(p2) - e3 = convert_to_cart(p4) - ang1 = spherical_angles(e1, e2, e3) - - e1 = convert_to_cart(p2) - e2 = convert_to_cart(p3) - e3 = convert_to_cart(p1) - ang2 = spherical_angles(e1, e2, e3) - - e1 = convert_to_cart(p3) - e2 = convert_to_cart(p4) - e3 = convert_to_cart(p2) - ang3 = spherical_angles(e1, e2, e3) - - e1 = convert_to_cart(p4) - e2 = convert_to_cart(p3) - e3 = convert_to_cart(p1) - ang4 = spherical_angles(e1, e2, e3) - - area = ang1 + ang2 + ang3 + ang4 - 2.0d0*MAPL_PI_R8 - - end function get_area_spherical_polygon - - subroutine get_points_in_spherical_domain(center_lons,center_lats,corner_lons,corner_lats,lons,lats,ii,jj,rc) - real(real64), intent(in) :: center_lats(:,:),center_lons(:,:) - real(real64), intent(in) :: corner_lats(:,:),corner_lons(:,:) - real(real64), intent(in) :: lons(:),lats(:) - integer, intent(out) :: ii(:),jj(:) - integer, intent(out), optional :: rc - - integer :: npts,i,n,niter,im,jm,ilb,jlb,iub,jub,ifound,jfound - integer :: lold,uold,lnew,unew - logical :: in_region,in_sub_region - - npts = size(lats) - - _ASSERT(size(lats)==size(lons),"lats and lons do not match") - _ASSERT(npts==size(ii),"size of ii does not match") - _ASSERT(npts==size(ii),"size of jj does not match") - - im=size(corner_lons,1)-1 - jm=size(corner_lons,2)-1 - niter = max(im,jm) - - do i=1,npts - ifound=-1 - jfound=-1 - ilb=1 - iub=im - jlb=1 - jub=jm - in_region = point_in_polygon([lons(i),lats(i)],[center_lons(ilb,jlb),center_lats(ilb,jlb)], & - [corner_lons(ilb,jlb),corner_lats(ilb,jlb)], & - [corner_lons(iub+1,jlb),corner_lats(iub+1,jlb)], & - [corner_lons(iub+1,jub+1),corner_lats(iub+1,jub+1)], & - [corner_lons(ilb,jub+1),corner_lats(ilb,jub+1)]) - if (in_region) then - ! bisect first dimension - lnew=ilb - unew=iub - do n = 1,niter - lold=lnew - uold=unew - unew=lold+(uold-lold)/2 - in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(lnew,jlb),center_lats(lnew,jlb)], & - [corner_lons(lnew,jlb),corner_lats(lnew,jlb)], & - [corner_lons(unew+1,jlb),corner_lats(unew+1,jlb)], & - [corner_lons(unew+1,jub+1),corner_lats(unew+1,jub+1)], & - [corner_lons(lnew,jub+1),corner_lats(lnew,jub+1)]) - if (in_sub_region) then - lnew=lold - unew=unew - else - lnew=unew+1 - unew=uold - end if - if (unew==lnew) then - ifound=unew - exit - end if - enddo - ! bisect 2nd dimension - lnew=jlb - unew=jub - do n = 1,niter - lold=lnew - uold=unew - unew=lold+(uold-lold)/2 - in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(ifound,lnew),center_lats(ifound,lnew)] , & - [corner_lons(ifound,lnew),corner_lats(ifound,lnew)], & - [corner_lons(ifound+1,lnew),corner_lats(ifound+1,lnew)], & - [corner_lons(ifound+1,unew+1),corner_lats(ifound+1,unew+1)], & - [corner_lons(ifound,unew+1),corner_lats(ifound,unew+1)]) - if (in_sub_region) then - lnew=lold - unew=unew - else - lnew=unew+1 - unew=uold - end if - if (unew==lnew) then - jfound=unew - exit - end if - enddo - end if - ii(i)=ifound - jj(i)=jfound - enddo - _RETURN(_SUCCESS) - - end subroutine get_points_in_spherical_domain - - function point_in_polygon(p0,pinside,a1,a2,a3,a4) result(in_poly) - real(real64), intent(in) :: p0(2),pinside(2),a1(2),a2(2),a3(2),a4(2) - logical :: in_poly - - real(real64) :: p1c(3),p2c(3),a1c(3),a2c(3),a3c(3),a4c(3) - logical :: intersect(4) - p1c=convert_to_cart(p0) - p2c=convert_to_cart(pinside) - a1c=convert_to_cart(a1) - a2c=convert_to_cart(a2) - a3c=convert_to_cart(a3) - a4c=convert_to_cart(a4) - - intersect(1) = lines_intersect(p1c,p2c,a1c,a2c) - intersect(2) = lines_intersect(p1c,p2c,a2c,a3c) - intersect(3) = lines_intersect(p1c,p2c,a3c,a4c) - intersect(4) = lines_intersect(p1c,p2c,a4c,a1c) - if (mod(count(intersect),2)==0) then - in_poly=.true. - else - in_poly=.false. - end if - - - end function point_in_polygon - -! it is claimed this should work but doesn't - !function point_in_polygon_crosprod(p1,a1,a2,a3,a4) result(in_poly) - !real(real64), intent(in) :: p1(2),a1(2),a2(2),a3(2),a4(2) - !logical :: in_poly - - !real(real64) :: p1c(3),a1c(3),a2c(3),a3c(3),a4c(3) - !real(real64) :: crs12(3),crs23(3),crs34(3),crs41(3) - !real(real64) :: d12,d23,d34,d41 - !logical :: signs(4) - !! a1 -> a2 -> a3 -> a4 so a4 connects to a1 - - !p1c=convert_to_cart(p1) - !a1c=convert_to_cart(a1) - !a2c=convert_to_cart(a2) - !a3c=convert_to_cart(a3) - !a4c=convert_to_cart(a4) - - !crs12 = cross_prod(a1c,a2c) - !crs23 = cross_prod(a2c,a3c) - !crs34 = cross_prod(a3c,a4c) - !crs41 = cross_prod(a4c,a1c) - !d12=dot_product(p1c,crs12) - !d23=dot_product(p1c,crs23) - !d34=dot_product(p1c,crs34) - !d41=dot_product(p1c,crs41) - !signs(1)= (d12<0.0) - !signs(2)= (d23<0.0) - !signs(3)= (d34<0.0) - !signs(4)= (d41<0.0) - !in_poly=( (count(signs)==0) .or. (count(signs)==4) ) - - !end function point_in_polygon_crossprod - - function lines_intersect(b0,b1,a0,a1) result(intersect) - real(real64), intent(in) :: b0(3),b1(3),a0(3),a1(3) - logical :: intersect - real(real64) :: p(3),q(3),t(3) - real(real64) :: s1,s2,s3,s4 - logical :: signs(4) - - intersect=.false. - q=cross_prod(b0,b1) - p=cross_prod(a0,a1) - t=normal_vect(cross_prod(p,q)) - - s1=dot_product(cross_prod(a0,p),t) - s2=dot_product(cross_prod(a1,p),t) - s3=dot_product(cross_prod(b0,q),t) - s4=dot_product(cross_prod(b1,q),t) - - signs(1) = -s1 <0.d0 - signs(2) = s2 <0.d0 - signs(3) = -s3 < 0.d0 - signs(4) = s4 < 0.d0 - - intersect = ((count(signs)==0) .or. (count(signs)==4)) - - end function lines_intersect - - function normal_vect(vin) result(vout) - real(real64), intent(in) :: vin(3) - real(real64) :: vout(3) - vout=vin/sqrt(vin(1)*vin(1)+vin(2)*vin(2)+vin(3)*vin(3)) - - end function normal_vect - - function cross_prod(v1,v2) result(vout) - real(real64), intent(in) :: v1(3),v2(3) - real(real64) :: vout(3) - vout(1)=v1(2)*v2(3)-v1(3)*v2(2) - vout(2)=v1(3)*v2(1)-v1(1)*v2(3) - vout(3)=v1(1)*v2(2)-v1(2)*v2(1) - end function cross_prod - - function convert_to_cart(v) result(xyz) - real(real64), intent(in) :: v(2) - real(real64) :: xyz(3) - - xyz(1)=cos(v(2))*cos(v(1)) - xyz(2)=cos(v(2))*sin(v(1)) - xyz(3)=sin(v(2)) - - end function convert_to_cart - -function vect_mag(v) result(mag) - real(real64), intent(in) :: v(3) - real :: mag - mag = sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) -end function vect_mag - -function spherical_angles(p1,p2,p3) result(spherical_angle) - real(real64) :: spherical_angle - real(real64), intent(in) :: p1(3),p2(3),p3(3) - - real (real64):: e1(3), e2(3), e3(3) - real (real64):: px, py, pz - real (real64):: qx, qy, qz - real (real64):: angle, ddd - integer n - - do n=1,3 - e1(n) = p1(n) - e2(n) = p2(n) - e3(n) = p3(n) - enddo - - !------------------------------------------------------------------- - ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry - !------------------------------------------------------------------- - ! Vector P: - px = e1(2)*e2(3) - e1(3)*e2(2) - py = e1(3)*e2(1) - e1(1)*e2(3) - pz = e1(1)*e2(2) - e1(2)*e2(1) - ! Vector Q: - qx = e1(2)*e3(3) - e1(3)*e3(2) - qy = e1(3)*e3(1) - e1(1)*e3(3) - qz = e1(1)*e3(2) - e1(2)*e3(1) - - ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz) - - if ( ddd <= 0.0d0 ) then - angle = 0.d0 - else - ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd) - if ( abs(ddd)>1.d0) then - angle = 0.5d0 * MAPL_PI_R8 - if (ddd < 0.d0) then - angle = MAPL_PI_R8 - else - angle = 0.d0 + ! get area of spherical rectangle given the four corners + ! p4 ------ p3 + ! | | + ! | | + ! | | + ! p1 ------ p2 + function get_area_spherical_polygon(p1,p4,p2,p3) result(area) + real(real64) :: area + real(real64), intent(in) :: p1(2),p2(2),p3(2),p4(2) + + real(real64) :: e1(3),e2(3),e3(3) + real(real64) :: ang1,ang2,ang3,ang4 + + e1 = convert_to_cart(p1) + e2 = convert_to_cart(p2) + e3 = convert_to_cart(p4) + ang1 = spherical_angles(e1, e2, e3) + + e1 = convert_to_cart(p2) + e2 = convert_to_cart(p3) + e3 = convert_to_cart(p1) + ang2 = spherical_angles(e1, e2, e3) + + e1 = convert_to_cart(p3) + e2 = convert_to_cart(p4) + e3 = convert_to_cart(p2) + ang3 = spherical_angles(e1, e2, e3) + + e1 = convert_to_cart(p4) + e2 = convert_to_cart(p3) + e3 = convert_to_cart(p1) + ang4 = spherical_angles(e1, e2, e3) + + area = ang1 + ang2 + ang3 + ang4 - 2.0d0*MAPL_PI_R8 + + end function get_area_spherical_polygon + + ! Computes distance between two points (in lat lon as radians), + ! and returns distance in radians (unit sphere) + ! Using formulae from: https://www.movable-type.co.uk/scripts/latlong.html + + elemental real function distance(lon1, lat1, lon2, lat2) result(d) + real, intent(in) :: lon1, lat1 + real, intent(in) :: lon2, lat2 + + associate(a => sin(lat2-lat1)**2 + cos(lat1)*cos(lat2)*sin((lon2-lon1)/2)**2) + d = 2*atan2(sqrt(a), sqrt(1-a)) + end associate + + end function distance + + subroutine get_points_in_spherical_domain(center_lons,center_lats,corner_lons,corner_lats,lons,lats,ii,jj,rc) + real(real64), intent(in) :: center_lats(:,:),center_lons(:,:) + real(real64), intent(in) :: corner_lats(:,:),corner_lons(:,:) + real(real64), intent(in) :: lons(:),lats(:) + integer, intent(out) :: ii(:),jj(:) + integer, intent(out), optional :: rc + + integer :: npts,i,n,niter,im,jm,ilb,jlb,iub,jub,ifound,jfound + integer :: lold,uold,lnew,unew + logical :: in_region,in_sub_region + + npts = size(lats) + + _ASSERT(size(lats)==size(lons),"lats and lons do not match") + _ASSERT(npts==size(ii),"size of ii does not match") + _ASSERT(npts==size(ii),"size of jj does not match") + + im=size(corner_lons,1)-1 + jm=size(corner_lons,2)-1 + niter = max(im,jm) + + do i=1,npts + ifound=-1 + jfound=-1 + ilb=1 + iub=im + jlb=1 + jub=jm + in_region = point_in_polygon([lons(i),lats(i)],[center_lons(ilb,jlb),center_lats(ilb,jlb)], & + [corner_lons(ilb,jlb),corner_lats(ilb,jlb)], & + [corner_lons(iub+1,jlb),corner_lats(iub+1,jlb)], & + [corner_lons(iub+1,jub+1),corner_lats(iub+1,jub+1)], & + [corner_lons(ilb,jub+1),corner_lats(ilb,jub+1)]) + if (in_region) then + ! bisect first dimension + lnew=ilb + unew=iub + do n = 1,niter + lold=lnew + uold=unew + unew=lold+(uold-lold)/2 + in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(lnew,jlb),center_lats(lnew,jlb)], & + [corner_lons(lnew,jlb),corner_lats(lnew,jlb)], & + [corner_lons(unew+1,jlb),corner_lats(unew+1,jlb)], & + [corner_lons(unew+1,jub+1),corner_lats(unew+1,jub+1)], & + [corner_lons(lnew,jub+1),corner_lats(lnew,jub+1)]) + if (in_sub_region) then + lnew=lold + unew=unew + else + lnew=unew+1 + unew=uold + end if + if (unew==lnew) then + ifound=unew + exit + end if + enddo + ! bisect 2nd dimension + lnew=jlb + unew=jub + do n = 1,niter + lold=lnew + uold=unew + unew=lold+(uold-lold)/2 + in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(ifound,lnew),center_lats(ifound,lnew)] , & + [corner_lons(ifound,lnew),corner_lats(ifound,lnew)], & + [corner_lons(ifound+1,lnew),corner_lats(ifound+1,lnew)], & + [corner_lons(ifound+1,unew+1),corner_lats(ifound+1,unew+1)], & + [corner_lons(ifound,unew+1),corner_lats(ifound,unew+1)]) + if (in_sub_region) then + lnew=lold + unew=unew + else + lnew=unew+1 + unew=uold + end if + if (unew==lnew) then + jfound=unew + exit + end if + enddo end if + ii(i)=ifound + jj(i)=jfound + enddo + _RETURN(_SUCCESS) + + end subroutine get_points_in_spherical_domain + + function point_in_polygon(p0,pinside,a1,a2,a3,a4) result(in_poly) + real(real64), intent(in) :: p0(2),pinside(2),a1(2),a2(2),a3(2),a4(2) + logical :: in_poly + + real(real64) :: p1c(3),p2c(3),a1c(3),a2c(3),a3c(3),a4c(3) + logical :: intersect(4) + p1c=convert_to_cart(p0) + p2c=convert_to_cart(pinside) + a1c=convert_to_cart(a1) + a2c=convert_to_cart(a2) + a3c=convert_to_cart(a3) + a4c=convert_to_cart(a4) + + intersect(1) = lines_intersect(p1c,p2c,a1c,a2c) + intersect(2) = lines_intersect(p1c,p2c,a2c,a3c) + intersect(3) = lines_intersect(p1c,p2c,a3c,a4c) + intersect(4) = lines_intersect(p1c,p2c,a4c,a1c) + if (mod(count(intersect),2)==0) then + in_poly=.true. + else + in_poly=.false. + end if + + + end function point_in_polygon + + ! it is claimed this should work but doesn't + !function point_in_polygon_crosprod(p1,a1,a2,a3,a4) result(in_poly) + !real(real64), intent(in) :: p1(2),a1(2),a2(2),a3(2),a4(2) + !logical :: in_poly + + !real(real64) :: p1c(3),a1c(3),a2c(3),a3c(3),a4c(3) + !real(real64) :: crs12(3),crs23(3),crs34(3),crs41(3) + !real(real64) :: d12,d23,d34,d41 + !logical :: signs(4) + !! a1 -> a2 -> a3 -> a4 so a4 connects to a1 + + !p1c=convert_to_cart(p1) + !a1c=convert_to_cart(a1) + !a2c=convert_to_cart(a2) + !a3c=convert_to_cart(a3) + !a4c=convert_to_cart(a4) + + !crs12 = cross_prod(a1c,a2c) + !crs23 = cross_prod(a2c,a3c) + !crs34 = cross_prod(a3c,a4c) + !crs41 = cross_prod(a4c,a1c) + !d12=dot_product(p1c,crs12) + !d23=dot_product(p1c,crs23) + !d34=dot_product(p1c,crs34) + !d41=dot_product(p1c,crs41) + !signs(1)= (d12<0.0) + !signs(2)= (d23<0.0) + !signs(3)= (d34<0.0) + !signs(4)= (d41<0.0) + !in_poly=( (count(signs)==0) .or. (count(signs)==4) ) + + !end function point_in_polygon_crossprod + + function lines_intersect(b0,b1,a0,a1) result(intersect) + real(real64), intent(in) :: b0(3),b1(3),a0(3),a1(3) + logical :: intersect + real(real64) :: p(3),q(3),t(3) + real(real64) :: s1,s2,s3,s4 + logical :: signs(4) + + intersect=.false. + q=cross_prod(b0,b1) + p=cross_prod(a0,a1) + t=normal_vect(cross_prod(p,q)) + + s1=dot_product(cross_prod(a0,p),t) + s2=dot_product(cross_prod(a1,p),t) + s3=dot_product(cross_prod(b0,q),t) + s4=dot_product(cross_prod(b1,q),t) + + signs(1) = -s1 <0.d0 + signs(2) = s2 <0.d0 + signs(3) = -s3 < 0.d0 + signs(4) = s4 < 0.d0 + + intersect = ((count(signs)==0) .or. (count(signs)==4)) + + end function lines_intersect + + function normal_vect(vin) result(vout) + real(real64), intent(in) :: vin(3) + real(real64) :: vout(3) + vout=vin/sqrt(vin(1)*vin(1)+vin(2)*vin(2)+vin(3)*vin(3)) + + end function normal_vect + + function cross_prod(v1,v2) result(vout) + real(real64), intent(in) :: v1(3),v2(3) + real(real64) :: vout(3) + vout(1)=v1(2)*v2(3)-v1(3)*v2(2) + vout(2)=v1(3)*v2(1)-v1(1)*v2(3) + vout(3)=v1(1)*v2(2)-v1(2)*v2(1) + end function cross_prod + + function convert_to_cart(v) result(xyz) + real(real64), intent(in) :: v(2) + real(real64) :: xyz(3) + + xyz(1)=cos(v(2))*cos(v(1)) + xyz(2)=cos(v(2))*sin(v(1)) + xyz(3)=sin(v(2)) + + end function convert_to_cart + + function spherical_angles(p1,p2,p3) result(spherical_angle) + real(real64) :: spherical_angle + real(real64), intent(in) :: p1(3),p2(3),p3(3) + + real (real64):: e1(3), e2(3), e3(3) + real (real64):: px, py, pz + real (real64):: qx, qy, qz + real (real64):: angle, ddd + integer n + + do n=1,3 + e1(n) = p1(n) + e2(n) = p2(n) + e3(n) = p3(n) + enddo + + !------------------------------------------------------------------- + ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry + !------------------------------------------------------------------- + ! Vector P: + px = e1(2)*e2(3) - e1(3)*e2(2) + py = e1(3)*e2(1) - e1(1)*e2(3) + pz = e1(1)*e2(2) - e1(2)*e2(1) + ! Vector Q: + qx = e1(2)*e3(3) - e1(3)*e3(2) + qy = e1(3)*e3(1) - e1(1)*e3(3) + qz = e1(1)*e3(2) - e1(2)*e3(1) + + ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz) + + if ( ddd <= 0.0d0 ) then + angle = 0.d0 else - angle = acos( ddd ) + ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd) + if ( abs(ddd)>1.d0) then + angle = 0.5d0 * MAPL_PI_R8 + if (ddd < 0.d0) then + angle = MAPL_PI_R8 + else + angle = 0.d0 + end if + else + angle = acos( ddd ) + endif endif - endif - spherical_angle = angle -end function + spherical_angle = angle + end function spherical_angles + + function vect_mag(v) result(mag) + real(real64), intent(in) :: v(3) + real :: mag + mag = sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) + end function vect_mag + +end module mapl_SphericalGeometry -end module MAPL_SphericalGeometry From 989b82677642fb6b53383bda9ac34f608312f9ce Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Thu, 18 Nov 2021 17:49:50 -0500 Subject: [PATCH 2/9] Fixes to support REAL64 precision. --- base/HorizontalFluxRegridder.F90 | 2 +- base/MAPL_SphericalGeometry.F90 | 26 +++++++++++++++++++++----- 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/base/HorizontalFluxRegridder.F90 b/base/HorizontalFluxRegridder.F90 index 5d355ed03795..59324165ea5e 100644 --- a/base/HorizontalFluxRegridder.F90 +++ b/base/HorizontalFluxRegridder.F90 @@ -74,7 +74,7 @@ subroutine initialize_subclass(this, unusable, rc) integer :: counts(5) integer :: status integer :: units ! unused - real, pointer :: lons(:,:), lats(:,:) + real(kind=ESMF_KIND_R8), pointer :: lons(:,:), lats(:,:) _UNUSED_DUMMY(unusable) spec = this%get_spec() diff --git a/base/MAPL_SphericalGeometry.F90 b/base/MAPL_SphericalGeometry.F90 index f901cf1383df..088e9de36f6c 100644 --- a/base/MAPL_SphericalGeometry.F90 +++ b/base/MAPL_SphericalGeometry.F90 @@ -4,7 +4,7 @@ module mapl_SphericalGeometry use mapl_ErrorHandlingMod use ESMF use mapl_Constants - use, intrinsic :: iso_fortran_env, only: REAL64,REAL32 + use, intrinsic :: iso_fortran_env, only: REAL64, REAL32 implicit none private @@ -12,6 +12,12 @@ module mapl_SphericalGeometry public get_area_spherical_polygon public :: distance + + interface distance + module procedure distance_r32 + module procedure distance_r64 + end interface distance + contains ! get area of spherical rectangle given the four corners @@ -55,15 +61,25 @@ end function get_area_spherical_polygon ! and returns distance in radians (unit sphere) ! Using formulae from: https://www.movable-type.co.uk/scripts/latlong.html - elemental real function distance(lon1, lat1, lon2, lat2) result(d) - real, intent(in) :: lon1, lat1 - real, intent(in) :: lon2, lat2 + elemental real(kind=REAL64) function distance_r64(lon1, lat1, lon2, lat2) result(d) + real(kind=REAL64), intent(in) :: lon1, lat1 + real(kind=REAL64), intent(in) :: lon2, lat2 + + associate(a => sin(lat2-lat1)**2 + cos(lat1)*cos(lat2)*sin((lon2-lon1)/2)**2) + d = 2*atan2(sqrt(a), sqrt(1-a)) + end associate + + end function distance_r64 + + elemental real(kind=REAL32) function distance_r32(lon1, lat1, lon2, lat2) result(d) + real(kind=REAL32), intent(in) :: lon1, lat1 + real(kind=REAL32), intent(in) :: lon2, lat2 associate(a => sin(lat2-lat1)**2 + cos(lat1)*cos(lat2)*sin((lon2-lon1)/2)**2) d = 2*atan2(sqrt(a), sqrt(1-a)) end associate - end function distance + end function distance_r32 subroutine get_points_in_spherical_domain(center_lons,center_lats,corner_lons,corner_lats,lons,lats,ii,jj,rc) real(real64), intent(in) :: center_lats(:,:),center_lons(:,:) From b5c282aa0a4ca97e17667bc99a071415ffae2c02 Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Thu, 2 Dec 2021 09:20:25 -0500 Subject: [PATCH 3/9] Changed default decomposition topology. --- CHANGELOG.md | 4 ++++ base/Base/Base_Base_implementation.F90 | 13 ++++++++++--- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 15ad87898131..f411561ac32b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -746,10 +746,14 @@ new interface with conventional ordering has been introduced. ### Changed + - Changed the way how a writer is chosen. Previously, a writing processor is chosen as long as it is idle. Now, an idle processor is chosen from a node with the most idle processors. - Changed error checking `_ASSERT` to use `__RC__` macro and `_VERIFY` for UserRC - Changed `_ASSERT` with `.and.` conditional to separate `_ASSERT` to improve error message +- Changed default decomposition algorithm in Base/base. + . More optimal for CS case + . Hopefully aligns with common choices for native decomp to reduce need for nontrivial regridding. - Changed usage of MAPL_IO subroutines in CubedSphere and LatLon Grid Factories to open command with newunit clause - Updated `components.yaml` - ESMA_env v3.7.0 (Use MPT 2.25 at NAS on TOSS) diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index ba3da7a2f233..469148e27800 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -833,7 +833,7 @@ module subroutine MAPL_MakeDecomposition(nx, ny, unusable, reduceFactor, rc) type (ESMF_VM) :: vm integer :: pet_count - + integer :: bias character(len=*), parameter :: Iam= __FILE__ // '::MAPL_MakeDecomposition()' integer :: status @@ -843,11 +843,18 @@ module subroutine MAPL_MakeDecomposition(nx, ny, unusable, reduceFactor, rc) _VERIFY(status) call ESMF_VMGet(vm, petCount=pet_count, rc=status) _VERIFY(status) - if (present(reduceFactor)) pet_count=pet_count/reduceFactor + if (present(reduceFactor)) then + pet_count=pet_count/reduceFactor + ! Assume CS + bias = 1 + else + ! Assume Lat-Lon + bias =2 + end if ! count down from sqrt(n) ! Note: inal iteration (nx=1) is guaranteed to succeed. - do nx = floor(sqrt(real(2*pet_count))), 1, -1 + do nx = nint(sqrt(real(bias*pet_count))), 1, -1 if (mod(pet_count, nx) == 0) then ! found a decomposition ny = pet_count / nx exit From 7dc923462e76a4f595f46286ab4b27ba5e6c3588 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Tue, 4 Apr 2023 12:43:57 -0400 Subject: [PATCH 4/9] Fix up changelog --- CHANGELOG.md | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f411561ac32b..5e386480a70a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,10 +11,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +- Changed default decomposition algorithm in Base/base. + - More optimal for CS case + - Hopefully aligns with common choices for native decomp to reduce need for nontrivial regridding. + ### Fixed + - Corrected bug in HorizontalFluxRegridder. Fluxes need to be multiplied by edge length for correct treatment. +### Removed + +### Deprecated ## [2.37.0] - 2023-04-03 @@ -746,14 +754,10 @@ new interface with conventional ordering has been introduced. ### Changed - - Changed the way how a writer is chosen. Previously, a writing processor is chosen as long as it is idle. Now, an idle processor is chosen from a node with the most idle processors. - Changed error checking `_ASSERT` to use `__RC__` macro and `_VERIFY` for UserRC - Changed `_ASSERT` with `.and.` conditional to separate `_ASSERT` to improve error message -- Changed default decomposition algorithm in Base/base. - . More optimal for CS case - . Hopefully aligns with common choices for native decomp to reduce need for nontrivial regridding. - Changed usage of MAPL_IO subroutines in CubedSphere and LatLon Grid Factories to open command with newunit clause - Updated `components.yaml` - ESMA_env v3.7.0 (Use MPT 2.25 at NAS on TOSS) From d2b0e12a2984eba1b7faeddca455b82aad36df05 Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Mon, 17 Jul 2023 16:44:59 -0400 Subject: [PATCH 5/9] Update HorizontalFluxRegridder.F90 Fail harder if hflux does not meet requirements. --- base/HorizontalFluxRegridder.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/base/HorizontalFluxRegridder.F90 b/base/HorizontalFluxRegridder.F90 index 59324165ea5e..0db9a1b839f3 100644 --- a/base/HorizontalFluxRegridder.F90 +++ b/base/HorizontalFluxRegridder.F90 @@ -57,6 +57,7 @@ logical function supports(spec, unusable, rc) call MAPL_GridGet(spec%grid_out, localCellCountPerDim=counts_out, _RC) supports = all(mod(counts_in(1:2), counts_out(1:2)) == 0) .or. all(mod(counts_out, counts_in) == 0) + _ASSERT(supports, "HFlux regridder requires local domains to be properly nested.") _RETURN(_SUCCESS) end function supports From 80473d506ed0647ca1b57d89f91ad1c3050aaa84 Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Tue, 18 Jul 2023 12:35:08 -0400 Subject: [PATCH 6/9] Fixed bad merge. --- base/MAPL_SphericalGeometry.F90 | 594 ++++++++++++++++---------------- 1 file changed, 297 insertions(+), 297 deletions(-) diff --git a/base/MAPL_SphericalGeometry.F90 b/base/MAPL_SphericalGeometry.F90 index 088e9de36f6c..0f95a131ae9b 100644 --- a/base/MAPL_SphericalGeometry.F90 +++ b/base/MAPL_SphericalGeometry.F90 @@ -1,13 +1,13 @@ #include "MAPL_ErrLog.h" -module mapl_SphericalGeometry - use mapl_KeywordEnforcerMod +module MAPL_SphericalGeometry + use MAPL_KeywordEnforcerMod use mapl_ErrorHandlingMod use ESMF - use mapl_Constants - use, intrinsic :: iso_fortran_env, only: REAL64, REAL32 + use MAPL_Constants + use, intrinsic :: iso_fortran_env, only: REAL64,REAL32 + implicit none private - public get_points_in_spherical_domain public get_area_spherical_polygon public :: distance @@ -18,44 +18,300 @@ module mapl_SphericalGeometry module procedure distance_r64 end interface distance -contains - - ! get area of spherical rectangle given the four corners - ! p4 ------ p3 - ! | | - ! | | - ! | | - ! p1 ------ p2 - function get_area_spherical_polygon(p1,p4,p2,p3) result(area) - real(real64) :: area - real(real64), intent(in) :: p1(2),p2(2),p3(2),p4(2) - - real(real64) :: e1(3),e2(3),e3(3) - real(real64) :: ang1,ang2,ang3,ang4 - - e1 = convert_to_cart(p1) - e2 = convert_to_cart(p2) - e3 = convert_to_cart(p4) - ang1 = spherical_angles(e1, e2, e3) - e1 = convert_to_cart(p2) - e2 = convert_to_cart(p3) - e3 = convert_to_cart(p1) - ang2 = spherical_angles(e1, e2, e3) - - e1 = convert_to_cart(p3) - e2 = convert_to_cart(p4) - e3 = convert_to_cart(p2) - ang3 = spherical_angles(e1, e2, e3) - - e1 = convert_to_cart(p4) - e2 = convert_to_cart(p3) - e3 = convert_to_cart(p1) - ang4 = spherical_angles(e1, e2, e3) +contains - area = ang1 + ang2 + ang3 + ang4 - 2.0d0*MAPL_PI_R8 + ! get area of spherical rectangle given the four corners + ! p4 ------ p3 + ! | | + ! | | + ! | | + ! p1 ------ p2 + function get_area_spherical_polygon(p1,p4,p2,p3) result(area) + real(real64) :: area + real(real64), intent(in) :: p1(2),p2(2),p3(2),p4(2) + + real(real64) :: e1(3),e2(3),e3(3) + real(real64) :: ang1,ang2,ang3,ang4 + + e1 = convert_to_cart(p1) + e2 = convert_to_cart(p2) + e3 = convert_to_cart(p4) + ang1 = spherical_angles(e1, e2, e3) + + e1 = convert_to_cart(p2) + e2 = convert_to_cart(p3) + e3 = convert_to_cart(p1) + ang2 = spherical_angles(e1, e2, e3) + + e1 = convert_to_cart(p3) + e2 = convert_to_cart(p4) + e3 = convert_to_cart(p2) + ang3 = spherical_angles(e1, e2, e3) + + e1 = convert_to_cart(p4) + e2 = convert_to_cart(p3) + e3 = convert_to_cart(p1) + ang4 = spherical_angles(e1, e2, e3) + + area = ang1 + ang2 + ang3 + ang4 - 2.0d0*MAPL_PI_R8 + + end function get_area_spherical_polygon + + subroutine get_points_in_spherical_domain(center_lons,center_lats,corner_lons,corner_lats,lons,lats,ii,jj,rc) + real(real64), intent(in) :: center_lats(:,:),center_lons(:,:) + real(real64), intent(in) :: corner_lats(:,:),corner_lons(:,:) + real(real64), intent(in) :: lons(:),lats(:) + integer, intent(out) :: ii(:),jj(:) + integer, intent(out), optional :: rc + + integer :: npts,i,n,niter,im,jm,ilb,jlb,iub,jub,ifound,jfound + integer :: lold,uold,lnew,unew + logical :: in_region,in_sub_region + + npts = size(lats) + + _ASSERT(size(lats)==size(lons),"lats and lons do not match") + _ASSERT(npts==size(ii),"size of ii does not match") + _ASSERT(npts==size(ii),"size of jj does not match") + + im=size(corner_lons,1)-1 + jm=size(corner_lons,2)-1 + niter = max(im,jm) + + do i=1,npts + ifound=-1 + jfound=-1 + ilb=1 + iub=im + jlb=1 + jub=jm + in_region = point_in_polygon([lons(i),lats(i)],[center_lons(ilb,jlb),center_lats(ilb,jlb)], & + [corner_lons(ilb,jlb),corner_lats(ilb,jlb)], & + [corner_lons(iub+1,jlb),corner_lats(iub+1,jlb)], & + [corner_lons(iub+1,jub+1),corner_lats(iub+1,jub+1)], & + [corner_lons(ilb,jub+1),corner_lats(ilb,jub+1)]) + if (in_region) then + ! bisect first dimension + lnew=ilb + unew=iub + do n = 1,niter + lold=lnew + uold=unew + unew=lold+(uold-lold)/2 + in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(lnew,jlb),center_lats(lnew,jlb)], & + [corner_lons(lnew,jlb),corner_lats(lnew,jlb)], & + [corner_lons(unew+1,jlb),corner_lats(unew+1,jlb)], & + [corner_lons(unew+1,jub+1),corner_lats(unew+1,jub+1)], & + [corner_lons(lnew,jub+1),corner_lats(lnew,jub+1)]) + if (in_sub_region) then + lnew=lold + unew=unew + else + lnew=unew+1 + unew=uold + end if + if (unew==lnew) then + ifound=unew + exit + end if + enddo + ! bisect 2nd dimension + lnew=jlb + unew=jub + do n = 1,niter + lold=lnew + uold=unew + unew=lold+(uold-lold)/2 + in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(ifound,lnew),center_lats(ifound,lnew)] , & + [corner_lons(ifound,lnew),corner_lats(ifound,lnew)], & + [corner_lons(ifound+1,lnew),corner_lats(ifound+1,lnew)], & + [corner_lons(ifound+1,unew+1),corner_lats(ifound+1,unew+1)], & + [corner_lons(ifound,unew+1),corner_lats(ifound,unew+1)]) + if (in_sub_region) then + lnew=lold + unew=unew + else + lnew=unew+1 + unew=uold + end if + if (unew==lnew) then + jfound=unew + exit + end if + enddo + end if + ii(i)=ifound + jj(i)=jfound + enddo + _RETURN(_SUCCESS) + + end subroutine get_points_in_spherical_domain + + function point_in_polygon(p0,pinside,a1,a2,a3,a4) result(in_poly) + real(real64), intent(in) :: p0(2),pinside(2),a1(2),a2(2),a3(2),a4(2) + logical :: in_poly + + real(real64) :: p1c(3),p2c(3),a1c(3),a2c(3),a3c(3),a4c(3) + logical :: intersect(4) + p1c=convert_to_cart(p0) + p2c=convert_to_cart(pinside) + a1c=convert_to_cart(a1) + a2c=convert_to_cart(a2) + a3c=convert_to_cart(a3) + a4c=convert_to_cart(a4) + + intersect(1) = lines_intersect(p1c,p2c,a1c,a2c) + intersect(2) = lines_intersect(p1c,p2c,a2c,a3c) + intersect(3) = lines_intersect(p1c,p2c,a3c,a4c) + intersect(4) = lines_intersect(p1c,p2c,a4c,a1c) + if (mod(count(intersect),2)==0) then + in_poly=.true. + else + in_poly=.false. + end if + + + end function point_in_polygon + +! it is claimed this should work but doesn't + !function point_in_polygon_crosprod(p1,a1,a2,a3,a4) result(in_poly) + !real(real64), intent(in) :: p1(2),a1(2),a2(2),a3(2),a4(2) + !logical :: in_poly + + !real(real64) :: p1c(3),a1c(3),a2c(3),a3c(3),a4c(3) + !real(real64) :: crs12(3),crs23(3),crs34(3),crs41(3) + !real(real64) :: d12,d23,d34,d41 + !logical :: signs(4) + ! a1 -> a2 -> a3 -> a4 so a4 connects to a1 + + !p1c=convert_to_cart(p1) + !a1c=convert_to_cart(a1) + !a2c=convert_to_cart(a2) + !a3c=convert_to_cart(a3) + !a4c=convert_to_cart(a4) + + !crs12 = cross_prod(a1c,a2c) + !crs23 = cross_prod(a2c,a3c) + !crs34 = cross_prod(a3c,a4c) + !crs41 = cross_prod(a4c,a1c) + !d12=dot_product(p1c,crs12) + !d23=dot_product(p1c,crs23) + !d34=dot_product(p1c,crs34) + !d41=dot_product(p1c,crs41) + !signs(1)= (d12<0.0) + !signs(2)= (d23<0.0) + !signs(3)= (d34<0.0) + !signs(4)= (d41<0.0) + !in_poly=( (count(signs)==0) .or. (count(signs)==4) ) + + !end function point_in_polygon_crossprod + + function lines_intersect(b0,b1,a0,a1) result(intersect) + real(real64), intent(in) :: b0(3),b1(3),a0(3),a1(3) + logical :: intersect + real(real64) :: p(3),q(3),t(3) + real(real64) :: s1,s2,s3,s4 + logical :: signs(4) + + intersect=.false. + q=cross_prod(b0,b1) + p=cross_prod(a0,a1) + t=normal_vect(cross_prod(p,q)) + + s1=dot_product(cross_prod(a0,p),t) + s2=dot_product(cross_prod(a1,p),t) + s3=dot_product(cross_prod(b0,q),t) + s4=dot_product(cross_prod(b1,q),t) + + signs(1) = -s1 <0.d0 + signs(2) = s2 <0.d0 + signs(3) = -s3 < 0.d0 + signs(4) = s4 < 0.d0 + + intersect = ((count(signs)==0) .or. (count(signs)==4)) + + end function lines_intersect + + function normal_vect(vin) result(vout) + real(real64), intent(in) :: vin(3) + real(real64) :: vout(3) + vout=vin/sqrt(vin(1)*vin(1)+vin(2)*vin(2)+vin(3)*vin(3)) + + end function normal_vect + + function cross_prod(v1,v2) result(vout) + real(real64), intent(in) :: v1(3),v2(3) + real(real64) :: vout(3) + vout(1)=v1(2)*v2(3)-v1(3)*v2(2) + vout(2)=v1(3)*v2(1)-v1(1)*v2(3) + vout(3)=v1(1)*v2(2)-v1(2)*v2(1) + end function cross_prod + + function convert_to_cart(v) result(xyz) + real(real64), intent(in) :: v(2) + real(real64) :: xyz(3) + + xyz(1)=cos(v(2))*cos(v(1)) + xyz(2)=cos(v(2))*sin(v(1)) + xyz(3)=sin(v(2)) + + end function convert_to_cart + +function vect_mag(v) result(mag) + real(real64), intent(in) :: v(3) + real :: mag + mag = sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) +end function vect_mag + +function spherical_angles(p1,p2,p3) result(spherical_angle) + real(real64) :: spherical_angle + real(real64), intent(in) :: p1(3),p2(3),p3(3) + + real (real64):: e1(3), e2(3), e3(3) + real (real64):: px, py, pz + real (real64):: qx, qy, qz + real (real64):: angle, ddd + integer n + + do n=1,3 + e1(n) = p1(n) + e2(n) = p2(n) + e3(n) = p3(n) + enddo + + !------------------------------------------------------------------- + ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry + !------------------------------------------------------------------- + ! Vector P: + px = e1(2)*e2(3) - e1(3)*e2(2) + py = e1(3)*e2(1) - e1(1)*e2(3) + pz = e1(1)*e2(2) - e1(2)*e2(1) + ! Vector Q: + qx = e1(2)*e3(3) - e1(3)*e3(2) + qy = e1(3)*e3(1) - e1(1)*e3(3) + qz = e1(1)*e3(2) - e1(2)*e3(1) + + ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz) + + if ( ddd <= 0.0d0 ) then + angle = 0.d0 + else + ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd) + if ( abs(ddd)>1.d0) then + angle = 0.5d0 * MAPL_PI_R8 + if (ddd < 0.d0) then + angle = MAPL_PI_R8 + else + angle = 0.d0 + end if + else + angle = acos( ddd ) + endif + endif - end function get_area_spherical_polygon + spherical_angle = angle +end function ! Computes distance between two points (in lat lon as radians), ! and returns distance in radians (unit sphere) @@ -81,260 +337,4 @@ elemental real(kind=REAL32) function distance_r32(lon1, lat1, lon2, lat2) result end function distance_r32 - subroutine get_points_in_spherical_domain(center_lons,center_lats,corner_lons,corner_lats,lons,lats,ii,jj,rc) - real(real64), intent(in) :: center_lats(:,:),center_lons(:,:) - real(real64), intent(in) :: corner_lats(:,:),corner_lons(:,:) - real(real64), intent(in) :: lons(:),lats(:) - integer, intent(out) :: ii(:),jj(:) - integer, intent(out), optional :: rc - - integer :: npts,i,n,niter,im,jm,ilb,jlb,iub,jub,ifound,jfound - integer :: lold,uold,lnew,unew - logical :: in_region,in_sub_region - - npts = size(lats) - - _ASSERT(size(lats)==size(lons),"lats and lons do not match") - _ASSERT(npts==size(ii),"size of ii does not match") - _ASSERT(npts==size(ii),"size of jj does not match") - - im=size(corner_lons,1)-1 - jm=size(corner_lons,2)-1 - niter = max(im,jm) - - do i=1,npts - ifound=-1 - jfound=-1 - ilb=1 - iub=im - jlb=1 - jub=jm - in_region = point_in_polygon([lons(i),lats(i)],[center_lons(ilb,jlb),center_lats(ilb,jlb)], & - [corner_lons(ilb,jlb),corner_lats(ilb,jlb)], & - [corner_lons(iub+1,jlb),corner_lats(iub+1,jlb)], & - [corner_lons(iub+1,jub+1),corner_lats(iub+1,jub+1)], & - [corner_lons(ilb,jub+1),corner_lats(ilb,jub+1)]) - if (in_region) then - ! bisect first dimension - lnew=ilb - unew=iub - do n = 1,niter - lold=lnew - uold=unew - unew=lold+(uold-lold)/2 - in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(lnew,jlb),center_lats(lnew,jlb)], & - [corner_lons(lnew,jlb),corner_lats(lnew,jlb)], & - [corner_lons(unew+1,jlb),corner_lats(unew+1,jlb)], & - [corner_lons(unew+1,jub+1),corner_lats(unew+1,jub+1)], & - [corner_lons(lnew,jub+1),corner_lats(lnew,jub+1)]) - if (in_sub_region) then - lnew=lold - unew=unew - else - lnew=unew+1 - unew=uold - end if - if (unew==lnew) then - ifound=unew - exit - end if - enddo - ! bisect 2nd dimension - lnew=jlb - unew=jub - do n = 1,niter - lold=lnew - uold=unew - unew=lold+(uold-lold)/2 - in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(ifound,lnew),center_lats(ifound,lnew)] , & - [corner_lons(ifound,lnew),corner_lats(ifound,lnew)], & - [corner_lons(ifound+1,lnew),corner_lats(ifound+1,lnew)], & - [corner_lons(ifound+1,unew+1),corner_lats(ifound+1,unew+1)], & - [corner_lons(ifound,unew+1),corner_lats(ifound,unew+1)]) - if (in_sub_region) then - lnew=lold - unew=unew - else - lnew=unew+1 - unew=uold - end if - if (unew==lnew) then - jfound=unew - exit - end if - enddo - end if - ii(i)=ifound - jj(i)=jfound - enddo - _RETURN(_SUCCESS) - - end subroutine get_points_in_spherical_domain - - function point_in_polygon(p0,pinside,a1,a2,a3,a4) result(in_poly) - real(real64), intent(in) :: p0(2),pinside(2),a1(2),a2(2),a3(2),a4(2) - logical :: in_poly - - real(real64) :: p1c(3),p2c(3),a1c(3),a2c(3),a3c(3),a4c(3) - logical :: intersect(4) - p1c=convert_to_cart(p0) - p2c=convert_to_cart(pinside) - a1c=convert_to_cart(a1) - a2c=convert_to_cart(a2) - a3c=convert_to_cart(a3) - a4c=convert_to_cart(a4) - - intersect(1) = lines_intersect(p1c,p2c,a1c,a2c) - intersect(2) = lines_intersect(p1c,p2c,a2c,a3c) - intersect(3) = lines_intersect(p1c,p2c,a3c,a4c) - intersect(4) = lines_intersect(p1c,p2c,a4c,a1c) - if (mod(count(intersect),2)==0) then - in_poly=.true. - else - in_poly=.false. - end if - - - end function point_in_polygon - - ! it is claimed this should work but doesn't - !function point_in_polygon_crosprod(p1,a1,a2,a3,a4) result(in_poly) - !real(real64), intent(in) :: p1(2),a1(2),a2(2),a3(2),a4(2) - !logical :: in_poly - - !real(real64) :: p1c(3),a1c(3),a2c(3),a3c(3),a4c(3) - !real(real64) :: crs12(3),crs23(3),crs34(3),crs41(3) - !real(real64) :: d12,d23,d34,d41 - !logical :: signs(4) - !! a1 -> a2 -> a3 -> a4 so a4 connects to a1 - - !p1c=convert_to_cart(p1) - !a1c=convert_to_cart(a1) - !a2c=convert_to_cart(a2) - !a3c=convert_to_cart(a3) - !a4c=convert_to_cart(a4) - - !crs12 = cross_prod(a1c,a2c) - !crs23 = cross_prod(a2c,a3c) - !crs34 = cross_prod(a3c,a4c) - !crs41 = cross_prod(a4c,a1c) - !d12=dot_product(p1c,crs12) - !d23=dot_product(p1c,crs23) - !d34=dot_product(p1c,crs34) - !d41=dot_product(p1c,crs41) - !signs(1)= (d12<0.0) - !signs(2)= (d23<0.0) - !signs(3)= (d34<0.0) - !signs(4)= (d41<0.0) - !in_poly=( (count(signs)==0) .or. (count(signs)==4) ) - - !end function point_in_polygon_crossprod - - function lines_intersect(b0,b1,a0,a1) result(intersect) - real(real64), intent(in) :: b0(3),b1(3),a0(3),a1(3) - logical :: intersect - real(real64) :: p(3),q(3),t(3) - real(real64) :: s1,s2,s3,s4 - logical :: signs(4) - - intersect=.false. - q=cross_prod(b0,b1) - p=cross_prod(a0,a1) - t=normal_vect(cross_prod(p,q)) - - s1=dot_product(cross_prod(a0,p),t) - s2=dot_product(cross_prod(a1,p),t) - s3=dot_product(cross_prod(b0,q),t) - s4=dot_product(cross_prod(b1,q),t) - - signs(1) = -s1 <0.d0 - signs(2) = s2 <0.d0 - signs(3) = -s3 < 0.d0 - signs(4) = s4 < 0.d0 - - intersect = ((count(signs)==0) .or. (count(signs)==4)) - - end function lines_intersect - - function normal_vect(vin) result(vout) - real(real64), intent(in) :: vin(3) - real(real64) :: vout(3) - vout=vin/sqrt(vin(1)*vin(1)+vin(2)*vin(2)+vin(3)*vin(3)) - - end function normal_vect - - function cross_prod(v1,v2) result(vout) - real(real64), intent(in) :: v1(3),v2(3) - real(real64) :: vout(3) - vout(1)=v1(2)*v2(3)-v1(3)*v2(2) - vout(2)=v1(3)*v2(1)-v1(1)*v2(3) - vout(3)=v1(1)*v2(2)-v1(2)*v2(1) - end function cross_prod - - function convert_to_cart(v) result(xyz) - real(real64), intent(in) :: v(2) - real(real64) :: xyz(3) - - xyz(1)=cos(v(2))*cos(v(1)) - xyz(2)=cos(v(2))*sin(v(1)) - xyz(3)=sin(v(2)) - - end function convert_to_cart - - function spherical_angles(p1,p2,p3) result(spherical_angle) - real(real64) :: spherical_angle - real(real64), intent(in) :: p1(3),p2(3),p3(3) - - real (real64):: e1(3), e2(3), e3(3) - real (real64):: px, py, pz - real (real64):: qx, qy, qz - real (real64):: angle, ddd - integer n - - do n=1,3 - e1(n) = p1(n) - e2(n) = p2(n) - e3(n) = p3(n) - enddo - - !------------------------------------------------------------------- - ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry - !------------------------------------------------------------------- - ! Vector P: - px = e1(2)*e2(3) - e1(3)*e2(2) - py = e1(3)*e2(1) - e1(1)*e2(3) - pz = e1(1)*e2(2) - e1(2)*e2(1) - ! Vector Q: - qx = e1(2)*e3(3) - e1(3)*e3(2) - qy = e1(3)*e3(1) - e1(1)*e3(3) - qz = e1(1)*e3(2) - e1(2)*e3(1) - - ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz) - - if ( ddd <= 0.0d0 ) then - angle = 0.d0 - else - ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd) - if ( abs(ddd)>1.d0) then - angle = 0.5d0 * MAPL_PI_R8 - if (ddd < 0.d0) then - angle = MAPL_PI_R8 - else - angle = 0.d0 - end if - else - angle = acos( ddd ) - endif - endif - - spherical_angle = angle - end function spherical_angles - - function vect_mag(v) result(mag) - real(real64), intent(in) :: v(3) - real :: mag - mag = sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) - end function vect_mag - -end module mapl_SphericalGeometry - +end module MAPL_SphericalGeometry From 0c281e3dae73b51206dd038ac17f70cca25ead78 Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Mon, 24 Jul 2023 08:37:48 -0400 Subject: [PATCH 7/9] Corners require 1 more pts than centers for a grid. --- base/HorizontalFluxRegridder.F90 | 52 ++++++++++++++++---------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/base/HorizontalFluxRegridder.F90 b/base/HorizontalFluxRegridder.F90 index 0db9a1b839f3..fe821725f9ed 100644 --- a/base/HorizontalFluxRegridder.F90 +++ b/base/HorizontalFluxRegridder.F90 @@ -75,7 +75,7 @@ subroutine initialize_subclass(this, unusable, rc) integer :: counts(5) integer :: status integer :: units ! unused - real(kind=ESMF_KIND_R8), pointer :: lons(:,:), lats(:,:) + real(kind=ESMF_KIND_R8), allocatable :: corner_lonsn(:,:), corner_lats(:,:) _UNUSED_DUMMY(unusable) spec = this%get_spec() @@ -98,31 +98,31 @@ subroutine initialize_subclass(this, unusable, rc) this%resolution_ratio = (IM_in / IM_out) - call ESMF_GridGetCoord(grid_in, coordDim=1, farrayPtr=lons, & - localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__) - call ESMF_GridGetCoord(grid_in, coordDim=2, farrayPtr=lats, & - localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__) - - this%dx_in = distance( & - lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & - lons(2:IM_in+1,1:JM_in), lats(2:IM_in+1,1:JM_in)) - - this%dy_in = distance( & - lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & - lons(1:IM_in,2:JM_in+1), lats(1:IM_in,2:JM_in+1)) - - call ESMF_GridGetCoord(grid_out, coordDim=1, farrayPtr=lons, & - localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__) - call ESMF_GridGetCoord(grid_out, coordDim=2, farrayPtr=lats, & - localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__) - - this%dx_out = distance( & - lons(1:IM_out,1:JM_out), lats(1:IM_out,1:JM_out), & - lons(2:IM_out+1,1:JM_out), lats(2:IM_out+1,1:JM_out)) - - this%dy_out = distance( & - lons(1:IM_out,1:JM_out), lats(1:IM_out,1:JM_out), & - lons(1:IM_out,2:JM_out+1), lats(1:IM_out,2:JM_out+1)) + allocate(corner_lons(IM_in+1,JM_in+1), corner_lats(IM_in+1,JM_in+1)) + associate(lons => corner_lons, lats => corner_lats) + call MAPL_GridGetCorners(gridCornerLons=lons, gridCornerLats=lats, _RC) + + this%dx_in = distance( & + lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & + lons(2:IM_in+1,1:JM_in), lats(2:IM_in+1,1:JM_in)) + + this%dy_in = distance( & + lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & + lons(1:IM_in,2:JM_in+1), lats(1:IM_in,2:JM_in+1)) + end associate + + allocate(corner_lons(IM_out+1,JM_out+1), corner_lats(IM_out+1,JM_out+1)) + associate(lons => corner_lons, lats => corner_lats) + call MAPL_GridGetCorners(gridCornerLons=lons, gridCornerLats=lats, _RC) + + this%dx_out = distance( & + lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & + lons(2:IM_in+1,1:JM_in), lats(2:IM_in+1,1:JM_in)) + + this%dy_out = distance( & + lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & + lons(1:IM_in,2:JM_in+1), lats(1:IM_in,2:JM_in+1)) + end associate end associate end associate From 4925e405fbc33867cda8ba48ec3ff52424350132 Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Mon, 24 Jul 2023 08:59:18 -0400 Subject: [PATCH 8/9] Minor corrections. --- base/HorizontalFluxRegridder.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/base/HorizontalFluxRegridder.F90 b/base/HorizontalFluxRegridder.F90 index fe821725f9ed..cfc0024d2808 100644 --- a/base/HorizontalFluxRegridder.F90 +++ b/base/HorizontalFluxRegridder.F90 @@ -9,7 +9,8 @@ module mapl_HorizontalFluxRegridder use mapl_RegridMethods use mapl_KeywordEnforcerMod use mapl_ErrorHandlingMod - use mapl_BaseMod + use mapl_MaplGrid + use mapl_Base use mapl_SphericalGeometry implicit none private @@ -75,7 +76,7 @@ subroutine initialize_subclass(this, unusable, rc) integer :: counts(5) integer :: status integer :: units ! unused - real(kind=ESMF_KIND_R8), allocatable :: corner_lonsn(:,:), corner_lats(:,:) + real(kind=ESMF_KIND_R8), allocatable :: corner_lons(:,:), corner_lats(:,:) _UNUSED_DUMMY(unusable) spec = this%get_spec() @@ -100,7 +101,7 @@ subroutine initialize_subclass(this, unusable, rc) allocate(corner_lons(IM_in+1,JM_in+1), corner_lats(IM_in+1,JM_in+1)) associate(lons => corner_lons, lats => corner_lats) - call MAPL_GridGetCorners(gridCornerLons=lons, gridCornerLats=lats, _RC) + call MAPL_GridGetCorners(grid_in, gridCornerLons=lons, gridCornerLats=lats, _RC) this%dx_in = distance( & lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & @@ -111,9 +112,10 @@ subroutine initialize_subclass(this, unusable, rc) lons(1:IM_in,2:JM_in+1), lats(1:IM_in,2:JM_in+1)) end associate + deallocate(corner_lons, corner_lats) allocate(corner_lons(IM_out+1,JM_out+1), corner_lats(IM_out+1,JM_out+1)) associate(lons => corner_lons, lats => corner_lats) - call MAPL_GridGetCorners(gridCornerLons=lons, gridCornerLats=lats, _RC) + call MAPL_GridGetCorners(grid_out, gridCornerLons=lons, gridCornerLats=lats, _RC) this%dx_out = distance( & lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), & From 66ba65bddca03110f440a85a37b2e62ed2e53088 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Sun, 20 Oct 2024 12:55:28 -0400 Subject: [PATCH 9/9] Fix changelog --- CHANGELOG.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 89f84ce42e44..443bb0054381 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added + - Allow update offsets of ±timestep in ExtData2G ### Changed @@ -17,6 +18,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed - Fixed issue of some Baselibs builds appearing to support zstandard. This is not possible due to Baselibs building HDF5 and netCDF as static libraries +- Corrected bug in HorizontalFluxRegridder. Fluxes need to be multiplied by edge length for correct treatment. ### Removed @@ -594,8 +596,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Created cubed-sphere grid factory with files split by face - Removed unneeded and confusing default in History Grid Comp (see #2081) - Fixes in CMake for fArgParse transition -- Corrected bug in HorizontalFluxRegridder. Fluxes need to be - multiplied by edge length for correct treatment. ### Deprecated