Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Horizontal flux regridding bug fix #36

Merged
merged 1 commit into from
Jul 10, 2024
Merged

Horizontal flux regridding bug fix #36

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 69 additions & 9 deletions base/HorizontalFluxRegridder.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -70,6 +75,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()
Expand All @@ -91,6 +98,34 @@ 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

Expand Down Expand Up @@ -129,9 +164,16 @@ 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
!ewl old: m_y = m_y + v_in(ii,jj)
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
!ewl old: v_out(i,j) = m_y
end do
end do

Expand All @@ -141,9 +183,15 @@ 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
!ewl old: m_x = m_x + u_in(ii,jj)
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
!ewl old: u_out(i,j) = m_x
end do
end do

Expand Down Expand Up @@ -186,9 +234,15 @@ 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
!ewl old: m_y = m_y + v_in(ii,jj)
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
!ewl old: v_out(i,j) = m_y
end do
end do

Expand All @@ -198,9 +252,15 @@ 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
!ewl old: m_x = m_x + u_in(ii,jj)
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
!ewl old: u_out(i,j) = m_x
end do
end do

Expand Down
31 changes: 31 additions & 0 deletions base/MAPL_SphericalGeometry.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,13 @@ module MAPL_SphericalGeometry
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

subroutine get_points_in_spherical_domain(center_lons,center_lats,corner_lons,corner_lats,lons,lats,ii,jj,rc)
Expand Down Expand Up @@ -217,4 +224,28 @@ function vect_mag(v) result(mag)
mag = sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
end function vect_mag

! 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
Loading