forked from MPAS-Dev/MPAS-Model
-
Notifications
You must be signed in to change notification settings - Fork 1
mpas_atm_core
David Fillmore edited this page Dec 12, 2025
·
1 revision
! Copyright (c) 2013, Los Alamos National Security, LLC (LANS)
! and the University Corporation for Atmospheric Research (UCAR).
!
! Unless noted otherwise source code is licensed under the BSD license.
! Additional copyright and license information can be found in the LICENSE file
! distributed with this code, or at http://mpas-dev.github.com/license.html
!
module atm_core
use mpas_derived_types
use mpas_pool_routines
use mpas_dmpar
use mpas_log, only : mpas_log_write, mpas_log_info
use mpas_io_units, only : mpas_new_unit, mpas_release_unit
! Provides definition of halo_exchange_routine
#include "mpas_halo_interface.inc"
type (MPAS_Clock_type), pointer :: clock
contains
function atm_core_init(domain, startTimeStamp) result(ierr)
use mpas_timekeeping
use mpas_kind_types
use mpas_stream_manager
use mpas_atm_dimensions, only : mpas_atm_set_dims
use mpas_atm_diagnostics_manager, only : mpas_atm_diag_setup
use mpas_atm_threading, only : mpas_atm_threading_init
use atm_time_integration, only : mpas_atm_dynamics_init
use mpas_timer, only : mpas_timer_start, mpas_timer_stop
use mpas_attlist, only : mpas_modify_att
use mpas_string_utils, only : mpas_string_replace
use mpas_atm_halos, only: atm_build_halo_groups, exchange_halo_group
#ifdef DO_CHEMISTRY
use mpas_atm_chemistry, only: chemistry_init
#endif
implicit none
type (domain_type), intent(inout) :: domain
character(len=*), intent(out) :: startTimeStamp
integer :: ierr
real (kind=RKIND), pointer :: dt
type (block_type), pointer :: block
logical, pointer :: config_do_restart
type (mpas_pool_type), pointer :: state
type (mpas_pool_type), pointer :: mesh
type (mpas_pool_type), pointer :: diag
type (field2DReal), pointer :: u_field, pv_edge_field, ru_field, rw_field
type (field0DReal), pointer :: Time_field
character (len=StrKIND), pointer :: xtime
character (len=StrKIND), pointer :: initial_time1, initial_time2
type (MPAS_Time_Type) :: startTime
real (kind=RKIND), pointer :: nominalMinDc
real (kind=RKIND), pointer :: config_len_disp
real (kind=RKIND), pointer :: Time
integer, pointer :: nVertLevels, maxEdges, maxEdges2, num_scalars
character (len=ShortStrKIND) :: init_stream_name
real (kind=R8KIND) :: input_start_time, input_stop_time
ierr = 0
!
! Setup threading
!
call mpas_atm_threading_init(domain % blocklist, ierr)
if ( ierr /= 0 ) then
call mpas_log_write('Threading setup failed for core '//trim(domain % core % coreName), messageType=MPAS_LOG_CRIT)
end if
!
! Set up inner dimensions used by arrays in optimized dynamics routines
!
call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state)
call mpas_pool_get_dimension(state, 'nVertLevels', nVertLevels)
call mpas_pool_get_dimension(state, 'maxEdges', maxEdges)
call mpas_pool_get_dimension(state, 'maxEdges2', maxEdges2)
call mpas_pool_get_dimension(state, 'num_scalars', num_scalars)
call mpas_atm_set_dims(nVertLevels, maxEdges, maxEdges2, num_scalars)
!
! Set "local" clock to point to the clock contained in the domain type
!
clock => domain % clock
mpas_log_info => domain % logInfo
!
! Build halo exchange groups and set method for exchanging halos in a group
!
call atm_build_halo_groups(domain, ierr)
if (ierr /= 0) then
call mpas_log_write('Failed to build halo exchange groups.', messageType=MPAS_LOG_ERR)
return
end if
call mpas_pool_get_config(domain % blocklist % configs, 'config_do_restart', config_do_restart)
call mpas_pool_get_config(domain % blocklist % configs, 'config_dt', dt)
!
! By default, the 'invariant' stream has an input_interval of "none", so
! the following stream read has no effect. However, in case the 'invariant'
! stream is defined with an input_interval of "initial_only", we read
! time-invariant fields first.
!
call MPAS_stream_mgr_read(domain % streamManager, streamID='invariant', whence=MPAS_STREAM_NEAREST, ierr=ierr)
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='invariant', direction=MPAS_STREAM_INPUT, ierr=ierr)
!
! If this is a restart run, read the restart stream, else read the input
! stream.
! Regardless of which stream we read for initial conditions, reset the
! input alarms for both input and restart before reading any remaining
! input streams.
!
call mpas_timer_start('read_ICs')
if (config_do_restart) then
init_stream_name = 'restart'
else
init_stream_name = 'input'
end if
call mpas_log_write('Reading initial state from '''//trim(init_stream_name)//''' stream')
call mpas_dmpar_get_time(input_start_time)
call MPAS_stream_mgr_read(domain % streamManager, streamID=trim(init_stream_name), ierr=ierr)
call mpas_dmpar_get_time(input_stop_time)
call mpas_timer_stop('read_ICs')
if (ierr /= MPAS_STREAM_MGR_NOERR) then
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
call mpas_log_write('Error reading initial conditions', messageType=MPAS_LOG_ERR)
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
end if
call mpas_log_write(' Timing for read of '''//trim(init_stream_name)//''' stream: $r s', &
realArgs=(/real(input_stop_time - input_start_time, kind=RKIND)/))
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='input', direction=MPAS_STREAM_INPUT, ierr=ierr)
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='restart', direction=MPAS_STREAM_INPUT, ierr=ierr)
call mpas_log_write(' ----- done reading initial state -----')
!
! Determine horizontal length scale used by horizontal diffusion and 3-d divergence damping
!
call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', mesh)
nullify(nominalMinDc)
call mpas_pool_get_array(mesh, 'nominalMinDc', nominalMinDc)
nullify(config_len_disp)
call mpas_pool_get_config(domain % blocklist % configs, 'config_len_disp', config_len_disp)
call mpas_log_write('')
!
! If config_len_disp was specified as a valid value, use that
!
if (config_len_disp > 0.0_RKIND) then
call mpas_log_write('Setting nominalMinDc to $r based on namelist option config_len_disp', realArgs=[config_len_disp])
!
! But if nominalMinDc was available in the input file and is different, print a warning
!
if (nominalMinDc > 0.0_RKIND .and. abs(nominalMinDc - config_len_disp) > 1.0e-6_RKIND * config_len_disp) then
call mpas_log_write('nominalMinDc was read from input file as a positive value ($r) that differs', &
realArgs=[nominalMinDc], messageType=MPAS_LOG_WARN)
call mpas_log_write('from the specified config_len_disp value ($r)', &
realArgs=[config_len_disp], messageType=MPAS_LOG_WARN)
end if
nominalMinDc = config_len_disp
!
! Otherwise, try to use nominalMinDc
!
else
if (nominalMinDc > 0.0_RKIND) then
call mpas_log_write('Setting config_len_disp to $r based on nominalMinDc value in input file', realArgs=[nominalMinDc])
config_len_disp = nominalMinDc
else
call mpas_log_write('Both config_len_disp and nominalMinDc are <= 0.0.', messageType=MPAS_LOG_ERR)
call mpas_log_write('Please either specify config_len_disp in the &nhyd_model namelist group,', &
messageType=MPAS_LOG_ERR)
call mpas_log_write('or use an input file that provides a valid value for the nominalMinDc variable.', &
messageType=MPAS_LOG_ERR)
ierr = 1
return
end if
end if
!
! Read all other inputs
! For now we don't do this here to match results with previous code; to match requires
! that we read in SST and seaice fields after the call to atm_mpas_init_block()
!
! call MPAS_stream_mgr_read(domain % streamManager, ierr=ierr)
! call MPAS_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_INPUT, ierr=ierr)
if (.not. config_do_restart) then
block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'state', state)
call mpas_pool_initialize_time_levels(state)
block => block % next
end do
end if
!
! Read a new data stream for Incremental Analysis Update (IAU), if config_IAU_option /= 'off' : HA (June-15-2016)
! FIXME: should I check xtime for the IAU fields? Maybe not.
! Note: Because the 'iau' stream has the 'iau' package attached to it, the MPAS_stream_mgr_read( )
! call here will actually try to read the stream only if IAU is being used in the run.
!
call MPAS_stream_mgr_read(domain % streamManager, streamID='iau', whence=MPAS_STREAM_NEAREST, ierr=ierr)
if (ierr /= MPAS_STREAM_MGR_NOERR) then
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
call mpas_log_write('Error reading IAU files', messageType=MPAS_LOG_ERR)
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
end if
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='iau', ierr=ierr)
!
! Set startTimeStamp based on the start time of the simulation clock
!
startTime = mpas_get_clock_time(clock, MPAS_START_TIME, ierr)
call mpas_get_time(startTime, dateTimeString=startTimeStamp)
call exchange_halo_group(domain, 'initialization:u')
!
! Perform basic compatibility checks among the fields that were read and the run-time options that were selected
!
call mpas_atm_run_compatibility(domain % dminfo, domain % blocklist, domain % streamManager, ierr)
if (ierr /= 0) then
call mpas_log_write('Please correct issues with the model input fields and/or namelist.', messageType=MPAS_LOG_ERR)
return
end if
block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'mesh', mesh)
call mpas_pool_get_subpool(block % structs, 'state', state)
call atm_mpas_init_block(domain % dminfo, domain % streamManager, block, mesh, dt)
call mpas_pool_get_array(state, 'xtime', xtime, 1)
xtime = startTimeStamp
! Initialize initial_time in second time level. We need to do this because initial state
! is read into time level 1, and if we write output from the set of state arrays that
! represent the original time level 2, the initial_time field will be invalid.
call mpas_pool_get_array(state, 'initial_time', initial_time1, 1)
call mpas_pool_get_array(state, 'initial_time', initial_time2, 2)
initial_time2 = initial_time1
! Set the units to be cf compliant 'seconds since <cf-timestamp>'
call mpas_pool_get_field(state, 'Time', Time_field)
call mpas_modify_att(Time_field % attLists(1) % attList, 'units', &
'seconds since ' // mpas_string_replace(initial_time1, '_', ' '))
block => block % next
end do
call exchange_halo_group(domain, 'initialization:pv_edge,ru,rw')
call mpas_atm_diag_setup(domain % streamManager, domain % blocklist % configs, &
domain % blocklist % structs, domain % clock, domain % dminfo)
!
! Prepare the dynamics for integration
!
call mpas_atm_dynamics_init(domain)
!
! Initialize the chemistry package
!
#ifdef DO_CHEMISTRY
call chemistry_init(domain)
#endif
end function atm_core_init
subroutine atm_simulation_clock_init(core_clock, configs, ierr)
use mpas_timekeeping
implicit none
type (MPAS_Clock_type), intent(inout) :: core_clock
type (mpas_pool_type), intent(inout) :: configs
integer, intent(out) :: ierr
type (MPAS_Time_Type) :: startTime, stopTime
type (MPAS_TimeInterval_type) :: runDuration, timeStep
integer :: local_err
real (kind=RKIND), pointer :: config_dt
character (len=StrKIND), pointer :: config_start_time
character (len=StrKIND), pointer :: config_restart_timestamp_name
character (len=StrKIND), pointer :: config_run_duration
character (len=StrKIND), pointer :: config_stop_time
character (len=StrKIND) :: startTimeStamp
integer :: iounit
ierr = 0
call mpas_pool_get_config(configs, 'config_dt', config_dt)
call mpas_pool_get_config(configs, 'config_start_time', config_start_time)
call mpas_pool_get_config(configs, 'config_restart_timestamp_name', config_restart_timestamp_name)
call mpas_pool_get_config(configs, 'config_run_duration', config_run_duration)
call mpas_pool_get_config(configs, 'config_stop_time', config_stop_time)
if(trim(config_start_time) == 'file') then
call mpas_new_unit(iounit)
open(iounit,file=trim(config_restart_timestamp_name),form='formatted',status='old')
read(iounit,*) startTimeStamp
close(iounit)
call mpas_release_unit(iounit)
else
startTimeStamp = config_start_time
end if
call mpas_set_time(curr_time=startTime, dateTimeString=startTimeStamp, ierr=local_err)
call mpas_set_timeInterval(timeStep, dt=config_dt, ierr=local_err)
if (trim(config_run_duration) /= "none") then
call mpas_set_timeInterval(runDuration, timeString=config_run_duration, ierr=local_err)
call mpas_create_clock(core_clock, startTime=startTime, timeStep=timeStep, runDuration=runDuration, ierr=local_err)
if (trim(config_stop_time) /= "none") then
call mpas_set_time(curr_time=stopTime, dateTimeString=config_stop_time, ierr=local_err)
if(startTime + runduration /= stopTime) then
call mpas_log_write('config_run_duration and config_stop_time are inconsistent: using config_run_duration.', messageType=MPAS_LOG_WARN)
end if
end if
else if (trim(config_stop_time) /= "none") then
call mpas_set_time(curr_time=stopTime, dateTimeString=config_stop_time, ierr=local_err)
call mpas_create_clock(core_clock, startTime=startTime, timeStep=timeStep, stopTime=stopTime, ierr=local_err)
else
call mpas_log_write('Neither config_run_duration nor config_stop_time were specified.', messageType=MPAS_LOG_ERR)
ierr = 1
end if
!TODO: set phyics alarms here...
!....
!....
end subroutine atm_simulation_clock_init
subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
use atm_time_integration
use mpas_rbf_interpolation
use mpas_vector_reconstruction
use mpas_stream_manager
use mpas_atm_boundaries, only : mpas_atm_setup_bdy_masks
#ifdef DO_PHYSICS
! use mpas_atmphys_aquaplanet
use mpas_atmphys_control
use mpas_atmphys_init
use mpas_atmphys_manager
#endif
implicit none
type (dm_info), intent(in) :: dminfo
type (MPAS_streamManager_type), intent(inout) :: stream_manager
type (block_type), intent(inout) :: block
type (mpas_pool_type), intent(inout) :: mesh !MGD does this need to be a pointer?
real (kind=RKIND), intent(in) :: dt
type (mpas_pool_type), pointer :: state
type (mpas_pool_type), pointer :: diag
type (mpas_pool_type), pointer :: tend
type (mpas_pool_type), pointer :: sfc_input
type (mpas_pool_type), pointer :: diag_physics
type (mpas_pool_type), pointer :: diag_physics_noahmp
type (mpas_pool_type), pointer :: ngw_input
type (mpas_pool_type), pointer :: atm_input
type (mpas_pool_type), pointer :: output_noahmp
integer :: iCell,iEdge,iVertex,k
real (kind=RKIND), dimension(:,:), pointer :: u, uReconstructX, uReconstructY, uReconstructZ, uReconstructZonal, uReconstructMeridional
real (kind=RKIND), dimension(:), pointer :: meshScalingDel2, meshScalingDel4
real (kind=RKIND), dimension(:), pointer :: areaCell, invAreaCell
real (kind=RKIND), dimension(:), pointer :: dvEdge, invDvEdge
real (kind=RKIND), dimension(:), pointer :: dcEdge, invDcEdge
real (kind=RKIND), dimension(:), pointer :: areaTriangle, invAreaTriangle
integer, pointer :: nCells_ptr, nEdges_ptr, nVertices_ptr, nVertLevels_ptr, nEdgesSolve
integer :: nCells, nEdges, nVertices, nVertLevels
integer :: thread
character(len=StrKIND), pointer :: mminlu
integer, pointer :: nThreads
integer, dimension(:), pointer :: cellThreadStart, cellThreadEnd
integer, dimension(:), pointer :: cellSolveThreadStart, cellSolveThreadEnd
integer, dimension(:), pointer :: edgeThreadStart, edgeThreadEnd
integer, dimension(:), pointer :: edgeSolveThreadStart, edgeSolveThreadEnd
integer, dimension(:), pointer :: vertexThreadStart, vertexThreadEnd
integer, dimension(:), pointer :: vertexSolveThreadStart, vertexSolveThreadEnd
logical, pointer :: config_do_restart, config_do_DAcycling
call atm_compute_signs(mesh)
call mpas_pool_get_subpool(block % structs, 'diag', diag)
call mpas_pool_get_subpool(block % structs, 'state', state)
call mpas_pool_get_subpool(block % structs, 'state', state)
call mpas_pool_get_config(block % configs, 'config_do_restart', config_do_restart)
call mpas_pool_get_config(block % configs, 'config_do_DAcycling', config_do_DAcycling)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! Compute inverses to avoid divides later
call mpas_pool_get_array(mesh, 'areaCell', areaCell)
call mpas_pool_get_array(mesh, 'invAreaCell', invAreaCell)
call mpas_pool_get_array(mesh, 'dvEdge', dvEdge)
call mpas_pool_get_array(mesh, 'dcEdge', dcEdge)
call mpas_pool_get_array(mesh, 'invDvEdge', invDvEdge)
call mpas_pool_get_array(mesh, 'invDcEdge', invDcEdge)
call mpas_pool_get_array(mesh, 'areaTriangle', areaTriangle)
call mpas_pool_get_array(mesh, 'invAreaTriangle', invAreaTriangle)
call mpas_pool_get_dimension(mesh, 'nCells', nCells_ptr)
call mpas_pool_get_dimension(mesh, 'nEdges', nEdges_ptr)
call mpas_pool_get_dimension(mesh, 'nVertices', nVertices_ptr)
nCells = nCells_ptr
nEdges = nEdges_ptr
nVertices = nVertices_ptr
do iCell=1,nCells
invAreaCell(iCell) = 1.0_RKIND / areaCell(iCell)
end do
do iEdge=1,nEdges
invDvEdge(iEdge) = 1.0_RKIND / dvEdge(iEdge)
end do
do iEdge=1,nEdges
invDcEdge(iEdge) = 1.0_RKIND / dcEdge(iEdge)
end do
do iVertex=1,nVertices
invAreaTriangle(iVertex) = 1.0_RKIND / areaTriangle(iVertex)
end do
!!!!! End compute inverses
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call atm_adv_coef_compression(mesh)
call atm_couple_coef_3rd_order(mesh, block % configs)
call mpas_pool_get_dimension(state, 'nVertices', nVertices_ptr)
call mpas_pool_get_dimension(state, 'nVertLevels', nVertLevels_ptr)
nVertices = nVertices_ptr
nVertLevels = nVertLevels_ptr
allocate(ke_vertex(nVertLevels,nVertices+1)) ! ke_vertex is a module variable defined in mpas_atm_time_integration.F
allocate(ke_edge(nVertLevels,nEdges+1)) ! ke_edge is a module variable defined in mpas_atm_time_integration.F
!$acc parallel default(present)
!$acc loop vector
do k = 1, nVertLevels
ke_vertex(k,nVertices+1) = 0.0_RKIND
ke_edge(k,nEdges+1) = 0.0_RKIND
end do
!$acc end parallel
call mpas_pool_get_dimension(block % dimensions, 'nThreads', nThreads)
call mpas_pool_get_dimension(block % dimensions, 'cellThreadStart', cellThreadStart)
call mpas_pool_get_dimension(block % dimensions, 'cellThreadEnd', cellThreadEnd)
call mpas_pool_get_dimension(block % dimensions, 'cellSolveThreadStart', cellSolveThreadStart)
call mpas_pool_get_dimension(block % dimensions, 'cellSolveThreadEnd', cellSolveThreadEnd)
call mpas_pool_get_dimension(block % dimensions, 'vertexThreadStart', vertexThreadStart)
call mpas_pool_get_dimension(block % dimensions, 'vertexThreadEnd', vertexThreadEnd)
call mpas_pool_get_dimension(block % dimensions, 'vertexSolveThreadStart', vertexSolveThreadStart)
call mpas_pool_get_dimension(block % dimensions, 'vertexSolveThreadEnd', vertexSolveThreadEnd)
call mpas_pool_get_dimension(block % dimensions, 'edgeThreadStart', edgeThreadStart)
call mpas_pool_get_dimension(block % dimensions, 'edgeThreadEnd', edgeThreadEnd)
call mpas_pool_get_dimension(block % dimensions, 'edgeSolveThreadStart', edgeSolveThreadStart)
call mpas_pool_get_dimension(block % dimensions, 'edgeSolveThreadEnd', edgeSolveThreadEnd)
!$OMP PARALLEL DO
do thread=1,nThreads
if (.not. config_do_restart .or. (config_do_restart .and. config_do_DAcycling)) then
call atm_init_coupled_diagnostics(state, 1, diag, mesh, block % configs, &
cellThreadStart(thread), cellThreadEnd(thread), &
vertexThreadStart(thread), vertexThreadEnd(thread), &
edgeThreadStart(thread), edgeThreadEnd(thread), &
cellSolveThreadStart(thread), cellSolveThreadEnd(thread), &
vertexSolveThreadStart(thread), vertexSolveThreadEnd(thread), &
edgeSolveThreadStart(thread), edgeSolveThreadEnd(thread))
end if
call atm_compute_solve_diagnostics(dt, state, 1, diag, mesh, block % configs, &
cellThreadStart(thread), cellThreadEnd(thread), &
vertexThreadStart(thread), vertexThreadEnd(thread), &
edgeThreadStart(thread), edgeThreadEnd(thread))
end do
!$OMP END PARALLEL DO
deallocate(ke_vertex)
deallocate(ke_edge)
call mpas_rbf_interp_initialize(mesh)
call mpas_init_reconstruct(mesh)
call mpas_pool_get_array(state, 'u', u, 1)
call mpas_pool_get_array(diag, 'uReconstructX', uReconstructX)
call mpas_pool_get_array(diag, 'uReconstructY', uReconstructY)
call mpas_pool_get_array(diag, 'uReconstructZ', uReconstructZ)
call mpas_pool_get_array(diag, 'uReconstructZonal', uReconstructZonal)
call mpas_pool_get_array(diag, 'uReconstructMeridional', uReconstructMeridional)
call mpas_reconstruct(mesh, u, &
uReconstructX, &
uReconstructY, &
uReconstructZ, &
uReconstructZonal, &
uReconstructMeridional &
)
#ifdef DO_PHYSICS
!proceed with initialization of physics parameterization if moist_physics is set to true:
call mpas_pool_get_subpool(block % structs, 'sfc_input', sfc_input)
! Before calling physics_init, ensure that mminlu contains the name of the land use dataset
call mpas_pool_get_array(sfc_input, 'mminlu', mminlu)
if (len_trim(mminlu) == 0) then
call mpas_log_write('****************************************************************')
call mpas_log_write('No information on land use dataset is available.')
call mpas_log_write('Assume that we are using ''USGS''.')
call mpas_log_write('****************************************************************')
write(mminlu,'(a)') 'USGS'
end if
if (moist_physics) then
!initialization of some input variables in registry:
call mpas_pool_get_subpool(block % structs, 'tend', tend)
call mpas_pool_get_subpool(block % structs, 'diag_physics', diag_physics)
call mpas_pool_get_subpool(block % structs, 'diag_physics_noahmp', diag_physics_noahmp)
call mpas_pool_get_subpool(block % structs, 'ngw_input', ngw_input)
call mpas_pool_get_subpool(block % structs, 'atm_input', atm_input)
call mpas_pool_get_subpool(block % structs, 'output_noahmp', output_noahmp)
call physics_tables_init(dminfo, block % configs)
call physics_registry_init(mesh, block % configs, sfc_input)
call physics_run_init(block % configs, mesh, state, clock, stream_manager)
!initialization of all physics:
call physics_init(dminfo, stream_manager, clock, block % configs, mesh, diag, tend, state, 1, &
diag_physics, diag_physics_noahmp, ngw_input, atm_input, sfc_input, output_noahmp)
endif
#endif
call atm_compute_mesh_scaling(mesh, block % configs)
call atm_compute_damping_coefs(mesh, block % configs)
!
! Set up mask fields used in limited-area simulations
!
call mpas_atm_setup_bdy_masks(mesh, block % configs)
call mpas_pool_get_dimension(mesh, 'nEdgesSolve', nEdgesSolve)
call mpas_pool_get_array(mesh, 'meshScalingDel2', meshScalingDel2)
call mpas_pool_get_array(mesh, 'meshScalingDel4', meshScalingDel4)
call mpas_log_write('min/max of meshScalingDel2 = $r $r', &
realArgs=(/minval(meshScalingDel2(1:nEdgesSolve)), maxval(meshScalingDel2(1:nEdgesSolve))/))
call mpas_log_write('min/max of meshScalingDel4 = $r $r', &
realArgs=(/minval(meshScalingDel4(1:nEdgesSolve)), maxval(meshScalingDel4(1:nEdgesSolve))/))
end subroutine atm_mpas_init_block
function atm_core_run(domain) result(ierr)
use mpas_timekeeping
use mpas_kind_types
use mpas_stream_manager
use mpas_derived_types, only : MPAS_STREAM_LATEST_BEFORE, MPAS_STREAM_INPUT, MPAS_STREAM_INPUT_OUTPUT
use mpas_timer, only : mpas_timer_start, mpas_timer_stop
use mpas_atm_boundaries, only : mpas_atm_update_bdy_tend
use mpas_atm_diagnostics_manager, only : mpas_atm_diag_update, mpas_atm_diag_compute, mpas_atm_diag_reset
implicit none
type (domain_type), intent(inout) :: domain
integer :: ierr
real (kind=RKIND) :: dt
logical, pointer :: config_do_restart
logical, pointer :: config_apply_lbcs
type (block_type), pointer :: block_ptr
type (MPAS_Time_Type) :: currTime
character(len=StrKIND) :: timeStamp
character (len=StrKIND), pointer :: config_restart_timestamp_name
integer :: itimestep
integer :: iounit
integer :: stream_dir
character(len=StrKIND) :: input_stream, read_time
type (mpas_pool_type), pointer :: state, diag, mesh, diag_physics, tend, tend_physics
! For timing information
real (kind=R8KIND) :: integ_start_time, integ_stop_time
real (kind=R8KIND) :: diag_start_time, diag_stop_time
real (kind=R8KIND) :: input_start_time, input_stop_time
real (kind=R8KIND) :: output_start_time, output_stop_time
ierr = 0
clock => domain % clock
mpas_log_info => domain % logInfo
call mpas_get_timeInterval(mpas_get_clock_timestep(clock, ierr), dt=dt)
call mpas_pool_get_config(domain % blocklist % configs, 'config_do_restart', config_do_restart)
call mpas_pool_get_config(domain % blocklist % configs, 'config_restart_timestamp_name', config_restart_timestamp_name)
! Avoid writing a restart file at the initial time
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='restart', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
! Also, for restart runs, avoid writing the initial history or diagnostics fields to avoid overwriting those from the preceding run
if (config_do_restart) then
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='output', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='diagnostics', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
end if
call mpas_dmpar_get_time(diag_start_time)
if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
block_ptr => domain % blocklist
do while (associated(block_ptr))
call mpas_pool_get_subpool(block_ptr % structs, 'state', state)
call mpas_pool_get_subpool(block_ptr % structs, 'diag', diag)
#ifdef DO_PHYSICS
call mpas_pool_get_subpool(block_ptr % structs, 'diag_physics', diag_physics)
#endif
call mpas_pool_get_subpool(block_ptr % structs, 'mesh', mesh)
call atm_compute_output_diagnostics(state, 1, diag, mesh)
block_ptr => block_ptr % next
end do
end if
call mpas_timer_start('diagnostic_fields')
call mpas_atm_diag_reset()
call mpas_atm_diag_update()
call mpas_atm_diag_compute()
call mpas_timer_stop('diagnostic_fields')
call mpas_dmpar_get_time(diag_stop_time)
call mpas_timer_start('stream_output')
call mpas_dmpar_get_time(output_start_time)
call mpas_stream_mgr_write(domain % streamManager, ierr=ierr)
call mpas_dmpar_get_time(output_stop_time)
call mpas_timer_stop('stream_output')
if (ierr /= MPAS_STREAM_MGR_NOERR .and. &
ierr /= MPAS_STREAM_MGR_ERR_CLOBBER_FILE .and. &
ierr /= MPAS_STREAM_MGR_ERR_CLOBBER_REC) then
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
call mpas_log_write('Error writing one or more output streams', messageType=MPAS_LOG_ERR)
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
end if
if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
call mpas_log_write('Timing for diagnostic computation: $r s', realArgs=(/real(diag_stop_time - diag_start_time, kind=RKIND)/))
call mpas_log_write('Timing for stream output: $r s', realArgs=(/real(output_stop_time - output_start_time, kind=RKIND)/))
end if
call mpas_timer_start('diagnostic_fields')
call mpas_atm_diag_reset()
call mpas_timer_stop('diagnostic_fields')
call mpas_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)
block_ptr => domain % blocklist
call mpas_pool_get_subpool(block_ptr % structs, 'mesh', mesh)
call mpas_pool_get_subpool(block_ptr % structs, 'state', state)
call mpas_pool_get_subpool(block_ptr % structs, 'diag', diag)
#ifdef DO_PHYSICS
call mpas_pool_get_subpool(block_ptr % structs, 'diag_physics', diag_physics)
#endif
call mpas_pool_get_config(domain % blocklist % configs, 'config_apply_lbcs', config_apply_lbcs)
!
! Read initial boundary state
!
if (config_apply_lbcs .and. &
MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='lbc_in', direction=MPAS_STREAM_INPUT, ierr=ierr)) then
block_ptr => domain % blocklist
do while (associated(block_ptr))
call mpas_atm_update_bdy_tend(clock, domain % streamManager, block_ptr, .true., ierr)
if (ierr /= 0) then
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp)
call mpas_log_write('Failed to process LBC data on or before '//trim(timeStamp), messageType=MPAS_LOG_ERR)
return
end if
block_ptr => block_ptr % next
end do
end if
! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
itimestep = 1
do while (.not. mpas_is_clock_stop_time(clock))
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
call mpas_log_write('')
call mpas_log_write('Begin timestep '//trim(timeStamp))
!
! Read future boundary state and compute boundary tendencies
!
if (config_apply_lbcs .and. &
MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='lbc_in', direction=MPAS_STREAM_INPUT, ierr=ierr)) then
block_ptr => domain % blocklist
do while (associated(block_ptr))
call mpas_atm_update_bdy_tend(clock, domain % streamManager, block_ptr, .false., ierr)
if (ierr /= 0) then
call mpas_log_write('Failed to process LBC data at next time after '//trim(timeStamp), messageType=MPAS_LOG_ERR)
return
end if
block_ptr => block_ptr % next
end do
end if
! Regardless of whether boundary tendencies were updated, above, we do not want to read the 'lbc_in' stream
! as a general input stream, below.
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='lbc_in', direction=MPAS_STREAM_INPUT, ierr=ierr)
!
! Read external field updates
!
call MPAS_stream_mgr_begin_iteration(domain % streamManager, ierr=ierr)
do while (MPAS_stream_mgr_get_next_stream(domain % streamManager, streamID=input_stream, directionProperty=stream_dir))
if (stream_dir == MPAS_STREAM_INPUT .or. stream_dir == MPAS_STREAM_INPUT_OUTPUT) then
if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID=input_stream, &
direction=MPAS_STREAM_INPUT, ierr=ierr)) then
call mpas_timer_start('stream_input')
call mpas_dmpar_get_time(input_start_time)
call MPAS_stream_mgr_read(domain % streamManager, streamID=input_stream, whence=MPAS_STREAM_LATEST_BEFORE, &
actualWhen=read_time, ierr=ierr)
call mpas_dmpar_get_time(input_stop_time)
call mpas_timer_stop('stream_input')
if (ierr /= MPAS_STREAM_MGR_NOERR) then
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
call mpas_log_write('Error reading input stream '//trim(input_stream), messageType=MPAS_LOG_ERR)
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
end if
call mpas_log_write('----------------------------------------------------------------------')
call mpas_log_write(' Read '''//trim(input_stream)//''' input stream valid at '//trim(read_time))
call mpas_log_write(' Timing for stream input: $r s', realArgs=(/real(input_stop_time - input_start_time, kind=RKIND)/))
call mpas_log_write('----------------------------------------------------------------------')
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID=input_stream, direction=MPAS_STREAM_INPUT, ierr=ierr)
end if
end if
end do
call mpas_timer_start("time integration")
call mpas_dmpar_get_time(integ_start_time)
call atm_do_timestep(domain, dt, itimestep)
call mpas_dmpar_get_time(integ_stop_time)
call mpas_timer_stop("time integration")
call mpas_log_write(' Timing for integration step: $r s', realArgs=(/real(integ_stop_time - integ_start_time, kind=RKIND)/))
! Move time level 2 fields back into time level 1 for next time step
call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state)
call mpas_pool_shift_time_levels(state)
! Advance clock before writing output
itimestep = itimestep + 1
call mpas_advance_clock(clock)
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
call mpas_dmpar_get_time(diag_start_time)
!
! Write any output streams that have alarms ringing, after computing diagnostics fields
!
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
block_ptr => domain % blocklist
do while (associated(block_ptr))
call mpas_pool_get_subpool(block_ptr % structs, 'state', state)
call mpas_pool_get_subpool(block_ptr % structs, 'diag', diag)
#ifdef DO_PHYSICS
call mpas_pool_get_subpool(block_ptr % structs, 'diag_physics', diag_physics)
call mpas_pool_get_subpool(block_ptr % structs, 'tend_physics', tend_physics)
#endif
call mpas_pool_get_subpool(block_ptr % structs, 'mesh', mesh)
call mpas_pool_get_subpool(block_ptr % structs, 'tend', tend)
call atm_compute_output_diagnostics(state, 1, diag, mesh)
block_ptr => block_ptr % next
end do
end if
call mpas_timer_start('diagnostic_fields')
call mpas_atm_diag_update()
call mpas_atm_diag_compute()
call mpas_timer_stop('diagnostic_fields')
call mpas_dmpar_get_time(diag_stop_time)
call mpas_timer_start('stream_output')
call mpas_dmpar_get_time(output_start_time)
call mpas_stream_mgr_write(domain % streamManager, ierr=ierr)
call mpas_dmpar_get_time(output_stop_time)
call mpas_timer_stop('stream_output')
if (ierr /= MPAS_STREAM_MGR_NOERR .and. &
ierr /= MPAS_STREAM_MGR_ERR_CLOBBER_FILE .and. &
ierr /= MPAS_STREAM_MGR_ERR_CLOBBER_REC) then
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
call mpas_log_write('Error writing one or more output streams', messageType=MPAS_LOG_ERR)
call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
end if
if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
call mpas_log_write('Timing for diagnostic computation: $r s', realArgs=(/real(diag_stop_time - diag_start_time, kind=RKIND)/))
call mpas_log_write('Timing for stream output: $r s', realArgs=(/real(output_stop_time - output_start_time, kind=RKIND)/))
end if
! reset any diagnostics here
if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='diagnostics', direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
block_ptr => domain % blocklist
do while (associated(block_ptr))
#ifdef DO_PHYSICS
call mpas_pool_get_subpool(block_ptr % structs, 'diag_physics', diag_physics)
#endif
call atm_reset_diagnostics(diag_physics)
block_ptr => block_ptr % next
end do
end if
! Only after we've successfully written the restart file should we we
! write the restart_timestamp file
if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='restart', direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
if (domain % dminfo % my_proc_id == 0) then
call mpas_new_unit(iounit)
open(iounit,file=trim(config_restart_timestamp_name),form='formatted',status='replace')
write(iounit,*) trim(timeStamp)
close(iounit)
call mpas_release_unit(iounit)
end if
end if
call mpas_timer_start('diagnostic_fields')
call mpas_atm_diag_reset()
call mpas_timer_stop('diagnostic_fields')
call mpas_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)
end do
end function atm_core_run
subroutine atm_compute_output_diagnostics(state, time_lev, diag, mesh)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields for a domain to be written to history files
!
! Input: state - contains model prognostic fields
! mesh - contains grid metadata
!
! Output: state - upon returning, diagnostic fields will have be computed
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use mpas_constants
implicit none
type (mpas_pool_type), intent(inout) :: state
integer, intent(in) :: time_lev ! which time level to use from state
type (mpas_pool_type), intent(inout) :: diag
type (mpas_pool_type), intent(in) :: mesh
integer :: iCell, k
integer, pointer :: nCells, nVertLevels, index_qv
real (kind=RKIND), dimension(:,:), pointer :: theta, rho, theta_m, rho_zz, zz
real (kind=RKIND), dimension(:,:), pointer :: pressure_base, pressure_p, pressure
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
call mpas_pool_get_dimension(mesh, 'nCells', nCells)
call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)
call mpas_pool_get_dimension(state, 'index_qv', index_qv)
call mpas_pool_get_array(state, 'theta_m', theta_m, time_lev)
call mpas_pool_get_array(state, 'rho_zz', rho_zz, time_lev)
call mpas_pool_get_array(state, 'scalars', scalars, time_lev)
call mpas_pool_get_array(diag, 'theta', theta)
call mpas_pool_get_array(diag, 'rho', rho)
call mpas_pool_get_array(diag, 'pressure_p', pressure_p)
call mpas_pool_get_array(diag, 'pressure_base', pressure_base)
call mpas_pool_get_array(diag, 'pressure', pressure)
call mpas_pool_get_array(mesh, 'zz', zz)
do iCell=1,nCells
do k=1,nVertLevels
theta(k,iCell) = theta_m(k,iCell) / (1._RKIND + rvord * scalars(index_qv,k,iCell))
rho(k,iCell) = rho_zz(k,iCell) * zz(k,iCell)
pressure(k,iCell) = pressure_base(k,iCell) + pressure_p(k,iCell)
end do
end do
end subroutine atm_compute_output_diagnostics
subroutine atm_reset_diagnostics(diag_physics)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! reset some diagnostics after output
!
! Input: diag_physics - contains physics diagnostic fields
!
! Output: whatever diagnostics need resetting after output
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
type (mpas_pool_type), pointer :: diag_physics
real (kind=RKIND), dimension(:), pointer :: refl10cm_1km_max
#ifdef DO_PHYSICS
call mpas_pool_get_array(diag_physics, 'refl10cm_1km_max', refl10cm_1km_max)
if(associated(refl10cm_1km_max)) then
refl10cm_1km_max(:) = 0.
endif
#endif
end subroutine atm_reset_diagnostics
subroutine atm_do_timestep(domain, dt, itimestep)
use mpas_timekeeping
use mpas_kind_types
use atm_time_integration
#ifdef DO_PHYSICS
use mpas_atmphys_control
use mpas_atmphys_driver
use mpas_atmphys_manager
use mpas_atmphys_update
#endif
#ifdef DO_CHEMISTRY
use mpas_atm_chemistry, only: chemistry_step
#endif
use mpas_atm_halos, only: exchange_halo_group
implicit none
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
integer, intent(in) :: itimestep
type (MPAS_Time_Type) :: startTime, currTime
type (MPAS_TimeInterval_Type) :: xtimeTime
integer :: s, s_n, s_d
real (kind=RKIND) :: xtime_s
integer :: ierr
clock => domain % clock
mpas_log_info => domain % logInfo
startTime = mpas_get_clock_time(clock, MPAS_START_TIME, ierr)
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
xtimeTime = currTime - startTime
call mpas_get_timeInterval(interval=xtimeTime, S=s, S_n=s_n, S_d=s_d, ierr=ierr)
xtime_s = (s + s_n / s_d)
#ifdef DO_PHYSICS
!proceed with physics if moist_physics is set to true:
if(moist_physics) then
call physics_timetracker(domain,dt,clock,itimestep,xtime_s)
call physics_driver(domain,itimestep,xtime_s)
endif
#endif
#ifdef DO_CHEMISTRY
call chemistry_step(domain, dt)
#endif
call atm_timestep(domain, dt, currTime, itimestep, exchange_halo_group)
end subroutine atm_do_timestep
function atm_core_finalize(domain) result(ierr)
use mpas_decomp
use mpas_timekeeping
use mpas_atm_diagnostics_manager, only : mpas_atm_diag_cleanup
use mpas_atm_threading, only : mpas_atm_threading_finalize
use atm_time_integration, only : mpas_atm_dynamics_finalize
use mpas_atm_halos, only: atm_destroy_halo_groups
#ifdef DO_PHYSICS
use mpas_atmphys_finalize
#endif
#ifdef DO_CHEMISTRY
use mpas_atm_chemistry, only: chemistry_finalize
#endif
implicit none
type (domain_type), intent(inout) :: domain
type (block_type), pointer :: block_ptr
integer :: ierr
ierr = 0
clock => domain % clock
mpas_log_info => domain % logInfo
!
! Finalize the dynamics
!
call mpas_atm_dynamics_finalize(domain)
!
! Finalize chemistry
!
#ifdef DO_CHEMISTRY
call chemistry_finalize()
#endif
call mpas_atm_diag_cleanup()
call mpas_destroy_clock(clock, ierr)
call mpas_decomp_destroy_decomp_list(domain % decompositions)
#ifdef DO_PHYSICS
block_ptr => domain % blocklist
do while (associated(block_ptr))
call atmphys_finalize(block_ptr%configs)
block_ptr => block_ptr%next
end do
#endif
!
! Destroy halo exchange groups
!
call atm_destroy_halo_groups(domain, ierr)
if (ierr /= 0) then
call mpas_log_write('Failed to destroy halo exchange groups.', messageType=MPAS_LOG_ERR)
end if
!
! Finalize threading
!
call mpas_atm_threading_finalize(domain % blocklist)
call mpas_log_write('')
call mpas_log_write('********************************************************')
call mpas_log_write(' Finished running the atmosphere core')
call mpas_log_write('********************************************************')
end function atm_core_finalize
subroutine atm_compute_mesh_scaling(mesh, configs)
implicit none
type (mpas_pool_type), intent(inout) :: mesh
type (mpas_pool_type), intent(in) :: configs
integer :: iCell,iEdge, cell1, cell2
integer, pointer :: nEdges, nCells
integer, dimension(:,:), pointer :: cellsOnEdge
real (kind=RKIND), dimension(:), pointer :: meshDensity, meshScalingDel2, meshScalingDel4
real (kind=RKIND), dimension(:), pointer :: meshScalingRegionalCell, meshScalingRegionalEdge
logical, pointer :: config_h_ScaleWithMesh
call mpas_pool_get_array(mesh, 'meshDensity', meshDensity)
call mpas_pool_get_array(mesh, 'meshScalingDel2', meshScalingDel2)
call mpas_pool_get_array(mesh, 'meshScalingDel4', meshScalingDel4)
call mpas_pool_get_array(mesh, 'meshScalingRegionalCell', meshScalingRegionalCell)
call mpas_pool_get_array(mesh, 'meshScalingRegionalEdge', meshScalingRegionalEdge)
call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_dimension(mesh, 'nEdges', nEdges)
call mpas_pool_get_dimension(mesh, 'nCells', nCells)
call mpas_pool_get_config(configs, 'config_h_ScaleWithMesh', config_h_ScaleWithMesh)
!
! Compute the scaling factors to be used in the del2 and del4 dissipation
!
meshScalingDel2(:) = 1.0
meshScalingDel4(:) = 1.0
if (config_h_ScaleWithMesh) then
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
meshScalingDel2(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**0.25
meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**0.75
end do
end if
!
! Compute the scaling factors to be used in relaxation zone of regional configuration
!
meshScalingRegionalCell(:) = 1.0
meshScalingRegionalEdge(:) = 1.0
if (config_h_ScaleWithMesh) then
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
meshScalingRegionalEdge(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**0.25
end do
do iCell=1,nCells
meshScalingRegionalCell(iCell) = 1.0 / (meshDensity(iCell))**0.25
end do
end if
end subroutine atm_compute_mesh_scaling
subroutine atm_compute_signs(mesh)
implicit none
type (mpas_pool_type), intent(inout) :: mesh
integer :: i, j, iCell, iVtx
integer, pointer :: nCells, nVertices, nEdges, vertexDegree
integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: edgesOnVertex, verticesOnEdge, cellsOnEdge, edgesOnCell
integer, dimension(:,:), pointer :: verticesOnCell, cellsOnVertex, kiteForCell
real (kind=RKIND), dimension(:,:), pointer :: edgesOnVertex_sign, edgesOnCell_sign
real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3, zb_cell, zb3_cell
call mpas_pool_get_array(mesh, 'edgesOnVertex', edgesOnVertex)
call mpas_pool_get_array(mesh, 'verticesOnEdge', verticesOnEdge)
call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(mesh, 'edgesOnCell', edgesOnCell)
call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(mesh, 'verticesOnCell', verticesOnCell)
call mpas_pool_get_array(mesh, 'cellsOnVertex', cellsOnVertex)
call mpas_pool_get_array(mesh, 'edgesOnVertex_sign', edgesOnVertex_sign)
call mpas_pool_get_array(mesh, 'edgesOnCell_sign', edgesOnCell_sign)
call mpas_pool_get_array(mesh, 'kiteForCell', kiteForCell)
call mpas_pool_get_array(mesh, 'zb', zb)
call mpas_pool_get_array(mesh, 'zb3', zb3)
call mpas_pool_get_array(mesh, 'zb_cell', zb_cell)
call mpas_pool_get_array(mesh, 'zb3_cell', zb3_cell)
call mpas_pool_get_dimension(mesh, 'nCells', nCells)
call mpas_pool_get_dimension(mesh, 'nVertices', nVertices)
call mpas_pool_get_dimension(mesh, 'nEdges', nEdges)
call mpas_pool_get_dimension(mesh, 'vertexDegree', vertexDegree)
do iVtx=1,nVertices
do i=1,vertexDegree
if (edgesOnVertex(i,iVtx) <= nEdges) then
if (iVtx == verticesOnEdge(2,edgesOnVertex(i,iVtx))) then
edgesOnVertex_sign(i,iVtx) = 1.0
else
edgesOnVertex_sign(i,iVtx) = -1.0
end if
else
edgesOnVertex_sign(i,iVtx) = 0.0
end if
end do
end do
do iCell=1,nCells
do i=1,nEdgesOnCell(iCell)
if (edgesOnCell(i,iCell) <= nEdges) then
if (iCell == cellsOnEdge(1,edgesOnCell(i,iCell))) then
edgesOnCell_sign(i,iCell) = 1.0
zb_cell(:,i,iCell) = zb(:,1,edgesOnCell(i,iCell))
zb3_cell(:,i,iCell) = zb3(:,1,edgesOnCell(i,iCell))
else
edgesOnCell_sign(i,iCell) = -1.0
zb_cell(:,i,iCell) = zb(:,2,edgesOnCell(i,iCell))
zb3_cell(:,i,iCell) = zb3(:,2,edgesOnCell(i,iCell))
end if
else
edgesOnCell_sign(i,iCell) = 0.0
end if
end do
end do
do iCell=1,nCells
do i=1,nEdgesOnCell(iCell)
iVtx = verticesOnCell(i,iCell)
if (iVtx <= nVertices) then
do j=1,vertexDegree
if (iCell == cellsOnVertex(j,iVtx)) then
kiteForCell(i,iCell) = j
exit
end if
end do
if (j > vertexDegree) then
call mpas_log_write('Unexpected error while identifying kiteForCell', messageType=MPAS_LOG_ERR)
end if
else
kiteForCell(i,iCell) = 1
end if
end do
end do
end subroutine atm_compute_signs
subroutine atm_compute_damping_coefs(mesh, configs)
implicit none
type (mpas_pool_type), intent(inout) :: mesh
type (mpas_pool_type), intent(in) :: configs
integer :: iCell, k
integer, pointer :: nCells, nVertLevels
real (kind=RKIND), pointer :: config_xnutr, config_zd
real (kind=RKIND) :: z, zt, m1, pii
real (kind=RKIND), dimension(:,:), pointer :: dss, zgrid
real (kind=RKIND), dimension(:), pointer :: meshDensity
real (kind=RKIND) :: dx_scale_power
m1 = -1.0
pii = acos(m1)
call mpas_pool_get_dimension(mesh, 'nCells', nCells)
call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)
call mpas_pool_get_array(mesh, 'meshDensity', meshDensity)
call mpas_pool_get_array(mesh, 'dss', dss)
call mpas_pool_get_array(mesh, 'zgrid', zgrid)
call mpas_pool_get_config(configs, 'config_zd', config_zd)
call mpas_pool_get_config(configs, 'config_xnutr', config_xnutr)
dx_scale_power = 1.0
dss(:,:) = 0.0
do iCell=1,nCells
zt = zgrid(nVertLevels+1,iCell)
do k=1,nVertLevels
z = 0.5*(zgrid(k,iCell) + zgrid(k+1,iCell))
if (z > config_zd) then
dss(k,iCell) = config_xnutr*sin(0.5*pii*(z-config_zd)/(zt-config_zd))**2.0
dss(k,iCell) = dss(k,iCell) / meshDensity(iCell)**(0.25*dx_scale_power)
end if
end do
end do
end subroutine atm_compute_damping_coefs
subroutine atm_adv_coef_compression( mesh )
implicit none
type (mpas_pool_type), intent(inout) :: mesh
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
integer, dimension(:,:), pointer :: cellsOnCell, cellsOnEdge, advCellsForEdge
integer, dimension(:), pointer :: nEdgesOnCell, nAdvCellsForEdge
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge
integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell
integer, pointer :: nCells, nEdges
integer :: cell_list(20)
logical :: addcell
call mpas_pool_get_array(mesh, 'deriv_two', deriv_two)
call mpas_pool_get_array(mesh, 'adv_coefs', adv_coefs)
call mpas_pool_get_array(mesh, 'adv_coefs_3rd', adv_coefs_3rd)
call mpas_pool_get_array(mesh, 'cellsOnCell', cellsOnCell)
call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(mesh, 'advCellsForEdge', advCellsForEdge)
call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(mesh, 'nAdvCellsForEdge', nAdvCellsForEdge)
call mpas_pool_get_array(mesh, 'dcEdge', dcEdge)
call mpas_pool_get_array(mesh, 'dvEdge', dvEdge)
call mpas_pool_get_dimension(mesh, 'nCells', nCells)
call mpas_pool_get_dimension(mesh, 'nEdges', nEdges)
do iEdge = 1, nEdges
nAdvCellsForEdge(iEdge) = 0
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
!
! do only if this edge flux is needed to update owned cells
!
if (cell1 <= nCells .or. cell2 <= nCells) then
cell_list(1) = cell1
cell_list(2) = cell2
n = 2
! add cells surrounding cell 1. n is number of cells currently in list
do i = 1, nEdgesOnCell(cell1)
if (cellsOnCell(i,cell1) /= cell2) then
n = n + 1
cell_list(n) = cellsOnCell(i,cell1)
end if
end do
! add cells surrounding cell 2 (brute force approach)
do iCell = 1, nEdgesOnCell(cell2)
addcell = .true.
do i=1,n
if (cell_list(i) == cellsOnCell(iCell,cell2)) addcell = .false.
end do
if (addcell) then
n = n+1
cell_list(n) = cellsOnCell(iCell,cell2)
end if
end do
nAdvCellsForEdge(iEdge) = n
do iCell = 1, nAdvCellsForEdge(iEdge)
advCellsForEdge(iCell,iEdge) = cell_list(iCell)
end do
! we have the ordered list, now construct coefficients
adv_coefs(:,iEdge) = 0.
adv_coefs_3rd(:,iEdge) = 0.
! pull together third and fourth order contributions to the flux
! first from cell1
j_in = 0
do j=1, n
if( cell_list(j) == cell1 ) j_in = j
end do
adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(1,1,iEdge)
adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(1,1,iEdge)
do iCell = 1, nEdgesOnCell(cell1)
j_in = 0
do j=1, n
if( cell_list(j) == cellsOnCell(iCell,cell1) ) j_in = j
end do
adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
end do
! pull together third and fourth order contributions to the flux
! now from cell2
j_in = 0
do j=1, n
if( cell_list(j) == cell2 ) j_in = j
end do
adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(1,2,iEdge)
adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(1,2,iEdge)
do iCell = 1, nEdgesOnCell(cell2)
j_in = 0
do j=1, n
if( cell_list(j) == cellsOnCell(iCell,cell2) ) j_in = j
end do
adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(iCell+1,2,iEdge)
adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(iCell+1,2,iEdge)
end do
do j = 1,n
adv_coefs (j,iEdge) = - (dcEdge(iEdge) **2) * adv_coefs (j,iEdge) / 12.
adv_coefs_3rd(j,iEdge) = - (dcEdge(iEdge) **2) * adv_coefs_3rd(j,iEdge) / 12.
end do
! 2nd order centered contribution - place this in the main flux weights
j_in = 0
do j=1, n
if( cell_list(j) == cell1 ) j_in = j
end do
adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
j_in = 0
do j=1, n
if( cell_list(j) == cell2 ) j_in = j
end do
adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
! multiply by edge length - thus the flux is just dt*ru times the results of the vector-vector multiply
do j=1,n
adv_coefs (j,iEdge) = dvEdge(iEdge) * adv_coefs (j,iEdge)
adv_coefs_3rd(j,iEdge) = dvEdge(iEdge) * adv_coefs_3rd(j,iEdge)
end do
end if ! only do for edges of owned-cells
end do ! end loop over edges
end subroutine atm_adv_coef_compression
subroutine atm_couple_coef_3rd_order(mesh, configs)
implicit none
type (mpas_pool_type), intent(inout) :: mesh
type (mpas_pool_type), intent(in) :: configs
real (kind=RKIND), dimension(:,:), pointer :: adv_coefs_3rd
real (kind=RKIND), dimension(:,:,:), pointer :: zb3_cell
real (kind=RKIND), pointer :: config_coef_3rd_order
call mpas_pool_get_array(mesh, 'zb3_cell', zb3_cell)
call mpas_pool_get_array(mesh, 'adv_coefs_3rd', adv_coefs_3rd)
call mpas_pool_get_config(configs, 'config_coef_3rd_order', config_coef_3rd_order)
adv_coefs_3rd(:,:) = config_coef_3rd_order * adv_coefs_3rd(:,:)
zb3_cell(:,:,:) = config_coef_3rd_order * zb3_cell(:,:,:)
end subroutine atm_couple_coef_3rd_order
!-----------------------------------------------------------------------
! routine mpas_atm_run_compatibility
!
!> \brief Checks input fields and options for compatibility
!> \author Michael Duda
!> \date 22 October 2018
!> \details
!> This routine checks the input fields and run-time options provided
!> by the user for compatibility. For example, two run-time options may
!> be mutually exclusive, or an option may require that a certain input
!> field is provided.
!>
!> A value of 0 is returned if there are no incompatibilities among
!> the provided input fields and run-time options, and a non-zero value
!> otherwise.
!
!-----------------------------------------------------------------------
subroutine mpas_atm_run_compatibility(dminfo, blockList, streamManager, ierr)
#ifdef DO_PHYSICS
use mpas_atmphys_control, only : physics_compatibility_check
#endif
use mpas_atm_boundaries, only : mpas_atm_bdy_checks
use atm_time_integration, only : mpas_atm_dynamics_checks
implicit none
type (dm_info), pointer :: dminfo
type (block_type), pointer :: blockList
type (MPAS_streamManager_type), pointer :: streamManager
integer, intent(out) :: ierr
integer :: local_ierr
ierr = 0
#ifdef DO_PHYSICS
!
! Physics specific checks found in physics/mpas_atmphys_control.F
!
call physics_compatibility_check(dminfo, blockList, streamManager, local_ierr)
ierr = ierr + local_ierr
#endif
!
! Checks for limited-area simulations
!
call mpas_atm_bdy_checks(dminfo, blockList, streamManager, local_ierr)
ierr = ierr + local_ierr
!
! Checks for dynamics options
!
call mpas_atm_dynamics_checks(dminfo, blockList, streamManager, local_ierr)
ierr = ierr + local_ierr
end subroutine mpas_atm_run_compatibility
end module atm_core