diff --git a/src/core_ocean/Registry.xml b/src/core_ocean/Registry.xml index b01dc27d45..d57eee1ae9 100644 --- a/src/core_ocean/Registry.xml +++ b/src/core_ocean/Registry.xml @@ -771,6 +771,32 @@ possible_values="Real Values greater than zero less than bottomDepth" /> + + + + + + + + + @@ -3003,6 +3030,18 @@ packages="landIceCouplingPKG" /> + + + + + + diff --git a/src/core_ocean/driver/mpas_ocn_core_interface.F b/src/core_ocean/driver/mpas_ocn_core_interface.F index 6eb6e04063..c9250acfd8 100644 --- a/src/core_ocean/driver/mpas_ocn_core_interface.F +++ b/src/core_ocean/driver/mpas_ocn_core_interface.F @@ -119,6 +119,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: landIceCouplingPKGActive logical, pointer :: thicknessBulkPKGActive logical, pointer :: frazilIceActive + logical, pointer :: tidalForcingActive logical, pointer :: inSituEOSActive logical, pointer :: variableShortwaveActive logical, pointer :: gmActive @@ -145,6 +146,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: config_use_freq_filtered_thickness logical, pointer :: config_use_frazil_ice_formation + logical, pointer :: config_use_tidal_forcing logical, pointer :: config_use_standardGM logical, pointer :: config_use_time_varying_atmospheric_forcing @@ -272,6 +274,15 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ frazilIceActive = .true. end if + ! + ! test for use of tidal forcing, tidalForcingActive + ! + call mpas_pool_get_package(packagePool, 'tidalForcingActive', tidalForcingActive) + call mpas_pool_get_config(configPool, 'config_use_tidal_forcing', config_use_tidal_forcing) + if (config_use_tidal_forcing) then + tidalForcingActive = .true. + end if + ! ! test for form of pressure gradient computation ! diff --git a/src/core_ocean/mode_forward/mpas_ocn_forward_mode.F b/src/core_ocean/mode_forward/mpas_ocn_forward_mode.F index bac96c6e62..5011df5122 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_forward_mode.F +++ b/src/core_ocean/mode_forward/mpas_ocn_forward_mode.F @@ -54,6 +54,7 @@ module ocn_forward_mode use ocn_surface_bulk_forcing use ocn_surface_land_ice_fluxes use ocn_frazil_forcing + use ocn_tidal_forcing use ocn_tracer_hmix use ocn_tracer_hmix_redi @@ -218,6 +219,8 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{ ierr = ior(ierr, err_tmp) call ocn_frazil_forcing_init(err_tmp) ierr = ior(ierr, err_tmp) + call ocn_tidal_forcing_init(err_tmp) + ierr = ior(ierr, err_tmp) call ocn_tracer_hmix_init(err_tmp) ierr = ior(ierr, err_tmp) @@ -583,6 +586,8 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{ call mpas_timer_stop("land_ice_build_arrays") call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, diagnosticsPool, statePool, err) + + call ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, diagnosticsPool, statePool, err) block_ptr => block_ptr % next end do diff --git a/src/core_ocean/mode_init/Makefile b/src/core_ocean/mode_init/Makefile index 4f435e7ab5..00bb43bcf0 100644 --- a/src/core_ocean/mode_init/Makefile +++ b/src/core_ocean/mode_init/Makefile @@ -22,7 +22,8 @@ TEST_CASES = mpas_ocn_init_baroclinic_channel.o \ mpas_ocn_init_sea_mount.o \ mpas_ocn_init_global_ocean.o \ mpas_ocn_init_isomip.o \ - mpas_ocn_init_isomip_plus.o + mpas_ocn_init_isomip_plus.o \ + mpas_ocn_init_tidal_boundary.o #mpas_ocn_init_TEMPLATE.o all: init_mode @@ -73,6 +74,8 @@ mpas_ocn_init_isomip_plus.o: $(UTILS) mpas_ocn_init_ziso.o: $(UTILS) +mpas_ocn_init_tidal_boundary.o: $(UTILS) + #mpas_ocn_init_TEMPLATE.o: $(UTILS) clean: diff --git a/src/core_ocean/mode_init/Registry.xml b/src/core_ocean/mode_init/Registry.xml index 2df1a0ffe0..6d827cbcae 100644 --- a/src/core_ocean/mode_init/Registry.xml +++ b/src/core_ocean/mode_init/Registry.xml @@ -13,4 +13,5 @@ #include "Registry_sea_mount.xml" #include "Registry_isomip.xml" #include "Registry_isomip_plus.xml" +#include "Registry_tidal_boundary.xml" // #include "Registry_TEMPLATE.xml" diff --git a/src/core_ocean/mode_init/Registry_tidal_boundary.xml b/src/core_ocean/mode_init/Registry_tidal_boundary.xml new file mode 100644 index 0000000000..56a2791d79 --- /dev/null +++ b/src/core_ocean/mode_init/Registry_tidal_boundary.xml @@ -0,0 +1,30 @@ + + + + + + + + + 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..5033d7c20c 100644 --- a/src/core_ocean/mode_init/mpas_ocn_init_mode.F +++ b/src/core_ocean/mode_init/mpas_ocn_init_mode.F @@ -55,6 +55,7 @@ module ocn_init_mode use ocn_init_sea_mount use ocn_init_isomip use ocn_init_isomip_plus + use ocn_init_tidal_boundary implicit none private @@ -280,6 +281,7 @@ function ocn_init_mode_run(domain) result(iErr)!{{{ call ocn_init_setup_sea_mount(domain, ierr) call ocn_init_setup_isomip(domain, ierr) call ocn_init_setup_isomip_plus(domain, ierr) + call ocn_init_setup_tidal_boundary(domain, ierr) !call ocn_init_setup_TEMPLATE(domain, ierr) call mpas_log_write( ' Completed setup of: ' // trim(config_init_configuration)) @@ -383,6 +385,8 @@ subroutine ocn_init_mode_validate_configuration(configPool, packagePool, ioconte iErr = ior(iErr, err_tmp) call ocn_init_validate_isomip_plus(configPool, packagePool, iocontext, iErr=err_tmp) iErr = ior(iErr, err_tmp) + call ocn_init_validate_tidal_boundary(configPool, packagePool, iocontext, iErr=err_tmp) + iErr = ior(iErr, err_tmp) ! call ocn_init_validate_TEMPLATE(configPool, packagePool, iocontext, iErr=err_tmp) ! iErr = ior(iErr, err_tmp) end subroutine ocn_init_mode_validate_configuration!}}} diff --git a/src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F b/src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F new file mode 100644 index 0000000000..ae0bc24c77 --- /dev/null +++ b/src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F @@ -0,0 +1,372 @@ +! 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_tidal_boundary +! +!> \brief MPAS ocean initialize case -- Tests tidal boundary +!> \author Phillip Wolfram +!> \date 04/05/2019 +!> \details +!> This module contains the routines for initializing the +!> the tidal_boundary test case +! +!----------------------------------------------------------------------- + +module ocn_init_tidal_boundary + + use mpas_kind_types + use mpas_io_units + use mpas_derived_types + use mpas_pool_routines + 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_tidal_boundary, & + ocn_init_validate_tidal_boundary + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_init_setup_tidal_boundary +! +!> \brief Setup for Tidal Boundary on drying slope +!> \author Phillip Wolfram +!> \date 04/05/2019 +!> \details +!> This routine sets up the initial conditions for the tidal_boundary test case. +! +!----------------------------------------------------------------------- + subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ + + !-------------------------------------------------------------------- + + implicit none + + type (domain_type), intent(inout) :: domain + integer, intent(out) :: iErr + real (kind=RKIND) :: yMin, yMax, dcEdgeMin + real (kind=RKIND) :: yMinGlobal, yMaxGlobal, dcEdgeMinGlobal + real (kind=RKIND) :: plugWidth + real (kind=RKIND) :: slopeCenter, slopeWidth, depCenter, depWidth + + type (block_type), pointer :: block_ptr + + type (mpas_pool_type), pointer :: meshPool + type (mpas_pool_type), pointer :: forcingPool + type (mpas_pool_type), pointer :: statePool + type (mpas_pool_type), pointer :: diagnosticsPool + type (mpas_pool_type), pointer :: verticalMeshPool + type (mpas_pool_type), pointer :: tracersPool + + integer :: iCell, k + + ! Define dimensions + integer, pointer :: nCellsSolve, nEdgesSolve, nVertLevels, nVertLevelsP1 + integer, pointer :: index_temperature, index_salinity, index_tracer1 + + ! Define arrays + integer, dimension(:), pointer :: maxLevelCell + real (kind=RKIND), dimension(:), pointer :: yCell, refBottomDepth, bottomDepth, vertCoordMovementWeights, dcEdge + real (kind=RKIND), dimension(:), pointer :: tidalInputMask + real (kind=RKIND), dimension(:), pointer :: ssh + real (kind=RKIND), dimension(:,:), pointer :: layerThickness, restingThickness, zMid + real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers, debugTracers + + ! Define configs + character (len=StrKIND), pointer :: config_init_configuration, & + config_vertical_grid, & + config_tidal_boundary_layer_type + logical, pointer :: config_tidal_boundary_use_distances, config_alter_ICs_for_pbcs + real (kind=RKIND), pointer :: config_tidal_boundary_right_bottom_depth, & + config_tidal_boundary_left_bottom_depth, & + config_tidal_boundary_plug_width_dist, config_tidal_boundary_plug_width_frac, & + config_tidal_boundary_plug_temperature, config_tidal_boundary_domain_temperature, & + config_tidal_boundary_salinity, config_tidal_boundary_isopycnal_min_thickness, & + config_tidal_boundary_water_thickness + + + real (kind=RKIND), dimension(:), pointer :: interfaceLocations + + iErr = 0 + + + call mpas_pool_get_config(ocnConfigs, 'config_init_configuration', config_init_configuration) + + call mpas_log_write('here') + if(config_init_configuration .ne. trim('tidal_boundary')) return + + call mpas_pool_get_config(ocnConfigs, 'config_vertical_grid', config_vertical_grid) + + call mpas_pool_get_config(ocnConfigs, 'config_tidal_boundary_right_bottom_depth', config_tidal_boundary_right_bottom_depth) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_boundary_left_bottom_depth', config_tidal_boundary_left_bottom_depth) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_boundary_plug_width_frac', config_tidal_boundary_plug_width_frac) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_boundary_plug_temperature', config_tidal_boundary_plug_temperature) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_boundary_domain_temperature', config_tidal_boundary_domain_temperature) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_boundary_salinity', config_tidal_boundary_salinity) + call mpas_pool_get_config(ocnConfigs, 'config_alter_ICs_for_pbcs', config_alter_ICs_for_pbcs) + + call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool) + call mpas_pool_get_dimension(meshPool, 'nVertLevelsP1', nVertLevelsP1) + + allocate(interfaceLocations(nVertLevelsP1)) + call ocn_generate_vertical_grid(config_vertical_grid, interfaceLocations) + + ! Initalize y values to large positive and negative values + yMin = 1.0E10_RKIND + yMax = -1.0E10_RKIND + dcEdgeMin = 1.0E10_RKIND + + ! Determine local min and max y value. + 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, 'yCell', yCell) + call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) + + yMin = min( yMin, minval(yCell(1:nCellssolve))) + yMax = max( yMax, maxval(yCell(1:nCellssolve))) + dcEdgeMin = min( dcEdgeMin, minval(dcEdge(1:nEdgesSolve))) + + block_ptr => block_ptr % next + end do + + ! Determine global min and max y value. This is so the domain + ! can be split into right and left. + 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, dcEdgeMin, dcEdgeMinGlobal) + + plugWidth = (yMaxGlobal - yMinGlobal) * config_tidal_boundary_plug_width_frac + + 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, 'forcing', forcingPool) + call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask) + + call mpas_pool_get_array(meshPool, 'yCell', yCell) + call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool) + call mpas_pool_get_subpool(block_ptr % structs, 'diagnostics', diagnosticsPool) + call mpas_pool_get_subpool(block_ptr % structs, 'verticalMesh', verticalMeshPool) + call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) + + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + + call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature) + call mpas_pool_get_dimension(tracersPool, 'index_salinity', index_salinity) + call mpas_pool_get_dimension(tracersPool, 'index_tracer1', index_tracer1) + + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'refBottomDepth', refBottomDepth) + call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_array(meshPool, 'vertCoordMovementWeights', vertCoordMovementWeights) + + call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1) + call mpas_pool_get_array(tracersPool, 'debugTracers', debugTracers, 1) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) + call mpas_pool_get_array(diagnosticsPool, 'zMid', zMid) + call mpas_pool_get_array(statePool, 'ssh', ssh, 1) + + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + + + if (config_tidal_boundary_right_bottom_depth < config_tidal_boundary_left_bottom_depth) then + call mpas_log_write('Right boundary must be deeper than left boundary!', MPAS_LOG_CRIT) + end if + + ! Set refBottomDepth, bottomDepth, and maxLevelCell + do k = 1, nVertLevels + refBottomDepth(k) = config_tidal_boundary_right_bottom_depth * interfaceLocations(k+1) + end do + + do iCell = 1, nCellsSolve + bottomDepth(iCell) = config_tidal_boundary_left_bottom_depth & + + (yCell(iCell) - yMin) / (yMax - yMin) * & + (config_tidal_boundary_right_bottom_depth - config_tidal_boundary_left_bottom_depth) + + if (.not. config_alter_ICs_for_pbcs) then + do k = 1,nVertLevels + if (refBottomDepth(k) > bottomDepth(iCell)) then + bottomDepth(iCell) = refBottomDepth(k) + exit + end if + end do + end if + + ! Determine maxLevelCell based on bottomDepth and refBottomDepth + ! Also set botomDepth based on refBottomDepth, since + ! above bottomDepth was set with continuous analytical functions, + ! and needs to be discrete + maxLevelCell(iCell) = nVertLevels + if (nVertLevels > 1) then + do k = 1, nVertLevels + if (bottomDepth(iCell) < refBottomDepth(k)) then + maxLevelCell(iCell) = k-1 + bottomDepth(iCell) = refBottomDepth(k-1) + exit + end if + end do + end if + + end do + + ! Set tidal boundary mask + do iCell = 1, nCellsSolve + tidalInputMask(iCell) = 0.0_RKIND + if (yCell(iCell) > (25.0e3 - dcEdgeMinGlobal/2.0_RKIND)) then + tidalInputMask(iCell) = 1.0_RKIND + ! spread it over multiple cells + !if (yCell(iCell) > (25.0e3 - 3*dcEdgeMinGlobal)) then + ! tidalInputMask(iCell) = exp(-((yCell(iCell)-25.0e3)/dcEdgeMinGlobal)**2.0) + end if + end do + + do iCell = 1, nCellsSolve + ssh(iCell) = 0.0_RKIND + end do + + do iCell = 1, nCellsSolve + ! Set temperature + if ( associated(activeTracers) ) then + do k = 1, maxLevelCell(iCell) + if(yCell(iCell) < yMinGlobal + plugWidth) then + activeTracers(index_temperature, k, iCell) = config_tidal_boundary_plug_temperature + else + activeTracers(index_temperature, k, iCell) = config_tidal_boundary_domain_temperature + end if + end do + ! Set debug tracer to match plug onshore + if ( associated(debugTracers) ) then + do k = 1, maxLevelCell(iCell) + if(yCell(iCell) < (yMinGlobal + plugWidth)) then + debugTracers(index_tracer1, k, iCell) = 1.0_RKIND + else + debugTracers(index_tracer1, k, iCell) = 0.0_RKIND + end if + end do + end if + end if + end do + + do iCell = 1, nCellsSolve + ! Set layerThickness and restingThickness + call ocn_compute_layerThickness_zMid_from_bottomDepth(layerThickness(:,iCell), zMid(:,iCell), & + refBottomDepth, bottomDepth(iCell), maxLevelCell(iCell), nVertLevels, iErr, & + restingThickness=restingThickness(:,iCell), & + ssh=ssh(iCell)) + ! also computes zMid + do k = 1, maxLevelCell(iCell) + restingThickness(k, iCell) = layerThickness(k, iCell) + end do + + ! Set salinity + if ( associated(activeTracers) ) then + activeTracers(index_salinity, :, iCell) = config_tidal_boundary_salinity + end if + + end do + + ! Set vertCoordMovementWeights + vertCoordMovementWeights(:) = 1.0_RKIND + + block_ptr => block_ptr % next + end do + + deallocate(interfaceLocations) + + !-------------------------------------------------------------------- + + end subroutine ocn_init_setup_tidal_boundary!}}} + + +!*********************************************************************** +! +! routine ocn_init_validate_tidal_boundary +! +!> \brief Validation for Wetting and drying by tidal boundary test case +!> \author Phillip Wolfram +!> \date 04/05/2019 +!> \details +!> This routine validates the configuration options for the tidal_boundary test case. +! +!----------------------------------------------------------------------- + subroutine ocn_init_validate_tidal_boundary(configPool, packagePool, iocontext, iErr)!{{{ + + !-------------------------------------------------------------------- + + type (mpas_pool_type), intent(inout) :: configPool, packagePool + type (mpas_io_context_type), intent(inout) :: iocontext + + integer, intent(out) :: iErr + + character (len=StrKIND), pointer :: config_init_configuration + integer, pointer :: config_tidal_boundary_vert_levels, config_vert_levels + + iErr = 0 + + call mpas_pool_get_config(configPool, 'config_init_configuration', config_init_configuration) + + if(config_init_configuration .ne. trim('tidal_boundary')) return + + call mpas_pool_get_config(configPool, 'config_vert_levels', config_vert_levels) + call mpas_pool_get_config(configPool, 'config_tidal_boundary_vert_levels', config_tidal_boundary_vert_levels) + + if(config_vert_levels <= 0 .and. config_tidal_boundary_vert_levels > 0) then + config_vert_levels = config_tidal_boundary_vert_levels + !call mpas_log_write( 'Using value of $i', intArgs=(/ config_vert_levels /)) + else if(config_vert_levels <= 0) then + call mpas_log_write( 'Validation failed for tidal_boundary test case. ' & + // 'Not given a usable value for vertical levels.', MPAS_LOG_CRIT) + iErr = 1 + end if + + !-------------------------------------------------------------------- + + end subroutine ocn_init_validate_tidal_boundary!}}} + + +!*********************************************************************** + +end module ocn_init_tidal_boundary + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/src/core_ocean/shared/Makefile b/src/core_ocean/shared/Makefile index 94a2b6409b..dffdbdc5de 100644 --- a/src/core_ocean/shared/Makefile +++ b/src/core_ocean/shared/Makefile @@ -61,6 +61,7 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_surface_land_ice_fluxes.o \ mpas_ocn_effective_density_in_land_ice.o \ mpas_ocn_frazil_forcing.o \ + mpas_ocn_tidal_forcing.o \ mpas_ocn_forcing_restoring.o \ mpas_ocn_time_average_coupled.o \ mpas_ocn_sea_ice.o \ @@ -71,7 +72,7 @@ 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_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_diagnostics.o mpas_ocn_diagnostics_routines.o: mpas_ocn_constants.o @@ -167,6 +168,8 @@ mpas_ocn_surface_land_ice_fluxes.o: mpas_ocn_constants.o mpas_ocn_equation_of_st mpas_ocn_frazil_forcing.o: mpas_ocn_constants.o mpas_ocn_equation_of_state.o +mpas_ocn_tidal_forcing.o: mpas_ocn_constants.o mpas_ocn_equation_of_state.o + mpas_ocn_effective_density_in_land_ice.o: mpas_ocn_constants.o mpas_ocn_forcing_restoring.o: mpas_ocn_constants.o diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 8a23c7252a..87b3abefd4 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -31,6 +31,7 @@ module ocn_tendency use ocn_surface_bulk_forcing use ocn_surface_land_ice_fluxes use ocn_frazil_forcing + use ocn_tidal_forcing use ocn_tracer_hmix use ocn_high_freq_thickness_hmix_del2 @@ -185,6 +186,11 @@ subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{ ! call ocn_frazil_forcing_layer_thickness(meshPool, forcingPool, tend_layerThickness, err) + ! + ! tidal boundary condition flux + ! + call ocn_tidal_forcing_layer_thickness(meshPool, forcingPool, tend_layerThickness, err) + call mpas_timer_stop("ocn_tend_thick") end subroutine ocn_tend_thick!}}} diff --git a/src/core_ocean/shared/mpas_ocn_tidal_forcing.F b/src/core_ocean/shared/mpas_ocn_tidal_forcing.F new file mode 100644 index 0000000000..7602e76938 --- /dev/null +++ b/src/core_ocean/shared/mpas_ocn_tidal_forcing.F @@ -0,0 +1,309 @@ +! 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_tidal_forcing +! +!> \brief MPAS ocean tidal forcing module +!> \author Phillip J. Wolfram +!> \date 04/05/2019 +!> \details +!> This module contains routines for the tidal forcing. +! +!----------------------------------------------------------------------- + +module ocn_tidal_forcing + + use mpas_kind_types + use mpas_constants + use mpas_derived_types + use mpas_pool_routines + use mpas_timekeeping + use mpas_timer + use ocn_constants + use ocn_equation_of_state + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_tidal_forcing_build_array, & + ocn_tidal_forcing_layer_thickness, & + ocn_tidal_forcing_init + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + + logical :: tidalFluxOn + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_tidal_forcing_layer_thickness +! +!> \brief Add tendency due to tidal forcing +!> \author Phillip J. Wolfram +!> \date 04/05/2019 +!> \details +!> This routine adds a tendency to layer thickness due to tidal forcing +! +!----------------------------------------------------------------------- + + subroutine ocn_tidal_forcing_layer_thickness(meshPool, forcingPool, layerThicknessTend, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information + type (mpas_pool_type), intent(in) :: forcingPool !< Input: Forcing information + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + real (kind=RKIND), intent(inout), dimension(:,:) :: layerThicknessTend + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: Error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, k, nCells + integer, dimension(:), pointer :: nCellsArray + integer, dimension(:), pointer :: maxLevelCell + real (kind=RKIND), dimension(:,:), pointer :: tidalLayerThicknessTendency + + err = 0 + + if ( .not. tidalfluxOn ) return + + call mpas_timer_start("tidal thickness tendency") + + call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(forcingPool, 'tidalLayerThicknessTendency', & + tidalLayerThicknessTendency) + + ! frazil fields are needed only over 0 and 1 halos + nCells = nCellsArray( 2 ) + + ! Build surface fluxes at cell centers + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k = 1, maxLevelCell(iCell) + layerThicknessTend(k,iCell) = layerThicknessTend(k,iCell) + & + tidalLayerThicknessTendency(k,iCell) + + end do + end do + !$omp end do + + call mpas_timer_stop("tidal thickness tendency") + + end subroutine ocn_tidal_forcing_layer_thickness !}}} + + +!*********************************************************************** +! +! routine ocn_tidal_forcing_build_array +! +!> \brief Computes tidal forcing functions. +!> \author Phillip J. Wolfram +!> \date 04/05/2019 +!> \details +!> ocn_tidal_forcing_build_array computes the tendencies to layer thickness +!> due to tidal forcing of layer thickness +!> +!> this routine should be call at the beginning of whatever time stepping method is utilized +!> and the tendencies should be retieved when building up the RHS of the thickess, temperature +!> and salinity equations. +!> +!> this routine is only applicable to the thickness fields +! +!----------------------------------------------------------------------- + + subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, diagnosticsPool, statePool, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), pointer, intent(in) :: meshPool !< Input: Mesh information + type (mpas_pool_type), pointer, intent(in) :: diagnosticsPool !< Input: Diagnostic information + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + type (domain_type), intent(inout) :: domain + type (mpas_pool_type), pointer, intent(inout) :: statePool !< Input: State information + type (mpas_pool_type), pointer, intent(inout) :: forcingPool !< Input: Forcing information + integer, intent(inout) :: err !< Error flag + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), pointer :: tidalLayerThicknessTendency + real (kind=RKIND), dimension(:), pointer :: tidalInputMask + real (kind=RKIND), pointer :: config_tidal_forcing_monochromatic_amp + real (kind=RKIND), pointer :: config_tidal_forcing_monochromatic_period + real (kind=RKIND), pointer :: config_tidal_forcing_monochromatic_baseline + real (kind=RKIND), pointer :: config_use_tidal_forcing_tau + character (len=StrKIND), pointer :: config_tidal_forcing_type + + integer :: iCell, k, nCells + integer, dimension(:), pointer :: nCellsArray + integer, pointer, dimension(:) :: maxLevelCell + integer, pointer :: nVertLevels + real (kind=RKIND) :: dt + + type (MPAS_time_type) :: currentTime + real (kind=RKIND), pointer, dimension(:) :: ssh + real (kind=RKIND), pointer, dimension(:,:) :: zMid + real (kind=RKIND), pointer, dimension(:,:) :: layerThickness + + ! local variables + real (kind=RKIND) :: totalDepth, tidalHeight + + character (len=StrKIND), pointer :: simulationStartTime, xtime + real (kind=RKIND), pointer :: daysSinceStartOfSim + type (MPAS_Time_type) :: startTime, xtime_timeType, simulationStartTime_timeType + + if ( .not. tidalfluxOn ) return + + call mpas_pool_get_array(statePool, 'ssh', ssh, 1) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) + call mpas_pool_get_array(forcingPool, 'tidalLayerThicknessTendency', & + tidalLayerThicknessTendency) + call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask) + + ! compute time since start of simulation, in days + call mpas_pool_get_array(diagnosticsPool, 'daysSinceStartOfSim', daysSinceStartOfSim) + + call mpas_pool_get_config(ocnConfigs, 'config_tidal_forcing_type', config_tidal_forcing_type) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_forcing_monochromatic_amp', & + config_tidal_forcing_monochromatic_amp) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_forcing_monochromatic_period', & + config_tidal_forcing_monochromatic_period) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_forcing_monochromatic_baseline', & + config_tidal_forcing_monochromatic_baseline) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_forcing_tau', & + config_use_tidal_forcing_tau) + + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + + ! tidal fields are needed only over 0 and 1 halos + nCells = nCellsArray( 2 ) + + ! loop over all columns + do iCell = 1, nCells + + ! compute the tidalHeight + if (trim(config_tidal_forcing_type) == 'monochromatic') then + tidalHeight = config_tidal_forcing_monochromatic_amp * & + SIN(2.0_RKIND*pii/config_tidal_forcing_monochromatic_period * daysSinceStartOfSim) - & + config_tidal_forcing_monochromatic_baseline + else + tidalHeight = 0.0_RKIND + end if + + ! compute total depth for relative thickness contribution + totalDepth = 0.0_RKIND + do k = 1, maxLevelCell(iCell) + totalDepth = totalDepth + layerThickness(k,iCell) + end do + + ! distribute tidal forcing tendency fractionally over water column + do k = 1, maxLevelCell(iCell) + tidalLayerThicknessTendency(k,iCell) = tidalInputMask(iCell) / config_use_tidal_forcing_tau & + * (layerThickness(k,iCell)/totalDepth) * (tidalHeight - ssh(iCell)) + end do + enddo ! do iCell = 1, nCells + + end subroutine ocn_tidal_forcing_build_array !}}} + + +!*********************************************************************** +! +! routine ocn_tidal_forcing_init +! +!> \brief Initializes ocean tidal forcing module. +!> \author Phillip J. Wolfram +!> \date 04/05/2019 +!> \details +!> This routine initializes the ocean tidal forcing module and variables. +! +!----------------------------------------------------------------------- + + subroutine ocn_tidal_forcing_init(err)!{{{ + + integer, intent(out) :: err !< Output: error flag + logical, pointer :: config_use_tidal_forcing + + err = 0 + + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_forcing', config_use_tidal_forcing) + + tidalFluxOn = .false. + + if(config_use_tidal_forcing) then + tidalFluxOn = .true. + end if + + end subroutine ocn_tidal_forcing_init!}}} + + +!*********************************************************************** + +end module ocn_tidal_forcing + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/testing_and_setup/compass/ocean/surface_waves/default/1km/config_analysis.xml b/testing_and_setup/compass/ocean/surface_waves/default/1km/config_analysis.xml new file mode 100644 index 0000000000..40e4aba8f1 --- /dev/null +++ b/testing_and_setup/compass/ocean/surface_waves/default/1km/config_analysis.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/surface_waves/default/1km/config_driver.xml b/testing_and_setup/compass/ocean/surface_waves/default/1km/config_driver.xml new file mode 100644 index 0000000000..f7ea893e5e --- /dev/null +++ b/testing_and_setup/compass/ocean/surface_waves/default/1km/config_driver.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/surface_waves/default/1km/config_forward.xml b/testing_and_setup/compass/ocean/surface_waves/default/1km/config_forward.xml new file mode 100644 index 0000000000..a6a1af1daa --- /dev/null +++ b/testing_and_setup/compass/ocean/surface_waves/default/1km/config_forward.xml @@ -0,0 +1,70 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + init.nc + + + init.nc + + + input + initial_only + forcing.nc + + + + + + output + truncate + output.nc + 0000-00-00_01:00:00 + + + + + + + + + + + + + + + + init.nc + + + 4 + + + + diff --git a/testing_and_setup/compass/ocean/surface_waves/default/1km/config_init1.xml b/testing_and_setup/compass/ocean/surface_waves/default/1km/config_init1.xml new file mode 100644 index 0000000000..79cd79c92e --- /dev/null +++ b/testing_and_setup/compass/ocean/surface_waves/default/1km/config_init1.xml @@ -0,0 +1,101 @@ + + + + + + + + + + + + + + + + + + + + + + + + mesh.nc + + + output + forcing.nc + truncate + 0000_00:00:01 + + + + + + output + 0000_00:01:00 + truncate + ocean.nc + + + + + + + + + + + + + + + + + + + + + + + + + namelist.input.default + namelist.input + + + + -i + s/nx.*=.*/nx = 6,/g + namelist.input + + + -i + s/ny.*=.*/ny = 30,/g + namelist.input + + + -i + s/dc.*=.*/dc = 1000.0,/g + namelist.input + + + + + --keepx + -f + grid.nc + + + grid.nc + culled_mesh.nc + + + culled_mesh.nc + mesh.nc + + + + + diff --git a/testing_and_setup/compass/ocean/surface_waves/default/analysis/comparison.py b/testing_and_setup/compass/ocean/surface_waves/default/analysis/comparison.py new file mode 100755 index 0000000000..3226122ac2 --- /dev/null +++ b/testing_and_setup/compass/ocean/surface_waves/default/analysis/comparison.py @@ -0,0 +1,31 @@ +#!/usr/bin/env python +""" + +Tidal channel comparison betewen MPAS-O and analytical forcing result. + +Phillip J. Wolfram +04/12/2019 + +""" + +import numpy as np +import xarray as xr +import matplotlib.pyplot as plt + +# render statically by default +plt.switch_backend('agg') + +# analytical case +x = np.linspace(0,24,100) +y = np.sin(x*2*np.pi/24) +plt.plot(x,y, lw=3, color='black', label='analytical') + +# data from MPAS-O on boundary +ds = xr.open_mfdataset('output.nc') +ds.ssh.where(ds.tidalInputMask).mean('nCells').plot(marker='o', label='MPAS-O') + +plt.legend() +plt.ylabel('Tidal amplitude (m)') +plt.xlabel('Time (hrs)') + +plt.savefig('tidalcomparison.png')