diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9b70b81415..e96a3807a7 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2378,6 +2378,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (associated(CS%visc%Kv_shear)) & call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) + if (associated(CS%visc%Kv_slow)) & + call pass_var(CS%visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) + call cpu_clock_end(id_clock_pass_init) call register_obsolete_diagnostics(param_file, CS%diag) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 09305eb9fb..02b0b622a3 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -233,6 +233,9 @@ module MOM_variables !! convection etc). TKE_turb => NULL() !< The turbulent kinetic energy per unit mass defined !! at the interfaces between each layer, in m2 s-2. + logical :: add_Kv_slow !< If True, adds Kv_slow when calculating the + !! 'coupling coefficient' (a[k]) at the interfaces. + !! This is done in find_coupling_coef. end type vertvisc_type !> The BT_cont_type structure contains information about the summed layer diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index 2be8beee4a..57b86c80ca 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -212,6 +212,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, CS, hbl) iFaceHeight(k+1) = iFaceHeight(k) - dh enddo + ! gets index of the level and interface above hbl kOBL = CVMix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) call CVMix_coeffs_conv(Mdiff_out=CS%kv_conv(i,j,:), & diff --git a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 new file mode 100644 index 0000000000..7137aabfa6 --- /dev/null +++ b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 @@ -0,0 +1,301 @@ +!> Interface to CVMix double diffusion scheme. +module MOM_CVMix_ddiff + +! 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_derivs +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_ddiff, only : cvmix_init_ddiff, CVMix_coeffs_ddiff +use cvmix_kpp, only : CVmix_kpp_compute_kOBL_depth +implicit none ; private + +#include + +public CVMix_ddiff_init, CVMix_ddiff_end, CVMix_ddiff_is_used, compute_ddiff_coeffs + +!> Control structure including parameters for CVMix double diffusion. +type, public :: CVMix_ddiff_cs + + ! Parameters + real :: strat_param_max !< maximum value for the stratification parameter (nondim) + real :: kappa_ddiff_s !< leading coefficient in formula for salt-fingering regime + !! for salinity diffusion (m^2/s) + real :: ddiff_exp1 !< interior exponent in salt-fingering regime formula (nondim) + real :: ddiff_exp2 !< exterior exponent in salt-fingering regime formula (nondim) + real :: mol_diff !< molecular diffusivity (m^2/s) + real :: kappa_ddiff_param1 !< exterior coefficient in diffusive convection regime (nondim) + real :: kappa_ddiff_param2 !< middle coefficient in diffusive convection regime (nondim) + real :: kappa_ddiff_param3 !< interior coefficient in diffusive convection regime (nondim) + real :: min_thickness !< Minimum thickness allowed (m) + character(len=4) :: diff_conv_type !< type of diffusive convection to use. Options are Marmorino & + !! Caldwell 1976 ("MC76"; default) and Kelley 1988, 1990 ("K90") + logical :: debug !< If true, turn on debugging + + ! Daignostic handles and pointers + type(diag_ctrl), pointer :: diag => NULL() + integer :: id_KT_extra = -1, id_KS_extra = -1, id_R_rho = -1 + + ! Diagnostics arrays + real, allocatable, dimension(:,:,:) :: KT_extra !< double diffusion diffusivity for temp (m2/s) + real, allocatable, dimension(:,:,:) :: KS_extra !< double diffusion diffusivity for salt (m2/s) + real, allocatable, dimension(:,:,:) :: R_rho !< double-diffusion density ratio (nondim) + +end type CVMix_ddiff_cs + +character(len=40) :: mdl = "MOM_CVMix_ddiff" !< This module's name. + +contains + +!> Initialized the CVMix double diffusion module. +logical function CVMix_ddiff_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_ddiff_cs), pointer :: CS !< This module's control structure. + +! This include declares and sets the variable "version". +#include "version_variable.h" + + if (associated(CS)) then + call MOM_error(WARNING, "CVMix_ddiff_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + ! Read parameters + call log_version(param_file, mdl, version, & + "Parameterization of mixing due to double diffusion processes via CVMix") + call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_init, & + "If true, turns on double diffusive processes via CVMix. \n"// & + "Note that double diffusive processes on viscosity are ignored \n"// & + "in CVMix, see http://cvmix.github.io/ for justification.",& + default=.false.) + + if (.not. CVMix_ddiff_init) return + + 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_DDIFF') + + call get_param(param_file, mdl, "STRAT_PARAM_MAX", CS%strat_param_max, & + "The maximum value for the double dissusion stratification parameter", & + units="nondim", default=2.55) + + call get_param(param_file, mdl, "KAPPA_DDIFF_S", CS%kappa_ddiff_s, & + "Leading coefficient in formula for salt-fingering regime \n"// & + "for salinity diffusion.", units="m2 s-1", default=1.0e-4) + + call get_param(param_file, mdl, "DDIFF_EXP1", CS%ddiff_exp1, & + "Interior exponent in salt-fingering regime formula.", & + units="nondim", default=1.0) + + call get_param(param_file, mdl, "DDIFF_EXP2", CS%ddiff_exp2, & + "Exterior exponent in salt-fingering regime formula.", & + units="nondim", default=3.0) + + call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM1", CS%kappa_ddiff_param1, & + "Exterior coefficient in diffusive convection regime.", & + units="nondim", default=0.909) + + call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM2", CS%kappa_ddiff_param2, & + "Middle coefficient in diffusive convection regime.", & + units="nondim", default=4.6) + + call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM3", CS%kappa_ddiff_param3, & + "Interior coefficient in diffusive convection regime.", & + units="nondim", default=-0.54) + + call get_param(param_file, mdl, "MOL_DIFF", CS%mol_diff, & + "Molecular diffusivity used in CVMix double diffusion.", & + units="m2 s-1", default=1.5e-6) + + call get_param(param_file, mdl, "DIFF_CONV_TYPE", CS%diff_conv_type, & + "type of diffusive convection to use. Options are Marmorino \n" //& + "and Caldwell 1976 (MC76) and Kelley 1988, 1990 (K90).", & + default="MC76") + + call closeParameterBlock(param_file) + + ! Register diagnostics + CS%diag => diag + + CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for temperature', 'm2 s-1') + + CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for salinity', 'm2 s-1') + + CS%id_R_rho = register_diag_field('ocean_model','R_rho',diag%axesTi,Time, & + 'Double-diffusion density ratio', 'nondim') + if (CS%id_R_rho > 0) & + allocate(CS%R_rho( SZI_(G), SZJ_(G), SZK_(G)+1)); CS%R_rho(:,:,:) = 0.0 + + call cvmix_init_ddiff(strat_param_max=CS%strat_param_max, & + kappa_ddiff_s=CS%kappa_ddiff_s, & + ddiff_exp1=CS%ddiff_exp1, & + ddiff_exp2=CS%ddiff_exp2, & + mol_diff=CS%mol_diff, & + kappa_ddiff_param1=CS%kappa_ddiff_param1, & + kappa_ddiff_param2=CS%kappa_ddiff_param2, & + kappa_ddiff_param3=CS%kappa_ddiff_param3, & + diff_conv_type=CS%diff_conv_type) + +end function CVMix_ddiff_init + +!> Subroutine for computing vertical diffusion coefficients for the +!! double diffusion mixing parameterization. +subroutine compute_ddiff_coeffs(h, tv, G, GV, j, Kd_T, Kd_S, 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),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_T !< Interface double diffusion diapycnal + !! diffusivity for temp (m2/sec). + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_S !< Interface double diffusion diapycnal + !! diffusivity for salt (m2/sec). + type(CVMix_ddiff_cs), pointer :: CS !< The control structure returned + !! by a previous call to CVMix_ddiff_init. + integer, intent(in) :: j !< Meridional grid indice. +! real, dimension(:,:), optional, pointer :: hbl !< Depth of ocean boundary layer (m) + + ! local variables + real, dimension(SZK_(G)) :: & + cellHeight, & !< Height of cell centers (m) + dRho_dT, & !< partial derivatives of density wrt temp (kg m-3 degC-1) + dRho_dS, & !< partial derivatives of density wrt saln (kg m-3 PPT-1) + pres_int, & !< pressure at each interface (Pa) + temp_int, & !< temp and at interfaces (degC) + salt_int, & !< salt at at interfaces + alpha_dT, & !< alpha*dT across interfaces + beta_dS, & !< beta*dS across interfaces + dT, & !< temp. difference between adjacent layers (degC) + dS !< salt difference between adjacent layers + + real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) + integer :: kOBL !< level of OBL extent + real :: pref, g_o_rho0, rhok, rhokm1, dz, dh, hcorr + integer :: i, k + + ! initialize dummy variables + pres_int(:) = 0.0; temp_int(:) = 0.0; salt_int(:) = 0.0 + alpha_dT(:) = 0.0; beta_dS(:) = 0.0; dRho_dT(:) = 0.0 + dRho_dS(:) = 0.0; dT(:) = 0.0; dS(:) = 0.0 + + ! set Kd_T and Kd_S to zero to avoid passing values from previous call + Kd_T(:,j,:) = 0.0; Kd_S(:,j,:) = 0.0 + + ! GMM, I am leaving some code commented below. We need to pass BLD to + ! this soubroutine to avoid adding diffusivity above that. This needs + ! to be done once we re-structure the order of the calls. + !if (.not. associated(hbl)) then + ! allocate(hbl(SZI_(G), SZJ_(G))); + ! hbl(:,:) = 0.0 + !endif + + do i = G%isc, G%iec + + ! skip calling at land points + if (G%mask2dT(i,j) == 0.) cycle + + pRef = 0. + pres_int(1) = pRef + ! we don't have SST and SSS, so let's use values at top-most layer + temp_int(1) = TV%T(i,j,1); salt_int(1) = TV%S(i,j,1) + do k=2,G%ke + ! pressure at interface + pres_int(k) = pRef + GV%H_to_Pa * h(i,j,k-1) + ! temp and salt at interface + ! for temp: (t1*h1 + t2*h2)/(h1+h2) + temp_int(k) = (TV%T(i,j,k-1)*h(i,j,k-1) + TV%T(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k)) + salt_int(k) = (TV%S(i,j,k-1)*h(i,j,k-1) + TV%S(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k)) + ! dT and dS + dT(k) = (TV%T(i,j,k-1)-TV%T(i,j,k)) + dS(k) = (TV%S(i,j,k-1)-TV%S(i,j,k)) + pRef = pRef + GV%H_to_Pa * h(i,j,k-1) + enddo ! k-loop finishes + + call calculate_density_derivs(temp_int(:), salt_int(:), pres_int(:), drho_dT(:), drho_dS(:), 1, G%ke, TV%EQN_OF_STATE) + + ! The "-1.0" below is needed so that the following criteria is satisfied: + ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then "salt finger" + ! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then "diffusive convection" + do k=1,G%ke + alpha_dT(k) = -1.0*drho_dT(k) * dT(k) + beta_dS(k) = drho_dS(k) * dS(k) + enddo + + if (CS%id_R_rho > 0.0) then + do k=1,G%ke + CS%R_rho(i,j,k) = alpha_dT(k)/beta_dS(k) + ! avoid NaN's + if(CS%R_rho(i,j,k) /= CS%R_rho(i,j,k)) CS%R_rho(i,j,k) = 0.0 + enddo + endif + + 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 + + ! gets index of the level and interface above hbl + !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) + + call CVMix_coeffs_ddiff(Tdiff_out=Kd_T(i,j,:), & + Sdiff_out=Kd_S(i,j,:), & + strat_param_num=alpha_dT(:), & + strat_param_denom=beta_dS(:), & + nlev=G%ke, & + max_nlev=G%ke) + + ! Do not apply mixing due to convection within the boundary layer + !do k=1,kOBL + ! Kd_T(i,j,k) = 0.0 + ! Kd_S(i,j,k) = 0.0 + !enddo + + enddo ! i-loop + +end subroutine compute_ddiff_coeffs + +!> Reads the parameter "USE_CVMIX_DDIFF" 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_ddiff_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_DDIFF", CVMix_ddiff_is_used, & + default=.false., do_not_log = .true.) + +end function CVMix_ddiff_is_used + +!> Clear pointers and dealocate memory +subroutine CVMix_ddiff_end(CS) + type(CVMix_ddiff_cs), pointer :: CS ! Control structure + + deallocate(CS) + +end subroutine CVMix_ddiff_end + + +end module MOM_CVMix_ddiff diff --git a/src/parameterizations/vertical/MOM_CVMix_shear.F90 b/src/parameterizations/vertical/MOM_CVMix_shear.F90 index 2635af7fb5..53339d3488 100644 --- a/src/parameterizations/vertical/MOM_CVMix_shear.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_shear.F90 @@ -30,14 +30,14 @@ module MOM_CVMix_shear !> Control structure including parameters for CVMix interior shear schemes. type, public :: CVMix_shear_cs logical :: use_LMD94, use_PP81 !< Flags for various schemes + logical :: smooth_ri !< If true, smooth Ri using a 1-2-1 filter real :: Ri_zero !< LMD94 critical Richardson number real :: Nu_zero !< LMD94 maximum interior diffusivity - real :: KPP_exp !< + real :: KPP_exp !< Exponent of unitless factor of diff. + !! for KPP internal shear mixing scheme. 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) ! Daignostic handles and pointers type(diag_ctrl), pointer :: diag => NULL() @@ -52,25 +52,26 @@ module MOM_CVMix_shear !> 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. + 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) :: kd !< The vertical diffusivity at each interface - !! (not layer!) in m2 s-1. - 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 - !! CVMix_shear_init. + 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) :: kd !< The vertical diffusivity at each interface + !! (not layer!) in m2 s-1. + 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 + !! CVMix_shear_init. ! Local variables integer :: i, j, k, kk, km1 - real :: gorho - real :: pref, DU, DV, DRHO, DZ, N2, S2 + real :: GoRho + real :: pref, DU, DV, DRHO, DZ, N2, S2, dummy real, dimension(2*(G%ke)) :: pres_1d, temp_1d, salt_1d, rho_1d real, dimension(G%ke+1) :: Ri_Grad !< Gradient Richardson number - + real, parameter :: epsln = 1.e-10 !< Threshold to identify + !! vanished layers ! some constants GoRho = GV%g_Earth / GV%Rho0 @@ -120,10 +121,30 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, & ! 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 + Ri_grad(G%ke+1) = Ri_grad(G%ke) + + if (CS%smooth_ri) then + ! 1) fill Ri_grad in vanished layers with adjacent value + do k = 2, G%ke + if (h(i,j,k) .le. epsln) Ri_grad(k) = Ri_grad(k-1) + enddo + + Ri_grad(G%ke+1) = Ri_grad(G%ke) + + ! 2) vertically smooth Ri with 1-2-1 filter + dummy = 0.25 * Ri_grad(1) + Ri_grad(G%ke+1) = Ri_grad(G%ke) + do k = 1, G%ke + Ri_Grad(k) = dummy + 0.5 * Ri_Grad(k) + 0.25 * Ri_grad(k+1) + dummy = 0.25 * Ri_grad(k) + enddo + endif + + if (CS%id_ri_grad > 0) CS%ri_grad(i,j,:) = Ri_Grad(:) + ! Call to CVMix wrapper for computing interior mixing coefficients. call CVMix_coeffs_shear(Mdiff_out=kv(i,j,:), & Tdiff_out=kd(i,j,:), & @@ -209,7 +230,11 @@ 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 get_param(param_file, mdl, "SMOOTH_RI", CS%smooth_ri, & + "If true, vertically smooth the Richardson"// & + "number by applying a 1-2-1 filter once.", & + default = .false.) + 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) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index f98185685a..79234c7e11 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -1573,7 +1573,7 @@ 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) + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD!< bnd. layer depth (m) ! Local variables integer :: i,j do j = G%jsc, G%jec ; do i = G%isc, G%iec diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 61c212db8b..bb1e0b11c1 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -408,7 +408,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, kd_lay, kv, j, G, GV, CS) 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%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 diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 6eb3b854f4..528dc33135 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -2,53 +2,6 @@ module MOM_diabatic_aux ! This file is part of MOM6. See LICENSE.md for the license. -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, April 1994 - July 2000 * -!* Alistair Adcroft, and Stephen Griffies * -!* * -!* This program contains the subroutine that, along with the * -!* subroutines that it calls, implements diapycnal mass and momentum * -!* fluxes and a bulk mixed layer. The diapycnal diffusion can be * -!* used without the bulk mixed layer. * -!* * -!* diabatic first determines the (diffusive) diapycnal mass fluxes * -!* based on the convergence of the buoyancy fluxes within each layer. * -!* The dual-stream entrainment scheme of MacDougall and Dewar (JPO, * -!* 1997) is used for combined diapycnal advection and diffusion, * -!* calculated implicitly and potentially with the Richardson number * -!* dependent mixing, as described by Hallberg (MWR, 2000). Diapycnal * -!* advection is fundamentally the residual of diapycnal diffusion, * -!* so the fully implicit upwind differencing scheme that is used is * -!* entirely appropriate. The downward buoyancy flux in each layer * -!* is determined from an implicit calculation based on the previously * -!* calculated flux of the layer above and an estimated flux in the * -!* layer below. This flux is subject to the following conditions: * -!* (1) the flux in the top and bottom layers are set by the boundary * -!* conditions, and (2) no layer may be driven below an Angstrom thick-* -!* ness. If there is a bulk mixed layer, the buffer layer is treat- * -!* ed as a fixed density layer with vanishingly small diffusivity. * -!* * -!* diabatic takes 5 arguments: the two velocities (u and v), the * -!* thicknesses (h), a structure containing the forcing fields, and * -!* the length of time over which to act (dt). The velocities and * -!* thickness are taken as inputs and modified within the subroutine. * -!* There is no limit on the time step. * -!* * -!* 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, T, S, buoy, ustar, 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 : post_data, register_diag_field, safe_alloc_ptr @@ -239,26 +192,19 @@ subroutine make_frazil(h, tv, G, GV, CS, p_surf) end subroutine make_frazil +!> Applies double diffusion to T & S, assuming no diapycal mass +!! fluxes, using a simple triadiagonal solver. subroutine differential_diffuse_T_S(h, tv, visc, dt, G, GV) 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(thermo_var_ptrs), intent(inout) :: tv - type(vertvisc_type), intent(in) :: visc - real, intent(in) :: dt - -! This subroutine applies double diffusion to T & S, assuming no diapycal mass -! fluxes, using a simple triadiagonal solver. - -! Arguments: h - Layer thickness, in m or kg m-2. -! (in) tv - A structure containing pointers to any available -! thermodynamic fields. Absent fields have NULL ptrs. -! (in) visc - A structure containing vertical viscosities, bottom boundary -! layer properies, and related fields. -! (in) dt - Time increment, in s. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. + type(thermo_var_ptrs), intent(inout) :: tv !< pointers to any available modynamic fields. + !! Absent fields have NULL ptrs. + type(vertvisc_type), intent(in) :: visc !< structure containing vertical viscosities, + !! layer properies, and related fields. + real, intent(in) :: dt !< Time increment, in s. + ! local variables real, dimension(SZI_(G)) :: & b1_T, b1_S, & ! Variables used by the tridiagonal solvers of T & S, in H. d1_T, d1_S ! Variables used by the tridiagonal solvers, nondim. @@ -345,30 +291,25 @@ subroutine differential_diffuse_T_S(h, tv, visc, dt, G, GV) S(i,j,k) = S(i,j,k) + c1_S(i,k+1)*S(i,j,k+1) enddo ; enddo enddo - end subroutine differential_diffuse_T_S +!> Keep salinity from falling below a small but positive threshold +!! This occurs when the ice model attempts to extract more salt then +!! is actually available to it from the ocean. subroutine adjust_salt(h, tv, G, GV, 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(thermo_var_ptrs), intent(inout) :: tv - type(diabatic_aux_CS), intent(in) :: CS - -! Keep salinity from falling below a small but positive threshold -! This occurs when the ice model attempts to extract more salt then -! is actually available to it from the ocean. - -! Arguments: h - Layer thickness, in m. -! (in/out) tv - A structure containing pointers to any available -! thermodynamic fields. Absent fields have NULL ptrs. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) CS - The control structure returned by a previous call to -! diabatic_driver_init. - real :: salt_add_col(SZI_(G),SZJ_(G)) ! The accumulated salt requirement - real :: S_min ! The minimum salinity - real :: mc ! A layer's mass kg m-2 . + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m + !! or kg m-2) + type(thermo_var_ptrs), intent(inout) :: tv !< structure containing pointers to any + !! available thermodynamic fields. + type(diabatic_aux_CS), intent(in) :: CS !< control structure returned by + !! a previous call to diabatic_driver_init. + + ! local variables + real :: salt_add_col(SZI_(G),SZJ_(G)) !< The accumulated salt requirement + real :: S_min !< The minimum salinity + real :: mc !< A layer's mass kg m-2 . integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -410,33 +351,29 @@ subroutine adjust_salt(h, tv, G, GV, CS) end subroutine adjust_salt +!> Insert salt from brine rejection into the first layer below +!! the mixed layer which both contains mass and in which the +!! change in layer density remains stable after the addition +!! of salt via brine rejection. subroutine insert_brine(h, tv, G, GV, fluxes, nkmb, CS, dt, id_brine_lay) 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(thermo_var_ptrs), intent(inout) :: tv - type(forcing), intent(in) :: fluxes - integer, intent(in) :: nkmb - type(diabatic_aux_CS), intent(in) :: CS - real, intent(in) :: dt + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m + !! or kg m-2) + type(thermo_var_ptrs), intent(inout) :: tv !< structure containing pointers to + !! any available hermodynamic fields. + type(forcing), intent(in) :: fluxes !< tructure containing pointers + !! any possible forcing fields + integer, intent(in) :: nkmb !< number of layers in the mixed and + !! buffer layers + type(diabatic_aux_CS), intent(in) :: CS !< control structure returned by a + !! previous call to diabatic_driver_init. + real, intent(in) :: dt !< time step between calls to this + !! function (s) ?? integer, intent(in) :: id_brine_lay -! Insert salt from brine rejection into the first layer below -! the mixed layer which both contains mass and in which the -! change in layer density remains stable after the addition -! of salt via brine rejection. - -! Arguments: h - Layer thickness, in m. -! (in/out) tv - A structure containing pointers to any available -! thermodynamic fields. Absent fields have NULL ptrs. -! (in) fluxes = A structure containing pointers to any possible -! forcing fields; unused fields have NULL ptrs. -! (in) nkmb - The number of layers in the mixed and buffer layers. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) CS - The control structure returned by a previous call to -! diabatic_driver_init. + ! local variables real :: salt(SZI_(G)) ! The amount of salt rejected from ! sea ice. [grams] real :: dzbr(SZI_(G)) ! cumulative depth over which brine is distributed @@ -539,10 +476,9 @@ subroutine insert_brine(h, tv, G, GV, fluxes, nkmb, CS, dt, id_brine_lay) end subroutine insert_brine +!> Simple tri-diagnonal solver for T and S. +!! "Simple" means it only uses arrays hold, ea and eb. subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) -! Simple tri-diagnonal solver for T and S -! "Simple" means it only uses arrays hold, ea and eb - ! Arguments type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure integer, intent(in) :: is, ie, js, je @@ -579,37 +515,22 @@ subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) enddo end subroutine triDiagTS - +!> Calculates u_h and v_h (velocities at thickness points), +!! optionally using the entrainments (in m) passed in as arguments. subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, ea, eb) 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(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< The zonal velocity, in m s-1 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< The meridional velocity, in m s-1 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(out) :: u_h - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: v_h - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: ea - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: eb -! This subroutine calculates u_h and v_h (velocities at thickness -! points), optionally using the entrainments (in m) passed in as arguments. - -! Arguments: u - Zonal velocity, in m s-1. -! (in) v - Meridional velocity, in m s-1. -! (in) h - Layer thickness, in m or kg m-2. -! (out) u_h - The zonal velocity at thickness points after -! entrainment, in m s-1. -! (out) v_h - The meridional velocity at thickness points after -! entrainment, in m s-1. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in, opt) ea - The amount of fluid entrained from the layer above within -! this time step, in units of m or kg m-2. Omitting ea is the -! same as setting it to 0. -! (in, opt) eb - The amount of fluid entrained from the layer below within -! this time step, in units of m or kg m-2. Omitting eb is the -! same as setting it to 0. ea and eb must either be both -! present or both absent. - + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: u_h, v_h !< zonal and meridional velocity at thickness + !! points entrainment, in m s-1. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in), optional :: ea, eb !< The amount of fluid entrained + !! from the layer above within this time step + !! , in units of m or kg m-2. Omitting ea is the + !! same as setting it to 0. + + ! local variables real :: b_denom_1 ! The first term in the denominator of b1 in m or kg m-2. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected, in m or kg m-2. @@ -1318,26 +1239,20 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, & end subroutine applyBoundaryFluxesInOut +!> Initializes this module. subroutine diabatic_aux_init(Time, G, GV, param_file, diag, CS, useALEalgorithm, use_ePBL) type(time_type), intent(in) :: Time type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - type(diag_ctrl), target, intent(inout) :: diag - type(diabatic_aux_CS), pointer :: CS - logical, intent(in) :: useALEalgorithm - logical, intent(in) :: use_ePBL - -! Arguments: -! (in) Time = current model time -! (in) G = ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) param_file = structure indicating the open file to parse for parameter values -! (in) diag = structure used to regulate diagnostic output -! (in/out) CS = pointer set to point to the control structure for this module -! (in) use_ePBL = If true, use the implicit energetics planetary boundary -! layer scheme to determine the diffusivity in the -! surface boundary layer. + type(diag_ctrl), target, intent(inout) :: diag !< structure used to regulate diagnostic output + type(diabatic_aux_CS), pointer :: CS !< pointer set to point to the ontrol structure for + !! this module + logical, intent(in) :: useALEalgorithm !< If True, uses ALE. + logical, intent(in) :: use_ePBL !< If true, use the implicit energetics + !! planetary boundary layer scheme to determine the + !! diffusivity in the surface boundary layer. + ! local variables type(vardesc) :: vd ! This "include" declares and sets the variable "version". @@ -1460,4 +1375,48 @@ subroutine diabatic_aux_end(CS) end subroutine diabatic_aux_end +!> \namespace MOM_diabatic_aux +!! +!! This module contains the subroutines that, along with the * +!! subroutines that it calls, implements diapycnal mass and momentum * +!! fluxes and a bulk mixed layer. The diapycnal diffusion can be * +!! used without the bulk mixed layer. * +!! * +!! diabatic first determines the (diffusive) diapycnal mass fluxes * +!! based on the convergence of the buoyancy fluxes within each layer. * +!! The dual-stream entrainment scheme of MacDougall and Dewar (JPO, * +!! 1997) is used for combined diapycnal advection and diffusion, * +!! calculated implicitly and potentially with the Richardson number * +!! dependent mixing, as described by Hallberg (MWR, 2000). Diapycnal * +!! advection is fundamentally the residual of diapycnal diffusion, * +!! so the fully implicit upwind differencing scheme that is used is * +!! entirely appropriate. The downward buoyancy flux in each layer * +!! is determined from an implicit calculation based on the previously * +!! calculated flux of the layer above and an estimated flux in the * +!! layer below. This flux is subject to the following conditions: * +!! (1) the flux in the top and bottom layers are set by the boundary * +!! conditions, and (2) no layer may be driven below an Angstrom thick-* +!! ness. If there is a bulk mixed layer, the buffer layer is treat- * +!! ed as a fixed density layer with vanishingly small diffusivity. * +!! * +!! diabatic takes 5 arguments: the two velocities (u and v), the * +!! thicknesses (h), a structure containing the forcing fields, and * +!! the length of time over which to act (dt). The velocities and * +!! thickness are taken as inputs and modified within the subroutine. * +!! There is no limit on the time step. * +!! * +!! 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, T, S, buoy, ustar, 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). * +!! * +!!********+*********+*********+*********+*********+*********+*********+** + end module MOM_diabatic_aux diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 6316fd40e6..698243a7f6 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -10,6 +10,7 @@ module MOM_diabatic_driver 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_ddiff, only : CVMix_ddiff_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 @@ -97,6 +98,7 @@ 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_ddiff !< If true, use the CVMix double diffusion module. 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. @@ -249,7 +251,7 @@ module MOM_diabatic_driver integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap -integer :: id_clock_kpp +integer :: id_clock_kpp, id_clock_CVMix_ddiff contains @@ -385,7 +387,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect*h_neglect Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 - if (nz == 1) return showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("diabatic(), MOM_diabatic_driver.F90") @@ -487,13 +488,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif if (CS%ML_mix_first > 0.0) then -! This subroutine -! (1) Cools the mixed layer. -! (2) Performs convective adjustment by mixed layer entrainment. -! (3) Heats the mixed layer and causes it to detrain to -! Monin-Obukhov depth or minimum mixed layer depth. -! (4) Uses any remaining TKE to drive mixed layer entrainment. -! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) + ! This subroutine: + ! (1) Cools the mixed layer. + ! (2) Performs convective adjustment by mixed layer entrainment. + ! (3) Heats the mixed layer and causes it to detrain to + ! Monin-Obukhov depth or minimum mixed layer depth. + ! (4) Uses any remaining TKE to drive mixed layer entrainment. + ! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) call find_uv_at_h(u, v, h, u_h, v_h, G, GV) call cpu_clock_begin(id_clock_mixedlayer) @@ -528,11 +529,12 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)") if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G) endif - endif + endif ! end CS%bulkmixedlayer 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%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) @@ -589,7 +591,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & CS%int_tide_input%tideamp, CS%int_tide_input%Nb, dt, G, GV, CS%int_tide_CSp) endif if (showCallTree) call callTree_waypoint("done with propagate_int_tide (diabatic)") - endif + endif ! end CS%use_int_tides call cpu_clock_begin(id_clock_set_diffusivity) ! Sets: Kd, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S @@ -730,10 +732,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! a diffusivity and happen before KPP. But generally in MOM, we do not match ! KPP boundary layer to interior, so this diffusivity can be computed when convenient. if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T)) then - call cpu_clock_begin(id_clock_differential_diff) + call cpu_clock_begin(id_clock_CVMix_ddiff) call differential_diffuse_T_S(h, tv, visc, dt, G, GV) - call cpu_clock_end(id_clock_differential_diff) + call cpu_clock_end(id_clock_CVMix_ddiff) if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)") if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G) @@ -746,7 +748,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo endif - endif @@ -1381,6 +1382,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! 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) + if (associated(visc%Kv_slow)) & + call pass_var(visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) + call cpu_clock_end(id_clock_pass) if (.not. CS%useALEalgorithm) then @@ -1885,7 +1889,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, real :: Kd integer :: num_mode - logical :: use_temperature, differentialDiffusion + logical :: use_temperature type(vardesc) :: vd ! This "include" declares and sets the variable "version". @@ -1936,11 +1940,10 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, the diffusivity from ePBL is added to all\n"//& "other diffusivities. Otherwise, the larger of kappa-\n"//& "shear and ePBL diffusivities are used.", default=.true.) - call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialDiffusion, & - "If true, apply parameterization of double-diffusion.", & - default=.false. ) + CS%use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) CS%use_kappa_shear = kappa_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"//& @@ -2399,8 +2402,8 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=CLOCK_MODULE) id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=CLOCK_ROUTINE) id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=CLOCK_ROUTINE) - id_clock_differential_diff = -1 ; if (differentialDiffusion) & - id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=CLOCK_ROUTINE) + id_clock_CVMix_ddiff = -1 ; if (CS%use_CVMix_ddiff) & + id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion)', grain=CLOCK_ROUTINE) ! initialize the auxiliary diabatic driver module call diabatic_aux_init(Time, G, GV, param_file, diag, CS%diabatic_aux_CSp, & diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 9906083597..903868795a 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -23,6 +23,8 @@ module MOM_set_diffusivity 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 : CVMix_shear_end +use MOM_CVMix_ddiff, only : CVMix_ddiff_init, CVMix_ddiff_end, CVMix_ddiff_cs +use MOM_CVMix_ddiff, only : compute_ddiff_coeffs 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 @@ -43,104 +45,101 @@ module MOM_set_diffusivity public set_diffusivity_end type, public :: set_diffusivity_CS ; private - logical :: debug ! If true, write verbose checksums for debugging. - - logical :: bulkmixedlayer ! If true, a refined bulk mixed layer is used with - ! GV%nk_rho_varies variable density mixed & buffer - ! layers. - real :: FluxRi_max ! The flux Richardson number where the stratification is - ! 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 :: 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 - ! from the BBL mixing and the other diffusivities. - ! Otherwise, diffusivities from the BBL_mixing is - ! added. - logical :: use_LOTW_BBL_diffusivity ! If true, use simpler/less precise, BBL diffusivity. - logical :: LOTW_BBL_use_omega ! If true, use simpler/less precise, BBL diffusivity. - real :: BBL_effic ! efficiency with which the energy extracted - ! by bottom drag drives BBL diffusion (nondim) - real :: cdrag ! quadratic drag coefficient (nondim) - real :: IMax_decay ! inverse of a maximum decay scale for - ! bottom-drag driven turbulence, (1/m) - - real :: Kd ! interior diapycnal diffusivity (m2/s) - real :: Kd_min ! minimum diapycnal diffusivity (m2/s) - real :: Kd_max ! maximum increment for diapycnal diffusivity (m2/s) - ! Set to a negative value to have no limit. - real :: Kd_add ! uniform diffusivity added everywhere without - ! filtering or scaling (m2/s) - real :: Kv ! interior vertical viscosity (m2/s) - real :: Kdml ! mixed layer diapycnal diffusivity (m2/s) - ! when bulkmixedlayer==.false. - real :: Hmix ! mixed layer thickness (meter) when - ! bulkmixedlayer==.false. + logical :: debug !< If true, write verbose checksums for debugging. + + logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with + !! GV%nk_rho_varies variable density mixed & buffer + !! layers. + real :: FluxRi_max !< The flux Richardson number where the stratification is + !! 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 :: 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 + !! from the BBL mixing and the other diffusivities. + !! Otherwise, diffusivities from the BBL_mixing is + !! added. + logical :: use_LOTW_BBL_diffusivity !< If true, use simpler/less precise, BBL diffusivity. + logical :: LOTW_BBL_use_omega !< If true, use simpler/less precise, BBL diffusivity. + real :: BBL_effic !< efficiency with which the energy extracted + !! by bottom drag drives BBL diffusion (nondim) + real :: cdrag !< quadratic drag coefficient (nondim) + real :: IMax_decay !< inverse of a maximum decay scale for + !! bottom-drag driven turbulence, (1/m) + real :: Kv !< The interior vertical viscosity (m2/s) + real :: Kd !< interior diapycnal diffusivity (m2/s) + real :: Kd_min !< minimum diapycnal diffusivity (m2/s) + real :: Kd_max !< maximum increment for diapycnal diffusivity (m2/s) + !! Set to a negative value to have no limit. + real :: Kd_add !< uniform diffusivity added everywhere without + !! filtering or scaling (m2/s) + real :: Kdml !< mixed layer diapycnal diffusivity (m2/s) + !! when bulkmixedlayer==.false. + real :: Hmix !< mixed layer thickness (meter) when + !! bulkmixedlayer==.false. type(diag_ctrl), pointer :: diag ! structure to regulate diagn output timing - logical :: limit_dissipation ! If enabled, dissipation is limited to be larger - ! than the following: - real :: dissip_min ! Minimum dissipation (W/m3) - real :: dissip_N0 ! Coefficient a in minimum dissipation = a+b*N (W/m3) - real :: dissip_N1 ! Coefficient b in minimum dissipation = a+b*N (J/m3) - real :: dissip_N2 ! Coefficient c in minimum dissipation = c*N2 (W m-3 s2) - real :: dissip_Kd_min ! Minimum Kd (m2/s) with dissipatio Rho0*Kd_min*N^2 - - 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) - 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 - ! (1) The depth of the mixed layer, or - ! (2) An Ekman length scale. - ! Energy availble to drive mixing below the mixed layer is - ! given by E = ML_RAD_COEFF*MSTAR*USTAR**3. Optionally, if - ! ML_rad_TKE_decay is true, this is further reduced by a factor - ! of exp(-h_ML*Idecay_len_TkE), where Idecay_len_TKE is - ! calculated the same way as in the mixed layer code. - ! The diapycnal diffusivity is KD(k) = E/(N2(k)+OMEGA2), - ! where N2 is the squared buoyancy frequency (s-2) and OMEGA2 - ! is the rotation rate of the earth squared. - real :: ML_rad_kd_max ! Maximum diapycnal diffusivity due to turbulence - ! radiated from the base of the mixed layer (m2/s) - real :: ML_rad_efold_coeff ! non-dim coefficient to scale penetration depth - real :: ML_rad_coeff ! coefficient, which scales MSTAR*USTAR^3 to - ! obtain energy available for mixing below - ! mixed layer base (nondimensional) - logical :: ML_rad_TKE_decay ! If true, apply same exponential decay - ! to ML_rad as applied to the other surface - ! sources of TKE in the mixed layer code. - real :: ustar_min ! A minimum value of ustar to avoid numerical - ! problems (m/s). If the value is small enough, - ! this parameter should not affect the solution. - real :: TKE_decay ! ratio of natural Ekman depth to TKE decay scale (nondim) - real :: mstar ! ratio of friction velocity cubed to - ! TKE input to the mixed layer (nondim) - logical :: ML_use_omega ! If true, use absolute rotation rate instead - ! of the vertical component of rotation when - ! setting the decay scale for mixed layer turbulence. - real :: ML_omega_frac ! When setting the decay scale for turbulence, use - ! this fraction of the absolute rotation rate blended - ! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2. - 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 :: 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 :: 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) + logical :: limit_dissipation !< If enabled, dissipation is limited to be larger + !! than the following: + real :: dissip_min !< Minimum dissipation (W/m3) + real :: dissip_N0 !< Coefficient a in minimum dissipation = a+b*N (W/m3) + real :: dissip_N1 !< Coefficient b in minimum dissipation = a+b*N (J/m3) + real :: dissip_N2 !< Coefficient c in minimum dissipation = c*N2 (W m-3 s2) + real :: dissip_Kd_min !< Minimum Kd (m2/s) with dissipatio Rho0*Kd_min*N^2 + + 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) + 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 + !! (1) The depth of the mixed layer, or + !! (2) An Ekman length scale. + !! Energy availble to drive mixing below the mixed layer is + !! given by E = ML_RAD_COEFF*MSTAR*USTAR**3. Optionally, if + !! ML_rad_TKE_decay is true, this is further reduced by a factor + !! of exp(-h_ML*Idecay_len_TkE), where Idecay_len_TKE is + !! calculated the same way as in the mixed layer code. + !! The diapycnal diffusivity is KD(k) = E/(N2(k)+OMEGA2), + !! where N2 is the squared buoyancy frequency (s-2) and OMEGA2 + !! is the rotation rate of the earth squared. + real :: ML_rad_kd_max !< Maximum diapycnal diffusivity due to turbulence + !! radiated from the base of the mixed layer (m2/s) + real :: ML_rad_efold_coeff !< non-dim coefficient to scale penetration depth + real :: ML_rad_coeff !< coefficient, which scales MSTAR*USTAR^3 to + !! obtain energy available for mixing below + !! mixed layer base (nondimensional) + logical :: ML_rad_TKE_decay !< If true, apply same exponential decay + !! to ML_rad as applied to the other surface + !! sources of TKE in the mixed layer code. + real :: ustar_min !< A minimum value of ustar to avoid numerical + !! problems (m/s). If the value is small enough, + !! this parameter should not affect the solution. + real :: TKE_decay !< ratio of natural Ekman depth to TKE decay scale (nondim) + real :: mstar !! ratio of friction velocity cubed to + !! TKE input to the mixed layer (nondim) + logical :: ML_use_omega !< If true, use absolute rotation rate instead + !! of the vertical component of rotation when + !! setting the decay scale for mixed layer turbulence. + real :: ML_omega_frac !< When setting the decay scale for turbulence, use + !! this fraction of the absolute rotation rate blended + !! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2. + 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 :: use_CVMix_shear !< If true, use one of the CVMix modules to find + !! shear-driven diapycnal diffusivity. + logical :: use_CVMix_ddiff !< If true, enable double-diffusive mixing via CVMix. + logical :: simple_TKE_to_Kd !< If true, uses a simple estimate of Kd/TKE that + !! does not rely on a layer-formulation. 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_ddiff_cs), pointer :: CVMix_ddiff_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() @@ -158,11 +157,6 @@ module MOM_set_diffusivity integer :: id_N2 = -1 integer :: id_N2_z = -1 - integer :: id_KT_extra = -1 - integer :: id_KS_extra = -1 - integer :: id_KT_extra_z = -1 - integer :: id_KS_extra_z = -1 - end type set_diffusivity_CS type diffusivity_diags @@ -172,12 +166,9 @@ module MOM_set_diffusivity Kd_BBL => NULL(),& ! BBL diffusivity at interfaces (m2/s) Kd_work => NULL(),& ! layer integrated work by diapycnal mixing (W/m2) maxTKE => NULL(),& ! energy required to entrain to h_max (m3/s3) - TKE_to_Kd => NULL(),& ! conversion rate (~1.0 / (G_Earth + dRho_lay)) + TKE_to_Kd => NULL() ! conversion rate (~1.0 / (G_Earth + dRho_lay)) ! between TKE dissipated within a layer and Kd ! in that layer, in m2 s-1 / m3 s-3 = s2 m-1 - KT_extra => NULL(),& ! double diffusion diffusivity for temp (m2/s) - KS_extra => NULL() ! double diffusion diffusivity for saln (m2/s) - end type diffusivity_diags ! Clocks @@ -185,6 +176,17 @@ module MOM_set_diffusivity contains +!> Sets the interior vertical diffusion of scalars due to the following processes: +!! 1) Shear-driven mixing: two options, Jackson et at. and KPP interior; +!! 2) Background mixing via CVMix (Bryan-Lewis profile) or the scheme described by +!! Harrison & Hallberg, JPO 2008; +!! 3) Double-diffusion aplpied via CVMix; +!! 4) Tidal mixing: many options available, see MOM_tidal_mixing.F90; +!! In addition, this subroutine has the option to set the interior vertical +!! viscosity associated with processes 2-4 listed above, which is stored in +!! visc%Kv_slow. Vertical viscosity due to shear-driven mixing is passed via +!! visc%Kv_shear +!! GMM, TODO: add contribution from tidal mixing into visc%Kv_slow subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & G, GV, CS, Kd, Kd_int) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -196,9 +198,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & 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(in) :: u_h + intent(in) :: u_h !< zonal thickness transport m^2/s. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: v_h + intent(in) :: v_h !< meridional thickness transport m^2/s. type(thermo_var_ptrs), intent(inout) :: tv !< Structure with pointers to thermodynamic !! fields. Out is for tv%TempxPmE. type(forcing), intent(in) :: fluxes !< Structure of surface fluxes that may be @@ -226,17 +228,15 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! massless layers filled vertically by diffusion. real, dimension(SZI_(G),SZK_(G)) :: & - N2_lay, & ! squared buoyancy frequency associated with layers (1/s2) - maxTKE, & ! energy required to entrain to h_max (m3/s3) - TKE_to_Kd ! conversion rate (~1.0 / (G_Earth + dRho_lay)) between - ! TKE dissipated within a layer and Kd in that layer, in - ! m2 s-1 / m3 s-3 = s2 m-1. + N2_lay, & !< squared buoyancy frequency associated with layers (1/s2) + maxTKE, & !< energy required to entrain to h_max (m3/s3) + TKE_to_Kd !< conversion rate (~1.0 / (G_Earth + dRho_lay)) between + !< TKE dissipated within a layer and Kd in that layer, in + !< m2 s-1 / m3 s-3 = s2 m-1. real, dimension(SZI_(G),SZK_(G)+1) :: & - N2_int, & ! squared buoyancy frequency associated at interfaces (1/s2) - dRho_int, & ! locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? - KT_extra, & ! double difusion diffusivity on temperature (m2/sec) - KS_extra ! double difusion diffusivity on salinity (m2/sec) + N2_int, & !< squared buoyancy frequency associated at interfaces (1/s2) + dRho_int !< locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? real :: I_Rho0 ! inverse of Boussinesq density (m3/kg) real :: dissip ! local variable for dissipation calculations (W/m3) @@ -271,10 +271,16 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & use_EOS = associated(tv%eqn_of_state) - if ((CS%double_diffusion) .and. & + if ((CS%use_CVMix_ddiff) .and. & .not.(associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S)) ) & call MOM_error(FATAL, "set_diffusivity: visc%Kd_extra_T and "//& - "visc%Kd_extra_S must be associated when DOUBLE_DIFFUSION is true.") + "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF is true.") + + ! Set Kd, Kd_int and Kv_slow to constant values. + ! If nothing else is specified, this will be the value used. + Kd(:,:,:) = CS%Kd + Kd_int(:,:,:) = CS%Kd + if (associated(visc%Kv_slow)) visc%Kv_slow(:,:,:) = CS%Kv ! Set up arrays for diagnostics. @@ -293,12 +299,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%id_TKE_to_Kd > 0) then allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz)) ; dd%TKE_to_Kd(:,:,:) = 0.0 endif - if ((CS%id_KT_extra > 0) .or. (CS%id_KT_extra_z > 0)) then - allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1)) ; dd%KT_extra(:,:,:) = 0.0 - endif - if ((CS%id_KS_extra > 0) .or. (CS%id_KS_extra_z > 0)) then - allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1)) ; dd%KS_extra(:,:,:) = 0.0 - endif if ((CS%id_Kd_BBL > 0) .or. (CS%id_Kd_BBL_z > 0)) then allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0 endif @@ -341,6 +341,10 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & 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_shear, visc%Kv_shear,G,GV,CS%CVMix_shear_csp) + if (CS%debug) then + call hchksum(visc%Kd_shear, "after CVMix_shear visc%Kd_shear",G%HI) + call hchksum(visc%Kv_shear, "after CVMix_shear visc%Kv_shear",G%HI) + endif elseif (associated(visc%Kv_shear)) then visc%Kv_shear(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled endif @@ -352,8 +356,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! set surface diffusivities (CS%bkgnd_mixing_csp%Kd_sfc) call sfc_bkgnd_mixing(G, CS%bkgnd_mixing_csp) -! GMM, fix OMP calls below - !$OMP parallel do default(none) shared(is,ie,js,je,nz,G,GV,CS,h,tv,T_f,S_f,fluxes,dd, & !$OMP Kd,visc, & !$OMP Kd_int,dt,u,v,Omega2) & @@ -370,35 +372,13 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo endif - ! add background mixing + ! 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 - if (KS_extra(i,K) > KT_extra(i,K)) then ! salt fingering - Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KT_extra(i,K) - Kd(i,j,k) = Kd(i,j,k) + 0.5*KT_extra(i,K) - visc%Kd_extra_S(i,j,k) = KS_extra(i,K)-KT_extra(i,K) - visc%Kd_extra_T(i,j,k) = 0.0 - elseif (KT_extra(i,K) > 0.0) then ! double-diffusive convection - Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KS_extra(i,K) - Kd(i,j,k) = Kd(i,j,k) + 0.5*KS_extra(i,K) - visc%Kd_extra_T(i,j,k) = KT_extra(i,K)-KS_extra(i,K) - visc%Kd_extra_S(i,j,k) = 0.0 - else ! There is no double diffusion at this interface. - visc%Kd_extra_T(i,j,k) = 0.0 - visc%Kd_extra_S(i,j,k) = 0.0 - endif - enddo ; enddo - if (associated(dd%KT_extra)) then ; do K=1,nz+1 ; do i=is,ie - dd%KT_extra(i,j,K) = KT_extra(i,K) - enddo ; enddo ; endif - - if (associated(dd%KS_extra)) then ; do K=1,nz+1 ; do i=is,ie - dd%KS_extra(i,j,K) = KS_extra(i,K) - enddo ; enddo ; endif + ! Apply double diffusion + ! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available. + if (CS%use_CVMix_ddiff) then + call compute_ddiff_coeffs(h, tv, G, GV, j, visc%Kd_extra_T, visc%Kd_extra_S, CS%CVMix_ddiff_csp) endif ! Add the input turbulent diffusivity. @@ -496,6 +476,11 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%useKappaShear) call hchksum(visc%Kd_shear,"Turbulent Kd",G%HI,haloshift=0) + if (CS%use_CVMix_ddiff) then + call hchksum(visc%Kd_extra_T, "MOM_set_diffusivity: Kd_extra_T",G%HI,haloshift=0) + call hchksum(visc%Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S",G%HI,haloshift=0) + endif + 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.) @@ -512,12 +497,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & 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) @@ -538,13 +517,28 @@ 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) + ! post diagnostics - num_z_diags = 0 + ! background mixing + 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) + ! double diffusive mixing + if (CS%CVMix_ddiff_csp%id_KT_extra > 0) & + call post_data(CS%CVMix_ddiff_csp%id_KT_extra, visc%Kd_extra_T, CS%CVMix_ddiff_csp%diag) + if (CS%CVMix_ddiff_csp%id_KS_extra > 0) & + call post_data(CS%CVMix_ddiff_csp%id_KS_extra, visc%Kd_extra_S, CS%CVMix_ddiff_csp%diag) + if (CS%CVMix_ddiff_csp%id_R_rho > 0) & + call post_data(CS%CVMix_ddiff_csp%id_R_rho, CS%CVMix_ddiff_csp%R_rho, CS%CVMix_ddiff_csp%diag) + + if (CS%id_Kd_layer > 0) call post_data(CS%id_Kd_layer, Kd, CS%diag) + + ! tidal mixing call post_tidal_diagnostics(G,GV,h,CS%tm_csp) + num_z_diags = 0 if (CS%tm_csp%Int_tide_dissipation .or. CS%tm_csp%Lee_wave_dissipation .or. & CS%tm_csp%Lowmode_itidal_dissipation) then @@ -568,26 +562,11 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif - if (CS%id_KT_extra > 0) call post_data(CS%id_KT_extra, dd%KT_extra, CS%diag) - if (CS%id_KS_extra > 0) call post_data(CS%id_KS_extra, dd%KS_extra, CS%diag) if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, dd%Kd_BBL, CS%diag) - if (CS%id_KT_extra_z > 0) then - num_z_diags = num_z_diags + 1 - z_ids(num_z_diags) = CS%id_KT_extra_z - z_ptrs(num_z_diags)%p => dd%KT_extra - endif - - if (CS%id_KS_extra_z > 0) then - num_z_diags = num_z_diags + 1 - z_ids(num_z_diags) = CS%id_KS_extra_z - z_ptrs(num_z_diags)%p => dd%KS_extra - endif - if (CS%id_Kd_BBL_z > 0) then num_z_diags = num_z_diags + 1 z_ids(num_z_diags) = CS%id_Kd_BBL_z - z_ptrs(num_z_diags)%p => dd%KS_extra endif if (num_z_diags > 0) & @@ -598,8 +577,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (associated(dd%Kd_user)) deallocate(dd%Kd_user) 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) - if (associated(dd%KS_extra)) deallocate(dd%KS_extra) if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL) if (showCallTree) call callTree_leave("set_diffusivity()") @@ -952,119 +929,6 @@ 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 -!! what was in Large et al. (1994). All the coefficients here should probably -!! be made run-time variables rather than hard-coded constants. -!! -!! \todo Find reference for NCAR tech note above. -subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, Kd_T_dd, Kd_S_dd) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available - !! thermodynamic fields; absent fields have NULL - !! ptrs. - 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(in) :: T_f !< layer temp in C with the values in massless layers - !! filled vertically by diffusion. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: S_f !< Layer salinities in PPT with values in massless - !! layers filled vertically by diffusion. - integer, intent(in) :: j !< Meridional index upon which to work. - type(set_diffusivity_CS), pointer :: CS !< Module control structure. - real, dimension(SZI_(G),SZK_(G)+1), & - intent(out) :: Kd_T_dd !< Interface double diffusion diapycnal - !! diffusivity for temp (m2/sec). - real, dimension(SZI_(G),SZK_(G)+1), & - intent(out) :: Kd_S_dd !< Interface double diffusion diapycnal - !! diffusivity for saln (m2/sec). - -! Arguments: -! (in) tv - structure containing pointers to any available -! thermodynamic fields; absent fields have NULL ptrs -! (in) h - layer thickness (m or kg m-2) -! (in) T_f - layer temp in C with the values in massless layers -! filled vertically by diffusion -! (in) S_f - layer salinities in PPT with values in massless layers -! filled vertically by diffusion -! (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_T_dd - interface double diffusion diapycnal diffusivity for temp (m2/sec) -! (out) Kd_S_dd - interface double diffusion diapycnal diffusivity for saln (m2/sec) - -! 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 -! what was in Large et al. (1994). All the coefficients here should probably -! be made run-time variables rather than hard-coded constants. - - real, dimension(SZI_(G)) :: & - dRho_dT, & ! partial derivatives of density wrt temp (kg m-3 degC-1) - dRho_dS, & ! partial derivatives of density wrt saln (kg m-3 PPT-1) - pres, & ! pressure at each interface (Pa) - Temp_int, & ! temp and saln at interfaces - Salin_int - - real :: alpha_dT ! density difference between layers due to temp diffs (kg/m3) - real :: beta_dS ! density difference between layers due to saln diffs (kg/m3) - - real :: Rrho ! vertical density ratio - real :: diff_dd ! factor for double-diffusion - real :: prandtl ! flux ratio for diffusive convection regime - - real, parameter :: Rrho0 = 1.9 ! limit for double-diffusive density ratio - real, parameter :: dsfmax = 1.e-4 ! max diffusivity in case of salt fingering - real, parameter :: Kv_molecular = 1.5e-6 ! molecular viscosity (m2/sec) - - integer :: i, k, is, ie, nz - is = G%isc ; ie = G%iec ; nz = G%ke - - if (associated(tv%eqn_of_state)) then - do i=is,ie - pres(i) = 0.0 ; Kd_T_dd(i,1) = 0.0 ; Kd_S_dd(i,1) = 0.0 - Kd_T_dd(i,nz+1) = 0.0 ; Kd_S_dd(i,nz+1) = 0.0 - enddo - do K=2,nz - do i=is,ie - pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1) - Temp_Int(i) = 0.5 * (T_f(i,j,k-1) + T_f(i,j,k)) - Salin_Int(i) = 0.5 * (S_f(i,j,k-1) + S_f(i,j,k)) - enddo - call calculate_density_derivs(Temp_int, Salin_int, pres, & - dRho_dT(:), dRho_dS(:), is, ie-is+1, tv%eqn_of_state) - - do i=is,ie - alpha_dT = -1.0*dRho_dT(i) * (T_f(i,j,k-1) - T_f(i,j,k)) - beta_dS = dRho_dS(i) * (S_f(i,j,k-1) - S_f(i,j,k)) - - if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then ! salt finger case - Rrho = min(alpha_dT/beta_dS,Rrho0) - diff_dd = 1.0 - ((RRho-1.0)/(RRho0-1.0)) - diff_dd = dsfmax*diff_dd*diff_dd*diff_dd - Kd_T_dd(i,K) = 0.7*diff_dd - Kd_S_dd(i,K) = diff_dd - elseif ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then ! diffusive convection - Rrho = alpha_dT/beta_dS - diff_dd = Kv_molecular*0.909*exp(4.6*exp(-0.54*(1/Rrho-1))) - prandtl = 0.15*Rrho - if (Rrho > 0.5) prandtl = (1.85-0.85/Rrho)*Rrho - Kd_T_dd(i,K) = diff_dd - Kd_S_dd(i,K) = prandtl*diff_dd - else - Kd_T_dd(i,K) = 0.0 ; Kd_S_dd(i,K) = 0.0 - endif - enddo - enddo - endif - -end subroutine double_diffusion !> This routine adds diffusion sustained by flow energy extracted by bottom drag. subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & maxTKE, kb, G, GV, CS, Kd, Kd_int, Kd_BBL) @@ -1974,6 +1838,11 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp ! 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, "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.) + 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"//& @@ -2076,45 +1945,6 @@ 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.", & - default=.false.) - if (CS%double_diffusion) then - call get_param(param_file, mdl, "MAX_RRHO_SALT_FINGERS", CS%Max_Rrho_salt_fingers, & - "Maximum density ratio for salt fingering regime.", & - default=2.55, units="nondim") - call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", CS%Max_salt_diff_salt_fingers, & - "Maximum salt diffusivity for salt fingering regime.", & - default=1.e-4, units="m2 s-1") - call get_param(param_file, mdl, "KV_MOLECULAR", CS%Kv_molecular, & - "Molecular viscosity for calculation of fluxes under \n"//& - "double-diffusive convection.", default=1.5e-6, units="m2 s-1") - ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults. - - CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for temperature', 'm2 s-1') - - CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for salinity', 'm2 s-1') - - if (associated(diag_to_Z_CSp)) then - vd = var_desc("KT_extra", "m2 s-1", & - "Double-Diffusive Temperature Diffusivity, interpolated to z", & - z_grid='z') - CS%id_KT_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("KS_extra", "m2 s-1", & - "Double-Diffusive Salinity Diffusivity, interpolated to z",& - z_grid='z') - CS%id_KS_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("Kd_BBL", "m2 s-1", & - "Bottom Boundary Layer Diffusivity", z_grid='z') - CS%id_Kd_BBL_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - endif - endif - if (CS%user_change_diff) then call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) endif @@ -2130,6 +1960,9 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp ! CVMix shear-driven mixing CS%use_CVMix_shear = CVMix_shear_init(Time, G, GV, param_file, CS%diag, CS%CVMix_shear_csp) + ! CVMix double diffusion mixing + CS%use_CVMix_ddiff = CVMix_ddiff_init(Time, G, GV, param_file, CS%diag, CS%CVMix_ddiff_csp) + end subroutine set_diffusivity_init !> Clear pointers and dealocate memory @@ -2146,6 +1979,9 @@ subroutine set_diffusivity_end(CS) if (CS%use_CVMix_shear) & call CVMix_shear_end(CS%CVMix_shear_csp) + if (CS%use_CVMix_ddiff) & + call CVMix_ddiff_end(CS%CVMix_ddiff_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 ec0b5a80b3..ec1b09a5ad 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -2,38 +2,6 @@ module MOM_set_visc ! This file is part of MOM6. See LICENSE.md for the license. -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, April 1994 - October 2006 * -!* Quadratic Bottom Drag by James Stephens and R. Hallberg. * -!* * -!* This file contains the subroutine that calculates various values * -!* related to the bottom boundary layer, such as the viscosity and * -!* thickness of the BBL (set_viscous_BBL). This would also be the * -!* module in which other viscous quantities that are flow-independent * -!* might be set. This information is transmitted to other modules * -!* via a vertvisc type structure. * -!* * -!* The same code is used for the two velocity components, by * -!* indirectly referencing the velocities and defining a handful of * -!* direction-specific defined variables. * -!* * -!* Macros written all in capital letters are defined in MOM_memory.h. * -!* * -!* A small fragment of the grid is shown below: * -!* * -!* j+1 x ^ x ^ x At x: q * -!* j+1 > o > o > At ^: v, frhatv, tauy * -!* j x ^ x ^ x At >: u, frhatu, taux * -!* j > o > o > At o: h * -!* 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_debugging, only : uvchksum, hchksum use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr @@ -44,8 +12,9 @@ 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_conv, only : CVMix_conv_is_used +use MOM_cvmix_shear, only : cvmix_shear_is_used +use MOM_cvmix_conv, only : cvmix_conv_is_used +use MOM_CVMix_ddiff, only : CVMix_ddiff_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 @@ -1791,8 +1760,10 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., & do_not_log=.true.) + 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_shear = CVMix_shear_is_used(param_file) @@ -1811,7 +1782,9 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) 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_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 + + ! MOM_bkgnd_mixing is always used, so always allocate visc%Kv_slow. GMM + allocate(visc%Kv_slow(isd:ied,jsd:jed,nz+1)) ; visc%Kv_slow(:,:,:) = 0.0 vd = var_desc("Kd_shear","m2 s-1","Shear-driven turbulent diffusivity at interfaces", & hor_grid='h', z_grid='i') @@ -1854,21 +1827,14 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) type(set_visc_CS), pointer :: CS !< A pointer that is set to point to the control !! structure for this module type(ocean_OBC_type), pointer :: OBC -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (out) visc - A structure containing vertical viscosities and related -! fields. Allocated here. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module + + ! local variables real :: Csmag_chan_dflt, smag_const1, TKE_decay_dflt, bulk_Ri_ML_dflt real :: Kv_background real :: omega_frac_dflt integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, i, j, n - logical :: use_kappa_shear, adiabatic, differential_diffusion, use_omega + logical :: use_kappa_shear, adiabatic, use_omega + logical :: use_CVMix_ddiff type(OBC_segment_type), pointer :: segment ! pointer to OBC segment type ! This include declares and sets the variable "version". #include "version_variable.h" @@ -1891,8 +1857,8 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) ! Set default, read and log parameters call log_version(param_file, mdl, version, "") - CS%RiNo_mix = .false. - use_kappa_shear = .false. ; differential_diffusion = .false. !; adiabatic = .false. ! Needed? -AJA + CS%RiNo_mix = .false. ; use_CVMix_ddiff = .false. + use_kappa_shear = .false. !; adiabatic = .false. ! Needed? -AJA call get_param(param_file, mdl, "BOTTOMDRAGLAW", CS%bottomdraglaw, & "If true, the bottom stress is calculated with a drag \n"//& "law of the form c_drag*|u|*u. The velocity magnitude \n"//& @@ -1919,11 +1885,9 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) if (.not.adiabatic) then use_kappa_shear = kappa_shear_is_used(param_file) CS%RiNo_mix = use_kappa_shear - call get_param(param_file, mdl, "DOUBLE_DIFFUSION", differential_diffusion, & - "If true, increase diffusivitives for temperature or salt \n"//& - "based on double-diffusive paramaterization from MOM4/KPP.", & - default=.false.) + use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) endif + call get_param(param_file, mdl, "PRANDTL_TURB", visc%Prandtl_turb, & "The turbulent Prandtl number applied to shear \n"//& "instability.", units="nondim", default=1.0) @@ -2016,6 +1980,15 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) "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.) + + call get_param(param_file, mdl, "ADD_KV_SLOW", visc%add_Kv_slow, & + "If true, the background vertical viscosity in the interior \n"//& + "(i.e., tidal + background + shear + convenction) is addded \n"// & + "when computing the coupling coefficient. The purpose of this \n"// & + "flag is to be able to recover previous answers and it will likely \n"// & + "be removed in the future since this option should always be true.", & + default=.false.) + call get_param(param_file, mdl, "KV_BBL_MIN", CS%KV_BBL_min, & "The minimum viscosities in the bottom boundary layer.", & units="m2 s-1", default=Kv_background) @@ -2065,7 +2038,7 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) Time, 'Rayleigh drag velocity at v points', 'm s-1') endif - if (differential_diffusion) then + if (use_CVMix_ddiff) then allocate(visc%Kd_extra_T(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_T = 0.0 allocate(visc%Kd_extra_S(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_S = 0.0 endif @@ -2113,4 +2086,37 @@ subroutine set_visc_end(visc, CS) deallocate(CS) end subroutine set_visc_end +!> \namespace MOM_set_visc +!!********+*********+*********+*********+*********+*********+*********+** +!!* * +!!* By Robert Hallberg, April 1994 - October 2006 * +!!* Quadratic Bottom Drag by James Stephens and R. Hallberg. * +!!* * +!!* This file contains the subroutine that calculates various values * +!!* related to the bottom boundary layer, such as the viscosity and * +!!* thickness of the BBL (set_viscous_BBL). This would also be the * +!!* module in which other viscous quantities that are flow-independent * +!!* might be set. This information is transmitted to other modules * +!!* via a vertvisc type structure. * +!!* * +!!* The same code is used for the two velocity components, by * +!!* indirectly referencing the velocities and defining a handful of * +!!* direction-specific defined variables. * +!!* * +!!* Macros written all in capital letters are defined in MOM_memory.h. * +!!* * +!!* A small fragment of the grid is shown below: * +!!* * +!!* j+1 x ^ x ^ x At x: q * +!!* j+1 > o > o > At ^: v, frhatv, tauy * +!!* j x ^ x ^ x At >: u, frhatu, taux * +!!* j > o > o > At o: h * +!!* 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). * +!!* * +!!********+*********+*********+*********+*********+*********+*********+** + end module MOM_set_visc diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 48a6380ead..bafbe5eb59 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -2,7 +2,7 @@ module MOM_vert_friction ! This file is part of MOM6. See LICENSE.md for the license. - +use MOM_domains, only : pass_var, To_All, Omit_corners use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl use MOM_debugging, only : uvchksum, hchksum @@ -116,6 +116,7 @@ module MOM_vert_friction integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 integer :: id_Ray_u = -1, id_Ray_v = -1, id_taux_bot = -1, id_tauy_bot = -1 + integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 !>@} type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() @@ -614,6 +615,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) ! based on harmonic mean thicknesses, in m or kg m-2. h_ml ! The mixed layer depth, in m or kg m-2. real, allocatable, dimension(:,:) :: hML_u, hML_v + real, allocatable, dimension(:,:,:) :: Kv_v, & !< Total vertical viscosity at u-points + Kv_u !< Total vertical viscosity at v-points real :: zcol(SZI_(G)) ! The height of an interface at h-points, in m or kg m-2. real :: botfn ! A function which goes from 1 at the bottom to 0 much more ! than Hbbl into the interior. @@ -646,6 +649,14 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) I_Hbbl(:) = 1.0 / (CS%Hbbl * GV%m_to_H + h_neglect) I_valBL = 0.0 ; if (CS%harm_BL_val > 0.0) I_valBL = 1.0 / CS%harm_BL_val + if (CS%id_Kv_u > 0) then + allocate(Kv_u(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) ; Kv_u(:,:,:) = 0.0 + endif + + if (CS%id_Kv_v > 0) then + allocate(Kv_v(G%isd:G%ied,G%JsdB:G%JedB,G%ke)) ; Kv_v(:,:,:) = 0.0 + endif + if (CS%debug .or. (CS%id_hML_u > 0)) then allocate(hML_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; hML_u(:,:) = 0.0 endif @@ -821,6 +832,13 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) CS%h_u(I,j,k) = hvel(I,k) ; enddo ; enddo endif + ! Diagnose total Kv at u-points + if (CS%id_Kv_u > 0) then + do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) Kv_u(I,j,k) = 0.5 * (CS%a_u(I,j,K)+CS%a_u(I,j,K+1)) * CS%h_u(I,j,k) + enddo ; enddo + endif + enddo @@ -984,6 +1002,14 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) do K=1,nz+1 ; do i=is,ie ; if (do_i(i)) CS%a_v(i,J,K) = a(i,K) ; enddo ; enddo do k=1,nz ; do i=is,ie ; if (do_i(i)) CS%h_v(i,J,k) = hvel(i,k) ; enddo ; enddo endif + + ! Diagnose total Kv at v-points + if (CS%id_Kv_v > 0) then + do k=1,nz ; do i=is,ie + if (do_i(I)) Kv_v(i,J,k) = 0.5 * (CS%a_v(i,J,K)+CS%a_v(i,J,K+1)) * CS%h_v(i,J,k) + enddo ; enddo + endif + enddo ! end of v-point j loop if (CS%debug) then @@ -997,6 +1023,9 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) endif ! Offer diagnostic fields for averaging. + if (CS%id_Kv_slow > 0) call post_data(CS%id_Kv_slow, visc%Kv_slow, CS%diag) + if (CS%id_Kv_u > 0) call post_data(CS%id_Kv_u, Kv_u, CS%diag) + if (CS%id_Kv_v > 0) call post_data(CS%id_Kv_v, Kv_v, CS%diag) if (CS%id_au_vv > 0) call post_data(CS%id_au_vv, CS%a_u, CS%diag) if (CS%id_av_vv > 0) call post_data(CS%id_av_vv, CS%a_v, CS%diag) if (CS%id_h_u > 0) call post_data(CS%id_h_u, CS%h_u, CS%diag) @@ -1165,6 +1194,44 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m endif endif + ! add "slow" varying vertical viscosity (e.g., from background, tidal etc) + if (associated(visc%Kv_slow) .and. (visc%add_Kv_slow)) then + ! GMM/ A factor of 2 is also needed here, see comment above from BGR. + if (work_on_u) then + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + Kv_add(i,K) = Kv_add(i,K) + 1.0 * (visc%Kv_slow(i,j,k) + visc%Kv_slow(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) = Kv_add(i,K) + 2. * visc%Kv_slow(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) = Kv_add(i,K) + 2. * visc%Kv_slow(i+1,j,k) ; enddo + endif + endif ; enddo + endif + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + a(i,K) = a(i,K) + Kv_add(i,K) + endif ; enddo ; enddo + else + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + Kv_add(i,K) = Kv_add(i,K) + 1.0*(visc%Kv_slow(i,j,k) + visc%Kv_slow(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) = Kv_add(i,K) + 2. * visc%Kv_slow(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) = Kv_add(i,K) + 2. * visc%Kv_slow(i,j+1,k) ; enddo + endif + endif ; enddo + endif + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + a(i,K) = a(i,K) + Kv_add(i,K) + endif ; enddo ; enddo + endif + endif + do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then ! botfn determines when a point is within the influence of the bottom ! boundary layer, going from 1 at the bottom to 0 in the interior. @@ -1671,17 +1738,30 @@ subroutine vertvisc_init(MIS, Time, G, GV, param_file, diag, ADp, dirs, & ALLOC_(CS%a_v(isd:ied,JsdB:JedB,nz+1)) ; CS%a_v(:,:,:) = 0.0 ALLOC_(CS%h_v(isd:ied,JsdB:JedB,nz)) ; CS%h_v(:,:,:) = 0.0 + CS%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, Time, & + 'Slow varying vertical viscosity', 'm2 s-1') + + CS%id_Kv_u = register_diag_field('ocean_model', 'Kv_u', diag%axesCuL, Time, & + 'Total vertical viscosity at u-points', 'm2 s-1') + + CS%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, Time, & + 'Total vertical viscosity at v-points', 'm2 s-1') + CS%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, Time, & 'Zonal Viscous Vertical Coupling Coefficient', 'm s-1') + CS%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, Time, & 'Meridional Viscous Vertical Coupling Coefficient', 'm s-1') CS%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, Time, & 'Thickness at Zonal Velocity Points for Viscosity', thickness_units) + CS%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, Time, & 'Thickness at Meridional Velocity Points for Viscosity', thickness_units) + CS%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, Time, & 'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', thickness_units) + CS%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, Time, & 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', thickness_units)