fixes #199 Bugfix/tclune/#199 hflux fix 3rd try

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
 

base/Base/Base_Base_implementation.F90

@@ -802,19 +802,26 @@ 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
 
     _UNUSED_DUMMY(unusable)
 
     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,15 +3440,15 @@ 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
      two_pi  = MAPL_PI_R8*2
 
      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)

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
 

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