subroutine check_topography_data(pio_file)
! Module(s) from CAM-SIMA.
use cam_abortutils, only: check_allocate, endrun
use cam_field_read, only: cam_read_field
use cam_logfile, only: debug_output, debugout_debug, debugout_none, debugout_verbose
use dynconst, only: constant_g => gravit
! Module(s) from CESM Share.
use shr_kind_mod, only: kind_r8 => shr_kind_r8
! Module(s) from external libraries.
use pio, only: file_desc_t, pio_file_is_open
type(file_desc_t), pointer, intent(in) :: pio_file
character(*), parameter :: subname = 'dyn_comp::check_topography_data'
integer :: ierr
logical :: success
real(kind_r8), parameter :: error_tolerance = 1.0E-3_kind_r8 ! Error tolerance for consistency check.
real(kind_r8), allocatable :: surface_geometric_height(:) ! Computed from topography file.
real(kind_r8), allocatable :: surface_geopotential(:) ! Read from topography file.
real(kind_dyn_mpas), pointer :: zgrid(:, :) ! From MPAS. Geometric height (m) at layer interfaces.
call dyn_debug_print(debugout_debug, subname // ' entered')
nullify(zgrid)
call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid')
if (associated(pio_file)) then
call dyn_debug_print(debugout_verbose, 'Topography file is used for consistency check')
if (.not. pio_file_is_open(pio_file)) then
call endrun('Invalid PIO file descriptor', subname, __LINE__)
end if
allocate(surface_geopotential(ncells_solve), stat=ierr)
call check_allocate(ierr, subname, 'surface_geopotential(ncells_solve)', 'dyn_comp', __LINE__)
allocate(surface_geometric_height(ncells_solve), stat=ierr)
call check_allocate(ierr, subname, 'surface_geometric_height(ncells_solve)', 'dyn_comp', __LINE__)
surface_geopotential(:) = 0.0_kind_r8
surface_geometric_height(:) = 0.0_kind_r8
call cam_read_field('PHIS', pio_file, surface_geopotential, success, &
gridname='cam_cell', timelevel=1, log_output=(debug_output > debugout_none))
if (.not. success) then
call endrun('Failed to find variable "PHIS"', subname, __LINE__)
end if
surface_geometric_height(:) = surface_geopotential(:) / constant_g
! Surface geometric height in MPAS should match the values in topography file.
if (any(abs(real(zgrid(1, 1:ncells_solve), kind_r8) - surface_geometric_height(:)) > error_tolerance)) then
call endrun('Surface geometric height in MPAS is not consistent with topography data', subname, __LINE__)
end if
deallocate(surface_geopotential)
deallocate(surface_geometric_height)
else
call dyn_debug_print(debugout_verbose, 'Topography file is not used for consistency check')
! Surface geometric height in MPAS should be zero.
if (any(abs(real(zgrid(1, 1:ncells_solve), kind_r8)) > error_tolerance)) then
call endrun('Surface geometric height in MPAS is not zero', subname, __LINE__)
end if
end if
nullify(zgrid)
call dyn_debug_print(debugout_debug, subname // ' completed')
end subroutine check_topography_data