From c9c7751e52da1e415bdd7e89b54ee3000b254659 Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Thu, 18 Nov 2021 16:28:00 -0500 Subject: [PATCH 1/3] Fixes #1199 --- CHANGELOG.md | 3 + base/HorizontalFluxRegridder.F90 | 66 ++++- base/MAPL_SphericalGeometry.F90 | 442 ++++++++++++++++--------------- 3 files changed, 290 insertions(+), 221 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 67a60a73330f..68fd30d7d6d8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Corrected bug in HorizontalFluxRegridder. Fluxes need to be + multiplied by edge length for correct treatment. + ### Added - Added ability to generate monthly checkpoints (fixes issue #1065) diff --git a/base/HorizontalFluxRegridder.F90 b/base/HorizontalFluxRegridder.F90 index b73b2bebdca4..712b2bfd626f 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 8cadc85c3d18..ff031a841cae 100644 --- a/base/MAPL_SphericalGeometry.F90 +++ b/base/MAPL_SphericalGeometry.F90 @@ -1,220 +1,236 @@ #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 -implicit none -private -public get_points_in_spherical_domain + public :: get_points_in_spherical_domain + public :: distance + contains - 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 - -end module MAPL_SphericalGeometry + ! 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 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 From 9c7d86962d6a97643977d95b05c4486266aa8cb3 Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Thu, 18 Nov 2021 17:49:50 -0500 Subject: [PATCH 2/3] Fixes to support REAL64 precision. --- base/HorizontalFluxRegridder.F90 | 2 +- base/MAPL_SphericalGeometry.F90 | 28 ++++++++++++++++++++++------ 2 files changed, 23 insertions(+), 7 deletions(-) diff --git a/base/HorizontalFluxRegridder.F90 b/base/HorizontalFluxRegridder.F90 index 712b2bfd626f..4a2887b6663f 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 ff031a841cae..fe799b489faf 100644 --- a/base/MAPL_SphericalGeometry.F90 +++ b/base/MAPL_SphericalGeometry.F90 @@ -4,28 +4,44 @@ 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 public :: get_points_in_spherical_domain public :: distance - + + + interface distance + module procedure distance_r32 + module procedure distance_r64 + end interface distance + contains ! 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 + 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 6e38aa0e6b6ae15ca746543b021b74f11815bd34 Mon Sep 17 00:00:00 2001 From: Tom Clune Date: Thu, 2 Dec 2021 09:20:25 -0500 Subject: [PATCH 3/3] Changed default decomposition topology. --- CHANGELOG.md | 3 +++ base/Base/Base_Base_implementation.F90 | 13 ++++++++++--- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 68fd30d7d6d8..8401dc46ae5c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,9 @@ 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. - Changed usage of MAPL_IO subroutines in CubedSphere and LatLon Grid Factories to open command with newunit clause - Updated `components.yaml` - ESMA_env v3.6.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 6d2cb9053397..d43dc716bba5 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -832,7 +832,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 @@ -842,11 +842,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