diff --git a/CHANGELOG.md b/CHANGELOG.md index faa2abd5873e..075de8c477b4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -36,6 +36,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 diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index b50e84e98897..98b675d3c501 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -802,7 +802,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 @@ -810,11 +810,18 @@ module subroutine MAPL_MakeDecomposition(nx, ny, unusable, reduceFactor, rc) call ESMF_VMGetCurrent(vm, _RC) call ESMF_VMGet(vm, petCount=pet_count, _RC) - 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 @@ -2772,8 +2779,8 @@ module subroutine MAPL_GetGlobalHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, _ASSERT( IM_WORLD*6 == JM_WORLD, "It only works for cubed-sphere grid") allocate(lons(npts),lats(npts)) - - call MAPL_Reverse_Schmidt(Grid, stretched, npts, lon=lon, lat=lat, lonR8=lonR8, latR8=latR8, lonRe=lons, latRe=lats, _RC) + + call MAPL_Reverse_Schmidt(Grid, stretched, npts, lon=lon, lat=lat, lonR8=lonR8, latR8=latR8, lonRe=lons, latRe=lats, _RC) dalpha = 2.0d0*alpha/IM_WORLD @@ -2887,7 +2894,7 @@ function grid_is_ok(grid) result(OK) if ( I1 == 1 .and. J1 == 1 ) then allocate(lonRe(j2-j1+1), latRe(j2-j1+1)) call MAPL_Reverse_Schmidt(grid, stretched, J2-J1+1, lonR8=corner_lons(1,:), & - latR8=corner_lats(1,:), lonRe=lonRe, latRe=latRe, _RC) + latR8=corner_lats(1,:), lonRe=lonRe, latRe=latRe, _RC) allocate(accurate_lon(j2-j1+1), accurate_lat(j2-j1+1)) @@ -3382,32 +3389,32 @@ module function MAPL_GetCorrectedPhase(gc,rc) result(phase) end function MAPL_GetCorrectedPhase module subroutine MAPL_Reverse_Schmidt(Grid, stretched, npts, lon, lat, lonR8, latR8, lonRe, latRe, rc) - type(ESMF_Grid), intent(inout) :: Grid + type(ESMF_Grid), intent(inout) :: Grid logical, intent(out) :: stretched integer, intent(in ) :: npts ! number of points in lat and lon arrays real, optional, intent(in ) :: lon(npts) ! array of longitudes in radians real, optional, intent(in ) :: lat(npts) ! array of latitudes in radians real(ESMF_KIND_R8), optional, intent(in ) :: lonR8(npts) ! array of longitudes in radians - real(ESMF_KIND_R8), optional, intent(in ) :: latR8(npts) ! - real(ESMF_KIND_R8), optional, intent(out ) :: lonRe(npts) ! - real(ESMF_KIND_R8), optional, intent(out ) :: latRe(npts) ! + real(ESMF_KIND_R8), optional, intent(in ) :: latR8(npts) ! + real(ESMF_KIND_R8), optional, intent(out ) :: lonRe(npts) ! + real(ESMF_KIND_R8), optional, intent(out ) :: latRe(npts) ! integer, optional, intent(out) :: rc logical :: factorPresent, lonPresent, latPresent integer :: status real(ESMF_KIND_R8) :: c2p1, c2m1, half_pi, two_pi, stretch_factor, target_lon, target_lat - real(ESMF_KIND_R8), dimension(npts) :: x,y,z, Xx, Yy, Zz + real(ESMF_KIND_R8), dimension(npts) :: x,y,z, Xx, Yy, Zz logical, dimension(npts) :: n_s _RETURN_IF( npts == 0 ) - + call ESMF_AttributeGet(grid, name='STRETCH_FACTOR', isPresent= factorPresent, _RC) call ESMF_AttributeGet(grid, name='TARGET_LON', isPresent= lonPresent, _RC) call ESMF_AttributeGet(grid, name='TARGET_LAT', isPresent= latPresent, _RC) if ( factorPresent .and. lonPresent .and. latPresent) then stretched = .true. - else + else stretched = .false. endif @@ -3433,7 +3440,7 @@ module subroutine MAPL_Reverse_Schmidt(Grid, stretched, npts, lon, lat, lonR8, l call ESMF_AttributeGet(grid, name='TARGET_LON', value=target_lon, _RC) call ESMF_AttributeGet(grid, name='TARGET_LAT', value=target_lat, _RC) - c2p1 = 1 + stretch_factor*stretch_factor + c2p1 = 1 + stretch_factor*stretch_factor c2m1 = 1 - stretch_factor*stretch_factor half_pi = MAPL_PI_R8/2 @@ -3441,7 +3448,7 @@ module subroutine MAPL_Reverse_Schmidt(Grid, stretched, npts, lon, lat, lonR8, l target_lon = target_lon*MAPL_DEGREES_TO_RADIANS_R8 target_lat = target_lat*MAPL_DEGREES_TO_RADIANS_R8 - + x = cos(latRe)*cos(lonRe - target_lon) y = cos(latRe)*sin(lonRe - target_lon) z = sin(latRe) diff --git a/base/HorizontalFluxRegridder.F90 b/base/HorizontalFluxRegridder.F90 index 44e2594dc51f..4c7cd572a6b1 100644 --- a/base/HorizontalFluxRegridder.F90 +++ b/base/HorizontalFluxRegridder.F90 @@ -9,7 +9,9 @@ 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 @@ -20,6 +22,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 @@ -54,6 +58,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 @@ -69,6 +74,8 @@ subroutine initialize_subclass(this, unusable, rc) integer :: counts(5) integer :: status + integer :: units ! unused + real(kind=ESMF_KIND_R8), allocatable :: corner_lons(:,:), corner_lats(:,:) _UNUSED_DUMMY(unusable) spec = this%get_spec() @@ -90,8 +97,37 @@ 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) + + 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(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), & + 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 + + 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(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), & + 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 + _RETURN(_SUCCESS) end subroutine initialize_subclass @@ -128,9 +164,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 @@ -140,9 +181,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 @@ -185,9 +230,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 @@ -197,9 +246,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 a6d0d09c22a7..0f95a131ae9b 100644 --- a/base/MAPL_SphericalGeometry.F90 +++ b/base/MAPL_SphericalGeometry.F90 @@ -6,10 +6,19 @@ module MAPL_SphericalGeometry 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 + implicit none + private + public get_points_in_spherical_domain + 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 @@ -304,4 +313,28 @@ function spherical_angles(p1,p2,p3) result(spherical_angle) spherical_angle = angle end function + ! 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(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_r32 + end module MAPL_SphericalGeometry