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

Skip static interpolation step if input grid isn't a unit sphere #1259

Open
wants to merge 4 commits into
base: develop
Choose a base branch
from
Open

Skip static interpolation step if input grid isn't a unit sphere #1259

Changes from 2 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
b3efd5a
Add init_atm function to check if grid is a unit sphere
gdicker1 Jan 3, 2025
74727af
Skip static_interp step if grid isn't a unit sphere
gdicker1 Jan 3, 2025
4c9a875
fixup! Add init_atm function to check if grid is a unit sphere
gdicker1 Jan 8, 2025
12ad874
fixup! Skip static_interp step if grid isn't a unit sphere
gdicker1 Jan 8, 2025
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
70 changes: 68 additions & 2 deletions src/core_init_atmosphere/mpas_init_atm_cases.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,6 @@ subroutine init_atm_setup_case(domain, stream_manager)
type (domain_type), intent(inout) :: domain
type (MPAS_streamManager_type), intent(inout) :: stream_manager


Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd recommend not changing whitespace here.

integer :: ierr
type (block_type), pointer :: block_ptr

Expand Down Expand Up @@ -211,7 +210,15 @@ subroutine init_atm_setup_case(domain, stream_manager)
end if

if (config_static_interp) then
call init_atm_static(mesh, block_ptr % dimensions, block_ptr % configs)
if ( on_unit_sphere(domain % dminfo, mesh) ) then
call init_atm_static(mesh, block_ptr % dimensions, block_ptr % configs)
else
call mpas_log_write('*************************************************************', messageType=MPAS_LOG_ERR)
call mpas_log_write('Input grid does not describe a unit sphere, static processing', messageType=MPAS_LOG_ERR)
call mpas_log_write('has likely already been applied. Repeating this step can ', messageType=MPAS_LOG_ERR)
call mpas_log_write('cause issues running the model and unrealistic results. ', messageType=MPAS_LOG_ERR)
call mpas_log_write('*************************************************************', messageType=MPAS_LOG_CRIT)
end if
end if

if (config_native_gwd_static) then
Expand Down Expand Up @@ -7140,4 +7147,63 @@ function read_text_array(dminfo, filename, xarray) result(ierr)
end function read_text_array


!-----------------------------------------------------------------------
! routine on_unit_sphere
!
!> \brief True if sphere surface area is close to unit sphere surface area (4*pii)
!> \author G. Dylan Dickerson
!> \date 03 Jan 2025
!> \details
!>
!> This routine will use the global sum of cell areas to
!> approximate the surface area of the sphere. If this surface area
!> is within a factor of 4*pii, the mesh is likely still a unit
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The statement "within a factor of 4pii" could be misinterpreted to mean that the surface area is +/- 4pii of something.

Perhaps it would be better to reframe this description slightly to say that the routine returns true if it determines that the mesh is defined on a unit sphere -- which is the critical point -- and to then follow up with some notes that due to inexactness of arithmetic, the routine actually checks that the total sphere area is not larger than twice the area of a unit sphere.

!> sphere and this returns true. Otherwise, returns false since
!> the mesh is likely not a unit sphere.
!>
!---------------------------------------------------------------------------
function on_unit_sphere(dminfo, mesh) result(unit_sphere)

use mpas_pool_routines, only : mpas_pool_get_array, mpas_pool_get_dimension
use mpas_dmpar, only : mpas_dmpar_sum_real
use mpas_constants, only : pii

implicit none

! Input vars
type (dm_info), intent(in) :: dminfo
type (mpas_pool_type), intent(inout) :: mesh

! Output vars
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd suggest Return value rather than Output vars.

logical unit_sphere

! Local vars
real (kind=RKIND) :: surfArea, g_surfArea
real (kind=RKIND) :: tol
integer :: iCell

real (kind=RKIND), dimension(:), pointer :: areaCell
integer, pointer :: nCellsSolve

unit_sphere = .false.

surfArea = 0.0_RKIND
g_surfArea = 0.0_RKIND
tol = 2.0_RKIND
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could be a Fortran parameter (i.e., a const).


call mpas_pool_get_array(mesh, 'areaCell', areaCell)
call mpas_pool_get_dimension(mesh, 'nCellsSolve', nCellsSolve)

! sum over all owned cells, don't double-count halo cells
do iCell = 1, nCellsSolve
surfArea = surfArea + areaCell(iCell)
end do
call mpas_dmpar_sum_real(dminfo, surfArea, g_surfArea)

if ( (g_surfArea / 4*pii) < tol ) then
unit_sphere = .true.
end if

end function on_unit_sphere

end module init_atm_cases