diff --git a/pkg/CVMix-src b/pkg/CVMix-src index d83f582714..c3a711e4e4 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit d83f582714e7f0f98d20efd8fac8fab01fa3bfe6 +Subproject commit c3a711e4e45f5ebcdc528f8ac690ad9a843375e9 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 60ce1d9c99..22dbb86b15 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1448,7 +1448,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. @@ -2288,8 +2289,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) @@ -2900,6 +2901,8 @@ subroutine MOM_end(CS) call tracer_registry_end(CS%tracer_Reg) call tracer_flow_control_end(CS%tracer_flow_CSp) + 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/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index df64bb7b21..87ce532a28 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -32,6 +32,7 @@ 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 @@ -936,7 +937,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, & diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 new file mode 100644 index 0000000000..14c6c3412e --- /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_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 new file mode 100644 index 0000000000..55e7d55d6e --- /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', 'conv_N2', 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', 'conv_kd', diag%axesTi, Time, & + 'Additional diffusivity added by MOM_cvmix_conv module', 'm2/s') + CS%id_kv_conv = register_diag_field('ocean_model', 'conv_kv', 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,NINT(hbl(i,j)) + 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 index 460dde7c47..345522126b 100644 --- a/src/parameterizations/vertical/MOM_cvmix_shear.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_shear.F90 @@ -25,37 +25,44 @@ module MOM_cvmix_shear #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 -character(len=40) :: mdl = "MOM_CVMix_shear" !< This module's name. +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,18 +116,31 @@ 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. @@ -133,7 +153,7 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) 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 @@ -193,9 +213,29 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) 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. + + ! 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 @@ -213,4 +253,15 @@ logical function cvmix_shear_is_used(param_file) 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 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 91b3c343e0..25c2464b56 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -9,7 +9,7 @@ module MOM_diabatic_driver use MOM_checksum_packages, only : MOM_state_chksum, MOM_state_stats 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_CVMix_shear, only : CVMix_shear_is_used +use MOM_CVMix_shear, only : cvmix_shear_is_used use MOM_diabatic_aux, only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_CS use MOM_diabatic_aux, only : make_frazil, adjust_salt, insert_brine, differential_diffuse_T_S, triDiagTS use MOM_diabatic_aux, only : find_uv_at_h, diagnoseMLDbyDensityDifference, applyBoundaryFluxesInOut @@ -22,8 +22,8 @@ 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_energetic_PBL, only : energetic_PBL, energetic_PBL_init @@ -47,7 +47,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 @@ -90,6 +90,8 @@ module MOM_diabatic_driver !! shear-driven diapycnal diffusivity. logical :: use_cvmix_shear !< If true, use the CVMix module to find the !! shear-driven diapycnal 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 @@ -149,8 +151,6 @@ module MOM_diabatic_driver logical :: useKPP !< 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 @@ -220,7 +220,7 @@ 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(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 @@ -522,7 +522,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G if (CS%debug) then call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, haloshift=0) endif - if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then + if (CS%use_kappa_shear .or. CS%use_cvmix_shear) then if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then call find_uv_at_h(u, v, h_orig, u_h, v_h, G, GV, eaml, ebml) if (CS%debug) then @@ -582,8 +582,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)") @@ -632,9 +632,14 @@ 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, CS%KPP_NLTscalar) + fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar) !$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 @@ -666,9 +671,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 @@ -813,10 +827,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) @@ -1349,9 +1363,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 @@ -1905,7 +1919,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"//& @@ -2324,9 +2338,9 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, endif - ! 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) @@ -2414,7 +2428,9 @@ subroutine diabatic_driver_end(CS) deallocate( CS%KPP_NLTscalar ) call KPP_end(CS%KPP_CSp) endif - if (CS%useConvection) call diffConvection_end(CS%Conv_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) & 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 ce32e25b49..6794d7b45b 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 15cbd6bb1c..eb12eb23d6 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -2,28 +2,6 @@ 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 @@ -41,12 +19,14 @@ module MOM_set_diffusivity 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 +50,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,21 +75,6 @@ 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) @@ -244,11 +171,9 @@ 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 @@ -266,7 +191,8 @@ module MOM_set_diffusivity 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() integer :: id_TKE_itidal = -1 @@ -380,27 +306,9 @@ 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 @@ -421,22 +329,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 +358,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) @@ -569,60 +461,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 +496,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 +536,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 @@ -799,21 +627,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,6 +673,7 @@ 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 @@ -1293,6 +1133,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 @@ -2282,6 +2124,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, endif ! Polzin end subroutine add_int_tide_diffusivity + !> This subroutine calculates several properties related to bottom !! boundary layer turbulence. subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, CS) @@ -2533,6 +2376,7 @@ subroutine set_density_ratios(h, tv, kb, G, GV, CS, j, ds_dsp1, rho_0) end subroutine set_density_ratios +!> Initialized the set_diffusivity module subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp, int_tide_CSp) type(time_type), intent(in) :: Time type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. @@ -2547,21 +2391,14 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp type(int_tide_CS), pointer :: int_tide_CSp !< pointer to the internal tides control !! structure (BDM) -! 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) - + ! local variables real :: decay_length, utide, zbot, hamp type(vardesc) :: vd logical :: read_tideamp, 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 @@ -2706,79 +2543,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"//& @@ -3131,6 +2897,8 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp 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.", & @@ -3168,31 +2936,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%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..3df3e7b780 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -44,7 +44,8 @@ module MOM_set_visc use MOM_grid, only : ocean_grid_type 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_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,18 +1784,19 @@ 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"// & "to calculate diffusivities and non-local transport in the OBL.", & @@ -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_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index fee1fb456a..ff14a698ed 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1094,21 +1094,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 @@ -1117,14 +1117,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