diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index d97cb96a40..a24dd03fd9 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -580,7 +580,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) if (debug .and. root_pe().eq.pe_here()) print *, "calling ocn_domain_mct" call ocn_domain_mct(lsize, MOM_MCT_gsmap, MOM_MCT_dom) - call ocn_domain_mct(lsize*km, MOM_MCT_gsmap3d, MOM_MCT_dom3d) !TODO: this is not used + !call ocn_domain_mct(lsize*km, MOM_MCT_gsmap3d, MOM_MCT_dom3d) !TODO: this is not used ! Inialize mct attribute vectors diff --git a/pkg/CVMix-src b/pkg/CVMix-src index d83f582714..534fc38e75 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit d83f582714e7f0f98d20efd8fac8fab01fa3bfe6 +Subproject commit 534fc38e759fcb8a2563fa0dc4c0bbf81f758db9 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 1fd75a71ef..5398da3e69 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1466,7 +1466,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & real, dimension(:,:), pointer :: shelf_area type(MOM_restart_CS), pointer :: restart_CSp_tmp => NULL() type(group_pass_type) :: tmp_pass_uv_T_S_h, pass_uv_T_S_h - type(group_pass_type) :: tmp_pass_Kv_turb + ! GMM, the following *is not* used. Should we delete it? + type(group_pass_type) :: tmp_pass_Kv_shear real :: default_val ! default value for a parameter logical :: write_geom_files ! If true, write out the grid geometry files. @@ -2308,8 +2309,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call do_group_pass(pass_uv_T_S_h, G%Domain) - if (associated(CS%visc%Kv_turb)) & - call pass_var(CS%visc%Kv_turb, G%Domain, To_All+Omit_Corners, halo=1) + if (associated(CS%visc%Kv_shear)) & + call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) call cpu_clock_end(id_clock_pass_init) @@ -2920,6 +2921,9 @@ subroutine MOM_end(CS) call tracer_registry_end(CS%tracer_Reg) call tracer_flow_control_end(CS%tracer_flow_CSp) + ! GMM, the following is commented because it fails on Travis. + !if (associated(CS%diabatic_CSp)) call diabatic_driver_end(CS%diabatic_CSp) + if (CS%offline_tracer_mode) call offline_transport_end(CS%offline_CSp) DEALLOC_(CS%uhtr) ; DEALLOC_(CS%vhtr) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index f0f3437f4b..f7fa45f12c 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -173,7 +173,7 @@ module MOM_variables !! coefficients, and related fields. type, public :: vertvisc_type real :: Prandtl_turb !< The Prandtl number for the turbulent diffusion - !! that is captured in Kd_turb. + !! that is captured in Kd_shear. real, pointer, dimension(:,:) :: & bbl_thick_u => NULL(), & !< The bottom boundary layer thickness at the !! u-points, in m. @@ -224,10 +224,13 @@ module MOM_variables ! Kd_extra_S is positive for salt fingering; Kd_extra_T ! is positive for double diffusive convection. These ! are only allocated if DOUBLE_DIFFUSION is true. - Kd_turb => NULL(), &!< The turbulent diapycnal diffusivity at the interfaces - !! between each layer, in m2 s-1. - Kv_turb => NULL(), &!< The turbulent vertical viscosity at the interfaces - !! between each layer, in m2 s-1. + Kd_shear => NULL(), &!< The shear-driven turbulent diapycnal diffusivity + !! at the interfaces between each layer, in m2 s-1. + Kv_shear => NULL(), &!< The shear-driven turbulent vertical viscosity + !! at the interfaces between each layer, in m2 s-1. + Kv_slow => NULL(), &!< The turbulent vertical viscosity component due to + !! "slow" processes (e.g., tidal, background, + !! convection etc). TKE_turb => NULL() !< The turbulent kinetic energy per unit mass defined !! at the interfaces between each layer, in m2 s-2. end type vertvisc_type diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 3bff792ad1..219a6b33fa 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -76,7 +76,7 @@ module MOM_state_initialization use Rossby_front_2d_initialization, only : Rossby_front_initialize_temperature_salinity use Rossby_front_2d_initialization, only : Rossby_front_initialize_velocity use SCM_idealized_hurricane, only : SCM_idealized_hurricane_TS_init -use SCM_CVmix_tests, only: SCM_CVmix_tests_TS_init +use SCM_CVMix_tests, only: SCM_CVMix_tests_TS_init use dyed_channel_initialization, only : dyed_channel_set_OBC_tracer_data use dyed_obcs_initialization, only : dyed_obcs_set_OBC_data use supercritical_initialization, only : supercritical_set_OBC_data @@ -338,7 +338,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & " \t dumbbell - sloshing channel ICs. \n"//& " \t rossby_front - a mixed layer front in thermal wind balance.\n"//& " \t SCM_ideal_hurr - used in the SCM idealized hurricane test.\n"//& - " \t SCM_CVmix_tests - used in the SCM CVmix tests.\n"//& + " \t SCM_CVMix_tests - used in the SCM CVMix tests.\n"//& " \t USER - call a user modified routine.", & fail_if_missing=new_sim, do_not_log=just_read) ! " \t baroclinic_zone - an analytic baroclinic zone. \n"//& @@ -371,7 +371,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & tv%S, h, G, GV, PF, eos, just_read_params=just_read) case ("SCM_ideal_hurr"); call SCM_idealized_hurricane_TS_init ( tv%T, & tv%S, h, G, GV, PF, just_read_params=just_read) - case ("SCM_CVmix_tests"); call SCM_CVmix_tests_TS_init (tv%T, & + case ("SCM_CVMix_tests"); call SCM_CVMix_tests_TS_init (tv%T, & tv%S, h, G, GV, PF, just_read_params=just_read) case ("dense"); call dense_water_initialize_TS(G, GV, PF, eos, tv%T, tv%S, & h, just_read_params=just_read) diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 new file mode 100644 index 0000000000..4b422ccf9a --- /dev/null +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -0,0 +1,270 @@ +!> Interface to CVMix convection scheme. +module MOM_CVMix_conv + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field +use MOM_diag_mediator, only : post_data +use MOM_EOS, only : calculate_density +use MOM_variables, only : thermo_var_ptrs +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_debugging, only : hchksum +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_file_parser, only : get_param, log_version, param_file_type +use CVMix_convection, only : CVMix_init_conv, CVMix_coeffs_conv +use CVMix_kpp, only : CVMix_kpp_compute_kOBL_depth + +implicit none ; private + +#include + +public CVMix_conv_init, calculate_CVMix_conv, CVMix_conv_end, CVMix_conv_is_used + +!> Control structure including parameters for CVMix convection. +type, public :: CVMix_conv_cs + + ! Parameters + real :: kd_conv_const !< diffusivity constant used in convective regime (m2/s) + real :: kv_conv_const !< viscosity constant used in convective regime (m2/s) + real :: bv_sqr_conv !< Threshold for squared buoyancy frequency + !! needed to trigger Brunt-Vaisala parameterization (1/s^2) + real :: min_thickness !< Minimum thickness allowed (m) + logical :: debug !< If true, turn on debugging + + ! Daignostic handles and pointers + type(diag_ctrl), pointer :: diag => NULL() + integer :: id_N2 = -1, id_kd_conv = -1, id_kv_conv = -1 + + ! Diagnostics arrays + real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) + real, allocatable, dimension(:,:,:) :: kd_conv !< Diffusivity added by convection (m2/s) + real, allocatable, dimension(:,:,:) :: kv_conv !< Viscosity added by convection (m2/s) + +end type CVMix_conv_cs + +character(len=40) :: mdl = "MOM_CVMix_conv" !< This module's name. + +contains + +!> Initialized the CVMix convection mixing routine. +logical function CVMix_conv_init(Time, G, GV, param_file, diag, CS) + + type(time_type), intent(in) :: Time !< The current time. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. + type(CVMix_conv_cs), pointer :: CS !< This module's control structure. + + ! Local variables + real :: prandtl_conv !< Turbulent Prandtl number used in convective instabilities. + logical :: useEPBL !< If True, use the ePBL boundary layer scheme. + +! This include declares and sets the variable "version". +#include "version_variable.h" + + if (associated(CS)) then + call MOM_error(WARNING, "CVMix_conv_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + ! Read parameters + call log_version(param_file, mdl, version, & + "Parameterization of enhanced mixing due to convection via CVMix") + call get_param(param_file, mdl, "USE_CVMix_CONVECTION", CVMix_conv_init, & + "If true, turns on the enhanced mixing due to convection \n"// & + "via CVMix. This scheme increases diapycnal diffs./viscs. \n"// & + " at statically unstable interfaces. Relevant parameters are \n"// & + "contained in the CVMix_CONVECTION% parameter block.", & + default=.false.) + + if (.not. CVMix_conv_init) return + + call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useEPBL, default=.false., & + do_not_log=.true.) + + ! Warn user if EPBL is being used, since in this case mixing due to convection will + ! be aplied in the boundary layer + if (useEPBL) then + call MOM_error(WARNING, 'MOM_CVMix_conv_init: '// & + 'CVMix convection may not be properly applied when ENERGETICS_SFC_PBL = True'//& + 'as convective mixing might occur in the boundary layer.') + endif + + call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) + + call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, do_not_log=.True.) + + call openParameterBlock(param_file,'CVMix_CONVECTION') + + call get_param(param_file, mdl, "PRANDTL_CONV", prandtl_conv, & + "The turbulent Prandtl number applied to convective \n"//& + "instabilities (i.e., used to convert KD_CONV into KV_CONV)", & + units="nondim", default=1.0) + + call get_param(param_file, mdl, 'KD_CONV', CS%kd_conv_const, & + "Diffusivity used in convective regime. Corresponding viscosity \n" // & + "(KV_CONV) will be set to KD_CONV * PRANDTL_TURB.", & + units='m2/s', default=1.00) + + call get_param(param_file, mdl, 'BV_SQR_CONV', CS%bv_sqr_conv, & + "Threshold for squared buoyancy frequency needed to trigger \n" // & + "Brunt-Vaisala parameterization.", & + units='1/s^2', default=0.0) + + call closeParameterBlock(param_file) + + ! set kv_conv_const based on kd_conv_const and prandtl_conv + CS%kv_conv_const = CS%kd_conv_const * prandtl_conv + + ! allocate arrays and set them to zero + allocate(CS%N2(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%N2(:,:,:) = 0. + allocate(CS%kd_conv(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kd_conv(:,:,:) = 0. + allocate(CS%kv_conv(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kv_conv(:,:,:) = 0. + + ! Register diagnostics + CS%diag => diag + CS%id_N2 = register_diag_field('ocean_model', 'N2_conv', diag%axesTi, Time, & + 'Square of Brunt-Vaisala frequency used by MOM_CVMix_conv module', '1/s2') + CS%id_kd_conv = register_diag_field('ocean_model', 'kd_conv', diag%axesTi, Time, & + 'Additional diffusivity added by MOM_CVMix_conv module', 'm2/s') + CS%id_kv_conv = register_diag_field('ocean_model', 'kv_conv', diag%axesTi, Time, & + 'Additional viscosity added by MOM_CVMix_conv module', 'm2/s') + + call CVMix_init_conv(convect_diff=CS%kd_conv_const, & + convect_visc=CS%kv_conv_const, & + lBruntVaisala=.true., & + BVsqr_convect=CS%bv_sqr_conv) + +end function CVMix_conv_init + +!> Subroutine for calculating enhanced diffusivity/viscosity +!! due to convection via CVMix +subroutine calculate_CVMix_conv(h, tv, G, GV, CS, hbl) + + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. + type(CVMix_conv_cs), pointer :: CS !< The control structure returned + !! by a previous call to CVMix_conv_init. + real, dimension(:,:), optional, pointer :: hbl!< Depth of ocean boundary layer (m) + + ! local variables + real, dimension(SZK_(G)) :: rho_lwr !< Adiabatic Water Density, this is a dummy + !! variable since here convection is always + !! computed based on Brunt Vaisala. + real, dimension(SZK_(G)) :: rho_1d !< water density in a column, this is also + !! a dummy variable, same reason as above. + real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) + real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) + integer :: kOBL !< level of OBL extent + real :: pref, g_o_rho0, rhok, rhokm1, dz, dh, hcorr + integer :: i, j, k + + g_o_rho0 = GV%g_Earth / GV%Rho0 + + ! initialize dummy variables + rho_lwr(:) = 0.0; rho_1d(:) = 0.0 + + if (.not. associated(hbl)) then + allocate(hbl(SZI_(G), SZJ_(G))); + hbl(:,:) = 0.0 + endif + + do j = G%jsc, G%jec + do i = G%isc, G%iec + + ! set N2 to zero at the top- and bottom-most interfaces + CS%N2(i,j,1) = 0. + CS%N2(i,j,G%ke+1) = 0. + + ! skip calling at land points + !if (G%mask2dT(i,j) == 0.) cycle + + pRef = 0. + ! Compute Brunt-Vaisala frequency (static stability) on interfaces + do k=2,G%ke + + ! pRef is pressure at interface between k and km1. + pRef = pRef + GV%H_to_Pa * h(i,j,k) + call calculate_density(tv%t(i,j,k), tv%s(i,j,k), pref, rhok, tv%eqn_of_state) + call calculate_density(tv%t(i,j,k-1), tv%s(i,j,k-1), pref, rhokm1, tv%eqn_of_state) + + dz = ((0.5*(h(i,j,k-1) + h(i,j,k))+GV%H_subroundoff)*GV%H_to_m) + CS%N2(i,j,k) = g_o_rho0 * (rhok - rhokm1) / dz ! Can be negative + + enddo + + iFaceHeight(1) = 0.0 ! BBL is all relative to the surface + hcorr = 0.0 + ! compute heights at cell center and interfaces + do k=1,G%ke + dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) + hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 + dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * dh + iFaceHeight(k+1) = iFaceHeight(k) - dh + enddo + + kOBL = CVMix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) + + call CVMix_coeffs_conv(Mdiff_out=CS%kv_conv(i,j,:), & + Tdiff_out=CS%kd_conv(i,j,:), & + Nsqr=CS%N2(i,j,:), & + dens=rho_1d(:), & + dens_lwr=rho_lwr(:), & + nlev=G%ke, & + max_nlev=G%ke, & + OBL_ind=kOBL) + + ! Do not apply mixing due to convection within the boundary layer + do k=1,kOBL + CS%kv_conv(i,j,k) = 0.0 + CS%kd_conv(i,j,k) = 0.0 + enddo + + enddo + enddo + + if (CS%debug) then + call hchksum(CS%N2, "MOM_CVMix_conv: N2",G%HI,haloshift=0) + call hchksum(CS%kd_conv, "MOM_CVMix_conv: kd_conv",G%HI,haloshift=0) + call hchksum(CS%kv_conv, "MOM_CVMix_conv: kv_conv",G%HI,haloshift=0) + endif + + ! send diagnostics to post_data + if (CS%id_N2 > 0) call post_data(CS%id_N2, CS%N2, CS%diag) + if (CS%id_kd_conv > 0) call post_data(CS%id_kd_conv, CS%kd_conv, CS%diag) + if (CS%id_kv_conv > 0) call post_data(CS%id_kv_conv, CS%kv_conv, CS%diag) + +end subroutine calculate_CVMix_conv + +!> Reads the parameter "USE_CVMix_CONVECTION" and returns state. +!! This function allows other modules to know whether this parameterization will +!! be used without needing to duplicate the log entry. +logical function CVMix_conv_is_used(param_file) + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + call get_param(param_file, mdl, "USE_CVMix_CONVECTION", CVMix_conv_is_used, & + default=.false., do_not_log = .true.) + +end function CVMix_conv_is_used + +!> Clear pointers and dealocate memory +subroutine CVMix_conv_end(CS) + type(CVMix_conv_cs), pointer :: CS ! Control structure + + deallocate(CS%N2) + deallocate(CS%kd_conv) + deallocate(CS%kv_conv) + deallocate(CS) + +end subroutine CVMix_conv_end + + +end module MOM_CVMix_conv diff --git a/src/parameterizations/vertical/MOM_cvmix_shear.F90 b/src/parameterizations/vertical/MOM_CVMix_shear.F90 similarity index 66% rename from src/parameterizations/vertical/MOM_cvmix_shear.F90 rename to src/parameterizations/vertical/MOM_CVMix_shear.F90 index 460dde7c47..f99a0d4dcb 100644 --- a/src/parameterizations/vertical/MOM_cvmix_shear.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_shear.F90 @@ -1,10 +1,10 @@ !> Interface to CVMix interior shear schemes -module MOM_cvmix_shear +module MOM_CVMix_shear ! This file is part of MOM6. See LICENSE.md for the license. !--------------------------------------------------- -! module MOM_cvmix_shear +! module MOM_CVMix_shear ! Author: Brandon Reichl ! Date: Aug 31, 2016 ! Purpose: Interface to CVMix interior shear schemes @@ -19,43 +19,50 @@ module MOM_cvmix_shear use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, EOS_type -use cvmix_shear, only : cvmix_init_shear, cvmix_coeffs_shear +use CVMix_shear, only : CVMix_init_shear, CVMix_coeffs_shear use MOM_kappa_shear, only : kappa_shear_is_used implicit none ; private #include -public Calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_is_used +public calculate_CVMix_shear, CVMix_shear_init, CVMix_shear_is_used, CVMix_shear_end !> Control structure including parameters for CVMix interior shear schemes. -type, public :: CVMix_shear_CS +type, public :: CVMix_shear_cs logical :: use_LMD94, use_PP81 !< Flags for various schemes real :: Ri_zero !< LMD94 critical Richardson number real :: Nu_zero !< LMD94 maximum interior diffusivity real :: KPP_exp !< real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) real, allocatable, dimension(:,:,:) :: S2 !< Squared shear frequency (1/s2) + real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number +! real, allocatable, dimension(:,:,:) :: kv !< vertical viscosity at interface (m2/s) +! real, allocatable, dimension(:,:,:) :: kd !< vertical diffusivity at interface (m2/s) character(10) :: Mix_Scheme !< Mixing scheme name (string) -end type CVMix_shear_CS + ! Daignostic handles and pointers + type(diag_ctrl), pointer :: diag => NULL() + integer :: id_N2 = -1, id_S2 = -1, id_ri_grad = -1, id_kv = -1, id_kd = -1 + +end type CVMix_shear_cs character(len=40) :: mdl = "MOM_CVMix_shear" !< This module's name. contains -!> Subroutine for calculating (internal) diffusivity -subroutine Calculate_cvmix_shear(u_H, v_H, h, tv, KH, & - KM, G, GV, CS ) +!> Subroutine for calculating (internal) vertical diffusivities/viscosities +subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, & + kv, G, GV, CS ) type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: u_H !< Initial zonal velocity on T points, in m s-1. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: v_H !< Initial meridional velocity on T points, in m s-1. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: KH !< The vertical viscosity at each interface + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd !< The vertical diffusivity at each interface !! (not layer!) in m2 s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: KM !< The vertical viscosity at each interface + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv !< The vertical viscosity at each interface !! (not layer!) in m2 s-1. - type(CVMix_shear_CS), pointer :: CS !< The control structure returned by a previous call to + type(CVMix_shear_cs), pointer :: CS !< The control structure returned by a previous call to !! CVMix_shear_init. ! Local variables integer :: i, j, k, kk, km1 @@ -109,31 +116,44 @@ subroutine Calculate_cvmix_shear(u_H, v_H, h, tv, KH, & N2 = DRHO/DZ S2 = (DU*DU+DV*DV)/(DZ*DZ) Ri_Grad(k) = max(0.,N2)/max(S2,1.e-16) + + ! fill 3d arrays, if user asks for diagsnostics + if (CS%id_N2 > 0) CS%N2(i,j,k) = N2 + if (CS%id_S2 > 0) CS%S2(i,j,k) = S2 + if (CS%id_ri_grad > 0) CS%ri_grad(i,j,k) = Ri_Grad(k) + enddo ! Call to CVMix wrapper for computing interior mixing coefficients. - call cvmix_coeffs_shear(Mdiff_out=KM(i,j,:), & - Tdiff_out=KH(i,j,:), & + call CVMix_coeffs_shear(Mdiff_out=kv(i,j,:), & + Tdiff_out=kd(i,j,:), & RICH=Ri_Grad, & nlev=G%ke, & max_nlev=G%ke) enddo enddo -end subroutine Calculate_cvmix_shear + ! write diagnostics + if (CS%id_kd > 0) call post_data(CS%id_kd,kd, CS%diag) + if (CS%id_kv > 0) call post_data(CS%id_kv,kv, CS%diag) + if (CS%id_N2 > 0) call post_data(CS%id_N2,CS%N2, CS%diag) + if (CS%id_S2 > 0) call post_data(CS%id_S2,CS%S2, CS%diag) + if (CS%id_ri_grad > 0) call post_data(CS%id_ri_grad,CS%ri_grad, CS%diag) +end subroutine calculate_CVMix_shear -!> Initialized the cvmix internal shear mixing routine. + +!> Initialized the CVMix internal shear mixing routine. !! \note *This is where we test to make sure multiple internal shear !! mixing routines (including JHL) are not enabled at the same time. -!! (returns) cvmix_shear_init - True if module is to be used, False otherwise -logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) +!! (returns) CVMix_shear_init - True if module is to be used, False otherwise +logical function CVMix_shear_init(Time, G, GV, param_file, diag, CS) type(time_type), intent(in) :: Time !< The current time. type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. - type(CVMix_shear_CS), pointer :: CS !< This module's control structure. + type(CVMix_shear_cs), pointer :: CS !< This module's control structure. ! Local variables integer :: NumberTrue=0 logical :: use_JHL @@ -141,7 +161,7 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) #include "version_variable.h" if (associated(CS)) then - call MOM_error(WARNING, "cvmix_shear_init called with an associated "// & + call MOM_error(WARNING, "CVMix_shear_init called with an associated "// & "control structure.") return endif @@ -169,14 +189,14 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) ! After testing for interior schemes, make sure only 0 or 1 are enabled. ! Otherwise, warn user and kill job. if ((NumberTrue).gt.1) then - call MOM_error(FATAL, 'MOM_cvmix_shear_init: '// & + call MOM_error(FATAL, 'MOM_CVMix_shear_init: '// & 'Multiple shear driven internal mixing schemes selected,'//& ' please disable all but one scheme to proceed.') endif - cvmix_shear_init=(CS%use_PP81.or.CS%use_LMD94) + CVMix_shear_init=(CS%use_PP81.or.CS%use_LMD94) ! Forego remainder of initialization if not using this scheme - if (.not. cvmix_shear_init) return + if (.not. CVMix_shear_init) return call get_param(param_file, mdl, "NU_ZERO", CS%Nu_Zero, & "Leading coefficient in KPP shear mixing.", & units="nondim", default=5.e-3) @@ -189,20 +209,40 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) "Exponent of unitless factor of diffusivities,"// & " for KPP internal shear mixing scheme." & ,units="nondim", default=3.0) - call cvmix_init_shear(mix_scheme=CS%mix_scheme, & + call CVMix_init_shear(mix_scheme=CS%mix_scheme, & KPP_nu_zero=CS%Nu_Zero, & KPP_Ri_zero=CS%Ri_zero, & KPP_exp=CS%KPP_exp) - ! Allocation and initialization - allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%N2(:,:,:) = 0. - allocate( CS%S2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%S2(:,:,:) = 0. -end function cvmix_shear_init + ! Register diagnostics; allocation and initialization + CS%diag => diag + + CS%id_N2 = register_diag_field('ocean_model', 'N2_shear', diag%axesTi, Time, & + 'Square of Brunt-Vaisala frequency used by MOM_CVMix_shear module', '1/s2') + if (CS%id_N2 > 0) & + allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%N2(:,:,:) = 0. + + CS%id_S2 = register_diag_field('ocean_model', 'S2_shear', diag%axesTi, Time, & + 'Square of vertical shear used by MOM_CVMix_shear module','1/s2') + if (CS%id_S2 > 0) & + allocate( CS%S2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%S2(:,:,:) = 0. + + CS%id_ri_grad = register_diag_field('ocean_model', 'ri_grad_shear', diag%axesTi, Time, & + 'Gradient Richarson number used by MOM_CVMix_shear module','nondim') + if (CS%id_ri_grad > 0) & !Initialize w/ large Richardson value + allocate( CS%ri_grad( SZI_(G), SZJ_(G), SZK_(G)+1 ));CS%ri_grad(:,:,:) = 1.e8 + + CS%id_kd = register_diag_field('ocean_model', 'kd_shear_CVMix', diag%axesTi, Time, & + 'Vertical diffusivity added by MOM_CVMix_shear module', 'm2/s') + CS%id_kv = register_diag_field('ocean_model', 'kv_shear_CVMix', diag%axesTi, Time, & + 'Vertical viscosity added by MOM_CVMix_shear module', 'm2/s') + +end function CVMix_shear_init !> Reads the parameters "LMD94" and "PP81" and returns state. !! This function allows other modules to know whether this parameterization will !! be used without needing to duplicate the log entry. -logical function cvmix_shear_is_used(param_file) +logical function CVMix_shear_is_used(param_file) type(param_file_type), intent(in) :: param_file !< Run-time parameter files handle. ! Local variables logical :: LMD94, PP81 @@ -210,7 +250,18 @@ logical function cvmix_shear_is_used(param_file) default=.false., do_not_log = .true.) call get_param(param_file, mdl, "Use_PP81", PP81, & default=.false., do_not_log = .true.) - cvmix_shear_is_used = (LMD94 .or. PP81) -end function cvmix_shear_is_used + CVMix_shear_is_used = (LMD94 .or. PP81) +end function CVMix_shear_is_used + +!> Clear pointers and dealocate memory +subroutine CVMix_shear_end(CS) + type(CVMix_shear_cs), pointer :: CS ! Control structure + + if (CS%id_N2 > 0) deallocate(CS%N2) + if (CS%id_S2 > 0) deallocate(CS%S2) + if (CS%id_ri_grad > 0) deallocate(CS%ri_grad) + deallocate(CS) + +end subroutine CVMix_shear_end -end module MOM_cvmix_shear +end module MOM_CVMix_shear diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 26dc241dd6..dd6c07813f 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -15,14 +15,14 @@ module MOM_KPP use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number -use CVmix_kpp, only : CVmix_init_kpp, CVmix_put_kpp, CVmix_get_kpp_real -use CVmix_kpp, only : CVmix_coeffs_kpp -use CVmix_kpp, only : CVmix_kpp_compute_OBL_depth -use CVmix_kpp, only : CVmix_kpp_compute_turbulent_scales -use CVmix_kpp, only : CVmix_kpp_compute_bulk_Richardson -use CVmix_kpp, only : CVmix_kpp_compute_unresolved_shear -use CVmix_kpp, only : CVmix_kpp_params_type -use CVmix_kpp, only : CVmix_kpp_compute_kOBL_depth +use CVMix_kpp, only : CVMix_init_kpp, CVMix_put_kpp, CVMix_get_kpp_real +use CVMix_kpp, only : CVMix_coeffs_kpp +use CVMix_kpp, only : CVMix_kpp_compute_OBL_depth +use CVMix_kpp, only : CVMix_kpp_compute_turbulent_scales +use CVMix_kpp, only : CVMix_kpp_compute_bulk_Richardson +use CVMix_kpp, only : CVMix_kpp_compute_unresolved_shear +use CVMix_kpp, only : CVMix_kpp_params_type +use CVMix_kpp, only : CVMix_kpp_compute_kOBL_depth implicit none ; private @@ -33,9 +33,10 @@ module MOM_KPP public :: KPP_end public :: KPP_NonLocalTransport_temp public :: KPP_NonLocalTransport_saln +public :: KPP_get_BLD ! Enumerated constants -integer, private, parameter :: NLT_SHAPE_CVMIX = 0 !< Use the CVmix profile +integer, private, parameter :: NLT_SHAPE_CVMix = 0 !< Use the CVMix profile integer, private, parameter :: NLT_SHAPE_LINEAR = 1 !< Linear, \f$ G(\sigma) = 1-\sigma \f$ integer, private, parameter :: NLT_SHAPE_PARABOLIC = 2 !< Parabolic, \f$ G(\sigma) = (1-\sigma)^2 \f$ integer, private, parameter :: NLT_SHAPE_CUBIC = 3 !< Cubic, \f$ G(\sigma) = 1 + (2\sigma-3) \sigma^2\f$ @@ -102,8 +103,8 @@ module MOM_KPP logical :: STOKES_MIXING !< Flag if model is mixing down Stokes gradient !! This is relavent for which current to use in RiB - !> CVmix parameters - type(CVmix_kpp_params_type), pointer :: KPP_params => NULL() + !> CVMix parameters + type(CVMix_kpp_params_type), pointer :: KPP_params => NULL() ! Diagnostic handles and pointers type(diag_ctrl), pointer :: diag => NULL() @@ -158,7 +159,7 @@ module MOM_KPP contains -!> Initialize the CVmix KPP module and set up diagnostics +!> Initialize the CVMix KPP module and set up diagnostics !! Returns True if KPP is to be used, False otherwise. logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) @@ -182,10 +183,10 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) allocate(CS) ! Read parameters - call log_version(paramFile, mdl, version, 'This is the MOM wrapper to CVmix:KPP\n' // & - 'See http://code.google.com/p/cvmix/') + call log_version(paramFile, mdl, version, 'This is the MOM wrapper to CVMix:KPP\n' // & + 'See http://cvmix.github.io/') call get_param(paramFile, mdl, "USE_KPP", KPP_init, & - "If true, turns on the [CVmix] KPP scheme of Large et al., 1994,\n"// & + "If true, turns on the [CVMix] KPP scheme of Large et al., 1994,\n"// & "to calculate diffusivities and non-local transport in the OBL.", & default=.false.) ! Forego remainder of initialization if not using this scheme @@ -274,14 +275,14 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) call get_param(paramFile, mdl, 'NLT_SHAPE', string, & 'MOM6 method to set nonlocal transport profile.\n'// & 'Over-rides the result from CVMix. Allowed values are: \n'// & - '\t CVMIX - Uses the profiles from CVmix specified by MATCH_TECHNIQUE\n'//& + '\t CVMix - Uses the profiles from CVMix specified by MATCH_TECHNIQUE\n'//& '\t LINEAR - A linear profile, 1-sigma\n'// & '\t PARABOLIC - A parablic profile, (1-sigma)^2\n'// & '\t CUBIC - A cubic profile, (1-sigma)^2(1+2*sigma)\n'// & '\t CUBIC_LMD - The original KPP profile', & - default='CVMIX') + default='CVMix') select case ( trim(string) ) - case ("CVMIX") ; CS%NLT_shape = NLT_SHAPE_CVMIX + case ("CVMix") ; CS%NLT_shape = NLT_SHAPE_CVMix case ("LINEAR") ; CS%NLT_shape = NLT_SHAPE_LINEAR case ("PARABOLIC") ; CS%NLT_shape = NLT_SHAPE_PARABOLIC case ("CUBIC") ; CS%NLT_shape = NLT_SHAPE_CUBIC @@ -299,7 +300,7 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) default='SimpleShapes') if (CS%MatchTechnique.eq.'ParabolicNonLocal') then ! This forces Cs2 (Cs in non-local computation) to equal 1 for parabolic non-local option. - ! May be used during CVmix initialization. + ! May be used during CVMix initialization. Cs_is_one=.true. endif call get_param(paramFile, mdl, 'KPP_ZERO_DIFFUSIVITY', CS%KPPzeroDiffusivity, & @@ -323,7 +324,7 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) case default ; call MOM_error(FATAL,"KPP_init: "// & "Unrecognized KPP_SHORTWAVE_METHOD option"//trim(string)) end select - call get_param(paramFile, mdl, 'CVMIX_ZERO_H_WORK_AROUND', CS%min_thickness, & + call get_param(paramFile, mdl, 'CVMix_ZERO_H_WORK_AROUND', CS%min_thickness, & 'A minimum thickness used to avoid division by small numbers in the vicinity\n'// & 'of vanished layers. This is independent of MIN_THICKNESS used in other parts of MOM.', & units='m', default=0.) @@ -400,7 +401,7 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) call closeParameterBlock(paramFile) call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) - call CVmix_init_kpp( Ri_crit=CS%Ri_crit, & + call CVMix_init_kpp( Ri_crit=CS%Ri_crit, & minOBLdepth=CS%minOBLdepth, & minVtsqr=CS%minVtsqr, & vonKarman=CS%vonKarman, & @@ -412,73 +413,73 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) MatchTechnique=CS%MatchTechnique, & lenhanced_diff=CS%enhance_diffusion,& lnonzero_surf_nonlocal=Cs_is_one ,& - CVmix_kpp_params_user=CS%KPP_params ) + CVMix_kpp_params_user=CS%KPP_params ) ! Register diagnostics CS%diag => diag CS%id_OBLdepth = register_diag_field('ocean_model', 'KPP_OBLdepth', diag%axesT1, Time, & - 'Thickness of the surface Ocean Boundary Layer calculated by [CVmix] KPP', 'meter', & + 'Thickness of the surface Ocean Boundary Layer calculated by [CVMix] KPP', 'meter', & cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', & cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme') ! CMOR names are placeholders; must be modified by time period ! for CMOR compliance. Diag manager will be used for omlmax and ! omldamax. CS%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, Time, & - 'Bulk difference in density used in Bulk Richardson number, as used by [CVmix] KPP', 'kg/m3') + 'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', 'kg/m3') CS%id_BulkUz2 = register_diag_field('ocean_model', 'KPP_BulkUz2', diag%axesTL, Time, & - 'Square of bulk difference in resolved velocity used in Bulk Richardson number via [CVmix] KPP', 'm2/s2') + 'Square of bulk difference in resolved velocity used in Bulk Richardson number via [CVMix] KPP', 'm2/s2') CS%id_BulkRi = register_diag_field('ocean_model', 'KPP_BulkRi', diag%axesTL, Time, & - 'Bulk Richardson number used to find the OBL depth used by [CVmix] KPP', 'nondim') + 'Bulk Richardson number used to find the OBL depth used by [CVMix] KPP', 'nondim') CS%id_Sigma = register_diag_field('ocean_model', 'KPP_sigma', diag%axesTi, Time, & - 'Sigma coordinate used by [CVmix] KPP', 'nondim') + 'Sigma coordinate used by [CVMix] KPP', 'nondim') CS%id_Ws = register_diag_field('ocean_model', 'KPP_Ws', diag%axesTL, Time, & - 'Turbulent vertical velocity scale for scalars used by [CVmix] KPP', 'm/s') + 'Turbulent vertical velocity scale for scalars used by [CVMix] KPP', 'm/s') CS%id_N = register_diag_field('ocean_model', 'KPP_N', diag%axesTi, Time, & - '(Adjusted) Brunt-Vaisala frequency used by [CVmix] KPP', '1/s') + '(Adjusted) Brunt-Vaisala frequency used by [CVMix] KPP', '1/s') CS%id_N2 = register_diag_field('ocean_model', 'KPP_N2', diag%axesTi, Time, & - 'Square of Brunt-Vaisala frequency used by [CVmix] KPP', '1/s2') + 'Square of Brunt-Vaisala frequency used by [CVMix] KPP', '1/s2') CS%id_Vt2 = register_diag_field('ocean_model', 'KPP_Vt2', diag%axesTL, Time, & - 'Unresolved shear turbulence used by [CVmix] KPP', 'm2/s2') + 'Unresolved shear turbulence used by [CVMix] KPP', 'm2/s2') CS%id_uStar = register_diag_field('ocean_model', 'KPP_uStar', diag%axesT1, Time, & - 'Friction velocity, u*, as used by [CVmix] KPP', 'm/s') + 'Friction velocity, u*, as used by [CVMix] KPP', 'm/s') CS%id_buoyFlux = register_diag_field('ocean_model', 'KPP_buoyFlux', diag%axesTi, Time, & - 'Surface (and penetrating) buoyancy flux, as used by [CVmix] KPP', 'm2/s3') + 'Surface (and penetrating) buoyancy flux, as used by [CVMix] KPP', 'm2/s3') CS%id_QminusSW = register_diag_field('ocean_model', 'KPP_QminusSW', diag%axesT1, Time, & - 'Net temperature flux ignoring short-wave, as used by [CVmix] KPP', 'K m/s') + 'Net temperature flux ignoring short-wave, as used by [CVMix] KPP', 'K m/s') CS%id_netS = register_diag_field('ocean_model', 'KPP_netSalt', diag%axesT1, Time, & - 'Effective net surface salt flux, as used by [CVmix] KPP', 'ppt m/s') + 'Effective net surface salt flux, as used by [CVMix] KPP', 'ppt m/s') CS%id_Kt_KPP = register_diag_field('ocean_model', 'KPP_Kheat', diag%axesTi, Time, & - 'Heat diffusivity due to KPP, as calculated by [CVmix] KPP', 'm2/s') + 'Heat diffusivity due to KPP, as calculated by [CVMix] KPP', 'm2/s') CS%id_Kd_in = register_diag_field('ocean_model', 'KPP_Kd_in', diag%axesTi, Time, & 'Diffusivity passed to KPP', 'm2/s') CS%id_Ks_KPP = register_diag_field('ocean_model', 'KPP_Ksalt', diag%axesTi, Time, & - 'Salt diffusivity due to KPP, as calculated by [CVmix] KPP', 'm2/s') + 'Salt diffusivity due to KPP, as calculated by [CVMix] KPP', 'm2/s') CS%id_Kv_KPP = register_diag_field('ocean_model', 'KPP_Kv', diag%axesTi, Time, & - 'Vertical viscosity due to KPP, as calculated by [CVmix] KPP', 'm2/s') + 'Vertical viscosity due to KPP, as calculated by [CVMix] KPP', 'm2/s') CS%id_NLTt = register_diag_field('ocean_model', 'KPP_NLtransport_heat', diag%axesTi, Time, & - 'Non-local transport (Cs*G(sigma)) for heat, as calculated by [CVmix] KPP', 'nondim') + 'Non-local transport (Cs*G(sigma)) for heat, as calculated by [CVMix] KPP', 'nondim') CS%id_NLTs = register_diag_field('ocean_model', 'KPP_NLtransport_salt', diag%axesTi, Time, & - 'Non-local tranpsort (Cs*G(sigma)) for scalars, as calculated by [CVmix] KPP', 'nondim') + 'Non-local tranpsort (Cs*G(sigma)) for scalars, as calculated by [CVMix] KPP', 'nondim') CS%id_NLT_dTdt = register_diag_field('ocean_model', 'KPP_NLT_dTdt', diag%axesTL, Time, & - 'Temperature tendency due to non-local transport of heat, as calculated by [CVmix] KPP', 'K/s') + 'Temperature tendency due to non-local transport of heat, as calculated by [CVMix] KPP', 'K/s') CS%id_NLT_dSdt = register_diag_field('ocean_model', 'KPP_NLT_dSdt', diag%axesTL, Time, & - 'Salinity tendency due to non-local transport of salt, as calculated by [CVmix] KPP', 'ppt/s') + 'Salinity tendency due to non-local transport of salt, as calculated by [CVMix] KPP', 'ppt/s') CS%id_NLT_temp_budget = register_diag_field('ocean_model', 'KPP_NLT_temp_budget', diag%axesTL, Time, & - 'Heat content change due to non-local transport, as calculated by [CVmix] KPP', 'W/m^2') + 'Heat content change due to non-local transport, as calculated by [CVMix] KPP', 'W/m^2') CS%id_NLT_saln_budget = register_diag_field('ocean_model', 'KPP_NLT_saln_budget', diag%axesTL, Time, & - 'Salt content change due to non-local transport, as calculated by [CVmix] KPP', 'kg/(sec*m^2)') + 'Salt content change due to non-local transport, as calculated by [CVMix] KPP', 'kg/(sec*m^2)') CS%id_Tsurf = register_diag_field('ocean_model', 'KPP_Tsurf', diag%axesT1, Time, & - 'Temperature of surface layer (10% of OBL depth) as passed to [CVmix] KPP', 'C') + 'Temperature of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'C') CS%id_Ssurf = register_diag_field('ocean_model', 'KPP_Ssurf', diag%axesT1, Time, & - 'Salinity of surface layer (10% of OBL depth) as passed to [CVmix] KPP', 'ppt') + 'Salinity of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'ppt') CS%id_Usurf = register_diag_field('ocean_model', 'KPP_Usurf', diag%axesCu1, Time, & - 'i-component flow of surface layer (10% of OBL depth) as passed to [CVmix] KPP', 'm/s') + 'i-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'm/s') CS%id_Vsurf = register_diag_field('ocean_model', 'KPP_Vsurf', diag%axesCv1, Time, & - 'j-component flow of surface layer (10% of OBL depth) as passed to [CVmix] KPP', 'm/s') + 'j-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'm/s') CS%id_EnhK = register_diag_field('ocean_model', 'EnhK', diag%axesTI, Time, & - 'Langmuir number enhancement to K as used by [CVmix] KPP','nondim') + 'Langmuir number enhancement to K as used by [CVMix] KPP','nondim') CS%id_EnhVt2 = register_diag_field('ocean_model', 'EnhVt2', diag%axesTL, Time, & - 'Langmuir number enhancement to Vt2 as used by [CVmix] KPP','nondim') + 'Langmuir number enhancement to Vt2 as used by [CVMix] KPP','nondim') allocate( CS%OBLdepthprev( SZI_(G), SZJ_(G) ) );CS%OBLdepthprev(:,:)=0.0; if (CS%id_OBLdepth > 0) allocate( CS%OBLdepth( SZI_(G), SZJ_(G) ) ) @@ -788,16 +789,16 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & N_1d(G%ke+1 ) = 0.0 ! turbulent velocity scales w_s and w_m computed at the cell centers. - ! Note that if sigma > CS%surf_layer_ext, then CVmix_kpp_compute_turbulent_scales + ! Note that if sigma > CS%surf_layer_ext, then CVMix_kpp_compute_turbulent_scales ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass ! sigma=CS%surf_layer_ext for this calculation. - call CVmix_kpp_compute_turbulent_scales( & + call CVMix_kpp_compute_turbulent_scales( & CS%surf_layer_ext, & ! (in) Normalized surface layer depth; sigma = CS%surf_layer_ext -cellHeight, & ! (in) Assume here that OBL depth (m) = -cellHeight(k) surfBuoyFlux2, & ! (in) Buoyancy flux at surface (m2/s3) surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) w_s=Ws_1d, & ! (out) Turbulent velocity scale profile (m/s) - CVmix_kpp_params_user=CS%KPP_params ) + CVMix_kpp_params_user=CS%KPP_params ) !Compute CVMix VT2 Vt2_1d(:) = CVmix_kpp_compute_unresolved_shear( & @@ -848,7 +849,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & zt_cntr = cellHeight(1:G%ke), & ! Depth of cell center (m) delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) (1/s) delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference (m2/s2) - Vt_sqr_cntr=VT2_1d, & + Vt_sqr_cntr=Vt2_1d, & ws_cntr=Ws_1d, & ! Turbulent velocity scale profile (m/s) N_iface=N_1d) ! Buoyancy frequency (1/s) @@ -856,7 +857,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & surfBuoyFlux = buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit ! h to Monin-Obukov (default is false, ie. not used) - call CVmix_kpp_compute_OBL_depth( & + call CVMix_kpp_compute_OBL_depth( & BulkRi_1d, & ! (in) Bulk Richardson number iFaceHeight, & ! (in) Height of interfaces (m) OBLdepth_0d, & ! (out) OBL depth (m) @@ -865,7 +866,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) Coriolis=Coriolis, & ! (in) Coriolis parameter (1/s) - CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters + CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters ! A hack to avoid KPP reaching the bottom. It was needed during development ! because KPP was unable to handle vanishingly small layers near the bottom. @@ -878,7 +879,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & if(CS%fixedOBLdepth) OBLdepth_0d = CS%fixedOBLdepth_value OBLdepth_0d = max( OBLdepth_0d, -iFaceHeight(2) ) ! no shallower than top layer OBLdepth_0d = min( OBLdepth_0d, -iFaceHeight(G%ke+1) ) ! no deeper than bottom - kOBL = CVmix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) + kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) !************************************************************************* ! smg: remove code below @@ -983,7 +984,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & Kviscosity(:)=Kv(i,j,:) endif - call cvmix_coeffs_kpp(Kviscosity, & ! (inout) Total viscosity (m2/s) + call CVMix_coeffs_kpp(Kviscosity, & ! (inout) Total viscosity (m2/s) Kdiffusivity(:,1), & ! (inout) Total heat diffusivity (m2/s) Kdiffusivity(:,2), & ! (inout) Total salt diffusivity (m2/s) iFaceHeight, & ! (in) Height of interfaces (m) @@ -999,7 +1000,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) G%ke, & ! (in) Number of levels to compute coeffs for G%ke, & ! (in) Number of levels in array shape - CVmix_kpp_params_user=CS%KPP_params ) + CVMix_kpp_params_user=CS%KPP_params ) IF (CS%LT_K_ENHANCEMENT) then if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then @@ -1093,7 +1094,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) w_s=Ws_1d, & ! (out) Turbulent velocity scale profile (m/s) - CVmix_kpp_params_user=CS%KPP_params) ! KPP parameters + CVmix_kpp_params_user=CS%KPP_params) ! KPP parameters CS%Ws(i,j,:) = Ws_1d(:) endif @@ -1190,7 +1191,18 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & end subroutine KPP_calculate - +!> Copies KPP surface boundary layer depth into BLD +subroutine KPP_get_BLD(CS, BLD, G) + type(KPP_CS), pointer :: CS !< Control structure for + !! this module + type(ocean_grid_type), intent(in) :: G !< Grid structure + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: BLD!< bnd. layer depth (m) + ! Local variables + integer :: i,j + do j = G%jsc, G%jec ; do i = G%isc, G%iec + BLD(i,j) = CS%OBLdepth(i,j) + enddo ; enddo +end subroutine KPP_get_BLD !> Apply KPP non-local transport of surface fluxes for temperature. subroutine KPP_NonLocalTransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, & @@ -1315,7 +1327,10 @@ end subroutine KPP_NonLocalTransport_saln subroutine KPP_end(CS) type(KPP_CS), pointer :: CS !< Control structure + if (.not.associated(CS)) return + deallocate(CS) + end subroutine KPP_end !> \namespace mom_kpp @@ -1323,11 +1338,11 @@ end subroutine KPP_end !! \section section_KPP The K-Profile Parameterization !! !! The K-Profile Parameterization (KPP) of Large et al., 1994, (http://dx.doi.org/10.1029/94RG01872) is -!! implemented via the Community Vertical Mixing package, [CVmix](https://code.google.com/p/cvmix), +!! implemented via the Community Vertical Mixing package, [CVMix](http://cvmix.github.io/), !! which is called directly by this module. !! !! The formulation and implementation of KPP is described in great detail in the -!! [CVMix manual](https://cvmix.googlecode.com/svn/trunk/manual/cvmix.pdf) (written by our own Stephen Griffies). +!! [CVMix manual](https://github.com/CVMix/CVMix-description/raw/master/cvmix.pdf) (written by our own Stephen Griffies). !! !! \subsection section_KPP_nutshell KPP in a nutshell !! @@ -1348,7 +1363,7 @@ end subroutine KPP_end !! In our implementation of KPP, we allow the shape functions used for \f$ K \f$ and for the non-local transport !! to be chosen independently. !! -!! [google_thread_NLT]: https://groups.google.com/forum/#!msg/cvmix-dev/i6rF-eHOtKI/Ti8BeyksrhAJ "Extreme values of non-local transport" +!! [google_thread_NLT]: https://groups.google.com/forum/#!msg/CVMix-dev/i6rF-eHOtKI/Ti8BeyksrhAJ "Extreme values of non-local transport" !! !! The particular shape function most widely used in the atmospheric community is !! \f[ G(\sigma) = \sigma (1-\sigma)^2 \f] diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 new file mode 100644 index 0000000000..e9441d36e5 --- /dev/null +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -0,0 +1,449 @@ +!> Interface to background mixing schemes, including the Bryan and Lewis (1979) +!! which is applied via CVMix. + +module MOM_bkgnd_mixing + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field +use MOM_diag_mediator, only : post_data +use MOM_EOS, only : calculate_density, calculate_density_derivs +use MOM_variables, only : thermo_var_ptrs +use MOM_forcing_type, only : forcing +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_error_handler, only : is_root_pe +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_debugging, only : hchksum +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_file_parser, only : get_param, log_version, param_file_type +use CVMix_background, only : CVMix_init_bkgnd, CVMix_coeffs_bkgnd +use MOM_variables, only : vertvisc_type +use MOM_intrinsic_functions, only : invcosh + +implicit none ; private + +#include + +public bkgnd_mixing_init +public bkgnd_mixing_end +public calculate_bkgnd_mixing +public sfc_bkgnd_mixing + +!> Control structure including parameters for this module. +type, public :: bkgnd_mixing_cs + + ! Parameters + real :: Bryan_Lewis_c1 !< The vertical diffusivity values for Bryan-Lewis profile + !! at |z|=D (m2/s) + real :: Bryan_Lewis_c2 !< The amplitude of variation in diffusivity for the + !! Bryan-Lewis diffusivity profile (m2/s) + real :: Bryan_Lewis_c3 !< The inverse length scale for transition region in the + !! Bryan-Lewis diffusivity profile (1/m) + real :: Bryan_Lewis_c4 !< The depth where diffusivity is Bryan_Lewis_bl1 in the + !! Bryan-Lewis profile (m) + real :: Kd_min !< minimum diapycnal diffusivity (m2/s) + real :: Kd !< interior diapycnal diffusivity (m2/s) + real :: N0_2Omega !< ratio of the typical Buoyancy frequency to + !! twice the Earth's rotation period, used with the + !! Henyey scaling from the mixing + real :: prandtl_bkgnd !< Turbulent Prandtl number used to convert + !! vertical background diffusivity into viscosity + real :: Kd_tanh_lat_scale !< A nondimensional scaling for the range of + !! diffusivities with Kd_tanh_lat_fn. Valid values + !! are in the range of -2 to 2; 0.4 reproduces CM2M. + real :: Kdml !< mixed layer diapycnal diffusivity (m2/s) + !! when bulkmixedlayer==.false. + real :: Hmix !< mixed layer thickness (meter) when + !! bulkmixedlayer==.false. + logical :: Kd_tanh_lat_fn !< If true, use the tanh dependence of Kd_sfc on + !! latitude, like GFDL CM2.1/CM2M. There is no + !! physical justification for this form, and it can + !! not be used with Henyey_IGW_background. + logical :: Bryan_Lewis_diffusivity!< If true, background vertical diffusivity + !! uses Bryan-Lewis (1979) like tanh profile. + logical :: Henyey_IGW_background !< If true, use a simplified variant of the + !! Henyey et al, JGR (1986) latitudinal scaling for the background diapycnal diffusivity, + !! which gives a marked decrease in the diffusivity near the equator. The simplification + !! here is to assume that the in-situ stratification is the same as the reference stratificaiton. + logical :: Henyey_IGW_background_new !< same as Henyey_IGW_background + !! but incorporate the effect of stratification on TKE dissipation, + !! e = f/f_0 * acosh(N/f) / acosh(N_0/f_0) * e_0 + !! where e is the TKE dissipation, and N_0 and f_0 + !! are the reference buoyancy frequency and inertial frequencies respectively. + !! e_0 is the reference dissipation at (N_0,f_0). In the previous version, N=N_0. + !! Additionally, the squared inverse relationship between diapycnal diffusivities + !! and stratification is included: + !! + !! kd = e/N^2 + !! + !! where kd is the diapycnal diffusivity. This approach assumes that work done + !! against gravity is uniformly distributed throughout the column. Whereas, kd=kd_0*e, + !! as in the original version, concentrates buoyancy work in regions of strong stratification. + logical :: bulkmixedlayer !< If true, a refined bulk mixed layer scheme is used + logical :: debug !< If true, turn on debugging in this module + ! Daignostic handles and pointers + type(diag_ctrl), pointer :: diag => NULL() + integer :: id_kd_bkgnd = -1, id_kv_bkgnd = -1 + + real, allocatable, dimension(:,:) :: Kd_sfc !< surface value of the diffusivity (m2/s) + ! Diagnostics arrays + real, allocatable, dimension(:,:,:) :: kd_bkgnd !< Background diffusivity (m2/s) + real, allocatable, dimension(:,:,:) :: kv_bkgnd !< Background viscosity (m2/s) + +end type bkgnd_mixing_cs + +character(len=40) :: mdl = "MOM_bkgnd_mixing" !< This module's name. + +contains + +!> Initialize the background mixing routine. +subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) + + type(time_type), intent(in) :: Time !< The current time. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. + type(bkgnd_mixing_cs), pointer :: CS !< This module's control structure. + + ! Local variables + +! This include declares and sets the variable "version". +#include "version_variable.h" + + if (associated(CS)) then + call MOM_error(WARNING, "bkgnd_mixing_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + ! Read parameters + call log_version(param_file, mdl, version, & + "Adding static vertical background mixing coefficients") + + call get_param(param_file, mdl, "KD", CS%Kd, & + "The background diapycnal diffusivity of density in the \n"//& + "interior. Zero or the molecular value, ~1e-7 m2 s-1, \n"//& + "may be used.", units="m2 s-1", fail_if_missing=.true.) + + call get_param(param_file, mdl, "KD_MIN", CS%Kd_min, & + "The minimum diapycnal diffusivity.", & + units="m2 s-1", default=0.01*CS%Kd) + + ! The following is needed to set one of the choices of vertical background mixing + + ! BULKMIXEDLAYER is not always defined (e.g., CM2G63L), so the following by pass + ! the need to include BULKMIXEDLAYER in MOM_input + CS%bulkmixedlayer = (GV%nkml > 0) + if (CS%bulkmixedlayer) then + ! Check that Kdml is not set when using bulk mixed layer + call get_param(param_file, mdl, "KDML", CS%Kdml, default=-1.) + if (CS%Kdml>0.) call MOM_error(FATAL, & + "bkgnd_mixing_init: KDML cannot be set when using"// & + "bulk mixed layer.") + CS%Kdml = CS%Kd ! This is not used with a bulk mixed layer, but also + ! cannot be a NaN. + else + call get_param(param_file, mdl, "KDML", CS%Kdml, & + "If BULKMIXEDLAYER is false, KDML is the elevated \n"//& + "diapycnal diffusivity in the topmost HMIX of fluid. \n"//& + "KDML is only used if BULKMIXEDLAYER is false.", & + units="m2 s-1", default=CS%Kd) + call get_param(param_file, mdl, "HMIX_FIXED", CS%Hmix, & + "The prescribed depth over which the near-surface \n"//& + "viscosity and diffusivity are elevated when the bulk \n"//& + "mixed layer is not used.", units="m", fail_if_missing=.true.) + endif + + call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) + + call get_param(param_file, mdl, "PRANDTL_BKGND", CS%prandtl_bkgnd, & + "Turbulent Prandtl number used to convert vertical \n"//& + "background diffusivities into viscosities.", & + units="nondim", default=1.0) + +! call openParameterBlock(param_file,'MOM_BACKGROUND_MIXING') + + call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", & + CS%Bryan_Lewis_diffusivity, & + "If true, use a Bryan & Lewis (JGR 1979) like tanh \n"//& + "profile of background diapycnal diffusivity with depth. \n"//& + "This is done via CVMix.", default=.false.) + + if (CS%Bryan_Lewis_diffusivity) then + + call get_param(param_file, mdl, "BRYAN_LEWIS_C1", & + CS%Bryan_Lewis_c1, & + "The vertical diffusivity values for Bryan-Lewis profile at |z|=D.", & + units="m2 s-1", fail_if_missing=.true.) + + call get_param(param_file, mdl, "BRYAN_LEWIS_C2", & + CS%Bryan_Lewis_c2, & + "The amplitude of variation in diffusivity for the Bryan-Lewis profile", & + units="m2 s-1", fail_if_missing=.true.) + + call get_param(param_file, mdl, "BRYAN_LEWIS_C3", & + CS%Bryan_Lewis_c3, & + "The inverse length scale for transition region in the Bryan-Lewis profile", & + units="m-1", fail_if_missing=.true.) + + call get_param(param_file, mdl, "BRYAN_LEWIS_C4", & + CS%Bryan_Lewis_c4, & + "The depth where diffusivity is BRYAN_LEWIS_C1 in the Bryan-Lewis profile",& + units="m", fail_if_missing=.true.) + + endif ! CS%Bryan_Lewis_diffusivity + + call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND", & + CS%Henyey_IGW_background, & + "If true, use a latitude-dependent scaling for the near \n"//& + "surface background diffusivity, as described in \n"//& + "Harrison & Hallberg, JPO 2008.", default=.false.) + + call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND_NEW", & + CS%Henyey_IGW_background_new, & + "If true, use a better latitude-dependent scaling for the\n"//& + "background diffusivity, as described in \n"//& + "Harrison & Hallberg, JPO 2008.", default=.false.) + + if (CS%Henyey_IGW_background .and. CS%Henyey_IGW_background_new) & + call MOM_error(FATAL, "set_diffusivity_init: HENYEY_IGW_BACKGROUND and \n"//& + "HENYEY_IGW_BACKGROUND_NEW are mutually exclusive. Set only one or none.") + + if (CS%Henyey_IGW_background) & + call get_param(param_file, mdl, "HENYEY_N0_2OMEGA", CS%N0_2Omega, & + "The ratio of the typical Buoyancy frequency to twice \n"//& + "the Earth's rotation period, used with the Henyey \n"//& + "scaling from the mixing.", units="nondim", default=20.0) + + call get_param(param_file, mdl, "KD_TANH_LAT_FN", & + CS%Kd_tanh_lat_fn, & + "If true, use a tanh dependence of Kd_sfc on latitude, \n"//& + "like CM2.1/CM2M. There is no physical justification \n"//& + "for this form, and it can not be used with \n"//& + "HENYEY_IGW_BACKGROUND.", default=.false.) + + if (CS%Kd_tanh_lat_fn) & + call get_param(param_file, mdl, "KD_TANH_LAT_SCALE", & + CS%Kd_tanh_lat_scale, & + "A nondimensional scaling for the range ofdiffusivities \n"//& + "with KD_TANH_LAT_FN. Valid values are in the range of \n"//& + "-2 to 2; 0.4 reproduces CM2M.", units="nondim", default=0.0) + + if (CS%Henyey_IGW_background .and. CS%Kd_tanh_lat_fn) call MOM_error(FATAL, & + "MOM_bkgnd_mixing: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.") + +! call closeParameterBlock(param_file) + + ! allocate arrays and set them to zero + allocate(CS%kd_bkgnd(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kd_bkgnd(:,:,:) = 0. + allocate(CS%kv_bkgnd(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kv_bkgnd(:,:,:) = 0. + allocate(CS%Kd_sfc(SZI_(G), SZJ_(G))); CS%Kd_sfc(:,:) = 0. + + ! Register diagnostics + CS%diag => diag + CS%id_kd_bkgnd = register_diag_field('ocean_model', 'Kd_bkgnd', diag%axesTi, Time, & + 'Background diffusivity added by MOM_bkgnd_mixing module', 'm2/s') + CS%id_kv_bkgnd = register_diag_field('ocean_model', 'Kv_bkgnd', diag%axesTi, Time, & + 'Background viscosity added by MOM_bkgnd_mixing module', 'm2/s') + +end subroutine bkgnd_mixing_init + +!> Get surface vertical background diffusivities/viscosities. +subroutine sfc_bkgnd_mixing(G, CS) + + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(bkgnd_mixing_cs), pointer, intent(inout) :: CS !< The control structure returned by + !! a previous call to bkgnd_mixing_init. + ! local variables + real :: I_x30 !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg))) + real :: deg_to_rad !< factor converting degrees to radians, pi/180. + real :: abs_sin !< absolute value of sine of latitude (nondim) + real :: epsilon + integer :: i, j, k, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + ! set some parameters + deg_to_rad = atan(1.0)/45.0 ! = PI/180 + epsilon = 1.e-10 + + + if (.not. CS%Bryan_Lewis_diffusivity) then +!$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc) + do j=js,je ; do i=is,ie + CS%Kd_sfc(i,j) = CS%Kd + enddo ; enddo + endif + + if (CS%Henyey_IGW_background) then + I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. +!$OMP parallel do default(none) +!shared(is,ie,js,je,Kd_sfc,CS,G,deg_to_rad,epsilon,I_x30) & +!$OMP private(abs_sin) + do j=js,je ; do i=is,ie + abs_sin = abs(sin(G%geoLatT(i,j)*deg_to_rad)) + CS%Kd_sfc(i,j) = max(CS%Kd_min, CS%Kd_sfc(i,j) * & + ((abs_sin * invcosh(CS%N0_2Omega/max(epsilon,abs_sin))) * I_x30) ) + enddo ; enddo + elseif (CS%Kd_tanh_lat_fn) then +!$OMP parallel do default(none) shared(is,ie,js,je,Kd_sfc,CS,G) + do j=js,je ; do i=is,ie + ! The transition latitude and latitude range are hard-scaled here, since + ! this is not really intended for wide-spread use, but rather for + ! comparison with CM2M / CM2.1 settings. + CS%Kd_sfc(i,j) = max(CS%Kd_min, CS%Kd_sfc(i,j) * (1.0 + & + CS%Kd_tanh_lat_scale * 0.5*tanh((abs(G%geoLatT(i,j)) - 35.0)/5.0) )) + enddo ; enddo + endif + + if (CS%debug) call hchksum(CS%Kd_sfc,"After sfc_bkgnd_mixing: Kd_sfc",G%HI,haloshift=0) + +end subroutine sfc_bkgnd_mixing + + +!> Calculates the vertical background diffusivities/viscosities +subroutine calculate_bkgnd_mixing(h, tv, N2_lay, kd_lay, kv, j, G, GV, CS) + + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. + real, dimension(SZI_(G),SZK_(G)), intent(in) :: N2_lay!< squared buoyancy frequency associated + !! with layers (1/s2) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: kd_lay!< Diapycnal diffusivity of each layer m2 s-1. + real, dimension(:,:,:), pointer :: kv !< The "slow" vertical viscosity at each interface + !! (not layer!) in m2 s-1. + integer, intent(in) :: j !< Meridional grid indice. + type(bkgnd_mixing_cs), pointer :: CS !< The control structure returned by + !! a previous call to bkgnd_mixing_init. + + ! local variables + real, dimension(SZI_(G), SZK_(G)+1) :: depth_2d !< distance from surface of an interface (m) + real, dimension(SZI_(G)) :: & + depth !< distance from surface of an interface (meter) + real :: depth_c !< depth of the center of a layer (meter) + real :: I_Hmix !< inverse of fixed mixed layer thickness (1/m) + real :: I_2Omega !< 1/(2 Omega) (sec) + real :: N_2Omega + real :: N02_N2 + real :: I_x30 !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg))) + real :: deg_to_rad !< factor converting degrees to radians, pi/180. + real :: abs_sin !< absolute value of sine of latitude (nondim) + real :: epsilon + integer :: i, k, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + + ! set some parameters + deg_to_rad = atan(1.0)/45.0 ! = PI/180 + epsilon = 1.e-10 + + depth_2d(:,:) = 0.0 + ! Set up the background diffusivity. + if (CS%Bryan_Lewis_diffusivity) then + + do i=is,ie + do k=2,nz+1 + depth_2d(i,k) = depth_2d(i,k-1) + GV%H_to_m*h(i,j,k-1) + enddo + ! if (is_root_pe()) write(*,*)'depth_3d(i,j,:)',depth_3d(i,j,:) + + call CVMix_init_bkgnd(max_nlev=nz, & + zw = depth_2d(i,:), & !< interface depth, must bepositive. + bl1 = CS%Bryan_Lewis_c1, & + bl2 = CS%Bryan_Lewis_c2, & + bl3 = CS%Bryan_Lewis_c3, & + bl4 = CS%Bryan_Lewis_c4, & + prandtl = CS%prandtl_bkgnd) + + call CVMix_coeffs_bkgnd(Mdiff_out=CS%kv_bkgnd(i,j,:), & + Tdiff_out=CS%kd_bkgnd(i,j,:), & + nlev=nz, & + max_nlev=nz) + + ! Update Kd + do k=1,nz + kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5*(CS%kd_bkgnd(i,j,K) + CS%kd_bkgnd(i,j,K+1)) + enddo + enddo ! i loop + + elseif ((.not. CS%Bryan_Lewis_diffusivity) .and. (.not.CS%bulkmixedlayer) .and. & + (CS%Kd/= CS%Kdml)) then + I_Hmix = 1.0 / CS%Hmix + do i=is,ie ; depth(i) = 0.0 ; enddo + do k=1,nz ; do i=is,ie + depth_c = depth(i) + 0.5*GV%H_to_m*h(i,j,k) + if (depth_c <= CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = CS%Kdml + elseif (depth_c >= 2.0*CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = CS%Kd_sfc(i,j) + else + kd_lay(i,j,k) = ((CS%Kd_sfc(i,j) - CS%Kdml) * I_Hmix) * depth_c + & + (2.0*CS%Kdml - CS%Kd_sfc(i,j)) + endif + + depth(i) = depth(i) + GV%H_to_m*h(i,j,k) + enddo ; enddo + + elseif (CS%Henyey_IGW_background_new) then + I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. + do k=1,nz ; do i=is,ie + abs_sin = max(epsilon,abs(sin(G%geoLatT(i,j)*deg_to_rad))) + N_2Omega = max(abs_sin,sqrt(N2_lay(i,k))*I_2Omega) + N02_N2 = (CS%N0_2Omega/N_2Omega)**2 + kd_lay(i,j,k) = max(CS%Kd_min, CS%Kd_sfc(i,j) * & + ((abs_sin * invcosh(N_2Omega/abs_sin)) * I_x30)*N02_N2) + enddo ; enddo + + else + do k=1,nz ; do i=is,ie + kd_lay(i,j,k) = CS%Kd_sfc(i,j) + enddo ; enddo + endif + + ! Update CS%kd_bkgnd and CS%kv_bkgnd for diagnostic purposes + if (.not. CS%Bryan_Lewis_diffusivity) then + do i=is,ie + CS%kd_bkgnd(i,j,1) = 0.0; CS%kv_bkgnd(i,j,1) = 0.0 + CS%kd_bkgnd(i,j,nz+1) = 0.0; CS%kv_bkgnd(i,j,nz+1) = 0.0 + do k=2,nz + CS%kd_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) + 0.5*(kd_lay(i,j,K-1) + kd_lay(i,j,K)) + CS%kv_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) * CS%prandtl_bkgnd + enddo + enddo + endif + + ! Update kv + if (associated(kv)) then + do i=is,ie + do k=1,nz+1 + kv(i,j,k) = kv(i,j,k) + CS%kv_bkgnd(i,j,k) + enddo + enddo + endif + +end subroutine calculate_bkgnd_mixing + +!> Reads the parameter "USE_CVMix_BACKGROUND" and returns state. +!! This function allows other modules to know whether this parameterization will +!! be used without needing to duplicate the log entry. +logical function CVMix_bkgnd_is_used(param_file) + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + call get_param(param_file, mdl, "USE_CVMix_BACKGROUND", CVMix_bkgnd_is_used, & + default=.false., do_not_log = .true.) + +end function CVMix_bkgnd_is_used + +!> Clear pointers and dealocate memory +subroutine bkgnd_mixing_end(CS) + type(bkgnd_mixing_cs), pointer :: CS ! Control structure + + deallocate(CS%kd_bkgnd) + deallocate(CS%kv_bkgnd) + deallocate(CS) + +end subroutine bkgnd_mixing_end + + +end module MOM_bkgnd_mixing diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 0239262ca6..b033f8bdfd 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -22,10 +22,12 @@ module MOM_diabatic_driver use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag, calc_Zint_diags use MOM_diapyc_energy_req, only : diapyc_energy_req_init, diapyc_energy_req_end use MOM_diapyc_energy_req, only : diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_CS -use MOM_diffConvection, only : diffConvection_CS, diffConvection_init -use MOM_diffConvection, only : diffConvection_calculate, diffConvection_end +use MOM_CVMix_conv, only : CVMix_conv_init, CVMix_conv_cs +use MOM_CVMix_conv, only : CVMix_conv_end, calculate_CVMix_conv use MOM_domains, only : pass_var, To_West, To_South, To_All, Omit_Corners use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type +use MOM_tidal_mixing, only : tidal_mixing_init, tidal_mixing_cs +use MOM_tidal_mixing, only : tidal_mixing_end use MOM_energetic_PBL, only : energetic_PBL, energetic_PBL_init use MOM_energetic_PBL, only : energetic_PBL_end, energetic_PBL_CS use MOM_energetic_PBL, only : energetic_PBL_get_MLD @@ -47,7 +49,7 @@ module MOM_diabatic_driver use MOM_internal_tides, only : propagate_int_tide use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS use MOM_kappa_shear, only : kappa_shear_is_used -use MOM_KPP, only : KPP_CS, KPP_init, KPP_calculate, KPP_end +use MOM_KPP, only : KPP_CS, KPP_init, KPP_calculate, KPP_end, KPP_get_BLD use MOM_KPP, only : KPP_NonLocalTransport_temp, KPP_NonLocalTransport_saln use MOM_opacity, only : opacity_init, set_opacity, opacity_end, opacity_CS use MOM_regularize_layers, only : regularize_layers, regularize_layers_init, regularize_layers_CS @@ -89,8 +91,11 @@ module MOM_diabatic_driver !! in the surface boundary layer. logical :: use_kappa_shear !< If true, use the kappa_shear module to find the !! shear-driven diapycnal diffusivity. - logical :: use_cvmix_shear !< If true, use the CVMix module to find the + logical :: use_CVMix_shear !< If true, use the CVMix module to find the !! shear-driven diapycnal diffusivity. + logical :: use_tidal_mixing !< If true, activate tidal mixing diffusivity. + logical :: use_CVMix_conv !< If true, use the CVMix module to get enhanced + !! mixing due to convection. logical :: use_sponge !< If true, sponges may be applied anywhere in the !! domain. The exact location and properties of !! those sponges are set by calls to @@ -147,11 +152,9 @@ module MOM_diabatic_driver real :: evap_CFL_limit = 0.8 !< The largest fraction of a layer that can be !! evaporated in one time-step (non-dim). - logical :: useKPP !< use CVmix/KPP diffusivities and non-local transport + logical :: useKPP = .false. !< use CVMix/KPP diffusivities and non-local transport logical :: salt_reject_below_ML !< If true, add salt below mixed layer (layer mode only) logical :: KPPisPassive !< If true, KPP is in passive mode, not changing answers. - logical :: useConvection !< If true, calculate large diffusivities when column - !! is statically unstable. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: debugConservation !< If true, monitor conservation and extrema. logical :: tracer_tridiag !< If true, use tracer_vertdiff instead of tridiagTS for @@ -221,7 +224,8 @@ module MOM_diabatic_driver type(optics_type), pointer :: optics => NULL() type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() type(KPP_CS), pointer :: KPP_CSp => NULL() - type(diffConvection_CS), pointer :: Conv_CSp => NULL() + type(tidal_mixing_cs), pointer :: tidal_mixing_csp => NULL() + type(CVMix_conv_cs), pointer :: CVMix_conv_csp => NULL() type(diapyc_energy_req_CS), pointer :: diapyc_en_rec_CSp => NULL() type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass @@ -584,8 +588,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G call cpu_clock_begin(id_clock_set_diffusivity) ! Sets: Kd, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S - ! Also changes: visc%Kd_turb, visc%TKE_turb (not clear that TKE_turb is used as input ????) - ! And sets visc%Kv_turb + ! Also changes: visc%Kd_shear, visc%TKE_turb (not clear that TKE_turb is used as input ???? + ! And sets visc%Kv_shear call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, CS%set_diff_CSp, Kd, Kd_int) call cpu_clock_end(id_clock_set_diffusivity) if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)") @@ -634,10 +638,15 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G !$OMP end parallel call KPP_calculate(CS%KPP_CSp, G, GV, h, tv%T, tv%S, u, v, tv%eqn_of_state, & - fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_turb, CS%KPP_NLTheat, & + fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, & CS%KPP_NLTscalar, Waves=Waves) !$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat) + if (associated(Hml)) then + call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G) + call pass_var(Hml, G%domain, halo=1) + endif + if (.not. CS%KPPisPassive) then !$OMP do do k=1,nz+1 ; do j=js,je ; do i=is,ie @@ -669,9 +678,18 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G endif ! endif for KPP - ! Check for static instabilities and increase Kd_int where unstable - if (CS%useConvection) call diffConvection_calculate(CS%Conv_CSp, & - G, GV, h, tv%T, tv%S, tv%eqn_of_state, Kd_int) + ! Add vertical diff./visc. due to convection (computed via CVMix) + if (CS%use_CVMix_conv) then + call calculate_CVMix_conv(h, tv, G, GV, CS%CVMix_conv_csp, Hml) + + !!!!!!!! GMM, the following needs to be checked !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + do k=1,nz ; do j=js,je ; do i=is,ie + Kd_int(i,j,k) = Kd_int(i,j,k) + CS%CVMix_conv_csp%kd_conv(i,j,k) + visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + CS%CVMix_conv_csp%kv_conv(i,j,k) + enddo ; enddo ; enddo + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + endif if (CS%useKPP) then @@ -816,10 +834,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G if (CS%ePBL_is_additive) then Kd_add_here = Kd_ePBL(i,j,K) - visc%Kv_turb(i,j,K) = visc%Kv_turb(i,j,K) + Kd_ePBL(i,j,K) + visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + Kd_ePBL(i,j,K) else - Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_turb(i,j,K), 0.0) - visc%Kv_turb(i,j,K) = max(visc%Kv_turb(i,j,K), Kd_ePBL(i,j,K)) + Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_shear(i,j,K), 0.0) + visc%Kv_shear(i,j,K) = max(visc%Kv_shear(i,j,K), Kd_ePBL(i,j,K)) endif Ent_int = Kd_add_here * (GV%m_to_H**2 * dt) / & (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) @@ -1352,9 +1370,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G call create_group_pass(CS%pass_hold_eb_ea, eb, G%Domain, dir_flag, halo=1) call create_group_pass(CS%pass_hold_eb_ea, ea, G%Domain, dir_flag, halo=1) call do_group_pass(CS%pass_hold_eb_ea, G%Domain) - ! visc%Kv_turb is not in the group pass because it has larger vertical extent. - if (associated(visc%Kv_turb)) & - call pass_var(visc%Kv_turb, G%Domain, To_All+Omit_Corners, halo=1) + ! visc%Kv_shear is not in the group pass because it has larger vertical extent. + if (associated(visc%Kv_shear)) & + call pass_var(visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) call cpu_clock_end(id_clock_pass) if (.not. CS%useALEalgorithm) then @@ -1908,7 +1926,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, apply parameterization of double-diffusion.", & default=.false. ) CS%use_kappa_shear = kappa_shear_is_used(param_file) - CS%use_CVMix_shear = cvmix_shear_is_used(param_file) + CS%use_CVMix_shear = CVMix_shear_is_used(param_file) if (CS%bulkmixedlayer) then call get_param(param_file, mod, "ML_MIX_FIRST", CS%ML_mix_first, & "The fraction of the mixed layer mixing that is applied \n"//& @@ -2328,10 +2346,12 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, allocate(CS%frazil_heat_diag(isd:ied,jsd:jed,nz) ) ; CS%frazil_heat_diag(:,:,:) = 0. endif + ! CS%use_tidal_mixing is set to True if an internal tidal dissipation scheme is to be used. + CS%use_tidal_mixing = tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, CS%tidal_mixing_CSp) - ! CS%useConvection is set to True IF convection will be used, otherwise False. - ! CS%Conv_CSp is allocated by diffConvection_init() - CS%useConvection = diffConvection_init(param_file, G, diag, Time, CS%Conv_CSp) + ! CS%use_CVMix_conv is set to True if CVMix convection will be used, otherwise + ! False. + CS%use_CVMix_conv = CVMix_conv_init(Time, G, GV, param_file, diag, CS%CVMix_conv_csp) call entrain_diffusive_init(Time, G, GV, param_file, diag, CS%entrain_diffusive_CSp) @@ -2347,7 +2367,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, endif ! initialize module for setting diffusivities - call set_diffusivity_init(Time, G, GV, param_file, diag, CS%set_diff_CSp, diag_to_Z_CSp, CS%int_tide_CSp) + call set_diffusivity_init(Time, G, GV, param_file, diag, CS%set_diff_CSp, diag_to_Z_CSp, CS%int_tide_CSp, CS%tidal_mixing_CSp) ! set up the clocks for this module @@ -2409,17 +2429,22 @@ subroutine diabatic_driver_end(CS) call entrain_diffusive_end(CS%entrain_diffusive_CSp) call set_diffusivity_end(CS%set_diff_CSp) - if (CS%useKPP) then + + if (CS%useKPP) then deallocate( CS%KPP_buoy_flux ) deallocate( CS%KPP_temp_flux ) deallocate( CS%KPP_salt_flux ) - endif + endif if (CS%useKPP) then deallocate( CS%KPP_NLTheat ) deallocate( CS%KPP_NLTscalar ) call KPP_end(CS%KPP_CSp) endif - if (CS%useConvection) call diffConvection_end(CS%Conv_CSp) + + if (CS%use_tidal_mixing) call tidal_mixing_end(CS%tidal_mixing_CSp) + + if (CS%use_CVMix_conv) call CVMix_conv_end(CS%CVMix_conv_csp) + if (CS%use_energetic_PBL) & call energetic_PBL_end(CS%energetic_PBL_CSp) if (CS%debug_energy_req) & @@ -2430,9 +2455,15 @@ subroutine diabatic_driver_end(CS) deallocate(CS%optics) endif - call diag_grid_storage_end(CS%diag_grids_prev) + ! GMM, the following is commented out because arrays in + ! CS%diag_grids_prev are neither pointers or allocatables + ! and, therefore, cannot be deallocated. + + !call diag_grid_storage_end(CS%diag_grids_prev) + if (associated(CS)) deallocate(CS) + end subroutine diabatic_driver_end diff --git a/src/parameterizations/vertical/MOM_diffConvection.F90 b/src/parameterizations/vertical/MOM_diffConvection.F90 deleted file mode 100644 index b63dbe4472..0000000000 --- a/src/parameterizations/vertical/MOM_diffConvection.F90 +++ /dev/null @@ -1,163 +0,0 @@ -module MOM_diffConvection - -! This file is part of MOM6. See LICENSE.md for the license. - -use MOM_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data -use MOM_diag_mediator, only : query_averaging_enabled, register_diag_field -use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_PE -use MOM_EOS, only : EOS_type, calculate_density -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_file_parser, only : openParameterBlock, closeParameterBlock -use MOM_grid, only : ocean_grid_type, isPointInCell -use MOM_verticalGrid, only : verticalGrid_type - -implicit none ; private - -#include "MOM_memory.h" - -public :: diffConvection_init, diffConvection_calculate, diffConvection_end - -! Control structure for containing KPP parameters/data -type, public :: diffConvection_CS ; private - - ! Parameters - real :: Kd_convection ! The value of diffusivity to add at statically unstable interfaces (m2/s) - logical :: debug ! If true, turn on debugging - logical :: passiveMode ! If true, make the motions but go nowhere - - ! Daignostic handles and pointers - type(diag_ctrl), pointer :: diag => NULL() - integer :: id_N2 = -1, id_Kd_conv = -1 - - ! Diagnostics arrays - real, allocatable, dimension(:,:,:) :: N2 ! Brunt-Vaisala frequency (1/s2) - real, allocatable, dimension(:,:,:) :: Kd_conv ! Diffusivity added by convection (m2/s) - -end type diffConvection_CS - -! Module data used for debugging only -logical, parameter :: verbose = .False. - -contains - -logical function diffConvection_init(paramFile, G, diag, Time, CS) -!< Initialize the CVmix KPP module and set up diagnostics -!! Returns True if module is to be used, otherwise returns False. - -! Arguments - type(param_file_type), intent(in) :: paramFile !< File parser - type(ocean_grid_type), intent(in) :: G !< Ocean grid - type(diag_ctrl), target, intent(in) :: diag !< Diagnostics - type(time_type), intent(in) :: Time !< Time - type(diffConvection_CS), pointer :: CS !< Control structure -! Local variables -#include "version_variable.h" - character(len=40) :: mdl = 'MOM_diffConvection' ! This module's name. - - if (associated(CS)) call MOM_error(FATAL, 'MOM_diffConvection, diffConvection_init: '// & - 'Control structure has already been initialized') - allocate(CS) - -! Read parameters - call log_version(paramFile, mdl, version, & - 'This module implements enhanced diffusivity as a\n' // & - 'function of static stability, N^2.') - call get_param(paramFile, mdl, "USE_CONVECTION", diffConvection_init, & - "If true, turns on the diffusive convection scheme that\n"// & - "increases diapycnal diffusivities at statically unstable\n"// & - "interfaces. Relevant parameters are contained in the\n"// & - "CONVECTION% parameter block.", & - default=.false.) - - call openParameterBlock(paramFile,'CONVECTION') - call get_param(paramFile, mdl, 'PASSIVE', CS%passiveMode, & - 'If True, puts KPP into a passive-diagnostic mode.', & - default=.False.) - call get_param(paramFile, mdl, 'KD_CONV', CS%Kd_convection, & - 'DIffusivity used in statically unstable regions of column.', & - units='m2/s', default=1.00) - call closeParameterBlock(paramFile) - call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) - -! Forego remainder of initialization if not using this scheme - if (.not. diffConvection_init) return - -! Register diagnostics - CS%diag => diag - CS%id_N2 = register_diag_field('ocean_model', 'Conv_N2', diag%axesTi, Time, & - 'Square of Brunt-Vaisala frequency used by diffConvection module', '1/s2') - if (CS%id_N2 > 0) allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) - CS%id_Kd_conv = register_diag_field('ocean_model', 'Conv_Kd', diag%axesTi, Time, & - 'Additional diffusivity added by diffConvection module', 'm2/s') - if (CS%id_Kd_conv > 0) allocate( CS%Kd_conv( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) - - if (CS%id_N2 > 0) CS%N2(:,:,:) = 0. - if (CS%id_Kd_conv > 0) CS%Kd_conv(:,:,:) = 0. - -end function diffConvection_init - - -subroutine diffConvection_calculate(CS, G, GV, h, Temp, Salt, EOS, Kd_int) -!< Calculates diffusivity and non-local transport for KPP parameterization - -! Arguments - type(diffConvection_CS), pointer :: CS !< Control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses (units of H) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Temp !< Pot. temperature (degrees C) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Salt !< Salinity (ppt) - type(EOS_type), pointer :: EOS !< Equation of state - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kd_int !< (in) Vertical diffusivity on interfaces (m2/s) - !! (out) Modified vertical diffusivity (m2/s) -! Local variables - integer :: i, j, k - real, dimension( G%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces (1/s2) - real, dimension( G%ke+1 ) :: Kd_1d ! Vertical diffusivity at interfaces (m2/s) - real :: GoRho, pRef, rhoK, rhoKm1 - - GoRho = GV%g_Earth / GV%Rho0 - - N2_1d( 1 ) = 0. - N2_1d( G%ke+1 ) = 0. - Kd_1d( 1 ) = 0. - Kd_1d( G%ke+1 ) = 0. - do j=G%jsc,G%jec ; do i=G%isc,G%iec - ! This k-loop calculates external quantities independent of any iterations - ! Start at bottom of top level - pRef = 0. ! Ignore atmospheric pressure - do K = 2, G%ke - ! Pressure at interface K is incremented by mass of level above - pRef = pRef + GV%g_Earth * GV%Rho0 * h(i,j,k-1) * GV%H_to_m ! Boussinesq approximation!!!! ????? - ! Compute Brunt-Vaisala frequency (static stability) on interfaces - call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS) - call calculate_density(Temp(i,j,k-1), Salt(i,j,k-1), pRef, rhoKm1, EOS) - N2_1d(K) = GoRho * (rhoK - rhoKm1) / & - (0.5*(h(i,j,k-1) + h(i,j,k)) + GV%H_subroundoff) ! Can be negative - Kd_1d(K) = 0. - if (N2_1d(K) < 0.) Kd_1d(K) = CS%Kd_convection - enddo ! k - - if (.not. CS%passiveMode) Kd_int(i,j,:) = Kd_int(i,j,:) + Kd_1d(:) - - if (CS%id_N2 > 0) CS%N2(i,j,:) = N2_1d(:) - if (CS%id_Kd_conv > 0) CS%Kd_conv(i,j,:) = Kd_1d(:) - - enddo ; enddo ! j - - if (CS%id_N2 > 0) call post_data(CS%id_N2, CS%N2, CS%diag) - if (CS%id_Kd_conv > 0) call post_data(CS%id_Kd_conv, CS%Kd_conv, CS%diag) - -end subroutine diffConvection_calculate - - -subroutine diffConvection_end(CS) -! Clear pointers, dealocate memory - type(diffConvection_CS), pointer :: CS ! Control structure - - if (CS%id_N2 > 0) deallocate(CS%N2, CS%diag) - if (CS%id_Kd_conv > 0) deallocate(CS%Kd_conv, CS%diag) - deallocate(CS) -end subroutine diffConvection_end - -end module MOM_diffConvection diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 96a150500f..3344f218bc 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -71,7 +71,7 @@ module MOM_kappa_shear real :: TKE_bg ! The background level of TKE, in m2 s-2. real :: kappa_0 ! The background diapycnal diffusivity, in m2 s-1. real :: kappa_tol_err ! The fractional error in kappa that is tolerated. - real :: Prandtl_turb ! Prandtl number used to convert Kd_turb into viscosity. + real :: Prandtl_turb ! Prandtl number used to convert Kd_shear into viscosity. integer :: nkml ! The number of layers in the mixed layer, as ! treated in this routine. If the pieces of the ! mixed layer are not to be treated collectively, @@ -130,7 +130,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & intent(inout) :: kv_io !< The vertical viscosity at each interface !! (not layer!) in m2 s-1. This discards any !! previous value i.e. intent(out) and simply - !! sets Kv = Prandtl * Kd_turb + !! sets Kv = Prandtl * Kd_shear real, intent(in) :: dt !< Time increment, in s. type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous !! call to kappa_shear_init. @@ -156,7 +156,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & ! the iteration toward convergence. ! (in/out) kv_io - The vertical viscosity at each interface ! (not layer!) in m2 s-1. This discards any previous value -! i.e. intent(out) and simply sets Kv = Prandtl * Kd_turb +! i.e. intent(out) and simply sets Kv = Prandtl * Kd_shear ! (in) dt - Time increment, in s. ! (in) G - The ocean's grid structure. ! (in) GV - The ocean's vertical grid structure. diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 0a707a5487..2e9eaca122 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -2,32 +2,10 @@ module MOM_set_diffusivity ! This file is part of MOM6. See LICENSE.md for the license. -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, September 1997 - June 2007 * -!* * -!* This file contains the subroutines that sets the diapycnal * -!* diffusivity, perhaps adding up pieces that are calculated in other * -!* files and passed in via the vertvisc type argument. * -!* * -!* A small fragment of the grid is shown below: * -!* * -!* j+1 x ^ x ^ x At x: q * -!* j+1 > o > o > At ^: v * -!* j x ^ x ^ x At >: u * -!* j > o > o > At o: h, buoy, ustar, T, S, Kd, ea, eb, etc. * -!* j-1 x ^ x ^ x * -!* i-1 i i+1 At x & ^: * -!* i i+1 At > & o: * -!* * -!* The boundaries always run through q grid points (x). * -!* * -!********+*********+*********+*********+*********+*********+*********+** - use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, time_type -use MOM_diag_mediator, only : safe_alloc_ptr, post_data, register_diag_field +use MOM_diag_mediator, only : post_data, register_diag_field use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag, calc_Zint_diags use MOM_debugging, only : hchksum, uvchksum use MOM_EOS, only : calculate_density, calculate_density_derivs @@ -38,15 +16,19 @@ module MOM_set_diffusivity use MOM_forcing_type, only : forcing, optics_type use MOM_grid, only : ocean_grid_type use MOM_internal_tides, only : int_tide_CS, get_lowmode_loss +use MOM_tidal_mixing, only : tidal_mixing_CS, calculate_tidal_mixing +use MOM_tidal_mixing, only : setup_tidal_diagnostics, post_tidal_diagnostics use MOM_intrinsic_functions, only : invcosh use MOM_io, only : slasher, vardesc, var_desc, MOM_read_data use MOM_kappa_shear, only : calculate_kappa_shear, kappa_shear_init, Kappa_shear_CS -use MOM_cvmix_shear, only : calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_CS +use MOM_CVMix_shear, only : calculate_CVMix_shear, CVMix_shear_init, CVMix_shear_cs +use MOM_CVMix_shear, only : CVMix_shear_end +use MOM_bkgnd_mixing, only : calculate_bkgnd_mixing, bkgnd_mixing_init, bkgnd_mixing_cs +use MOM_bkgnd_mixing, only : bkgnd_mixing_end, sfc_bkgnd_mixing use MOM_string_functions, only : uppercase use MOM_thickness_diffuse, only : vert_fill_TS use MOM_variables, only : thermo_var_ptrs, vertvisc_type, p3d -use MOM_verticalGrid, only : verticalGrid_type - +use MOM_verticalGrid, only : verticalGrid_type use user_change_diffusivity, only : user_change_diff, user_change_diff_init use user_change_diffusivity, only : user_change_diff_end, user_change_diff_CS @@ -70,44 +52,6 @@ module MOM_set_diffusivity ! large enough that N2 > omega2. The full expression for ! the Flux Richardson number is usually ! FLUX_RI_MAX*N2/(N2+OMEGA2). The default is 0.2. - logical :: Henyey_IGW_background ! If true, use a simplified variant of the - ! Henyey et al, JGR (1986) latitudinal scaling for - ! the background diapycnal diffusivity, which gives - ! a marked decrease in the diffusivity near the - ! equator. The simplification here is to assume - ! that the in-situ stratification is the same as - ! the reference stratificaiton. - logical :: Henyey_IGW_background_new ! same as Henyey_IGW_background - ! but incorporate the effect of - ! stratification on TKE dissipation, - ! - ! e = f/f_0 * acosh(N/f) / acosh(N_0/f_0) * e_0 - ! - ! where e is the TKE dissipation, and N_0 and f_0 are the - ! reference buoyancy frequency and inertial frequencies respectively. - ! e_0 is the reference dissipation at (N_0,f_0). In the - ! previous version, N=N_0. - ! - ! Additionally, the squared inverse relationship between - ! diapycnal diffusivities and stratification is included - ! - ! kd = e/N^2 - ! - ! where kd is the diapycnal diffusivity. - ! This approach assumes that work done - ! against gravity is uniformly distributed - ! throughout the column. Whereas, kd=kd_0*e, - ! as in the original version, concentrates buoyancy - ! work in regions of strong stratification. - - logical :: Kd_tanh_lat_fn ! If true, use the tanh dependence of Kd_sfc on - ! latitude, like GFDL CM2.1/CM2M. There is no physical - ! justification for this form, and it can not be - ! used with Henyey_IGW_background. - real :: Kd_tanh_lat_scale ! A nondimensional scaling for the range of - ! diffusivities with Kd_tanh_lat_fn. Valid values - ! are in the range of -2 to 2; 0.4 reproduces CM2M. - logical :: bottomdraglaw ! If true, the bottom stress is calculated with a ! drag law c_drag*|u|*u. logical :: BBL_mixing_as_max ! If true, take the maximum of the diffusivity @@ -133,65 +77,8 @@ module MOM_set_diffusivity ! when bulkmixedlayer==.false. real :: Hmix ! mixed layer thickness (meter) when ! bulkmixedlayer==.false. - - logical :: Bryan_Lewis_diffusivity ! If true, background vertical diffusivity - ! uses Bryan-Lewis (1979) like tanh profile. - real :: Kd_Bryan_Lewis_deep ! abyssal value of Bryan-Lewis profile (m2/s) - real :: Kd_Bryan_Lewis_surface ! surface value of Bryan-Lewis profile (m2/s) - real :: Bryan_Lewis_depth_cent ! center of transition depth in Bryan-Lewis (meter) - real :: Bryan_Lewis_width_trans ! width of transition for Bryan-Lewis (meter) - - real :: N0_2Omega ! ratio of the typical Buoyancy frequency to - ! twice the Earth's rotation period, used with the - ! Henyey scaling from the mixing - real :: N2_FLOOR_IOMEGA2 ! floor applied to N2(k) scaled by Omega^2 - ! If =0., N2(k) is positive definite - ! If =1., N2(k) > Omega^2 everywhere - type(diag_ctrl), pointer :: diag ! structure to regulate diagn output timing - real :: Int_tide_decay_scale ! decay scale for internal wave TKE (meter) - real :: Mu_itides ! efficiency for conversion of dissipation - ! to potential energy (nondimensional) - real :: Gamma_itides ! fraction of local dissipation (nondimensional) - real :: Gamma_lee ! fraction of local dissipation for lee waves - ! (Nikurashin's energy input) (nondimensional) - real :: Decay_scale_factor_lee ! Scaling factor for the decay scale of lee - ! wave energy dissipation (nondimensional) - real :: min_zbot_itides ! minimum depth for internal tide conversion (meter) - logical :: Int_tide_dissipation ! Internal tide conversion (from barotropic) with - ! the schemes of St Laurent et al (2002)/ - ! Simmons et al (2004) - logical :: Lowmode_itidal_dissipation ! Internal tide conversion (from low modes) with - ! the schemes of St Laurent et al (2002)/ - ! Simmons et al (2004) !BDM - integer :: Int_tide_profile ! A coded integer indicating the vertical profile - ! for dissipation of the internal waves. Schemes that - ! are currently encoded are St Laurent et al (2002) and - ! Polzin (2009). - real :: Nu_Polzin ! The non-dimensional constant used in Polzin form of - ! the vertical scale of decay of tidal dissipation - real :: Nbotref_Polzin ! Reference value for the buoyancy frequency at the - ! ocean bottom used in Polzin formulation of the - ! vertical scale of decay of tidal dissipation (1/s) - real :: Polzin_decay_scale_factor ! Scaling factor for the decay length scale - ! of the tidal dissipation profile in Polzin - ! (nondimensional) - real :: Polzin_decay_scale_max_factor ! The decay length scale of tidal - ! dissipation profile in Polzin formulation should not - ! exceed Polzin_decay_scale_max_factor * depth of the - ! ocean (nondimensional). - real :: Polzin_min_decay_scale ! minimum decay scale of the tidal dissipation - ! profile in Polzin formulation (meter) - logical :: Lee_wave_dissipation ! Enable lee-wave driven mixing, following - ! Nikurashin (2010), with a vertical energy - ! deposition profile specified by Lee_wave_profile. - ! St Laurent et al (2002) or - ! Simmons et al (2004) scheme - integer :: Lee_wave_profile ! A coded integer indicating the vertical profile - ! for dissipation of the lee waves. Schemes that are - ! currently encoded are St Laurent et al (2002) and - ! Polzin (2009). logical :: limit_dissipation ! If enabled, dissipation is limited to be larger ! than the following: real :: dissip_min ! Minimum dissipation (W/m3) @@ -203,10 +90,6 @@ module MOM_set_diffusivity real :: TKE_itide_max ! maximum internal tide conversion (W m-2) ! available to mix above the BBL real :: omega ! Earth's rotation frequency (s-1) - real :: utide ! constant tidal amplitude (m s-1) used if - ! tidal amplitude file is not present - real :: kappa_itides ! topographic wavenumber and non-dimensional scaling - real :: kappa_h2_factor ! factor for the product of wavenumber * rms sgs height logical :: ML_radiation ! allow a fraction of TKE available from wind work ! to penetrate below mixed layer base with a vertical ! decay scale determined by the minimum of @@ -244,62 +127,36 @@ module MOM_set_diffusivity logical :: user_change_diff ! If true, call user-defined code to change diffusivity. logical :: useKappaShear ! If true, use the kappa_shear module to find the ! shear-driven diapycnal diffusivity. - - logical :: useCVmix ! If true, use one of the CVMix modules to find + logical :: use_CVMix_shear ! If true, use one of the CVMix modules to find ! shear-driven diapycnal diffusivity. - - logical :: double_diffusion ! If true, enable double-diffusive mixing. + logical :: double_diffusion ! If true, enable double-diffusive mixing. logical :: simple_TKE_to_Kd ! If true, uses a simple estimate of Kd/TKE that ! does not rely on a layer-formulation. real :: Max_Rrho_salt_fingers ! max density ratio for salt fingering real :: Max_salt_diff_salt_fingers ! max salt diffusivity for salt fingers (m2/s) real :: Kv_molecular ! molecular visc for double diff convect (m2/s) - real, pointer, dimension(:,:) :: TKE_Niku => NULL() - real, pointer, dimension(:,:) :: TKE_itidal => NULL() - real, pointer, dimension(:,:) :: Nb => NULL() - real, pointer, dimension(:,:) :: mask_itidal => NULL() - real, pointer, dimension(:,:) :: h2 => NULL() - real, pointer, dimension(:,:) :: tideamp => NULL() ! RMS tidal amplitude (m/s) - character(len=200) :: inputdir type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() type(Kappa_shear_CS), pointer :: kappaShear_CSp => NULL() - type(CVMix_shear_CS), pointer :: CVMix_Shear_CSp => NULL() + type(CVMix_shear_cs), pointer :: CVMix_shear_csp => NULL() + type(bkgnd_mixing_cs), pointer :: bkgnd_mixing_csp => NULL() type(int_tide_CS), pointer :: int_tide_CSp => NULL() + type(tidal_mixing_cs), pointer :: tm_csp => NULL() - integer :: id_TKE_itidal = -1 - integer :: id_TKE_leewave = -1 integer :: id_maxTKE = -1 integer :: id_TKE_to_Kd = -1 - integer :: id_Kd_itidal = -1 - integer :: id_Kd_Niku = -1 - integer :: id_Kd_lowmode = -1 integer :: id_Kd_user = -1 integer :: id_Kd_layer = -1 integer :: id_Kd_BBL = -1 integer :: id_Kd_BBL_z = -1 - integer :: id_Kd_itidal_z = -1 - integer :: id_Kd_Niku_z = -1 - integer :: id_Kd_lowmode_z = -1 integer :: id_Kd_user_z = -1 integer :: id_Kd_Work = -1 - integer :: id_Kd_Itidal_Work = -1 - integer :: id_Kd_Niku_Work = -1 - integer :: id_Kd_Lowmode_Work= -1 - - integer :: id_Fl_itidal = -1 - integer :: id_Fl_lowmode = -1 - integer :: id_Polzin_decay_scale = -1 - integer :: id_Polzin_decay_scale_scaled = -1 - integer :: id_Nb = -1 integer :: id_N2 = -1 integer :: id_N2_z = -1 - integer :: id_N2_bot = -1 - integer :: id_N2_meanz = -1 integer :: id_KT_extra = -1 integer :: id_KS_extra = -1 @@ -311,19 +168,9 @@ module MOM_set_diffusivity type diffusivity_diags real, pointer, dimension(:,:,:) :: & N2_3d => NULL(),& ! squared buoyancy frequency at interfaces (1/s2) - Kd_itidal => NULL(),& ! internal tide diffusivity at interfaces (m2/s) - Fl_itidal => NULL(),& ! vertical flux of tidal turbulent dissipation (m3/s3) - Kd_lowmode => NULL(),& ! internal tide diffusivity at interfaces - ! due to propagating low modes (m2/s) (BDM) - Fl_lowmode => NULL(),& ! vertical flux of tidal turbulent dissipation - ! due to propagating low modes (m3/s3) (BDM) - Kd_Niku => NULL(),& ! lee-wave diffusivity at interfaces (m2/s) Kd_user => NULL(),& ! user-added diffusivity at interfaces (m2/s) Kd_BBL => NULL(),& ! BBL diffusivity at interfaces (m2/s) Kd_work => NULL(),& ! layer integrated work by diapycnal mixing (W/m2) - Kd_Niku_work => NULL(),& ! layer integrated work by lee-wave driven mixing (W/m2) - Kd_Itidal_Work => NULL(),& ! layer integrated work by int tide driven mixing (W/m2) - Kd_Lowmode_Work=> NULL(),& ! layer integrated work by low mode driven mixing (W/m2) BDM maxTKE => NULL(),& ! energy required to entrain to h_max (m3/s3) TKE_to_Kd => NULL(),& ! conversion rate (~1.0 / (G_Earth + dRho_lay)) ! between TKE dissipated within a layer and Kd @@ -331,20 +178,8 @@ module MOM_set_diffusivity KT_extra => NULL(),& ! double diffusion diffusivity for temp (m2/s) KS_extra => NULL() ! double diffusion diffusivity for saln (m2/s) - real, pointer, dimension(:,:) :: & - TKE_itidal_used => NULL(),& ! internal tide TKE input at ocean bottom (W/m2) - N2_bot => NULL(),& ! bottom squared buoyancy frequency (1/s2) - N2_meanz => NULL(),& ! vertically averaged buoyancy frequency (1/s2) - Polzin_decay_scale_scaled => NULL(),& ! vertical scale of decay for tidal dissipation - Polzin_decay_scale => NULL() ! vertical decay scale for tidal diss with Polzin (meter) - end type diffusivity_diags -character*(20), parameter :: STLAURENT_PROFILE_STRING = "STLAURENT_02" -character*(20), parameter :: POLZIN_PROFILE_STRING = "POLZIN_09" -integer, parameter :: STLAURENT_02 = 1 -integer, parameter :: POLZIN_09 = 2 - ! Clocks integer :: id_clock_kappaShear @@ -380,29 +215,11 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & optional, intent(out) :: Kd_int !< Diapycnal diffusivity at each interface !! (m2/sec). -! Arguments: -! (in) u - zonal velocity (m/s) -! (in) v - meridional velocity (m/s) -! (in) h - Layer thickness (m or kg/m2) -! (in) tv - structure with pointers to thermodynamic fields -! (in) fluxes - structure of surface fluxes that may be used -! (in) visc - structure containing vertical viscosities, bottom boundary -! layer properies, and related fields -! (in) dt - time increment (sec) -! (in) G - ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) CS - module control structure -! (in) j - meridional index upon which to work -! (out) Kd - diapycnal diffusivity of each layer (m2/sec) -! (out,opt) Kd_int - diapycnal diffusivity at each interface (m2/sec) - + ! local variables real, dimension(SZI_(G)) :: & - depth, & ! distance from surface of an interface (meter) N2_bot ! bottom squared buoyancy frequency (1/s2) - real, dimension(SZI_(G), SZJ_(G)) :: & - Kd_sfc ! surface value of the diffusivity (m2/s) - type(diffusivity_diags) :: dd ! structure w/ arrays of pointers to avail diags + type(diffusivity_diags) :: dd ! structure w/ arrays of pointers to avail diags real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & T_f, S_f ! temperature and salinity (deg C and ppt); @@ -421,22 +238,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & KT_extra, & ! double difusion diffusivity on temperature (m2/sec) KS_extra ! double difusion diffusivity on salinity (m2/sec) - real :: I_trans ! inverse of the transitional for Bryan-Lewis (1/m) - real :: depth_c ! depth of the center of a layer (meter) - real :: I_Hmix ! inverse of fixed mixed layer thickness (1/m) real :: I_Rho0 ! inverse of Boussinesq density (m3/kg) - real :: I_x30 ! 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg))) - real :: abs_sin ! absolute value of sine of latitude (nondim) - real :: atan_fn_sfc ! surface value of Bryan-Lewis profile (nondim) - real :: atan_fn_lay ! value of Bryan-Lewis profile in layer middle (nondim) - real :: I_atan_fn ! inverse of change in Bryan-Lewis profile from surface to infinite depth (nondim) - real :: deg_to_rad ! factor converting degrees to radians, pi/180. real :: dissip ! local variable for dissipation calculations (W/m3) real :: Omega2 ! squared absolute rotation rate (1/s2) - real :: I_2Omega ! 1/(2 Omega) (sec) - real :: N_2Omega - real :: N02_N2 - real :: epsilon logical :: use_EOS ! If true, compute density from T/S using equation of state. type(p3d) :: z_ptrs(6) ! pointers to diagns to be interpolated into depth space @@ -463,10 +267,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & I_Rho0 = 1.0/GV%Rho0 kappa_fill = 1.e-3 ! m2 s-1 dt_fill = 7200. - deg_to_rad = atan(1.0)/45.0 ! = PI/180 Omega2 = CS%Omega*CS%Omega - I_2Omega = 0.5/CS%Omega - epsilon = 1.e-10 use_EOS = associated(tv%eqn_of_state) @@ -480,58 +281,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if ((CS%id_N2 > 0) .or. (CS%id_N2_z > 0)) then allocate(dd%N2_3d(isd:ied,jsd:jed,nz+1)) ; dd%N2_3d(:,:,:) = 0.0 endif - if ((CS%id_Kd_itidal > 0) .or. (CS%id_Kd_itidal_z > 0) .or. & - (CS%id_Kd_Itidal_work > 0)) then - allocate(dd%Kd_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Kd_itidal(:,:,:) = 0.0 - endif - if ((CS%id_Kd_lowmode > 0) .or. (CS%id_Kd_lowmode_z > 0) .or. & - (CS%id_Kd_lowmode_work > 0)) then - allocate(dd%Kd_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Kd_lowmode(:,:,:) = 0.0 - endif - if ( (CS%id_Fl_itidal > 0) ) then - allocate(dd%Fl_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Fl_itidal(:,:,:) = 0.0 - endif - if ( (CS%id_Fl_lowmode > 0) ) then - allocate(dd%Fl_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Fl_lowmode(:,:,:) = 0.0 - endif - if ( (CS%id_Polzin_decay_scale > 0) ) then - allocate(dd%Polzin_decay_scale(isd:ied,jsd:jed)) - dd%Polzin_decay_scale(:,:) = 0.0 - endif - if ( (CS%id_Polzin_decay_scale_scaled > 0) ) then - allocate(dd%Polzin_decay_scale_scaled(isd:ied,jsd:jed)) - dd%Polzin_decay_scale_scaled(:,:) = 0.0 - endif - if ( (CS%id_N2_bot > 0) ) then - allocate(dd%N2_bot(isd:ied,jsd:jed)) ; dd%N2_bot(:,:) = 0.0 - endif - if ( (CS%id_N2_meanz > 0) ) then - allocate(dd%N2_meanz(isd:ied,jsd:jed)) ; dd%N2_meanz(:,:) = 0.0 - endif - if ((CS%id_Kd_Niku > 0) .or. (CS%id_Kd_Niku_z > 0) .or. & - (CS%id_Kd_Niku_work > 0)) then - allocate(dd%Kd_Niku(isd:ied,jsd:jed,nz+1)) ; dd%Kd_Niku(:,:,:) = 0.0 - endif if ((CS%id_Kd_user > 0) .or. (CS%id_Kd_user_z > 0)) then allocate(dd%Kd_user(isd:ied,jsd:jed,nz+1)) ; dd%Kd_user(:,:,:) = 0.0 endif if (CS%id_Kd_work > 0) then allocate(dd%Kd_work(isd:ied,jsd:jed,nz)) ; dd%Kd_work(:,:,:) = 0.0 endif - if (CS%id_Kd_Niku_work > 0) then - allocate(dd%Kd_Niku_work(isd:ied,jsd:jed,nz)) ; dd%Kd_Niku_work(:,:,:) = 0.0 - endif - if (CS%id_Kd_Itidal_work > 0) then - allocate(dd%Kd_Itidal_work(isd:ied,jsd:jed,nz)) - dd%Kd_Itidal_work(:,:,:) = 0.0 - endif - if (CS%id_Kd_Lowmode_Work > 0) then - allocate(dd%Kd_Lowmode_Work(isd:ied,jsd:jed,nz)) - dd%Kd_Lowmode_Work(:,:,:) = 0.0 - endif - if (CS%id_TKE_itidal > 0) then - allocate(dd%TKE_Itidal_used(isd:ied,jsd:jed)) ; dd%TKE_Itidal_used(:,:) = 0. - endif if (CS%id_maxTKE > 0) then allocate(dd%maxTKE(isd:ied,jsd:jed,nz)) ; dd%maxTKE(:,:,:) = 0.0 endif @@ -548,6 +303,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0 endif + ! set up arrays for tidal mixing diagnostics + call setup_tidal_diagnostics(G,CS%tm_csp) + ! Smooth the properties through massless layers. if (use_EOS) then if (CS%debug) then @@ -569,60 +327,33 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & call hchksum(v_h, "before calc_KS v_h",G%HI) endif call cpu_clock_begin(id_clock_kappaShear) - ! Changes: visc%Kd_turb, visc%TKE_turb (not clear that TKE_turb is used as input ????) - ! Sets visc%Kv_turb - call calculate_kappa_shear(u_h, v_h, h, tv, fluxes%p_surf, visc%Kd_turb, visc%TKE_turb, & - visc%Kv_turb, dt, G, GV, CS%kappaShear_CSp) + ! Changes: visc%Kd_shear, visc%TKE_turb (not clear that TKE_turb is used as input ????) + ! Sets visc%Kv_shear + call calculate_kappa_shear(u_h, v_h, h, tv, fluxes%p_surf, visc%Kd_shear, visc%TKE_turb, & + visc%Kv_shear, dt, G, GV, CS%kappaShear_CSp) call cpu_clock_end(id_clock_kappaShear) if (CS%debug) then - call hchksum(visc%Kd_turb, "after calc_KS visc%Kd_turb",G%HI) - call hchksum(visc%Kv_turb, "after calc_KS visc%Kv_turb",G%HI) + call hchksum(visc%Kd_shear, "after calc_KS visc%Kd_shear",G%HI) + call hchksum(visc%Kv_shear, "after calc_KS visc%Kv_shear",G%HI) call hchksum(visc%TKE_turb, "after calc_KS visc%TKE_turb",G%HI) endif if (showCallTree) call callTree_waypoint("done with calculate_kappa_shear (set_diffusivity)") - elseif (CS%useCVMix) then + elseif (CS%use_CVMix_shear) then !NOTE{BGR}: this needs to be cleaned up. It works in 1D case, but has not been tested outside. - call calculate_cvmix_shear(u_h, v_h, h, tv, visc%Kd_turb, visc%Kv_turb,G,GV,CS%CVMix_shear_CSp) - elseif (associated(visc%Kv_turb)) then - visc%Kv_turb(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled + call calculate_CVMix_shear(u_h, v_h, h, tv, visc%Kd_shear, visc%Kv_shear,G,GV,CS%CVMix_shear_csp) + elseif (associated(visc%Kv_shear)) then + visc%Kv_shear(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled endif -! Calculate the diffusivity, Kd, for each layer. This would be -! the appropriate place to add a depth-dependent parameterization or -! another explicit parameterization of Kd. + ! Calculate the diffusivity, Kd, for each layer. This would be + ! the appropriate place to add a depth-dependent parameterization or + ! another explicit parameterization of Kd. - if (CS%Bryan_Lewis_diffusivity) then -!$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc) - do j=js,je ; do i=is,ie - Kd_sfc(i,j) = CS%Kd_Bryan_Lewis_surface - enddo ; enddo - else -!$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc) - do j=js,je ; do i=is,ie - Kd_sfc(i,j) = CS%Kd - enddo ; enddo - endif - if (CS%Henyey_IGW_background) then - I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. -!$OMP parallel do default(none) shared(is,ie,js,je,Kd_sfc,CS,G,deg_to_rad,epsilon,I_x30) & -!$OMP private(abs_sin) - do j=js,je ; do i=is,ie - abs_sin = abs(sin(G%geoLatT(i,j)*deg_to_rad)) - Kd_sfc(i,j) = max(CS%Kd_min, Kd_sfc(i,j) * & - ((abs_sin * invcosh(CS%N0_2Omega/max(epsilon,abs_sin))) * I_x30) ) - enddo ; enddo - elseif (CS%Kd_tanh_lat_fn) then -!$OMP parallel do default(none) shared(is,ie,js,je,Kd_sfc,CS,G) - do j=js,je ; do i=is,ie - ! The transition latitude and latitude range are hard-scaled here, since - ! this is not really intended for wide-spread use, but rather for - ! comparison with CM2M / CM2.1 settings. - Kd_sfc(i,j) = max(CS%Kd_min, Kd_sfc(i,j) * (1.0 + & - CS%Kd_tanh_lat_scale * 0.5*tanh((abs(G%geoLatT(i,j)) - 35.0)/5.0) )) - enddo ; enddo - endif + ! set surface diffusivities (CS%bkgnd_mixing_csp%Kd_sfc) + call sfc_bkgnd_mixing(G, CS%bkgnd_mixing_csp) + +! GMM, fix OMP calls below - if (CS%debug) call hchksum(Kd_sfc,"Kd_sfc",G%HI,haloshift=0) !$OMP parallel do default(none) shared(is,ie,js,je,nz,G,GV,CS,h,tv,T_f,S_f,fluxes,dd, & !$OMP Kd,Kd_sfc,epsilon,deg_to_rad,I_2Omega,visc, & !$OMP Kd_int,dt,u,v,Omega2) & @@ -631,55 +362,18 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & !$OMP I_x30,abs_sin,N_2Omega,N02_N2,KT_extra, KS_extra, & !$OMP TKE_to_Kd,maxTKE,dissip,kb) do j=js,je + ! Set up variables related to the stratification. call find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, N2_lay, N2_int, N2_bot) + if (associated(dd%N2_3d)) then do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo endif - ! Set up the background diffusivity. - if (CS%Bryan_Lewis_diffusivity) then - I_trans = 1.0 / CS%Bryan_Lewis_width_trans - atan_fn_sfc = atan(CS%Bryan_Lewis_depth_cent*I_trans) - I_atan_fn = 1.0 / (2.0*atan(1.0) + atan_fn_sfc) - do i=is,ie ; depth(i) = 0.0 ; enddo - do k=1,nz ; do i=is,ie - atan_fn_lay = atan((CS%Bryan_Lewis_depth_cent - & - (depth(i)+0.5*GV%H_to_m*h(i,j,k)))*I_trans) - Kd(i,j,k) = Kd_sfc(i,j) + (CS%Kd_Bryan_Lewis_deep - Kd_sfc(i,j)) * & - (atan_fn_sfc - atan_fn_lay) * I_atan_fn - depth(i) = depth(i) + GV%H_to_m*h(i,j,k) - enddo ; enddo - elseif ((.not.CS%bulkmixedlayer) .and. (CS%Kd /= CS%Kdml)) then - I_Hmix = 1.0 / CS%Hmix - do i=is,ie ; depth(i) = 0.0 ; enddo - do k=1,nz ; do i=is,ie - depth_c = depth(i) + 0.5*GV%H_to_m*h(i,j,k) - - if (depth_c <= CS%Hmix) then ; Kd(i,j,k) = CS%Kdml - elseif (depth_c >= 2.0*CS%Hmix) then ; Kd(i,j,k) = Kd_sfc(i,j) - else - Kd(i,j,k) = ((Kd_sfc(i,j) - CS%Kdml) * I_Hmix) * depth_c + & - (2.0*CS%Kdml - Kd_sfc(i,j)) - endif - - depth(i) = depth(i) + GV%H_to_m*h(i,j,k) - enddo ; enddo - elseif (CS%Henyey_IGW_background_new) then - I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. - do k=1,nz ; do i=is,ie - abs_sin = max(epsilon,abs(sin(G%geoLatT(i,j)*deg_to_rad))) - N_2Omega = max(abs_sin,sqrt(N2_lay(i,k))*I_2Omega) - N02_N2 = (CS%N0_2Omega/N_2Omega)**2 - Kd(i,j,k) = max(CS%Kd_min, Kd_sfc(i,j) * & - ((abs_sin * invcosh(N_2Omega/abs_sin)) * I_x30)*N02_N2) - enddo ; enddo - else - do k=1,nz ; do i=is,ie - Kd(i,j,k) = Kd_sfc(i,j) - enddo ; enddo - endif + ! add background mixing + call calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc%Kv_slow, j, G, GV, CS%bkgnd_mixing_csp) + ! GMM, the following will go into the MOM_CVMix_double_diffusion module if (CS%double_diffusion) then call double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, KT_extra, KS_extra) do K=2,nz ; do i=is,ie @@ -708,18 +402,18 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif ! Add the input turbulent diffusivity. - if (CS%useKappaShear .or. CS%useCVMix) then + if (CS%useKappaShear .or. CS%use_CVMix_shear) then if (present(Kd_int)) then do K=2,nz ; do i=is,ie - Kd_int(i,j,K) = visc%Kd_turb(i,j,K) + 0.5*(Kd(i,j,k-1) + Kd(i,j,k)) + Kd_int(i,j,K) = visc%Kd_shear(i,j,K) + 0.5*(Kd(i,j,k-1) + Kd(i,j,k)) enddo ; enddo do i=is,ie - Kd_int(i,j,1) = visc%Kd_turb(i,j,1) ! This isn't actually used. It could be 0. + Kd_int(i,j,1) = visc%Kd_shear(i,j,1) ! This isn't actually used. It could be 0. Kd_int(i,j,nz+1) = 0.0 enddo endif do k=1,nz ; do i=is,ie - Kd(i,j,k) = Kd(i,j,k) + 0.5*(visc%Kd_turb(i,j,K) + visc%Kd_turb(i,j,K+1)) + Kd(i,j,k) = Kd(i,j,k) + 0.5*(visc%Kd_shear(i,j,K) + visc%Kd_shear(i,j,K+1)) enddo ; enddo else if (present(Kd_int)) then @@ -746,9 +440,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & call add_MLrad_diffusivity(h, fluxes, j, G, GV, CS, Kd, TKE_to_Kd, Kd_int) ! Add the Nikurashin and / or tidal bottom-driven mixing - if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) & - call add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, maxTKE, G, GV, CS, & - dd, N2_lay, Kd, Kd_int) + call calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, maxTKE, G, GV, CS%tm_csp, & + N2_lay, N2_int, Kd, Kd_int, CS%Kd_max) ! This adds the diffusion sustained by the energy extracted from the flow ! by the bottom drag. @@ -799,21 +492,32 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & enddo ! j-loop if (CS%debug) then - call hchksum(Kd,"BBL Kd",G%HI,haloshift=0) - if (CS%useKappaShear) call hchksum(visc%Kd_turb,"Turbulent Kd",G%HI,haloshift=0) + call hchksum(Kd ,"Kd",G%HI,haloshift=0) + + if (CS%useKappaShear) call hchksum(visc%Kd_shear,"Turbulent Kd",G%HI,haloshift=0) + if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, & G%HI, 0, symmetric=.true.) endif + if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) then call uvchksum("BBL bbl_thick_[uv]", visc%bbl_thick_u, & visc%bbl_thick_v, G%HI, 0, symmetric=.true.) endif + if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) then call uvchksum("Ray_[uv]", visc%Ray_u, visc%Ray_v, G%HI, 0, symmetric=.true.) endif + endif + ! send bkgnd_mixing diagnostics to post_data + if (CS%bkgnd_mixing_csp%id_kd_bkgnd > 0) & + call post_data(CS%bkgnd_mixing_csp%id_kd_bkgnd, CS%bkgnd_mixing_csp%kd_bkgnd, CS%bkgnd_mixing_csp%diag) + if (CS%bkgnd_mixing_csp%id_kv_bkgnd > 0) & + call post_data(CS%bkgnd_mixing_csp%id_kv_bkgnd, CS%bkgnd_mixing_csp%kv_bkgnd, CS%bkgnd_mixing_csp%diag) + if (CS%Kd_add > 0.0) then if (present(Kd_int)) then !$OMP parallel do default(none) shared(is,ie,js,je,nz,Kd_int,CS,Kd) @@ -834,54 +538,21 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & T_f, S_f, dd%Kd_user) endif + ! GMM, post diags... if (CS%id_Kd_layer > 0) call post_data(CS%id_Kd_layer, Kd, CS%diag) num_z_diags = 0 - if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) then - if (CS%id_TKE_itidal > 0) call post_data(CS%id_TKE_itidal, dd%TKE_itidal_used, CS%diag) - if (CS%id_TKE_leewave > 0) call post_data(CS%id_TKE_leewave, CS%TKE_Niku, CS%diag) - if (CS%id_Nb > 0) call post_data(CS%id_Nb, CS%Nb, CS%diag) - if (CS%id_N2 > 0) call post_data(CS%id_N2, dd%N2_3d, CS%diag) - if (CS%id_N2_bot > 0) call post_data(CS%id_N2_bot, dd%N2_bot, CS%diag) - if (CS%id_N2_meanz > 0) call post_data(CS%id_N2_meanz,dd%N2_meanz,CS%diag) - - if (CS%id_Fl_itidal > 0) call post_data(CS%id_Fl_itidal, dd%Fl_itidal, CS%diag) - if (CS%id_Kd_itidal > 0) call post_data(CS%id_Kd_itidal, dd%Kd_itidal, CS%diag) - if (CS%id_Kd_Niku > 0) call post_data(CS%id_Kd_Niku, dd%Kd_Niku, CS%diag) - if (CS%id_Kd_lowmode> 0) call post_data(CS%id_Kd_lowmode, dd%Kd_lowmode, CS%diag) - if (CS%id_Fl_lowmode> 0) call post_data(CS%id_Fl_lowmode, dd%Fl_lowmode, CS%diag) - if (CS%id_Kd_user > 0) call post_data(CS%id_Kd_user, dd%Kd_user, CS%diag) - if (CS%id_Kd_Work > 0) call post_data(CS%id_Kd_Work, dd%Kd_Work, CS%diag) - if (CS%id_Kd_Itidal_Work > 0) & - call post_data(CS%id_Kd_Itidal_Work, dd%Kd_Itidal_Work, CS%diag) - if (CS%id_Kd_Niku_Work > 0) call post_data(CS%id_Kd_Niku_Work, dd%Kd_Niku_Work, CS%diag) - if (CS%id_Kd_Lowmode_Work > 0) & - call post_data(CS%id_Kd_Lowmode_Work, dd%Kd_Lowmode_Work, CS%diag) - if (CS%id_maxTKE > 0) call post_data(CS%id_maxTKE, dd%maxTKE, CS%diag) - if (CS%id_TKE_to_Kd > 0) call post_data(CS%id_TKE_to_Kd, dd%TKE_to_Kd, CS%diag) - - if (CS%id_Polzin_decay_scale > 0 ) & - call post_data(CS%id_Polzin_decay_scale, dd%Polzin_decay_scale, CS%diag) - if (CS%id_Polzin_decay_scale_scaled > 0 ) & - call post_data(CS%id_Polzin_decay_scale_scaled, dd%Polzin_decay_scale_scaled, CS%diag) - - if (CS%id_Kd_itidal_z > 0) then - num_z_diags = num_z_diags + 1 - z_ids(num_z_diags) = CS%id_Kd_itidal_z - z_ptrs(num_z_diags)%p => dd%Kd_itidal - endif - if (CS%id_Kd_Niku_z > 0) then - num_z_diags = num_z_diags + 1 - z_ids(num_z_diags) = CS%id_Kd_Niku_z - z_ptrs(num_z_diags)%p => dd%Kd_Niku - endif + call post_tidal_diagnostics(G,GV,h,CS%tm_csp) - if (CS%id_Kd_lowmode_z > 0) then - num_z_diags = num_z_diags + 1 - z_ids(num_z_diags) = CS%id_Kd_lowmode_z - z_ptrs(num_z_diags)%p => dd%Kd_lowmode - endif + if (CS%tm_csp%Int_tide_dissipation .or. CS%tm_csp%Lee_wave_dissipation .or. & + CS%tm_csp%Lowmode_itidal_dissipation) then + + if (CS%id_N2 > 0) call post_data(CS%id_N2, dd%N2_3d, CS%diag) + if (CS%id_Kd_user > 0) call post_data(CS%id_Kd_user, dd%Kd_user, CS%diag) + if (CS%id_Kd_Work > 0) call post_data(CS%id_Kd_Work, dd%Kd_Work, CS%diag) + if (CS%id_maxTKE > 0) call post_data(CS%id_maxTKE, dd%maxTKE, CS%diag) + if (CS%id_TKE_to_Kd > 0) call post_data(CS%id_TKE_to_Kd, dd%TKE_to_Kd, CS%diag) if (CS%id_N2_z > 0) then num_z_diags = num_z_diags + 1 @@ -923,21 +594,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & call calc_Zint_diags(h, z_ptrs, z_ids, num_z_diags, G, GV, CS%diag_to_Z_CSp) if (associated(dd%N2_3d)) deallocate(dd%N2_3d) - if (associated(dd%Kd_itidal)) deallocate(dd%Kd_itidal) - if (associated(dd%Kd_lowmode)) deallocate(dd%Kd_lowmode) - if (associated(dd%Fl_itidal)) deallocate(dd%Fl_itidal) - if (associated(dd%Fl_lowmode)) deallocate(dd%Fl_lowmode) - if (associated(dd%Polzin_decay_scale)) deallocate(dd%Polzin_decay_scale) - if (associated(dd%Polzin_decay_scale_scaled)) deallocate(dd%Polzin_decay_scale_scaled) - if (associated(dd%N2_bot)) deallocate(dd%N2_bot) - if (associated(dd%N2_meanz)) deallocate(dd%N2_meanz) if (associated(dd%Kd_work)) deallocate(dd%Kd_work) if (associated(dd%Kd_user)) deallocate(dd%Kd_user) - if (associated(dd%Kd_Niku)) deallocate(dd%Kd_Niku) - if (associated(dd%Kd_Niku_work)) deallocate(dd%Kd_Niku_work) - if (associated(dd%Kd_Itidal_Work)) deallocate(dd%Kd_Itidal_Work) - if (associated(dd%Kd_Lowmode_Work)) deallocate(dd%Kd_Lowmode_Work) - if (associated(dd%TKE_itidal_used)) deallocate(dd%TKE_itidal_used) if (associated(dd%maxTKE)) deallocate(dd%maxTKE) if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd) if (associated(dd%KT_extra)) deallocate(dd%KT_extra) @@ -1228,8 +886,9 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) do_i(i) = (G%mask2dT(i,j) > 0.5) - if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation) then - h_amp(i) = sqrt(CS%h2(i,j)) ! for computing Nb + if ( (CS%tm_csp%Int_tide_dissipation .or. CS%tm_csp%Lee_wave_dissipation) .and. & + .not. CS%tm_csp%use_CVMix_tidal ) then + h_amp(i) = sqrt(CS%tm_csp%h2(i,j)) ! for computing Nb else h_amp(i) = 0.0 endif @@ -1293,6 +952,8 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & end subroutine find_N2 +! GMM, the following will be moved to a new module + !> This subroutine sets the additional diffusivities of temperature and !! salinity due to double diffusion, using the same functional form as is !! used in MOM4.1, and taken from an NCAR technical note (REF?) that updates @@ -1892,396 +1553,6 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, CS, Kd, TKE_to_Kd, Kd_int) end subroutine add_MLrad_diffusivity - !> This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities. - !! The mechanisms considered are (1) local dissipation of internal waves generated by the - !! barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating - !! low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves. - !! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction, - !! Froude-number-depending breaking, PSI, etc.). -subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, & - dd, N2_lay, Kd, Kd_int ) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) - real, dimension(SZI_(G)), intent(in) :: N2_bot - real, dimension(SZI_(G),SZK_(G)), intent(in) :: N2_lay - integer, intent(in) :: j - real, dimension(SZI_(G),SZK_(G)), intent(in) :: TKE_to_Kd, max_TKE - type(set_diffusivity_CS), pointer :: CS - type(diffusivity_diags), intent(inout) :: dd - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int - - ! This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities. - ! The mechanisms considered are (1) local dissipation of internal waves generated by the - ! barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating - ! low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves. - ! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction, - ! Froude-number-depending breaking, PSI, etc.). - - real, dimension(SZI_(G)) :: & - htot, & ! total thickness above or below a layer, or the - ! integrated thickness in the BBL (meter) - htot_WKB, & ! distance from top to bottom (meter) WKB scaled - TKE_itidal_bot, & ! internal tide TKE at ocean bottom (m3/s3) - TKE_Niku_bot, & ! lee-wave TKE at ocean bottom (m3/s3) - TKE_lowmode_bot, & ! internal tide TKE at ocean bottom lost from all remote low modes (m3/s3) (BDM) - Inv_int, & ! inverse of TKE decay for int tide over the depth of the ocean (nondim) - Inv_int_lee, & ! inverse of TKE decay for lee waves over the depth of the ocean (nondim) - Inv_int_low, & ! inverse of TKE decay for low modes over the depth of the ocean (nondim) (BDM) - z0_Polzin, & ! TKE decay scale in Polzin formulation (meter) - z0_Polzin_scaled, & ! TKE decay scale in Polzin formulation (meter) - ! multiplied by N2_bot/N2_meanz to be coherent with the WKB scaled z - ! z*=int(N2/N2_bot) * N2_bot/N2_meanz = int(N2/N2_meanz) - ! z0_Polzin_scaled = z0_Polzin * N2_bot/N2_meanz - N2_meanz, & ! vertically averaged squared buoyancy frequency (1/s2) for WKB scaling - TKE_itidal_rem, & ! remaining internal tide TKE (from barotropic source) - TKE_Niku_rem, & ! remaining lee-wave TKE - TKE_lowmode_rem, & ! remaining internal tide TKE (from propagating low mode source) (BDM) - TKE_frac_top, & ! fraction of bottom TKE that should appear at top of a layer (nondim) - TKE_frac_top_lee, & ! fraction of bottom TKE that should appear at top of a layer (nondim) - TKE_frac_top_lowmode, & - ! fraction of bottom TKE that should appear at top of a layer (nondim) (BDM) - z_from_bot, & ! distance from bottom (meter) - z_from_bot_WKB ! distance from bottom (meter), WKB scaled - - real :: I_rho0 ! 1 / RHO0, (m3/kg) - real :: Kd_add ! diffusivity to add in a layer (m2/sec) - real :: TKE_itide_lay ! internal tide TKE imparted to a layer (from barotropic) (m3/s3) - real :: TKE_Niku_lay ! lee-wave TKE imparted to a layer (m3/s3) - real :: TKE_lowmode_lay ! internal tide TKE imparted to a layer (from low mode) (m3/s3) (BDM) - real :: frac_used ! fraction of TKE that can be used in a layer (nondim) - real :: Izeta ! inverse of TKE decay scale (1/meter) - real :: Izeta_lee ! inverse of TKE decay scale for lee waves (1/meter) - real :: z0_psl ! temporary variable with units of meter - real :: TKE_lowmode_tot ! TKE from all low modes (W/m2) (BDM) - - - logical :: use_Polzin, use_Simmons - integer :: i, k, is, ie, nz - integer :: a, fr, m - is = G%isc ; ie = G%iec ; nz = G%ke - - if (.not.(CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation)) return - - do i=is,ie ; htot(i) = 0.0 ; Inv_int(i) = 0.0 ; Inv_int_lee(i) = 0.0 ; Inv_int_low(i) = 0.0 ;enddo - do k=1,nz ; do i=is,ie - htot(i) = htot(i) + GV%H_to_m*h(i,j,k) - enddo ; enddo - - I_Rho0 = 1.0/GV%Rho0 - - use_Polzin = ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09)) .or. & - (CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09)) .or. & - (CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == POLZIN_09))) - use_Simmons = ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == STLAURENT_02)) .or. & - (CS%lee_wave_dissipation .and. (CS%lee_wave_profile == STLAURENT_02)) .or. & - (CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == STLAURENT_02))) - - ! Calculate parameters for vertical structure of dissipation - ! Simmons: - if ( use_Simmons ) then - Izeta = 1.0 / max(CS%Int_tide_decay_scale, GV%H_subroundoff*GV%H_to_m) - Izeta_lee = 1.0 / max(CS%Int_tide_decay_scale*CS%Decay_scale_factor_lee, & - GV%H_subroundoff*GV%H_to_m) - do i=is,ie - CS%Nb(i,j) = sqrt(N2_bot(i)) - if (associated(dd%N2_bot)) dd%N2_bot(i,j) = N2_bot(i) - if ( CS%Int_tide_dissipation ) then - if (Izeta*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule. - Inv_int(i) = 1.0 / (1.0 - exp(-Izeta*htot(i))) - endif - endif - if ( CS%Lee_wave_dissipation ) then - if (Izeta_lee*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule. - Inv_int_lee(i) = 1.0 / (1.0 - exp(-Izeta_lee*htot(i))) - endif - endif - if ( CS%Lowmode_itidal_dissipation) then - if (Izeta*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule. - Inv_int_low(i) = 1.0 / (1.0 - exp(-Izeta*htot(i))) - endif - endif - z_from_bot(i) = GV%H_to_m*h(i,j,nz) - enddo - endif ! Simmons - - ! Polzin: - if ( use_Polzin ) then - ! WKB scaling of the vertical coordinate - do i=is,ie ; N2_meanz(i)=0.0 ; enddo - do k=1,nz ; do i=is,ie - N2_meanz(i) = N2_meanz(i) + N2_lay(i,k)*GV%H_to_m*h(i,j,k) - enddo ; enddo - do i=is,ie - N2_meanz(i) = N2_meanz(i) / (htot(i) + GV%H_subroundoff*GV%H_to_m) - if (associated(dd%N2_meanz)) dd%N2_meanz(i,j) = N2_meanz(i) - enddo - - ! WKB scaled z*(z=H) z* at the surface using the modified Polzin WKB scaling - do i=is,ie ; htot_WKB(i) = htot(i) ; enddo -! do i=is,ie ; htot_WKB(i) = 0.0 ; enddo -! do k=1,nz ; do i=is,ie -! htot_WKB(i) = htot_WKB(i) + GV%H_to_m*h(i,j,k)*N2_lay(i,k) / N2_meanz(i) -! enddo ; enddo - ! htot_WKB(i) = htot(i) ! Nearly equivalent and simpler - - do i=is,ie - CS%Nb(i,j) = sqrt(N2_bot(i)) - if ((CS%tideamp(i,j) > 0.0) .and. & - (CS%kappa_itides**2 * CS%h2(i,j) * CS%Nb(i,j)**3 > 1.0e-14) ) then - z0_polzin(i) = CS%Polzin_decay_scale_factor * CS%Nu_Polzin * & - CS%Nbotref_Polzin**2 * CS%tideamp(i,j) / & - ( CS%kappa_itides**2 * CS%h2(i,j) * CS%Nb(i,j)**3 ) - if (z0_polzin(i) < CS%Polzin_min_decay_scale) & - z0_polzin(i) = CS%Polzin_min_decay_scale - if (N2_meanz(i) > 1.0e-14 ) then - z0_polzin_scaled(i) = z0_polzin(i)*CS%Nb(i,j)**2 / N2_meanz(i) - else - z0_polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * htot(i) - endif - if (z0_polzin_scaled(i) > (CS%Polzin_decay_scale_max_factor * htot(i)) ) & - z0_polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * htot(i) - else - z0_polzin(i) = CS%Polzin_decay_scale_max_factor * htot(i) - z0_polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * htot(i) - endif - - if (associated(dd%Polzin_decay_scale)) & - dd%Polzin_decay_scale(i,j) = z0_polzin(i) - if (associated(dd%Polzin_decay_scale_scaled)) & - dd%Polzin_decay_scale_scaled(i,j) = z0_polzin_scaled(i) - if (associated(dd%N2_bot)) dd%N2_bot(i,j) = CS%Nb(i,j)*CS%Nb(i,j) - - if ( CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09) ) then - ! For the Polzin formulation, this if loop prevents the vertical - ! flux of energy dissipation from having NaN values - if (htot_WKB(i) > 1.0e-14) then - Inv_int(i) = ( z0_polzin_scaled(i) / htot_WKB(i) ) + 1 - endif - endif - if ( CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09) ) then - ! For the Polzin formulation, this if loop prevents the vertical - ! flux of energy dissipation from having NaN values - if (htot_WKB(i) > 1.0e-14) then - Inv_int_lee(i) = ( z0_polzin_scaled(i)*CS%Decay_scale_factor_lee / htot_WKB(i) ) + 1 - endif - endif - if ( CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == POLZIN_09) ) then - ! For the Polzin formulation, this if loop prevents the vertical - ! flux of energy dissipation from having NaN values - if (htot_WKB(i) > 1.0e-14) then - Inv_int_low(i) = ( z0_polzin_scaled(i) / htot_WKB(i) ) + 1 - endif - endif - - z_from_bot(i) = GV%H_to_m*h(i,j,nz) - ! Use the new formulation for WKB scaling. N2 is referenced to its - ! vertical mean. - if (N2_meanz(i) > 1.0e-14 ) then - z_from_bot_WKB(i) = GV%H_to_m*h(i,j,nz)*N2_lay(i,nz) / N2_meanz(i) - else ; z_from_bot_WKB(i) = 0 ; endif - enddo - endif ! Polzin - - ! Calculate/get dissipation values at bottom - ! Both Polzin and Simmons: - do i=is,ie - ! Dissipation of locally trapped internal tide (non-propagating high modes) - TKE_itidal_bot(i) = min(CS%TKE_itidal(i,j)*CS%Nb(i,j),CS%TKE_itide_max) - if (associated(dd%TKE_itidal_used)) & - dd%TKE_itidal_used(i,j) = TKE_itidal_bot(i) - TKE_itidal_bot(i) = (I_rho0 * CS%Mu_itides * CS%Gamma_itides) * TKE_itidal_bot(i) - ! Dissipation of locally trapped lee waves - TKE_Niku_bot(i) = 0.0 - if (CS%Lee_wave_dissipation) then - TKE_Niku_bot(i) = (I_rho0 * CS%Mu_itides * CS%Gamma_lee) * CS%TKE_Niku(i,j) - endif - ! Dissipation of propagating internal tide (baroclinic low modes; rays) (BDM) - TKE_lowmode_tot = 0.0 - TKE_lowmode_bot(i) = 0.0 - if (CS%Lowmode_itidal_dissipation) then - ! get loss rate due to wave drag on low modes (already multiplied by q) - call get_lowmode_loss(i,j,G,CS%int_tide_CSp,"WaveDrag",TKE_lowmode_tot) - TKE_lowmode_bot(i) = CS%Mu_itides * I_rho0 * TKE_lowmode_tot - endif - ! Vertical energy flux at bottom - TKE_itidal_rem(i) = Inv_int(i) * TKE_itidal_bot(i) - TKE_Niku_rem(i) = Inv_int_lee(i) * TKE_Niku_bot(i) - TKE_lowmode_rem(i) = Inv_int_low(i) * TKE_lowmode_bot(i) - - if (associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,nz) = TKE_itidal_rem(i) !why is this here? BDM - enddo - - ! Estimate the work that would be done by mixing in each layer. - ! Simmons: - if ( use_Simmons ) then - do k=nz-1,2,-1 ; do i=is,ie - if (max_TKE(i,k) <= 0.0) cycle - z_from_bot(i) = z_from_bot(i) + GV%H_to_m*h(i,j,k) - - ! Fraction of bottom flux predicted to reach top of this layer - TKE_frac_top(i) = Inv_int(i) * exp(-Izeta * z_from_bot(i)) - TKE_frac_top_lee(i) = Inv_int_lee(i) * exp(-Izeta_lee * z_from_bot(i)) - TKE_frac_top_lowmode(i) = Inv_int_low(i) * exp(-Izeta * z_from_bot(i)) - - ! Actual influx at bottom of layer minus predicted outflux at top of layer to give - ! predicted power expended - TKE_itide_lay = TKE_itidal_rem(i) - TKE_itidal_bot(i) * TKE_frac_top(i) - TKE_Niku_lay = TKE_Niku_rem(i) - TKE_Niku_bot(i) * TKE_frac_top_lee(i) - TKE_lowmode_lay = TKE_lowmode_rem(i) - TKE_lowmode_bot(i)* TKE_frac_top_lowmode(i) - - ! Actual power expended may be less than predicted if stratification is weak; adjust - if (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay > max_TKE(i,k)) then - frac_used = max_TKE(i,k) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) - TKE_itide_lay = frac_used * TKE_itide_lay - TKE_Niku_lay = frac_used * TKE_Niku_lay - TKE_lowmode_lay = frac_used * TKE_lowmode_lay - endif - - ! Calculate vertical flux available to bottom of layer above - TKE_itidal_rem(i) = TKE_itidal_rem(i) - TKE_itide_lay - TKE_Niku_rem(i) = TKE_Niku_rem(i) - TKE_Niku_lay - TKE_lowmode_rem(i) = TKE_lowmode_rem(i) - TKE_lowmode_lay - - ! Convert power to diffusivity - Kd_add = TKE_to_Kd(i,k) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) - - if (CS%Kd_max >= 0.0) Kd_add = min(Kd_add, CS%Kd_max) - Kd(i,j,k) = Kd(i,j,k) + Kd_add - - if (present(Kd_int)) then - Kd_int(i,j,K) = Kd_int(i,j,K) + 0.5*Kd_add - Kd_int(i,j,K+1) = Kd_int(i,j,K+1) + 0.5*Kd_add - endif - - ! diagnostics - if (associated(dd%Kd_itidal)) then - ! If at layers, dd%Kd_itidal is just TKE_to_Kd(i,k) * TKE_itide_lay - ! The following sets the interface diagnostics. - Kd_add = TKE_to_Kd(i,k) * TKE_itide_lay - if (CS%Kd_max >= 0.0) Kd_add = min(Kd_add, CS%Kd_max) - if (k>1) dd%Kd_itidal(i,j,K) = dd%Kd_itidal(i,j,K) + 0.5*Kd_add - if (k= 0.0) Kd_add = min(Kd_add, CS%Kd_max) - if (k>1) dd%Kd_Niku(i,j,K) = dd%Kd_Niku(i,j,K) + 0.5*Kd_add - if (k= 0.0) Kd_add = min(Kd_add, CS%Kd_max) - if (k>1) dd%Kd_lowmode(i,j,K) = dd%Kd_lowmode(i,j,K) + 0.5*Kd_add - if (k 1.0e-14 ) then - z_from_bot_WKB(i) = z_from_bot_WKB(i) + GV%H_to_m*h(i,j,k)*N2_lay(i,k)/N2_meanz(i) - else ; z_from_bot_WKB(i) = 0 ; endif - - ! Fraction of bottom flux predicted to reach top of this layer - TKE_frac_top(i) = ( Inv_int(i) * z0_polzin_scaled(i) ) / & - ( z0_polzin_scaled(i) + z_from_bot_WKB(i) ) - z0_psl = z0_polzin_scaled(i)*CS%Decay_scale_factor_lee - TKE_frac_top_lee(i) = (Inv_int_lee(i) * z0_psl) / (z0_psl + z_from_bot_WKB(i)) - TKE_frac_top_lowmode(i) = ( Inv_int_low(i) * z0_polzin_scaled(i) ) / & - ( z0_polzin_scaled(i) + z_from_bot_WKB(i) ) - - ! Actual influx at bottom of layer minus predicted outflux at top of layer to give - ! predicted power expended - TKE_itide_lay = TKE_itidal_rem(i) - TKE_itidal_bot(i) *TKE_frac_top(i) - TKE_Niku_lay = TKE_Niku_rem(i) - TKE_Niku_bot(i) * TKE_frac_top_lee(i) - TKE_lowmode_lay = TKE_lowmode_rem(i) - TKE_lowmode_bot(i)*TKE_frac_top_lowmode(i) - - ! Actual power expended may be less than predicted if stratification is weak; adjust - if (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay > max_TKE(i,k)) then - frac_used = max_TKE(i,k) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) - TKE_itide_lay = frac_used * TKE_itide_lay - TKE_Niku_lay = frac_used * TKE_Niku_lay - TKE_lowmode_lay = frac_used * TKE_lowmode_lay - endif - - ! Calculate vertical flux available to bottom of layer above - TKE_itidal_rem(i) = TKE_itidal_rem(i) - TKE_itide_lay - TKE_Niku_rem(i) = TKE_Niku_rem(i) - TKE_Niku_lay - TKE_lowmode_rem(i) = TKE_lowmode_rem(i) - TKE_lowmode_lay - - ! Convert power to diffusivity - Kd_add = TKE_to_Kd(i,k) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) - - if (CS%Kd_max >= 0.0) Kd_add = min(Kd_add, CS%Kd_max) - Kd(i,j,k) = Kd(i,j,k) + Kd_add - - if (present(Kd_int)) then - Kd_int(i,j,K) = Kd_int(i,j,K) + 0.5*Kd_add - Kd_int(i,j,K+1) = Kd_int(i,j,K+1) + 0.5*Kd_add - endif - - ! diagnostics - if (associated(dd%Kd_itidal)) then - ! If at layers, this is just dd%Kd_itidal(i,j,K) = TKE_to_Kd(i,k) * TKE_itide_lay - ! The following sets the interface diagnostics. - Kd_add = TKE_to_Kd(i,k) * TKE_itide_lay - if (CS%Kd_max >= 0.0) Kd_add = min(Kd_add, CS%Kd_max) - if (k>1) dd%Kd_itidal(i,j,K) = dd%Kd_itidal(i,j,K) + 0.5*Kd_add - if (k= 0.0) Kd_add = min(Kd_add, CS%Kd_max) - if (k>1) dd%Kd_Niku(i,j,K) = dd%Kd_Niku(i,j,K) + 0.5*Kd_add - if (k= 0.0) Kd_add = min(Kd_add, CS%Kd_max) - if (k>1) dd%Kd_lowmode(i,j,K) = dd%Kd_lowmode(i,j,K) + 0.5*Kd_add - if (k This subroutine calculates several properties related to bottom !! boundary layer turbulence. subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, CS) @@ -2533,7 +1804,8 @@ subroutine set_density_ratios(h, tv, kb, G, GV, CS, j, ds_dsp1, rho_0) end subroutine set_density_ratios -subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp, int_tide_CSp) +subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp, int_tide_CSp, & + tm_CSp) type(time_type), intent(in) :: Time type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -2546,26 +1818,18 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp !! structure. type(int_tide_CS), pointer :: int_tide_CSp !< pointer to the internal tides control !! structure (BDM) + type(tidal_mixing_cs), pointer :: tm_csp !< pointer to tidal mixing control + !! structure -! Arguments: -! (in) Time - current model time -! (in) G - ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) param_file - structure indicating open file to parse for params -! (in) diag - structure used to regulate diagnostic output -! (in/out) CS - pointer set to point to the module control structure -! (in) diag_to_Z_CSp - pointer to the Z-diagnostics control structure -! (in) int_tide_CSp - pointer to the internal tides control structure (BDM) - - real :: decay_length, utide, zbot, hamp + ! local variables + real :: decay_length type(vardesc) :: vd - logical :: read_tideamp, ML_use_omega + logical :: ML_use_omega + ! This include declares and sets the variable "version". #include "version_variable.h" + character(len=40) :: mdl = "MOM_set_diffusivity" ! This module's name. - character(len=20) :: tmpstr - character(len=200) :: filename, tideamp_file, h2_file, Niku_TKE_input_file - real :: Niku_scale ! local variable for scaling the Nikurashin TKE flux data real :: omega_frac_dflt integer :: i, j, is, ie, js, je integer :: isd, ied, jsd, jed @@ -2582,6 +1846,7 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp CS%diag => diag if (associated(int_tide_CSp)) CS%int_tide_CSp => int_tide_CSp + if (associated(tm_csp)) CS%tm_csp => tm_csp if (associated(diag_to_Z_CSp)) CS%diag_to_Z_CSp => diag_to_Z_CSp ! These default values always need to be set. @@ -2706,79 +1971,8 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp "for an isopycnal layer-formulation.", & default=.false.) - call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", & - CS%Bryan_Lewis_diffusivity, & - "If true, use a Bryan & Lewis (JGR 1979) like tanh \n"//& - "profile of background diapycnal diffusivity with depth.", & - default=.false.) - if (CS%Bryan_Lewis_diffusivity) then - call get_param(param_file, mdl, "KD_BRYAN_LEWIS_DEEP", & - CS%Kd_Bryan_Lewis_deep, & - "The abyssal value of a Bryan-Lewis diffusivity profile. \n"//& - "KD_BRYAN_LEWIS_DEEP is only used if \n"//& - "BRYAN_LEWIS_DIFFUSIVITY is true.", units="m2 s-1", & - fail_if_missing=.true.) - call get_param(param_file, mdl, "KD_BRYAN_LEWIS_SURFACE", & - CS%Kd_Bryan_Lewis_surface, & - "The surface value of a Bryan-Lewis diffusivity profile. \n"//& - "KD_BRYAN_LEWIS_SURFACE is only used if \n"//& - "BRYAN_LEWIS_DIFFUSIVITY is true.", units="m2 s-1", & - fail_if_missing=.true.) - call get_param(param_file, mdl, "BRYAN_LEWIS_DEPTH_CENT", & - CS%Bryan_Lewis_depth_cent, & - "The depth about which the transition in the Bryan-Lewis \n"//& - "profile is centered. BRYAN_LEWIS_DEPTH_CENT is only \n"//& - "used if BRYAN_LEWIS_DIFFUSIVITY is true.", units="m", & - fail_if_missing=.true.) - call get_param(param_file, mdl, "BRYAN_LEWIS_WIDTH_TRANS", & - CS%Bryan_Lewis_width_trans, & - "The width of the transition in the Bryan-Lewis \n"//& - "profile. BRYAN_LEWIS_WIDTH_TRANS is only \n"//& - "used if BRYAN_LEWIS_DIFFUSIVITY is true.", units="m", & - fail_if_missing=.true.) - endif - - call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND", & - CS%Henyey_IGW_background, & - "If true, use a latitude-dependent scaling for the near \n"//& - "surface background diffusivity, as described in \n"//& - "Harrison & Hallberg, JPO 2008.", default=.false.) - call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND_NEW", & - CS%Henyey_IGW_background_new, & - "If true, use a better latitude-dependent scaling for the\n"//& - "background diffusivity, as described in \n"//& - "Harrison & Hallberg, JPO 2008.", default=.false.) - if (CS%Henyey_IGW_background .and. CS%Henyey_IGW_background_new) call MOM_error(FATAL, & - "set_diffusivity_init: HENYEY_IGW_BACKGROUND and HENYEY_IGW_BACKGROUND_NEW "// & - "are mutually exclusive. Set only one or none.") - if (CS%Henyey_IGW_background) & - call get_param(param_file, mdl, "HENYEY_N0_2OMEGA", CS%N0_2Omega, & - "The ratio of the typical Buoyancy frequency to twice \n"//& - "the Earth's rotation period, used with the Henyey \n"//& - "scaling from the mixing.", units="nondim", default=20.0) - call get_param(param_file, mdl, "N2_FLOOR_IOMEGA2", CS%N2_FLOOR_IOMEGA2, & - "The floor applied to N2(k) scaled by Omega^2:\n"//& - "\tIf =0., N2(k) is simply positive definite.\n"//& - "\tIf =1., N2(k) > Omega^2 everywhere.", units="nondim", & - default=1.0) - - call get_param(param_file, mdl, "KD_TANH_LAT_FN", & - CS%Kd_tanh_lat_fn, & - "If true, use a tanh dependence of Kd_sfc on latitude, \n"//& - "like CM2.1/CM2M. There is no physical justification \n"//& - "for this form, and it can not be used with \n"//& - "HENYEY_IGW_BACKGROUND.", default=.false.) - if (CS%Kd_tanh_lat_fn) & - call get_param(param_file, mdl, "KD_TANH_LAT_SCALE", & - CS%Kd_tanh_lat_scale, & - "A nondimensional scaling for the range ofdiffusivities \n"//& - "with KD_TANH_LAT_FN. Valid values are in the range of \n"//& - "-2 to 2; 0.4 reproduces CM2M.", units="nondim", default=0.0) - - call get_param(param_file, mdl, "KV", CS%Kv, & - "The background kinematic viscosity in the interior. \n"//& - "The molecular value, ~1e-6 m2 s-1, may be used.", & - units="m2 s-1", fail_if_missing=.true.) + ! set params releted to the background mixing + call bkgnd_mixing_init(Time, G, GV, param_file, CS%diag, CS%bkgnd_mixing_csp) call get_param(param_file, mdl, "KD", CS%Kd, & "The background diapycnal diffusivity of density in the \n"//& @@ -2824,90 +2018,6 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) - call get_param(param_file, mdl, "INT_TIDE_DISSIPATION", CS%Int_tide_dissipation, & - "If true, use an internal tidal dissipation scheme to \n"//& - "drive diapycnal mixing, along the lines of St. Laurent \n"//& - "et al. (2002) and Simmons et al. (2004).", default=.false.) - if (CS%Int_tide_dissipation) then - call get_param(param_file, mdl, "INT_TIDE_PROFILE", tmpstr, & - "INT_TIDE_PROFILE selects the vertical profile of energy \n"//& - "dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//& - "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//& - "\t decay profile.\n"//& - "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//& - "\t decay profile.", & - default=STLAURENT_PROFILE_STRING) - tmpstr = uppercase(tmpstr) - select case (tmpstr) - case (STLAURENT_PROFILE_STRING) ; CS%int_tide_profile = STLAURENT_02 - case (POLZIN_PROFILE_STRING) ; CS%int_tide_profile = POLZIN_09 - case default - call MOM_error(FATAL, "set_diffusivity_init: Unrecognized setting "// & - "#define INT_TIDE_PROFILE "//trim(tmpstr)//" found in input file.") - end select - endif - - call get_param(param_file, mdl, "LEE_WAVE_DISSIPATION", CS%Lee_wave_dissipation, & - "If true, use an lee wave driven dissipation scheme to \n"//& - "drive diapycnal mixing, along the lines of Nikurashin \n"//& - "(2010) and using the St. Laurent et al. (2002) \n"//& - "and Simmons et al. (2004) vertical profile", default=.false.) - if (CS%lee_wave_dissipation) then - call get_param(param_file, mdl, "LEE_WAVE_PROFILE", tmpstr, & - "LEE_WAVE_PROFILE selects the vertical profile of energy \n"//& - "dissipation with LEE_WAVE_DISSIPATION. Valid values are:\n"//& - "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//& - "\t decay profile.\n"//& - "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//& - "\t decay profile.", & - default=STLAURENT_PROFILE_STRING) - tmpstr = uppercase(tmpstr) - select case (tmpstr) - case (STLAURENT_PROFILE_STRING) ; CS%lee_wave_profile = STLAURENT_02 - case (POLZIN_PROFILE_STRING) ; CS%lee_wave_profile = POLZIN_09 - case default - call MOM_error(FATAL, "set_diffusivity_init: Unrecognized setting "// & - "#define LEE_WAVE_PROFILE "//trim(tmpstr)//" found in input file.") - end select - endif - - call get_param(param_file, mdl, "INT_TIDE_LOWMODE_DISSIPATION", CS%Lowmode_itidal_dissipation, & - "If true, consider mixing due to breaking low modes that \n"//& - "have been remotely generated; as with itidal drag on the \n"//& - "barotropic tide, use an internal tidal dissipation scheme to \n"//& - "drive diapycnal mixing, along the lines of St. Laurent \n"//& - "et al. (2002) and Simmons et al. (2004).", default=.false.) - - if ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09)) .or. & - (CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09))) then - call get_param(param_file, mdl, "NU_POLZIN", CS%Nu_Polzin, & - "When the Polzin decay profile is used, this is a \n"//& - "non-dimensional constant in the expression for the \n"//& - "vertical scale of decay for the tidal energy dissipation.", & - units="nondim", default=0.0697) - call get_param(param_file, mdl, "NBOTREF_POLZIN", CS%Nbotref_Polzin, & - "When the Polzin decay profile is used, this is the \n"//& - "Rreference value of the buoyancy frequency at the ocean \n"//& - "bottom in the Polzin formulation for the vertical \n"//& - "scale of decay for the tidal energy dissipation.", & - units="s-1", default=9.61e-4) - call get_param(param_file, mdl, "POLZIN_DECAY_SCALE_FACTOR", & - CS%Polzin_decay_scale_factor, & - "When the Polzin decay profile is used, this is a \n"//& - "scale factor for the vertical scale of decay of the tidal \n"//& - "energy dissipation.", default=1.0, units="nondim") - call get_param(param_file, mdl, "POLZIN_SCALE_MAX_FACTOR", & - CS%Polzin_decay_scale_max_factor, & - "When the Polzin decay profile is used, this is a factor \n"//& - "to limit the vertical scale of decay of the tidal \n"//& - "energy dissipation to POLZIN_DECAY_SCALE_MAX_FACTOR \n"//& - "times the depth of the ocean.", units="nondim", default=1.0) - call get_param(param_file, mdl, "POLZIN_MIN_DECAY_SCALE", CS%Polzin_min_decay_scale, & - "When the Polzin decay profile is used, this is the \n"//& - "minimum vertical decay scale for the vertical profile\n"//& - "of internal tide dissipation with the Polzin (2009) formulation", & - units="m", default=0.0) - endif call get_param(param_file, mdl, "USER_CHANGE_DIFFUSIVITY", CS%user_change_diff, & "If true, call user-defined code to change the diffusivity.", & default=.false.) @@ -2935,171 +2045,19 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp if (CS%FluxRi_max > 0.0) & CS%dissip_N2 = CS%dissip_Kd_min * GV%Rho0 / CS%FluxRi_max - if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation) then - call get_param(param_file, mdl, "INT_TIDE_DECAY_SCALE", CS%Int_tide_decay_scale, & - "The decay scale away from the bottom for tidal TKE with \n"//& - "the new coding when INT_TIDE_DISSIPATION is used.", & - units="m", default=0.0) - call get_param(param_file, mdl, "MU_ITIDES", CS%Mu_itides, & - "A dimensionless turbulent mixing efficiency used with \n"//& - "INT_TIDE_DISSIPATION, often 0.2.", units="nondim", default=0.2) - call get_param(param_file, mdl, "GAMMA_ITIDES", CS%Gamma_itides, & - "The fraction of the internal tidal energy that is \n"//& - "dissipated locally with INT_TIDE_DISSIPATION. \n"//& - "THIS NAME COULD BE BETTER.", & - units="nondim", default=0.3333) - call get_param(param_file, mdl, "MIN_ZBOT_ITIDES", CS%min_zbot_itides, & - "Turn off internal tidal dissipation when the total \n"//& - "ocean depth is less than this value.", units="m", default=0.0) - - call safe_alloc_ptr(CS%Nb,isd,ied,jsd,jed) - call safe_alloc_ptr(CS%h2,isd,ied,jsd,jed) - call safe_alloc_ptr(CS%TKE_itidal,isd,ied,jsd,jed) - call safe_alloc_ptr(CS%mask_itidal,isd,ied,jsd,jed) ; CS%mask_itidal(:,:) = 1.0 - - call get_param(param_file, mdl, "KAPPA_ITIDES", CS%kappa_itides, & - "A topographic wavenumber used with INT_TIDE_DISSIPATION. \n"//& - "The default is 2pi/10 km, as in St.Laurent et al. 2002.", & - units="m-1", default=8.e-4*atan(1.0)) - - call get_param(param_file, mdl, "UTIDE", CS%utide, & - "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & - units="m s-1", default=0.0) - call safe_alloc_ptr(CS%tideamp,is,ie,js,je) ; CS%tideamp(:,:) = CS%utide - - call get_param(param_file, mdl, "KAPPA_H2_FACTOR", CS%kappa_h2_factor, & - "A scaling factor for the roughness amplitude with n"//& - "INT_TIDE_DISSIPATION.", units="nondim", default=1.0) - call get_param(param_file, mdl, "TKE_ITIDE_MAX", CS%TKE_itide_max, & - "The maximum internal tide energy source availble to mix \n"//& - "above the bottom boundary layer with INT_TIDE_DISSIPATION.", & - units="W m-2", default=1.0e3) - - call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, & - "If true, read a file (given by TIDEAMP_FILE) containing \n"//& - "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.) - if (read_tideamp) then - call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, & - "The path to the file containing the spatially varying \n"//& - "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc") - filename = trim(CS%inputdir) // trim(tideamp_file) - call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename) - call MOM_read_data(filename, 'tideamp', CS%tideamp, G%domain, timelevel=1) - endif - - call get_param(param_file, mdl, "H2_FILE", h2_file, & - "The path to the file containing the sub-grid-scale \n"//& - "topographic roughness amplitude with INT_TIDE_DISSIPATION.", & - fail_if_missing=.true.) - filename = trim(CS%inputdir) // trim(h2_file) - call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename) - call MOM_read_data(filename, 'h2', CS%h2, G%domain, timelevel=1) - - do j=js,je ; do i=is,ie - if (G%bathyT(i,j) < CS%min_zbot_itides) CS%mask_itidal(i,j) = 0.0 - CS%tideamp(i,j) = CS%tideamp(i,j) * CS%mask_itidal(i,j) * G%mask2dT(i,j) - - ! Restrict rms topo to 10 percent of column depth. - zbot = G%bathyT(i,j) - hamp = sqrt(CS%h2(i,j)) - hamp = min(0.1*zbot,hamp) - CS%h2(i,j) = hamp*hamp - - utide = CS%tideamp(i,j) - ! Compute the fixed part of internal tidal forcing; units are [kg s-2] here. - CS%TKE_itidal(i,j) = 0.5*CS%kappa_h2_factor*GV%Rho0*& - CS%kappa_itides*CS%h2(i,j)*utide*utide - enddo; enddo - - endif - CS%id_Kd_layer = register_diag_field('ocean_model', 'Kd_layer', diag%axesTL, Time, & 'Diapycnal diffusivity of layers (as set)', 'm2 s-1') - if (CS%Lee_wave_dissipation) then - - call get_param(param_file, mdl, "NIKURASHIN_TKE_INPUT_FILE",Niku_TKE_input_file, & - "The path to the file containing the TKE input from lee \n"//& - "wave driven mixing. Used with LEE_WAVE_DISSIPATION.", & - fail_if_missing=.true.) - call get_param(param_file, mdl, "NIKURASHIN_SCALE",Niku_scale, & - "A non-dimensional factor by which to scale the lee-wave \n"//& - "driven TKE input. Used with LEE_WAVE_DISSIPATION.", & - units="nondim", default=1.0) - - filename = trim(CS%inputdir) // trim(Niku_TKE_input_file) - call log_param(param_file, mdl, "INPUTDIR/NIKURASHIN_TKE_INPUT_FILE", & - filename) - call safe_alloc_ptr(CS%TKE_Niku,is,ie,js,je); CS%TKE_Niku(:,:) = 0.0 - call MOM_read_data(filename, 'TKE_input', CS%TKE_Niku, G%domain, timelevel=1 ) ! ??? timelevel -aja - CS%TKE_Niku(:,:) = Niku_scale * CS%TKE_Niku(:,:) - - call get_param(param_file, mdl, "GAMMA_NIKURASHIN",CS%Gamma_lee, & - "The fraction of the lee wave energy that is dissipated \n"//& - "locally with LEE_WAVE_DISSIPATION.", units="nondim", & - default=0.3333) - call get_param(param_file, mdl, "DECAY_SCALE_FACTOR_LEE",CS%Decay_scale_factor_lee, & - "Scaling for the vertical decay scaleof the local \n"//& - "dissipation of lee waves dissipation.", units="nondim", & - default=1.0) - - CS%id_TKE_leewave = register_diag_field('ocean_model','TKE_leewave',diag%axesT1,Time, & - 'Lee wave Driven Turbulent Kinetic Energy', 'W m-2') - CS%id_Kd_Niku = register_diag_field('ocean_model','Kd_Nikurashin',diag%axesTi,Time, & - 'Lee Wave Driven Diffusivity', 'm2 s-1') - else - CS%Decay_scale_factor_lee = -9.e99 ! This should never be used if CS%Lee_wave_dissipation = False - endif - if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. & - CS%Lowmode_itidal_dissipation) then + if (CS%tm_csp%Int_tide_dissipation .or. CS%tm_csp%Lee_wave_dissipation .or. & + CS%tm_csp%Lowmode_itidal_dissipation) then - CS%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal',diag%axesT1,Time, & - 'Internal Tide Driven Turbulent Kinetic Energy', 'W m-2') + CS%id_Kd_Work = register_diag_field('ocean_model','Kd_Work',diag%axesTL,Time, & + 'Work done by Diapycnal Mixing', 'W m-2') CS%id_maxTKE = register_diag_field('ocean_model','maxTKE',diag%axesTL,Time, & 'Maximum layer TKE', 'm3 s-3') CS%id_TKE_to_Kd = register_diag_field('ocean_model','TKE_to_Kd',diag%axesTL,Time, & 'Convert TKE to Kd', 's2 m') - - CS%id_Nb = register_diag_field('ocean_model','Nb',diag%axesT1,Time, & - 'Bottom Buoyancy Frequency', 's-1') - - CS%id_Kd_itidal = register_diag_field('ocean_model','Kd_itides',diag%axesTi,Time, & - 'Internal Tide Driven Diffusivity', 'm2 s-1') - - CS%id_Kd_lowmode = register_diag_field('ocean_model','Kd_lowmode',diag%axesTi,Time, & - 'Internal Tide Driven Diffusivity (from propagating low modes)', 'm2 s-1') - - CS%id_Fl_itidal = register_diag_field('ocean_model','Fl_itides',diag%axesTi,Time, & - 'Vertical flux of tidal turbulent dissipation', 'm3 s-3') - - CS%id_Fl_lowmode = register_diag_field('ocean_model','Fl_lowmode',diag%axesTi,Time, & - 'Vertical flux of tidal turbulent dissipation (from propagating low modes)', 'm3 s-3') - - CS%id_Polzin_decay_scale = register_diag_field('ocean_model','Polzin_decay_scale',diag%axesT1,Time, & - 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', 'm') - - CS%id_Polzin_decay_scale_scaled = register_diag_field('ocean_model','Polzin_decay_scale_scaled',diag%axesT1,Time, & - 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme, scaled by N2_bot/N2_meanz', 'm') - - CS%id_N2_bot = register_diag_field('ocean_model','N2_b',diag%axesT1,Time, & - 'Bottom Buoyancy frequency squared', 's-2') - - CS%id_N2_meanz = register_diag_field('ocean_model','N2_meanz',diag%axesT1,Time, & - 'Buoyancy frequency squared averaged over the water column', 's-2') - - CS%id_Kd_Work = register_diag_field('ocean_model','Kd_Work',diag%axesTL,Time, & - 'Work done by Diapycnal Mixing', 'W m-2') - - CS%id_Kd_Itidal_Work = register_diag_field('ocean_model','Kd_Itidal_Work',diag%axesTL,Time, & - 'Work done by Internal Tide Diapycnal Mixing', 'W m-2') - - CS%id_Kd_Niku_Work = register_diag_field('ocean_model','Kd_Nikurashin_Work',diag%axesTL,Time, & - 'Work done by Nikurashin Lee Wave Drag Scheme', 'W m-2') - - CS%id_Kd_Lowmode_Work = register_diag_field('ocean_model','Kd_Lowmode_Work',diag%axesTL,Time, & - 'Work done by Internal Tide Diapycnal Mixing (low modes)', 'W m-2') - CS%id_N2 = register_diag_field('ocean_model','N2',diag%axesTi,Time, & 'Buoyancy frequency squared', 's-2', cmor_field_name='obvfsq', & cmor_long_name='Square of seawater buoyancy frequency',& @@ -3113,25 +2071,13 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp vd = var_desc("N2", "s-2",& "Buoyancy frequency, interpolated to z", z_grid='z') CS%id_N2_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("Kd_itides","m2 s-1", & - "Internal Tide Driven Diffusivity, interpolated to z", z_grid='z') - CS%id_Kd_itidal_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - if (CS%Lee_wave_dissipation) then - vd = var_desc("Kd_Nikurashin", "m2 s-1", & - "Lee Wave Driven Diffusivity, interpolated to z", z_grid='z') - CS%id_Kd_Niku_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - endif - if (CS%Lowmode_itidal_dissipation) then - vd = var_desc("Kd_lowmode","m2 s-1", & - "Internal Tide Driven Diffusivity (from low modes), interpolated to z",& - z_grid='z') - CS%id_Kd_lowmode_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - endif if (CS%user_change_diff) & CS%id_Kd_user_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) endif endif + + ! GMM, the following should be moved to the DD module call get_param(param_file, mdl, "DOUBLE_DIFFUSION", CS%double_diffusion, & "If true, increase diffusivitives for temperature or salt \n"//& "based on double-diffusive paramaterization from MOM4/KPP.", & @@ -3169,31 +2115,37 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp endif endif - if (CS%Int_tide_dissipation .and. CS%Bryan_Lewis_diffusivity) & - call MOM_error(FATAL,"MOM_Set_Diffusivity: "// & - "Bryan-Lewis and internal tidal dissipation are both enabled. Choose one.") - if (CS%Henyey_IGW_background .and. CS%Kd_tanh_lat_fn) call MOM_error(FATAL, & - "Set_diffusivity: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.") - if (CS%user_change_diff) then call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) endif + if (CS%tm_csp%Int_tide_dissipation .and. CS%bkgnd_mixing_csp%Bryan_Lewis_diffusivity) & + call MOM_error(FATAL,"MOM_Set_Diffusivity: "// & + "Bryan-Lewis and internal tidal dissipation are both enabled. Choose one.") + CS%useKappaShear = kappa_shear_init(Time, G, GV, param_file, CS%diag, CS%kappaShear_CSp) if (CS%useKappaShear) & id_clock_kappaShear = cpu_clock_id('(Ocean kappa_shear)', grain=CLOCK_MODULE) - CS%useCVMix = CVMix_shear_init(Time, G, GV, param_file, CS%diag, CS%CVMix_shear_CSp) - + ! CVMix shear-driven mixing + CS%use_CVMix_shear = CVMix_shear_init(Time, G, GV, param_file, CS%diag, CS%CVMix_shear_csp) end subroutine set_diffusivity_init +!> Clear pointers and dealocate memory subroutine set_diffusivity_end(CS) - type(set_diffusivity_CS), pointer :: CS + type(set_diffusivity_CS), pointer :: CS !< Control structure for this module + + if (.not.associated(CS)) return + + call bkgnd_mixing_end(CS%bkgnd_mixing_csp) if (CS%user_change_diff) & call user_change_diff_end(CS%user_change_diff_CSp) + if (CS%use_CVMix_shear) & + call CVMix_shear_end(CS%CVMix_shear_csp) + if (associated(CS)) deallocate(CS) end subroutine set_diffusivity_end diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 394b17dbd2..90401313dc 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -45,6 +45,7 @@ module MOM_set_visc use MOM_hor_index, only : hor_index_type use MOM_kappa_shear, only : kappa_shear_is_used use MOM_CVMix_shear, only : CVMix_shear_is_used +use MOM_CVMix_conv, only : CVMix_conv_is_used use MOM_io, only : vardesc, var_desc use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_variables, only : thermo_var_ptrs @@ -1783,20 +1784,21 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) ! (in) restart_CS - A pointer to the restart control structure. type(vardesc) :: vd logical :: use_kappa_shear, adiabatic, useKPP, useEPBL - logical :: use_CVMix, MLE_use_PBL_MLD + logical :: use_CVMix_shear, MLE_use_PBL_MLD, use_CVMix_conv integer :: isd, ied, jsd, jed, nz character(len=40) :: mdl = "MOM_set_visc" ! This module's name. isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., & do_not_log=.true.) - use_kappa_shear = .false. ; use_CVMix = .false. ; - useKPP = .false. ; useEPBL = .false. + use_kappa_shear = .false. ; use_CVMix_shear = .false. ; + useKPP = .false. ; useEPBL = .false. ; use_CVMix_conv = .false. ; if (.not.adiabatic) then use_kappa_shear = kappa_shear_is_used(param_file) - use_CVMix = CVMix_shear_is_used(param_file) + use_CVMix_shear = CVMix_shear_is_used(param_file) + use_CVMix_conv = CVMix_conv_is_used(param_file) call get_param(param_file, mdl, "USE_KPP", useKPP, & - "If true, turns on the [CVmix] KPP scheme of Large et al., 1984,\n"// & + "If true, turns on the [CVMix] KPP scheme of Large et al., 1984,\n"// & "to calculate diffusivities and non-local transport in the OBL.", & default=.false., do_not_log=.true.) call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useEPBL, & @@ -1805,21 +1807,26 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) "in the surface boundary layer.", default=.false., do_not_log=.true.) endif - if (use_kappa_shear .or. useKPP .or. useEPBL .or. use_CVMix) then - allocate(visc%Kd_turb(isd:ied,jsd:jed,nz+1)) ; visc%Kd_turb(:,:,:) = 0.0 + if (use_kappa_shear .or. useKPP .or. useEPBL .or. use_CVMix_shear .or. use_CVMix_conv) then + allocate(visc%Kd_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kd_shear(:,:,:) = 0.0 allocate(visc%TKE_turb(isd:ied,jsd:jed,nz+1)) ; visc%TKE_turb(:,:,:) = 0.0 - allocate(visc%Kv_turb(isd:ied,jsd:jed,nz+1)) ; visc%Kv_turb(:,:,:) = 0.0 + allocate(visc%Kv_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kv_shear(:,:,:) = 0.0 + allocate(visc%Kv_slow(isd:ied,jsd:jed,nz+1)) ; visc%Kv_slow(:,:,:) = 0.0 - vd = var_desc("Kd_turb","m2 s-1","Turbulent diffusivity at interfaces", & + vd = var_desc("Kd_shear","m2 s-1","Shear-driven turbulent diffusivity at interfaces", & hor_grid='h', z_grid='i') - call register_restart_field(visc%Kd_turb, vd, .false., restart_CS) + call register_restart_field(visc%Kd_shear, vd, .false., restart_CS) vd = var_desc("TKE_turb","m2 s-2","Turbulent kinetic energy per unit mass at interfaces", & hor_grid='h', z_grid='i') call register_restart_field(visc%TKE_turb, vd, .false., restart_CS) - vd = var_desc("Kv_turb","m2 s-1","Turbulent viscosity at interfaces", & + vd = var_desc("Kv_shear","m2 s-1","Shear-driven turbulent viscosity at interfaces", & hor_grid='h', z_grid='i') - call register_restart_field(visc%Kv_turb, vd, .false., restart_CS) + call register_restart_field(visc%Kv_shear, vd, .false., restart_CS) + vd = var_desc("Kv_slow","m2 s-1","Vertical turbulent viscosity at interfaces due \n" // & + " to slow processes", hor_grid='h', z_grid='i') + call register_restart_field(visc%Kv_slow, vd, .false., restart_CS) + endif ! visc%MLD is used to communicate the state of the (e)PBL to the rest of the model @@ -2090,9 +2097,10 @@ subroutine set_visc_end(visc, CS) if (CS%dynamic_viscous_ML) then deallocate(visc%nkml_visc_u) ; deallocate(visc%nkml_visc_v) endif - if (associated(visc%Kd_turb)) deallocate(visc%Kd_turb) + if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear) + if (associated(visc%Kv_slow)) deallocate(visc%Kv_slow) if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb) - if (associated(visc%Kv_turb)) deallocate(visc%Kv_turb) + if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear) if (associated(visc%ustar_bbl)) deallocate(visc%ustar_bbl) if (associated(visc%TKE_bbl)) deallocate(visc%TKE_bbl) if (associated(visc%taux_shelf)) deallocate(visc%taux_shelf) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 new file mode 100644 index 0000000000..5524ef074a --- /dev/null +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -0,0 +1,1330 @@ +!> Interface to vertical tidal mixing schemes including CVMix tidal mixing. +module MOM_tidal_mixing + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field +use MOM_diag_mediator, only : safe_alloc_ptr, post_data +use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag +use MOM_diag_to_Z, only : calc_Zint_diags +use MOM_EOS, only : calculate_density +use MOM_variables, only : thermo_var_ptrs, p3d +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_debugging, only : hchksum +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_string_functions, only : uppercase, lowercase +use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc +use CVMix_tidal, only : CVMix_init_tidal, CVMix_compute_Simmons_invariant +use CVMix_tidal, only : CVMix_coeffs_tidal, CVMix_tidal_params_type +use CVMix_kinds_and_types, only : CVMix_global_params_type +use CVMix_put_get, only : CVMix_put + +implicit none ; private + +#include + +public tidal_mixing_init +public setup_tidal_diagnostics +public calculate_tidal_mixing +public post_tidal_diagnostics +public tidal_mixing_end + +!> Containers for tidal mixing diagnostics +type, public :: tidal_mixing_diags + real, pointer, dimension(:,:,:) :: & + Kd_itidal => NULL(),& ! internal tide diffusivity at interfaces (m2/s) + Fl_itidal => NULL(),& ! vertical flux of tidal turbulent dissipation (m3/s3) + Kd_lowmode => NULL(),& ! internal tide diffusivity at interfaces + ! due to propagating low modes (m2/s) (BDM) + Fl_lowmode => NULL(),& ! vertical flux of tidal turbulent dissipation + ! due to propagating low modes (m3/s3) (BDM) + Kd_Niku => NULL(),& ! lee-wave diffusivity at interfaces (m2/s) + Kd_Niku_work => NULL(),& ! layer integrated work by lee-wave driven mixing (W/m2) + Kd_Itidal_Work => NULL(),& ! layer integrated work by int tide driven mixing (W/m2) + Kd_Lowmode_Work => NULL(),& ! layer integrated work by low mode driven mixing (W/m2) BDM + N2_int => NULL(),& + vert_dep_3d => NULL() + + real, pointer, dimension(:,:) :: & + TKE_itidal_used => NULL(),& ! internal tide TKE input at ocean bottom (W/m2) + N2_bot => NULL(),& ! bottom squared buoyancy frequency (1/s2) + N2_meanz => NULL(),& ! vertically averaged buoyancy frequency (1/s2) + Polzin_decay_scale_scaled => NULL(),& ! vertical scale of decay for tidal dissipation + Polzin_decay_scale => NULL(),& ! vertical decay scale for tidal diss with Polzin (meter) + Simmons_coeff_2d => NULL() + +end type + +!> Control structure for tidal mixing module. +type, public :: tidal_mixing_cs + logical :: debug = .true. + + ! Parameters + logical :: int_tide_dissipation = .false. ! Internal tide conversion (from barotropic) + ! with the schemes of St Laurent et al (2002)/ + ! Simmons et al (2004) + + integer :: Int_tide_profile ! A coded integer indicating the vertical profile + ! for dissipation of the internal waves. Schemes that + ! are currently encoded are St Laurent et al (2002) and + ! Polzin (2009). + logical :: Lee_wave_dissipation = .false. ! Enable lee-wave driven mixing, following + ! Nikurashin (2010), with a vertical energy + ! deposition profile specified by Lee_wave_profile. + ! St Laurent et al (2002) or + ! Simmons et al (2004) scheme + + integer :: Lee_wave_profile ! A coded integer indicating the vertical profile + ! for dissipation of the lee waves. Schemes that are + ! currently encoded are St Laurent et al (2002) and + ! Polzin (2009). + real :: Int_tide_decay_scale ! decay scale for internal wave TKE (meter) + + real :: Mu_itides ! efficiency for conversion of dissipation + ! to potential energy (nondimensional) + + real :: Gamma_itides ! fraction of local dissipation (nondimensional) + + real :: Gamma_lee ! fraction of local dissipation for lee waves + ! (Nikurashin's energy input) (nondimensional) + real :: Decay_scale_factor_lee ! Scaling factor for the decay scale of lee + ! wave energy dissipation (nondimensional) + + real :: min_zbot_itides ! minimum depth for internal tide conversion (meter) + logical :: Lowmode_itidal_dissipation = .false. ! Internal tide conversion (from low modes) + ! with the schemes of St Laurent et al (2002)/ + ! Simmons et al (2004) !BDM + + real :: Nu_Polzin ! The non-dimensional constant used in Polzin form of + ! the vertical scale of decay of tidal dissipation + + real :: Nbotref_Polzin ! Reference value for the buoyancy frequency at the + ! ocean bottom used in Polzin formulation of the + ! vertical scale of decay of tidal dissipation (1/s) + real :: Polzin_decay_scale_factor ! Scaling factor for the decay length scale + ! of the tidal dissipation profile in Polzin + ! (nondimensional) + real :: Polzin_decay_scale_max_factor ! The decay length scale of tidal + ! dissipation profile in Polzin formulation should not + ! exceed Polzin_decay_scale_max_factor * depth of the + ! ocean (nondimensional). + real :: Polzin_min_decay_scale ! minimum decay scale of the tidal dissipation + ! profile in Polzin formulation (meter) + + real :: TKE_itide_max ! maximum internal tide conversion (W m-2) + ! available to mix above the BBL + + real :: utide ! constant tidal amplitude (m s-1) used if + real :: kappa_itides ! topographic wavenumber and non-dimensional scaling + real :: kappa_h2_factor ! factor for the product of wavenumber * rms sgs height + character(len=200) :: inputdir + + logical :: use_CVMix_tidal = .false. ! true if CVMix is to be used for determining + ! diffusivity due to tidal mixing + + real :: min_thickness ! Minimum thickness allowed [m] + + ! CVMix-specific parameters + type(CVMix_tidal_params_type) :: CVMix_tidal_params + type(CVMix_global_params_type) :: CVMix_glb_params ! for Prandtl number only + real :: tidal_max_coef ! maximum allowable tidal diffusivity. [m^2/s] + + ! Data containers + real, pointer, dimension(:,:) :: TKE_Niku => NULL() + real, pointer, dimension(:,:) :: TKE_itidal => NULL() + real, pointer, dimension(:,:) :: Nb => NULL() + real, pointer, dimension(:,:) :: mask_itidal => NULL() + real, pointer, dimension(:,:) :: h2 => NULL() + real, pointer, dimension(:,:) :: tideamp => NULL() ! RMS tidal amplitude (m/s) + real, allocatable,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) + + ! Diagnostics + type(diag_ctrl), pointer :: diag => NULL() ! structure to regulate diagn output timing + type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() + type(tidal_mixing_diags), pointer :: dd => NULL() + + ! Diagnostic identifiers + integer :: id_TKE_itidal = -1 + integer :: id_TKE_leewave = -1 + integer :: id_Kd_itidal = -1 + integer :: id_Kd_Niku = -1 + integer :: id_Kd_lowmode = -1 + integer :: id_Kd_itidal_z = -1 + integer :: id_Kd_Niku_z = -1 + integer :: id_Kd_lowmode_z = -1 + integer :: id_Kd_Itidal_Work = -1 + integer :: id_Kd_Niku_Work = -1 + integer :: id_Kd_Lowmode_Work = -1 + integer :: id_Nb = -1 + integer :: id_N2_bot = -1 + integer :: id_N2_meanz = -1 + integer :: id_Fl_itidal = -1 + integer :: id_Fl_lowmode = -1 + integer :: id_Polzin_decay_scale = -1 + integer :: id_Polzin_decay_scale_scaled = -1 + integer :: id_N2_int = -1 + integer :: id_Simmons_coeff = -1 + integer :: id_vert_dep = -1 + +end type tidal_mixing_cs + +character(len=40) :: mdl = "MOM_tidal_mixing" !< This module's name. +character*(20), parameter :: STLAURENT_PROFILE_STRING = "STLAURENT_02" +character*(20), parameter :: POLZIN_PROFILE_STRING = "POLZIN_09" +character*(20), parameter :: SIMMONS_PROFILE_STRING = "SIMMONS" +character*(20), parameter :: SCHMITTNER_PROFILE_STRING = "SCHMITTNER" +integer, parameter :: STLAURENT_02 = 1 +integer, parameter :: POLZIN_09 = 2 +integer, parameter :: SIMMONS_04 = 3 +integer, parameter :: SCHMITTNER = 4 + +contains + +!> Initializes internal tidal dissipation scheme for diapycnal mixing +logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, CS) + + type(time_type), intent(in) :: Time !< The current time. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. + type(diag_to_Z_CS), pointer :: diag_to_Z_CSp !< pointer to the Z-diagnostics control + type(tidal_mixing_cs), pointer :: CS !< This module's control structure. + + ! Local variables + logical :: read_tideamp + character(len=20) :: tmpstr, int_tide_profile_str + character(len=20) :: default_profile_string, tidal_energy_type + character(len=200) :: filename, h2_file, Niku_TKE_input_file + character(len=200) :: tidal_energy_file, tideamp_file + type(vardesc) :: vd + real :: utide, zbot, hamp, prandtl_tidal + real :: Niku_scale ! local variable for scaling the Nikurashin TKE flux data + integer :: i, j, is, ie, js, je + integer :: isd, ied, jsd, jed + +! This include declares and sets the variable "version". +#include "version_variable.h" + + if (associated(CS)) then + call MOM_error(WARNING, "tidal_mixing_init called when control structure "// & + "is already associated.") + return + endif + allocate(CS) + allocate(CS%dd) + + CS%debug = CS%debug.and.is_root_pe() + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + + CS%diag => diag + if (associated(diag_to_Z_CSp)) CS%diag_to_Z_CSp => diag_to_Z_CSp + + ! Read parameters + call log_version(param_file, mdl, version, & + "Vertical Tidal Mixing Parameterization") + call get_param(param_file, mdl, "USE_CVMix_TIDAL", CS%use_CVMix_tidal, & + "If true, turns on tidal mixing via CVMix", & + default=.false.) + + call get_param(param_file, mdl, "INPUTDIR", CS%inputdir, default=".",do_not_log=.true.) + CS%inputdir = slasher(CS%inputdir) + call get_param(param_file, mdl, "INT_TIDE_DISSIPATION", CS%int_tide_dissipation, & + "If true, use an internal tidal dissipation scheme to \n"//& + "drive diapycnal mixing, along the lines of St. Laurent \n"//& + "et al. (2002) and Simmons et al. (2004).", default=CS%use_CVMix_tidal) + + ! return if tidal mixing is inactive + tidal_mixing_init = CS%int_tide_dissipation + if (.not. tidal_mixing_init) return + + if (CS%int_tide_dissipation) then + default_profile_string = STLAURENT_PROFILE_STRING + if (CS%use_CVMix_tidal) default_profile_string = SIMMONS_PROFILE_STRING + call get_param(param_file, mdl, "INT_TIDE_PROFILE", int_tide_profile_str, & + "INT_TIDE_PROFILE selects the vertical profile of energy \n"//& + "dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//& + "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//& + "\t decay profile.\n"//& + "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//& + "\t decay profile.", & + default=default_profile_string) + ! TODO: list the newly available profile selections + int_tide_profile_str = uppercase(int_tide_profile_str) + select case (int_tide_profile_str) + case (STLAURENT_PROFILE_STRING) ; CS%int_tide_profile = STLAURENT_02 + case (POLZIN_PROFILE_STRING) ; CS%int_tide_profile = POLZIN_09 + case (SIMMONS_PROFILE_STRING) ; CS%int_tide_profile = SIMMONS_04 + case (SCHMITTNER_PROFILE_STRING) ; CS%int_tide_profile = SCHMITTNER + case default + call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & + "#define INT_TIDE_PROFILE "//trim(int_tide_profile_str)//" found in input file.") + end select + + ! Check profile consistency + if (CS%use_CVMix_tidal .and. (CS%int_tide_profile.eq.STLAURENT_02 .or. & + CS%int_tide_profile.eq.POLZIN_09)) then + call MOM_error(FATAL, "tidal_mixing_init: Tidal mixing profile"// & + " "//trim(int_tide_profile_str)//" unavailable in CVMix. Available "//& + "profiles in CVMix are "//trim(SIMMONS_PROFILE_STRING)//" and "//& + trim(SCHMITTNER_PROFILE_STRING)//".") + else if (.not.CS%use_CVMix_tidal .and. (CS%int_tide_profile.eq.SIMMONS_04.or. & + CS%int_tide_profile.eq.SCHMITTNER)) then + call MOM_error(FATAL, "tidal_mixing_init: Tidal mixing profiles "// & + trim(SIMMONS_PROFILE_STRING)//" and "//trim(SCHMITTNER_PROFILE_STRING)//& + " are available only when USE_CVMix_TIDAL is True.") + endif + + else if (CS%use_CVMix_tidal) then + call MOM_error(FATAL, "tidal_mixing_init: Cannot set INT_TIDE_DISSIPATION to False "// & + "when USE_CVMix_TIDAL is set to True.") + endif + + call get_param(param_file, mdl, "LEE_WAVE_DISSIPATION", CS%Lee_wave_dissipation, & + "If true, use an lee wave driven dissipation scheme to \n"//& + "drive diapycnal mixing, along the lines of Nikurashin \n"//& + "(2010) and using the St. Laurent et al. (2002) \n"//& + "and Simmons et al. (2004) vertical profile", default=.false.) + if (CS%lee_wave_dissipation) then + if (CS%use_CVMix_tidal) then + call MOM_error(FATAL, "tidal_mixing_init: Lee wave driven dissipation scheme cannot "// & + "be used when CVMix tidal mixing scheme is active.") + end if + call get_param(param_file, mdl, "LEE_WAVE_PROFILE", tmpstr, & + "LEE_WAVE_PROFILE selects the vertical profile of energy \n"//& + "dissipation with LEE_WAVE_DISSIPATION. Valid values are:\n"//& + "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//& + "\t decay profile.\n"//& + "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//& + "\t decay profile.", & + default=STLAURENT_PROFILE_STRING) + tmpstr = uppercase(tmpstr) + select case (tmpstr) + case (STLAURENT_PROFILE_STRING) ; CS%lee_wave_profile = STLAURENT_02 + case (POLZIN_PROFILE_STRING) ; CS%lee_wave_profile = POLZIN_09 + case default + call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & + "#define LEE_WAVE_PROFILE "//trim(tmpstr)//" found in input file.") + end select + endif + + call get_param(param_file, mdl, "INT_TIDE_LOWMODE_DISSIPATION", CS%Lowmode_itidal_dissipation, & + "If true, consider mixing due to breaking low modes that \n"//& + "have been remotely generated; as with itidal drag on the \n"//& + "barotropic tide, use an internal tidal dissipation scheme to \n"//& + "drive diapycnal mixing, along the lines of St. Laurent \n"//& + "et al. (2002) and Simmons et al. (2004).", default=.false.) + + if ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09)) .or. & + (CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09))) then + if (CS%use_CVMix_tidal) then + call MOM_error(FATAL, "tidal_mixing_init: Polzin scheme cannot "// & + "be used when CVMix tidal mixing scheme is active.") + end if + call get_param(param_file, mdl, "NU_POLZIN", CS%Nu_Polzin, & + "When the Polzin decay profile is used, this is a \n"//& + "non-dimensional constant in the expression for the \n"//& + "vertical scale of decay for the tidal energy dissipation.", & + units="nondim", default=0.0697) + call get_param(param_file, mdl, "NBOTREF_POLZIN", CS%Nbotref_Polzin, & + "When the Polzin decay profile is used, this is the \n"//& + "Rreference value of the buoyancy frequency at the ocean \n"//& + "bottom in the Polzin formulation for the vertical \n"//& + "scale of decay for the tidal energy dissipation.", & + units="s-1", default=9.61e-4) + call get_param(param_file, mdl, "POLZIN_DECAY_SCALE_FACTOR", & + CS%Polzin_decay_scale_factor, & + "When the Polzin decay profile is used, this is a \n"//& + "scale factor for the vertical scale of decay of the tidal \n"//& + "energy dissipation.", default=1.0, units="nondim") + call get_param(param_file, mdl, "POLZIN_SCALE_MAX_FACTOR", & + CS%Polzin_decay_scale_max_factor, & + "When the Polzin decay profile is used, this is a factor \n"//& + "to limit the vertical scale of decay of the tidal \n"//& + "energy dissipation to POLZIN_DECAY_SCALE_MAX_FACTOR \n"//& + "times the depth of the ocean.", units="nondim", default=1.0) + call get_param(param_file, mdl, "POLZIN_MIN_DECAY_SCALE", CS%Polzin_min_decay_scale, & + "When the Polzin decay profile is used, this is the \n"//& + "minimum vertical decay scale for the vertical profile\n"//& + "of internal tide dissipation with the Polzin (2009) formulation", & + units="m", default=0.0) + endif + + if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation) then + call get_param(param_file, mdl, "INT_TIDE_DECAY_SCALE", CS%Int_tide_decay_scale, & + "The decay scale away from the bottom for tidal TKE with \n"//& + "the new coding when INT_TIDE_DISSIPATION is used.", & + !units="m", default=0.0) + units="m", default=500.0) ! TODO: confirm this new default + call get_param(param_file, mdl, "MU_ITIDES", CS%Mu_itides, & + "A dimensionless turbulent mixing efficiency used with \n"//& + "INT_TIDE_DISSIPATION, often 0.2.", units="nondim", default=0.2) + call get_param(param_file, mdl, "GAMMA_ITIDES", CS%Gamma_itides, & + "The fraction of the internal tidal energy that is \n"//& + "dissipated locally with INT_TIDE_DISSIPATION. \n"//& + "THIS NAME COULD BE BETTER.", & + units="nondim", default=0.3333) + call get_param(param_file, mdl, "MIN_ZBOT_ITIDES", CS%min_zbot_itides, & + "Turn off internal tidal dissipation when the total \n"//& + "ocean depth is less than this value.", units="m", default=0.0) + endif + + if ( (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation) .and. & + .not. CS%use_CVMix_tidal) then + + call safe_alloc_ptr(CS%Nb,isd,ied,jsd,jed) + call safe_alloc_ptr(CS%h2,isd,ied,jsd,jed) + call safe_alloc_ptr(CS%TKE_itidal,isd,ied,jsd,jed) + call safe_alloc_ptr(CS%mask_itidal,isd,ied,jsd,jed) ; CS%mask_itidal(:,:) = 1.0 + + call get_param(param_file, mdl, "KAPPA_ITIDES", CS%kappa_itides, & + "A topographic wavenumber used with INT_TIDE_DISSIPATION. \n"//& + "The default is 2pi/10 km, as in St.Laurent et al. 2002.", & + units="m-1", default=8.e-4*atan(1.0)) + + call get_param(param_file, mdl, "UTIDE", CS%utide, & + "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & + units="m s-1", default=0.0) + call safe_alloc_ptr(CS%tideamp,is,ie,js,je) ; CS%tideamp(:,:) = CS%utide + + call get_param(param_file, mdl, "KAPPA_H2_FACTOR", CS%kappa_h2_factor, & + "A scaling factor for the roughness amplitude with n"//& + "INT_TIDE_DISSIPATION.", units="nondim", default=1.0) + call get_param(param_file, mdl, "TKE_ITIDE_MAX", CS%TKE_itide_max, & + "The maximum internal tide energy source availble to mix \n"//& + "above the bottom boundary layer with INT_TIDE_DISSIPATION.", & + units="W m-2", default=1.0e3) + + call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, & + "If true, read a file (given by TIDEAMP_FILE) containing \n"//& + "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.) + if (read_tideamp) then + if (CS%use_CVMix_tidal) then + call MOM_error(FATAL, "tidal_mixing_init: Tidal amplitude files are "// & + "not compatible with CVMix tidal mixing. ") + end if + call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, & + "The path to the file containing the spatially varying \n"//& + "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc") + filename = trim(CS%inputdir) // trim(tideamp_file) + call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename) + call MOM_read_data(filename, 'tideamp', CS%tideamp, G%domain, timelevel=1) + endif + + call get_param(param_file, mdl, "H2_FILE", h2_file, & + "The path to the file containing the sub-grid-scale \n"//& + "topographic roughness amplitude with INT_TIDE_DISSIPATION.", & + fail_if_missing=(.not.CS%use_CVMix_tidal)) + filename = trim(CS%inputdir) // trim(h2_file) + call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename) + call MOM_read_data(filename, 'h2', CS%h2, G%domain, timelevel=1) + + do j=js,je ; do i=is,ie + if (G%bathyT(i,j) < CS%min_zbot_itides) CS%mask_itidal(i,j) = 0.0 + CS%tideamp(i,j) = CS%tideamp(i,j) * CS%mask_itidal(i,j) * G%mask2dT(i,j) + + ! Restrict rms topo to 10 percent of column depth. + zbot = G%bathyT(i,j) + hamp = sqrt(CS%h2(i,j)) + hamp = min(0.1*zbot,hamp) + CS%h2(i,j) = hamp*hamp + + utide = CS%tideamp(i,j) + ! Compute the fixed part of internal tidal forcing; units are [kg s-2] here. + CS%TKE_itidal(i,j) = 0.5*CS%kappa_h2_factor*GV%Rho0*& + CS%kappa_itides*CS%h2(i,j)*utide*utide + enddo; enddo + + endif + + if (CS%Lee_wave_dissipation) then + + call get_param(param_file, mdl, "NIKURASHIN_TKE_INPUT_FILE",Niku_TKE_input_file, & + "The path to the file containing the TKE input from lee \n"//& + "wave driven mixing. Used with LEE_WAVE_DISSIPATION.", & + fail_if_missing=.true.) + call get_param(param_file, mdl, "NIKURASHIN_SCALE",Niku_scale, & + "A non-dimensional factor by which to scale the lee-wave \n"//& + "driven TKE input. Used with LEE_WAVE_DISSIPATION.", & + units="nondim", default=1.0) + + filename = trim(CS%inputdir) // trim(Niku_TKE_input_file) + call log_param(param_file, mdl, "INPUTDIR/NIKURASHIN_TKE_INPUT_FILE", & + filename) + call safe_alloc_ptr(CS%TKE_Niku,is,ie,js,je); CS%TKE_Niku(:,:) = 0.0 + call MOM_read_data(filename, 'TKE_input', CS%TKE_Niku, G%domain, timelevel=1 ) ! ??? timelevel -aja + CS%TKE_Niku(:,:) = Niku_scale * CS%TKE_Niku(:,:) + + call get_param(param_file, mdl, "GAMMA_NIKURASHIN",CS%Gamma_lee, & + "The fraction of the lee wave energy that is dissipated \n"//& + "locally with LEE_WAVE_DISSIPATION.", units="nondim", & + default=0.3333) + call get_param(param_file, mdl, "DECAY_SCALE_FACTOR_LEE",CS%Decay_scale_factor_lee, & + "Scaling for the vertical decay scaleof the local \n"//& + "dissipation of lee waves dissipation.", units="nondim", & + default=1.0) + else + CS%Decay_scale_factor_lee = -9.e99 ! This should never be used if CS%Lee_wave_dissipation = False + endif + + ! Configure CVMix + if (CS%use_CVMix_tidal) then + + ! Read in CVMix params + !call openParameterBlock(param_file,'CVMix_TIDAL') + call get_param(param_file, mdl, "TIDAL_MAX_COEF", CS%tidal_max_coef, & + "largest acceptable value for tidal diffusivity", & + units="m^2/s", default=50e-4) ! the default is 50e-4 in CVMix, 100e-4 in POP. + call get_param(param_file, mdl, "TIDAL_ENERGY_FILE",tidal_energy_file, & + "The path to the file containing tidal energy \n"//& + "dissipation. Used with CVMix tidal mixing schemes.", & + fail_if_missing=.true.) + tidal_energy_file = trim(CS%inputdir) // trim(tidal_energy_file) + call get_param(param_file, mdl, "TIDAL_ENERGY_TYPE",tidal_energy_type, & + "The type of input tidal energy flux dataset.",& + fail_if_missing=.true.) + ! TODO: list all available tidal energy types here + call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, & + do_not_log=.True.) + call get_param(param_file, mdl, "PRANDTL_TIDAL", prandtl_tidal, & + "Prandtl number used by CVMix tidal mixing schemes \n"//& + "to convert vertical diffusivities into viscosities.", & + units="nondim", default=1.0, & + do_not_log=.true.) + call CVMix_put(CS%CVMix_glb_params,'Prandtl',prandtl_tidal) + + int_tide_profile_str = lowercase(int_tide_profile_str) + + + ! TODO: check parameter consistency. (see POP::tidal_mixing.F90::tidal_check) + + ! Set up CVMix + call CVMix_init_tidal(CVMix_tidal_params_user = CS%CVMix_tidal_params, & + mix_scheme = int_tide_profile_str, & + efficiency = CS%Mu_itides, & + vertical_decay_scale = CS%int_tide_decay_scale, & + max_coefficient = CS%tidal_max_coef, & + local_mixing_frac = CS%Gamma_itides, & + depth_cutoff = CS%min_zbot_itides) + + call read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) + + !call closeParameterBlock(param_file) + + endif ! CVMix on + + ! Register Diagnostics fields + + if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. & + CS%Lowmode_itidal_dissipation) then + + CS%id_Kd_itidal = register_diag_field('ocean_model','Kd_itides',diag%axesTi,Time, & + 'Internal Tide Driven Diffusivity', 'm2 s-1') + + if (CS%use_CVMix_tidal) then + CS%id_N2_int = register_diag_field('ocean_model','N2_int',diag%axesTi,Time, & + 'Bouyancy frequency squared, at interfaces', 's-2') + ! TODO: add units + CS%id_Simmons_coeff = register_diag_field('ocean_model','Simmons_coeff',diag%axesT1,Time, & + 'time-invariant portion of the tidal mixing coefficient using the Simmons', '') + CS%id_vert_dep = register_diag_field('ocean_model','vert_dep',diag%axesTi,Time, & + 'vertical deposition function needed for Simmons et al tidal mixing', '') + + else + CS%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal',diag%axesT1,Time, & + 'Internal Tide Driven Turbulent Kinetic Energy', 'W m-2') + CS%id_Nb = register_diag_field('ocean_model','Nb',diag%axesT1,Time, & + 'Bottom Buoyancy Frequency', 's-1') + + CS%id_Kd_lowmode = register_diag_field('ocean_model','Kd_lowmode',diag%axesTi,Time, & + 'Internal Tide Driven Diffusivity (from propagating low modes)', 'm2 s-1') + + CS%id_Fl_itidal = register_diag_field('ocean_model','Fl_itides',diag%axesTi,Time, & + 'Vertical flux of tidal turbulent dissipation', 'm3 s-3') + + CS%id_Fl_lowmode = register_diag_field('ocean_model','Fl_lowmode',diag%axesTi,Time, & + 'Vertical flux of tidal turbulent dissipation (from propagating low modes)', 'm3 s-3') + + CS%id_Polzin_decay_scale = register_diag_field('ocean_model','Polzin_decay_scale',diag%axesT1,Time, & + 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', 'm') + + CS%id_Polzin_decay_scale_scaled = register_diag_field('ocean_model','Polzin_decay_scale_scaled',diag%axesT1,Time, & + 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme, scaled by N2_bot/N2_meanz', 'm') + + CS%id_N2_bot = register_diag_field('ocean_model','N2_b',diag%axesT1,Time, & + 'Bottom Buoyancy frequency squared', 's-2') + + CS%id_N2_meanz = register_diag_field('ocean_model','N2_meanz',diag%axesT1,Time, & + 'Buoyancy frequency squared averaged over the water column', 's-2') + + CS%id_Kd_Itidal_Work = register_diag_field('ocean_model','Kd_Itidal_Work',diag%axesTL,Time, & + 'Work done by Internal Tide Diapycnal Mixing', 'W m-2') + + CS%id_Kd_Niku_Work = register_diag_field('ocean_model','Kd_Nikurashin_Work',diag%axesTL,Time, & + 'Work done by Nikurashin Lee Wave Drag Scheme', 'W m-2') + + CS%id_Kd_Lowmode_Work = register_diag_field('ocean_model','Kd_Lowmode_Work',diag%axesTL,Time, & + 'Work done by Internal Tide Diapycnal Mixing (low modes)', 'W m-2') + + if (CS%Lee_wave_dissipation) then + CS%id_TKE_leewave = register_diag_field('ocean_model','TKE_leewave',diag%axesT1,Time, & + 'Lee wave Driven Turbulent Kinetic Energy', 'W m-2') + CS%id_Kd_Niku = register_diag_field('ocean_model','Kd_Nikurashin',diag%axesTi,Time, & + 'Lee Wave Driven Diffusivity', 'm2 s-1') + endif + endif ! S%use_CVMix_tidal + + if (associated(CS%diag_to_Z_CSp)) then + vd = var_desc("Kd_itides","m2 s-1", & + "Internal Tide Driven Diffusivity, interpolated to z", z_grid='z') + CS%id_Kd_itidal_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + if (CS%Lee_wave_dissipation) then + vd = var_desc("Kd_Nikurashin", "m2 s-1", & + "Lee Wave Driven Diffusivity, interpolated to z", z_grid='z') + CS%id_Kd_Niku_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + endif + if (CS%Lowmode_itidal_dissipation) then + vd = var_desc("Kd_lowmode","m2 s-1", & + "Internal Tide Driven Diffusivity (from low modes), interpolated to z",& + z_grid='z') + CS%id_Kd_lowmode_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + endif + endif + + endif + +end function tidal_mixing_init + + +!> Depending on whether or not CVMix is active, calls the associated subroutine to compute internal +!! tidal dissipation and to add the effect of internal-tide-driven mixing to the layer or interface +!! diffusivities. +subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, & + N2_lay, N2_int, Kd, Kd_int, Kd_max) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) + real, dimension(SZI_(G)), intent(in) :: N2_bot + real, dimension(SZI_(G),SZK_(G)), intent(in) :: N2_lay + real, dimension(SZI_(G),SZK_(G)+1), intent(in) :: N2_int + integer, intent(in) :: j + real, dimension(SZI_(G),SZK_(G)), intent(in) :: TKE_to_Kd, max_TKE + type(tidal_mixing_cs), pointer :: CS + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int + real, intent(inout) :: Kd_max + + if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) then + if (CS%use_CVMix_tidal) then + call calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) + else + call add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, & + N2_lay, Kd, Kd_int, Kd_max) + endif + endif +end subroutine + + +!> Calls the CVMix routines to compute tidal dissipation and to add the effect of internal-tide-driven +!! mixing to the interface diffusivities. +subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) + integer, intent(in) :: j + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(tidal_mixing_cs), pointer :: CS !< This module's control structure. + real, dimension(SZI_(G),SZK_(G)+1), intent(in) :: N2_int + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd + + ! local + real, dimension(SZK_(G)+1) :: Kd_tidal !< tidal diffusivity [m2/s] + real, dimension(SZK_(G)+1) :: Kv_tidal !< tidal viscosity [m2/s] + real, dimension(SZK_(G)+1) :: vert_dep !< vertical deposition needed for Simmons tidal mixing. + real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) + real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) + integer :: i, k, is, ie + real :: dh, hcorr, Simmons_coeff + real, parameter :: rho_fw = 1000.0 ! fresh water density [kg/m^3] ! TODO: when coupled, get this from CESM (SHR_CONST_RHOFW) + type(tidal_mixing_diags), pointer :: dd + + is = G%isc ; ie = G%iec + dd => CS%dd + + select case (CS%int_tide_profile) + case (SIMMONS_04) + do i=is,ie + + if (G%mask2dT(i,j)<1) cycle + + iFaceHeight = 0.0 ! BBL is all relative to the surface + hcorr = 0.0 + do k=1,G%ke + ! cell center and cell bottom in meters (negative values in the ocean) + dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) + hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 + dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * dh + iFaceHeight(k+1) = iFaceHeight(k) - dh + enddo + + call CVMix_compute_Simmons_invariant( nlev = G%ke, & + energy_flux = CS%tidal_qe_2d(i,j), & + rho = rho_fw, & + SimmonsCoeff = Simmons_coeff, & + VertDep = vert_dep, & + zw = iFaceHeight, & + zt = cellHeight, & + CVMix_tidal_params_user = CS%CVMix_tidal_params) + + ! Since we pass tidal_qe_2d=(CS%Gamma_itides)*tidal_energy_flux_2d, and not tidal_energy_flux_2d in + ! above subroutine call, we divide Simmons_coeff by CS%Gamma_itides as a corrective step: + ! TODO: (CS%Gamma_itides)*tidal_energy_flux_2d is unnecessary, directly use tidal_energy_flux_2d + Simmons_coeff = Simmons_coeff / CS%Gamma_itides + + + call CVMix_coeffs_tidal( Mdiff_out = Kv_tidal, & + Tdiff_out = Kd_tidal, & + Nsqr = N2_int(i,:), & + OceanDepth = -iFaceHeight(G%ke+1),& + SimmonsCoeff = Simmons_coeff, & + vert_dep = vert_dep, & + nlev = G%ke, & + max_nlev = G%ke, & + CVMix_params = CS%CVMix_glb_params, & + CVMix_tidal_params_user = CS%CVMix_tidal_params) + + do k=1,G%ke + Kd(i,j,k) = Kd(i,j,k) + 0.5*(Kd_tidal(k) + Kd_tidal(k+1) ) + !TODO: Kv(i,j,k) = ???????????? + enddo + + ! diagnostics + if (associated(dd%Kd_itidal)) then + dd%Kd_itidal(i,j,:) = Kd_tidal(:) + endif + if (associated(dd%N2_int)) then + dd%N2_int(i,j,:) = N2_int(i,:) + endif + if (associated(dd%Simmons_coeff_2d)) then + dd%Simmons_coeff_2d(i,j) = Simmons_coeff + endif + if (associated(dd%vert_dep_3d)) then + dd%vert_dep_3d(i,j,:) = vert_dep(:) + endif + + enddo ! i=is,ie + + ! TODO: case (SCHMITTNER) + case default + call MOM_error(FATAL, "tidal_mixing_init: The selected"// & + " INT_TIDE_PROFILE is unavailable in CVMix") + end select + +end subroutine calculate_CVMix_tidal + + +!> This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities. +!! The mechanisms considered are (1) local dissipation of internal waves generated by the +!! barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating +!! low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves. +!! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction, +!! Froude-number-depending breaking, PSI, etc.). +subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, & + N2_lay, Kd, Kd_int, Kd_max) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) + real, dimension(SZI_(G)), intent(in) :: N2_bot + real, dimension(SZI_(G),SZK_(G)), intent(in) :: N2_lay + integer, intent(in) :: j + real, dimension(SZI_(G),SZK_(G)), intent(in) :: TKE_to_Kd, max_TKE + type(tidal_mixing_cs), pointer :: CS + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int + real, intent(inout) :: Kd_max + + ! This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities. + ! The mechanisms considered are (1) local dissipation of internal waves generated by the + ! barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating + ! low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves. + ! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction, + ! Froude-number-depending breaking, PSI, etc.). + + real, dimension(SZI_(G)) :: & + htot, & ! total thickness above or below a layer, or the + ! integrated thickness in the BBL (meter) + htot_WKB, & ! distance from top to bottom (meter) WKB scaled + TKE_itidal_bot, & ! internal tide TKE at ocean bottom (m3/s3) + TKE_Niku_bot, & ! lee-wave TKE at ocean bottom (m3/s3) + TKE_lowmode_bot, & ! internal tide TKE at ocean bottom lost from all remote low modes (m3/s3) (BDM) + Inv_int, & ! inverse of TKE decay for int tide over the depth of the ocean (nondim) + Inv_int_lee, & ! inverse of TKE decay for lee waves over the depth of the ocean (nondim) + Inv_int_low, & ! inverse of TKE decay for low modes over the depth of the ocean (nondim) (BDM) + z0_Polzin, & ! TKE decay scale in Polzin formulation (meter) + z0_Polzin_scaled, & ! TKE decay scale in Polzin formulation (meter) + ! multiplied by N2_bot/N2_meanz to be coherent with the WKB scaled z + ! z*=int(N2/N2_bot) * N2_bot/N2_meanz = int(N2/N2_meanz) + ! z0_Polzin_scaled = z0_Polzin * N2_bot/N2_meanz + N2_meanz, & ! vertically averaged squared buoyancy frequency (1/s2) for WKB scaling + TKE_itidal_rem, & ! remaining internal tide TKE (from barotropic source) + TKE_Niku_rem, & ! remaining lee-wave TKE + TKE_lowmode_rem, & ! remaining internal tide TKE (from propagating low mode source) (BDM) + TKE_frac_top, & ! fraction of bottom TKE that should appear at top of a layer (nondim) + TKE_frac_top_lee, & ! fraction of bottom TKE that should appear at top of a layer (nondim) + TKE_frac_top_lowmode, & + ! fraction of bottom TKE that should appear at top of a layer (nondim) (BDM) + z_from_bot, & ! distance from bottom (meter) + z_from_bot_WKB ! distance from bottom (meter), WKB scaled + + real :: I_rho0 ! 1 / RHO0, (m3/kg) + real :: Kd_add ! diffusivity to add in a layer (m2/sec) + real :: TKE_itide_lay ! internal tide TKE imparted to a layer (from barotropic) (m3/s3) + real :: TKE_Niku_lay ! lee-wave TKE imparted to a layer (m3/s3) + real :: TKE_lowmode_lay ! internal tide TKE imparted to a layer (from low mode) (m3/s3) (BDM) + real :: frac_used ! fraction of TKE that can be used in a layer (nondim) + real :: Izeta ! inverse of TKE decay scale (1/meter) + real :: Izeta_lee ! inverse of TKE decay scale for lee waves (1/meter) + real :: z0_psl ! temporary variable with units of meter + real :: TKE_lowmode_tot ! TKE from all low modes (W/m2) (BDM) + + logical :: use_Polzin, use_Simmons + integer :: i, k, is, ie, nz + integer :: a, fr, m + type(tidal_mixing_diags), pointer :: dd + + is = G%isc ; ie = G%iec ; nz = G%ke + dd => CS%dd + + if (.not.(CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation)) return + + do i=is,ie ; htot(i) = 0.0 ; Inv_int(i) = 0.0 ; Inv_int_lee(i) = 0.0 ; Inv_int_low(i) = 0.0 ;enddo + do k=1,nz ; do i=is,ie + htot(i) = htot(i) + GV%H_to_m*h(i,j,k) + enddo ; enddo + + I_Rho0 = 1.0/GV%Rho0 + + use_Polzin = ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09)) .or. & + (CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09)) .or. & + (CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == POLZIN_09))) + use_Simmons = ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == STLAURENT_02)) .or. & + (CS%lee_wave_dissipation .and. (CS%lee_wave_profile == STLAURENT_02)) .or. & + (CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == STLAURENT_02))) + + ! Calculate parameters for vertical structure of dissipation + ! Simmons: + if ( use_Simmons ) then + Izeta = 1.0 / max(CS%Int_tide_decay_scale, GV%H_subroundoff*GV%H_to_m) + Izeta_lee = 1.0 / max(CS%Int_tide_decay_scale*CS%Decay_scale_factor_lee, & + GV%H_subroundoff*GV%H_to_m) + do i=is,ie + CS%Nb(i,j) = sqrt(N2_bot(i)) + if (associated(dd%N2_bot)) dd%N2_bot(i,j) = N2_bot(i) + if ( CS%Int_tide_dissipation ) then + if (Izeta*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule. + Inv_int(i) = 1.0 / (1.0 - exp(-Izeta*htot(i))) + endif + endif + if ( CS%Lee_wave_dissipation ) then + if (Izeta_lee*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule. + Inv_int_lee(i) = 1.0 / (1.0 - exp(-Izeta_lee*htot(i))) + endif + endif + if ( CS%Lowmode_itidal_dissipation) then + if (Izeta*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule. + Inv_int_low(i) = 1.0 / (1.0 - exp(-Izeta*htot(i))) + endif + endif + z_from_bot(i) = GV%H_to_m*h(i,j,nz) + enddo + endif ! Simmons + + ! Polzin: + if ( use_Polzin ) then + ! WKB scaling of the vertical coordinate + do i=is,ie ; N2_meanz(i)=0.0 ; enddo + do k=1,nz ; do i=is,ie + N2_meanz(i) = N2_meanz(i) + N2_lay(i,k)*GV%H_to_m*h(i,j,k) + enddo ; enddo + do i=is,ie + N2_meanz(i) = N2_meanz(i) / (htot(i) + GV%H_subroundoff*GV%H_to_m) + if (associated(dd%N2_meanz)) dd%N2_meanz(i,j) = N2_meanz(i) + enddo + + ! WKB scaled z*(z=H) z* at the surface using the modified Polzin WKB scaling + do i=is,ie ; htot_WKB(i) = htot(i) ; enddo +! do i=is,ie ; htot_WKB(i) = 0.0 ; enddo +! do k=1,nz ; do i=is,ie +! htot_WKB(i) = htot_WKB(i) + GV%H_to_m*h(i,j,k)*N2_lay(i,k) / N2_meanz(i) +! enddo ; enddo + ! htot_WKB(i) = htot(i) ! Nearly equivalent and simpler + + do i=is,ie + CS%Nb(i,j) = sqrt(N2_bot(i)) + if ((CS%tideamp(i,j) > 0.0) .and. & + (CS%kappa_itides**2 * CS%h2(i,j) * CS%Nb(i,j)**3 > 1.0e-14) ) then + z0_polzin(i) = CS%Polzin_decay_scale_factor * CS%Nu_Polzin * & + CS%Nbotref_Polzin**2 * CS%tideamp(i,j) / & + ( CS%kappa_itides**2 * CS%h2(i,j) * CS%Nb(i,j)**3 ) + if (z0_polzin(i) < CS%Polzin_min_decay_scale) & + z0_polzin(i) = CS%Polzin_min_decay_scale + if (N2_meanz(i) > 1.0e-14 ) then + z0_polzin_scaled(i) = z0_polzin(i)*CS%Nb(i,j)**2 / N2_meanz(i) + else + z0_polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * htot(i) + endif + if (z0_polzin_scaled(i) > (CS%Polzin_decay_scale_max_factor * htot(i)) ) & + z0_polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * htot(i) + else + z0_polzin(i) = CS%Polzin_decay_scale_max_factor * htot(i) + z0_polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * htot(i) + endif + + if (associated(dd%Polzin_decay_scale)) & + dd%Polzin_decay_scale(i,j) = z0_polzin(i) + if (associated(dd%Polzin_decay_scale_scaled)) & + dd%Polzin_decay_scale_scaled(i,j) = z0_polzin_scaled(i) + if (associated(dd%N2_bot)) dd%N2_bot(i,j) = CS%Nb(i,j)*CS%Nb(i,j) + + if ( CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09) ) then + ! For the Polzin formulation, this if loop prevents the vertical + ! flux of energy dissipation from having NaN values + if (htot_WKB(i) > 1.0e-14) then + Inv_int(i) = ( z0_polzin_scaled(i) / htot_WKB(i) ) + 1 + endif + endif + if ( CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09) ) then + ! For the Polzin formulation, this if loop prevents the vertical + ! flux of energy dissipation from having NaN values + if (htot_WKB(i) > 1.0e-14) then + Inv_int_lee(i) = ( z0_polzin_scaled(i)*CS%Decay_scale_factor_lee / htot_WKB(i) ) + 1 + endif + endif + if ( CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == POLZIN_09) ) then + ! For the Polzin formulation, this if loop prevents the vertical + ! flux of energy dissipation from having NaN values + if (htot_WKB(i) > 1.0e-14) then + Inv_int_low(i) = ( z0_polzin_scaled(i) / htot_WKB(i) ) + 1 + endif + endif + + z_from_bot(i) = GV%H_to_m*h(i,j,nz) + ! Use the new formulation for WKB scaling. N2 is referenced to its + ! vertical mean. + if (N2_meanz(i) > 1.0e-14 ) then + z_from_bot_WKB(i) = GV%H_to_m*h(i,j,nz)*N2_lay(i,nz) / N2_meanz(i) + else ; z_from_bot_WKB(i) = 0 ; endif + enddo + endif ! Polzin + + ! Calculate/get dissipation values at bottom + ! Both Polzin and Simmons: + do i=is,ie + ! Dissipation of locally trapped internal tide (non-propagating high modes) + TKE_itidal_bot(i) = min(CS%TKE_itidal(i,j)*CS%Nb(i,j),CS%TKE_itide_max) + if (associated(dd%TKE_itidal_used)) & + dd%TKE_itidal_used(i,j) = TKE_itidal_bot(i) + TKE_itidal_bot(i) = (I_rho0 * CS%Mu_itides * CS%Gamma_itides) * TKE_itidal_bot(i) + ! Dissipation of locally trapped lee waves + TKE_Niku_bot(i) = 0.0 + if (CS%Lee_wave_dissipation) then + TKE_Niku_bot(i) = (I_rho0 * CS%Mu_itides * CS%Gamma_lee) * CS%TKE_Niku(i,j) + endif + ! Dissipation of propagating internal tide (baroclinic low modes; rays) (BDM) + TKE_lowmode_tot = 0.0 + TKE_lowmode_bot(i) = 0.0 + if (CS%Lowmode_itidal_dissipation) then + ! get loss rate due to wave drag on low modes (already multiplied by q) + + ! TODO: uncomment the following call and fix it + !call get_lowmode_loss(i,j,G,CS%int_tide_CSp,"WaveDrag",TKE_lowmode_tot) + print *, "========", __FILE__, __LINE__ + call MOM_error(FATAL,"this block not supported yet. (aa)") + + TKE_lowmode_bot(i) = CS%Mu_itides * I_rho0 * TKE_lowmode_tot + endif + ! Vertical energy flux at bottom + TKE_itidal_rem(i) = Inv_int(i) * TKE_itidal_bot(i) + TKE_Niku_rem(i) = Inv_int_lee(i) * TKE_Niku_bot(i) + TKE_lowmode_rem(i) = Inv_int_low(i) * TKE_lowmode_bot(i) + + if (associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,nz) = TKE_itidal_rem(i) !why is this here? BDM + enddo + + ! Estimate the work that would be done by mixing in each layer. + ! Simmons: + if ( use_Simmons ) then + do k=nz-1,2,-1 ; do i=is,ie + if (max_TKE(i,k) <= 0.0) cycle + z_from_bot(i) = z_from_bot(i) + GV%H_to_m*h(i,j,k) + + ! Fraction of bottom flux predicted to reach top of this layer + TKE_frac_top(i) = Inv_int(i) * exp(-Izeta * z_from_bot(i)) + TKE_frac_top_lee(i) = Inv_int_lee(i) * exp(-Izeta_lee * z_from_bot(i)) + TKE_frac_top_lowmode(i) = Inv_int_low(i) * exp(-Izeta * z_from_bot(i)) + + ! Actual influx at bottom of layer minus predicted outflux at top of layer to give + ! predicted power expended + TKE_itide_lay = TKE_itidal_rem(i) - TKE_itidal_bot(i) * TKE_frac_top(i) + TKE_Niku_lay = TKE_Niku_rem(i) - TKE_Niku_bot(i) * TKE_frac_top_lee(i) + TKE_lowmode_lay = TKE_lowmode_rem(i) - TKE_lowmode_bot(i)* TKE_frac_top_lowmode(i) + + ! Actual power expended may be less than predicted if stratification is weak; adjust + if (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay > max_TKE(i,k)) then + frac_used = max_TKE(i,k) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) + TKE_itide_lay = frac_used * TKE_itide_lay + TKE_Niku_lay = frac_used * TKE_Niku_lay + TKE_lowmode_lay = frac_used * TKE_lowmode_lay + endif + + ! Calculate vertical flux available to bottom of layer above + TKE_itidal_rem(i) = TKE_itidal_rem(i) - TKE_itide_lay + TKE_Niku_rem(i) = TKE_Niku_rem(i) - TKE_Niku_lay + TKE_lowmode_rem(i) = TKE_lowmode_rem(i) - TKE_lowmode_lay + + ! Convert power to diffusivity + Kd_add = TKE_to_Kd(i,k) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) + + if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max) + Kd(i,j,k) = Kd(i,j,k) + Kd_add + + if (present(Kd_int)) then + Kd_int(i,j,K) = Kd_int(i,j,K) + 0.5*Kd_add + Kd_int(i,j,K+1) = Kd_int(i,j,K+1) + 0.5*Kd_add + endif + + ! diagnostics + if (associated(dd%Kd_itidal)) then + ! If at layers, dd%Kd_itidal is just TKE_to_Kd(i,k) * TKE_itide_lay + ! The following sets the interface diagnostics. + Kd_add = TKE_to_Kd(i,k) * TKE_itide_lay + if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max) + if (k>1) dd%Kd_itidal(i,j,K) = dd%Kd_itidal(i,j,K) + 0.5*Kd_add + if (k= 0.0) Kd_add = min(Kd_add, Kd_max) + if (k>1) dd%Kd_Niku(i,j,K) = dd%Kd_Niku(i,j,K) + 0.5*Kd_add + if (k= 0.0) Kd_add = min(Kd_add, Kd_max) + if (k>1) dd%Kd_lowmode(i,j,K) = dd%Kd_lowmode(i,j,K) + 0.5*Kd_add + if (k 1.0e-14 ) then + z_from_bot_WKB(i) = z_from_bot_WKB(i) + GV%H_to_m*h(i,j,k)*N2_lay(i,k)/N2_meanz(i) + else ; z_from_bot_WKB(i) = 0 ; endif + + ! Fraction of bottom flux predicted to reach top of this layer + TKE_frac_top(i) = ( Inv_int(i) * z0_polzin_scaled(i) ) / & + ( z0_polzin_scaled(i) + z_from_bot_WKB(i) ) + z0_psl = z0_polzin_scaled(i)*CS%Decay_scale_factor_lee + TKE_frac_top_lee(i) = (Inv_int_lee(i) * z0_psl) / (z0_psl + z_from_bot_WKB(i)) + TKE_frac_top_lowmode(i) = ( Inv_int_low(i) * z0_polzin_scaled(i) ) / & + ( z0_polzin_scaled(i) + z_from_bot_WKB(i) ) + + ! Actual influx at bottom of layer minus predicted outflux at top of layer to give + ! predicted power expended + TKE_itide_lay = TKE_itidal_rem(i) - TKE_itidal_bot(i) *TKE_frac_top(i) + TKE_Niku_lay = TKE_Niku_rem(i) - TKE_Niku_bot(i) * TKE_frac_top_lee(i) + TKE_lowmode_lay = TKE_lowmode_rem(i) - TKE_lowmode_bot(i)*TKE_frac_top_lowmode(i) + + ! Actual power expended may be less than predicted if stratification is weak; adjust + if (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay > max_TKE(i,k)) then + frac_used = max_TKE(i,k) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) + TKE_itide_lay = frac_used * TKE_itide_lay + TKE_Niku_lay = frac_used * TKE_Niku_lay + TKE_lowmode_lay = frac_used * TKE_lowmode_lay + endif + + ! Calculate vertical flux available to bottom of layer above + TKE_itidal_rem(i) = TKE_itidal_rem(i) - TKE_itide_lay + TKE_Niku_rem(i) = TKE_Niku_rem(i) - TKE_Niku_lay + TKE_lowmode_rem(i) = TKE_lowmode_rem(i) - TKE_lowmode_lay + + ! Convert power to diffusivity + Kd_add = TKE_to_Kd(i,k) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) + + if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max) + Kd(i,j,k) = Kd(i,j,k) + Kd_add + + if (present(Kd_int)) then + Kd_int(i,j,K) = Kd_int(i,j,K) + 0.5*Kd_add + Kd_int(i,j,K+1) = Kd_int(i,j,K+1) + 0.5*Kd_add + endif + + ! diagnostics + if (associated(dd%Kd_itidal)) then + ! If at layers, this is just dd%Kd_itidal(i,j,K) = TKE_to_Kd(i,k) * TKE_itide_lay + ! The following sets the interface diagnostics. + Kd_add = TKE_to_Kd(i,k) * TKE_itide_lay + if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max) + if (k>1) dd%Kd_itidal(i,j,K) = dd%Kd_itidal(i,j,K) + 0.5*Kd_add + if (k= 0.0) Kd_add = min(Kd_add, Kd_max) + if (k>1) dd%Kd_Niku(i,j,K) = dd%Kd_Niku(i,j,K) + 0.5*Kd_add + if (k= 0.0) Kd_add = min(Kd_add, Kd_max) + if (k>1) dd%Kd_lowmode(i,j,K) = dd%Kd_lowmode(i,j,K) + 0.5*Kd_add + if (k Sets up diagnostics arrays for tidal mixing. +subroutine setup_tidal_diagnostics(G,CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(tidal_mixing_cs), pointer :: CS + + ! local + integer :: isd, ied, jsd, jed, nz + type(tidal_mixing_diags), pointer :: dd + + isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed; nz = G%ke + dd => CS%dd + + if ((CS%id_Kd_itidal > 0) .or. (CS%id_Kd_itidal_z > 0) .or. & + (CS%id_Kd_Itidal_work > 0)) then + allocate(dd%Kd_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Kd_itidal(:,:,:) = 0.0 + endif + if ((CS%id_Kd_lowmode > 0) .or. (CS%id_Kd_lowmode_z > 0) .or. & + (CS%id_Kd_lowmode_work > 0)) then + allocate(dd%Kd_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Kd_lowmode(:,:,:) = 0.0 + endif + if ( (CS%id_Fl_itidal > 0) ) then + allocate(dd%Fl_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Fl_itidal(:,:,:) = 0.0 + endif + if ( (CS%id_Fl_lowmode > 0) ) then + allocate(dd%Fl_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Fl_lowmode(:,:,:) = 0.0 + endif + if ( (CS%id_Polzin_decay_scale > 0) ) then + allocate(dd%Polzin_decay_scale(isd:ied,jsd:jed)) + dd%Polzin_decay_scale(:,:) = 0.0 + endif + if ( (CS%id_N2_bot > 0) ) then + allocate(dd%N2_bot(isd:ied,jsd:jed)) ; dd%N2_bot(:,:) = 0.0 + endif + if ( (CS%id_N2_meanz > 0) ) then + allocate(dd%N2_meanz(isd:ied,jsd:jed)) ; dd%N2_meanz(:,:) = 0.0 + endif + if ( (CS%id_Polzin_decay_scale_scaled > 0) ) then + allocate(dd%Polzin_decay_scale_scaled(isd:ied,jsd:jed)) + dd%Polzin_decay_scale_scaled(:,:) = 0.0 + endif + if ((CS%id_Kd_Niku > 0) .or. (CS%id_Kd_Niku_z > 0) .or. & + (CS%id_Kd_Niku_work > 0)) then + allocate(dd%Kd_Niku(isd:ied,jsd:jed,nz+1)) ; dd%Kd_Niku(:,:,:) = 0.0 + endif + if (CS%id_Kd_Niku_work > 0) then + allocate(dd%Kd_Niku_work(isd:ied,jsd:jed,nz)) ; dd%Kd_Niku_work(:,:,:) = 0.0 + endif + if (CS%id_Kd_Itidal_work > 0) then + allocate(dd%Kd_Itidal_work(isd:ied,jsd:jed,nz)) + dd%Kd_Itidal_work(:,:,:) = 0.0 + endif + if (CS%id_Kd_Lowmode_Work > 0) then + allocate(dd%Kd_Lowmode_Work(isd:ied,jsd:jed,nz)) + dd%Kd_Lowmode_Work(:,:,:) = 0.0 + endif + if (CS%id_TKE_itidal > 0) then + allocate(dd%TKE_Itidal_used(isd:ied,jsd:jed)) ; dd%TKE_Itidal_used(:,:) = 0. + endif + ! additional diags for CVMix + if (CS%id_N2_int > 0) then + allocate(dd%N2_int(isd:ied,jsd:jed,nz+1)) ; dd%N2_int(:,:,:) = 0.0 + endif + if (CS%id_Simmons_coeff > 0) then + allocate(dd%Simmons_coeff_2d(isd:ied,jsd:jed)) ; dd%Simmons_coeff_2d(:,:) = 0.0 + endif + if (CS%id_vert_dep > 0) then + allocate(dd%vert_dep_3d(isd:ied,jsd:jed,nz+1)) ; dd%vert_dep_3d(:,:,:) = 0.0 + endif +end subroutine setup_tidal_diagnostics + +subroutine post_tidal_diagnostics(G,GV,h,CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). + type(tidal_mixing_cs), pointer :: CS + + ! local + integer :: num_z_diags + integer :: z_ids(6) ! id numbers of diagns to be interpolated to depth space + type(p3d) :: z_ptrs(6) ! pointers to diagns to be interpolated into depth space + type(tidal_mixing_diags), pointer :: dd + + num_z_diags = 0 + dd => CS%dd + + if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) then + if (CS%id_TKE_itidal > 0) call post_data(CS%id_TKE_itidal, dd%TKE_itidal_used, CS%diag) + if (CS%id_TKE_leewave > 0) call post_data(CS%id_TKE_leewave, CS%TKE_Niku, CS%diag) + if (CS%id_Nb > 0) call post_data(CS%id_Nb, CS%Nb, CS%diag) + if (CS%id_N2_bot > 0) call post_data(CS%id_N2_bot, dd%N2_bot, CS%diag) + if (CS%id_N2_meanz > 0) call post_data(CS%id_N2_meanz,dd%N2_meanz,CS%diag) + + if (CS%id_Fl_itidal > 0) call post_data(CS%id_Fl_itidal, dd%Fl_itidal, CS%diag) + if (CS%id_Kd_itidal > 0) call post_data(CS%id_Kd_itidal, dd%Kd_itidal, CS%diag) + if (CS%id_Kd_Niku > 0) call post_data(CS%id_Kd_Niku, dd%Kd_Niku, CS%diag) + if (CS%id_Kd_lowmode> 0) call post_data(CS%id_Kd_lowmode, dd%Kd_lowmode, CS%diag) + if (CS%id_Fl_lowmode> 0) call post_data(CS%id_Fl_lowmode, dd%Fl_lowmode, CS%diag) + + if (CS%id_N2_int> 0) call post_data(CS%id_N2_int, dd%N2_int, CS%diag) + if (CS%id_vert_dep> 0) call post_data(CS%id_vert_dep, dd%vert_dep_3d, CS%diag) + if (CS%id_Simmons_coeff> 0) call post_data(CS%id_Simmons_coeff, dd%Simmons_coeff_2d, CS%diag) + + if (CS%id_Kd_Itidal_Work > 0) & + call post_data(CS%id_Kd_Itidal_Work, dd%Kd_Itidal_Work, CS%diag) + if (CS%id_Kd_Niku_Work > 0) call post_data(CS%id_Kd_Niku_Work, dd%Kd_Niku_Work, CS%diag) + if (CS%id_Kd_Lowmode_Work > 0) & + call post_data(CS%id_Kd_Lowmode_Work, dd%Kd_Lowmode_Work, CS%diag) + + if (CS%id_Polzin_decay_scale > 0 ) & + call post_data(CS%id_Polzin_decay_scale, dd%Polzin_decay_scale, CS%diag) + if (CS%id_Polzin_decay_scale_scaled > 0 ) & + call post_data(CS%id_Polzin_decay_scale_scaled, dd%Polzin_decay_scale_scaled, CS%diag) + + if (CS%id_Kd_itidal_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_Kd_itidal_z + z_ptrs(num_z_diags)%p => dd%Kd_itidal + endif + + if (CS%id_Kd_Niku_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_Kd_Niku_z + z_ptrs(num_z_diags)%p => dd%Kd_Niku + endif + + if (CS%id_Kd_lowmode_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_Kd_lowmode_z + z_ptrs(num_z_diags)%p => dd%Kd_lowmode + endif + + endif + + if (num_z_diags > 0) & + call calc_Zint_diags(h, z_ptrs, z_ids, num_z_diags, G, GV, CS%diag_to_Z_CSp) + + if (associated(dd%Kd_itidal)) deallocate(dd%Kd_itidal) + if (associated(dd%Kd_lowmode)) deallocate(dd%Kd_lowmode) + if (associated(dd%Fl_itidal)) deallocate(dd%Fl_itidal) + if (associated(dd%Fl_lowmode)) deallocate(dd%Fl_lowmode) + if (associated(dd%Polzin_decay_scale)) deallocate(dd%Polzin_decay_scale) + if (associated(dd%Polzin_decay_scale_scaled)) deallocate(dd%Polzin_decay_scale_scaled) + if (associated(dd%N2_bot)) deallocate(dd%N2_bot) + if (associated(dd%N2_meanz)) deallocate(dd%N2_meanz) + if (associated(dd%Kd_Niku)) deallocate(dd%Kd_Niku) + if (associated(dd%Kd_Niku_work)) deallocate(dd%Kd_Niku_work) + if (associated(dd%Kd_Itidal_Work)) deallocate(dd%Kd_Itidal_Work) + if (associated(dd%Kd_Lowmode_Work)) deallocate(dd%Kd_Lowmode_Work) + if (associated(dd%TKE_itidal_used)) deallocate(dd%TKE_itidal_used) + if (associated(dd%N2_int)) deallocate(dd%N2_int) + if (associated(dd%vert_dep_3d)) deallocate(dd%vert_dep_3d) + if (associated(dd%Simmons_coeff_2d)) deallocate(dd%Simmons_coeff_2d) + +end subroutine post_tidal_diagnostics + +! TODO: move this subroutine to MOM_internal_tide_input module (?) +subroutine read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + character(len=20), intent(in) :: tidal_energy_type + character(len=200), intent(in) :: tidal_energy_file + type(tidal_mixing_cs), pointer :: CS + ! local + integer :: isd, ied, jsd, jed + real, allocatable, dimension(:,:) :: tidal_energy_flux_2d ! input tidal energy flux at T-grid points (W/m^2) + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + + if (.not. allocated(CS%tidal_qe_2d)) allocate(CS%tidal_qe_2d(isd:ied,jsd:jed)) + allocate(tidal_energy_flux_2d(isd:ied,jsd:jed)) + + select case (uppercase(tidal_energy_type(1:4))) + case ('JAYN') ! Jayne 2009 input tidal energy flux + call MOM_read_data(tidal_energy_file,'wave_dissipation',tidal_energy_flux_2d, G%domain) + CS%tidal_qe_2d = (CS%Gamma_itides) * tidal_energy_flux_2d + case default + call MOM_error(FATAL, "read_tidal_energy: Unknown tidal energy file type.") + ! TODO: add more tidal energy file types, e.g., Arbic, ER03, GN13, LGM0, etc. + ! see POP::tidal_mixing.F90 + end select + + deallocate(tidal_energy_flux_2d) + +end subroutine read_tidal_energy + + +!> Clear pointers and deallocate memory +subroutine tidal_mixing_end(CS) + type(tidal_mixing_cs), pointer :: CS ! This module's control structure + + if (.not.associated(CS)) return + + !TODO deallocate all the dynamically allocated members here ... + if (allocated(CS%tidal_qe_2d)) deallocate(CS%tidal_qe_2d) + deallocate(CS%dd) + deallocate(CS) + +end subroutine tidal_mixing_end + + +end module MOM_tidal_mixing diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 2bb0b30206..3dfb82acd1 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1162,21 +1162,21 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m endif endif ; enddo - if (associated(visc%Kv_turb)) then - ! BGR/ Add factor of 2. * the averaged Kv_turb. + if (associated(visc%Kv_shear)) then + ! BGR/ Add factor of 2. * the averaged Kv_shear. ! this is needed to reproduce the analytical solution to ! a simple diffusion problem, likely due to h_shear being ! equal to 2 x \delta z if (work_on_u) then do K=2,nz ; do i=is,ie ; if (do_i(i)) then - Kv_add(i,K) = (2.*0.5)*(visc%Kv_turb(i,j,k) + visc%Kv_turb(i+1,j,k)) + Kv_add(i,K) = (2.*0.5)*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k)) endif ; enddo ; enddo if (do_OBCs) then do I=is,ie ; if (do_i(I) .and. (OBC%segnum_u(I,j) /= OBC_NONE)) then if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then - do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_turb(i,j,k) ; enddo + do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_shear(i,j,k) ; enddo elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then - do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_turb(i+1,j,k) ; enddo + do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_shear(i+1,j,k) ; enddo endif endif ; enddo endif @@ -1185,14 +1185,14 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m endif ; enddo ; enddo else do K=2,nz ; do i=is,ie ; if (do_i(i)) then - Kv_add(i,K) = (2.*0.5)*(visc%Kv_turb(i,j,k) + visc%Kv_turb(i,j+1,k)) + Kv_add(i,K) = (2.*0.5)*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i,j+1,k)) endif ; enddo ; enddo if (do_OBCs) then do i=is,ie ; if (do_i(i) .and. (OBC%segnum_v(i,J) /= OBC_NONE)) then if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then - do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_turb(i,j,k) ; enddo + do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_shear(i,j,k) ; enddo elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then - do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_turb(i,j+1,k) ; enddo + do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_shear(i,j+1,k) ; enddo endif endif ; enddo endif diff --git a/src/user/SCM_CVmix_tests.F90 b/src/user/SCM_CVmix_tests.F90 index 01ee725626..8795dab494 100644 --- a/src/user/SCM_CVmix_tests.F90 +++ b/src/user/SCM_CVmix_tests.F90 @@ -1,6 +1,6 @@ -!> Initial conditions and forcing for the single column model (SCM) CVmix +!> Initial conditions and forcing for the single column model (SCM) CVMix !! test set. -module SCM_CVmix_tests +module SCM_CVMix_tests ! This file is part of MOM6. See LICENSE.md for the license. @@ -18,14 +18,14 @@ module SCM_CVmix_tests #include -public SCM_CVmix_tests_TS_init -public SCM_CVmix_tests_surface_forcing_init -public SCM_CVmix_tests_wind_forcing -public SCM_CVmix_tests_buoyancy_forcing -public SCM_CVmix_tests_CS +public SCM_CVMix_tests_TS_init +public SCM_CVMix_tests_surface_forcing_init +public SCM_CVMix_tests_wind_forcing +public SCM_CVMix_tests_buoyancy_forcing +public SCM_CVMix_tests_CS !> Container for surface forcing parameters -type SCM_CVmix_tests_CS ; +type SCM_CVMix_tests_CS ; private logical :: UseWindStress !< True to use wind stress logical :: UseHeatFlux !< True to use heat flux @@ -42,12 +42,12 @@ module SCM_CVmix_tests ! This include declares and sets the variable "version". #include "version_variable.h" -character(len=40) :: mdl = "SCM_CVmix_tests" ! This module's name. +character(len=40) :: mdl = "SCM_CVMix_tests" ! This module's name. contains -!> Initializes temperature and salinity for the SCM CVmix test example -subroutine SCM_CVmix_tests_TS_init(T, S, h, G, GV, param_file, just_read_params) +!> Initializes temperature and salinity for the SCM CVMix test example +subroutine SCM_CVMix_tests_TS_init(T, S, h, G, GV, param_file, just_read_params) real, dimension(NIMEM_,NJMEM_, NKMEM_), intent(out) :: T !< Potential temperature (degC) real, dimension(NIMEM_,NJMEM_, NKMEM_), intent(out) :: S !< Salinity (psu) real, dimension(NIMEM_,NJMEM_, NKMEM_), intent(in) :: h !< Layer thickness (m or Pa) @@ -118,20 +118,20 @@ subroutine SCM_CVmix_tests_TS_init(T, S, h, G, GV, param_file, just_read_params) enddo ! k enddo ; enddo -end subroutine SCM_CVmix_tests_TS_init +end subroutine SCM_CVMix_tests_TS_init -!> Initializes surface forcing for the CVmix test case suite -subroutine SCM_CVmix_tests_surface_forcing_init(Time, G, param_file, CS) +!> Initializes surface forcing for the CVMix test case suite +subroutine SCM_CVMix_tests_surface_forcing_init(Time, G, param_file, CS) type(time_type), intent(in) :: Time !< Time type(ocean_grid_type), intent(in) :: G !< Grid structure type(param_file_type), intent(in) :: param_file !< Input parameter structure - type(SCM_CVmix_tests_CS), pointer :: CS !< Parameter container + type(SCM_CVMix_tests_CS), pointer :: CS !< Parameter container ! This include declares and sets the variable "version". #include "version_variable.h" if (associated(CS)) then - call MOM_error(FATAL, "SCM_CVmix_tests_surface_forcing_init called with an associated "// & + call MOM_error(FATAL, "SCM_CVMix_tests_surface_forcing_init called with an associated "// & "control structure.") return endif @@ -141,57 +141,57 @@ subroutine SCM_CVmix_tests_surface_forcing_init(Time, G, param_file, CS) call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "SCM_USE_WIND_STRESS", & CS%UseWindStress, "Wind Stress switch "// & - "used in the SCM CVmix surface forcing.", & + "used in the SCM CVMix surface forcing.", & units='', default=.false.) call get_param(param_file, mdl, "SCM_USE_HEAT_FLUX", & CS%UseHeatFlux, "Heat flux switch "// & - "used in the SCM CVmix test surface forcing.", & + "used in the SCM CVMix test surface forcing.", & units='', default=.false.) call get_param(param_file, mdl, "SCM_USE_EVAPORATION", & CS%UseEvaporation, "Evaporation switch "// & - "used in the SCM CVmix test surface forcing.", & + "used in the SCM CVMix test surface forcing.", & units='', default=.false.) call get_param(param_file, mdl, "SCM_USE_DIURNAL_SW", & CS%UseDiurnalSW, "Diurnal sw radation switch "// & - "used in the SCM CVmix test surface forcing.", & + "used in the SCM CVMix test surface forcing.", & units='', default=.false.) if (CS%UseWindStress) then call get_param(param_file, mdl, "SCM_TAU_X", & CS%tau_x, "Constant X-dir wind stress "// & - "used in the SCM CVmix test surface forcing.", & + "used in the SCM CVMix test surface forcing.", & units='N/m2', fail_if_missing=.true.) call get_param(param_file, mdl, "SCM_TAU_Y", & CS%tau_y, "Constant y-dir wind stress "// & - "used in the SCM CVmix test surface forcing.", & + "used in the SCM CVMix test surface forcing.", & units='N/m2', fail_if_missing=.true.) endif if (CS%UseHeatFlux) then call get_param(param_file, mdl, "SCM_HEAT_FLUX", & CS%surf_HF, "Constant surface heat flux "// & - "used in the SCM CVmix test surface forcing.", & + "used in the SCM CVMix test surface forcing.", & units='m K/s', fail_if_missing=.true.) endif if (CS%UseEvaporation) then call get_param(param_file, mdl, "SCM_EVAPORATION", & CS%surf_evap, "Constant surface evaporation "// & - "used in the SCM CVmix test surface forcing.", & + "used in the SCM CVMix test surface forcing.", & units='m/s', fail_if_missing=.true.) endif if (CS%UseDiurnalSW) then call get_param(param_file, mdl, "SCM_DIURNAL_SW_MAX", & CS%Max_sw, "Maximum diurnal sw radiation "// & - "used in the SCM CVmix test surface forcing.", & + "used in the SCM CVMix test surface forcing.", & units='m K/s', fail_if_missing=.true.) endif -end subroutine SCM_CVmix_tests_surface_forcing_init +end subroutine SCM_CVMix_tests_surface_forcing_init -subroutine SCM_CVmix_tests_wind_forcing(state, forces, day, G, CS) +subroutine SCM_CVMix_tests_wind_forcing(state, forces, day, G, CS) type(surface), intent(in) :: state !< Surface state structure type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces type(time_type), intent(in) :: day !< Time in days type(ocean_grid_type), intent(inout) :: G !< Grid structure - type(SCM_CVmix_tests_CS), pointer :: CS !< Container for SCM parameters + type(SCM_CVMix_tests_CS), pointer :: CS !< Container for SCM parameters ! Local variables integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -216,16 +216,15 @@ subroutine SCM_CVmix_tests_wind_forcing(state, forces, day, G, CS) forces%ustar(i,j) = sqrt( mag_tau / CS%Rho0 ) enddo ; enddo ; endif +end subroutine SCM_CVMix_tests_wind_forcing -end subroutine SCM_CVmix_tests_wind_forcing - -subroutine SCM_CVmix_tests_buoyancy_forcing(state, fluxes, day, G, CS) +subroutine SCM_CVMix_tests_buoyancy_forcing(state, fluxes, day, G, CS) type(surface), intent(in) :: state !< Surface state structure type(forcing), intent(inout) :: fluxes !< Surface fluxes structure type(time_type), intent(in) :: day !< Time in days (seconds?) type(ocean_grid_type), intent(inout) :: G !< Grid structure - type(SCM_CVmix_tests_CS), pointer :: CS !< Container for SCM parameters + type(SCM_CVMix_tests_CS), pointer :: CS !< Container for SCM parameters ! Local variables integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq @@ -241,7 +240,7 @@ subroutine SCM_CVmix_tests_buoyancy_forcing(state, fluxes, day, G, CS) IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB if (CS%UseHeatFlux) then - ! Note CVmix test inputs give Heat flux in [m K/s] + ! Note CVMix test inputs give Heat flux in [m K/s] ! therefore must convert to W/m2 by multiplying ! by Rho0*Cp do J=Jsq,Jeq ; do i=is,ie @@ -251,7 +250,7 @@ subroutine SCM_CVmix_tests_buoyancy_forcing(state, fluxes, day, G, CS) if (CS%UseEvaporation) then do J=Jsq,Jeq ; do i=is,ie - ! Note CVmix test inputs give evaporation in m/s + ! Note CVMix test inputs give evaporation in m/s ! This therefore must be converted to mass flux ! by multiplying by density fluxes%evap(i,J) = CS%surf_evap * CS%Rho0 @@ -260,7 +259,7 @@ subroutine SCM_CVmix_tests_buoyancy_forcing(state, fluxes, day, G, CS) if (CS%UseDiurnalSW) then do J=Jsq,Jeq ; do i=is,ie - ! Note CVmix test inputs give max sw rad in [m K/s] + ! Note CVMix test inputs give max sw rad in [m K/s] ! therefore must convert to W/m2 by multiplying ! by Rho0*Cp ! Note diurnal cycle peaks at Noon. @@ -269,6 +268,6 @@ subroutine SCM_CVmix_tests_buoyancy_forcing(state, fluxes, day, G, CS) enddo; enddo endif -end subroutine SCM_CVmix_tests_buoyancy_forcing +end subroutine SCM_CVMix_tests_buoyancy_forcing -end module SCM_CVmix_tests +end module SCM_CVMix_tests