From edd16c99aa5f9318b551350a0a2cbeb11fa8d1a7 Mon Sep 17 00:00:00 2001 From: "Phillip J. Wolfram" Date: Thu, 4 Apr 2019 08:16:23 -0600 Subject: [PATCH 1/2] Adds wetting and drying capability to RK4 --- src/core_ocean/Registry.xml | 30 ++ .../mpas_ocn_time_integration_rk4.F | 63 ++- .../mpas_ocn_time_integration_split.F | 2 +- src/core_ocean/shared/Makefile | 7 +- src/core_ocean/shared/mpas_ocn_diagnostics.F | 56 ++- src/core_ocean/shared/mpas_ocn_tendency.F | 25 +- .../shared/mpas_ocn_wetting_drying.F | 461 ++++++++++++++++++ 7 files changed, 625 insertions(+), 19 deletions(-) create mode 100644 src/core_ocean/shared/mpas_ocn_wetting_drying.F diff --git a/src/core_ocean/Registry.xml b/src/core_ocean/Registry.xml index b01dc27d45..a0ae283fa6 100644 --- a/src/core_ocean/Registry.xml +++ b/src/core_ocean/Registry.xml @@ -916,6 +916,32 @@ possible_values="any positive real, typically 1.0e-3" /> + + + + + + + + + domain % blocklist + do while (associated(block)) + call ocn_prevent_drying_rk4(block, dt, rk_substep_weights(rk_step), config_zero_drying_velocity, err) + block => block % next + end do + ! exchange fields for parallelization + call mpas_pool_get_field(statePool, 'normalVelocity', normalVelocityField, 1) + call mpas_dmpar_exch_halo_field(normalVelocityField) + call mpas_pool_get_field(diagnosticsPool, 'wettingVelocity', wettingVelocityField) + call mpas_dmpar_exch_halo_field(wettingVelocityField) + end if + + call mpas_timer_stop("RK4-prevent drying") + call mpas_threading_barrier() + end if ! Compute tendencies for velocity, thickness, and tracers. ! In RK4 notation, we are computing the right hand side f(t,y), @@ -533,6 +566,23 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ ! END RK loop !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! verify that cells are not dry at conclusion of time step + if (config_use_wetting_drying) then + call mpas_timer_start("RK4- check wet cells") + + ! ensure existing layerThickness is valid + if (config_verify_not_dry) then + block => domain % blocklist + do while (associated(block)) + call ocn_wetting_drying_verify(block, config_drying_min_cell_height, err) + block => block % next + end do + end if + + call mpas_timer_stop("RK4- check wet cells") + call mpas_threading_barrier() + end if + call mpas_timer_stop("RK4-main loop") call mpas_threading_barrier() @@ -776,7 +826,7 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(block, dt, rkWeight, err)!{ endif call mpas_threading_barrier() - call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1) + call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1, dt) call mpas_threading_barrier() end subroutine ocn_time_integrator_rk4_compute_vel_tends!}}} @@ -977,7 +1027,7 @@ subroutine ocn_time_integrator_rk4_compute_tends(block, dt, rkWeight, err)!{{{ endif call mpas_threading_barrier() - call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1) + call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1, dt) call mpas_threading_barrier() if (associated(highFreqThicknessProvis)) then @@ -1025,6 +1075,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkWeight, err)!{ real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, layerThicknessProvis, layerThicknessTend real (kind=RKIND), dimension(:, :), pointer :: lowFreqDivergenceCur, lowFreqDivergenceProvis, lowFreqDivergenceTend real (kind=RKIND), dimension(:, :), pointer :: normalTransportVelocity, normalGMBolusVelocity + real (kind=RKIND), dimension(:, :), pointer :: wettingVelocity real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupCur, tracersGroupProvis, tracersGroupTend @@ -1075,6 +1126,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkWeight, err)!{ call mpas_pool_get_array(diagnosticsPool, 'normalTransportVelocity', normalTransportVelocity) call mpas_pool_get_array(diagnosticsPool, 'normalGMBolusVelocity', normalGMBolusVelocity) + call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', wettingVelocity) call mpas_threading_barrier() !$omp do schedule(runtime) private(k) @@ -1082,6 +1134,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkWeight, err)!{ do k = 1, maxLevelEdgeTop(iEdge) normalVelocityProvis(k, iEdge) = normalVelocityCur(k, iEdge) + rkWeight & * normalVelocityTend(k, iEdge) + normalVelocityProvis(k, iEdge) = normalVelocityProvis(k, iEdge) * (1.0_RKIND - wettingVelocity(k, iEdge)) end do end do !$omp end do @@ -1191,13 +1244,14 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ integer, pointer :: nCells, nEdges integer :: iCell, iEdge, k - type (mpas_pool_type), pointer :: statePool, tendPool, meshPool + type (mpas_pool_type), pointer :: statePool, tendPool, meshPool, diagnosticsPool type (mpas_pool_type), pointer :: tracersPool, tracersTendPool real (kind=RKIND), dimension(:, :), pointer :: normalVelocityNew, normalVelocityTend real (kind=RKIND), dimension(:, :), pointer :: layerThicknessNew, layerThicknessTend real (kind=RKIND), dimension(:, :), pointer :: highFreqThicknessNew, highFreqThicknessTend real (kind=RKIND), dimension(:, :), pointer :: lowFreqDivergenceNew, lowFreqDivergenceTend + real (kind=RKIND), dimension(:, :), pointer :: wettingVelocity real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupNew, tracersGroupTend @@ -1216,6 +1270,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ call mpas_pool_get_subpool(block % structs, 'state', statePool) call mpas_pool_get_subpool(block % structs, 'tend', tendPool) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool) call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) call mpas_pool_get_subpool(tendPool, 'tracersTend', tracersTendPool) @@ -1235,6 +1290,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ call mpas_pool_get_array(tendPool, 'highFreqThickness', highFreqThicknessTend) call mpas_pool_get_array(tendPool, 'lowFreqDivergence', lowFreqDivergenceTend) + call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', wettingVelocity) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) @@ -1249,6 +1305,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ !$omp do schedule(runtime) do iEdge = 1, nEdges normalVelocityNew(:, iEdge) = normalVelocityNew(:, iEdge) + rkWeight * normalVelocityTend(:, iEdge) + normalVelocityNew(:, iEdge) = normalVelocityNew(:, iEdge) * (1.0_RKIND - wettingVelocity(:, iEdge)) end do !$omp end do diff --git a/src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F b/src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F index 52f679279e..98e668d2e6 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F +++ b/src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F @@ -506,7 +506,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_pool_get_array(diagnosticsPool, 'layerThicknessEdge', layerThicknessEdge) - call ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshPool, scratchPool, stage1_tend_time) + call ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshPool, scratchPool, stage1_tend_time, dt) block => block % next end do diff --git a/src/core_ocean/shared/Makefile b/src/core_ocean/shared/Makefile index 94a2b6409b..ef7b6130da 100644 --- a/src/core_ocean/shared/Makefile +++ b/src/core_ocean/shared/Makefile @@ -65,13 +65,14 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_time_average_coupled.o \ mpas_ocn_sea_ice.o \ mpas_ocn_framework_forcing.o \ - mpas_ocn_time_varying_forcing.o + mpas_ocn_time_varying_forcing.o \ + mpas_ocn_wetting_drying.o all: $(OBJS) mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_diagnostics.o mpas_ocn_gm.o mpas_ocn_forcing.o mpas_ocn_surface_land_ice_fluxes.o -mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_frazil_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_diagnostics.o +mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_frazil_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_diagnostics_routines.o: mpas_ocn_constants.o @@ -197,6 +198,8 @@ mpas_ocn_framework_forcing.o: mpas_ocn_time_varying_forcing.o: +mpas_ocn_wetting_drying.o: mpas_ocn_diagnostics.o mpas_ocn_gm.o + clean: $(RM) *.o *.i *.mod *.f90 diff --git a/src/core_ocean/shared/mpas_ocn_diagnostics.F b/src/core_ocean/shared/mpas_ocn_diagnostics.F index 0a9555b818..84e2e49fe4 100644 --- a/src/core_ocean/shared/mpas_ocn_diagnostics.F +++ b/src/core_ocean/shared/mpas_ocn_diagnostics.F @@ -59,7 +59,8 @@ module ocn_diagnostics ocn_reconstruct_gm_vectors, & ocn_diagnostics_init, & ocn_compute_kpp_input_fields, & - ocn_validate_state + ocn_validate_state, & + ocn_build_log_filename !-------------------------------------------------------------------- ! @@ -145,8 +146,10 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic integer :: timeLevel integer, pointer :: indexTemperature, indexSalinity logical, pointer :: config_use_cvmix_kpp + logical, pointer :: config_use_wetting_drying real (kind=RKIND), pointer :: config_apvm_scale_factor, config_coef_3rd_order, config_cvmix_kpp_surface_layer_averaging character (len=StrKIND), pointer :: config_pressure_gradient_type + character (len=StrKIND), pointer :: config_thickness_flux_type real (kind=RKIND), pointer :: config_flux_attenuation_coefficient real (kind=RKIND), pointer :: config_flux_attenuation_coefficient_runoff @@ -173,6 +176,8 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient', config_flux_attenuation_coefficient) call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient_runoff', & config_flux_attenuation_coefficient_runoff) + call mpas_pool_get_config(ocnConfigs, 'config_use_wetting_drying', config_use_wetting_drying) + call mpas_pool_get_config(ocnConfigs, 'config_thickness_flux_type', config_thickness_flux_type) call mpas_pool_get_dimension(tracersPool, 'index_temperature', indexTemperature) call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity) @@ -265,18 +270,45 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic nCells = nCellsArray( size(nCellsArray) ) nVertices = nVerticesArray( size(nVerticesArray) ) - - !$omp do schedule(runtime) private(cell1, cell2, k) - do iEdge = 1, nEdges - ! initialize layerThicknessEdge to avoid divide by zero and NaN problems. - layerThicknessEdge(:, iEdge) = -1.0e34_RKIND - cell1 = cellsOnEdge(1,iEdge) - cell2 = cellsOnEdge(2,iEdge) - do k = 1, maxLevelEdgeTop(iEdge) + if (.not. config_use_wetting_drying .or. (config_use_wetting_drying .and. (trim(config_thickness_flux_type) .eq. 'centered'))) then + !$omp do schedule(runtime) private(cell1, cell2, k) + do iEdge = 1, nEdges + ! initialize layerThicknessEdge to avoid divide by zero and NaN problems. + layerThicknessEdge(:, iEdge) = -1.0e34_RKIND + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, maxLevelEdgeTop(iEdge) + ! central differenced layerThicknessEdge(k,iEdge) = 0.5_RKIND * (layerThickness(k,cell1) + layerThickness(k,cell2)) - end do - end do - !$omp end do + end do + end do + !$omp end do + else + if (config_use_wetting_drying .and. (trim(config_thickness_flux_type) .ne. 'centered')) then + if ( trim(config_thickness_flux_type) .eq. 'upwind') then + !$omp do schedule(runtime) private(cell1, cell2, k) + do iEdge = 1, nEdges + ! initialize layerThicknessEdge to avoid divide by zero and NaN problems. + layerThicknessEdge(:, iEdge) = -1.0e34_RKIND + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, maxLevelEdgeTop(iEdge) + ! upwind + if (normalVelocity(k, iEdge) > 0.0_RKIND) then + layerThicknessEdge(k,iEdge) = layerThickness(k,cell1) + elseif (normalVelocity(k, iEdge) < 0.0_RKIND) then + layerThicknessEdge(k,iEdge) = layerThickness(k,cell2) + else + layerThicknessEdge(k,iEdge) = max(layerThickness(k,cell1), layerThickness(k,cell2)) + end if + end do + end do + !$omp end do + end if + else + call mpas_log_write('Thickness flux option of ' // trim(config_thickness_flux_type) // 'is not known', MPAS_LOG_CRIT) + end if + end if coef_3rd_order = config_coef_3rd_order diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 8a23c7252a..4fa5dedc16 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -57,6 +57,7 @@ module ocn_tendency use ocn_vel_hmix use ocn_vel_forcing use ocn_vmix + use ocn_wetting_drying implicit none private @@ -201,7 +202,7 @@ end subroutine ocn_tend_thick!}}} ! !----------------------------------------------------------------------- - subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshPool, scratchPool, timeLevelIn)!{{{ + subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshPool, scratchPool, timeLevelIn, dt)!{{{ implicit none type (mpas_pool_type), intent(inout) :: tendPool !< Input/Output: Tendency structure @@ -209,6 +210,7 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP type (mpas_pool_type), intent(in) :: forcingPool !< Input: Forcing information type (mpas_pool_type), intent(in) :: diagnosticsPool !< Input: Diagnostic information type (mpas_pool_type), intent(in) :: meshPool !< Input: Mesh information + real (kind=RKIND), intent(in) :: dt type (mpas_pool_type), intent(inout) :: scratchPool !< Input: Scratch structure integer, intent(in), optional :: timeLevelIn !< Input: Time level for state fields @@ -217,10 +219,12 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP real (kind=RKIND), dimension(:), pointer :: surfaceStress, surfaceStressMagnitude, surfaceFluxAttenuationCoefficient real (kind=RKIND), dimension(:,:), pointer :: & layerThicknessEdge, normalVelocity, tangentialVelocity, density, potentialDensity, zMid, pressure, & + layerThickness, & tend_normalVelocity, circulation, relativeVorticity, viscosity, kineticEnergyCell, & normalizedRelativeVorticityEdge, normalizedPlanetaryVorticityEdge, & montgomeryPotential, vertAleTransportTop, divergence, vertViscTopOfEdge, & inSituThermalExpansionCoeff, inSituSalineContractionCoeff + real (kind=RKIND), dimension(:,:), pointer :: wettingVelocity real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers integer :: timeLevel @@ -229,6 +233,8 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP integer, pointer :: indexTemperature, indexSalinity, nEdges, nCells logical, pointer :: config_disable_vel_all_tend + logical, pointer :: config_use_wetting_drying + logical, pointer :: config_prevent_drying character (len=StrKIND), pointer :: config_pressure_gradient_type if (present(timeLevelIn)) then @@ -239,6 +245,8 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP call mpas_pool_get_config(ocnConfigs, 'config_disable_vel_all_tend', config_disable_vel_all_tend) call mpas_pool_get_config(ocnConfigs, 'config_pressure_gradient_type', config_pressure_gradient_type) + call mpas_pool_get_config(ocnConfigs, 'config_use_wetting_drying', config_use_wetting_drying) + call mpas_pool_get_config(ocnConfigs, 'config_prevent_drying', config_prevent_drying) call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) @@ -246,6 +254,7 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevel) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevel) call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, timeLevel) call mpas_pool_get_dimension(tracersPool, 'index_temperature', indexTemperature) call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity) @@ -272,6 +281,7 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP call mpas_pool_get_array(forcingPool, 'surfaceStress', surfaceStress) call mpas_pool_get_array(forcingPool, 'surfaceStressMagnitude', surfaceStressMagnitude) + call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', wettingVelocity) ! ! velocity tendency: start accumulating tendency terms ! @@ -291,6 +301,7 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP if(config_disable_vel_all_tend) return call mpas_timer_start("ocn_tend_vel") + call mpas_log_write("ocn_tend_vel start") ! Build bulk forcing surface stress call ocn_surface_bulk_forcing_vel(meshPool, forcingPool, surfaceStress, surfaceStressMagnitude, err) @@ -333,19 +344,31 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP call ocn_vel_hmix_tend(meshPool, scratchPool, divergence, relativeVorticity, normalVelocity, tangentialVelocity, viscosity, & tend_normalVelocity, err) + call mpas_log_write("error start") ! ! velocity tendency: forcing and bottom drag ! call ocn_vel_forcing_tend(meshPool, normalVelocity, surfaceFluxAttenuationCoefficient, & surfaceStress, kineticEnergyCell, layerThicknessEdge, & tend_normalVelocity, err) + call mpas_log_write("error end") + ! + ! velocity tendency: zero if drying + ! + !$omp do schedule(runtime) + do iEdge = 1, nEdges + tend_normalVelocity(:, iEdge) = tend_normalVelocity(:, iEdge) * (1.0_RKIND - wettingVelocity(:, iEdge)) + end do + !$omp end do ! ! velocity tendency: vertical mixing d/dz( nu_v du/dz)) ! call mpas_timer_stop("ocn_tend_vel") + call mpas_log_write("ocn_tend_vel end") call mpas_threading_barrier() + call mpas_log_write("ocn_tend_vel 2") end subroutine ocn_tend_vel!}}} diff --git a/src/core_ocean/shared/mpas_ocn_wetting_drying.F b/src/core_ocean/shared/mpas_ocn_wetting_drying.F new file mode 100644 index 0000000000..f2dad6d1f1 --- /dev/null +++ b/src/core_ocean/shared/mpas_ocn_wetting_drying.F @@ -0,0 +1,461 @@ +! 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_wetting_drying +! +!> \brief MPAS ocean wetting and drying +!> \author Phillip J. Wolfram +!> \date 03/09/2018 +!> \details +!> This module contains the wetting and drying algorithms +!> to prevent cell thicknesses from becoming too small. +!> +! +!----------------------------------------------------------------------- + +module ocn_wetting_drying + + use mpas_kind_types + use mpas_derived_types + use mpas_pool_routines + use mpas_timer + use mpas_io_units + + use ocn_constants + use ocn_diagnostics + use ocn_gm + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_wetting_drying_verify, ocn_prevent_drying_rk4 + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + real(kind=RKIND), parameter :: eps = 1.0e-12_RKIND + + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_wetting_drying_verify +! +!> \brief Verifies that cells are not too dry. +!> \author Phillip J. Wolfram +!> \date 03/09/2018 +!> \details +!> This routine checks that the minimum thickness in a cell is not +!> too small. +! +!----------------------------------------------------------------------- + + subroutine ocn_wetting_drying_verify( block , minHeight, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (block_type), intent(in) :: block !< block for computation + real (kind=RKIND), intent(in) :: minHeight + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), pointer :: statePool, meshPool, tendPool + integer, dimension(:), pointer :: maxLevelCell + real (kind=RKIND), dimension(:), pointer :: sshSubcycleNew + real (kind=RKIND), dimension(:), pointer :: bottomDepth + integer, pointer :: nCellsSolve + integer :: iCell, k + integer :: debugUnit + real (kind=RKIND), dimension(:,:), pointer :: layerThicknessCur + real (kind=RKIND), dimension(:,:), pointer :: layerThicknessNew + real (kind=RKIND), dimension(:,:), pointer :: layerThicknessTend + real (kind=RKIND) :: minThickness, layerThick + character (len=StrKIND) :: debugFilename + character (len=StrKIND), pointer :: config_time_integrator + + !----------------------------------------------------------------- + ! + ! call relevant routines for computing coefficients + ! + !----------------------------------------------------------------- + + call mpas_pool_get_subpool(block % structs, 'state', statePool) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'tend', tendPool) + + call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, timeLevel=1) + call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessNew, timeLevel=2) + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) + call mpas_pool_get_array(statePool, 'sshSubcycle', sshSubcycleNew, 2) + call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_config(ocnConfigs, 'config_time_integrator', config_time_integrator) + + err = 0 + + call mpas_log_write( 'Verifying that cells are not dry...') + + ! check to make sure that there is no layer that is too dry + minThickness = +1.0E34 + do iCell = 1, nCellsSolve + do k = 1, maxLevelCell(iCell) + ! use ssh as a proxy too for baroclinic mode + if (trim(config_time_integrator) == 'split_explicit') then + layerThick = min(layerThicknessNew(k, iCell), (sshSubcycleNew(iCell)+bottomDepth(iCell))/maxLevelCell(iCell)) + else + layerThick = layerThicknessNew(k, iCell) + end if + minThickness = min(minThickness, layerThick) + if ( layerThick < minHeight ) then + ! report error if layerThickness constraint is violated + if (err == 0) then + debugFilename = ocn_build_log_filename('mpas_ocean_layerThickness_block_stats_', block % blockID) + call mpas_new_unit(debugUnit) + open(unit=debugUnit, file=debugFilename, form='formatted', status='unknown') + end if + err = iCell + write(debugUnit, '(A, I5, A, I5, A, ES14.7, A)') 'ERROR: layerThickness too small at iCell=', iCell, & + ', k=', k, ' with thickness of ', layerThick , '.' + end if + end do + end do + + call mpas_log_write('Minimum thickness is $r.', realArgs=(/minThickness/)) + + if ( err > 0) then + ! end the simulation if layerThickness constraint is violated + call mpas_log_write( 'Done, some cells are have dried.') + flush(debugUnit) + close(debugUnit) + call mpas_release_unit(debugUnit) + call mpas_log_write( 'ERROR: Layer thickness smaller than $r, see ' // debugFilename , & + realArgs=(/ minHeight /), messageType=MPAS_LOG_CRIT) + else + call mpas_log_write( 'Done verifying that cells are wet.') + end if + + end subroutine ocn_wetting_drying_verify !}}} + + +!*********************************************************************** +! +! routine ocn_prevent_drying_rk4 +! +!> \brief Prevents velocity tendency from causing cell drying +!> \author Phillip J. Wolfram +!> \date 03/20/2018 +!> \details +!> This routine modifies velocity tendency to prevent cell drying. +! +!----------------------------------------------------------------------- + + subroutine ocn_prevent_drying_rk4(block, dt, rkWeight, config_zero_drying_velocity, err) !{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (block_type), intent(in) :: block + real (kind=RKIND), intent(in) :: dt + real (kind=RKIND), intent(in) :: rkWeight + logical, pointer :: config_zero_drying_velocity + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + + type (mpas_pool_type), pointer :: tendPool + type (mpas_pool_type), pointer :: meshPool + type (mpas_pool_type), pointer :: statePool + type (mpas_pool_type), pointer :: provisStatePool + type (mpas_pool_type), pointer :: diagnosticsPool + real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur + real (kind=RKIND), dimension(:, :), pointer :: layerThicknessProvis + real (kind=RKIND), dimension(:, :), pointer :: layerThicknessEdge + real (kind=RKIND), dimension(:, :), pointer :: normalTransportVelocity + real (kind=RKIND), dimension(:, :), pointer :: normalVelocity + real (kind=RKIND), dimension(:, :), pointer :: wettingVelocity + type (field2DReal), pointer :: wettingVelocityField + + integer, dimension(:), pointer :: maxLevelEdgeTop + integer, dimension(:), pointer :: maxLevelEdgeBot + integer, pointer :: nEdges + integer :: iEdge, k + + err = 0 + + call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool) + call mpas_pool_get_subpool(block % structs, 'tend', tendPool) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'state', statePool) + call mpas_pool_get_subpool(block % structs, 'provis_state', provisStatePool) + + call mpas_pool_get_array(diagnosticsPool, 'normalTransportVelocity', normalTransportVelocity) + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 1) + call mpas_pool_get_array(diagnosticsPool, 'layerThicknessEdge', layerThicknessEdge) + call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', wettingVelocity) + ! use thickness at n because constraint is h_n + dt*T_h > h_min + call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) + call mpas_pool_get_array(provisStatePool, 'layerThickness', layerThicknessProvis, 1) + + call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) + call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot) + + + !$omp do schedule(runtime) + do iEdge = 1, nEdges + wettingVelocity(:, iEdge) = 0.0_RKIND + end do + !$omp end do + call mpas_threading_barrier() + + ! ensure cells stay wet by selectively damping cells with a damping tendency to make sure tendency doesn't dry cells + + call ocn_wetting_drying_wettingVelocity(meshPool, layerThicknessEdge, layerThicknessCur, layerThicknessProvis, & + normalTransportVelocity, rkWeight, wettingVelocity, err) + + ! prevent drying from happening with selective wettingVelocity + !$omp do schedule(runtime) + do iEdge = 1, nEdges + do k = 1, maxLevelEdgeBot(iEdge) + if (abs(normalTransportVelocity(k,iEdge) + wettingVelocity(k,iEdge)) < eps) then + ! prevent spurious flux for close to zero values + normalTransportVelocity(k, iEdge) = 0.0_RKIND + normalVelocity(k, iEdge) = 0.0_RKIND + else if (abs(normalTransportVelocity(k,iEdge) + wettingVelocity(k,iEdge)) <= abs(normalTransportVelocity(k,iEdge))) then + normalTransportVelocity(k, iEdge) = normalTransportVelocity(k, iEdge) + wettingVelocity(k, iEdge) + normalVelocity(k, iEdge) = normalVelocity(k, iEdge) + wettingVelocity(k, iEdge) + end if + + if (abs(wettingVelocity(k, iEdge)) > 0.0_RKIND .and. config_zero_drying_velocity) then + normalTransportVelocity(k, iEdge) = 0.0_RKIND + normalVelocity(k, iEdge) = 0.0_RKIND + end if + + end do + end do + !$omp end do + call mpas_threading_barrier() + + end subroutine ocn_prevent_drying_rk4 !}}} + + +!*********************************************************************** +! +! routine ocn_wetting_drying_wettingVelocity +! +!> \brief Computes velocity to prevent cell drying +!> \author Phillip J. Wolfram +!> \date 03/19/2018 +!> \details +!> This routine adds wetting velocity opposed to drying motion +!> to prevent cells from drying. +! +!----------------------------------------------------------------------- + subroutine ocn_wetting_drying_wettingVelocity(meshPool, layerThicknessEdge, layerThicknessCur, layerThicknessProvis, & + + normalVelocity, dt, wettingVelocity, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: & + meshPool !< Input: horizonal mesh information + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThicknessCur !< Input: layer thickness at old time + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThicknessProvis !< Input: provisional layer thickness + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThicknessEdge !< Input: layerThickness interpolated to an edge + + real (kind=RKIND), dimension(:,:), intent(in) :: & + normalVelocity !< Input: transport + + real (kind=RKIND), intent(in) :: & + dt !< Input: time step + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + wettingVelocity !< Input/Output: velocity wettingVelocityency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iEdge, iCell, k, i + integer, pointer :: nVertLevels, nCells + integer, dimension(:), pointer :: nEdgesOnCell + integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: maxLevelEdgeBot + integer, dimension(:,:), pointer :: edgesOnCell + integer, dimension(:,:), pointer :: edgeSignOnCell + + real (kind=RKIND) :: divOutFlux + real (kind=RKIND) :: invAreaCell + real (kind=RKIND) :: layerThickness + real (kind=RKIND), dimension(:), pointer :: dvEdge + real (kind=RKIND), dimension(:), pointer :: areaCell + + real (kind=RKIND), pointer :: config_drying_min_cell_height + logical, pointer :: config_zero_drying_velocity + + err = 0 + + call mpas_pool_get_config(ocnConfigs, 'config_drying_min_cell_height', config_drying_min_cell_height) + call mpas_pool_get_config(ocnConfigs, 'config_zero_drying_velocity', config_zero_drying_velocity) + + call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) + call mpas_pool_get_array(meshPool, 'areaCell', areaCell) + call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) + call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot) + call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) + call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + + call mpas_threading_barrier() + + ! need predicted transport velocity to limit drying flux + !$omp do schedule(runtime) private(invAreaCell, i, iEdge, k) + do iCell = 1, nCells + invAreaCell = 1.0_RKIND / areaCell(iCell) + ! can switch with maxLevelEdgeBot(iEdge) + do k = 1, maxLevelCell(iCell) + divOutFlux = 0.0_RKIND + layerThickness = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell)) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + if (k <= maxLevelEdgeBot(iEdge)) then + ! only consider divergence flux leaving the cell + if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) < 0.0_RKIND ) then + divOutFlux = divOutFlux + normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) & + * layerThicknessEdge(k, iEdge) * dvEdge(iEdge) * invAreaCell + end if + end if + end do + + ! if layer thickness is too small, limit divergence flux outwards with opposite velocity + if ((layerThickness + dt*divOutFlux ) <= (config_drying_min_cell_height)) then + ! limit divOutFlux out of cell to keep it wet + divOutFlux = abs(divOutFlux) + divOutFlux = (layerThickness - (config_drying_min_cell_height + eps)) / (dt*divOutFlux + eps) + + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + if (k <= maxLevelEdgeBot(iEdge)) then + if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then + ! each outgoing velocity is penalized (but not the incoming, wetting velocities) + ! square the fractional term to make values near zero go to zero much quicker (to prevent threshold from being hit) + wettingVelocity(k, iEdge) = - (min(max(0.0_RKIND, 1.0_RKIND - (divOutFlux*divOutFlux)), 1.0_RKIND)) * normalVelocity(k, iEdge) + ! just go with simple boolean approach for zero wetting velocity for debugging purposes + if (config_zero_drying_velocity) then + wettingVelocity(k, iEdge) = 1.0_RKIND + end if + end if + end if + end do + end if + + end do + end do + !$omp end do + + call mpas_threading_barrier() + + end subroutine ocn_wetting_drying_wettingVelocity !}}} + + +end module ocn_wetting_drying + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| + +! vim: foldmethod=marker From 4ac7cecf61974fe1e0665f07f3ad861652715853 Mon Sep 17 00:00:00 2001 From: ZHENDONG CAO Date: Wed, 13 Feb 2019 11:54:22 -0700 Subject: [PATCH 2/2] Adds dam break with comparisons against obs & ROMS Digitized and script written by Zhendong Cao from Warner, J. C., Defne, Z., Haas, K., & Arango, H. G. (2013). A wetting and drying scheme for ROMS. Computers & geosciences, 58, 54-61. Comparisons are made against laboratory observations and ROMS output detailed in Warner et al (2013). --- src/core_ocean/mode_init/Makefile | 1 + src/core_ocean/mode_init/Registry.xml | 1 + .../mode_init/Registry_dam_break.xml | 30 ++ .../mode_init/mpas_ocn_init_dam_break.F | 325 ++++++++++++++++++ src/core_ocean/mode_init/mpas_ocn_init_mode.F | 4 + src/core_ocean/shared/mpas_ocn_tendency.F | 5 - .../default/004m/config_analysis.xml | 12 + .../dam_break/default/004m/config_driver.xml | 14 + .../dam_break/default/004m/config_forward.xml | 62 ++++ .../dam_break/default/004m/config_init1.xml | 79 +++++ .../dam_break/default/004m/config_init2.xml | 51 +++ .../default/012m/config_analysis.xml | 12 + .../dam_break/default/012m/config_driver.xml | 14 + .../dam_break/default/012m/config_forward.xml | 62 ++++ .../dam_break/default/012m/config_init1.xml | 79 +++++ .../dam_break/default/012m/config_init2.xml | 51 +++ .../dam_break/default/analysis/comparison.py | 95 +++++ .../dam_break/default/analysis/dam_break.png | Bin 0 -> 9973 bytes .../default/analysis/stations/-5A-sim.csv | 26 ++ .../default/analysis/stations/0-sim.csv | 27 ++ .../default/analysis/stations/4-sim.csv | 24 ++ .../default/analysis/stations/8A-sim.csv | 19 + .../default/analysis/stations/C-sim.csv | 25 ++ .../default/analysis/stations/Station-5A.csv | 18 + .../default/analysis/stations/Station0.csv | 22 ++ .../default/analysis/stations/Station4.csv | 29 ++ .../default/analysis/stations/Station8A.csv | 18 + .../default/analysis/stations/StationC.csv | 19 + .../analysis/stations/stationCoords.csv | 5 + 29 files changed, 1124 insertions(+), 5 deletions(-) create mode 100644 src/core_ocean/mode_init/Registry_dam_break.xml create mode 100644 src/core_ocean/mode_init/mpas_ocn_init_dam_break.F create mode 100644 testing_and_setup/compass/ocean/dam_break/default/004m/config_analysis.xml create mode 100644 testing_and_setup/compass/ocean/dam_break/default/004m/config_driver.xml create mode 100644 testing_and_setup/compass/ocean/dam_break/default/004m/config_forward.xml create mode 100644 testing_and_setup/compass/ocean/dam_break/default/004m/config_init1.xml create mode 100644 testing_and_setup/compass/ocean/dam_break/default/004m/config_init2.xml create mode 100644 testing_and_setup/compass/ocean/dam_break/default/012m/config_analysis.xml create mode 100644 testing_and_setup/compass/ocean/dam_break/default/012m/config_driver.xml create mode 100644 testing_and_setup/compass/ocean/dam_break/default/012m/config_forward.xml create mode 100644 testing_and_setup/compass/ocean/dam_break/default/012m/config_init1.xml create mode 100644 testing_and_setup/compass/ocean/dam_break/default/012m/config_init2.xml create mode 100755 testing_and_setup/compass/ocean/dam_break/default/analysis/comparison.py create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/dam_break.png create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/-5A-sim.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/0-sim.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/4-sim.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/8A-sim.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/C-sim.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/Station-5A.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/Station0.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/Station4.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/Station8A.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/StationC.csv create mode 100644 testing_and_setup/compass/ocean/dam_break/default/analysis/stations/stationCoords.csv diff --git a/src/core_ocean/mode_init/Makefile b/src/core_ocean/mode_init/Makefile index 4f435e7ab5..43ac23c4a6 100644 --- a/src/core_ocean/mode_init/Makefile +++ b/src/core_ocean/mode_init/Makefile @@ -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 \ diff --git a/src/core_ocean/mode_init/Registry.xml b/src/core_ocean/mode_init/Registry.xml index 2df1a0ffe0..ca89781f4d 100644 --- a/src/core_ocean/mode_init/Registry.xml +++ b/src/core_ocean/mode_init/Registry.xml @@ -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" diff --git a/src/core_ocean/mode_init/Registry_dam_break.xml b/src/core_ocean/mode_init/Registry_dam_break.xml new file mode 100644 index 0000000000..26eccadaed --- /dev/null +++ b/src/core_ocean/mode_init/Registry_dam_break.xml @@ -0,0 +1,30 @@ + + + + + + + + + diff --git a/src/core_ocean/mode_init/mpas_ocn_init_dam_break.F b/src/core_ocean/mode_init/mpas_ocn_init_dam_break.F new file mode 100644 index 0000000000..9249df17f2 --- /dev/null +++ b/src/core_ocean/mode_init/mpas_ocn_init_dam_break.F @@ -0,0 +1,325 @@ +! 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-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 + 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 diff --git a/src/core_ocean/mode_init/mpas_ocn_init_mode.F b/src/core_ocean/mode_init/mpas_ocn_init_mode.F index 129ea41a13..2cbdaec308 100644 --- a/src/core_ocean/mode_init/mpas_ocn_init_mode.F +++ b/src/core_ocean/mode_init/mpas_ocn_init_mode.F @@ -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 @@ -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) @@ -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) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 4fa5dedc16..e03552ba6d 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -301,7 +301,6 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP if(config_disable_vel_all_tend) return call mpas_timer_start("ocn_tend_vel") - call mpas_log_write("ocn_tend_vel start") ! Build bulk forcing surface stress call ocn_surface_bulk_forcing_vel(meshPool, forcingPool, surfaceStress, surfaceStressMagnitude, err) @@ -344,14 +343,12 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP call ocn_vel_hmix_tend(meshPool, scratchPool, divergence, relativeVorticity, normalVelocity, tangentialVelocity, viscosity, & tend_normalVelocity, err) - call mpas_log_write("error start") ! ! velocity tendency: forcing and bottom drag ! call ocn_vel_forcing_tend(meshPool, normalVelocity, surfaceFluxAttenuationCoefficient, & surfaceStress, kineticEnergyCell, layerThicknessEdge, & tend_normalVelocity, err) - call mpas_log_write("error end") ! ! velocity tendency: zero if drying @@ -366,9 +363,7 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP ! velocity tendency: vertical mixing d/dz( nu_v du/dz)) ! call mpas_timer_stop("ocn_tend_vel") - call mpas_log_write("ocn_tend_vel end") call mpas_threading_barrier() - call mpas_log_write("ocn_tend_vel 2") end subroutine ocn_tend_vel!}}} diff --git a/testing_and_setup/compass/ocean/dam_break/default/004m/config_analysis.xml b/testing_and_setup/compass/ocean/dam_break/default/004m/config_analysis.xml new file mode 100644 index 0000000000..d01ff823ee --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/004m/config_analysis.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/004m/config_driver.xml b/testing_and_setup/compass/ocean/dam_break/default/004m/config_driver.xml new file mode 100644 index 0000000000..4dd53b2570 --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/004m/config_driver.xml @@ -0,0 +1,14 @@ + + + + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/004m/config_forward.xml b/testing_and_setup/compass/ocean/dam_break/default/004m/config_forward.xml new file mode 100644 index 0000000000..3f862cd9ff --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/004m/config_forward.xml @@ -0,0 +1,62 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + init.nc + + + init.nc + + + output + truncate + output.nc + 0000-00-00_00:00:00.3 + + + + + + + + + + + + + init.nc + + + 6 + + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/004m/config_init1.xml b/testing_and_setup/compass/ocean/dam_break/default/004m/config_init1.xml new file mode 100644 index 0000000000..921a3748a5 --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/004m/config_init1.xml @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + mesh.nc + + + output + truncate + 0000_00:00:01 + ocean.nc + + + + + + + + + + + + + + + + + namelist.input.default + namelist.input + + + + -i + s/nx.*=.*/nx = 325,/g + namelist.input + + + -i + s/ny.*=.*/ny = 700,/g + namelist.input + + + -i + s/dc.*=.*/dc = 0.04,/g + namelist.input + + + + + grid.nc + + + + ocean.nc + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/004m/config_init2.xml b/testing_and_setup/compass/ocean/dam_break/default/004m/config_init2.xml new file mode 100644 index 0000000000..d0b4ed2630 --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/004m/config_init2.xml @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + init.nc + + + output + 0000_00:00:01 + truncate + ocean.nc + + + + + + + + + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/012m/config_analysis.xml b/testing_and_setup/compass/ocean/dam_break/default/012m/config_analysis.xml new file mode 100644 index 0000000000..d01ff823ee --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/012m/config_analysis.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/012m/config_driver.xml b/testing_and_setup/compass/ocean/dam_break/default/012m/config_driver.xml new file mode 100644 index 0000000000..4dd53b2570 --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/012m/config_driver.xml @@ -0,0 +1,14 @@ + + + + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/012m/config_forward.xml b/testing_and_setup/compass/ocean/dam_break/default/012m/config_forward.xml new file mode 100644 index 0000000000..ff2c715d82 --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/012m/config_forward.xml @@ -0,0 +1,62 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + init.nc + + + init.nc + + + output + truncate + output.nc + 0000-00-00_00:00:00.3 + + + + + + + + + + + + + init.nc + + + 6 + + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/012m/config_init1.xml b/testing_and_setup/compass/ocean/dam_break/default/012m/config_init1.xml new file mode 100644 index 0000000000..a92495a289 --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/012m/config_init1.xml @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + mesh.nc + + + output + truncate + 0000_00:00:01 + ocean.nc + + + + + + + + + + + + + + + + + namelist.input.default + namelist.input + + + + -i + s/nx.*=.*/nx = 108,/g + namelist.input + + + -i + s/ny.*=.*/ny = 232,/g + namelist.input + + + -i + s/dc.*=.*/dc = 0.12,/g + namelist.input + + + + + grid.nc + + + + ocean.nc + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/012m/config_init2.xml b/testing_and_setup/compass/ocean/dam_break/default/012m/config_init2.xml new file mode 100644 index 0000000000..d0b4ed2630 --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/012m/config_init2.xml @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + init.nc + + + output + 0000_00:00:01 + truncate + ocean.nc + + + + + + + + + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/dam_break/default/analysis/comparison.py b/testing_and_setup/compass/ocean/dam_break/default/analysis/comparison.py new file mode 100755 index 0000000000..74d81a6b5a --- /dev/null +++ b/testing_and_setup/compass/ocean/dam_break/default/analysis/comparison.py @@ -0,0 +1,95 @@ +#!/usr/bin/env python +""" + +Dam break comparison betewen MPAS-O, data, and ROMS results. + +Zhendong Cao and Phillip J. Wolfram +04/05/2019 + +""" + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.animation as animation +import xarray +import pandas as pd +from scipy import spatial +from PIL import Image +from collections import OrderedDict + +# render statically by default +plt.switch_backend('agg') + +## read output.nc +data = xarray.open_dataset('output.nc') +xCell = data.xCell.values +yCell = data.yCell.values +ssh = data.ssh.values +ssh = ssh+0.6 + +## read station coordinates +StaData = pd.read_csv('stations/stationCoords.csv',header=None) +StaName = StaData.iloc[:,0].values +StaCoord = StaData.iloc[:,1:] + +## find the nearest cell of each station +# coordinate shift from MPAS-O grid to dam break case +xCell = 13 - xCell +yCell = 13 - yCell + +matrix = np.array([xCell,yCell]) +tree = spatial.KDTree(list(zip(*matrix))) +StaCell = tree.query(StaCoord)[1] + +## cells representing station locations +Station = OrderedDict(list(zip(StaName,StaCell))) + +ii=0 + +## plot +fig = plt.figure() +for cell in Station: + ii+=1 + plt.subplot(3,2,ii+1) + + plt.tight_layout() + + # MPAS-O simulation results + mpaso=plt.plot(np.linspace(0,100,ssh.shape[0]),ssh[:,Station[cell]],color='#228B22',linewidth=2,alpha=0.6) + + # Measured data + data = pd.read_csv('stations/'+'Station'+cell+'.csv',header=None) + measured=plt.scatter(data[0]*10,data[1],4,marker='o',color='k') + + # ROMS simulation results (Warner et al., 2013) + sim = pd.read_csv('stations/'+cell+'-sim.csv',header=None) + roms=plt.scatter(sim[0]*10,sim[1],4,marker='v',color='b') + + plt.xlim(0,100) + plt.xticks(np.arange(0,101,20), np.arange(0,11,2)) + plt.ylim(0,0.7) + plt.yticks(np.arange(0,0.7,0.2)) + plt.text(35,0.5,'Station '+cell) + + if ii%2==0: + plt.ylabel('h (m)') + if ii>=4: + plt.xlabel('time (s)') + + plt.legend([mpaso[0],measured,roms],['MPAS-O','Measured','ROMS'], + fontsize='xx-small',frameon=False) + +#station location map +im = Image.open('dam_break.png') +im2 = im.resize((650, 300)) +plt.subplot(3,2,1) +plt.imshow(im2,interpolation='bicubic') +#plt.axis('off') +plt.xlabel('x (m)') +plt.ylabel('y (m)') +plt.locator_params(axis='x',nbins=3) +plt.xticks([0,325,650],[4,2,0]) +plt.locator_params(axis='y',nbins=5) +plt.yticks([0,75,150,220,300],reversed([2,1.5,1,0.5,0])) +#plt.show() +plt.savefig('dam_break_comparison.pdf',dpi=600) diff --git a/testing_and_setup/compass/ocean/dam_break/default/analysis/dam_break.png b/testing_and_setup/compass/ocean/dam_break/default/analysis/dam_break.png new file mode 100644 index 0000000000000000000000000000000000000000..a6ee7d969eda22c8e3b9fdc57fea2db44297ffae GIT binary patch literal 9973 zcmeHsRa9J0v+v+;L4#`u?(P=c-Q8`lFoeO~LkRA!fgk~byGwu|3GOnuJ9qeh-;s6J zx^MU8u65>NudbS^>Z-2p>RrDbrKTc_fl7=D001!L<)k$L02mN-oQ{kLebSv%@c;nm zvi4F^YVuN2lxps-Hug@|0DxRnvM!RI<}g8yp@zI&5V9Ofdz40zB1*mt_GiWgWO>T` zzBFv7to0Xf1oALG+5s!%g zg`Jf#K%QoAl7ff2G2oo`uI@GxdoJ`hH$gcrB>>4v9DYT(!)PQvCI(Rr@gyL%i7$`m zYbC2fb-=`n?i4QaHW=`kKavaH-JiVg0j5jwRzqJDAWc!$l#|G&y&;cHGGrEvt(eTO z7^nD^Un-7Cn>Xr=G}0SZyyOg2#|&u5%+Jy8nhX~H+eG;>8M*Hbz;g&2>0Cw_rY{jY zVQ4h7r0*9dzZJMQ>56Y6mP&=3%@MiMQ5{7Y9&#uITDgGZA`50 z9QwmYCR>BAZ)PZu6RrCCB&)E_d~3z6u{ad)8K{$fUzahMP6Ua-DMcD1wEcP?sZ52L zx-}r}c=FdnHgx_?<9 zem110^8WmyocNwYCY$R63pV|rcohbRW=;}zabverv~*UQ1W7uzNqQ$npTDHMbYuoa z<`)LEg!zTh^R198ix_DOmsmvI?U*+-Ji^Y#QmHgN&N@eDPVw&kH8K;ty=HEG7unjXVs5043Z;50{8 z(vhD52McPr-=^Lp`sfRWr1Y^{yyy5rq~h79d@8;~kq|y_#=qo7vElb#vrULZ~n505R_kjn*0U*rs17;rA-k`xXIjNRy z+i~jZyolVaI^32(h=rn=5Q@yg~P!2y=Y+m7_<7k~*!{O`h;t3$=y;YXKe#|Z)3FOm0LR~hB^0pe-9 zod?k;aM+&(+1Y`6m8`@ar4P9<;zxHC^)*O={_uhIVyF0Zus<+$O^^doJ=aLA(Z_;# zfPq)daBOQ>!4gD7XciLGatNAcsWL#3(Jpj$5@tYg9w{@Dqd5KoGH=L{8Q}t2Mi&*3 zGz*ci%Nj_w6joV^CXCd+#2o_q2%@3r{l+c zL^kSymn1Mp%>Sw;v7o>&&X2HK!rEg9Y-2^o&RW<$bJ93&vT%jmG z+Bd^#TS};u{a-*uU=0R)GJE>{_xq^jF_$@2GL4k2Nxj4NRwzHyj3ibhL?t^U{S}T& z_cWXs{J-i&2XZFW#ny#QOa7GJFSyjaD!cv-aG>VNYzZ`wch1WxSSjq7o1U(lKA$D7 zp_|L7(JucX3#pou`i8~rmw_86l_{rLw67j-{hh+xru<-{wTh zd^MX?6Y_UU2rE-do1|SMRs`B@?`i#np4k0@{G*<2o^bE8(PPcW*Vu1k!(v%tHDg)f z3WC~U<|!BRmdsKmG9_yAJJ#)bt$HPU9eR7i$&_O42e&vq(nr&L)A`aUBGNG|v2ZaB zg{#aLlFd(I!$}ZGP~t{nKav>p-P?_EVscGzn3ykFfmh-A+&JI5a`sw<3^OF?e9Vg7TtZoF=FWQB3KV!@$l*3YrxM}k+$l>Bd{{eBgY%pyA}_B-7~VJh{+OU*Duy`!H2Ke4O5%C^FwM)UooM1Gh{yWZok6 zt(Ql(pfJHQe)H4ds7hHwwBQ{!p%kT>)V(&K&_A?+Y{a4kQ9gBKR7EJ_Ay2mQV;2Vb^u zV1f_%I60fFPy#{x3j2KueXmE>ZbHttoZR-db}JJMljRlOmHic-6+S<7@CG>LHtUM+ zHvD$tLS*ZqdvcU*)wO2T^qKWp;92M~H?Ed~hP2*V$u>|Mz-Fc?E=$S8E^@lt6;C-z zX`<%*tI!N$CUH8qR_Y}Ck=K|0K_=2Ml6gyfNGP7d7IRncJpnO}0l{B_5{eW&9sFy` zB}QSUHqkb|9xe-=zFIj+E3K{Qp-dw#)?CBEBQl7{$BBg`^@g(0vWhahSuh|wO zgh$gsYYBttN$x34he>DR%C$zY`ZjVvfDpwfEn*ShtnkPw#OZ5!-I}$KP2v zBQuw=g2ZTi20A3$CFR%q(^a?#C&31ICYPhkV1fW0h|}mFG#jz3-kf@SR!%-k7l<&# zhwh>Uy2DpH0$V6;yHEP7`k?a@JUrq~NA1rmX7WG*~<+v?t3j z3F*ce&NgTKXmPsAnCUe-?^p3g+qed^e8t$(puTx$ z`BV2j6G<>1w-fiGUdP$y~S-CDLx zRv{jpltu8bx6-GoUC`Xw#9a9|Z&f9Jx1+7wB>=b2^HNiLKEKjsQovUq!QhGfFx}h*;9=vQylvcf@um8({p{h00xX(!mv(+;{}QWOq?(*d zCe{Rwx_5V_bTtznNjlRqgh9vT=E4osU0f5T!UqgkHwUDZiSCxetLG|fBR|qtiK12X zp#>s>G+pKe1_29j5zF8qjND|yt3I>9Pat%RxFIf)rCB3x&<0f-1rf*@5m1g0I9V4I z^n3*em)q23oQ*looug;#sut{oL-Ov0I+EGT6&&(%1VNk zuFkCHR<0J-tRQC~wB-N*gh7JPp|iE8IVH&1$;CquBtrd_wy+uNJfn~T-e-Ik3* zKtO@u2yak^itGZS7&{ZV&XdcXgqBwQFwS>g6dyP5nC2 zfBydEr>DKmf6e6L@$ayp0kXaR!^XkN&i0?Sp{l~KTtN+M4_7CzSM%B~_MW1g!v6^W zPxgP!_Af0dS7)HRwTA~(O_bx`D*v7LKjl^3?X98ly-xF=^8cOpKjqczTs>W(Ywd1t zDevNG?GByn-`fA-;Qw9X-?oI=UYF~?E#JR#^A8uA98pwZw*O4DC@Pi=V-f&Bq#-XY zp#_3DT1GL#*Iq?BUCl#o25cZY6zp8*eX&omS1mD*PF4)hC{Ig~vA2u)mY|n%)!6K+ zcibG5lJQXeGleM{wK zdM^%xx>OoSD}R#?Csy4cZZ8#nNE!lf$^hl_hX*EEa)*4&10X=f9AOYajyNuwfv|W$ zY?w{#I4@s*ZdaBbXH+#xHW7Ec$820j5qkEM8uqJCqpF=cIpC&v;`NIIrdXof8l0 z`d3`DK;4Ni0_srpUta-mZ6z%s$@$Qb^0J_Wd&*eUt9l7M(HhSzSrz)LFWXQex!y|U zRh^rz3s6&jY-auH%X@k#kwxcG@T$%LT@Yn2Bqpg>Uj(4UHuL|6^8aO4@)}AUdSCY! zosKxNx&`F841X7=t?-(8Q3eK8_!PRcupSctFsScO? zP&nlX%iqzs%^+2AZ$CbrIQs-aAz8{YYu`3JA4EY&!p=`W9Fh{*wS-H@$Xr2O#Tu)v z^bA*xE)y7ea=uaFwa>J}8(gHS^=oU02$ZbK5^C-0av5z{LP&FA+kA-Fx;yJ}vZlRb z0S_Cy&H6-IKff9~P_G+$I-H{LF3r-{M3!s@46gqgSer^U24guld@mPYRMa(2Q|@t54;z5U$F-2p8h%aCbbf}E$y1>C>4Z<=Q5)gfncg|v#waK?|{ z;rmj24Y?KaW#{OE`=UjB+}pSj?@W8lQUnd1YW=nrNqG%Xg;XE}J!>#rbwz3+K1Oab zxeI}HVs7x}C#}k2sVIvFlqgm>PRGD4CEOfs9G~I_!RpYiaPsQVev{m4v691pJ3F_; zkA3hkom$@s{~0NF5I9K_zQJHCH3F)NUI_0*&&0gT2yWi*MmCtisi%i8$dYQTy(Q@f zm_GV-5r0JvLbbcMpcl5_9~%&h%;iI0el#q6+FQlW{7zc<$v-f%!S?hz)_iZzW)P>1 zKU^k>cU|Yi>=KJKmP&-tTG5J<#esZ>mj58I-`Hr81J+NJIV69yQ|czeTA@!>T)v7H z4lIAjGxw7Sp0-)^nX1F5cb@<27yYPA+3#LfgxnP~WQ9qe%nm24(=I)EGA%FXk|+3u zk!G513kLH5T0v|%5Q@ME^r&MODUdYx@!LkSJ?&I}Tq$2?nOXgrv#Ub=0DJ+p3+`FG ze^shNdzfLuAlyc?r_%8?`@(_^PrQ!1Q+lUn8&$oI*dvbu;8x-avOQC7T5?Ugo2wW$ zTn2Hxjb`HMjAWHO)2;RuZ$(MA@&4kOWp3ha{i@a&O2sHsSagf+WvcuE;ZM?-0Da%P z0BD7$1#0%aaO%N>gwRso8A}5RO{pjX48X2|8a>6^wif-fttAP63|b$P_BMAI43uK< z1vej>q?eipLxu5jvrB9Chgqa;V?%hDHUAj~XJ}x%WIQUu_sv1Dj-ee$$5~ z*#?|lZ2JN7SjOF@cyJRo+yc{s}_4N_G^q2DX)Tc}F>gP^5XUFVEy!;r+BjODMrydw;5PD)8t=drPISMx^L!$w?r!D<-)FKL6pRB6(&@;Ig zw)xq9IQp20*T!6le?iLZ&G`L;?ca|CCi$7t75bPJH3t*X;}o}Os$X%nTNCh^CQ62S z6})8}dfYjxeZACCXG2$!AbmcSo!p+OHc>%dSz%ig?;4jHKWScl;5P{>$;BC6aPG*J zJ*`6|MBQQm66&%&#f@$5F(#|^gUXBUN@<&KJA&B;E+yffT=>|uuD&O;!&C;CgebC2 z4CrV7qLANyEENhq1@tDbpVNh1?3lZCT+9q5j&AARpk1&QX(p73MqBVOWm>M|q)8r2|l20yi`hM2m2HC*4pLU*Z9vlL7s4Wq~ z34gN=jBN5oA9@s5*VgTR$lG)n%60&lTP0S{$Wc@KZd& zwX|^k-omPuI_Q|%n{$S3EKVc6criL=D9dTZSj`u23=FM1LRwdTqKn<#YSrs4bXzP@ z=EP*VNANAS*}IrhwEIjGSW5Xxp7!H!NVvpN-1_1Nm2@j>+4G7ZVpai4=4oIrB@ivjSI7M6{>Q?{iEGs~s^a#V@LWpY87?veawQl-AHg8& z)y-K!U7*>m$0^n`4p$RdPwnQQ+xZsE;gj>5=bD!!JHn4xFdl^Kr2r{FYJ|!c3%3nv zZcNvE0#oHU*+1`0n!VPDEZ$UalWO4bCE}5UEO?08&HNd(W5HzWxwW3sr?<3B0?s(~ zNx8cPRdN*93%|{%l$=j3rL>ucxyvzr+6Z^xVH`GJ`NV)m;S(_DjuUXG+CW-Tt5Tis zocIZs!6sgyQfJ7#=5qYD#TwympjlndG>_)&+)mhuVtRVcTx+vVpi(nUBwQ2g+_qj7 z)4HAb`XU(?j~!)CBR-M`HzIl-E;DB5hf#e8S20+6d=OS96=z0T>(Hwn^7ktDLAJ`T zGP>V;^w!wdA(-E0M9`5%L?g7>(4R1FA7fSh1{=rT|C#scyWaT(O+ky6`tQT%LcP>g zp4&O=ofc;kZ|`+J(SZ*NGuH9(TN@y0!`@%}j(b1fQ9Dn=Ndi*=zZ zwOxfv@=rVf!^^7jusK zg%)4x?!&TGxjbbonz{EA*ydR4d`+lAa*l7Q$NoCBxJop6px97wIFn45UyT52Qn*L+yWrkY zY6l>x51)^v35kk3>c%00Zdy)V_i_Y+vbL7-#&ZJTdUp@|%HcfN>$J67-^FoxL9CrJ zsl-0j*6~%mt^PcIi|b)?_?co;Ll_MXeAZkK?LNfOV4w5#K<9{CzVmsC(F6I~Tta}885DHMD34KmRnloOzMa_w|pt2rp(y82;-PMiUZ-0L4_(I=Lin6 z>4;vGtHloh)_pl6ehiQ7R;nUh!ZyK%>_^Y#L%PlfBJ};}wPctMWVLiU1X8Iy;kL@f z3?hN1xu5wOSB&Qf$#iz|t$|M-5rqNbn}-4)84D@`d>aC6sS7q-tV|}XbRD`USc#Tm%}|9OEA_Xr2U+%OVS$_vg+Rd^HlBf(lKu0LUqSo zGt-9zVY+d=FSc+ThZ*a{N`-TM1qLn{LPCK_gUL3~aCDH+$JLd2gF6qsb-7^ZK_9^g z?TohBl;h1{{Q0q!af9~wPSzruZ#pgsdA1c>a{1$AD(?-Kv7ANY`f|}^^H!5Ofb)Es z=L3aQs}|Bi0!p@|gDvwWy$=uPnaiAx{`qK=h@0PR!UaOxPW4_^Y$pGRL|Jq3l|N@0 zBZ8ctFKc!1euW3h>bBk3SB}x^Eicv)cbe2??S3#ah}bsN+-kYc>$jCUt&^2UVVF4G z7oZU8#1z+tpn;42f-X@`m;>$@)BfP{ZM{vC>3I_>n<5+0J!z<Fh|obJgrHg?yOOmc@(B|jC+1)gxw z_pLgOz19?*;<={h#g5t_kWL*JcAIs*5_$bWHW1TtzPE ze4FlAt6%&4Uq#{{z<3>|*hvtSscDvnjm=$%buK_O`Ynv=+W5G=yP-n97C6KOFb8Hf-rwk&ejxECGHr3{B9j6F=t^&&FhCy_UBPKe%q286fQC=#8Cz8sL~XFWk2 z65FZ|OQRty1CVidZO>#W?r3*TE@&-yfb}}s6RhGkxJB`L1 zJ}~%z4lYFP9B4FbNcAKC2thWBMZYO;9ySu6Org|VbCA_m-enn+%r@cA_G z8!a+fHE)0{wx@Xthu>)5TuWXkDujM={hSPfTdO|_8yR(n#Iq=D-_onw!D(9!#-s4t z2nB2Kk`rn;iw?`+IK_7}^$aONA(t=opQA4pZOwIVVWph~BGWpz^H^0ybzY+APdy_j z(GLkD=T=H&%l8&bI}J(|cHb|JO)X#|$C!qxUn$Ytby-)T*~fP6r+O6e9+k!9@jT$X zwaZ%Qy3`SMypHYW>e>w%YI?-4p zKuOsPl7EsWN$QN z|8S%@XLWU4s)2rbQ#d%O0Yx`VC~4e;bgE24OmA`q+hH3%6!op$$XK2ZuZhER-$5N-0I-8aarC#EZ1%D`u9QIxIN{c^gjg5tkXB30-Yxb2Xg9S4=5pwLysmB=HNFCv{?(Cnqj*K~ zQ1wvv-t=;mZ33mO5Fssge3ozL&@V$OtvTyiY(YOzdxKX4zw(+|^kUREX&G!2wJH2! zhQf*OyJ1@j)A+V#4%aaGe_-^7!uTL~OBvENbNdq(MsSNW^s?fczFpJqvr&geFOP>U zw|Fnoqw}`9?yP&J##>giwXl@X0P*Gn3B!eFX$9b;uKKJ&E&D&%Bu;a%2vgHWR;?(; zuL>XS%@h$Pjb+Tem1-Q~hp+YYKP~L}5o?I(xS&UeR6{(gyIbTXEkmb+oOqtq@^7Q< zH!TwM#Yd#@KS?^>+`_@S*h%`(YKi$CctbkTx5B|QR*r6@u-M`-@Ku`05&hI)A7$<} zMS(rggNGWzDFnS#y5X}N^IfFssopM=6|HK%gxZ^HyCG^kTASy zve;A)XmH+keu==78wS(>sOzD~E2&XeIm5jwc}Ry>78gmyq9{zhM?8&V^$_@Pbxe29 z4gFd|4i)m5$0&Vmohw%t6J&S^nsFq2;JF_z*L*^6-oQPjAEF`BjreY@5D@w)nck@X z#}P8c0{+P(S#W=ni6l#MP4$-37tH(KAscmZXPCiUWlqK4JD%MV*;BRFxVJp#GG4om za?vpoW5@VuVw(Vs?LV}=orf`VCL=4l4hej4ZGTr`Q)Ho&ey3midbl6VBBWA8+M8mn zSnGBpdbpY!xm9WFcDiKTe0yK`+>3Drf?OjX&f|&@gYw>EE*+iO1U4OVxj+*bh6lC4 zuMx;bn`A*zKWO-;^DSPb;6IxHT5a>TB9B