diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e96a3807a7..9b70b81415 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2378,9 +2378,6 @@ 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_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 38eb78b89a..ef40f0170c 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1814,10 +1814,10 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) J=segment%HI%JsdB allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) do k=1,nz - rx_tangential(segment%HI%IsdB,J,k) = segment%ry_normal(segment%HI%isd,J,k) - rx_tangential(segment%HI%IedB,J,k) = segment%ry_normal(segment%HI%ied,J,k) + rx_tangential(segment%HI%IsdB,J,k) = segment%rx_normal(segment%HI%isd,J,k) + rx_tangential(segment%HI%IedB,J,k) = segment%rx_normal(segment%HI%ied,J,k) do I=segment%HI%IsdB+1,segment%HI%IedB-1 - rx_tangential(I,J,k) = 0.5*(segment%ry_normal(i,J,k) + segment%ry_normal(i+1,J,k)) + rx_tangential(I,J,k) = 0.5*(segment%rx_normal(i,J,k) + segment%rx_normal(i+1,J,k)) enddo enddo if (segment%radiation_tan) then @@ -1925,10 +1925,10 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) J=segment%HI%JsdB allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) do k=1,nz - rx_tangential(segment%HI%IsdB,J,k) = segment%ry_normal(segment%HI%isd,J,k) - rx_tangential(segment%HI%IedB,J,k) = segment%ry_normal(segment%HI%ied,J,k) + rx_tangential(segment%HI%IsdB,J,k) = segment%rx_normal(segment%HI%isd,J,k) + rx_tangential(segment%HI%IedB,J,k) = segment%rx_normal(segment%HI%ied,J,k) do I=segment%HI%IsdB+1,segment%HI%IedB-1 - rx_tangential(I,J,k) = 0.5*(segment%ry_normal(i,J,k) + segment%ry_normal(i+1,J,k)) + rx_tangential(I,J,k) = 0.5*(segment%rx_normal(i,J,k) + segment%rx_normal(i+1,J,k)) enddo enddo if (segment%radiation_tan) then diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 02b0b622a3..09305eb9fb 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -233,9 +233,6 @@ 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 57b86c80ca..2be8beee4a 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -212,7 +212,6 @@ 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 deleted file mode 100644 index 7137aabfa6..0000000000 --- a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 +++ /dev/null @@ -1,301 +0,0 @@ -!> 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 53339d3488..2635af7fb5 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 !< Exponent of unitless factor of diff. - !! for KPP internal shear mixing scheme. + real :: KPP_exp !< real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) real, allocatable, dimension(:,:,:) :: S2 !< Squared shear frequency (1/s2) real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number +! real, allocatable, dimension(:,:,:) :: kv !< vertical viscosity at interface (m2/s) +! real, allocatable, dimension(:,:,:) :: kd !< vertical diffusivity at interface (m2/s) character(10) :: Mix_Scheme !< Mixing scheme name (string) ! Daignostic handles and pointers type(diag_ctrl), pointer :: diag => NULL() @@ -52,26 +52,25 @@ 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, dummy + real :: gorho + real :: pref, DU, DV, DRHO, DZ, N2, S2 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 @@ -121,30 +120,10 @@ 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,:), & @@ -230,11 +209,7 @@ 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 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, & + 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 79234c7e11..f98185685a 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(inout) :: BLD!< bnd. layer depth (m) + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: BLD!< bnd. layer depth (m) ! Local variables integer :: i,j do j = G%jsc, G%jec ; do i = G%isc, G%iec diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index bb1e0b11c1..61c212db8b 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) = 0.5*(kd_lay(i,j,K-1) + kd_lay(i,j,K)) + CS%kd_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) + 0.5*(kd_lay(i,j,K-1) + kd_lay(i,j,K)) CS%kv_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) * CS%prandtl_bkgnd enddo enddo diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 528dc33135..6eb3b854f4 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -2,6 +2,53 @@ 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 @@ -192,19 +239,26 @@ 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 !< 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. + 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. - ! 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. @@ -291,25 +345,30 @@ 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 !< 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 . + 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 . integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -351,29 +410,33 @@ 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 !< 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) ?? + 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 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 @@ -476,9 +539,10 @@ 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 @@ -515,22 +579,37 @@ 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, 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, 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 :: 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. @@ -1239,20 +1318,26 @@ 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 !< 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(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(vardesc) :: vd ! This "include" declares and sets the variable "version". @@ -1375,48 +1460,4 @@ 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 698243a7f6..6316fd40e6 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -10,7 +10,6 @@ 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 @@ -98,7 +97,6 @@ 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. @@ -251,7 +249,7 @@ module MOM_diabatic_driver integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap -integer :: id_clock_kpp, id_clock_CVMix_ddiff +integer :: id_clock_kpp contains @@ -387,6 +385,7 @@ 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") @@ -488,13 +487,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif if (CS%ML_mix_first > 0.0) then - ! This subroutine: - ! (1) Cools the mixed layer. - ! (2) Performs convective adjustment by mixed layer entrainment. - ! (3) Heats the mixed layer and causes it to detrain to - ! Monin-Obukhov depth or minimum mixed layer depth. - ! (4) Uses any remaining TKE to drive mixed layer entrainment. - ! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) +! This subroutine +! (1) Cools the mixed layer. +! (2) Performs convective adjustment by mixed layer entrainment. +! (3) Heats the mixed layer and causes it to detrain to +! Monin-Obukhov depth or minimum mixed layer depth. +! (4) Uses any remaining TKE to drive mixed layer entrainment. +! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) call find_uv_at_h(u, v, h, u_h, v_h, G, GV) call cpu_clock_begin(id_clock_mixedlayer) @@ -529,12 +528,11 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)") if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G) endif - endif ! end CS%bulkmixedlayer + endif if (CS%debug) then call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, haloshift=0) endif - if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then call find_uv_at_h(u, v, h_orig, u_h, v_h, G, GV, eaml, ebml) @@ -591,7 +589,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 ! end CS%use_int_tides + endif call cpu_clock_begin(id_clock_set_diffusivity) ! Sets: Kd, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S @@ -732,10 +730,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! a diffusivity and happen before KPP. But generally in MOM, we do not match ! KPP boundary layer to interior, so this diffusivity can be computed when convenient. if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T)) then - call cpu_clock_begin(id_clock_CVMix_ddiff) + 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_CVMix_ddiff) + 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) @@ -748,6 +746,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo endif + endif @@ -1382,9 +1381,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! visc%Kv_shear is not in the group pass because it has larger vertical extent. if (associated(visc%Kv_shear)) & call pass_var(visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) - if (associated(visc%Kv_slow)) & - call pass_var(visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) - call cpu_clock_end(id_clock_pass) if (.not. CS%useALEalgorithm) then @@ -1889,7 +1885,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, real :: Kd integer :: num_mode - logical :: use_temperature + logical :: use_temperature, differentialDiffusion type(vardesc) :: vd ! This "include" declares and sets the variable "version". @@ -1940,10 +1936,11 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, the diffusivity from ePBL is added to all\n"//& "other diffusivities. Otherwise, the larger of kappa-\n"//& "shear and ePBL diffusivities are used.", default=.true.) - CS%use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) + call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialDiffusion, & + "If true, apply parameterization of double-diffusion.", & + default=.false. ) CS%use_kappa_shear = kappa_shear_is_used(param_file) CS%use_CVMix_shear = CVMix_shear_is_used(param_file) - 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"//& @@ -2402,8 +2399,8 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=CLOCK_MODULE) id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=CLOCK_ROUTINE) id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=CLOCK_ROUTINE) - id_clock_CVMix_ddiff = -1 ; if (CS%use_CVMix_ddiff) & - id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion)', grain=CLOCK_ROUTINE) + id_clock_differential_diff = -1 ; if (differentialDiffusion) & + id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=CLOCK_ROUTINE) ! initialize the auxiliary diabatic driver module call diabatic_aux_init(Time, G, GV, param_file, diag, CS%diabatic_aux_CSp, & diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 903868795a..9906083597 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -23,8 +23,6 @@ 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 @@ -45,101 +43,104 @@ 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 :: 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. + 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. 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 :: 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. + 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) 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,6 +158,11 @@ module MOM_set_diffusivity integer :: id_N2 = -1 integer :: id_N2_z = -1 + integer :: id_KT_extra = -1 + integer :: id_KS_extra = -1 + integer :: id_KT_extra_z = -1 + integer :: id_KS_extra_z = -1 + end type set_diffusivity_CS type diffusivity_diags @@ -166,9 +172,12 @@ module MOM_set_diffusivity Kd_BBL => NULL(),& ! BBL diffusivity at interfaces (m2/s) Kd_work => NULL(),& ! layer integrated work by diapycnal mixing (W/m2) maxTKE => NULL(),& ! energy required to entrain to h_max (m3/s3) - TKE_to_Kd => NULL() ! conversion rate (~1.0 / (G_Earth + dRho_lay)) + TKE_to_Kd => NULL(),& ! conversion rate (~1.0 / (G_Earth + dRho_lay)) ! between TKE dissipated within a layer and Kd ! in that layer, in m2 s-1 / m3 s-3 = s2 m-1 + KT_extra => NULL(),& ! double diffusion diffusivity for temp (m2/s) + KS_extra => NULL() ! double diffusion diffusivity for saln (m2/s) + end type diffusivity_diags ! Clocks @@ -176,17 +185,6 @@ module MOM_set_diffusivity contains -!> Sets the interior vertical diffusion of scalars due to the following processes: -!! 1) Shear-driven mixing: two options, Jackson et at. and KPP interior; -!! 2) Background mixing via CVMix (Bryan-Lewis profile) or the scheme described by -!! Harrison & Hallberg, JPO 2008; -!! 3) Double-diffusion 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. @@ -198,9 +196,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 !< zonal thickness transport m^2/s. + intent(in) :: u_h real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: v_h !< meridional thickness transport m^2/s. + intent(in) :: v_h 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 @@ -228,15 +226,17 @@ 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? + 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) real :: dissip ! local variable for dissipation calculations (W/m3) @@ -271,16 +271,10 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & use_EOS = associated(tv%eqn_of_state) - if ((CS%use_CVMix_ddiff) .and. & + 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 USE_CVMIX_DDIFF is true.") - - ! Set Kd, Kd_int and Kv_slow to constant values. - ! If nothing else is specified, this will be the value used. - Kd(:,:,:) = CS%Kd - Kd_int(:,:,:) = CS%Kd - if (associated(visc%Kv_slow)) visc%Kv_slow(:,:,:) = CS%Kv + "visc%Kd_extra_S must be associated when DOUBLE_DIFFUSION is true.") ! Set up arrays for diagnostics. @@ -299,6 +293,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%id_TKE_to_Kd > 0) then allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz)) ; dd%TKE_to_Kd(:,:,:) = 0.0 endif + if ((CS%id_KT_extra > 0) .or. (CS%id_KT_extra_z > 0)) then + allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1)) ; dd%KT_extra(:,:,:) = 0.0 + endif + if ((CS%id_KS_extra > 0) .or. (CS%id_KS_extra_z > 0)) then + allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1)) ; dd%KS_extra(:,:,:) = 0.0 + endif if ((CS%id_Kd_BBL > 0) .or. (CS%id_Kd_BBL_z > 0)) then allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0 endif @@ -341,10 +341,6 @@ 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 @@ -356,6 +352,8 @@ 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) & @@ -372,13 +370,35 @@ 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) - ! Apply double diffusion - ! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available. - if (CS%use_CVMix_ddiff) then - call compute_ddiff_coeffs(h, tv, G, GV, j, visc%Kd_extra_T, visc%Kd_extra_S, CS%CVMix_ddiff_csp) + ! GMM, the following will go into the MOM_CVMix_double_diffusion module + if (CS%double_diffusion) then + call double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, KT_extra, KS_extra) + do K=2,nz ; do i=is,ie + if (KS_extra(i,K) > KT_extra(i,K)) then ! salt fingering + Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KT_extra(i,K) + Kd(i,j,k) = Kd(i,j,k) + 0.5*KT_extra(i,K) + visc%Kd_extra_S(i,j,k) = KS_extra(i,K)-KT_extra(i,K) + visc%Kd_extra_T(i,j,k) = 0.0 + elseif (KT_extra(i,K) > 0.0) then ! double-diffusive convection + Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KS_extra(i,K) + Kd(i,j,k) = Kd(i,j,k) + 0.5*KS_extra(i,K) + visc%Kd_extra_T(i,j,k) = KT_extra(i,K)-KS_extra(i,K) + visc%Kd_extra_S(i,j,k) = 0.0 + else ! There is no double diffusion at this interface. + visc%Kd_extra_T(i,j,k) = 0.0 + visc%Kd_extra_S(i,j,k) = 0.0 + endif + enddo ; enddo + if (associated(dd%KT_extra)) then ; do K=1,nz+1 ; do i=is,ie + dd%KT_extra(i,j,K) = KT_extra(i,K) + enddo ; enddo ; endif + + if (associated(dd%KS_extra)) then ; do K=1,nz+1 ; do i=is,ie + dd%KS_extra(i,j,K) = KS_extra(i,K) + enddo ; enddo ; endif endif ! Add the input turbulent diffusivity. @@ -476,11 +496,6 @@ 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.) @@ -497,6 +512,12 @@ 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) @@ -517,28 +538,13 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & T_f, S_f, dd%Kd_user) endif - ! post diagnostics - - ! 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) - + ! GMM, post diags... if (CS%id_Kd_layer > 0) call post_data(CS%id_Kd_layer, Kd, CS%diag) - ! tidal mixing + num_z_diags = 0 + 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 @@ -562,11 +568,26 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif + if (CS%id_KT_extra > 0) call post_data(CS%id_KT_extra, dd%KT_extra, CS%diag) + if (CS%id_KS_extra > 0) call post_data(CS%id_KS_extra, dd%KS_extra, CS%diag) if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, dd%Kd_BBL, CS%diag) + if (CS%id_KT_extra_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_KT_extra_z + z_ptrs(num_z_diags)%p => dd%KT_extra + endif + + if (CS%id_KS_extra_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_KS_extra_z + z_ptrs(num_z_diags)%p => dd%KS_extra + endif + if (CS%id_Kd_BBL_z > 0) then num_z_diags = num_z_diags + 1 z_ids(num_z_diags) = CS%id_Kd_BBL_z + z_ptrs(num_z_diags)%p => dd%KS_extra endif if (num_z_diags > 0) & @@ -577,6 +598,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (associated(dd%Kd_user)) deallocate(dd%Kd_user) if (associated(dd%maxTKE)) deallocate(dd%maxTKE) if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd) + if (associated(dd%KT_extra)) deallocate(dd%KT_extra) + if (associated(dd%KS_extra)) deallocate(dd%KS_extra) if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL) if (showCallTree) call callTree_leave("set_diffusivity()") @@ -929,6 +952,119 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & end subroutine find_N2 +! GMM, the following will be moved to a new module + +!> This subroutine sets the additional diffusivities of temperature and +!! salinity due to double diffusion, using the same functional form as is +!! used in MOM4.1, and taken from an NCAR technical note (REF?) that updates +!! what was in Large et al. (1994). All the coefficients here should probably +!! be made run-time variables rather than hard-coded constants. +!! +!! \todo Find reference for NCAR tech note above. +subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, Kd_T_dd, Kd_S_dd) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available + !! thermodynamic fields; absent fields have NULL + !! ptrs. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: T_f !< layer temp in C with the values in massless layers + !! filled vertically by diffusion. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: S_f !< Layer salinities in PPT with values in massless + !! layers filled vertically by diffusion. + integer, intent(in) :: j !< Meridional index upon which to work. + type(set_diffusivity_CS), pointer :: CS !< Module control structure. + real, dimension(SZI_(G),SZK_(G)+1), & + intent(out) :: Kd_T_dd !< Interface double diffusion diapycnal + !! diffusivity for temp (m2/sec). + real, dimension(SZI_(G),SZK_(G)+1), & + intent(out) :: Kd_S_dd !< Interface double diffusion diapycnal + !! diffusivity for saln (m2/sec). + +! Arguments: +! (in) tv - structure containing pointers to any available +! thermodynamic fields; absent fields have NULL ptrs +! (in) h - layer thickness (m or kg m-2) +! (in) T_f - layer temp in C with the values in massless layers +! filled vertically by diffusion +! (in) S_f - layer salinities in PPT with values in massless layers +! filled vertically by diffusion +! (in) G - ocean grid structure +! (in) GV - The ocean's vertical grid structure. +! (in) CS - module control structure +! (in) j - meridional index upon which to work +! (out) Kd_T_dd - interface double diffusion diapycnal diffusivity for temp (m2/sec) +! (out) Kd_S_dd - interface double diffusion diapycnal diffusivity for saln (m2/sec) + +! This subroutine sets the additional diffusivities of temperature and +! salinity due to double diffusion, using the same functional form as is +! used in MOM4.1, and taken from an NCAR technical note (###REF?) that updates +! what was in Large et al. (1994). All the coefficients here should probably +! be made run-time variables rather than hard-coded constants. + + real, dimension(SZI_(G)) :: & + dRho_dT, & ! partial derivatives of density wrt temp (kg m-3 degC-1) + dRho_dS, & ! partial derivatives of density wrt saln (kg m-3 PPT-1) + pres, & ! pressure at each interface (Pa) + Temp_int, & ! temp and saln at interfaces + Salin_int + + real :: alpha_dT ! density difference between layers due to temp diffs (kg/m3) + real :: beta_dS ! density difference between layers due to saln diffs (kg/m3) + + real :: Rrho ! vertical density ratio + real :: diff_dd ! factor for double-diffusion + real :: prandtl ! flux ratio for diffusive convection regime + + real, parameter :: Rrho0 = 1.9 ! limit for double-diffusive density ratio + real, parameter :: dsfmax = 1.e-4 ! max diffusivity in case of salt fingering + real, parameter :: Kv_molecular = 1.5e-6 ! molecular viscosity (m2/sec) + + integer :: i, k, is, ie, nz + is = G%isc ; ie = G%iec ; nz = G%ke + + if (associated(tv%eqn_of_state)) then + do i=is,ie + pres(i) = 0.0 ; Kd_T_dd(i,1) = 0.0 ; Kd_S_dd(i,1) = 0.0 + Kd_T_dd(i,nz+1) = 0.0 ; Kd_S_dd(i,nz+1) = 0.0 + enddo + do K=2,nz + do i=is,ie + pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1) + Temp_Int(i) = 0.5 * (T_f(i,j,k-1) + T_f(i,j,k)) + Salin_Int(i) = 0.5 * (S_f(i,j,k-1) + S_f(i,j,k)) + enddo + call calculate_density_derivs(Temp_int, Salin_int, pres, & + dRho_dT(:), dRho_dS(:), is, ie-is+1, tv%eqn_of_state) + + do i=is,ie + alpha_dT = -1.0*dRho_dT(i) * (T_f(i,j,k-1) - T_f(i,j,k)) + beta_dS = dRho_dS(i) * (S_f(i,j,k-1) - S_f(i,j,k)) + + if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then ! salt finger case + Rrho = min(alpha_dT/beta_dS,Rrho0) + diff_dd = 1.0 - ((RRho-1.0)/(RRho0-1.0)) + diff_dd = dsfmax*diff_dd*diff_dd*diff_dd + Kd_T_dd(i,K) = 0.7*diff_dd + Kd_S_dd(i,K) = diff_dd + elseif ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then ! diffusive convection + Rrho = alpha_dT/beta_dS + diff_dd = Kv_molecular*0.909*exp(4.6*exp(-0.54*(1/Rrho-1))) + prandtl = 0.15*Rrho + if (Rrho > 0.5) prandtl = (1.85-0.85/Rrho)*Rrho + Kd_T_dd(i,K) = diff_dd + Kd_S_dd(i,K) = prandtl*diff_dd + else + Kd_T_dd(i,K) = 0.0 ; Kd_S_dd(i,K) = 0.0 + endif + enddo + enddo + endif + +end subroutine double_diffusion !> This routine adds diffusion sustained by flow energy extracted by bottom drag. subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & maxTKE, kb, G, GV, CS, Kd, Kd_int, Kd_BBL) @@ -1838,11 +1974,6 @@ 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"//& @@ -1945,6 +2076,45 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp endif endif + + ! GMM, the following should be moved to the DD module + call get_param(param_file, mdl, "DOUBLE_DIFFUSION", CS%double_diffusion, & + "If true, increase diffusivitives for temperature or salt \n"//& + "based on double-diffusive paramaterization from MOM4/KPP.", & + default=.false.) + if (CS%double_diffusion) then + call get_param(param_file, mdl, "MAX_RRHO_SALT_FINGERS", CS%Max_Rrho_salt_fingers, & + "Maximum density ratio for salt fingering regime.", & + default=2.55, units="nondim") + call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", CS%Max_salt_diff_salt_fingers, & + "Maximum salt diffusivity for salt fingering regime.", & + default=1.e-4, units="m2 s-1") + call get_param(param_file, mdl, "KV_MOLECULAR", CS%Kv_molecular, & + "Molecular viscosity for calculation of fluxes under \n"//& + "double-diffusive convection.", default=1.5e-6, units="m2 s-1") + ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults. + + CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for temperature', 'm2 s-1') + + CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for salinity', 'm2 s-1') + + if (associated(diag_to_Z_CSp)) then + vd = var_desc("KT_extra", "m2 s-1", & + "Double-Diffusive Temperature Diffusivity, interpolated to z", & + z_grid='z') + CS%id_KT_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + vd = var_desc("KS_extra", "m2 s-1", & + "Double-Diffusive Salinity Diffusivity, interpolated to z",& + z_grid='z') + CS%id_KS_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + vd = var_desc("Kd_BBL", "m2 s-1", & + "Bottom Boundary Layer Diffusivity", z_grid='z') + CS%id_Kd_BBL_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + endif + endif + if (CS%user_change_diff) then call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) endif @@ -1960,9 +2130,6 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp ! CVMix shear-driven mixing CS%use_CVMix_shear = CVMix_shear_init(Time, G, GV, param_file, CS%diag, CS%CVMix_shear_csp) - ! CVMix double diffusion mixing - CS%use_CVMix_ddiff = CVMix_ddiff_init(Time, G, GV, param_file, CS%diag, CS%CVMix_ddiff_csp) - end subroutine set_diffusivity_init !> Clear pointers and dealocate memory @@ -1979,9 +2146,6 @@ 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 ec1b09a5ad..ec0b5a80b3 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -2,6 +2,38 @@ 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 @@ -12,9 +44,8 @@ module MOM_set_visc use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_kappa_shear, only : kappa_shear_is_used -use MOM_cvmix_shear, only : cvmix_shear_is_used -use MOM_cvmix_conv, only : cvmix_conv_is_used -use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used +use MOM_CVMix_shear, only : CVMix_shear_is_used +use MOM_CVMix_conv, only : CVMix_conv_is_used use MOM_io, only : vardesc, var_desc use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_variables, only : thermo_var_ptrs @@ -1760,10 +1791,8 @@ 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) @@ -1782,9 +1811,7 @@ 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 - - ! 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 + 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') @@ -1827,14 +1854,21 @@ 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 - - ! local variables +! 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 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, use_omega - logical :: use_CVMix_ddiff + logical :: use_kappa_shear, adiabatic, differential_diffusion, use_omega type(OBC_segment_type), pointer :: segment ! pointer to OBC segment type ! This include declares and sets the variable "version". #include "version_variable.h" @@ -1857,8 +1891,8 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) ! Set default, read and log parameters call log_version(param_file, mdl, version, "") - CS%RiNo_mix = .false. ; use_CVMix_ddiff = .false. - use_kappa_shear = .false. !; adiabatic = .false. ! Needed? -AJA + CS%RiNo_mix = .false. + use_kappa_shear = .false. ; differential_diffusion = .false. !; adiabatic = .false. ! Needed? -AJA call get_param(param_file, mdl, "BOTTOMDRAGLAW", CS%bottomdraglaw, & "If true, the bottom stress is calculated with a drag \n"//& "law of the form c_drag*|u|*u. The velocity magnitude \n"//& @@ -1885,9 +1919,11 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) if (.not.adiabatic) then use_kappa_shear = kappa_shear_is_used(param_file) CS%RiNo_mix = use_kappa_shear - use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) + call get_param(param_file, mdl, "DOUBLE_DIFFUSION", differential_diffusion, & + "If true, increase diffusivitives for temperature or salt \n"//& + "based on double-diffusive paramaterization from MOM4/KPP.", & + default=.false.) 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) @@ -1980,15 +2016,6 @@ 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) @@ -2038,7 +2065,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 (use_CVMix_ddiff) then + if (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 @@ -2086,37 +2113,4 @@ 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 bafbe5eb59..48a6380ead 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,7 +116,6 @@ 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() @@ -615,8 +614,6 @@ 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. @@ -649,14 +646,6 @@ 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 @@ -832,13 +821,6 @@ 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 @@ -1002,14 +984,6 @@ 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 @@ -1023,9 +997,6 @@ 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) @@ -1194,44 +1165,6 @@ 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. @@ -1738,30 +1671,17 @@ 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)