Skip to content

Commit

Permalink
Adds dam break case
Browse files Browse the repository at this point in the history
  • Loading branch information
@caozd999
caozd999 committed Apr 1, 2019
1 parent 01b9180 commit 5f5f4f4
Show file tree
Hide file tree
Showing 13 changed files with 727 additions and 0 deletions.
1 change: 1 addition & 0 deletions src/core_ocean/mode_init/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ UTILS = mpas_ocn_init_spherical_utils.o \

TEST_CASES = mpas_ocn_init_baroclinic_channel.o \
mpas_ocn_init_lock_exchange.o \
mpas_ocn_init_dam_break.o \
mpas_ocn_init_internal_waves.o \
mpas_ocn_init_overflow.o \
mpas_ocn_init_cvmix_WSwSBF.o \
Expand Down
1 change: 1 addition & 0 deletions src/core_ocean/mode_init/Registry.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#include "Registry_lock_exchange.xml"
#include "Registry_internal_waves.xml"
#include "Registry_overflow.xml"
#include "Registry_dam_break.xml"
#include "Registry_global_ocean.xml"
#include "Registry_cvmix_WSwSBF.xml"
#include "Registry_iso.xml"
Expand Down
30 changes: 30 additions & 0 deletions src/core_ocean/mode_init/Registry_dam_break.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
<nml_record name="dam_break" mode="init" configuration="dam_break">
<nml_option name="config_dam_break_vert_levels" type="integer" default_value="1" units="unitless"
description="Number of vertical levels in dam_break case. Default value is 1."
possible_values="Any positive integer number greater than 0."
/>
<nml_option name="config_dam_break_eta0" type="real" default_value="0.6" units="m"
description="Depth of the domain ($H$)."
possible_values="Any real number larger than zero."
/>
<nml_option name="config_dam_break_dc" type="real" default_value="0.04" units="m"
description="gird resolution in meter ($dc$)."
possible_values="Any real number larger than zero."
/>
<nml_option name="config_dam_break_R0" type="real" default_value="24.2" units="m"
description="max wave propagation radius in 10.0s."
possible_values="sqrt(eta*9.8)*10.0."
/>
<nml_option name="config_dam_break_Xl" type="real" default_value="1.0" units="m"
description="The length of dam along the X-direction."
possible_values="Any real number larger than or equal to zero."
/>
<nml_option name="config_dam_break_Yl" type="real" default_value="2.0" units="m"
description="The length of dam along the Y-direction."
possible_values="Any real number larger than or equal to zero."
/>
<nml_option name="config_dam_break_Inlet" type="real" default_value="0.4" units="m"
description="The width of inlet (dam mouth)."
possible_values="Any real number larger than or equal to zero."
/>
</nml_record>
326 changes: 326 additions & 0 deletions src/core_ocean/mode_init/mpas_ocn_init_dam_break.F
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,326 @@
! 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
!
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
! ocn_init_dam_break
!
!> \brief MPAS ocean initialize case -- Zonally periodic Idealized Southern Ocean (dam_break)
!> \author Phillip J. Wolfram, Luke Van Roekel, Todd Ringler
!> \date 09/14/2015
!> \details
!> This module contains the routines for initializing the
!> dam_break initial condition.
!
!-----------------------------------------------------------------------

module ocn_init_dam_break

use mpas_kind_types
use mpas_io_units
use mpas_derived_types
use mpas_pool_routines
use mpas_constants
use mpas_stream_manager
use mpas_dmpar

use ocn_constants
use ocn_init_vertical_grids
use ocn_init_cell_markers

implicit none
private
save

!--------------------------------------------------------------------
!
! Public parameters
!
!--------------------------------------------------------------------

!--------------------------------------------------------------------
!
! Public member functions
!
!--------------------------------------------------------------------

public :: ocn_init_setup_dam_break, &
ocn_init_validate_dam_break
!--------------------------------------------------------------------
!
! Private module variables
!
!--------------------------------------------------------------------

!***********************************************************************

contains
!***********************************************************************
!
! routine ocn_init_setup_dam_break
!
!> \brief Setup for this initial condition
!> \author Zhendong Cao
!> \date 01/17/2019
!> \details
!> This routine sets up the initial conditions for the dam_break configuration.
!
!-----------------------------------------------------------------------

