diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index 7d43b5ee83..ba0ac56c8d 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -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 @@ -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