From 2bfa4bce5aaecdaacc8606485b0c2fa164efb8f1 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 10 Jun 2021 15:13:34 -0800 Subject: [PATCH 01/15] Adding a pass_var to surface h This is required for the u,v sponges to be invariant of tiling. I don't know why, but the problem only showed up for me in a narrow channel in the Bering domain. --- src/parameterizations/vertical/MOM_ALE_sponge.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 5d9af389c9..84e91aaffc 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -1096,6 +1096,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (CS%id_sp_u_tendency > 0) then allocate(tmp_u(G%isdB:G%iedB,G%jsd:G%jed,nz));tmp_u(:,:,:)=0.0 endif + call pass_var(h(:,:,1),G%Domain) ! u points do c=1,CS%num_col_u I = CS%col_i_u(c) ; j = CS%col_j_u(c) From 7bbc3096e3b38334e7720602ff0f8b13fa79677e Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 16 Jun 2021 14:49:33 -0400 Subject: [PATCH 02/15] Replace MOM_control_struct pointers as locals The MOM_control_struct is declared and passed as a pointer to a CS (usually allocated anonymously) and treated as if it were a pointer, even though there is currently no real advantage to doing so. After the FMS update, the deallocation of this CS was causing a segmentation fault in the PGI compilers. While the underlying cause was never determined, it is likely due to some automated deallocation of the CS contents, whose addressing became scrambled. This problem can be resolved by moving all of the CS contents to stack, so that the contents are automatically removed upon exiting whatever function it was instantiated. Subsequent calls can reference the local (or parent) stack contents. --- .../drivers/FMS_cap/ocean_model_MOM.F90 | 4 +-- .../drivers/mct_cap/mom_ocean_model_mct.F90 | 4 +-- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 4 +-- config_src/drivers/solo_driver/MOM_driver.F90 | 3 +- src/core/MOM.F90 | 34 ++++++++----------- 5 files changed, 21 insertions(+), 28 deletions(-) diff --git a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 index c3e13329f2..50ea6c943d 100644 --- a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 +++ b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 @@ -195,8 +195,8 @@ module ocean_model_mod type(unit_scale_type), pointer :: & US => NULL() !< A pointer to a structure containing dimensional !! unit scaling factors. - type(MOM_control_struct), pointer :: & - MOM_CSp => NULL() !< A pointer to the MOM control structure + type(MOM_control_struct) :: MOM_CSp + !< MOM control structure type(ice_shelf_CS), pointer :: & Ice_shelf_CSp => NULL() !< A pointer to the control structure for the !! ice shelf model that couples with MOM6. This diff --git a/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 b/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 index 9b40a9e7b4..3bd0e1e28d 100644 --- a/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 +++ b/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 @@ -195,8 +195,8 @@ module MOM_ocean_model_mct !! about the vertical grid. type(unit_scale_type), pointer :: US => NULL() !< A pointer to a structure containing !! dimensional unit scaling factors. - type(MOM_control_struct), pointer :: & - MOM_CSp => NULL() !< A pointer to the MOM control structure + type(MOM_control_struct) :: MOM_CSp + !< MOM control structure type(ice_shelf_CS), pointer :: & Ice_shelf_CSp => NULL() !< A pointer to the control structure for the !! ice shelf model that couples with MOM6. This diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 1d5de0dd3e..7f91e15a69 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -197,8 +197,8 @@ module MOM_ocean_model_nuopc !! about the vertical grid. type(unit_scale_type), pointer :: US => NULL() !< A pointer to a structure containing !! dimensional unit scaling factors. - type(MOM_control_struct), pointer :: & - MOM_CSp => NULL() !< A pointer to the MOM control structure + type(MOM_control_struct) :: MOM_CSp + !< MOM control structure type(ice_shelf_CS), pointer :: & Ice_shelf_CSp => NULL() !< A pointer to the control structure for the !! ice shelf model that couples with MOM6. This diff --git a/config_src/drivers/solo_driver/MOM_driver.F90 b/config_src/drivers/solo_driver/MOM_driver.F90 index ebb953be93..7dfce01f68 100644 --- a/config_src/drivers/solo_driver/MOM_driver.F90 +++ b/config_src/drivers/solo_driver/MOM_driver.F90 @@ -180,8 +180,7 @@ program MOM_main ! and diffusion equation are read in from files stored from ! a previous integration of the prognostic model - type(MOM_control_struct), pointer :: MOM_CSp => NULL() - !> A pointer to the tracer flow control structure. + type(MOM_control_struct) :: MOM_CSp !> MOM control structure type(tracer_flow_control_CS), pointer :: & tracer_flow_CSp => NULL() !< A pointer to the tracer flow control structure type(surface_forcing_CS), pointer :: surface_forcing_CSp => NULL() diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 05c6fe6c43..91e2ea6afd 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -434,7 +434,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS type(surface), target, intent(inout) :: sfc_state !< surface ocean state type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type real, intent(in) :: time_int_in !< time interval covered by this run segment [s]. - type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM + type(MOM_control_struct), intent(inout), target :: CS !< control structure from initialize_MOM type(Wave_parameters_CS), & optional, pointer :: Waves !< An optional pointer to a wave property CS logical, optional, intent(in) :: do_dynamics !< Present and false, do not do updates due @@ -981,7 +981,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & real, intent(in) :: bbl_time_int !< time interval over which updates to the !! bottom boundary layer properties will apply [T ~> s], !! or zero not to update the properties. - type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM + type(MOM_control_struct), intent(inout), target :: CS !< control structure from initialize_MOM type(time_type), intent(in) :: Time_local !< End time of a segment, as a time type type(wave_parameters_CS), & optional, pointer :: Waves !< Container for wave related parameters; the @@ -1432,7 +1432,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS type(surface), intent(inout) :: sfc_state !< surface ocean state type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type real, intent(in) :: time_interval !< time interval - type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM + type(MOM_control_struct), intent(inout) :: CS !< control structure from initialize_MOM ! Local pointers type(ocean_grid_type), pointer :: G => NULL() ! Pointer to a structure containing @@ -1630,7 +1630,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar type(param_file_type), intent(out) :: param_file !< structure indicating parameter file to parse type(directories), intent(out) :: dirs !< structure with directory paths - type(MOM_control_struct), pointer :: CS !< pointer set in this routine to MOM control structure + type(MOM_control_struct), intent(inout), target :: CS !< pointer set in this routine to MOM control structure type(MOM_restart_CS), pointer :: restart_CSp !< pointer set in this routine to the !! restart control structure that will !! be used for MOM. @@ -1730,13 +1730,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & type(ocean_internal_state) :: MOM_internal_state character(len=200) :: area_varname, ice_shelf_file, inputdir, filename - if (associated(CS)) then - call MOM_error(WARNING, "initialize_MOM called with an associated "// & - "control structure.") - return - endif - allocate(CS) - CS%Time => Time id_clock_init = cpu_clock_id('Ocean Initialization', grain=CLOCK_SUBCOMPONENT) @@ -2818,7 +2811,7 @@ end subroutine initialize_MOM subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) type(time_type), intent(in) :: Time !< model time, used in this routine type(directories), intent(in) :: dirs !< structure with directory paths - type(MOM_control_struct), pointer :: CS !< pointer to MOM control structure + type(MOM_control_struct), intent(inout) :: CS !< MOM control structure type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control !! structure that will be used for MOM. ! Local variables @@ -3044,7 +3037,7 @@ end subroutine adjust_ssh_for_p_atm !! setting the appropriate fields in sfc_state. Unused fields !! are set to NULL or are unallocated. subroutine extract_surface_state(CS, sfc_state_in) - type(MOM_control_struct), pointer :: CS !< Master MOM control structure + type(MOM_control_struct), intent(inout), target :: CS !< Master MOM control structure type(surface), target, intent(inout) :: sfc_state_in !< transparent ocean surface state !! structure shared with the calling routine !! data in this structure is intent out. @@ -3471,7 +3464,7 @@ end subroutine rotate_initial_state !> Return true if all phases of step_MOM are at the same point in time. function MOM_state_is_synchronized(CS, adv_dyn) result(in_synch) - type(MOM_control_struct), pointer :: CS !< MOM control structure + type(MOM_control_struct), intent(inout) :: CS !< MOM control structure logical, optional, intent(in) :: adv_dyn !< If present and true, only check !! whether the advection is up-to-date with !! the dynamics. @@ -3492,7 +3485,7 @@ end function MOM_state_is_synchronized !> This subroutine offers access to values or pointers to other types from within !! the MOM_control_struct, allowing the MOM_control_struct to be opaque. subroutine get_MOM_state_elements(CS, G, GV, US, C_p, C_p_scaled, use_temp) - type(MOM_control_struct), pointer :: CS !< MOM control structure + type(MOM_control_struct), intent(inout), target :: CS !< MOM control structure type(ocean_grid_type), optional, pointer :: G !< structure containing metrics and grid info type(verticalGrid_type), optional, pointer :: GV !< structure containing vertical grid info type(unit_scale_type), optional, pointer :: US !< A dimensional unit scaling type @@ -3511,7 +3504,7 @@ end subroutine get_MOM_state_elements !> Find the global integrals of various quantities. subroutine get_ocean_stocks(CS, mass, heat, salt, on_PE_only) - type(MOM_control_struct), pointer :: CS !< MOM control structure + type(MOM_control_struct), intent(inout) :: CS !< MOM control structure real, optional, intent(out) :: heat !< The globally integrated integrated ocean heat [J]. real, optional, intent(out) :: salt !< The globally integrated integrated ocean salt [kg]. real, optional, intent(out) :: mass !< The globally integrated integrated ocean mass [kg]. @@ -3528,7 +3521,7 @@ end subroutine get_ocean_stocks !> End of ocean model, including memory deallocation subroutine MOM_end(CS) - type(MOM_control_struct), pointer :: CS !< MOM control structure + type(MOM_control_struct), intent(inout) :: CS !< MOM control structure call MOM_sum_output_end(CS%sum_output_CSp) @@ -3604,7 +3597,6 @@ subroutine MOM_end(CS) if (associated(CS%update_OBC_CSp)) call OBC_register_end(CS%update_OBC_CSp) call verticalGridEnd(CS%GV) - call unit_scaling_end(CS%US) call MOM_grid_end(CS%G) if (CS%debug .or. CS%G%symmetric) & @@ -3613,9 +3605,11 @@ subroutine MOM_end(CS) if (CS%rotate_index) & call deallocate_MOM_domain(CS%G%Domain) - call deallocate_MOM_domain(CS%G_in%domain) + ! The MPP domains may be needed by an external coupler, so use `cursory`. + ! TODO: This may create a domain memory leak, and needs investigation. + call deallocate_MOM_domain(CS%G_in%domain, cursory=.true.) - deallocate(CS) + call unit_scaling_end(CS%US) end subroutine MOM_end !> \namespace mom From 95a770ac76882b5c6da94e5ad3c3e5bdd986c7a7 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 22 Jun 2021 12:29:49 -0800 Subject: [PATCH 03/15] Matt's suggestion for not sponging at the grid edge. --- src/parameterizations/vertical/MOM_ALE_sponge.F90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 84e91aaffc..a84e45c1b0 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -16,7 +16,7 @@ module MOM_ALE_sponge use MOM_coms, only : sum_across_PEs use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field use MOM_diag_mediator, only : diag_ctrl -use MOM_domains, only : pass_var +use MOM_domains, only : pass_var use MOM_error_handler, only : MOM_error, FATAL, NOTE, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type @@ -1008,10 +1008,11 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,& answers_2018=CS%hor_regrid_answers_2018) - call pass_var(sp_val,G%Domain) + call pass_var(sp_val, G%Domain) + call pass_var(mask_z, G%Domain) do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB sp_val_u(I,j,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i+1,j,1:nz_data)) - mask_u(I,j,1:nz_data) = max(mask_z(i,j,1:nz_data),mask_z(i+1,j,1:nz_data)) + mask_u(I,j,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i+1,j,1:nz_data)) enddo ; enddo allocate( hsrc(nz_data) ) @@ -1055,10 +1056,11 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., & spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,& answers_2018=CS%hor_regrid_answers_2018) - call pass_var(sp_val,G%Domain) + call pass_var(sp_val, G%Domain) + call pass_var(mask_z, G%Domain) do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec sp_val_v(i,J,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i,j+1,1:nz_data)) - mask_v(i,J,1:nz_data) = max(mask_z(i,j,1:nz_data),mask_z(i,j+1,1:nz_data)) + mask_v(i,J,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i,j+1,1:nz_data)) enddo ; enddo !call pass_var(mask_z,G%Domain) allocate( hsrc(nz_data) ) From 53dfdc71dbb317317a393509c1c4250b8fb3a5a8 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 23 Jun 2021 18:53:06 -0800 Subject: [PATCH 04/15] Finish the masking out of edge in uv sponge. - Without this change, the edges don't reproduce on restart due to the h values outside being nonsense. --- .../vertical/MOM_ALE_sponge.F90 | 32 +++++++++---------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index d19de4ce46..e122452368 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -709,7 +709,6 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, integer :: nz_data !< the number of vertical levels in this input field character(len=256) :: mesg ! String for error messages ! Local variables for ALE remapping - real, dimension(:), allocatable :: tmpT1d real :: zTopOfCell, zBottomOfCell ! Heights [Z ~> m]. type(remapping_CS) :: remapCS ! Remapping parameters and work arrays @@ -841,7 +840,7 @@ subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename else CS%Ref_val_v%id = init_external_field(filename_v, fieldname_v) endif - fld_sz(1:4)=-1 fld_sz(1:4)=-1 + fld_sz(1:4)=-1 call get_external_field_info(CS%Ref_val_v%id, size=fld_sz) CS%Ref_val_v%nz_data = fld_sz(3) CS%Ref_val_v%num_tlevs = fld_sz(4) @@ -1027,23 +1026,23 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) enddo ; enddo allocate( hsrc(nz_data) ) - allocate( tmpT1d(nz_data) ) do c=1,CS%num_col_u ! c is an index for the next 3 lines but a multiplier for the rest of the loop ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i_u(c) ; j = CS%col_j_u(c) - CS%Ref_val_u%p(1:nz_data,c) = sp_val_u(i,j,1:nz_data) + if (mask_u(i,j,1) == 1.0) then + CS%Ref_val_u%p(1:nz_data,c) = sp_val_u(i,j,1:nz_data) + else + CS%Ref_val_u%p(1:nz_data,c) = 0.0 + endif ! Build the source grid - zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0; tmpT1d(:) = -99.9 + zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0 do k=1,nz_data if (mask_u(i,j,k) == 1.0) then zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(i,j) ) - tmpT1d(k) = sp_val_u(i,j,k) elseif (k>1) then zBottomOfCell = -G%bathyT(i,j) - tmpT1d(k) = tmpT1d(k-1) else ! This next block should only ever be reached over land - tmpT1d(k) = -99.9 endif hsrc(k) = zTopOfCell - zBottomOfCell if (hsrc(k)>0.) nPoints = nPoints + 1 @@ -1053,7 +1052,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(i,j) ) CS%Ref_val_u%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) enddo - deallocate(sp_val, sp_val_u, mask_u, mask_z, hsrc, tmpT1d) + deallocate(sp_val, sp_val_u, mask_u, mask_z, hsrc) nz_data = CS%Ref_val_v%nz_data allocate(sp_val( G%isd:G%ied,G%jsd:G%jed,1:nz_data)) allocate(sp_val_v(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data)) @@ -1075,23 +1074,23 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) enddo ; enddo !call pass_var(mask_z,G%Domain) allocate( hsrc(nz_data) ) - allocate( tmpT1d(nz_data) ) do c=1,CS%num_col_v ! c is an index for the next 3 lines but a multiplier for the rest of the loop ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i_v(c) ; j = CS%col_j_v(c) - CS%Ref_val_v%p(1:nz_data,c) = sp_val_v(i,j,1:nz_data) + if (mask_v(i,j,1) == 1.0) then + CS%Ref_val_v%p(1:nz_data,c) = sp_val_v(i,j,1:nz_data) + else + CS%Ref_val_v%p(1:nz_data,c) = 0.0 + endif ! Build the source grid - zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0; tmpT1d(:) = -99.9 + zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0 do k=1,nz_data if (mask_v(i,j,k) == 1.0) then zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(i,j) ) - tmpT1d(k) = sp_val_v(i,j,k) elseif (k>1) then zBottomOfCell = -G%bathyT(i,j) - tmpT1d(k) = tmpT1d(k-1) else ! This next block should only ever be reached over land - tmpT1d(k) = -99.9 endif hsrc(k) = zTopOfCell - zBottomOfCell if (hsrc(k)>0.) nPoints = nPoints + 1 @@ -1101,7 +1100,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(i,j) ) CS%Ref_val_v%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) enddo - deallocate(sp_val, sp_val_v, mask_v, mask_z, hsrc, tmpT1d) + deallocate(sp_val, sp_val_v, mask_v, mask_z, hsrc) endif call pass_var(h,G%Domain) @@ -1110,7 +1109,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (CS%id_sp_u_tendency > 0) then allocate(tmp_u(G%isdB:G%iedB,G%jsd:G%jed,nz));tmp_u(:,:,:)=0.0 endif - call pass_var(h(:,:,1),G%Domain) ! u points do c=1,CS%num_col_u I = CS%col_i_u(c) ; j = CS%col_j_u(c) From b7f2cea3ebc1ee52fdc5fdf41990c9e0971580f7 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Mon, 19 Jul 2021 13:29:21 -0800 Subject: [PATCH 05/15] Add location to ssh warning when dry. --- src/core/MOM_barotropic.F90 | 9 +++++++-- src/core/MOM_boundary_update.F90 | 3 --- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 51f9a5cb85..1fa0ae5bcd 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2356,8 +2356,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (GV%Boussinesq) then do j=js,je ; do i=is,ie - if (eta(i,j) < -GV%Z_to_H*G%bathyT(i,j)) & - call MOM_error(WARNING, "btstep: eta has dropped below bathyT.") + if (eta(i,j) < -GV%Z_to_H*G%bathyT(i,j)) then + ioff = G%idg_offset + i + joff = G%jdg_offset + j + write(mesg,"('btstep: eta has dropped below bathyT at ', i5, i5)") ioff, joff + call MOM_error(WARNING, trim(mesg)) +! call MOM_error(WARNING, "btstep: eta has dropped below bathyT.") + endif enddo ; enddo else do j=js,je ; do i=is,ie diff --git a/src/core/MOM_boundary_update.F90 b/src/core/MOM_boundary_update.F90 index 2e25af2460..53ca87ea2c 100644 --- a/src/core/MOM_boundary_update.F90 +++ b/src/core/MOM_boundary_update.F90 @@ -121,9 +121,6 @@ subroutine update_OBC_data(OBC, G, GV, US, tv, h, CS, Time) type(update_OBC_CS), pointer :: CS !< Control structure for OBCs type(time_type), intent(in) :: Time !< Model time -! Something here... with CS%file_OBC_CSp? -! if (CS%use_files) & -! call update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (CS%use_tidal_bay) & call tidal_bay_set_OBC_data(OBC, CS%tidal_bay_OBC_CSp, G, GV, h, Time) if (CS%use_Kelvin) & From 7e9fd3d84cc262b8c06a17d1fd10f0de2c5c96c1 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 18 Aug 2021 09:45:31 -0800 Subject: [PATCH 06/15] Tiny typo. --- src/parameterizations/lateral/MOM_internal_tides.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 8c08691675..68b0594be4 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -2297,7 +2297,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) "The default is 2pi/10 km, as in St.Laurent et al. 2002.", & units="m-1", default=8.e-4*atan(1.0), scale=US%L_to_m) call get_param(param_file, mdl, "KAPPA_H2_FACTOR", kappa_h2_factor, & - "A scaling factor for the roughness amplitude with n"//& + "A scaling factor for the roughness amplitude with "//& "INT_TIDE_DISSIPATION.", units="nondim", default=1.0) ! Allocate various arrays needed for loss rates From dabf8eaafa88e17ea0623cedd6b068a30a83032b Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 19 Aug 2021 15:54:28 -0800 Subject: [PATCH 07/15] Documented internal tidal mixing. - Also updated frazil ice documentation. --- docs/discrete_space.rst | 1 - docs/parameterizations_vertical.rst | 2 + docs/tracers.rst | 1 + docs/zotero.bib | 29 +++ src/core/_Finite_difference.dox | 5 - src/core/_Sea_ice.dox | 8 +- src/equation_of_state/_Equation_of_State.dox | 2 +- .../vertical/MOM_diabatic_aux.F90 | 2 +- src/parameterizations/vertical/_Frazil.dox | 33 +++ .../vertical/_Internal_tides.dox | 215 ++++++++++++++++++ 10 files changed, 289 insertions(+), 9 deletions(-) delete mode 100644 src/core/_Finite_difference.dox create mode 100644 src/parameterizations/vertical/_Frazil.dox create mode 100644 src/parameterizations/vertical/_Internal_tides.dox diff --git a/docs/discrete_space.rst b/docs/discrete_space.rst index b954915256..08a41a5f2d 100644 --- a/docs/discrete_space.rst +++ b/docs/discrete_space.rst @@ -13,7 +13,6 @@ algorithm. :maxdepth: 2 api/generated/pages/Discrete_Grids - api/generated/pages/Finite_Difference_Operators api/generated/pages/PPM api/generated/pages/Discrete_Coriolis api/generated/pages/Discrete_PG diff --git a/docs/parameterizations_vertical.rst b/docs/parameterizations_vertical.rst index 27285034d7..0d22787294 100644 --- a/docs/parameterizations_vertical.rst +++ b/docs/parameterizations_vertical.rst @@ -26,6 +26,8 @@ Kappa-shear Internal-tide driven mixing The schemes of :cite:`st_laurent2002`, :cite:`polzin2009`, and :cite:`melet2012`, are all implemented through MOM_set_diffusivity and MOM_diabatic_driver. + :ref:`Internal_Tidal_Mixing` + Vertical friction ----------------- diff --git a/docs/tracers.rst b/docs/tracers.rst index 6190fe096d..8b5a21ee12 100644 --- a/docs/tracers.rst +++ b/docs/tracers.rst @@ -9,4 +9,5 @@ Tracers in MOM6 api/generated/pages/Horizontal_Diffusion.rst api/generated/pages/Vertical_Diffusion.rst api/generated/pages/Passive_Tracers + api/generated/pages/Frazil_Ice diff --git a/docs/zotero.bib b/docs/zotero.bib index f0e1a3b44d..957097f217 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -847,6 +847,16 @@ @article{melet2012 pages = {602--615} } +@article{simmons2004, + title = {Tidally driven mixing in a numerical model of the ocean general circulation}, + volume = {6}, + author = {Simmons, H. L. and S. R. Jayne and L. C. St. Laurent and A. J. Weaver}, + year = {2004}, + journal = {Ocean Modell.}, + pages = {245--263}, + doi = {10.1016/S1463-5003(03)00011-8} +} + @article{polzin2009, title = {An abyssal recipe}, volume = {30}, @@ -863,6 +873,15 @@ @article{polzin2009 pages = {298--309} } +@article{polzin2004, + title = {Idealized solutions for the energy balance of the finescale internal wave field}, + volume = {34}, + journal = {J. Phys. Oceanogr.}, + author = {Polzin, Kurt L.}, + year = {2004}, + pages = {231--246} +} + @article{white2009, title = {High-order regridding-remapping schemes for continuous isopycnal and generalized coordinates in ocean models}, volume = {228}, @@ -2524,3 +2543,13 @@ @article{hallberg2005 volume = {8}, doi = {10.1016/j.ocemod.2004.01.001} } + +@article{bell1975, + author = {T. H. Bell}, + year = {1975}, + title = {Lee wavews in stratified flows with simple harmonic time dependence"}, + journal = {J. Fluid Mech.}, + volume = {67}, + pages = {705--722} +} + diff --git a/src/core/_Finite_difference.dox b/src/core/_Finite_difference.dox deleted file mode 100644 index ecbd37d8b7..0000000000 --- a/src/core/_Finite_difference.dox +++ /dev/null @@ -1,5 +0,0 @@ -/*! \page Finite_Difference_Operators Finite Difference Operators - -\brief Finite Difference Operators - -*/ diff --git a/src/core/_Sea_ice.dox b/src/core/_Sea_ice.dox index bec05af17c..232bac1bb8 100644 --- a/src/core/_Sea_ice.dox +++ b/src/core/_Sea_ice.dox @@ -1,5 +1,11 @@ /*! \page Sea_Ice Sea Ice Considerations -\section Frazil Ice Formation +\section section_seaice Sea Ice Considerations + +For realistic domains, it is assumed that MOM6 will be run in a coupled mode, such that either the +sea-ice model or the coupler will be computing atmospheric bulk fluxes and passing them to the ocean. +Likewise, MOM6 can compute the frazil ice formation as described in \ref section_frazil, which it +then passes to the sea-ice model, expecting to get back the rejected brine or melted fresh water in +return. */ diff --git a/src/equation_of_state/_Equation_of_State.dox b/src/equation_of_state/_Equation_of_State.dox index 2e18b49f54..791c7001b1 100644 --- a/src/equation_of_state/_Equation_of_State.dox +++ b/src/equation_of_state/_Equation_of_State.dox @@ -36,7 +36,7 @@ Compute the required quantities using the equation of state from \cite jackett19 Compute the required quantities using the equation of state from [TEOS-10](http://www.teos-10.org/). -\section TFREEZE Freezing Temperature of Sea Water +\section section_TFREEZE Freezing Temperature of Sea Water There are three choices for computing the freezing point of sea water: diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index ee27c6c5df..072bc1445e 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -45,7 +45,7 @@ module MOM_diabatic_aux real :: rivermix_depth = 0.0 !< The depth to which rivers are mixed if do_rivermix = T [Z ~> m]. logical :: reclaim_frazil !< If true, try to use any frazil heat deficit to !! to cool the topmost layer down to the freezing - !! point. The default is false. + !! point. The default is true. logical :: pressure_dependent_frazil !< If true, use a pressure dependent !! freezing temperature when making frazil. The !! default is false, which will be faster but is diff --git a/src/parameterizations/vertical/_Frazil.dox b/src/parameterizations/vertical/_Frazil.dox new file mode 100644 index 0000000000..06321231a3 --- /dev/null +++ b/src/parameterizations/vertical/_Frazil.dox @@ -0,0 +1,33 @@ +/*! \page Frazil_Ice Frazil Ice Formation + +\section section_frazil Frazil Ice Formation + +Frazil ice forms in the model when the in situ temperature drops below +the local freezing point, taking into account the in situ salinity and +pressure. Starting at the bottom and working up through the water column, +if the water is below freezing, set it to freezing and add the heat +required to the heat deficit. If the water above is warmer than freezing, +use that heat to take away the heat deficit and to cool the water. If +you get all the way to the surface with a heat deficit, that quantity +is passed to the ice model as a heat flux it will need to provide to +the ocean. + +The local freezing point code is provided by the equation of state being +used by MOM6. See \ref section_TFREEZE for the MOM6 options. + +The salinity is adjusted only at the surface when frazil ice is +formed. This happens when the ice model creates ice with the heat deficit, +taking salt out of the surface waters. We inherit this behavior from +older versions of MOM, but the effect of not adjusting the in situ +salinity is thought to be small. + +Note that versions simply whisking all the heat deficit to the surface +without checking for warm water above tended to produce rapidly-melting +ice floes in warm waters. This was deemed unphysical and was corrected. + +A similar process that we are also omitting is the formation of salt +crystals when the salinity becomes too high. The salt crystals should +form and sink, leaving a layer on the bed that will be diluted when the +salinity drops again. This process can be seen in a lake in Death Valley. + +*/ diff --git a/src/parameterizations/vertical/_Internal_tides.dox b/src/parameterizations/vertical/_Internal_tides.dox new file mode 100644 index 0000000000..16be5a695c --- /dev/null +++ b/src/parameterizations/vertical/_Internal_tides.dox @@ -0,0 +1,215 @@ +/*! \page Internal_Tidal_Mixing Internal Tidal Mixing + +Two parameterizations of vertical mixing due to internal tides are +available with the option INT_TIDE_DISSIPATION. The first is that of +\cite st_laurent2002 while the second is that of \cite polzin2009. Choose +between them with the INT_TIDE_PROFILE option. There are other relevant +paramters which can be seen in MOM_parameter_doc.all once the main tidal +dissipation switch is turned on. + +\section section_st_laurent St Laurent et al. + +The estimated turbulent dissipation rate of +internal tide energy \f$\epsilon\f$ is: + +\f[ + \epsilon = \frac{q E(x,y)}{\rho} F(z). +\f] + +where \f$\rho\f$ is the reference density of seawater, \f$E(x,y)\f$ is +the energy flux per unit area transferred from barotropic to baroclinic +tides, \f$q\f$ is the fraction of the internal-tide energy dissipated +locally, and \f$F(z)\f$ is the vertical structure of the dissipation. +This \f$q\f$ is estimated to be roughly 0.3 based on observations. The +term \f$E(x,y)\f$ is given by \cite st_laurent2002 as: + +\f[ + E(x,y) \simeq \frac{1}{2} \rho N_b \kappa h^2 \langle U^2 \rangle +\f] + +where \f$\rho\f$ is the reference density of seawater, \f$N_b\f$ is +the buoyancy frequency along the seafloor, and \f$(\kappa, h)\f$ are +the wavenumber and amplitude scales for the topographic roughness, and +\f$\langle U^2 \rangle\f$ is the barotropic tide variance. It is assumed +that the model will read in topographic roughness squared \f$h^2\f$ +from a file (the variable must be named "h2"). + +To convert from energy dissipation to vertical diffusion \f$K_d\f$, +the simple estimate is: + +\f[ + K_d \approx \frac{\Gamma q E(x,y) F(z)}{\rho N^2} +\f] + +where \f$\Gamma\f$ is the mixing efficiency, generally set to 0.2 +and \f$F(z)\f$ is a vertical structure function with exponential decay +from the bottom: + +\f[ + F(z) = \frac{e^{-(H+z)/\zeta}}{\zeta (1 - e^{H/\zeta}}. +\f] + +Here, \f$\zeta\f$ is a vertical decay scale with a default of 500 meters. +One change in MOM6 from the St. Laurent scheme is to use this form of \f$\Gamma\f$: + +\f[ + \Gamma = 0.2 \frac{N^2}{N^2 + \Omega^2} +\f] + +instead of \f$\Gamma = 0.2\f$, where \f$\Omega\f$ is the angular velocity +of the Earth. This allows the buoyancy fluxes to tend to zero in regions +of very weak stratification, allowing a no-flux bottom boundary condition +to be satisfied. + +This \f$K_d\f$ gets added to the diffusivity due to the background and +other contributions unless you set BBL_MIXING_AS_MAX to True, in which +case the maximum of all the contributions is used. + +\section section_polzin Polzin + +The vertical diffusion profile of \cite polzin2009 is a WKB-stretched +algebraic decay profile. It is based on a radiation balance equation, +which links the dissipation profile associtated with internal breaking to +the finescale internal wave shear producing that dissipation. The vertical +profile of internal-tide driven energy dissipation can then vary in time +and space, and evolve in a changing climate (\cite melet2012). \cite +melet2012 describes how the Polzin scheme is implemented in MOM6, +copied here. + +The parameterization of \cite polzin2009 links the energy dissipation +profile to the finescale internal wave shear producing that +dissipation, using an idealized vertical wavenumber energy spectrum +to identify analytic solutions to a radiation balance equation +(\cite polzin2004). These solutions yield a dissipation profile +\f$\epsilon(z)\f$: + +\f[ + \epsilon = \frac{\epsilon_0}{[1 + (z/z_p)]^2}, +\f] + +where the magnitude \f$\epsilon_0\f$ and scale height \f$z_p\f$ can be expressed in terms of the +spectral amplitude and bandwidth of the idealized vertical wavenumber energy spectrum in uniform +stratification (\cite polzin2009). + +To take into account the nonuniform stratification, \cite polzin2009 applied a buoyancy scaling +using the Wentzel-Kramers-Brillouin (WKB) approximation. As a result, the vertical wavenumber of a +wave packet varies in proportion to the buoyancy frequency \f$N\f$, which in turn implies an +additional transport of energy to smaller scales, and thus a possible enhanced mixing in regions of +strong stratification. Such effects can be described by buoyancy scaling the vertical coordinate +\f$z\f$ as + +\f[ + z^{\ast}(z) = \int_{0}^{z} \left[ \frac{N^2 (z^\prime )}{N_b^2} \right] dz^{\prime} , +\f] + +with \f$z^\prime\f$ being positive upward relative to the bottom of the ocean. The turbulent +dissipation rate then becomes + +\f[ + \epsilon = \frac{\epsilon_0}{[1 + (z^{\ast} /z_p)]^2} \frac{N^2(z)}{N_b^2} . +\f] + +The spectral amplitude and bandwidth of the idealized vertical wavenumber +energy spectrum are identified after WKB scaling using a quasi-linear +spectral model of internal-tide generation that incorporates horizontal +advection of the barotropic tide into the momentum equation (\cite bell1975). +As a result, Polzin's formulation leads to an expression for +the spatially and temporally varying dissipation of internal tide energy +at the bottom \f$\epsilon_0\f$, and the vertical scale of decay for the +dissipation of internal tide energy \f$z_p\f$. + +\subsection subsection_energy_conserving Energy-conserving form + +To satisfy energy conservation (the integral of the vertical structure for the turbulent dissipation +over depth should be unity), the dissipation is rewritten as + +\f[ + \epsilon = \frac{\epsilon_0 z_p}{1 + (z^\ast/z_p)]^2} \frac{N^2(z)}{N^2_b} \left[ + \frac{1}{z^{\ast(z=H)}} + \frac{1}{z_p} \right] . +\f] + +In the MOM6 implementation, we use the \cite st_laurent2002 template for the vertical flux of energy +at the ocean floor, so that in both formulations: + +\f[ + \int_{0}^{H} \epsilon (z) dz = \frac{qE}{\rho} . +\f] + +Whereas \cite polzin2009 assumed tthat the total dissipation was locally in balance with the +barotropic to baroclinic energy conversion rate \f$(q=1)\f$, here we use the \cite simmons2004 value +of \f$q=1/3\f$ to retain as much consistency as passible between both parameterizations. + +\subsection subsection_vertical_decay_scale Vertical decay-scale reformulation + +We follow the \cite polzin2009 prescription for the vertical scale of decay for the dissipation of +internal-tide energy. However, we assume that the topographic power law, denoted as \f$\nu\f$ in \cite +polzin2009, is equal to 1 (instead of 0.9) and we reformulated the expression of \f$z_p\f$ to put it +in a more readable form: + +\f[ + z_p = \frac{z_p^\mbox{ref} (\kappa^\mbox{ref})^2 (h^\mbox{ref})^2 (N_b^\mbox{ref})^3} {U^\mbox{ref}} + \frac{U}{h^2 \kappa^2 N_b^3} . +\f] + +The superscript ref refers to reference values of the various parameters, as given by +observations from the Brazil basin. Therefore, the above can be rewritten as + +\f[ + z_p = \mu (N_b^\mbox{ref} )^2 + \frac{U}{h^2 \kappa^2 N_b^3} . +\f] + +where \f$\mu\f$ is a nondimensional constant \f$(\mu = 0.06970)\f$ and \f$N_b^\mbox{ref} = 9.6 \times +10^{-4} s^{-1}\f$. Finally, a minimum decay scale of \f$z_p = 100 m\f$ is imposed in our +implementation. + +\subsection subsection_reformulation_WKB Reformulation of the WKB scaling + +Since the dissipation is expressed as a function of the ratio \f$z^\ast / z_p\f$, a different WKB +scaling can be used so long as we modify \f$z_p\f$ accordingly. In the implemented parameterization, +we define the scaled height coordinate \f$z^\ast\f$ by + +\f[ + z^\ast (z) = \frac{1}{\overline{N^2 (z)}^z} \int_{0}^{z} N^2(z^\prime ) dz ^\prime , +\f] + +with \f$z^\prime\f$ defined to be the height above the ocean bottom. By normalizing \f$N^2\f$ by its +vertical mean \f$\overline{N^2}^z\f$, \f$z^\ast\f$ ranges from \f$0\f$ to \f$H\f$, the depth of the +ocean. + +The WKB-scaled vertical decay scale for the Polzin formulation becomes + +\f[ + z^\ast_p = \mu(N_b^\mbox{ref})^2 \frac{U}{h^2 \kappa^2 N_b \overline{N^2}^z} . +\f] + +Unlike the \cite st_laurent2002 parameterization, the vertical decay scale now depends on physical +variables and can evolve with a changing climate. + +Finally, the Polzin vertical profile of dissipation implemented in the model is given by + +\f[ + \epsilon = \frac{qE(x,y)}{\rho [1 + (z^\ast/z_p^\ast)]^2} \frac{N^2(z)}{\overline{N^2}^z} + \left( \frac{1}{H} + \frac{1}{z_p^\ast} \right) . +\f] + +In both parameterizations, turbulent diapycnal diffusivities are inferred from the dissipation +\f$\epsilon\f$ by: + +\f[ + K_d = \frac{\Gamma \epsilon}{N^2} +\f] + +and using this form of \f$\Gamma\f$: + +\f[ + \Gamma = 0.2 \frac{N^2}{N^2 + \Omega^2} +\f] + +instead of \f$\Gamma = 0.2\f$, where \f$\Omega\f$ is the angular velocity +of the Earth. This allows the buoyancy fluxes to tend to zero in regions +of very weak stratification, allowing a no-flux bottom boundary condition +to be satisfied. + +*/ + From c9c4cc9b37b842017c10a82bc8b1154b2cbb5895 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 19 Aug 2021 16:29:09 -0800 Subject: [PATCH 08/15] Fixed the \cite arguments? --- src/parameterizations/vertical/_Internal_tides.dox | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/parameterizations/vertical/_Internal_tides.dox b/src/parameterizations/vertical/_Internal_tides.dox index 16be5a695c..882b73dd1b 100644 --- a/src/parameterizations/vertical/_Internal_tides.dox +++ b/src/parameterizations/vertical/_Internal_tides.dox @@ -72,8 +72,8 @@ algebraic decay profile. It is based on a radiation balance equation, which links the dissipation profile associtated with internal breaking to the finescale internal wave shear producing that dissipation. The vertical profile of internal-tide driven energy dissipation can then vary in time -and space, and evolve in a changing climate (\cite melet2012). \cite -melet2012 describes how the Polzin scheme is implemented in MOM6, +and space, and evolve in a changing climate (\cite melet2012). \cite melet2012 +describes how the Polzin scheme is implemented in MOM6, copied here. The parameterization of \cite polzin2009 links the energy dissipation @@ -141,10 +141,11 @@ of \f$q=1/3\f$ to retain as much consistency as passible between both parameteri \subsection subsection_vertical_decay_scale Vertical decay-scale reformulation -We follow the \cite polzin2009 prescription for the vertical scale of decay for the dissipation of -internal-tide energy. However, we assume that the topographic power law, denoted as \f$\nu\f$ in \cite -polzin2009, is equal to 1 (instead of 0.9) and we reformulated the expression of \f$z_p\f$ to put it -in a more readable form: +We follow the \cite polzin2009 prescription for the vertical scale of +decay for the dissipation of internal-tide energy. However, we assume +that the topographic power law, denoted as \f$\nu\f$ in \cite polzin2009, +is equal to 1 (instead of 0.9) and we reformulated the expression of +\f$z_p\f$ to put it in a more readable form: \f[ z_p = \frac{z_p^\mbox{ref} (\kappa^\mbox{ref})^2 (h^\mbox{ref})^2 (N_b^\mbox{ref})^3} {U^\mbox{ref}} From 72a525ad8cf1369a1fed06a3b27b627b96deff2e Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 27 Aug 2021 17:49:49 -0800 Subject: [PATCH 09/15] Go back to gfdl verison of MOM_boundary_update --- src/core/MOM_boundary_update.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/core/MOM_boundary_update.F90 b/src/core/MOM_boundary_update.F90 index 599e4b92ca..dc89f3f92c 100644 --- a/src/core/MOM_boundary_update.F90 +++ b/src/core/MOM_boundary_update.F90 @@ -143,6 +143,9 @@ subroutine update_OBC_data(OBC, G, GV, US, tv, h, CS, Time) type(update_OBC_CS), pointer :: CS !< Control structure for OBCs type(time_type), intent(in) :: Time !< Model time +! Something here... with CS%file_OBC_CSp? +! if (CS%use_files) & +! call update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (CS%use_tidal_bay) & call tidal_bay_set_OBC_data(OBC, CS%tidal_bay_OBC_CSp, G, GV, h, Time) if (CS%use_Kelvin) & From 710b087c7a4a11c4ed5a4780540b2e9d193c756e Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 27 Aug 2021 17:52:15 -0800 Subject: [PATCH 10/15] More documentation - Jackson start. --- docs/parameterizations_vertical.rst | 7 +- docs/zotero.bib | 107 ++++++++++++- .../vertical/MOM_set_diffusivity.F90 | 10 -- .../vertical/MOM_set_viscosity.F90 | 74 --------- .../vertical/_Internal_tides.dox | 11 ++ .../vertical/_V_diffusivity.dox | 150 ++++++++++++++++++ .../vertical/_V_viscosity.dox | 64 ++++++++ 7 files changed, 333 insertions(+), 90 deletions(-) create mode 100644 src/parameterizations/vertical/_V_diffusivity.dox create mode 100644 src/parameterizations/vertical/_V_viscosity.dox diff --git a/docs/parameterizations_vertical.rst b/docs/parameterizations_vertical.rst index 0d22787294..c9404c5088 100644 --- a/docs/parameterizations_vertical.rst +++ b/docs/parameterizations_vertical.rst @@ -21,9 +21,12 @@ Interior and bottom-driven mixing --------------------------------- Kappa-shear - MOM_kappa_shear implement the shear-driven mixing of :cite:`jackson2008`. + MOM_kappa_shear implements the shear-driven mixing of :cite:`jackson2008`. + + :ref:`Internal_Shear_Mixing` Internal-tide driven mixing + The schemes of :cite:`st_laurent2002`, :cite:`polzin2009`, and :cite:`melet2012`, are all implemented through MOM_set_diffusivity and MOM_diabatic_driver. :ref:`Internal_Tidal_Mixing` @@ -33,6 +36,8 @@ Vertical friction Vertical viscosity is implemented in MOM_vert_frict and coefficient computed in MOM_set_viscosity, although contributions to viscosity from other parameterizations are calculated in those respective modules (e.g. MOM_kappa_shear, MOM_KPP, MOM_energetic_PBL). + :ref:`Vertical_Viscosity` + Vertical diffusion ------------------ diff --git a/docs/zotero.bib b/docs/zotero.bib index 957097f217..8a5c14dcfa 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -655,6 +655,30 @@ @article{killworth1992 pages = {1379--1387} } +@article{killworth1999, + doi = {10.1175/1520-0485(1999)029<1221:atbblc>2.0.co;2}, + year = 1999, + publisher = {American Meteorological Society}, + volume = {29}, + number = {6}, + pages = {1221--1238}, + author = {P. D. Killworth and N. R. Edwards}, + title = {A Turbulent Bottom Boundary Layer Code for Use in Numerical Ocean Models}, + journal = {J. Phys. Oceanography} +} + +@article{zilitinkevich1996, + doi = {10.1007/bf02430334}, + year = 1996, + publisher = {Springer Science and Business Media {LLC}}, + volume = {81}, + number = {3-4}, + pages = {325--351}, + author = {S. Zilitinkevich and D. V. Mironov}, + title = {A multi-limit formulation for the equilibrium depth of a stably stratified boundary layer}, + journal = {Boundary-Layer Meteorology} +} + @article{gent1995, title = {Parameterizing {Eddy}-{Induced} {Tracer} {Transports} in {Ocean} {Circulation} {Models}}, volume = {25}, @@ -800,6 +824,18 @@ @article{jackson2008 pages = {1033--1053} } +@article{turner1986, + doi = {10.1017/s0022112086001222}, + year = 1986, + publisher = {Cambridge University Press ({CUP})}, + volume = {173}, + pages = {431--471}, + author = {J. S. Turner}, + title = {Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows}, + journal = {J. Fluid Mech.} +} + + @article{reichl2018, title = {A simplified energetics based planetary boundary layer ({ePBL}) approach for ocean climate simulations.}, volume = {132}, @@ -1761,6 +1797,18 @@ @article{large1994 pages = {363--403} } +@article{pacanowski1981, + doi = {10.1175/1520-0485(1981)011<1443:povmin>2.0.co;2}, + year = 1981, + publisher = {American Meteorological Society}, + volume = {11}, + number = {11}, + pages = {1443--1451}, + author = {R. C. Pacanowski and S. G. H. Philander}, + title = {Parameterization of Vertical Mixing in Numerical Models of Tropical Oceans}, + journal = {J. Phys. Oceanography} +} + @article{van_roekel2018, title = {The {KPP} {Boundary} {Layer} {Scheme} for the {Ocean}: {Revisiting} {Its} {Formulation} and {Benchmarking} {One}-{Dimensional} {Simulations} {Relative} to {LES}}, volume = {10}, @@ -2343,6 +2391,19 @@ @article{hallberg2000 pages = {1402--1419} } +@article{umlauf2005, + doi = {10.1016/j.csr.2004.08.004}, + year = 2005, + publisher = {Elsevier {BV}}, + volume = {25}, + number = {7-8}, + pages = {795--827}, + author = {L. Umlauf and H. Burchard}, + title = {Second-order turbulence closure models for geophysical boundary layers. A review of recent work}, + journal = {Continental Shelf Res.} +} + + @article{easter1993, title = {Two Modified Versions of Bott's Positive-Definite Numerical Advection Scheme}, @@ -2545,11 +2606,47 @@ @article{hallberg2005 } @article{bell1975, - author = {T. H. Bell}, - year = {1975}, - title = {Lee wavews in stratified flows with simple harmonic time dependence"}, - journal = {J. Fluid Mech.}, + doi = {10.1017/s0022112075000560}, + year = 1975, + publisher = {Cambridge University Press ({CUP})}, volume = {67}, - pages = {705--722} + number = {4}, + pages = {705--722}, + author = {T. H. Bell}, + title = {Lee waves in stratified flows with simple harmonic time dependence}, + journal = {J. Fluid Mech.} +} + +@article{nikurashin2010a, + doi = {10.1175/2009jpo4199.1}, + year = 2010, + publisher = {American Meteorological Society}, + volume = {40}, + number = {5}, + pages = {1055--1074}, + author = {M. Nikurashin and R. Ferrari}, + title = {Radiation and Dissipation of Internal Waves Generated by Geostrophic Motions Impinging on Small-Scale Topography: Theory}, + journal = {J. Phys. Oceanography} +} + +@article{nikurashin2010b, + doi = {10.1175/2010jpo4315.1}, + year = 2010, + publisher = {American Meteorological Society}, + volume = {40}, + number = {9}, + pages = {2025--2042}, + author = {M. Nikurashin and R. Ferrari}, + title = {Radiation and Dissipation of Internal Waves Generated by Geostrophic Motions Impinging on Small-Scale Topography: Application to the Southern Ocean}, + journal = {J. Phys. Oceanography} } +@article{miles1961, + title = {On the stability of heterogeneous shear flows}, + author = {JW Miles}, + year = {1961}, + journal = {J. of Fluid Mech.}, + volume = {10}, + pages = {496--508}, + doi = {10.1017/S0022112061000305} +} diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index f4874252f4..0d07f0fea4 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -198,16 +198,6 @@ module MOM_set_diffusivity contains -!> Sets the interior vertical diffusion of scalars due to the following processes: -!! 1. Shear-driven mixing: two options, Jackson et at. and KPP interior; -!! 2. Background mixing via CVMix (Bryan-Lewis profile) or the scheme described by -!! Harrison & Hallberg, JPO 2008; -!! 3. Double-diffusion, old method and new method via CVMix; -!! 4. Tidal mixing: many options available, see MOM_tidal_mixing.F90; -!! In addition, this subroutine has the option to set the interior vertical -!! viscosity associated with processes 1,2 and 4 listed above, which is stored in -!! visc%Kv_slow. Vertical viscosity due to shear-driven mixing is passed via -!! visc%Kv_shear subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & G, GV, US, CS, Kd_lay, Kd_int, Kd_extra_T, Kd_extra_S) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 99bd91d8f8..9a2680ecc1 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -115,80 +115,6 @@ module MOM_set_visc contains !> Calculates the thickness of the bottom boundary layer and the viscosity within that layer. -!! -!! A drag law is used, either linearized about an assumed bottom velocity or using the -!! actual near-bottom velocities combined with an assumed unresolved velocity. The bottom -!! boundary layer thickness is limited by a combination of stratification and rotation, as -!! in the paper of Killworth and Edwards, JPO 1999. It is not necessary to calculate the -!! thickness and viscosity every time step; instead previous values may be used. -!! -!! \section set_viscous_BBL Viscous Bottom Boundary Layer -!! -!! If set_visc_cs.bottomdraglaw is True then a bottom boundary layer viscosity and thickness -!! are calculated so that the bottom stress is -!! \f[ -!! \mathbf{\tau}_b = C_d | U_{bbl} | \mathbf{u}_{bbl} -!! \f] -!! If set_visc_cs.bottomdraglaw is True then the term \f$|U_{bbl}|\f$ is set equal to the -!! value in set_visc_cs.drag_bg_vel so that \f$C_d |U_{bbl}|\f$ becomes a Rayleigh bottom drag. -!! Otherwise \f$|U_{bbl}|\f$ is found by averaging the flow over the bottom set_visc_cs.hbbl -!! of the model, adding the amplitude of tides set_visc_cs.tideamp and a constant -!! set_visc_cs.drag_bg_vel. For these calculations the vertical grid at the velocity -!! component locations is found by -!! \f[ -!! \begin{array}{ll} -!! \frac{2 h^- h^+}{h^- + h^+} & u \left( h^+ - h^-\right) >= 0 -!! \\ -!! \frac{1}{2} \left( h^- + h^+ \right) & u \left( h^+ - h^-\right) < 0 -!! \end{array} -!! \f] -!! which biases towards the thin cell if the thin cell is upwind. Biasing the grid toward -!! thin upwind cells helps increase the effect of viscosity and inhibits flow out of these -!! thin cells. -!! -!! After diagnosing \f$|U_{bbl}|\f$ over a fixed depth an active viscous boundary layer -!! thickness is found using the ideas of Killworth and Edwards, 1999 (hereafter KW99). -!! KW99 solve the equation -!! \f[ -!! \left( \frac{h_{bbl}}{h_f} \right)^2 + \frac{h_{bbl}}{h_N} = 1 -!! \f] -!! for the boundary layer depth \f$h_{bbl}\f$. Here -!! \f[ -!! h_f = \frac{C_n u_*}{f} -!! \f] -!! is the rotation controlled boundary layer depth in the absence of stratification. -!! \f$u_*\f$ is the surface friction speed given by -!! \f[ -!! u_*^2 = C_d |U_{bbl}|^2 -!! \f] -!! and is a function of near bottom model flow. -!! \f[ -!! h_N = \frac{C_i u_*}{N} = \frac{ (C_i u_* )^2 }{g^\prime} -!! \f] -!! is the stratification controlled boundary layer depth. The non-dimensional parameters -!! \f$C_n=0.5\f$ and \f$C_i=20\f$ are suggested by Zilitinkevich and Mironov, 1996. -!! -!! If a Richardson number dependent mixing scheme is being used, as indicated by -!! set_visc_cs.rino_mix, then the boundary layer thickness is bounded to be no larger -!! than a half of set_visc_cs.hbbl . -!! -!! \todo Channel drag needs to be explained -!! -!! A BBL viscosity is calculated so that the no-slip boundary condition in the vertical -!! viscosity solver implies the stress \f$\mathbf{\tau}_b\f$. -!! -!! \subsection set_viscous_BBL_ref References -!! -!! \arg Killworth, P. D., and N. R. Edwards, 1999: -!! A Turbulent Bottom Boundary Layer Code for Use in Numerical Ocean Models. -!! J. Phys. Oceanogr., 29, 1221-1238, -!! doi:10.1175/1520-0485(1999)029<1221:ATBBLC>2.0.CO;2 -!! \arg Zilitinkevich, S., Mironov, D.V., 1996: -!! A multi-limit formulation for the equilibrium depth of a stably stratified boundary layer. -!! Boundary-Layer Meteorology 81, 325-351. -!! doi:10.1007/BF02430334 -!! subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. diff --git a/src/parameterizations/vertical/_Internal_tides.dox b/src/parameterizations/vertical/_Internal_tides.dox index 882b73dd1b..bf6e16ea5c 100644 --- a/src/parameterizations/vertical/_Internal_tides.dox +++ b/src/parameterizations/vertical/_Internal_tides.dox @@ -212,5 +212,16 @@ of the Earth. This allows the buoyancy fluxes to tend to zero in regions of very weak stratification, allowing a no-flux bottom boundary condition to be satisfied. +\section Nikurashin Lee Wave Mixing + +If one has the INT_TIDE_DISSIPATION flag on, there is an option to also turn on +LEE_WAVE_DISSIPATION. The theory is presented in \cite nikurashin2010a +while the application of it is presented in \cite nikurashin2010b. For +the implementation in MOM6, it is required that you provide an estimate +of the TKE loss due to the Lee waves which is then applied with either +the St. Laurent or the Polzin vertical profile. + +IS THERE A SCRIPT to produce this somewhere or what??? + */ diff --git a/src/parameterizations/vertical/_V_diffusivity.dox b/src/parameterizations/vertical/_V_diffusivity.dox new file mode 100644 index 0000000000..5c40768eaf --- /dev/null +++ b/src/parameterizations/vertical/_V_diffusivity.dox @@ -0,0 +1,150 @@ +/*! \page Internal_Shear_Mixing Internal Vertical Mixing + +Sets the interior vertical diffusion of scalars due to the following processes: + +-# Shear-driven mixing: two options, \cite jackson2008 and KPP interior; +-# Background mixing via CVMix (Bryan-Lewis profile) or the scheme described by + \cite harrison2008. +-# Double-diffusion, old method and new method via CVMix; +-# Tidal mixing: many options available, see \ref Internal_Tidal_Mixing. + +In addition, the MOM_set_diffusivity has the option to set the interior vertical +viscosity associated with processes 1,2 and 4 listed above, which is stored in +visc\%Kv\_slow. Vertical viscosity due to shear-driven mixing is passed via +visc\%Kv\_shear + +The resulting diffusivity, \f$K_d\f$, is the sum of all the contributions +unless you set BBL_MIXING_AS_MAX to True, in which case the maximum of +all the contributions is used. + +In addition, \f$K_d\f$ is multiplied by the term: + +\f[ + \frac{N^2}{N^2 + \Omega^2} +\f] + +where \f$N\f$ is the buoyancy frequency and \f$\Omega\f$ is the angular velocity +of the Earth. This allows the buoyancy fluxes to tend to zero in regions +of very weak stratification, allowing a no-flux bottom boundary condition +to be satisfied. + +\section section_Shear Shear-driven Mixing + +Below the surface mixed layer, there are places in the world's oceans +where shear mixing is known to take place. This shear-driven mixing can +be represented in MOM6 through either CVMix or the parameterization of +\cite jackson2008. + +\subsection subsection_CVMix_shear Shear-driven mixing in CVMix + +The community vertical mixing (CVMix) code contains options for shear +mixing from either \cite large1994 or from \cite pacanowski1981. In MOM6, +CVMix is included via a git submodule which loads the external CVMix +package. The shear mixing routine in CVMix was developed to reproduce the +observed mixing of the equatorial undercurrent in the Pacific. + +We first compute the gradient Richardson number \f$\mbox{Ri} = N^2 / S^2\f$, +where \f$S\f$ is the vertical shear (\f$S = ||\bf{u}_z ||\f$) and \f$N\f$ +is the buoyancy frequency (\f$N^2 = -g \rho_z / \rho_0\f$). The +parameterization of \cite large1994 is as follows, where the diffusivity \f$\kappa\f$ +is given by + +\f[ + \kappa = \kappa_0 \left[ 1 - \min \left( 1, \frac{\mbox{Ri}}{\mbox{Ri}_c} \right) ^2 \right] ^3 , +\f] + +with \f$\kappa_0 = 5 \times 10^{-3}\, \mbox{m}^2 \,\mbox{s}^{-1}\f$ and \f$\mbox{Ri}_c = 0.7\f$. + +\subsection subsection_kappa_shear Shear-driven mixing in Jackson + +While the above parameterization works well enough in the equatorial +Pacific, another place one can expect shear-mixing to matter is +in overflows of dense water. \cite jackson2008 proposes a new shear +parameterization with the goal of working in both the equatorial undercurrent +and for overflows, also to have smooth transitions between unstable and +stable regions. Their scheme looks like: + +\f{eqnarray} + \frac{\partial^2 \kappa}{\partial z^2} - \frac{\kappa}{L^2_d} &= - 2 SF(\mbox{Ri}) . + \label{eq:Jackson_10} +\f} + +This is similar to the locally constant stratification limit of +\cite turner1986, but with the addition of a decay length scale +\f$L_d = \lambda L_b\f$. Here \f$L_b = Q^{1/2} / N\f$ is the buoyancy +length scale where \f$Q\f$ is the turbulent kinetic energy (TKE) per +unit mass, and \f$\lambda\f$ is a nondimensional constant. The function +\f$F(\mbox{Ri})\f$ is a function of the Richardson number that remains +to be determined. As in \cite turner1986, there must be a critical +value of \f$\mbox{Ri}\f$ above which \f$F(\mbox{Ri}) = 0\f$. There +are two length scales: the width of the low Richardson number region +as in \cite turner1986, and the buoyancy length scale, which is the +length scale over which the TKE is affected by the stratification (see +\cite jackson2008 for more details). In particular, the inclusion of a +decay length scale means that the diffusivity decays exponentially away +from the mixing region with a length scale of \f$L_d\f$. This is important +since turbulent eddies generated in the low \f$\mbox{Ri}\f$ layer can +be vertically self-advected and mix nearby regions. This method yields +a smoother diffusivity than that in \cite hallberg2000, especially in +areas where the Richardson number is noisy. + +This parameterization predicts the turbulent eddy diffusivity in terms +of the vertical profiles of velocity and density, providing that the +TKE is known. To complete the parameterization we use a TKE \f$Q\f$ +budget such as that used in second-order turbulence closure models +(\cite umlauf2005). We make a few additional assumptions, however, +and use the simplified form + +\f{eqnarray} + \frac{\partial}{\partial z} \left[ (\kappa + \nu_0) \frac{\partial Q} + {\partial z} \right] + \kappa (S^2 - N^2) - Q(c_N N + c_S S) &= 0. + \label{eq:Jackson_11} +\f} + +The system is therefore in balance between a vertical diffusion of +TKE caused by both the eddy and molecular viscosity \f$(\nu_0)\f$, +the production of TKE by shear, a sink due to stratification, and the +dissipation. Note that we are assuming a Prandtl number of 1, although a +parameterization for the Prandtl number could be added. We have assumed +that the TKE reaches a quasi-steady state faster than the flow is evolving +and faster than it can be affected by mean-flow advection so that \f$DQ/Dt = +0\f$. Since this parameterization is meant to be used in climate models +with low horizontal resolution and large time steps compared to the +mixing time scales, this is a reasonable assumtion. The most tenuous +assumption is in the form of the dissipation \f$\epsilon = Q(C_N N + +c_S S)\f$ (where \f$c_N\f$ and \f$c_S\f$ are to be determined), +which is assumed to be dependent on the buoyancy frequeny (through loss +of energy to internal waves) and the velocity shear (through the energy +cascade to smaller scales). + +We can rewrite \eqref{eq:Jackson_10} as the steady "transport" equation +for the turbulent diffusivity (i.e., with \f$D\kappa/Dt = 0\f$), + +\f[ + \frac{\partial}{\partial z} \left( \kappa \frac{\partial \kappa}{\partial z} + \right) + 2\kappa SF(\mbox{Ri}) - \left( \frac{\kappa}{L_d} \right)^2 - + \left( \frac{\partial \kappa}{\partial z} \right) ^2 = 0 . +\f] + +The first term on the left can be regarded as a vertical transport of +diffusivity, the second term as a source, and the final two as sinks. +This equation with \eqref{eq:Jackson_11} are simple enough to solve quickly +using an iterative technique. + +We also need boundary contitions for \eqref{eq:Jackson_10} +and \eqref{eq:Jackson_11}. For the turbulent diffusivity we use +\f$\kappa = 0\f$ since our diffusivity is numerically defined on +layer interfaces. This ensures that there is no turbulent flux across +boundaries. For the TKE we use boundary conditions of \f$Q = Q_0\f$ where +\f$Q_0\f$ is a constant value of TKE, used to prevent a singularity +in \eqref{eq:Jackson_10}, that is chosen to be small enough to not +influence results. Note that the value of \f$\kappa\f$ calculated here +reflects shear-driven turbulent mixing only; the total diffusivity would +be this value plus any diffusivities due to other turbulent processes +or a background value. + +\section section_Background Background Mixing + +\section section_Double_Diff Double Diffusion + +*/ diff --git a/src/parameterizations/vertical/_V_viscosity.dox b/src/parameterizations/vertical/_V_viscosity.dox new file mode 100644 index 0000000000..cc59e83457 --- /dev/null +++ b/src/parameterizations/vertical/_V_viscosity.dox @@ -0,0 +1,64 @@ +/*! \page Vertical_Viscosity Viscous Bottom Boundary Layer + +A drag law is used, either linearized about an assumed bottom velocity or using the +actual near-bottom velocities combined with an assumed unresolved velocity. The bottom +boundary layer thickness is limited by a combination of stratification and rotation, as +in the paper of \cite killworth1999. It is not necessary to calculate the +thickness and viscosity every time step; instead previous values may be used. + +\section set_viscous_BBL Viscous Bottom Boundary Layer + +If set_visc_CS\%bottomdraglaw is True then a bottom boundary layer viscosity and thickness +are calculated so that the bottom stress is +\f[ +\mathbf{\tau}_b = C_d | U_{bbl} | \mathbf{u}_{bbl} +\f] +If set_visc_CS\%bottomdraglaw is True then the term \f$|U_{bbl}|\f$ is set equal to the +value in set_visc_CS.drag_bg_vel so that \f$C_d |U_{bbl}|\f$ becomes a Rayleigh bottom drag. +Otherwise \f$|U_{bbl}|\f$ is found by averaging the flow over the bottom set_visc_CS\%hbbl +of the model, adding the amplitude of tides set_visc_CS\%tideamp and a constant +set_visc_CS\%drag_bg_vel. For these calculations the vertical grid at the velocity +component locations is found by +\f[ +\begin{array}{ll} +\frac{2 h^- h^+}{h^- + h^+} & u \left( h^+ - h^-\right) >= 0 +\\ +\frac{1}{2} \left( h^- + h^+ \right) & u \left( h^+ - h^-\right) < 0 +\end{array} +\f] +which biases towards the thin cell if the thin cell is upwind. Biasing the grid toward +thin upwind cells helps increase the effect of viscosity and inhibits flow out of these +thin cells. + +After diagnosing \f$|U_{bbl}|\f$ over a fixed depth an active viscous boundary layer +thickness is found using the ideas of Killworth and Edwards, 1999 (hereafter KW99). +KW99 solve the equation +\f[ +\left( \frac{h_{bbl}}{h_f} \right)^2 + \frac{h_{bbl}}{h_N} = 1 +\f] +for the boundary layer depth \f$h_{bbl}\f$. Here +\f[ +h_f = \frac{C_n u_*}{f} +\f] +is the rotation controlled boundary layer depth in the absence of stratification. +\f$u_*\f$ is the surface friction speed given by +\f[ +u_*^2 = C_d |U_{bbl}|^2 +\f] +and is a function of near bottom model flow. +\f[ +h_N = \frac{C_i u_*}{N} = \frac{ (C_i u_* )^2 }{g^\prime} +\f] +is the stratification controlled boundary layer depth. The non-dimensional parameters +\f$C_n=0.5\f$ and \f$C_i=20\f$ are suggested by \cite zilitinkevich1996. + +If a Richardson number dependent mixing scheme is being used, as indicated by +set_visc_CS\%rino_mix, then the boundary layer thickness is bounded to be no larger +than a half of set_visc_CS\%hbbl . + +\todo Channel drag needs to be explained + +A BBL viscosity is calculated so that the no-slip boundary condition in the vertical +viscosity solver implies the stress \f$\mathbf{\tau}_b\f$. + +*/ From 7f76390159b0fd4d5221ac3257c695bc7ba1ab5f Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 31 Aug 2021 17:17:14 -0800 Subject: [PATCH 11/15] Some background diffusivity text --- docs/zotero.bib | 25 ++++ .../vertical/_V_diffusivity.dox | 107 +++++++++++++++++- 2 files changed, 130 insertions(+), 2 deletions(-) diff --git a/docs/zotero.bib b/docs/zotero.bib index 8a5c14dcfa..a00fe569bd 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -1462,6 +1462,18 @@ @article{harrison2008 pages = {1894--1912} } +@article{danabasoglu2012, + doi = {10.1175/jcli-d-11-00091.1}, + year = 2012, + publisher = {American Meteorological Society}, + volume = {25}, + number = {5}, + pages = {1361--1389}, + author = {G. Danabasoglu and S. C. Bates and B. P. Briegleb and S. R. Jayne and M. Jochum and W. G. Large and S. Peacock and S. G. Yeager}, + title = {The {CCSM}4 Ocean Component}, + journal = {J. Climate} +} + @article{henyey1986, title = {Energy and action flow through the internal wave field: {An} eikonal approach}, volume = {91}, @@ -2650,3 +2662,16 @@ @article{miles1961 pages = {496--508}, doi = {10.1017/S0022112061000305} } + +@article{bryan1979, + doi = {10.1029/jc084ic05p02503}, + year = 1979, + publisher = {American Geophysical Union ({AGU})}, + volume = {84}, + number = {C5}, + pages = {2503}, + author = {K. Bryan and L. J. Lewis}, + title = {A water mass model of the World Ocean}, + journal = {J. Geophys. Res.} +} + diff --git a/src/parameterizations/vertical/_V_diffusivity.dox b/src/parameterizations/vertical/_V_diffusivity.dox index 5c40768eaf..1d79f58997 100644 --- a/src/parameterizations/vertical/_V_diffusivity.dox +++ b/src/parameterizations/vertical/_V_diffusivity.dox @@ -76,8 +76,15 @@ length scale where \f$Q\f$ is the turbulent kinetic energy (TKE) per unit mass, and \f$\lambda\f$ is a nondimensional constant. The function \f$F(\mbox{Ri})\f$ is a function of the Richardson number that remains to be determined. As in \cite turner1986, there must be a critical -value of \f$\mbox{Ri}\f$ above which \f$F(\mbox{Ri}) = 0\f$. There -are two length scales: the width of the low Richardson number region +value of \f$\mbox{Ri}\f$ above which \f$F(\mbox{Ri}) = 0\f$. +For better agreement with observations in a law-of-the-wall configuration, +we modify \f$L_d\f$ to be \f$\min (\lambda L_b, L_z)\f$, where \f$L_z\f$ +is the distance to the nearest solid boundary. This can be understood by +considering \f$L_d\f$ to be the size of the largest turbulent eddies, +whether they are constrained by the stratification (through \f$L_b\f$) +or through the geometry (through \f$L_z\f$). + +There are two length scales: the width of the low Richardson number region as in \cite turner1986, and the buoyancy length scale, which is the length scale over which the TKE is affected by the stratification (see \cite jackson2008 for more details). In particular, the inclusion of a @@ -143,8 +150,104 @@ reflects shear-driven turbulent mixing only; the total diffusivity would be this value plus any diffusivities due to other turbulent processes or a background value. +Based on \cite turner1986, we choose \f$F(\mbox{Ri})\f$ of the form + +\f[ + F(\mbox{Ri}) = F_0 \left( \frac{1 - \mbox{Ri} / \mbox{Ri}_c} + {1 + \alpha \mbox{Ri} / \mbox{Ri}_c} \right) , +\f] + +where \f$\alpha\f$ is the curvature parameter. This table shows the default +values of the relevant parameters: + + + +
Shear mixing parameters
Parameter Default value MOM6 parameter +
\f$\mbox{Ri}_c\f$ 0.25 RINO_CRIT +
\f$\nu_0\f$ \f$1.5 \times 10^{-5}\f$ KD_KAPPA_SHEAR_0 +
\f$F_0\f$ 0.089 SHEARMIX_RATE +
\f$\alpha\f$ -0.97 FRI_CURVATURE +
\f$\lambda\f$ 0.82 KAPPA_BUOY_SCALE_COEF +
\f$c_N\f$ 0.24 TKE_N_DECAY_CONST +
\f$c_S\f$ 0.14 TKE_SHEAR_DECAY_CONST +
+ +These can all be adjusted at run time, plus some other parameters such as the maximum number of iterations +to perform. + \section section_Background Background Mixing +There are three choices for the vertical background mixing: that in +CVMix (\cite bryan1979), that in \cite harrison2008, and that in +\cite danabasoglu2012. + +\subsection subsection_bryan_lewis CVMix background mixing + +The background vertical mixing in \cite bryan1979 is of the form: + +\f[ + \kappa = C_1 + C_2 \mbox{atan} [ C_3 ( |z| - C_4 )] +\f] + +where the contants are runtime parameters as shown here: + + + +
Bryan Lewis parameters
Parameter Units MOM6 parameter +
\f$C_1\f$ m2 s-1 BRYAN_LEWIS_C1 +
\f$C_2\f$ m2 s-1 BRYAN_LEWIS_C2 +
\f$C_3\f$ m-1 BRYAN_LEWIS_C3 +
\f$C_4\f$ m BRYAN_LEWIS_C4 +
+ +\subsection subsection_henyey Henyey IGW background mixing + +\cite harrison2008 choose a vertical background mixing with a latitudinal +dependence based on \cite henyey1986. Specifically, theory predicts +a minimum in mixing due to wave-wave interactions at the equator and +observations support that theory. In this option, the surface background +diffusivity is + +\f[ + \kappa_s (\phi) = \max \left[ 10^{-7}, \kappa_0 \left| \frac{f}{f_{30}} \right| + \frac{ \cosh^{-1} (1/f) }{ \cosh^{-1} (1/f_{30})} \right] , +\f] + +where \f$f_{30}\f$ is the Coriolis frequency at \f$30^\circ\f$ latitude. The two-dimensional equation for +the diffusivity is + +\f[ + \kappa(\phi, z) = \kappa_s + \Gamma \mbox{atan} \left( \frac{H_t}{\delta_t} \right) + + \Gamma \mbox{atan} \left( \frac{z - H_t}{\delta_t} \right) , +\f] +\f[ + \Gamma = \frac{(\kappa_d - \kappa_s) }{\left[ 0.5 \pi + \mbox{atan} \left( \frac{H_t}{\delta_t} \right) + \right] }, +\f] + +where \f$H_t = 2500\, \mbox{m}\f$, \f$\delta_t = 222\, \mbox{m}\f$, and +\f$\kappa_d\f$ is the deep ocean diffusivity of \f$10^{-4}\, \mbox{m}^2 +\, \mbox{s}^{-1}\f$. + +There is also a "new" Henyey version, taking into account the effect of stratification on +TKE dissipation, + +\f[ + \epsilon = \epsilon_0 \frac{f}{f_0} \frac{\mbox{acosh} (N/f)}{\mbox{acosh} (N_0 / f_0)} +\f] + +where \f$N_0\f$ and \f$f_0\f$ are the reference buoyancy frequency and inertial frequencies, respectively +and \f$\epsilon_0\f$ is the reference dissipation at \f$(N_0, f_0)\f$. In the previous version, \f$N = +N_0\f$. Additionally, the relationship between diapycnal diffusivities and stratification is included: + +\f[ + \kappa = \frac{\epsilon}{N^2} +\f] +This approach assumes that work done against gravity is uniformly distributed throughout the water column. +The original version concentrates buoyancy work in regions of strong stratification. + +\subsection subsection_danabasoglu_back Danabasoglu background mixing + \section section_Double_Diff Double Diffusion */ From 95cecb61280a123f975b403a64558df1b3cb2750 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 17 Sep 2021 11:57:47 -0800 Subject: [PATCH 12/15] Done with internal diffusion? --- docs/parameterizations_vertical.rst | 5 +- docs/zotero.bib | 9 + .../vertical/_Internal_tides.dox | 227 ------------ .../vertical/_V_diffusivity.dox | 346 +++++++++++++++++- 4 files changed, 347 insertions(+), 240 deletions(-) delete mode 100644 src/parameterizations/vertical/_Internal_tides.dox diff --git a/docs/parameterizations_vertical.rst b/docs/parameterizations_vertical.rst index c9404c5088..4705cf6c48 100644 --- a/docs/parameterizations_vertical.rst +++ b/docs/parameterizations_vertical.rst @@ -23,13 +23,12 @@ Interior and bottom-driven mixing Kappa-shear MOM_kappa_shear implements the shear-driven mixing of :cite:`jackson2008`. - :ref:`Internal_Shear_Mixing` - Internal-tide driven mixing The schemes of :cite:`st_laurent2002`, :cite:`polzin2009`, and :cite:`melet2012`, are all implemented through MOM_set_diffusivity and MOM_diabatic_driver. - :ref:`Internal_Tidal_Mixing` + :ref:`Internal_Vert_Mixing` + Vertical friction ----------------- diff --git a/docs/zotero.bib b/docs/zotero.bib index a00fe569bd..13ed6fd6eb 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2675,3 +2675,12 @@ @article{bryan1979 journal = {J. Geophys. Res.} } +@techreport{griffies2015a, + author = {S. M. Griffies and M. Levy and A. J. Adcroft and G. Danabasoglu and R. + W. Hallberg and D. Jacobsen and W. Large and T. Ringler}, + title = {Theory and Numerics of the Community Ocean Vertical Mixing (CVMix) + Project}, + year = {2015}, + pages = {98 pp}, + institution = {NOAA GFDL} +} diff --git a/src/parameterizations/vertical/_Internal_tides.dox b/src/parameterizations/vertical/_Internal_tides.dox deleted file mode 100644 index a07663d4a1..0000000000 --- a/src/parameterizations/vertical/_Internal_tides.dox +++ /dev/null @@ -1,227 +0,0 @@ -/*! \page Internal_Tidal_Mixing Internal Tidal Mixing - -Two parameterizations of vertical mixing due to internal tides are -available with the option INT_TIDE_DISSIPATION. The first is that of -\cite st_laurent2002 while the second is that of \cite polzin2009. Choose -between them with the INT_TIDE_PROFILE option. There are other relevant -parameters which can be seen in MOM_parameter_doc.all once the main tidal -dissipation switch is turned on. - -\section section_st_laurent St Laurent et al. - -The estimated turbulent dissipation rate of -internal tide energy \f$\epsilon\f$ is: - -\f[ - \epsilon = \frac{q E(x,y)}{\rho} F(z). -\f] - -where \f$\rho\f$ is the reference density of seawater, \f$E(x,y)\f$ is -the energy flux per unit area transferred from barotropic to baroclinic -tides, \f$q\f$ is the fraction of the internal-tide energy dissipated -locally, and \f$F(z)\f$ is the vertical structure of the dissipation. -This \f$q\f$ is estimated to be roughly 0.3 based on observations. The -term \f$E(x,y)\f$ is given by \cite st_laurent2002 as: - -\f[ - E(x,y) \simeq \frac{1}{2} \rho N_b \kappa h^2 \langle U^2 \rangle -\f] - -where \f$\rho\f$ is the reference density of seawater, \f$N_b\f$ is -the buoyancy frequency along the seafloor, and \f$(\kappa, h)\f$ are -the wavenumber and amplitude scales for the topographic roughness, and -\f$\langle U^2 \rangle\f$ is the barotropic tide variance. It is assumed -that the model will read in topographic roughness squared \f$h^2\f$ -from a file (the variable must be named "h2"). - -To convert from energy dissipation to vertical diffusion \f$K_d\f$, -the simple estimate is: - -\f[ - K_d \approx \frac{\Gamma q E(x,y) F(z)}{\rho N^2} -\f] - -where \f$\Gamma\f$ is the mixing efficiency, generally set to 0.2 -and \f$F(z)\f$ is a vertical structure function with exponential decay -from the bottom: - -\f[ - F(z) = \frac{e^{-(H+z)/\zeta}}{\zeta (1 - e^{H/\zeta}}. -\f] - -Here, \f$\zeta\f$ is a vertical decay scale with a default of 500 meters. -One change in MOM6 from the St. Laurent scheme is to use this form of \f$\Gamma\f$: - -\f[ - \Gamma = 0.2 \frac{N^2}{N^2 + \Omega^2} -\f] - -instead of \f$\Gamma = 0.2\f$, where \f$\Omega\f$ is the angular velocity -of the Earth. This allows the buoyancy fluxes to tend to zero in regions -of very weak stratification, allowing a no-flux bottom boundary condition -to be satisfied. - -This \f$K_d\f$ gets added to the diffusivity due to the background and -other contributions unless you set BBL_MIXING_AS_MAX to True, in which -case the maximum of all the contributions is used. - -\section section_polzin Polzin - -The vertical diffusion profile of \cite polzin2009 is a WKB-stretched -algebraic decay profile. It is based on a radiation balance equation, -which links the dissipation profile associated with internal breaking to -the finescale internal wave shear producing that dissipation. The vertical -profile of internal-tide driven energy dissipation can then vary in time -and space, and evolve in a changing climate (\cite melet2012). \cite melet2012 -describes how the Polzin scheme is implemented in MOM6, -copied here. - -The parameterization of \cite polzin2009 links the energy dissipation -profile to the finescale internal wave shear producing that -dissipation, using an idealized vertical wavenumber energy spectrum -to identify analytic solutions to a radiation balance equation -(\cite polzin2004). These solutions yield a dissipation profile -\f$\epsilon(z)\f$: - -\f[ - \epsilon = \frac{\epsilon_0}{[1 + (z/z_p)]^2}, -\f] - -where the magnitude \f$\epsilon_0\f$ and scale height \f$z_p\f$ can be expressed in terms of the -spectral amplitude and bandwidth of the idealized vertical wavenumber energy spectrum in uniform -stratification (\cite polzin2009). - -To take into account the nonuniform stratification, \cite polzin2009 applied a buoyancy scaling -using the Wentzel-Kramers-Brillouin (WKB) approximation. As a result, the vertical wavenumber of a -wave packet varies in proportion to the buoyancy frequency \f$N\f$, which in turn implies an -additional transport of energy to smaller scales, and thus a possible enhanced mixing in regions of -strong stratification. Such effects can be described by buoyancy scaling the vertical coordinate -\f$z\f$ as - -\f[ - z^{\ast}(z) = \int_{0}^{z} \left[ \frac{N^2 (z^\prime )}{N_b^2} \right] dz^{\prime} , -\f] - -with \f$z^\prime\f$ being positive upward relative to the bottom of the ocean. The turbulent -dissipation rate then becomes - -\f[ - \epsilon = \frac{\epsilon_0}{[1 + (z^{\ast} /z_p)]^2} \frac{N^2(z)}{N_b^2} . -\f] - -The spectral amplitude and bandwidth of the idealized vertical wavenumber -energy spectrum are identified after WKB scaling using a quasi-linear -spectral model of internal-tide generation that incorporates horizontal -advection of the barotropic tide into the momentum equation (\cite bell1975). -As a result, Polzin's formulation leads to an expression for -the spatially and temporally varying dissipation of internal tide energy -at the bottom \f$\epsilon_0\f$, and the vertical scale of decay for the -dissipation of internal tide energy \f$z_p\f$. - -\subsection subsection_energy_conserving Energy-conserving form - -To satisfy energy conservation (the integral of the vertical structure for the turbulent dissipation -over depth should be unity), the dissipation is rewritten as - -\f[ - \epsilon = \frac{\epsilon_0 z_p}{1 + (z^\ast/z_p)]^2} \frac{N^2(z)}{N^2_b} \left[ - \frac{1}{z^{\ast(z=H)}} + \frac{1}{z_p} \right] . -\f] - -In the MOM6 implementation, we use the \cite st_laurent2002 template for the vertical flux of energy -at the ocean floor, so that in both formulations: - -\f[ - \int_{0}^{H} \epsilon (z) dz = \frac{qE}{\rho} . -\f] - -Whereas \cite polzin2009 assumed that the total dissipation was locally in balance with the -barotropic to baroclinic energy conversion rate \f$(q=1)\f$, here we use the \cite simmons2004 value -of \f$q=1/3\f$ to retain as much consistency as possible between both parameterizations. - -\subsection subsection_vertical_decay_scale Vertical decay-scale reformulation - -We follow the \cite polzin2009 prescription for the vertical scale of -decay for the dissipation of internal-tide energy. However, we assume -that the topographic power law, denoted as \f$\nu\f$ in \cite polzin2009, -is equal to 1 (instead of 0.9) and we reformulated the expression of -\f$z_p\f$ to put it in a more readable form: - -\f[ - z_p = \frac{z_p^\mbox{ref} (\kappa^\mbox{ref})^2 (h^\mbox{ref})^2 (N_b^\mbox{ref})^3} {U^\mbox{ref}} - \frac{U}{h^2 \kappa^2 N_b^3} . -\f] - -The superscript ref refers to reference values of the various parameters, as given by -observations from the Brazil basin. Therefore, the above can be rewritten as - -\f[ - z_p = \mu (N_b^\mbox{ref} )^2 - \frac{U}{h^2 \kappa^2 N_b^3} . -\f] - -where \f$\mu\f$ is a nondimensional constant \f$(\mu = 0.06970)\f$ and \f$N_b^\mbox{ref} = 9.6 \times -10^{-4} s^{-1}\f$. Finally, a minimum decay scale of \f$z_p = 100 m\f$ is imposed in our -implementation. - -\subsection subsection_reformulation_WKB Reformulation of the WKB scaling - -Since the dissipation is expressed as a function of the ratio \f$z^\ast / z_p\f$, a different WKB -scaling can be used so long as we modify \f$z_p\f$ accordingly. In the implemented parameterization, -we define the scaled height coordinate \f$z^\ast\f$ by - -\f[ - z^\ast (z) = \frac{1}{\overline{N^2 (z)}^z} \int_{0}^{z} N^2(z^\prime ) dz ^\prime , -\f] - -with \f$z^\prime\f$ defined to be the height above the ocean bottom. By normalizing \f$N^2\f$ by its -vertical mean \f$\overline{N^2}^z\f$, \f$z^\ast\f$ ranges from \f$0\f$ to \f$H\f$, the depth of the -ocean. - -The WKB-scaled vertical decay scale for the Polzin formulation becomes - -\f[ - z^\ast_p = \mu(N_b^\mbox{ref})^2 \frac{U}{h^2 \kappa^2 N_b \overline{N^2}^z} . -\f] - -Unlike the \cite st_laurent2002 parameterization, the vertical decay scale now depends on physical -variables and can evolve with a changing climate. - -Finally, the Polzin vertical profile of dissipation implemented in the model is given by - -\f[ - \epsilon = \frac{qE(x,y)}{\rho [1 + (z^\ast/z_p^\ast)]^2} \frac{N^2(z)}{\overline{N^2}^z} - \left( \frac{1}{H} + \frac{1}{z_p^\ast} \right) . -\f] - -In both parameterizations, turbulent diapycnal diffusivities are inferred from the dissipation -\f$\epsilon\f$ by: - -\f[ - K_d = \frac{\Gamma \epsilon}{N^2} -\f] - -and using this form of \f$\Gamma\f$: - -\f[ - \Gamma = 0.2 \frac{N^2}{N^2 + \Omega^2} -\f] - -instead of \f$\Gamma = 0.2\f$, where \f$\Omega\f$ is the angular velocity -of the Earth. This allows the buoyancy fluxes to tend to zero in regions -of very weak stratification, allowing a no-flux bottom boundary condition -to be satisfied. - -\section Nikurashin Lee Wave Mixing - -If one has the INT_TIDE_DISSIPATION flag on, there is an option to also turn on -LEE_WAVE_DISSIPATION. The theory is presented in \cite nikurashin2010a -while the application of it is presented in \cite nikurashin2010b. For -the implementation in MOM6, it is required that you provide an estimate -of the TKE loss due to the Lee waves which is then applied with either -the St. Laurent or the Polzin vertical profile. - -\todo Is there a script to produce this somewhere or what??? - -*/ - diff --git a/src/parameterizations/vertical/_V_diffusivity.dox b/src/parameterizations/vertical/_V_diffusivity.dox index 4d671fec88..8c4c8ce7aa 100644 --- a/src/parameterizations/vertical/_V_diffusivity.dox +++ b/src/parameterizations/vertical/_V_diffusivity.dox @@ -1,12 +1,12 @@ -/*! \page Internal_Shear_Mixing Internal Vertical Mixing +/*! \page Internal_Vert_Mixing Internal Vertical Mixing Sets the interior vertical diffusion of scalars due to the following processes: --# Shear-driven mixing: two options, \cite jackson2008 and KPP interior; --# Background mixing via CVMix (Bryan-Lewis profile), the scheme described by - \cite harrison2008, or that in \cite danabasoglu2012. --# Double-diffusion, old method and new method via CVMix; --# Tidal mixing: many options available, see \ref Internal_Tidal_Mixing. +-# Shear-driven mixing (\ref section_Shear): \cite jackson2008 or KPP interior; +-# Background mixing (\ref section_Background): via CVMix (Bryan-Lewis profile), + the scheme described by \cite harrison2008, or that in \cite danabasoglu2012. +-# Double-diffusion (\ref section_Double_Diff): old method or new method via CVMix; +-# Tidal mixing: many options available, see \ref section_Internal_Tidal_Mixing. In addition, the MOM_set_diffusivity has the option to set the interior vertical viscosity associated with processes 1,2 and 4 listed above, which is stored in @@ -50,13 +50,14 @@ parameterization of \cite large1994 is as follows, where the diffusivity \f$\kap is given by \f[ - \kappa = \kappa_0 \left[ 1 - \min \left( 1, \frac{\mbox{Ri}}{\mbox{Ri}_c} \right) ^2 \right] ^3 ,\ + \kappa = \kappa_0 \left[ 1 - \min \left( 1, \frac{\mbox{Ri}}{\mbox{Ri}_c} \right) + ^2 \right] ^3 , \f] with \f$\kappa_0 = 5 \times 10^{-3}\, \mbox{m}^2 \,\mbox{s}^{-1}\f$ and \f$\mbox{Ri}_c = 0.7\f$. One can instead select the \cite pacanowski1981 scheme within CVMix. Unlike -the \cite large1994 scheme, they propose that the\ vertical shear +the \cite large1994 scheme, they propose that the vertical shear viscosity \f$\nu_{\mbox{shear}}\f$ be different from the vertical shear diffusivity \f$\kappa_{\mbox{shear}}\f$. For gravitationally stable profiles (i.e., \f$N^2 > 0\f$), they chose @@ -248,13 +249,16 @@ the diffusivity is where \f$H_t = 2500\, \mbox{m}\f$, \f$\delta_t = 222\, \mbox{m}\f$, and \f$\kappa_d\f$ is the deep ocean diffusivity of \f$10^{-4}\, \mbox{m}^2 \, \mbox{s}^{-1}\f$. Note that this is the vertical structure described -in \cite harrison2008, but that isn't what is in the code. Instead, the surface +in \cite harrison2008, but that isn't what is in the MOM6 code. Instead, the surface value is propagated down, with the assumption that the tidal mixing parameterization -will provide the deep mixing: \ref Internal_Tidal_Mixing. +will provide the deep mixing: \ref section_Internal_Tidal_Mixing. There is also a "new" Henyey version, taking into account the effect of stratification on TKE dissipation, +\todo Harrison (personal communication) recommends that this option be made obsolete and +eventually removed. + \f[ \epsilon = \epsilon_0 \frac{f}{f_0} \frac{\mbox{acosh} (N/f)}{\mbox{acosh} (N_0 / f_0)} \f] @@ -281,4 +285,326 @@ Some parameters of this curve are set in the input file, some are hard-coded in \section section_Double_Diff Double Diffusion +From \cite large1994, \cite griffies2015a, double-diffusive mixing +can occur when the vertical gradient of density is stable but the +vertical gradient of either salinity or temperature is unstable in its +contribution to density. The key stratification parameter for double +diffusive processes is + +\f[ + R_\rho = \frac{\alpha}{\beta} \left( \frac{\partial \Theta / \partial z}{\partial S / + \partial z} \right) , +\f] + +where the thermal expansion coefficient is given by + +\f[ + \alpha = - \frac{1}{\rho} \left( \frac{\partial \rho}{\partial \Theta} \right) , +\f] + +and the haline contraction coefficient is + +\f[ + \beta = \frac{1}{\rho} \left( \frac{\partial \rho}{\partial S} \right) . +\f] + +Note that the effects from double diffusive processes on viscosity are not well known and +are ignored in MOM6. + +In MOM6, there are two choices for the implementation of double +diffusion. The older DOUBLE_DIFFUSION option, with reference to an +unknown tech report from NCAR, aims to match the scheme used in MOM4, an update on +\cite large1994. The newer option is to call the routines from CVMix (USE_CVMIX_DDIFF). + +There are two regimes of double diffusive processes, salt fingering and diffusive +convective, with differing parameterizations in the two regimes. + +\subsection subsection_salt_finger Salt fingering regime + +The salt fingering regime occurs when salinity is destabilizing the water column (salty, +above fresh water) and when the stratification parameter \f$R_\rho\f$ is within a +particular range: + +\f[ + \frac{\partial S}{\partial z} > 0 +\f] +\f[ + 1 < R_\rho < R_\rho^0. +\f] + +The value of the cutoff \f$R_\rho\f$ is 1.9 in the old code, 2.55 in CVMix. + +The form of the diffusivity for both is + +\f{eqnarray}{ + \kappa_d =& \kappa_d^0 \left[ 1 - \left( \frac{R_\rho - 1}{R_\rho^0 - 1} \right) + \right]^3 & \mbox{for } 1 < R_\rho < R_\rho^0 \\ + \kappa_d =& 0 & \mbox{otherwise.} +\f} + +The default values of \f$\kappa_d^0\f$ are + +\f{eqnarray}{ + \kappa_d^0 =& 1 \times 10^{-4} & \mbox{for salinity and other tracers} \\ + \kappa_d^0 =& 0.7 \times 10^{-4} & \mbox{for temperature.} +\f} + +Note that the form in \cite large1994 is slightly different. + +\subsection subsection_diffusive_convective Diffusive convective regime + +Both implementations of the diffusive convective double diffusion have the same form +(\cite large1994) and are active when + +\f[ + \frac{\partial \Theta}{\partial z} < 0 +\f] +\f[ + 0 < R_\rho < 1. +\f] + +For temperature, the vertical diffusivity is given by + +\f[ + \kappa_d = \nu_\mbox{molecular} \times 0.909 \exp \left( 4.6 \exp \left[ -.54 + \left( R_\rho^{-1} - 1 \right) \right] \right) , +\f] + +where + +\f[ + \nu_\mbox{molecular} = 1.5 \times 10^{-6} \mbox{m}^2 \mbox{s}^{-1} +\f] + +is the molecular viscosity of water. Multiplying the diffusivity by the Prandtl number +\f$Pr\f$ + +\f{eqnarray}{ + Pr = \left\{ \begin{matrix} (1.85 - 0.85 R_\rho^{-1} ) R_\rho & 0.5 \leq R_\rho < 1 \\ + 0.15 R_\rho & R_\rho < 0.5 , \end{matrix} \right. +\f} + +gives the diffusivity for salinity and other tracers. + +\section section_Internal_Tidal_Mixing Internal Tidal Mixing + +Two parameterizations of vertical mixing due to internal tides are +available with the option INT_TIDE_DISSIPATION. The first is that of +\cite st_laurent2002 while the second is that of \cite polzin2009. Choose +between them with the INT_TIDE_PROFILE option. There are other relevant +parameters which can be seen in MOM_parameter_doc.all once the main tidal +dissipation switch is turned on. + +\subsection subsection_st_laurent St Laurent et al. + +The estimated turbulent dissipation rate of +internal tide energy \f$\epsilon\f$ is: + +\f[ + \epsilon = \frac{q E(x,y)}{\rho} F(z). +\f] + +where \f$\rho\f$ is the reference density of seawater, \f$E(x,y)\f$ is +the energy flux per unit area transferred from barotropic to baroclinic +tides, \f$q\f$ is the fraction of the internal-tide energy dissipated +locally, and \f$F(z)\f$ is the vertical structure of the dissipation. +This \f$q\f$ is estimated to be roughly 0.3 based on observations. The +term \f$E(x,y)\f$ is given by \cite st_laurent2002 as: + +\f[ + E(x,y) \simeq \frac{1}{2} \rho N_b \kappa h^2 \langle U^2 \rangle +\f] + +where \f$\rho\f$ is the reference density of seawater, \f$N_b\f$ is +the buoyancy frequency along the seafloor, and \f$(\kappa, h)\f$ are +the wavenumber and amplitude scales for the topographic roughness, and +\f$\langle U^2 \rangle\f$ is the barotropic tide variance. It is assumed +that the model will read in topographic roughness squared \f$h^2\f$ +from a file (the variable must be named "h2"). + +To convert from energy dissipation to vertical diffusion \f$K_d\f$, +the simple estimate is: + +\f[ + K_d \approx \frac{\Gamma q E(x,y) F(z)}{\rho N^2} +\f] + +where \f$\Gamma\f$ is the mixing efficiency, generally set to 0.2 +and \f$F(z)\f$ is a vertical structure function with exponential decay +from the bottom: + +\f[ + F(z) = \frac{e^{-(H+z)/\zeta}}{\zeta (1 - e^{H/\zeta}}. +\f] + +Here, \f$\zeta\f$ is a vertical decay scale with a default of 500 meters. +One change in MOM6 from the St. Laurent scheme is to use this form of \f$\Gamma\f$: + +\f[ + \Gamma = 0.2 \frac{N^2}{N^2 + \Omega^2} +\f] + +instead of \f$\Gamma = 0.2\f$, where \f$\Omega\f$ is the angular velocity +of the Earth. This allows the buoyancy fluxes to tend to zero in regions +of very weak stratification, allowing a no-flux bottom boundary condition +to be satisfied. + +\subsection subsection_polzin Polzin + +The vertical diffusion profile of \cite polzin2009 is a WKB-stretched +algebraic decay profile. It is based on a radiation balance equation, +which links the dissipation profile associated with internal breaking to +the finescale internal wave shear producing that dissipation. The vertical +profile of internal-tide driven energy dissipation can then vary in time +and space, and evolve in a changing climate (\cite melet2012). \cite melet2012 +describes how the Polzin scheme is implemented in MOM6, +copied here. + +The parameterization of \cite polzin2009 links the energy dissipation +profile to the finescale internal wave shear producing that +dissipation, using an idealized vertical wavenumber energy spectrum +to identify analytic solutions to a radiation balance equation +(\cite polzin2004). These solutions yield a dissipation profile +\f$\epsilon(z)\f$: + +\f[ + \epsilon = \frac{\epsilon_0}{[1 + (z/z_p)]^2}, +\f] + +where the magnitude \f$\epsilon_0\f$ and scale height \f$z_p\f$ can be expressed in terms of the +spectral amplitude and bandwidth of the idealized vertical wavenumber energy spectrum in uniform +stratification (\cite polzin2009). + +To take into account the nonuniform stratification, \cite polzin2009 applied a buoyancy scaling +using the Wentzel-Kramers-Brillouin (WKB) approximation. As a result, the vertical wavenumber of a +wave packet varies in proportion to the buoyancy frequency \f$N\f$, which in turn implies an +additional transport of energy to smaller scales, and thus a possible enhanced mixing in regions of +strong stratification. Such effects can be described by buoyancy scaling the vertical coordinate +\f$z\f$ as + +\f[ + z^{\ast}(z) = \int_{0}^{z} \left[ \frac{N^2 (z^\prime )}{N_b^2} \right] dz^{\prime} , +\f] + +with \f$z^\prime\f$ being positive upward relative to the bottom of the ocean. The turbulent +dissipation rate then becomes + +\f[ + \epsilon = \frac{\epsilon_0}{[1 + (z^{\ast} /z_p)]^2} \frac{N^2(z)}{N_b^2} . +\f] + +The spectral amplitude and bandwidth of the idealized vertical wavenumber +energy spectrum are identified after WKB scaling using a quasi-linear +spectral model of internal-tide generation that incorporates horizontal +advection of the barotropic tide into the momentum equation (\cite bell1975). +As a result, Polzin's formulation leads to an expression for +the spatially and temporally varying dissipation of internal tide energy +at the bottom \f$\epsilon_0\f$, and the vertical scale of decay for the +dissipation of internal tide energy \f$z_p\f$. + +\subsubsection subsection_energy_conserving Energy-conserving form + +To satisfy energy conservation (the integral of the vertical structure for the turbulent dissipation +over depth should be unity), the dissipation is rewritten as + +\f[ + \epsilon = \frac{\epsilon_0 z_p}{1 + (z^\ast/z_p)]^2} \frac{N^2(z)}{N^2_b} \left[ + \frac{1}{z^{\ast(z=H)}} + \frac{1}{z_p} \right] . +\f] + +In the MOM6 implementation, we use the \cite st_laurent2002 template for the vertical flux of energy +at the ocean floor, so that in both formulations: + +\f[ + \int_{0}^{H} \epsilon (z) dz = \frac{qE}{\rho} . +\f] + +Whereas \cite polzin2009 assumed that the total dissipation was locally in balance with the +barotropic to baroclinic energy conversion rate \f$(q=1)\f$, here we use the \cite simmons2004 value +of \f$q=1/3\f$ to retain as much consistency as possible between both parameterizations. + +\subsubsection subsection_vertical_decay_scale Vertical decay-scale reformulation + +We follow the \cite polzin2009 prescription for the vertical scale of +decay for the dissipation of internal-tide energy. However, we assume +that the topographic power law, denoted as \f$\nu\f$ in \cite polzin2009, +is equal to 1 (instead of 0.9) and we reformulated the expression of +\f$z_p\f$ to put it in a more readable form: + +\f[ + z_p = \frac{z_p^\mbox{ref} (\kappa^\mbox{ref})^2 (h^\mbox{ref})^2 (N_b^\mbox{ref})^3} {U^\mbox{ref}} + \frac{U}{h^2 \kappa^2 N_b^3} . +\f] + +The superscript ref refers to reference values of the various parameters, as given by +observations from the Brazil basin. Therefore, the above can be rewritten as + +\f[ + z_p = \mu (N_b^\mbox{ref} )^2 + \frac{U}{h^2 \kappa^2 N_b^3} . +\f] + +where \f$\mu\f$ is a nondimensional constant \f$(\mu = 0.06970)\f$ and \f$N_b^\mbox{ref} = 9.6 \times +10^{-4} s^{-1}\f$. Finally, a minimum decay scale of \f$z_p = 100 m\f$ is imposed in our +implementation. + +\subsubsection subsection_reformulation_WKB Reformulation of the WKB scaling + +Since the dissipation is expressed as a function of the ratio \f$z^\ast / z_p\f$, a different WKB +scaling can be used so long as we modify \f$z_p\f$ accordingly. In the implemented parameterization, +we define the scaled height coordinate \f$z^\ast\f$ by + +\f[ + z^\ast (z) = \frac{1}{\overline{N^2 (z)}^z} \int_{0}^{z} N^2(z^\prime ) dz ^\prime , +\f] + +with \f$z^\prime\f$ defined to be the height above the ocean bottom. By normalizing \f$N^2\f$ by its +vertical mean \f$\overline{N^2}^z\f$, \f$z^\ast\f$ ranges from \f$0\f$ to \f$H\f$, the depth of the +ocean. + +The WKB-scaled vertical decay scale for the Polzin formulation becomes + +\f[ + z^\ast_p = \mu(N_b^\mbox{ref})^2 \frac{U}{h^2 \kappa^2 N_b \overline{N^2}^z} . +\f] + +Unlike the \cite st_laurent2002 parameterization, the vertical decay scale now depends on physical +variables and can evolve with a changing climate. + +Finally, the Polzin vertical profile of dissipation implemented in the model is given by + +\f[ + \epsilon = \frac{qE(x,y)}{\rho [1 + (z^\ast/z_p^\ast)]^2} \frac{N^2(z)}{\overline{N^2}^z} + \left( \frac{1}{H} + \frac{1}{z_p^\ast} \right) . +\f] + +In both parameterizations, turbulent diapycnal diffusivities are inferred from the dissipation +\f$\epsilon\f$ by: + +\f[ + K_d = \frac{\Gamma \epsilon}{N^2} +\f] + +and using this form of \f$\Gamma\f$: + +\f[ + \Gamma = 0.2 \frac{N^2}{N^2 + \Omega^2} +\f] + +instead of \f$\Gamma = 0.2\f$, where \f$\Omega\f$ is the angular velocity +of the Earth. This allows the buoyancy fluxes to tend to zero in regions +of very weak stratification, allowing a no-flux bottom boundary condition +to be satisfied. + +\subsection subsection_Lee_waves Nikurashin Lee Wave Mixing + +If one has the INT_TIDE_DISSIPATION flag on, there is an option to also turn on +LEE_WAVE_DISSIPATION. The theory is presented in \cite nikurashin2010a +while the application of it is presented in \cite nikurashin2010b. For +the implementation in MOM6, it is required that you provide an estimate +of the TKE loss due to the Lee waves which is then applied with either +the St. Laurent or the Polzin vertical profile. + +\todo Is there a script to produce this somewhere or what??? + */ From b723674fbb799e36aedd8b0461845ad095fe7268 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Mon, 20 Sep 2021 13:53:31 -0800 Subject: [PATCH 13/15] Adding information to MOM6 warning. - where is the water depth negative? --- src/core/MOM_barotropic.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 471999c60c..611cf8aea6 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2359,7 +2359,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (GV%Boussinesq) then do j=js,je ; do i=is,ie if ((eta(i,j) < -GV%Z_to_H*G%bathyT(i,j)) .and. (G%mask2dT(i,j) > 0.0)) then - write(mesg,'(ES24.16," vs. ",ES24.16)') GV%H_to_m*eta(i,j), -US%Z_to_m*G%bathyT(i,j) + write(mesg,'(ES24.16," vs. ",ES24.16, " at ", i7, i7)') GV%H_to_m*eta(i,j), -US%Z_to_m*G%bathyT(i,j), & + i + G%isd_global, j + G%jsd_global if (err_count < 2) & call MOM_error(WARNING, "btstep: eta has dropped below bathyT: "//trim(mesg), all_print=.true.) err_count = err_count + 1 From e368bfe959afc1331d2b1917612b8c45fb748344 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 21 Sep 2021 19:22:49 -0800 Subject: [PATCH 14/15] Several small things, including fix to sponge verbosity. - Also, added a link to OBC wiki page. - Working around bibtex hashing issue I don't understand, renaming some unused tags. - Making MOM_barotropic when the water depth goes negative. --- docs/parameterizations_lateral.rst | 2 ++ docs/zotero.bib | 4 ++-- src/core/MOM_barotropic.F90 | 4 ++-- src/framework/MOM_horizontal_regridding.F90 | 8 ++++++-- 4 files changed, 12 insertions(+), 6 deletions(-) diff --git a/docs/parameterizations_lateral.rst b/docs/parameterizations_lateral.rst index 102090b7a4..3a3266a2bb 100644 --- a/docs/parameterizations_lateral.rst +++ b/docs/parameterizations_lateral.rst @@ -43,4 +43,6 @@ Tidal forcing ------------- Astronomical tidal forcings and self-attraction and loading are implement in MOM_tidal_forcing. +Tides can also be added via an open boundary tidal specification, +see [OBC wiki page](https://github.com/NOAA-GFDL/MOM6-examples/wiki/Open-Boundary-Conditions). diff --git a/docs/zotero.bib b/docs/zotero.bib index 13ed6fd6eb..bb400542b8 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -1728,7 +1728,7 @@ @article{adcroft2006 pages = {224--233} } -@article{adcroft2004, +@article{adcroft2004-1, title = {Rescaled height coordinates for accurate representation of free-surface flows in ocean circulation models}, volume = {7}, issn = {1463-5003}, @@ -2308,7 +2308,7 @@ @article{carpenter1990 doi = {https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2} } -@article{kasahara1974, +@article{kasahara1974-1, title = {Various {Vertical} {Coordinate} {Systems} {Used} for {Numerical} {Weather} {Prediction}}, volume = {102}, issn = {0027-0644}, diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 611cf8aea6..f3e37e73de 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2359,8 +2359,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (GV%Boussinesq) then do j=js,je ; do i=is,ie if ((eta(i,j) < -GV%Z_to_H*G%bathyT(i,j)) .and. (G%mask2dT(i,j) > 0.0)) then - write(mesg,'(ES24.16," vs. ",ES24.16, " at ", i7, i7)') GV%H_to_m*eta(i,j), -US%Z_to_m*G%bathyT(i,j), & - i + G%isd_global, j + G%jsd_global + write(mesg,'(ES24.16," vs. ",ES24.16, " at ", ES12.4, ES12.4, i7, i7)') GV%H_to_m*eta(i,j), & + -US%Z_to_m*G%bathyT(i,j), G%geoLonT(i,j), G%geoLatT(i,j), i + G%isd_global, j + G%jsd_global if (err_count < 2) & call MOM_error(WARNING, "btstep: eta has dropped below bathyT: "//trim(mesg), all_print=.true.) err_count = err_count + 1 diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index de2a76a746..cfdc8bc273 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -9,6 +9,7 @@ module MOM_horizontal_regridding use MOM_domains, only : pass_var use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint +use MOM_error_handler, only : MOM_get_verbosity use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_interpolate, only : time_interp_external, horiz_interp_init @@ -676,6 +677,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t real, dimension(SZI_(G),SZJ_(G)) :: good2 ! 1 where the data is valid after Ice-9 real, dimension(SZI_(G),SZJ_(G)) :: fill2 ! 1 for points that still need to be filled after Ice-9 integer :: turns + integer :: verbosity turns = G%HI%turns @@ -696,6 +698,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_value) + verbosity = MOM_get_verbosity() + id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3) spongeDataOngrid = .false. @@ -764,7 +768,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t if (.not.spongeDataOngrid) then if (is_root_pe()) & - call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) + call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>2), turns=turns) ! Loop through each data level and interpolate to model grid. ! After interpolating, fill in points which will be needed to define the layers. do k=1,kd @@ -880,7 +884,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t enddo ! kd else - call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) + call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>2), turns=turns) do k=1,kd do j=js,je do i=is,ie From a3a4d8efec70595b25478a4e4683dbfff200b0ff Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 22 Sep 2021 10:58:49 -0800 Subject: [PATCH 15/15] Changing verbosity in MOM_horizontal_regridding.F90 --- src/framework/MOM_horizontal_regridding.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index cfdc8bc273..ba59f50abe 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -768,7 +768,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t if (.not.spongeDataOngrid) then if (is_root_pe()) & - call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>2), turns=turns) + call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>5), turns=turns) ! Loop through each data level and interpolate to model grid. ! After interpolating, fill in points which will be needed to define the layers. do k=1,kd @@ -884,7 +884,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t enddo ! kd else - call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>2), turns=turns) + call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>5), turns=turns) do k=1,kd do j=js,je do i=is,ie