subroutine ocn_init_setup_dam_break(domain, iErr)!{{{

!--------------------------------------------------------------------

type (domain_type), intent(inout) :: domain
integer, intent(out) :: iErr

type (block_type), pointer :: block_ptr
type (mpas_pool_type), pointer :: meshPool
type (mpas_pool_type), pointer :: statePool
type (mpas_pool_type), pointer :: tracersPool
type (mpas_pool_type), pointer :: verticalMeshPool

! local variables
integer :: iCell, k, idx
real (kind=RKIND) :: yMin, yMax, xMin, xMax, dcEdgeMin, dcEdgeMinGlobal
real (kind=RKIND) :: yMinGlobal, yMaxGlobal, yMidGlobal, xMinGlobal, xMaxGlobal
real (kind=RKIND) :: localVar1, localVar2
real (kind=RKIND), dimension(:), pointer :: interfaceLocations
real (kind=RKIND) :: Tscale, dx, dy

! Define config variable pointers
character (len=StrKIND), pointer :: config_init_configuration, config_vertical_grid
real (kind=RKIND), pointer :: config_dam_break_eta0
real (kind=RKIND), pointer :: config_dam_break_R0
real (kind=RKIND), pointer :: config_drying_min_cell_height
real (kind=RKIND), pointer :: config_dam_break_Xl,config_dam_break_Yl,config_dam_break_Inlet
real (kind=RKIND), pointer :: config_dam_break_dc
integer, pointer :: config_dam_break_vert_levels
logical, pointer :: config_write_cull_cell_mask
logical, pointer :: config_use_wetting_drying
real (kind=RKIND), parameter :: eps=1.0e-12

! Define dimension pointers
integer, pointer :: nCellsSolve, nEdgesSolve, nVertLevels, nVertLevelsP1
integer, pointer :: index_temperature, index_salinity

! Define variable pointers
logical, pointer :: on_a_sphere
integer, dimension(:), pointer :: maxLevelCell
real (kind=RKIND), dimension(:), pointer :: ssh
real (kind=RKIND), dimension(:), pointer :: xCell, yCell,refBottomDepth, refZMid, &
vertCoordMovementWeights, bottomDepth, &
fCell, fEdge, fVertex, dcEdge
real (kind=RKIND), dimension(:,:), pointer :: layerThickness, restingThickness
real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers
integer, dimension(:), pointer :: cullCell

iErr = 0

call mpas_pool_get_config(ocnConfigs, 'config_init_configuration', config_init_configuration)

if(config_init_configuration .ne. trim('dam_break')) return

! Get config flag settings
call mpas_pool_get_config(domain % configs, 'config_write_cull_cell_mask', config_write_cull_cell_mask)
call mpas_pool_get_config(ocnConfigs, 'config_vertical_grid', config_vertical_grid)
call mpas_pool_get_config(ocnConfigs, 'config_dam_break_eta0', config_dam_break_eta0)
call mpas_pool_get_config(ocnConfigs, 'config_dam_break_R0', config_dam_break_R0)
call mpas_pool_get_config(ocnConfigs, 'config_drying_min_cell_height', config_drying_min_cell_height)
call mpas_pool_get_config(ocnConfigs, 'config_dam_break_Xl', config_dam_break_Xl)
call mpas_pool_get_config(ocnConfigs, 'config_dam_break_Yl', config_dam_break_Yl)
call mpas_pool_get_config(ocnConfigs, 'config_dam_break_Inlet', config_dam_break_Inlet)
call mpas_pool_get_config(ocnConfigs, 'config_dam_break_dc', config_dam_break_dc)
call mpas_pool_get_config(ocnConfigs, 'config_dam_break_vert_levels', config_dam_break_vert_levels)

! Determine vertical grid for configuration
call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)
call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
call mpas_pool_get_dimension(meshPool, 'nVertLevelsP1', nVertLevelsP1)
call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere)

! you may restrict your case geometry as follows:
if ( on_a_sphere ) call mpas_log_write('The dam_break configuration can only be applied ' &
// 'to a planar mesh. Exiting...', MPAS_LOG_CRIT)

nVertLevels = config_dam_break_vert_levels
nVertLevelsP1 = nVertLevels + 1
allocate(interfaceLocations(nVertLevelsP1))
call ocn_generate_vertical_grid( config_vertical_grid, interfaceLocations, ocnConfigs)

! Intialize the min/max values to large postive and negative values
yMin = 1.0E10_RKIND
yMax = -1.0E10_RKIND
xMin = 1.0E10_RKIND
xMax = -1.0E10_RKIND
dcEdgeMin = 1.0E10_RKIND

! Determine local min and max values.
block_ptr => domain % blocklist
do while (associated(block_ptr))
call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool)
call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve)
call mpas_pool_get_array(meshPool, 'xCell', xCell)
call mpas_pool_get_array(meshPool, 'yCell', yCell)
call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)

