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..963547f3c0 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. @@ -277,6 +279,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & type(diabatic_CS), pointer :: CS !< module control structure type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves + ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & ea, & ! amount of fluid entrained from the layer above within ! one time step (m for Bouss, kg/m^2 for non-Bouss) @@ -385,11 +388,12 @@ 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") + if (.not. (CS%useALEalgorithm)) call MOM_error(FATAL, "MOM_diabatic_driver: "// & + "The ALE algorithm must be enabled when using MOM_diabatic_driver.") ! Offer diagnostics of various state varables at the start of diabatic ! these are mostly for debugging purposes. @@ -451,7 +455,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%id_frazil_h > 0) call post_data(CS%id_frazil_h, h, CS%diag) endif call disable_averaging(CS%diag) - endif + endif !associated(tv%T) .AND. associated(tv%frazil) + ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep call enable_averaging(dt, Time_end, CS%diag) if (CS%debugConservation) call MOM_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, G) @@ -481,58 +486,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (associated(CS%optics)) & call set_opacity(CS%optics, fluxes, G, GV, CS%opacity_CSp) - if (CS%bulkmixedlayer) then - if (CS%debug) then - call MOM_forcing_chksum("Before mixedlayer", fluxes, G, haloshift=0) - 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) - call find_uv_at_h(u, v, h, u_h, v_h, G, GV) - - call cpu_clock_begin(id_clock_mixedlayer) - if (CS%ML_mix_first < 1.0) then - ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???) - call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*CS%ML_mix_first, & - eaml,ebml, G, GV, CS%bulkmixedlayer_CSp, CS%optics, & - Hml, CS%aggregate_FW_forcing, dt, last_call=.false.) - if (CS%salt_reject_below_ML) & - call insert_brine(h, tv, G, GV, fluxes, nkmb, CS%diabatic_aux_CSp, & - dt*CS%ML_mix_first, CS%id_brine_lay) - else - ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???) - call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, & - G, GV, CS%bulkmixedlayer_CSp, CS%optics, & - Hml, CS%aggregate_FW_forcing, dt, last_call=.true.) - endif - - ! Keep salinity from falling below a small but positive threshold. - ! This constraint is needed for SIS1 ice model, which can extract - ! more salt than is present in the ocean. SIS2 does not suffer - ! from this limitation, in which case we can let salinity=0 and still - ! have salt conserved with SIS2 ice. So for SIS2, we can run with - ! BOUND_SALINITY=False in MOM.F90. - if (associated(tv%S) .and. associated(tv%salt_deficit)) & - call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - call cpu_clock_end(id_clock_mixedlayer) - if (CS%debug) then - call MOM_state_chksum("After mixedlayer ", u, v, h, G, GV, haloshift=0) - call MOM_forcing_chksum("After mixedlayer", fluxes, G, haloshift=0) - endif - 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 - - if (CS%debug) then + if (CS%debug) & 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 +545,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 @@ -607,7 +563,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call hchksum(Kd_Int, "after set_diffusivity Kd_Int",G%HI,haloshift=0) endif - if (CS%useKPP) then call cpu_clock_begin(id_clock_kpp) ! KPP needs the surface buoyancy flux but does not update state variables. @@ -729,11 +684,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! If using matching within the KPP scheme, then this step needs to provide ! 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) + if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T) .and. .not. & + CS%use_CVMix_ddiff) then + call cpu_clock_begin(id_clock_differential_diff) call differential_diffuse_T_S(h, tv, visc, dt, G, GV) call cpu_clock_end(id_clock_differential_diff) + 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,52 +703,26 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo endif - endif ! This block sets ea, eb from Kd or Kd_int. - ! If using ALE algorithm, set ea=eb=Kd_int on interfaces for - ! use in the tri-diagonal solver. - ! Otherwise, call entrainment_diffusive() which sets ea and eb - ! based on KD and target densities (ie. does remapping as well). - if (CS%useALEalgorithm) then + ! set ea=eb=Kd_int on interfaces for use in the tri-diagonal solver. - do j=js,je ; do i=is,ie - ea(i,j,1) = 0. - enddo ; enddo + do j=js,je ; do i=is,ie + ea(i,j,1) = 0. + enddo ; enddo !$OMP parallel do default(none) shared(is,ie,js,je,nz,h_neglect,h,ea,GV,dt,Kd_int,eb) & !$OMP private(hval) - do k=2,nz ; do j=js,je ; do i=is,ie - hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) - ea(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_int(i,j,k) - eb(i,j,k-1) = ea(i,j,k) - enddo ; enddo ; enddo - do j=js,je ; do i=is,ie - eb(i,j,nz) = 0. - enddo ; enddo - if (showCallTree) call callTree_waypoint("done setting ea,eb from Kd_int (diabatic)") - - else ! .not. CS%useALEalgorithm - ! When not using ALE, calculate layer entrainments/detrainments from - ! diffusivities and differences between layer and target densities - call cpu_clock_begin(id_clock_entrain) - ! Calculate appropriately limited diapycnal mass fluxes to account - ! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb - call Entrainment_diffusive(u, v, h, tv, fluxes, dt, G, GV, CS%entrain_diffusive_CSp, & - ea, eb, kb, Kd_Lay=Kd, Kd_int=Kd_int) - call cpu_clock_end(id_clock_entrain) - if (showCallTree) call callTree_waypoint("done with Entrainment_diffusive (diabatic)") - - endif ! endif for (CS%useALEalgorithm) - - if (CS%debug) then - call MOM_forcing_chksum("after calc_entrain ", fluxes, G, haloshift=0) - call MOM_thermovar_chksum("after calc_entrain ", tv, G) - call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, haloshift=0) - call hchksum(ea, "after calc_entrain ea", G%HI, haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after calc_entrain eb", G%HI, haloshift=0, scale=GV%H_to_m) - endif + do k=2,nz ; do j=js,je ; do i=is,ie + hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) + ea(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_int(i,j,k) + eb(i,j,k-1) = ea(i,j,k) + enddo ; enddo ; enddo + do j=js,je ; do i=is,ie + eb(i,j,nz) = 0. + enddo ; enddo + if (showCallTree) call callTree_waypoint("done setting ea,eb from Kd_int (diabatic)") ! Save fields before boundary forcing is applied for tendency diagnostics if (CS%boundary_forcing_tendency_diag) then @@ -802,96 +733,90 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo endif - ! Apply forcing when using the ALE algorithm - if (CS%useALEalgorithm) then - call cpu_clock_begin(id_clock_remap) - - ! Changes made to following fields: h, tv%T and tv%S. - - do k=1,nz ; do j=js,je ; do i=is,ie - h_prebound(i,j,k) = h(i,j,k) - enddo ; enddo ; enddo - if (CS%use_energetic_PBL) then - - skinbuoyflux(:,:) = 0.0 - call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & - h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) - - if (CS%debug) then - call hchksum(ea, "after applyBoundaryFluxes ea",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after applyBoundaryFluxes eb",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(cTKE, "after applyBoundaryFluxes cTKE",G%HI,haloshift=0) - call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT",G%HI,haloshift=0) - call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS",G%HI,haloshift=0) - endif + ! Apply forcing + call cpu_clock_begin(id_clock_remap) - call find_uv_at_h(u, v, h, u_h, v_h, G, GV) - call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, & - CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) - - ! If visc%MLD exists, copy the ePBL's MLD into it - if (associated(visc%MLD)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G) - call pass_var(visc%MLD, G%domain, halo=1) - Hml(:,:) = visc%MLD(:,:) - endif + ! Changes made to following fields: h, tv%T and tv%S. + do k=1,nz ; do j=js,je ; do i=is,ie + h_prebound(i,j,k) = h(i,j,k) + enddo ; enddo ; enddo + if (CS%use_energetic_PBL) then - ! Augment the diffusivities due to those diagnosed in energetic_PBL. - do K=2,nz ; do j=js,je ; do i=is,ie + skinbuoyflux(:,:) = 0.0 + call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & + h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) - if (CS%ePBL_is_additive) then - Kd_add_here = 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_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) - eb(i,j,k-1) = eb(i,j,k-1) + Ent_int - ea(i,j,k) = ea(i,j,k) + Ent_int - Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here + if (CS%debug) then + call hchksum(ea, "after applyBoundaryFluxes ea",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "after applyBoundaryFluxes eb",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(cTKE, "after applyBoundaryFluxes cTKE",G%HI,haloshift=0) + call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT",G%HI,haloshift=0) + call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS",G%HI,haloshift=0) + endif - ! for diagnostics - Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_int(i,j,K) - Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_int(i,j,K) + call find_uv_at_h(u, v, h, u_h, v_h, G, GV) + call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, & + CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) - enddo ; enddo ; enddo + ! If visc%MLD exists, copy the ePBL's MLD into it + if (associated(visc%MLD)) then + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G) + call pass_var(visc%MLD, G%domain, halo=1) + Hml(:,:) = visc%MLD(:,:) + endif - if (CS%debug) then - call hchksum(ea, "after ePBL ea",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after ePBL eb",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(Kd_ePBL, "after ePBL Kd_ePBL",G%HI,haloshift=0) + ! Augment the diffusivities due to those diagnosed in energetic_PBL. + do K=2,nz ; do j=js,je ; do i=is,ie + if (CS%ePBL_is_additive) then + Kd_add_here = 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_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) + eb(i,j,k-1) = eb(i,j,k-1) + Ent_int + ea(i,j,k) = ea(i,j,k) + Ent_int + Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here - else - call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & - h, tv, CS%aggregate_FW_forcing, & - CS%evap_CFL_limit, CS%minimum_forcing_depth) - - endif ! endif for CS%use_energetic_PBL - - ! diagnose the tendencies due to boundary forcing - ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme - ! so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards - if (CS%boundary_forcing_tendency_diag) then - call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, G, GV, CS) - if (CS%id_boundary_forcing_h > 0) call post_data(CS%id_boundary_forcing_h, h, CS%diag, alt_h = h_diag) - endif - ! Boundary fluxes may have changed T, S, and h - call diag_update_remap_grids(CS%diag) + ! for diagnostics + Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_int(i,j,K) + Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_int(i,j,K) + + enddo ; enddo ; enddo - call cpu_clock_end(id_clock_remap) if (CS%debug) then - call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, haloshift=0) - call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G) - call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, haloshift=0) + call hchksum(ea, "after ePBL ea",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "after ePBL eb",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(Kd_ePBL, "after ePBL Kd_ePBL",G%HI,haloshift=0) endif - if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)") - if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, G) - endif ! endif for (CS%useALEalgorithm) + else + call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & + h, tv, CS%aggregate_FW_forcing, & + CS%evap_CFL_limit, CS%minimum_forcing_depth) + + endif ! endif for CS%use_energetic_PBL + + ! diagnose the tendencies due to boundary forcing + ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme + ! so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards + if (CS%boundary_forcing_tendency_diag) then + call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, G, GV, CS) + if (CS%id_boundary_forcing_h > 0) call post_data(CS%id_boundary_forcing_h, h, CS%diag, alt_h = h_diag) + endif + ! Boundary fluxes may have changed T, S, and h + call diag_update_remap_grids(CS%diag) + call cpu_clock_end(id_clock_remap) + if (CS%debug) then + call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, haloshift=0) + call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G) + call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, haloshift=0) + endif + if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)") + if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, G) ! Update h according to divergence of the difference between ! ea and eb. We keep a record of the original h in hold. @@ -901,6 +826,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! enough iterations are permitted in Calculate_Entrainment. ! Even if too few iterations are allowed, it is still guarded ! against. In other words the checks are probably unnecessary. + + ! GMM, should the code below be deleted? eb(i,j,k-1) = ea(i,j,k), + ! see above, so h should not change. + !$OMP parallel do default(shared) do j=js,je do i=is,ie @@ -936,204 +865,55 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G) - ! Here, T and S are updated according to ea and eb. - ! If using the bulk mixed layer, T and S are also updated - ! by surface fluxes (in fluxes%*). - ! This is a very long block. - if (CS%bulkmixedlayer) then + ! calculate change in temperature & salinity due to dia-coordinate surface diffusion + if (associated(tv%T)) then - if (associated(tv%T)) then - call cpu_clock_begin(id_clock_tridiag) - ! Temperature and salinity (as state variables) are treated - ! differently from other tracers to insure massless layers that - ! are lighter than the mixed layer have temperatures and salinities - ! that correspond to their prescribed densities. - if (CS%massless_match_targets) then - !$OMP parallel do default (shared) private(h_tr,b1,d1,c1,b_denom_1) - do j=js,je - do i=is,ie - h_tr = hold(i,j,1) + h_neglect - b1(i) = 1.0 / (h_tr + eb(i,j,1)) - d1(i) = h_tr * b1(i) - tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1)) - tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1)) - enddo - do k=2,nkmb ; do i=is,ie - c1(i,k) = eb(i,j,k-1) * b1(i) - h_tr = hold(i,j,k) + h_neglect - b_denom_1 = h_tr + d1(i)*ea(i,j,k) - b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) - if (k kb(i,j)) then - c1(i,k) = eb(i,j,k-1) * b1(i) - h_tr = hold(i,j,k) + h_neglect - b_denom_1 = h_tr + d1(i)*ea(i,j,k) - b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) - d1(i) = b_denom_1 * b1(i) - tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1)) - tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1)) - elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j)) - ! The bottommost buffer layer might entrain all the mass from some - ! of the interior layers that are thin and lighter in the coordinate - ! density than that buffer layer. The T and S of these newly - ! massless interior layers are unchanged. - tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k) - tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k) - endif - enddo ; enddo - - do k=nz-1,nkmb,-1 ; do i=is,ie - if (k >= kb(i,j)) then - tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) - tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) - endif - enddo ; enddo - do i=is,ie ; if (kb(i,j) <= nz) then - tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j)) - tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j)) - endif ; enddo - do k=nkmb-1,1,-1 ; do i=is,ie - tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) - tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) - enddo ; enddo - enddo ! end of j loop - else ! .not. massless_match_targets - ! This simpler form allows T & S to be too dense for the layers - ! between the buffer layers and the interior. - ! Changes: T, S - if (CS%tracer_tridiag) then - call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV) - call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV) - else - call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) - endif - endif ! massless_match_targets - call cpu_clock_end(id_clock_tridiag) + if (CS%debug) then + call hchksum(ea, "before triDiagTS ea ",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "before triDiagTS eb ",G%HI,haloshift=0, scale=GV%H_to_m) + endif + call cpu_clock_begin(id_clock_tridiag) - endif ! endif for associated(T) - if (CS%debugConservation) call MOM_state_stats('BML tridiag', u, v, h, tv%T, tv%S, G) + ! Keep salinity from falling below a small but positive threshold. + ! This constraint is needed for SIS1 ice model, which can extract + ! more salt than is present in the ocean. SIS2 does not suffer + ! from this limitation, in which case we can let salinity=0 and still + ! have salt conserved with SIS2 ice. So for SIS2, we can run with + ! BOUND_SALINITY=False in MOM.F90. + if (associated(tv%S) .and. associated(tv%salt_deficit)) & + call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then - ! The mixed layer code has already been called, but there is some needed - ! bookkeeping. - !$OMP parallel do default(shared) + if (CS%diabatic_diff_tendency_diag) then do k=1,nz ; do j=js,je ; do i=is,ie - hold(i,j,k) = h_orig(i,j,k) - ea(i,j,k) = ea(i,j,k) + eaml(i,j,k) - eb(i,j,k) = eb(i,j,k) + ebml(i,j,k) + temp_diag(i,j,k) = tv%T(i,j,k) + saln_diag(i,j,k) = tv%S(i,j,k) enddo ; enddo ; enddo - if (CS%debug) then - call hchksum(ea, "after ea = ea + eaml",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after eb = eb + ebml",G%HI,haloshift=0, scale=GV%H_to_m) - endif - endif - - if (CS%ML_mix_first < 1.0) then - ! Call the mixed layer code now, perhaps for a second time. - ! 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 the buffer layer into two isopycnal layers. - - call find_uv_at_h(u, v, hold, u_h, v_h, G, GV, ea, eb) - if (CS%debug) call MOM_state_chksum("find_uv_at_h1 ", u, v, h, G, GV, haloshift=0) - - dt_mix = min(dt,dt*(1.0 - CS%ML_mix_first)) - call cpu_clock_begin(id_clock_mixedlayer) - ! Changes: h, tv%T, tv%S, ea and eb (G is also inout???) - call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, & - G, GV, CS%bulkmixedlayer_CSp, CS%optics, & - Hml, CS%aggregate_FW_forcing, dt, last_call=.true.) - - if (CS%salt_reject_below_ML) & - call insert_brine(h, tv, G, GV, fluxes, nkmb, CS%diabatic_aux_CSp, dt_mix, & - CS%id_brine_lay) - - ! Keep salinity from falling below a small but positive threshold. - ! This constraint is needed for SIS1 ice model, which can extract - ! more salt than is present in the ocean. SIS2 does not suffer - ! from this limitation, in which case we can let salinity=0 and still - ! have salt conserved with SIS2 ice. So for SIS2, we can run with - ! BOUND_SALINITY=False in MOM.F90. - if (associated(tv%S) .and. associated(tv%salt_deficit)) & - call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - - call cpu_clock_end(id_clock_mixedlayer) - if (showCallTree) call callTree_waypoint("done with 2nd bulkmixedlayer (diabatic)") - if (CS%debugConservation) call MOM_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, G) endif - else ! following block for when NOT using BULKMIXEDLAYER - - - ! calculate change in temperature & salinity due to dia-coordinate surface diffusion - if (associated(tv%T)) then - - if (CS%debug) then - call hchksum(ea, "before triDiagTS ea ",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "before triDiagTS eb ",G%HI,haloshift=0, scale=GV%H_to_m) - endif - call cpu_clock_begin(id_clock_tridiag) - - ! Keep salinity from falling below a small but positive threshold. - ! This constraint is needed for SIS1 ice model, which can extract - ! more salt than is present in the ocean. SIS2 does not suffer - ! from this limitation, in which case we can let salinity=0 and still - ! have salt conserved with SIS2 ice. So for SIS2, we can run with - ! BOUND_SALINITY=False in MOM.F90. - if (associated(tv%S) .and. associated(tv%salt_deficit)) & - call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - - if (CS%diabatic_diff_tendency_diag) then - do k=1,nz ; do j=js,je ; do i=is,ie - temp_diag(i,j,k) = tv%T(i,j,k) - saln_diag(i,j,k) = tv%S(i,j,k) - enddo ; enddo ; enddo - endif - - ! Changes T and S via the tridiagonal solver; no change to h - if (CS%tracer_tridiag) then - call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV) - call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV) - else - call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) - endif - - ! diagnose temperature, salinity, heat, and salt tendencies - ! Note: hold here refers to the thicknesses from before the dual-entraintment when using - ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will have changed - ! In either case, tendencies should be posted on hold - if (CS%diabatic_diff_tendency_diag) then - call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, G, GV, CS) - if (CS%id_diabatic_diff_h > 0) call post_data(CS%id_diabatic_diff_h, hold, CS%diag, alt_h = hold) - endif + ! Changes T and S via the tridiagonal solver; no change to h + if (CS%tracer_tridiag) then + call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV) + call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV) + else - call cpu_clock_end(id_clock_tridiag) - if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)") + call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) + endif - endif ! endif corresponding to if (associated(tv%T)) - if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G) + ! diagnose temperature, salinity, heat, and salt tendencies + ! Note: hold here refers to the thicknesses from before the dual-entraintment when using + ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will have changed + ! In either case, tendencies should be posted on hold + if (CS%diabatic_diff_tendency_diag) then + call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, G, GV, CS) + if (CS%id_diabatic_diff_h > 0) call post_data(CS%id_diabatic_diff_h, hold, CS%diag, alt_h = hold) + endif + call cpu_clock_end(id_clock_tridiag) + if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)") - endif ! endif for the BULKMIXEDLAYER block + endif ! endif corresponding to if (associated(tv%T)) + if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G) if (CS%debug) then call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, haloshift=0) @@ -1142,14 +922,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call hchksum(eb, "after mixed layer eb", G%HI, scale=GV%H_to_m) endif - if (.not. CS%useALEalgorithm) then - call cpu_clock_begin(id_clock_remap) - call regularize_layers(h, tv, dt, ea, eb, G, GV, CS%regularize_layers_CSp) - call cpu_clock_end(id_clock_remap) - if (showCallTree) call callTree_waypoint("done with regularize_layers (diabatic)") - if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, tv%T, tv%S, G) - endif - ! Whenever thickness changes let the diag manager know, as the ! target grids for vertical remapping may need to be regenerated. call diag_update_remap_grids(CS%diag) @@ -1186,6 +958,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! mixing of passive tracers from massless boundary layers to interior call cpu_clock_begin(id_clock_tracers) + if (CS%mix_boundary_tracers) then Tr_ea_BBL = sqrt(dt*CS%Kd_BBL_tr) !$OMP parallel do default(shared) private(htot,in_boundary,add_ent) @@ -1234,17 +1007,12 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo - if (CS%useALEalgorithm) then ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied ! so hold should be h_orig - call call_tracer_column_fns(h_prebound, h, ea, eb, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug, & - evap_CFL_limit = CS%evap_CFL_limit, & - minimum_forcing_depth = CS%minimum_forcing_depth) - else - call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug) - endif + call call_tracer_column_fns(h_prebound, h, ea, eb, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug, & + evap_CFL_limit = CS%evap_CFL_limit, & + minimum_forcing_depth = CS%minimum_forcing_depth) elseif (associated(visc%Kd_extra_S)) then ! extra diffusivity for passive tracers @@ -1264,56 +1032,31 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & eatr(i,j,k) = ea(i,j,k) + add_ent enddo ; enddo ; enddo - if (CS%useALEalgorithm) then ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied - call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug,& - evap_CFL_limit = CS%evap_CFL_limit, & - minimum_forcing_depth = CS%minimum_forcing_depth) - else - call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug) - endif + call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug,& + evap_CFL_limit = CS%evap_CFL_limit, & + minimum_forcing_depth = CS%minimum_forcing_depth) else - if (CS%useALEalgorithm) then ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied - call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug, & - evap_CFL_limit = CS%evap_CFL_limit, & - minimum_forcing_depth = CS%minimum_forcing_depth) - else - call call_tracer_column_fns(hold, h, ea, eb, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug) - endif + call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug, & + evap_CFL_limit = CS%evap_CFL_limit, & + minimum_forcing_depth = CS%minimum_forcing_depth) endif ! (CS%mix_boundary_tracers) - - call cpu_clock_end(id_clock_tracers) - ! sponges if (CS%use_sponge) then call cpu_clock_begin(id_clock_sponge) if (associated(CS%ALE_sponge_CSp)) then ! ALE sponge call apply_ALE_sponge(h, dt, G, CS%ALE_sponge_CSp, CS%Time) - else - ! Layer mode sponge - if (CS%bulkmixedlayer .and. associated(tv%eqn_of_state)) then - do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo - !$OMP parallel do default(shared) - do j=js,je - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), & - is, ie-is+1, tv%eqn_of_state) - enddo - call apply_sponge(h, dt, G, GV, ea, eb, CS%sponge_CSp, Rcv_ml) - else - call apply_sponge(h, dt, G, GV, ea, eb, CS%sponge_CSp) - endif endif + call cpu_clock_end(id_clock_sponge) if (CS%debug) then call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, haloshift=0) @@ -1336,29 +1079,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo endif -! For momentum, it is only the net flux that homogenizes within -! the mixed layer. Vertical viscosity that is proportional to the -! mixed layer turbulence is applied elsewhere. - if (CS%bulkmixedlayer) then - if (CS%debug) then - call hchksum(ea, "before net flux rearrangement ea",G%HI, scale=GV%H_to_m) - call hchksum(eb, "before net flux rearrangement eb",G%HI, scale=GV%H_to_m) - endif - !$OMP parallel do default(shared) private(net_ent) - do j=js,je - do K=2,GV%nkml ; do i=is,ie - net_ent = ea(i,j,k) - eb(i,j,k-1) - ea(i,j,k) = max(net_ent, 0.0) - eb(i,j,k-1) = max(-net_ent, 0.0) - enddo ; enddo - enddo - if (CS%debug) then - call hchksum(ea, "after net flux rearrangement ea",G%HI, scale=GV%H_to_m) - call hchksum(eb, "after net flux rearrangement eb",G%HI, scale=GV%H_to_m) - endif - endif - -! Initialize halo regions of ea, eb, and hold to default values. + ! Initialize halo regions of ea, eb, and hold to default values. !$OMP parallel do default(shared) do k=1,nz do i=is-1,ie+1 @@ -1381,85 +1102,10 @@ 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) - call cpu_clock_end(id_clock_pass) + if (associated(visc%Kv_slow)) & + call pass_var(visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) - if (.not. CS%useALEalgorithm) then - ! Use a tridiagonal solver to determine effect of the diapycnal - ! advection on velocity field. It is assumed that water leaves - ! or enters the ocean with the surface velocity. - if (CS%debug) then - call MOM_state_chksum("before u/v tridiag ", u, v, h, G, GV, haloshift=0) - call hchksum(ea, "before u/v tridiag ea",G%HI, scale=GV%H_to_m) - call hchksum(eb, "before u/v tridiag eb",G%HI, scale=GV%H_to_m) - call hchksum(hold, "before u/v tridiag hold",G%HI, scale=GV%H_to_m) - endif - call cpu_clock_begin(id_clock_tridiag) - !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval) - do j=js,je - do I=Isq,Ieq - if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,1) = u(I,j,1) - hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect - b1(I) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1))) - d1(I) = hval * b1(I) - u(I,j,1) = b1(I) * (hval * u(I,j,1)) - enddo - do k=2,nz ; do I=Isq,Ieq - if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,k) = u(I,j,k) - c1(I,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(I) - eaval = ea(i,j,k) + ea(i+1,j,k) - hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect - b1(I) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(I)*eaval)) - d1(I) = (hval + d1(I)*eaval) * b1(I) - u(I,j,k) = (hval*u(I,j,k) + eaval*u(I,j,k-1))*b1(I) - enddo ; enddo - do k=nz-1,1,-1 ; do I=Isq,Ieq - u(I,j,k) = u(I,j,k) + c1(I,k+1)*u(I,j,k+1) - if (associated(ADp%du_dt_dia)) & - ADp%du_dt_dia(I,j,k) = (u(I,j,k) - ADp%du_dt_dia(I,j,k)) * Idt - enddo ; enddo - if (associated(ADp%du_dt_dia)) then - do I=Isq,Ieq - ADp%du_dt_dia(I,j,nz) = (u(I,j,nz)-ADp%du_dt_dia(I,j,nz)) * Idt - enddo - endif - enddo - if (CS%debug) then - call MOM_state_chksum("aft 1st loop tridiag ", u, v, h, G, GV, haloshift=0) - endif - !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval) - do J=Jsq,Jeq - do i=is,ie - if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,1) = v(i,J,1) - hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect - b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1))) - d1(I) = hval * b1(I) - v(i,J,1) = b1(i) * (hval * v(i,J,1)) - enddo - do k=2,nz ; do i=is,ie - if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,k) = v(i,J,k) - c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i) - eaval = ea(i,j,k) + ea(i,j+1,k) - hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect - b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval)) - d1(i) = (hval + d1(i)*eaval) * b1(i) - v(i,J,k) = (hval*v(i,J,k) + eaval*v(i,J,k-1))*b1(i) - enddo ; enddo - do k=nz-1,1,-1 ; do i=is,ie - v(i,J,k) = v(i,J,k) + c1(i,k+1)*v(i,J,k+1) - if (associated(ADp%dv_dt_dia)) & - ADp%dv_dt_dia(i,J,k) = (v(i,J,k) - ADp%dv_dt_dia(i,J,k)) * Idt - enddo ; enddo - if (associated(ADp%dv_dt_dia)) then - do i=is,ie - ADp%dv_dt_dia(i,J,nz) = (v(i,J,nz)-ADp%dv_dt_dia(i,J,nz)) * Idt - enddo - endif - enddo - call cpu_clock_end(id_clock_tridiag) - if (CS%debug) then - call MOM_state_chksum("after u/v tridiag ", u, v, h, G, GV, haloshift=0) - endif - endif ! useALEalgorithm + call cpu_clock_end(id_clock_pass) call disable_averaging(CS%diag) ! Frazil formation keeps temperature above the freezing point. @@ -1939,8 +1585,18 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, 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) + + if (CS%use_CVMix_ddiff .and. differentialDiffusion) then + call MOM_error(FATAL, 'diabatic_driver_init: '// & + 'Multiple double-diffusion options selected (DOUBLE_DIFFUSION and'//& + 'USE_CVMIX_DDIFF), please disable all but one option to proceed.') + endif + 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"//& diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 9906083597..c40b3b0a2b 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,105 @@ 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 :: double_diffusion !< If true, enable double-diffusive mixing using an old method. + 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. + 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) 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() @@ -157,7 +160,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 @@ -181,10 +183,21 @@ module MOM_set_diffusivity end type diffusivity_diags ! Clocks -integer :: id_clock_kappaShear +integer :: id_clock_kappaShear, id_clock_CVMix_ddiff 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 +209,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,16 +239,16 @@ 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) + 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) real :: I_Rho0 ! inverse of Boussinesq density (m3/kg) @@ -271,10 +284,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. & - .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.") + if ((CS%use_CVMix_ddiff .or. CS%double_diffusion) .and. .not. & + (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S))) & + call MOM_error(FATAL, "set_diffusivity: both visc%Kd_extra_T and "//& + "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF or DOUBLE_DIFFUSION are 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. @@ -341,6 +360,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 +375,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,10 +391,10 @@ 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 + ! Double-diffusion (old method) 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 @@ -401,6 +422,14 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & enddo ; enddo ; endif endif + ! Apply double diffusion via CVMix + ! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available. + if (CS%use_CVMix_ddiff) then + call cpu_clock_begin(id_clock_CVMix_ddiff) + call compute_ddiff_coeffs(h, tv, G, GV, j, visc%Kd_extra_T, visc%Kd_extra_S, CS%CVMix_ddiff_csp) + call cpu_clock_end(id_clock_CVMix_ddiff) + endif + ! Add the input turbulent diffusivity. if (CS%useKappaShear .or. CS%use_CVMix_shear) then if (present(Kd_int)) then @@ -496,6 +525,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 +546,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 +566,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 @@ -952,8 +995,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 @@ -984,27 +1025,6 @@ subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, Kd_T_dd, Kd_S_dd) 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) @@ -1065,6 +1085,7 @@ subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, Kd_T_dd, Kd_S_dd) 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 +1995,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,12 +2102,11 @@ 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.", & @@ -2130,6 +2155,11 @@ 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) + if (CS%use_CVMix_ddiff) & + id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=CLOCK_MODULE) + end subroutine set_diffusivity_init !> Clear pointers and dealocate memory @@ -2146,6 +2176,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..33a7fbaa4f 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, differential_diffusion 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,9 @@ 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 + differential_diffusion = .false. 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"//& @@ -1923,7 +1890,9 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) "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 +1985,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 +2043,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 .or. differential_diffusion) 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 +2091,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)