diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 4061d22fccc0..a77a75315dac 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -992,6 +992,15 @@ add_default($nl, 'config_transport_tests_flow_id'); add_default($nl, 'config_vert_levels'); +######################################### +# Namelist group: manufactured_solution # +######################################### + +add_default($nl, 'config_use_manufactured_solution'); +add_default($nl, 'config_manufactured_solution_wavelength_x'); +add_default($nl, 'config_manufactured_solution_wavelength_y'); +add_default($nl, 'config_manufactured_solution_amplitude'); + ################################################ # Namelist group: tracer_forcing_activeTracers # ################################################ @@ -1761,6 +1770,7 @@ my @groups = qw(run_modes testing transport_tests init_mode_vert_levels + manufactured_solution tracer_forcing_activetracers tracer_forcing_debugtracers tracer_forcing_ecosystracers diff --git a/components/mpas-ocean/bld/build-namelist-group-list b/components/mpas-ocean/bld/build-namelist-group-list index 7d5b8cdd473f..11864939025d 100644 --- a/components/mpas-ocean/bld/build-namelist-group-list +++ b/components/mpas-ocean/bld/build-namelist-group-list @@ -39,6 +39,7 @@ my @groups = qw(run_modes testing transport_tests init_mode_vert_levels + manufactured_solution tracer_forcing_activetracers tracer_forcing_debugtracers tracer_forcing_ecosystracers diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index 1bcfc1484d09..465bf12d8eb2 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -509,6 +509,15 @@ add_default($nl, 'config_transport_tests_flow_id'); add_default($nl, 'config_vert_levels'); +######################################### +# Namelist group: manufactured_solution # +######################################### + +add_default($nl, 'config_use_manufactured_solution'); +add_default($nl, 'config_manufactured_solution_wavelength_x'); +add_default($nl, 'config_manufactured_solution_wavelength_y'); +add_default($nl, 'config_manufactured_solution_amplitude'); + ################################################ # Namelist group: tracer_forcing_activeTracers # ################################################ diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index 17e0a56ddeb8..aa5b8ed0f068 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -534,6 +534,12 @@ -1 + +.false. +2000000.0 +2000000.0 +1 + .true. .true. diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index 72369e414459..b8af04867437 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -2573,6 +2573,41 @@ Default: Defined in namelist_defaults.xml + + + +This flag includes additional thickness and velocity tendencies necessary for testing with a manufactured solution. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +Wavelength of manufactured solution in the x direction + +Valid values: Any positive real number +Default: Defined in namelist_defaults.xml + + + +Wavelength of manufactured solution in the y direction + +Valid values: Any positive real number +Default: Defined in namelist_defaults.xml + + + +Amplitude of the manufactured solution + +Valid values: Any positive real number +Default: Defined in namelist_defaults.xml + + + + + + + + + diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F index 0e5235791f14..3fe3395d13c2 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F @@ -78,6 +78,7 @@ module ocn_forward_mode use ocn_gm use ocn_submesoscale_eddies use ocn_stokes_drift + use ocn_manufactured_solution use ocn_high_freq_thickness_hmix_del2 @@ -499,6 +500,8 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{ call ocn_vertical_remap_init(err_tmp) call ocn_stokes_drift_init(err_tmp) ierr = ior(ierr,err_tmp) + call ocn_manufactured_solution_init(domain, err_tmp) + ierr = ior(ierr,err_tmp) if (ierr /= 0) then call mpas_log_write( & diff --git a/components/mpas-ocean/src/ocean.cmake b/components/mpas-ocean/src/ocean.cmake index c1b71dcd5806..6bc035bb5925 100644 --- a/components/mpas-ocean/src/ocean.cmake +++ b/components/mpas-ocean/src/ocean.cmake @@ -105,6 +105,7 @@ list(APPEND RAW_SOURCES core_ocean/shared/mpas_ocn_wetting_drying.F core_ocean/shared/mpas_ocn_vel_tidal_potential.F core_ocean/shared/mpas_ocn_stokes_drift.F + core_ocean/shared/mpas_ocn_manufactured_solution.F ) set(OCEAN_DRIVER diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index a9100a8bf9fe..5c5bc21cff12 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -75,13 +75,14 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_transport_tests.o \ mpas_ocn_vel_self_attraction_loading.o \ mpas_ocn_vertical_advection.o \ - mpas_ocn_stokes_drift.o + mpas_ocn_stokes_drift.o \ + mpas_ocn_manufactured_solution.o all: $(OBJS) mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics.o mpas_ocn_diagnostics_variables.o mpas_ocn_gm.o mpas_ocn_submesoscale_eddies.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_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_tracer_CFC.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_vel_self_attraction_loading.o mpas_ocn_vel_tidal_potential.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_thick_hadv.o mpas_ocn_thick_vadv.o mpas_ocn_vel_hadv_coriolis.o mpas_ocn_vel_pressure_grad.o mpas_ocn_vel_vadv.o mpas_ocn_vel_hmix.o mpas_ocn_vel_forcing.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_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_tracer_CFC.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_vel_self_attraction_loading.o mpas_ocn_vel_tidal_potential.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_thick_hadv.o mpas_ocn_thick_vadv.o mpas_ocn_vel_hadv_coriolis.o mpas_ocn_vel_pressure_grad.o mpas_ocn_vel_vadv.o mpas_ocn_vel_hmix.o mpas_ocn_vel_forcing.o mpas_ocn_manufactured_solution.o mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_vertical_advection.o mpas_ocn_submesoscale_eddies.o @@ -235,6 +236,8 @@ mpas_ocn_vertical_remap.o: mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_diagn mpas_ocn_stokes_drift.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_diagnostics_variables.o mpas_ocn_mesh.o +mpas_ocn_manufactured_solution.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o + clean: $(RM) *.o *.i *.mod *.f90 diff --git a/components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F b/components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F new file mode 100644 index 000000000000..e9ad6786103c --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F @@ -0,0 +1,251 @@ +! 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.io/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_manufactured_solution +! +!> \brief Computes tendency terms for a manufactured solution +!> \author Steven Brus, Carolyn Begeman +!> \date June 2023 +!> \details +!> This module contains the routines for computing the thickness and +!> normal velocity tendencies for a manufactured solution case. See +!> Bishnu et al. 2022 (https://doi.org/10.22541/essoar.167100170.03833124/v1) +! +!----------------------------------------------------------------------- + +module ocn_manufactured_solution + + use mpas_constants + use mpas_timer + use mpas_timekeeping + use ocn_constants + use ocn_config + use ocn_mesh + use ocn_diagnostics_variables + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_manufactured_solution_tend_thick, & + ocn_manufactured_solution_tend_vel, & + ocn_manufactured_solution_init + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + + real (kind=RKIND) :: kx, ky + real (kind=RKIND) :: ang_freq + real (kind=RKIND) :: eta0 + real (kind=RKIND) :: H0 + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_manufactured_solution_tend_thick +! +!> \brief Computes manufactured solution thickness tendency +!> \author Steven Brus, Carolyn Begeman +!> \date June 2023 +!> \details +!> This routine computes the thickness tendency for the manufactured solution +! +!----------------------------------------------------------------------- + + subroutine ocn_manufactured_solution_tend_thick(tend, err)!{{{ + + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: thickness tendency + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: iCell, kmin, kmax, k + real (kind=RKIND) :: phase, time + + ! End preamble + !------------- + ! Begin code + + if (.not. config_use_manufactured_solution) return + + time = daysSinceStartOfSim*86400.0_RKIND + + do iCell = 1,nCellsOwned + + kmin = minLevelCell(iCell) + kmax = maxLevelCell(iCell) + + phase = kx*xCell(iCell) + ky*yCell(iCell) - ang_freq*time + do k = kmin, kmax + tend(k,iCell) = tend(k,iCell) + eta0*(-H0*(kx + ky)*sin(phase) & + - ang_freq*cos(phase) & + + eta0*(kx + ky)*cos(2.0_RKIND*phase)) + enddo + + enddo + + err = 0 + + !-------------------------------------------------------------------- + + end subroutine ocn_manufactured_solution_tend_thick!}}} + +!*********************************************************************** +! +! routine ocn_manufactured_solution_tend_vel +! +!> \brief Computes manufactured solution velocity tendency +!> \author Steven Brus, Carolyn Begeman +!> \date June 2023 +!> \details +!> This routine computes the velocity tendency for the manufactured solution +! +!----------------------------------------------------------------------- + subroutine ocn_manufactured_solution_tend_vel(tend, err)!{{{ + + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: velocity tendency + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: iEdge, kmin, kmax, k + real (kind=RKIND) :: phase, u, v, time + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + if (.not. config_use_manufactured_solution) return + + time = daysSinceStartOfSim*86400.0_RKIND + + do iEdge = 1, nEdgesOwned + + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + + phase = kx*xEdge(iEdge) + ky*yEdge(iEdge) - ang_freq*time + + do k = kmin, kmax + + u = eta0*((-fEdge(iEdge) + gravity*kx)*cos(phase) & + + ang_freq*sin(phase) & + - 0.5_RKIND*eta0*(kx + ky)*sin(2.0_RKIND*(phase))) + v = eta0*((fEdge(iEdge) + gravity*ky)*cos(phase) & + + ang_freq*sin(phase) & + - 0.5_RKIND*eta0*(kx + ky)*sin(2.0_RKIND*(phase))) + + tend(k,iEdge) = tend(k,iEdge) + u*cos(angleEdge(iEdge)) + v*sin(angleEdge(iEdge)) + enddo + + enddo + + err = 0 + + !-------------------------------------------------------------------- + + end subroutine ocn_manufactured_solution_tend_vel!}}} + +!*********************************************************************** +! +! routine ocn_manufactured_solution_init +! +!> \brief Initialize the manufactured solution tendencies +!> \author Steven Brus, Carolyn Begeman +!> \date June 2023 +!> \details +!> This routine initializes constants related to the manufactured +!> solution tendencies +! +!----------------------------------------------------------------------- + subroutine ocn_manufactured_solution_init(domain, err)!{{{ + + type (domain_type), intent(inout) :: domain + integer, intent(out) :: err !< Output: Error flag + + type (block_type), pointer :: block + type (mpas_pool_type), pointer :: verticalMeshPool + real (kind=RKIND), dimension(:,:), pointer :: restingThickness + + if (.not. config_use_manufactured_solution) return + + kx = 2.0_RKIND*pii / config_manufactured_solution_wavelength_x + ky = 2.0_RKIND*pii / config_manufactured_solution_wavelength_y + + eta0 = config_manufactured_solution_amplitude + + + ! This test case assumes that the restingThickness is horizontally uniform + ! and that only one vertical level is used so only one set of indices is + ! used here. + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) + + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + + H0 = restingThickness(1,1) + + block => block % next + enddo + + ang_freq = sqrt(H0*gravity * (kx**2+ky**2)) + + err = 0 + + end subroutine ocn_manufactured_solution_init!}}} + +!*********************************************************************** + +end module ocn_manufactured_solution + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index 2341ea5c28da..fddabab62039 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F @@ -67,6 +67,7 @@ module ocn_tendency use ocn_wetting_drying use ocn_vel_tidal_potential use ocn_vel_self_attraction_loading + use ocn_manufactured_solution implicit none private @@ -225,6 +226,9 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{ call ocn_tidal_forcing_layer_thickness(forcingPool, & tendThick, err) + ! Compute thickness tendency to manufactured solution + call ocn_manufactured_solution_tend_thick(tendThick, err) + #ifdef MPAS_OPENACC !$acc exit data copyout(tendThick, surfaceThicknessFlux, surfaceThicknessFluxRunoff) #endif @@ -447,6 +451,9 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, & sfcStress, kineticEnergyCell, & layerThickEdgeDrag, & layerThickEdgeMean, tendVel, err) + + ! Compute tendency term for manufactured solution + call ocn_manufactured_solution_tend_vel(tendVel, err) ! vertical mixing treated implicitly in a later routine ! adjust total velocity tendency based on wetting/drying