yMin = min( yMin, minval(yCell(1:nCellsSolve)))
yMax = max( yMax, maxval(yCell(1:nCellsSolve)))
xMin = min( xMin, minval(xCell(1:nCellsSolve)))
xMax = max( xMax, maxval(xCell(1:nCellsSolve)))
dcEdgeMin = min( dcEdgeMin, minval(dcEdge(1:nEdgesSolve)))

block_ptr => block_ptr % next
end do ! do while(associated(block_ptr))

!--------------------------------------------------------------------
! Use this section to set initial values
!--------------------------------------------------------------------

block_ptr => domain % blocklist
do while(associated(block_ptr))
call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool)
call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool)
call mpas_pool_get_subpool(block_ptr % structs, 'verticalMesh', verticalMeshPool)
call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)

call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
call mpas_pool_get_array(meshPool, 'xCell', xCell)
call mpas_pool_get_array(meshPool, 'yCell', yCell)
call mpas_pool_get_array(meshPool,'bottomDepth',bottomDepth)
call mpas_pool_get_array(meshPool,'vertCoordMovementWeights',vertCoordMovementWeights)
call mpas_pool_get_array(meshPool,'refBottomDepth',refBottomDepth)
call mpas_pool_get_array(meshPool,'maxLevelCell',maxLevelCell)
call mpas_pool_get_array(meshPool,'cullCell',cullCell)

call mpas_pool_get_array(statePool,'layerThickness',layerThickness,1)
call mpas_pool_get_array(statePool,'ssh',ssh,1)
call mpas_pool_get_array(verticalMeshPool,'restingThickness',restingThickness)
call mpas_pool_get_array(verticalMeshPool,'refZMid',refZMid)

!set refBottomDepth and refZMid
do k = 1, nVertLevels
refBottomDepth(k) = config_dam_break_eta0 * interfaceLocations(k+1)
refZMid(k) = -0.5_RKIND * (interfaceLocations(k+1) + interfaceLocations(k)) * config_dam_break_eta0
enddo

!set vertCoordMovementWeights
vertCoordMovementWeights(:) = 1.0_RKIND

!cullCell dam
do iCell = 1, nCellsSolve
cullCell(iCell)=0
if (xCell(iCell).ge.config_dam_break_R0-2.0_RKIND*config_dam_break_dc) then
!! if (yCell(iCell).le.(config_dam_break_R0-0.5_RKIND*config_dam_break_Yl) .or. &
!! yCell(iCell).ge.(config_dam_break_R0+0.5_RKIND*config_dam_break_Yl)) then
if (yCell(iCell).le.(config_dam_break_R0-0.5_RKIND*config_dam_break_Yl-3.0_RKIND*config_dam_break_dc) .or. &
yCell(iCell).ge.(config_dam_break_R0+0.5_RKIND*config_dam_break_Yl+3.0_RKIND*config_dam_break_dc)) then
cullCell(iCell)=1
endif
endif
!cullCell dam mouth
if (xCell(iCell).ge.(config_dam_break_R0-11.0_RKIND/4.0_RKIND*config_dam_break_dc) .AND. &
xCell(iCell).le.(config_dam_break_R0-5.0_RKIND/4.0_RKIND*config_dam_break_dc)) then
if ( yCell(iCell) .le.(config_dam_break_R0-0.5_RKIND*config_dam_break_Inlet) .OR. &
yCell(iCell) .ge. (config_dam_break_R0+0.5_RKIND*config_dam_break_Inlet+config_dam_break_dc)) then
cullCell(iCell)=1
endif
endif

!initial bathymetry !!! need modified, caozd999 01252019
bottomDepth(iCell) = config_dam_break_eta0
if (xCell(iCell).ge.(config_dam_break_R0-2.0_RKIND*config_dam_break_dc)) then
ssh(iCell) = config_dam_break_eta0
else
ssh(iCell) = config_dam_break_vert_levels*config_drying_min_cell_height + eps
endif

! reorient for coordinate system in vertical
ssh(iCell) = -bottomDepth(iCell) + ssh(iCell)

! set maxLevelCell
maxLevelCell(iCell) = config_dam_break_vert_levels

!layerThickness
do k=1, maxLevelCell(iCell)
layerThickness(k,iCell) = max(config_drying_min_cell_height + eps, &
1.0_RKIND/float(maxLevelCell(iCell))*(ssh(iCell)+bottomDepth(iCell)))
restingThickness(k,iCell) = layerThickness(k,iCell)
enddo
enddo

