diff --git a/src/core_ocean/Registry.xml b/src/core_ocean/Registry.xml index 92fd636456..ee4f4323f0 100644 --- a/src/core_ocean/Registry.xml +++ b/src/core_ocean/Registry.xml @@ -744,6 +744,10 @@ /> + + + + @@ -1557,6 +1570,7 @@ + @@ -1965,6 +1979,7 @@ + + diff --git a/src/core_ocean/driver/mpas_ocn_core_interface.F b/src/core_ocean/driver/mpas_ocn_core_interface.F index 082ece0852..9cd46cf5cb 100644 --- a/src/core_ocean/driver/mpas_ocn_core_interface.F +++ b/src/core_ocean/driver/mpas_ocn_core_interface.F @@ -113,6 +113,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: thicknessFilterActive logical, pointer :: splitTimeIntegratorActive logical, pointer :: windStressBulkPKGActive + logical, pointer :: variableBottomDragPKGActive logical, pointer :: tracerBudgetActive logical, pointer :: landIcePressurePKGActive logical, pointer :: landIceFluxesPKGActive @@ -159,6 +160,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ character (len=StrKIND), pointer :: config_pressure_gradient_type character (len=StrKIND), pointer :: config_sw_absorption_type + logical, pointer :: config_use_variable_drag logical, pointer :: config_use_bulk_wind_stress logical, pointer :: config_use_bulk_thickness_flux logical, pointer :: config_compute_active_tracer_budgets @@ -230,6 +232,14 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ windStressBulkPKGActive = .true. end if + ! + ! test for variable bottom drag of momentum, variableBottomDragPKG + call mpas_pool_get_package(packagePool, 'variableBottomDragPKGActive', variableBottomDragPKGActive) + call mpas_pool_get_config(configPool, 'config_use_variable_drag', config_use_variable_drag) + if ( config_use_variable_drag ) then + variableBottomDragPKGActive = .true. + end if + ! ! test for tracer budget ! diff --git a/src/core_ocean/mode_init/Registry_tidal_boundary.xml b/src/core_ocean/mode_init/Registry_tidal_boundary.xml index 5e2134cb7b..edf9e3aefb 100644 --- a/src/core_ocean/mode_init/Registry_tidal_boundary.xml +++ b/src/core_ocean/mode_init/Registry_tidal_boundary.xml @@ -51,4 +51,12 @@ description="Fraction of the domain the plug should take up initially. Only in the y direction." possible_values="Any real number between 0 and 1." /> + + 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 index 3a925500a9..8a025fe688 100644 --- a/src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F +++ b/src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F @@ -102,6 +102,7 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ 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 :: bottomDrag real (kind=RKIND), dimension(:), pointer :: ssh real (kind=RKIND), dimension(:,:), pointer :: layerThickness, restingThickness, zMid real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers, debugTracers @@ -112,6 +113,7 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ config_tidal_boundary_layer_type logical, pointer :: config_tidal_boundary_use_distances, config_alter_ICs_for_pbcs logical, pointer :: config_use_wetting_drying, config_tidal_start_dry + logical, pointer :: config_use_variable_drag 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, & @@ -119,7 +121,10 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ config_tidal_boundary_salinity, config_tidal_boundary_isopycnal_min_thickness, & config_tidal_boundary_left_value, & config_tidal_boundary_right_value, & - config_drying_min_cell_height + config_drying_min_cell_height, & + config_tidal_boundary_water_thickness, & + config_tidal_forcing_left_Cd_or_n, & + config_tidal_forcing_right_Cd_or_n integer, pointer :: config_tidal_boundary_min_vert_levels @@ -147,6 +152,9 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ 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_config(ocnConfigs, 'config_tidal_forcing_left_Cd_or_n', config_tidal_forcing_left_Cd_or_n) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_forcing_right_Cd_or_n', config_tidal_forcing_right_Cd_or_n) + call mpas_pool_get_config(ocnConfigs, 'config_use_variable_drag', config_use_variable_drag) call mpas_pool_get_config(ocnConfigs, 'config_tidal_start_dry', config_tidal_start_dry) call mpas_pool_get_config(ocnConfigs, 'config_use_wetting_drying', config_use_wetting_drying) @@ -244,6 +252,15 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ (config_tidal_boundary_right_bottom_depth - config_tidal_boundary_left_bottom_depth) end do + if (config_use_variable_drag) then + call mpas_pool_get_array(forcingPool, 'bottomDrag', bottomDrag) + do iCell = 1, nCellsSolve + bottomDrag(iCell) = config_tidal_forcing_left_Cd_or_n & + + (yCell(iCell) - yMin) / (yMax - yMin) * & + (config_tidal_forcing_right_Cd_or_n - config_tidal_forcing_left_Cd_or_n) + end do + end if + if (config_use_wetting_drying .and. config_tidal_start_dry .and. & trim(config_tidal_boundary_layer_type) == 'zstar') then do iCell = 1, nCellsSolve diff --git a/src/core_ocean/shared/mpas_ocn_vmix.F b/src/core_ocean/shared/mpas_ocn_vmix.F index d928363293..bdabf8b663 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix.F @@ -26,6 +26,7 @@ module ocn_vmix use mpas_pool_routines use mpas_timer + use mpas_constants use ocn_constants use ocn_vmix_coefs_const use ocn_vmix_coefs_tanh @@ -481,6 +482,322 @@ subroutine ocn_vel_vmix_tend_implicit(meshPool, dt, kineticEnergyCell, vertViscT end subroutine ocn_vel_vmix_tend_implicit!}}} +!*********************************************************************** +! +! routine ocn_vel_vmix_tend_implicit_variable +! +!> \brief Computes tendencies for implicit momentum vertical mixing +!> with spatially-variable bottom drag +!> \author Mark Petersen, Phillip J. Wolfram +!> \date September 2011, September 2019 +!> \details +!> This routine computes the tendencies for implicit vertical mixing for momentum +!> using computed coefficients from spatially-variable bottom drag. +!> Except for bottom drag coefficient, routine should be identifcal to +!> ocn_vel_vmix_tend_implicit above. +! +!----------------------------------------------------------------------- + + subroutine ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, dt, kineticEnergyCell, & !{{{ + vertViscTopOfEdge, layerThickness, & + layerThicknessEdge, normalVelocity, err) + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: & + meshPool !< Input: mesh information + + real (kind=RKIND), dimension(:,:), intent(in) :: & + kineticEnergyCell !< Input: kinetic energy at cell + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertViscTopOfEdge !< Input: vertical mixing coefficients + + real (kind=RKIND), intent(in) :: & + dt !< Input: time step + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< Input: thickness at cell center + + real (kind=RKIND), dimension(:), intent(in) :: & + bottomDrag !< Input: bottomDrag at cell centeres + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + normalVelocity !< Input: velocity + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + layerThicknessEdge !< Input: thickness at edge + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iEdge, k, cell1, cell2, N, nEdges + integer, pointer :: nVertLevels + real (kind=RKIND) :: implicitCd + integer, dimension(:), pointer :: nEdgesArray + + integer, dimension(:), pointer :: maxLevelEdgeTop + + integer, dimension(:,:), pointer :: cellsOnEdge + + real (kind=RKIND), dimension(:), allocatable :: A, B, C, velTemp + + err = 0 + + if(.not.velVmixOn) return + + call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + + call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + + nEdges = nEdgesArray( 1 ) + + allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),velTemp(nVertLevels)) + A(1)=0.0_RKIND + + !$omp do schedule(runtime) + do iEdge = 1, nEdges + N = maxLevelEdgeTop(iEdge) + if (N .gt. 0) then + + ! Compute A(k), B(k), C(k) + ! layerThicknessEdge is computed in compute_solve_diag, and is not available yet, + ! so recompute layerThicknessEdge here. + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, N + layerThicknessEdge(k,iEdge) = 0.5_RKIND * (layerThickness(k,cell1) + layerThickness(k,cell2)) + end do + + ! average cell-based implicit bottom drag to edges + implicitCd = 0.5_RKIND*(bottomDrag(cell1) + bottomDrag(cell2)) + + ! A is lower diagonal term + do k = 2, N + A(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k,iEdge) & + / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) & + / layerThicknessEdge(k,iEdge) + enddo + + ! C is upper diagonal term + do k = 1, N-1 + C(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k+1,iEdge) & + / (layerThicknessEdge(k,iEdge) + layerThicknessEdge(k+1,iEdge)) & + / layerThicknessEdge(k,iEdge) + enddo + + ! B is diagonal term + B(1) = 1.0_RKIND - C(1) + do k = 2, N-1 + B(k) = 1.0_RKIND - A(k) - C(k) + enddo + + ! Apply bottom drag boundary condition on the viscous term + ! second line uses sqrt(2.0*kineticEnergyEdge(k,iEdge)) + ! use implicitCd from spatially variable bottom drag + B(N) = 1.0_RKIND - A(N) + dt*implicitCd & + * sqrt(kineticEnergyCell(k,cell1) + kineticEnergyCell(k,cell2)) / layerThicknessEdge(k,iEdge) + + call tridiagonal_solve(A(2:N),B,C(1:N-1),normalVelocity(:,iEdge),velTemp,N) + + normalVelocity(1:N,iEdge) = velTemp(1:N) + normalVelocity(N+1:nVertLevels,iEdge) = 0.0_RKIND + + end if + end do + !$omp end do + + deallocate(A,B,C,velTemp) + + !-------------------------------------------------------------------- + + end subroutine ocn_vel_vmix_tend_implicit_spatially_variable!}}} + + +!*********************************************************************** +! +! routine ocn_vel_vmix_tend_implicit_variable_mannings +! +!> \brief Computes tendencies for implicit momentum vertical mixing +!> with spatially-variable bottom drag using Mannings friction +!> \author Mark Petersen, Phillip J. Wolfram +!> \date September 2011, September 2019 +!> \details +!> This routine computes the tendencies for implicit vertical mixing for momentum +!> using computed coefficients from spatially-variable bottom drag. +!> Except for bottom drag coefficient, routine should be identifcal to +!> ocn_vel_vmix_tend_implicit above. Cd uses Mannings' n values for the +!> Cd=g*n^2*h^{-1/3}. +! +!----------------------------------------------------------------------- + + subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, bottomDrag, dt, kineticEnergyCell, & !{{{ + vertViscTopOfEdge, layerThickness, & + layerThicknessEdge, normalVelocity, ssh, bottomDepth, err) + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: & + meshPool !< Input: mesh information + + real (kind=RKIND), dimension(:,:), intent(in) :: & + kineticEnergyCell !< Input: kinetic energy at cell + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertViscTopOfEdge !< Input: vertical mixing coefficients + + real (kind=RKIND), intent(in) :: & + dt !< Input: time step + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< Input: thickness at cell center + + real (kind=RKIND), dimension(:), intent(in) :: & + bottomDrag !< Input: bottomDrag at cell centeres + + real (kind=RKIND), dimension(:), pointer :: ssh + + real (kind=RKIND), dimension(:), pointer :: bottomDepth + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + normalVelocity !< Input: velocity + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + layerThicknessEdge !< Input: thickness at edge + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iEdge, k, cell1, cell2, N, nEdges + integer, pointer :: nVertLevels + real (kind=RKIND) :: implicitCd + integer, dimension(:), pointer :: nEdgesArray + + integer, dimension(:), pointer :: maxLevelEdgeTop + + integer, dimension(:,:), pointer :: cellsOnEdge + + real (kind=RKIND), dimension(:), allocatable :: A, B, C, velTemp + + err = 0 + + if(.not.velVmixOn) return + + call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + + call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + + nEdges = nEdgesArray( 1 ) + + allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),velTemp(nVertLevels)) + A(1)=0.0_RKIND + + !$omp do schedule(runtime) + do iEdge = 1, nEdges + N = maxLevelEdgeTop(iEdge) + if (N .gt. 0) then + + ! Compute A(k), B(k), C(k) + ! layerThicknessEdge is computed in compute_solve_diag, and is not available yet, + ! so recompute layerThicknessEdge here. + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, N + layerThicknessEdge(k,iEdge) = 0.5_RKIND * (layerThickness(k,cell1) + layerThickness(k,cell2)) + end do + + ! average cell-based implicit bottom drag to edges and convert Mannings n to Cd + implicitCd = gravity*(0.5_RKIND*(bottomDrag(cell1) + bottomDrag(cell2)))**2.0 * & + 0.5_RKIND * (ssh(cell1) + ssh(cell2) + bottomDepth(cell1) + bottomDepth(cell2))**(-1.0_RKIND/3.0_RKIND) + + ! A is lower diagonal term + do k = 2, N + A(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k,iEdge) & + / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) & + / layerThicknessEdge(k,iEdge) + enddo + + ! C is upper diagonal term + do k = 1, N-1 + C(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k+1,iEdge) & + / (layerThicknessEdge(k,iEdge) + layerThicknessEdge(k+1,iEdge)) & + / layerThicknessEdge(k,iEdge) + enddo + + ! B is diagonal term + B(1) = 1.0_RKIND - C(1) + do k = 2, N-1 + B(k) = 1.0_RKIND - A(k) - C(k) + enddo + + ! Apply bottom drag boundary condition on the viscous term + ! second line uses sqrt(2.0*kineticEnergyEdge(k,iEdge)) + ! use implicitCd from spatially variable bottom drag + B(N) = 1.0_RKIND - A(N) + dt*implicitCd & + * sqrt(kineticEnergyCell(k,cell1) + kineticEnergyCell(k,cell2)) / layerThicknessEdge(k,iEdge) + + call tridiagonal_solve(A(2:N),B,C(1:N-1),normalVelocity(:,iEdge),velTemp,N) + + normalVelocity(1:N,iEdge) = velTemp(1:N) + normalVelocity(N+1:nVertLevels,iEdge) = 0.0_RKIND + + end if + end do + !$omp end do + + deallocate(A,B,C,velTemp) + + !-------------------------------------------------------------------- + + end subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings!}}} + + !*********************************************************************** ! ! routine ocn_tracer_vmix_tend_implicit @@ -644,6 +961,7 @@ subroutine ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPo integer :: iCell, timeLevel, k, cell1, cell2, iEdge, nCells, nEdges integer, dimension(:), pointer :: nCellsArray, nEdgesArray + real (kind=RKIND), dimension(:), pointer :: bottomDrag, ssh, bottomDepth ! needed for depth-variable computation real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThickness, layerThicknessEdge, vertViscTopOfEdge, & vertDiffTopOfCell, kineticEnergyCell real (kind=RKIND), dimension(:,:), pointer :: vertViscTopOfCell, nonLocalSurfaceTracerFlux, tracerGroupSurfaceFlux @@ -656,6 +974,7 @@ subroutine ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPo logical, pointer :: config_cvmix_kpp_nonlocal_with_implicit_mix logical, pointer :: config_Rayleigh_friction, config_Rayleigh_bottom_friction, & config_Rayleigh_damping_depth_variable + logical, pointer :: config_use_implicit_bottom_drag_variable, config_use_implicit_bottom_drag_variable_mannings type (mpas_pool_iterator_type) :: groupItr @@ -685,8 +1004,15 @@ subroutine ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPo config_Rayleigh_bottom_friction) call mpas_pool_get_config(ocnConfigs, 'config_Rayleigh_damping_depth_variable', & config_Rayleigh_damping_depth_variable) + call mpas_pool_get_config(ocnConfigs, 'config_Rayleigh_damping_depth_variable', & + config_Rayleigh_damping_depth_variable) + call mpas_pool_get_config(ocnConfigs, 'config_use_implicit_bottom_drag_variable', & + config_use_implicit_bottom_drag_variable) + call mpas_pool_get_config(ocnConfigs, 'config_use_implicit_bottom_drag_variable_mannings', & + config_use_implicit_bottom_drag_variable_mannings) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevel) call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevel) + call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel) call mpas_pool_get_array(diagnosticsPool, 'kineticEnergyCell', kineticEnergyCell) call mpas_pool_get_array(diagnosticsPool, 'layerThicknessEdge', layerThicknessEdge) @@ -699,11 +1025,14 @@ subroutine ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPo call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) + call mpas_pool_get_array(forcingPool, 'bottomDrag', bottomDrag) + call mpas_timer_start('vmix coefs', .false.) call ocn_vmix_coefs(meshPool, statePool, forcingPool, diagnosticsPool, scratchPool, err, timeLevel) call mpas_timer_stop('vmix coefs') @@ -734,12 +1063,25 @@ subroutine ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPo call mpas_timer_start('vmix solve momentum', .false.) if (.not. (config_Rayleigh_friction .or. & config_Rayleigh_bottom_friction .or. & - config_Rayleigh_damping_depth_variable)) then + config_Rayleigh_damping_depth_variable .or. & + config_use_implicit_bottom_drag_variable .or. & + config_use_implicit_bottom_drag_variable_mannings)) then call ocn_vel_vmix_tend_implicit(meshPool, dt, kineticEnergyCell, & vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, err) else - call ocn_vel_vmix_tend_implicit_rayleigh(meshPool, dt, kineticEnergyCell, & - vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, err) + if (config_use_implicit_bottom_drag_variable) then + call ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, dt, kineticEnergyCell, & + vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, err) + else if (config_use_implicit_bottom_drag_variable_mannings) then + ! update bottomDrag via Cd=g*n^2*h^(-1/3) + call ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, bottomDrag, & + dt, kineticEnergyCell, & + vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, & + ssh, bottomDepth, err) + else if (config_Rayleigh_damping_depth_variable) then + call ocn_vel_vmix_tend_implicit_rayleigh(meshPool, dt, kineticEnergyCell, & + vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, err) + end if end if call mpas_timer_stop('vmix solve momentum') diff --git a/testing_and_setup/compass/ocean/drying_slope/forward_template.xml b/testing_and_setup/compass/ocean/drying_slope/forward_template.xml index 29edbf55df..031256db38 100644 --- a/testing_and_setup/compass/ocean/drying_slope/forward_template.xml +++ b/testing_and_setup/compass/ocean/drying_slope/forward_template.xml @@ -23,6 +23,7 @@ + @@ -44,6 +45,7 @@ forcing.nc + diff --git a/testing_and_setup/compass/ocean/drying_slope/init_template.xml b/testing_and_setup/compass/ocean/drying_slope/init_template.xml index ed4c8e56f2..41e3b30ada 100644 --- a/testing_and_setup/compass/ocean/drying_slope/init_template.xml +++ b/testing_and_setup/compass/ocean/drying_slope/init_template.xml @@ -11,6 +11,7 @@ + @@ -24,6 +25,7 @@ 0000_00:00:01 + diff --git a/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/analysis/comparison.py b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/analysis/comparison.py new file mode 120000 index 0000000000..b40167a8a1 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/analysis/comparison.py @@ -0,0 +1 @@ +../comparison.py \ No newline at end of file diff --git a/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/comparison.py b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/comparison.py new file mode 100755 index 0000000000..0a4cb15021 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/comparison.py @@ -0,0 +1,148 @@ +#!/usr/bin/env python +""" + +Drying slope comparison betewen MPAS-O, analytical, and ROMS results from + +Warner, J. C., Defne, Z., Haas, K., & Arango, H. G. (2013). A wetting and +drying scheme for ROMS. Computers & geosciences, 58, 54-61. + +Phillip J. Wolfram and Zhendong Cao +04/30/2019 + +""" + +import numpy as np +import xarray as xr +import matplotlib.pyplot as plt +import pandas as pd + +# render statically by default +plt.switch_backend('agg') + +def setup_fig(): + fig, ax = plt.subplots(nrows=4,ncols=1, sharex=True, sharey=True, figsize=(6,8)) + fig.text(0.02, 0.5, 'Channel depth (m)', va='center', rotation='vertical') + #fig.text(0.5, 0.02, 'Along channel distance (km)', ha='center') + +def setup_subplot(): + plt.xlim(0,25) + plt.ylim(-1, 11) + ax = plt.gca() + ax.invert_yaxis() + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + + x = np.linspace(0,25,100) + y = 10.0/25.0*x + plt.plot(x, y, 'k-', lw=3) + + +def upper_plot(): + plt.subplot(2,1,1) + plt.gca().set_xticklabels([]) + setup_subplot() + + +def lower_plot(): + plt.subplot(2,1,2) + setup_subplot() + + +def plot_data(rval='0.0025', dtime='0.05', datatype='analytical', *args, **kwargs): + datafile = 'data/r' + rval + 'd' + dtime + '-' + datatype + '.csv' + data = pd.read_csv(datafile, header=None) + measured=plt.scatter(data[0], data[1], *args, **kwargs) + + +def plot_datasets(rval, times, fileno, plotdata=True): + for ii, dtime in enumerate(times): + if plotdata: + plot_data(rval=rval, dtime = dtime, datatype = 'analytical', + marker = '.', color = 'b', label='analytical') + plot_data(rval=rval, dtime = dtime, datatype = 'roms', + marker = '.', color = 'g', label='ROMS') + plot_MPASO([dtime], fileno, 'k-', lw=0.5, label='MPAS-O') + + if ii == 0: + plt.legend(frameon=False, loc='lower left') + place_time_labels(times) + plt.text(0.5, 5, 'Cd = ' + str(rval)) + + +def place_time_labels(times): + locs = [9.3, 7.2, 4.2, 2.2, 1.2, 0.2] + for atime, ay in zip(times, locs): + plt.text(25.2, ay, atime + ' days', size=8) + +def plot_MPASO(times, fileno, *args, **kwargs): + for atime in times: + # factor of 1e-16 needed to account for annoying round-off issue to get right time slices + plottime = int((float(atime)/0.2 + 1e-16)*24.0) + ds = xr.open_mfdataset('output'+ fileno + '.nc') + print('{} {} {}'.format(atime, plottime, ds.isel(Time=plottime).xtime.values)) + ds = ds.drop(np.setdiff1d([i for i in ds.variables], ['yCell','ssh'])) + ymean = ds.isel(Time=plottime).groupby('yCell').mean(dim=xr.ALL_DIMS) + x = ymean.yCell.values/1000.0 + y = ymean.ssh.values + #print('ymin={} ymax={}\n{}\n{}'.format(y.min(), y.max(),x, y)) + plt.plot(x, -y, *args, **kwargs) + ds.close() + + +def plot_tidal_forcing_comparison(): + # data from MPAS-O on boundary + for fileno in ['1','2']: + ds = xr.open_mfdataset('output' + fileno +'.nc') + ympas = ds.ssh.where(ds.tidalInputMask).mean('nCells').values + x = np.linspace(0, 1.0, len(ds.xtime))*12.0 + plt.plot(x, ympas, marker='o', label='MPAS-O forward' + fileno) + + # analytical case + x = np.linspace(0,12.0,100) + y = 10.0*np.sin(x*np.pi/12.0) - 10.0 + plt.plot(x, y, lw=3, color='black', label='analytical') + + plt.legend(frameon=False) + plt.ylabel('Tidal amplitude (m)') + plt.xlabel('Time (hrs)') + plt.suptitle('Tidal amplitude forcing (right side) for MPAS-O and analytical') + plt.savefig('tidalcomparison.png') + + +def main(): + ################ plot tidal forcing comparison ############### + #plot_tidal_forcing_comparison() + ############################################################## + + ################ plot drying front comparison ############### + times = ['0.50', '0.05', '0.40', '0.15', '0.30', '0.25'] + setup_fig() + + ############################################################## + plt.subplot(4,1,1) + plt.title('MPAS-O Drying slope comparison for variable Cd') + plt.gca().set_xticklabels([]) + setup_subplot() + plot_datasets(rval='1e-1 to 1e-2', times=times, fileno='4', plotdata=False) + + plt.subplot(4,1,2) + plt.gca().set_xticklabels([]) + setup_subplot() + plot_datasets(rval='1e-2', times=times, fileno='2', plotdata=False) + + plt.subplot(4,1,3) + plt.gca().set_xticklabels([]) + setup_subplot() + plot_datasets(rval='1e-2 to 1e-3', times=times, fileno='3', plotdata=False) + + # subplot Cd=1e-2 + plt.subplot(4,1,4) + setup_subplot() + plot_datasets(rval='1e-3', times=times, fileno='1', plotdata=False) + plt.xlabel('Along channel distance (km)') + + for outtype in ['.png','.pdf']: + plt.savefig('dryingslopecomparison' + outtype, tight_layout=True) + +if __name__ == "__main__": + main() diff --git a/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/config_analysis.xml b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/config_analysis.xml new file mode 100644 index 0000000000..43c1057ba4 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/config_analysis.xml @@ -0,0 +1,30 @@ + + + + + + + + + + + + + + output1.nc + vtk_output1 + allOnCells + maxEdges=0 + nVertLevels=0:10 + + + + output2.nc + vtk_output2 + allOnCells + maxEdges=0 + nVertLevels=0:10 + + + + diff --git a/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/config_driver.xml b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/config_driver.xml new file mode 100644 index 0000000000..e0545f3143 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/config_driver.xml @@ -0,0 +1,26 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/config_forward1.xml b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/config_forward1.xml new file mode 100644 index 0000000000..7d0b0cfebe --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/zstar_variableCd/1km/config_forward1.xml @@ -0,0 +1,25 @@ + + + + + + + +