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.

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
Changes from all commits
Commits
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
66 changes: 65 additions & 1 deletion src/core_init_atmosphere/mpas_init_atm_cases.F
Original file line number Diff line number Diff line change
Expand Up @@ -211,7 +211,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 +7148,60 @@ function read_text_array(dminfo, filename, xarray) result(ierr)
end function read_text_array


!-----------------------------------------------------------------------
! routine on_unit_sphere
!
!> \brief If mesh describes a unit sphere, return true
!> \author G. Dylan Dickerson
!> \date 03 Jan 2025
!> \details
!>
!> This routine actually checks that the total mesh surface
!> area is less than twice that of a unit sphere to account for
!> potential numerical errors.
!>
!---------------------------------------------------------------------------
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

! Return value
logical unit_sphere

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

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

unit_sphere = .false.

surfArea = 0.0_RKIND
g_surfArea = 0.0_RKIND

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