! Determine global min and max values.
call mpas_dmpar_min_real(domain % dminfo, yMin, yMinGlobal)
call mpas_dmpar_max_real(domain % dminfo, yMax, yMaxGlobal)
call mpas_dmpar_min_real(domain % dminfo, xMin, xMinGlobal)
call mpas_dmpar_max_real(domain % dminfo, xMax, xMaxGlobal)
call mpas_dmpar_min_real(domain % dminfo, dcEdgeMin, dcEdgeMinGlobal)

!mark periodic boundaries
if(config_write_cull_cell_mask) then
call ocn_mark_north_boundary(meshPool, yMaxGlobal, dcEdgeMinGlobal, iErr)
call ocn_mark_south_boundary(meshPool, yMinGlobal, dcEdgeMinGlobal, iErr)
call ocn_mark_east_boundary(meshPool, xMaxGlobal, dcEdgeMinGlobal, iErr)
call ocn_mark_west_boundary(meshPool, xMinGlobal, dcEdgeMinGlobal, iErr)
endif
block_ptr => block_ptr % next
end do !!!(do while(associated(block_ptr)))

deallocate(interfaceLocations)

end subroutine ocn_init_setup_dam_break!}}}

!***********************************************************************
!
! routine ocn_init_validate_dam_break
!
!> \brief Validation for this initial condition
!> \author Phillip J. Wolfram, Luke Van Roekel, Todd Ringler
!> \date 09/14/2015
!> \details
!> This routine validates the configuration options for this case.
!
!-----------------------------------------------------------------------

subroutine ocn_init_validate_dam_break(configPool, packagePool, iErr)!{{{

!--------------------------------------------------------------------

type (mpas_pool_type), intent(in) :: configPool, packagePool
integer, intent(out) :: iErr

character (len=StrKIND), pointer :: config_init_configuration
integer, pointer :: config_vert_levels, config_dam_break_vert_levels

iErr = 0

call mpas_pool_get_config(configPool, 'config_init_configuration', config_init_configuration)
if(config_init_configuration .ne. trim('dam_break')) return

call mpas_pool_get_config(configPool, 'config_vert_levels', config_vert_levels)
call mpas_pool_get_config(configPool, 'config_dam_break_vert_levels', config_dam_break_vert_levels)

if(config_vert_levels <= 0 .and. config_dam_break_vert_levels > 0) then
config_vert_levels = config_dam_break_vert_levels
else if (config_vert_levels <= 0) then
call mpas_log_write( 'Validation failed for dam_break. Not given a usable value for vertical levels.', MPAS_LOG_CRIT)
iErr = 1
end if

!--------------------------------------------------------------------

end subroutine ocn_init_validate_dam_break!}}}


!***********************************************************************

end module ocn_init_dam_break

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
! vim: foldmethod=marker
4 changes: 4 additions & 0 deletions src/core_ocean/mode_init/mpas_ocn_init_mode.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ module ocn_init_mode
use ocn_init_lock_exchange
use ocn_init_internal_waves
use ocn_init_overflow
use ocn_init_dam_break
use ocn_init_global_ocean
use ocn_init_cvmix_WSwSBF
use ocn_init_iso
Expand Down Expand Up @@ -269,6 +270,7 @@ function ocn_init_mode_run(domain) result(iErr)!{{{
call ocn_init_setup_lock_exchange(domain, ierr)
call ocn_init_setup_internal_waves(domain, ierr)
call ocn_init_setup_overflow(domain, ierr)
call ocn_init_setup_dam_break(domain, ierr)
call ocn_init_setup_global_ocean(domain, ierr)
call ocn_init_setup_cvmix_WSwSBF(domain, ierr)
call ocn_init_setup_iso(domain, ierr)
Expand Down Expand Up @@ -361,6 +363,8 @@ subroutine ocn_init_mode_validate_configuration(configPool, packagePool, ioconte
iErr = ior(iErr, err_tmp)
call ocn_init_validate_overflow(configPool, packagePool, iocontext, iErr=err_tmp)
iErr = ior(iErr, err_tmp)
call ocn_init_validate_dam_break(configPool, packagePool, iErr=err_tmp)
iErr = ior(iErr, err_tmp)
call ocn_init_validate_global_ocean(configPool, packagePool, iocontext, iErr=err_tmp)
iErr = ior(iErr, err_tmp)
call ocn_init_validate_cvmix_WSwSBF(configPool, packagePool, iocontext, iErr=err_tmp)
Expand Down
Loading

0 comments on commit 5f5f4f4

Please sign in to comment.