diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 8ee4113f71..e96a3807a7 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -52,6 +52,7 @@ module MOM use MOM_ALE, only : ALE_getCoordinate, ALE_getCoordinateUnits, ALE_writeCoordinateFile use MOM_ALE, only : ALE_updateVerticalGridType, ALE_remap_init_conds, ALE_register_diags use MOM_boundary_update, only : call_OBC_register, OBC_register_end, update_OBC_CS +use MOM_legacy_diabatic_driver,only : legacy_diabatic use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS use MOM_diabatic_driver, only : adiabatic, adiabatic_driver_init, diabatic_driver_end use MOM_diagnostics, only : calculate_diagnostic_fields, MOM_diagnostics_init @@ -197,6 +198,9 @@ module MOM logical :: adiabatic !< If true, there are no diapycnal mass fluxes, and no calls !! to routines to calculate or apply diapycnal fluxes. + logical :: use_legacy_diabatic_driver!< If true (default), use the a legacy version of the + !! diabatic subroutine. This is temporary and is needed + !! to avoid change in answers. logical :: diabatic_first !< If true, apply diabatic and thermodynamic !! processes before time stepping the dynamics. logical :: use_ALE_algorithm !< If true, use the ALE algorithm rather than layered @@ -1151,8 +1155,14 @@ subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, & endif call cpu_clock_begin(id_clock_diabatic) - call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, & - dtdia, Time_end_thermo, G, GV, CS%diabatic_CSp, Waves=Waves) + if (CS%use_legacy_diabatic_driver) then + ! the following subroutine is legacy and will be deleted in the near future. + call legacy_diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, & + dtdia, Time_end_thermo, G, GV, CS%diabatic_CSp, Waves=Waves) + else + call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, & + dtdia, Time_end_thermo, G, GV, CS%diabatic_CSp, Waves=Waves) + endif fluxes%fluxes_used = .true. call cpu_clock_end(id_clock_diabatic) @@ -1627,6 +1637,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "true. This assumes that KD = KDML = 0.0 and that \n"//& "there is no buoyancy forcing, but makes the model \n"//& "faster by eliminating subroutine calls.", default=.false.) + call get_param(param_file, "MOM", "USE_LEGACY_DIABATIC_DRIVER", CS%use_legacy_diabatic_driver, & + "If true, use the a legacy version of the diabatic subroutine. \n"//& + "This is temporary and is needed avoid change in answers.", & + default=.true.) call get_param(param_file, "MOM", "DO_DYNAMICS", CS%do_dynamics, & "If False, skips the dynamics calls that update u & v, as well as \n"//& "the gravity wave adjustment to h. This is a fragile feature and \n"//& @@ -2364,6 +2378,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (associated(CS%visc%Kv_shear)) & call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) + if (associated(CS%visc%Kv_slow)) & + call pass_var(CS%visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) + call cpu_clock_end(id_clock_pass_init) call register_obsolete_diagnostics(param_file, CS%diag) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 09305eb9fb..02b0b622a3 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -233,6 +233,9 @@ module MOM_variables !! convection etc). TKE_turb => NULL() !< The turbulent kinetic energy per unit mass defined !! at the interfaces between each layer, in m2 s-2. + logical :: add_Kv_slow !< If True, adds Kv_slow when calculating the + !! 'coupling coefficient' (a[k]) at the interfaces. + !! This is done in find_coupling_coef. end type vertvisc_type !> The BT_cont_type structure contains information about the summed layer diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index 2be8beee4a..57b86c80ca 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -212,6 +212,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, CS, hbl) iFaceHeight(k+1) = iFaceHeight(k) - dh enddo + ! gets index of the level and interface above hbl kOBL = CVMix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) call CVMix_coeffs_conv(Mdiff_out=CS%kv_conv(i,j,:), & diff --git a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 new file mode 100644 index 0000000000..7137aabfa6 --- /dev/null +++ b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 @@ -0,0 +1,301 @@ +!> Interface to CVMix double diffusion scheme. +module MOM_CVMix_ddiff + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field +use MOM_diag_mediator, only : post_data +use MOM_EOS, only : calculate_density_derivs +use MOM_variables, only : thermo_var_ptrs +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_debugging, only : hchksum +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_file_parser, only : get_param, log_version, param_file_type +use cvmix_ddiff, only : cvmix_init_ddiff, CVMix_coeffs_ddiff +use cvmix_kpp, only : CVmix_kpp_compute_kOBL_depth +implicit none ; private + +#include + +public CVMix_ddiff_init, CVMix_ddiff_end, CVMix_ddiff_is_used, compute_ddiff_coeffs + +!> Control structure including parameters for CVMix double diffusion. +type, public :: CVMix_ddiff_cs + + ! Parameters + real :: strat_param_max !< maximum value for the stratification parameter (nondim) + real :: kappa_ddiff_s !< leading coefficient in formula for salt-fingering regime + !! for salinity diffusion (m^2/s) + real :: ddiff_exp1 !< interior exponent in salt-fingering regime formula (nondim) + real :: ddiff_exp2 !< exterior exponent in salt-fingering regime formula (nondim) + real :: mol_diff !< molecular diffusivity (m^2/s) + real :: kappa_ddiff_param1 !< exterior coefficient in diffusive convection regime (nondim) + real :: kappa_ddiff_param2 !< middle coefficient in diffusive convection regime (nondim) + real :: kappa_ddiff_param3 !< interior coefficient in diffusive convection regime (nondim) + real :: min_thickness !< Minimum thickness allowed (m) + character(len=4) :: diff_conv_type !< type of diffusive convection to use. Options are Marmorino & + !! Caldwell 1976 ("MC76"; default) and Kelley 1988, 1990 ("K90") + logical :: debug !< If true, turn on debugging + + ! Daignostic handles and pointers + type(diag_ctrl), pointer :: diag => NULL() + integer :: id_KT_extra = -1, id_KS_extra = -1, id_R_rho = -1 + + ! Diagnostics arrays + real, allocatable, dimension(:,:,:) :: KT_extra !< double diffusion diffusivity for temp (m2/s) + real, allocatable, dimension(:,:,:) :: KS_extra !< double diffusion diffusivity for salt (m2/s) + real, allocatable, dimension(:,:,:) :: R_rho !< double-diffusion density ratio (nondim) + +end type CVMix_ddiff_cs + +character(len=40) :: mdl = "MOM_CVMix_ddiff" !< This module's name. + +contains + +!> Initialized the CVMix double diffusion module. +logical function CVMix_ddiff_init(Time, G, GV, param_file, diag, CS) + + type(time_type), intent(in) :: Time !< The current time. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. + type(CVMix_ddiff_cs), pointer :: CS !< This module's control structure. + +! This include declares and sets the variable "version". +#include "version_variable.h" + + if (associated(CS)) then + call MOM_error(WARNING, "CVMix_ddiff_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + ! Read parameters + call log_version(param_file, mdl, version, & + "Parameterization of mixing due to double diffusion processes via CVMix") + call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_init, & + "If true, turns on double diffusive processes via CVMix. \n"// & + "Note that double diffusive processes on viscosity are ignored \n"// & + "in CVMix, see http://cvmix.github.io/ for justification.",& + default=.false.) + + if (.not. CVMix_ddiff_init) return + + call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) + + call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, do_not_log=.True.) + + call openParameterBlock(param_file,'CVMIX_DDIFF') + + call get_param(param_file, mdl, "STRAT_PARAM_MAX", CS%strat_param_max, & + "The maximum value for the double dissusion stratification parameter", & + units="nondim", default=2.55) + + call get_param(param_file, mdl, "KAPPA_DDIFF_S", CS%kappa_ddiff_s, & + "Leading coefficient in formula for salt-fingering regime \n"// & + "for salinity diffusion.", units="m2 s-1", default=1.0e-4) + + call get_param(param_file, mdl, "DDIFF_EXP1", CS%ddiff_exp1, & + "Interior exponent in salt-fingering regime formula.", & + units="nondim", default=1.0) + + call get_param(param_file, mdl, "DDIFF_EXP2", CS%ddiff_exp2, & + "Exterior exponent in salt-fingering regime formula.", & + units="nondim", default=3.0) + + call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM1", CS%kappa_ddiff_param1, & + "Exterior coefficient in diffusive convection regime.", & + units="nondim", default=0.909) + + call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM2", CS%kappa_ddiff_param2, & + "Middle coefficient in diffusive convection regime.", & + units="nondim", default=4.6) + + call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM3", CS%kappa_ddiff_param3, & + "Interior coefficient in diffusive convection regime.", & + units="nondim", default=-0.54) + + call get_param(param_file, mdl, "MOL_DIFF", CS%mol_diff, & + "Molecular diffusivity used in CVMix double diffusion.", & + units="m2 s-1", default=1.5e-6) + + call get_param(param_file, mdl, "DIFF_CONV_TYPE", CS%diff_conv_type, & + "type of diffusive convection to use. Options are Marmorino \n" //& + "and Caldwell 1976 (MC76) and Kelley 1988, 1990 (K90).", & + default="MC76") + + call closeParameterBlock(param_file) + + ! Register diagnostics + CS%diag => diag + + CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for temperature', 'm2 s-1') + + CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for salinity', 'm2 s-1') + + CS%id_R_rho = register_diag_field('ocean_model','R_rho',diag%axesTi,Time, & + 'Double-diffusion density ratio', 'nondim') + if (CS%id_R_rho > 0) & + allocate(CS%R_rho( SZI_(G), SZJ_(G), SZK_(G)+1)); CS%R_rho(:,:,:) = 0.0 + + call cvmix_init_ddiff(strat_param_max=CS%strat_param_max, & + kappa_ddiff_s=CS%kappa_ddiff_s, & + ddiff_exp1=CS%ddiff_exp1, & + ddiff_exp2=CS%ddiff_exp2, & + mol_diff=CS%mol_diff, & + kappa_ddiff_param1=CS%kappa_ddiff_param1, & + kappa_ddiff_param2=CS%kappa_ddiff_param2, & + kappa_ddiff_param3=CS%kappa_ddiff_param3, & + diff_conv_type=CS%diff_conv_type) + +end function CVMix_ddiff_init + +!> Subroutine for computing vertical diffusion coefficients for the +!! double diffusion mixing parameterization. +subroutine compute_ddiff_coeffs(h, tv, G, GV, j, Kd_T, Kd_S, CS) + + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_T !< Interface double diffusion diapycnal + !! diffusivity for temp (m2/sec). + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_S !< Interface double diffusion diapycnal + !! diffusivity for salt (m2/sec). + type(CVMix_ddiff_cs), pointer :: CS !< The control structure returned + !! by a previous call to CVMix_ddiff_init. + integer, intent(in) :: j !< Meridional grid indice. +! real, dimension(:,:), optional, pointer :: hbl !< Depth of ocean boundary layer (m) + + ! local variables + real, dimension(SZK_(G)) :: & + cellHeight, & !< Height of cell centers (m) + dRho_dT, & !< partial derivatives of density wrt temp (kg m-3 degC-1) + dRho_dS, & !< partial derivatives of density wrt saln (kg m-3 PPT-1) + pres_int, & !< pressure at each interface (Pa) + temp_int, & !< temp and at interfaces (degC) + salt_int, & !< salt at at interfaces + alpha_dT, & !< alpha*dT across interfaces + beta_dS, & !< beta*dS across interfaces + dT, & !< temp. difference between adjacent layers (degC) + dS !< salt difference between adjacent layers + + real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) + integer :: kOBL !< level of OBL extent + real :: pref, g_o_rho0, rhok, rhokm1, dz, dh, hcorr + integer :: i, k + + ! initialize dummy variables + pres_int(:) = 0.0; temp_int(:) = 0.0; salt_int(:) = 0.0 + alpha_dT(:) = 0.0; beta_dS(:) = 0.0; dRho_dT(:) = 0.0 + dRho_dS(:) = 0.0; dT(:) = 0.0; dS(:) = 0.0 + + ! set Kd_T and Kd_S to zero to avoid passing values from previous call + Kd_T(:,j,:) = 0.0; Kd_S(:,j,:) = 0.0 + + ! GMM, I am leaving some code commented below. We need to pass BLD to + ! this soubroutine to avoid adding diffusivity above that. This needs + ! to be done once we re-structure the order of the calls. + !if (.not. associated(hbl)) then + ! allocate(hbl(SZI_(G), SZJ_(G))); + ! hbl(:,:) = 0.0 + !endif + + do i = G%isc, G%iec + + ! skip calling at land points + if (G%mask2dT(i,j) == 0.) cycle + + pRef = 0. + pres_int(1) = pRef + ! we don't have SST and SSS, so let's use values at top-most layer + temp_int(1) = TV%T(i,j,1); salt_int(1) = TV%S(i,j,1) + do k=2,G%ke + ! pressure at interface + pres_int(k) = pRef + GV%H_to_Pa * h(i,j,k-1) + ! temp and salt at interface + ! for temp: (t1*h1 + t2*h2)/(h1+h2) + temp_int(k) = (TV%T(i,j,k-1)*h(i,j,k-1) + TV%T(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k)) + salt_int(k) = (TV%S(i,j,k-1)*h(i,j,k-1) + TV%S(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k)) + ! dT and dS + dT(k) = (TV%T(i,j,k-1)-TV%T(i,j,k)) + dS(k) = (TV%S(i,j,k-1)-TV%S(i,j,k)) + pRef = pRef + GV%H_to_Pa * h(i,j,k-1) + enddo ! k-loop finishes + + call calculate_density_derivs(temp_int(:), salt_int(:), pres_int(:), drho_dT(:), drho_dS(:), 1, G%ke, TV%EQN_OF_STATE) + + ! The "-1.0" below is needed so that the following criteria is satisfied: + ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then "salt finger" + ! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then "diffusive convection" + do k=1,G%ke + alpha_dT(k) = -1.0*drho_dT(k) * dT(k) + beta_dS(k) = drho_dS(k) * dS(k) + enddo + + if (CS%id_R_rho > 0.0) then + do k=1,G%ke + CS%R_rho(i,j,k) = alpha_dT(k)/beta_dS(k) + ! avoid NaN's + if(CS%R_rho(i,j,k) /= CS%R_rho(i,j,k)) CS%R_rho(i,j,k) = 0.0 + enddo + endif + + iFaceHeight(1) = 0.0 ! BBL is all relative to the surface + hcorr = 0.0 + ! compute heights at cell center and interfaces + do k=1,G%ke + dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) + hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 + dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * dh + iFaceHeight(k+1) = iFaceHeight(k) - dh + enddo + + ! gets index of the level and interface above hbl + !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) + + call CVMix_coeffs_ddiff(Tdiff_out=Kd_T(i,j,:), & + Sdiff_out=Kd_S(i,j,:), & + strat_param_num=alpha_dT(:), & + strat_param_denom=beta_dS(:), & + nlev=G%ke, & + max_nlev=G%ke) + + ! Do not apply mixing due to convection within the boundary layer + !do k=1,kOBL + ! Kd_T(i,j,k) = 0.0 + ! Kd_S(i,j,k) = 0.0 + !enddo + + enddo ! i-loop + +end subroutine compute_ddiff_coeffs + +!> Reads the parameter "USE_CVMIX_DDIFF" and returns state. +!! This function allows other modules to know whether this parameterization will +!! be used without needing to duplicate the log entry. +logical function CVMix_ddiff_is_used(param_file) + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_is_used, & + default=.false., do_not_log = .true.) + +end function CVMix_ddiff_is_used + +!> Clear pointers and dealocate memory +subroutine CVMix_ddiff_end(CS) + type(CVMix_ddiff_cs), pointer :: CS ! Control structure + + deallocate(CS) + +end subroutine CVMix_ddiff_end + + +end module MOM_CVMix_ddiff diff --git a/src/parameterizations/vertical/MOM_CVMix_shear.F90 b/src/parameterizations/vertical/MOM_CVMix_shear.F90 index 2635af7fb5..53339d3488 100644 --- a/src/parameterizations/vertical/MOM_CVMix_shear.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_shear.F90 @@ -30,14 +30,14 @@ module MOM_CVMix_shear !> Control structure including parameters for CVMix interior shear schemes. type, public :: CVMix_shear_cs logical :: use_LMD94, use_PP81 !< Flags for various schemes + logical :: smooth_ri !< If true, smooth Ri using a 1-2-1 filter real :: Ri_zero !< LMD94 critical Richardson number real :: Nu_zero !< LMD94 maximum interior diffusivity - real :: KPP_exp !< + real :: KPP_exp !< Exponent of unitless factor of diff. + !! for KPP internal shear mixing scheme. real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) real, allocatable, dimension(:,:,:) :: S2 !< Squared shear frequency (1/s2) real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number -! real, allocatable, dimension(:,:,:) :: kv !< vertical viscosity at interface (m2/s) -! real, allocatable, dimension(:,:,:) :: kd !< vertical diffusivity at interface (m2/s) character(10) :: Mix_Scheme !< Mixing scheme name (string) ! Daignostic handles and pointers type(diag_ctrl), pointer :: diag => NULL() @@ -52,25 +52,26 @@ module MOM_CVMix_shear !> Subroutine for calculating (internal) vertical diffusivities/viscosities subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, & kv, G, GV, CS ) - type(ocean_grid_type), intent(in) :: G !< Grid structure. - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: u_H !< Initial zonal velocity on T points, in m s-1. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: v_H !< Initial meridional velocity on T points, in m s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. - type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd !< The vertical diffusivity at each interface - !! (not layer!) in m2 s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv !< The vertical viscosity at each interface - !! (not layer!) in m2 s-1. - type(CVMix_shear_cs), pointer :: CS !< The control structure returned by a previous call to - !! CVMix_shear_init. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd !< The vertical diffusivity at each interface + !! (not layer!) in m2 s-1. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv !< The vertical viscosity at each interface + !! (not layer!) in m2 s-1. + type(CVMix_shear_cs), pointer :: CS !< The control structure returned by a previous call to + !! CVMix_shear_init. ! Local variables integer :: i, j, k, kk, km1 - real :: gorho - real :: pref, DU, DV, DRHO, DZ, N2, S2 + real :: GoRho + real :: pref, DU, DV, DRHO, DZ, N2, S2, dummy real, dimension(2*(G%ke)) :: pres_1d, temp_1d, salt_1d, rho_1d real, dimension(G%ke+1) :: Ri_Grad !< Gradient Richardson number - + real, parameter :: epsln = 1.e-10 !< Threshold to identify + !! vanished layers ! some constants GoRho = GV%g_Earth / GV%Rho0 @@ -120,10 +121,30 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, & ! fill 3d arrays, if user asks for diagsnostics if (CS%id_N2 > 0) CS%N2(i,j,k) = N2 if (CS%id_S2 > 0) CS%S2(i,j,k) = S2 - if (CS%id_ri_grad > 0) CS%ri_grad(i,j,k) = Ri_Grad(k) enddo + Ri_grad(G%ke+1) = Ri_grad(G%ke) + + if (CS%smooth_ri) then + ! 1) fill Ri_grad in vanished layers with adjacent value + do k = 2, G%ke + if (h(i,j,k) .le. epsln) Ri_grad(k) = Ri_grad(k-1) + enddo + + Ri_grad(G%ke+1) = Ri_grad(G%ke) + + ! 2) vertically smooth Ri with 1-2-1 filter + dummy = 0.25 * Ri_grad(1) + Ri_grad(G%ke+1) = Ri_grad(G%ke) + do k = 1, G%ke + Ri_Grad(k) = dummy + 0.5 * Ri_Grad(k) + 0.25 * Ri_grad(k+1) + dummy = 0.25 * Ri_grad(k) + enddo + endif + + if (CS%id_ri_grad > 0) CS%ri_grad(i,j,:) = Ri_Grad(:) + ! Call to CVMix wrapper for computing interior mixing coefficients. call CVMix_coeffs_shear(Mdiff_out=kv(i,j,:), & Tdiff_out=kd(i,j,:), & @@ -209,7 +230,11 @@ logical function CVMix_shear_init(Time, G, GV, param_file, diag, CS) "Exponent of unitless factor of diffusivities,"// & " for KPP internal shear mixing scheme." & ,units="nondim", default=3.0) - call CVMix_init_shear(mix_scheme=CS%mix_scheme, & + call get_param(param_file, mdl, "SMOOTH_RI", CS%smooth_ri, & + "If true, vertically smooth the Richardson"// & + "number by applying a 1-2-1 filter once.", & + default = .false.) + call cvmix_init_shear(mix_scheme=CS%mix_scheme, & KPP_nu_zero=CS%Nu_Zero, & KPP_Ri_zero=CS%Ri_zero, & KPP_exp=CS%KPP_exp) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 2861218128..79234c7e11 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -14,6 +14,7 @@ module MOM_KPP use MOM_grid, only : ocean_grid_type, isPointInCell use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number +use MOM_domains, only : pass_var use CVMix_kpp, only : CVMix_init_kpp, CVMix_put_kpp, CVMix_get_kpp_real use CVMix_kpp, only : CVMix_coeffs_kpp @@ -29,6 +30,7 @@ module MOM_KPP #include "MOM_memory.h" public :: KPP_init +public :: KPP_compute_BLD public :: KPP_calculate public :: KPP_end public :: KPP_NonLocalTransport_temp @@ -86,6 +88,7 @@ module MOM_KPP character(len=30) :: MatchTechnique !< Method used in CVMix for setting diffusivity and NLT profile functions integer :: NLT_shape !< MOM6 over-ride of CVMix NLT shape function logical :: applyNonLocalTrans !< If True, apply non-local transport to heat and scalars + logical :: smoothBLD !< If True, apply a 1-1-4-1-1 Laplacian filter one time on HBLT. logical :: KPPzeroDiffusivity !< If True, will set diffusivity and viscosity from KPP to zero !! for testing purposes. logical :: KPPisAdditive !< If True, will add KPP diffusivity to initial diffusivity. @@ -135,6 +138,7 @@ module MOM_KPP ! Diagnostics arrays real, allocatable, dimension(:,:) :: OBLdepth !< Depth (positive) of OBL (m) + real, allocatable, dimension(:,:) :: kOBL !< Level (+fraction) of OBL extent real, allocatable, dimension(:,:) :: OBLdepthprev !< previous Depth (positive) of OBL (m) real, allocatable, dimension(:,:,:) :: dRho !< Bulk difference in density (kg/m3) real, allocatable, dimension(:,:,:) :: Uz2 !< Square of bulk difference in resolved velocity (m2/s2) @@ -210,6 +214,10 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) 'If False, calculates the non-local transport and tendencies but\n'//& 'purely for diagnostic purposes.', & default=.not. CS%passiveMode) + call get_param(paramFile, mdl, 'SMOOTH_BLD', CS%smoothBLD, & + 'If True, applies a 1-1-4-1-1 Laplacian filter one time on HBLT.\n'// & + 'computed via CVMix to reduce any horizontal two-grid-point noise.', & + default=.false.) call get_param(paramFile, mdl, 'RI_CRIT', CS%Ri_crit, & 'Critical bulk Richardson number used to define depth of the\n'// & 'surface Ocean Boundary Layer (OBL).', & @@ -486,9 +494,12 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) CS%id_EnhVt2 = register_diag_field('ocean_model', 'EnhVt2', diag%axesTL, Time, & 'Langmuir number enhancement to Vt2 as used by [CVMix] KPP','nondim') + allocate( CS%OBLdepth( SZI_(G), SZJ_(G) ) ) + CS%OBLdepth(:,:) = 0. + allocate( CS%kOBL( SZI_(G), SZJ_(G) ) ) + CS%kOBL(:,:) = 0. + allocate( CS%OBLdepthprev( SZI_(G), SZJ_(G) ) );CS%OBLdepthprev(:,:)=0.0 - if (CS%id_OBLdepth > 0) allocate( CS%OBLdepth( SZI_(G), SZJ_(G) ) ) - if (CS%id_OBLdepth > 0) CS%OBLdepth(:,:) = 0. if (CS%id_BulkDrho > 0) allocate( CS%dRho( SZI_(G), SZJ_(G), SZK_(G) ) ) if (CS%id_BulkDrho > 0) CS%dRho(:,:,:) = 0. if (CS%id_BulkUz2 > 0) allocate( CS%Uz2( SZI_(G), SZJ_(G), SZK_(G) ) ) @@ -527,8 +538,6 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) end function KPP_init - - !> KPP vertical diffusivity/viscosity and non-local tracer transport subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,& @@ -556,7 +565,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransHeat !< Temp non-local transport (m/s) real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local transport (m/s) - ! Local variables +! Local variables integer :: i, j, k, km1,kp1 ! Loop indices real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface (m) (negative in ocean) real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface (m) (negative in ocean) @@ -579,10 +588,9 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & real, dimension( 3*G%ke ) :: Temp_1D real, dimension( 3*G%ke ) :: Salt_1D - real :: kOBL, OBLdepth_0d, surfFricVel, surfBuoyFlux, Coriolis + real :: surfFricVel, surfBuoyFlux real :: GoRho, pRef, rho1, rhoK, rhoKm1, Uk, Vk, sigma, sigmaRatio - real :: zBottomMinusOffset ! Height of bottom plus a little bit (m) real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth. real :: hTot ! Running sum of thickness used in the surface layer average (m) real :: delH ! Thickness of a layer (m) @@ -643,8 +651,6 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & enddo ! things independent of position within the column - Coriolis = 0.25*( (G%CoriolisBu(i,j) + G%CoriolisBu(i-1,j-1)) & - +(G%CoriolisBu(i-1,j) + G%CoriolisBu(i,j-1)) ) surfFricVel = uStar(i,j) ! Bullk Richardson number computed for each cell in a column, @@ -683,6 +689,8 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & surfHu =0.0 surfHv =0.0 surfHuS =0.0 + surfHuS =0.0 + surfHvS =0.0 surfHvS =0.0 hTot =0.0 do ktmp = 1,ksfc @@ -851,109 +859,6 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & surfBuoyFlux = buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit ! h to Monin-Obukov (default is false, ie. not used) - call CVMix_kpp_compute_OBL_depth( & - BulkRi_1d, & ! (in) Bulk Richardson number - iFaceHeight, & ! (in) Height of interfaces (m) - OBLdepth_0d, & ! (out) OBL depth (m) - kOBL, & ! (out) level (+fraction) of OBL extent - zt_cntr=cellHeight, & ! (in) Height of cell centers (m) - surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) - surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) - Coriolis=Coriolis, & ! (in) Coriolis parameter (1/s) - CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters - - ! A hack to avoid KPP reaching the bottom. It was needed during development - ! because KPP was unable to handle vanishingly small layers near the bottom. - if (CS%deepOBLoffset>0.) then - zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) - OBLdepth_0d = min( OBLdepth_0d, -zBottomMinusOffset ) - endif - - ! apply some constraints on OBLdepth - if (CS%fixedOBLdepth) OBLdepth_0d = CS%fixedOBLdepth_value - OBLdepth_0d = max( OBLdepth_0d, -iFaceHeight(2) ) ! no shallower than top layer - OBLdepth_0d = min( OBLdepth_0d, -iFaceHeight(G%ke+1) ) ! no deeper than bottom - kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) - -!************************************************************************* -! smg: remove code below - -! Following "correction" step has been found to be unnecessary. -! Code should be removed after further testing. -! BGR: 03/15/2018-> Restructured code (Vt2 changed to compute from call in MOM_KPP now) -! I have not taken this restructuring into account here. -! Do we ever run with correctSurfLayerAvg? -! smg's suggested testing and removal is advised, in the meantime -! I have added warning if correctSurfLayerAvg is attempted. - ! if (CS%correctSurfLayerAvg) then - ! SLdepth_0d = CS%surf_layer_ext * OBLdepth_0d - ! hTot = h(i,j,1) - ! surfTemp = Temp(i,j,1) ; surfHtemp = surfTemp * hTot - ! surfSalt = Salt(i,j,1) ; surfHsalt = surfSalt * hTot - ! surfU = 0.5*(u(i,j,1)+u(i-1,j,1)) ; surfHu = surfU * hTot - ! surfV = 0.5*(v(i,j,1)+v(i,j-1,1)) ; surfHv = surfV * hTot - ! pRef = 0.0 - - ! do k = 2, G%ke - - ! ! Recalculate differences with surface layer - ! Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU - ! Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV - ! deltaU2(k) = Uk**2 + Vk**2 - ! pRef = pRef + GV%H_to_Pa * h(i,j,k) - ! call calculate_density(surfTemp, surfSalt, pRef, rho1, EOS) - ! call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS) - ! deltaRho(k) = rhoK - rho1 - - ! ! Surface layer averaging (needed for next k+1 iteration of this loop) - ! if (hTot < SLdepth_0d) then - ! delH = min( max(0., SLdepth_0d - hTot), h(i,j,k)*GV%H_to_m ) - ! hTot = hTot + delH - ! surfHtemp = surfHtemp + Temp(i,j,k) * delH ; surfTemp = surfHtemp / hTot - ! surfHsalt = surfHsalt + Salt(i,j,k) * delH ; surfSalt = surfHsalt / hTot - ! surfHu = surfHu + 0.5*(u(i,j,k)+u(i-1,j,k)) * delH ; surfU = surfHu / hTot - ! surfHv = surfHv + 0.5*(v(i,j,k)+v(i,j-1,k)) * delH ; surfV = surfHv / hTot - ! endif - - ! enddo - - ! BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( & - ! cellHeight(1:G%ke), & ! Depth of cell center (m) - ! GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) (1/s) - ! deltaU2, & ! Square of resolved velocity difference (m2/s2) - ! ws_cntr=Ws_1d, & ! Turbulent velocity scale profile (m/s) - ! N_iface=N_1d ) ! Buoyancy frequency (1/s) - - ! surfBuoyFlux = buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit - ! ! h to Monin-Obukov (default is false, ie. not used) - - ! call CVmix_kpp_compute_OBL_depth( & - ! BulkRi_1d, & ! (in) Bulk Richardson number - ! iFaceHeight, & ! (in) Height of interfaces (m) - ! OBLdepth_0d, & ! (out) OBL depth (m) - ! kOBL, & ! (out) level (+fraction) of OBL extent - ! zt_cntr=cellHeight, & ! (in) Height of cell centers (m) - ! surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) - ! surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) - ! Coriolis=Coriolis, & ! (in) Coriolis parameter (1/s) - ! CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters - - ! if (CS%deepOBLoffset>0.) then - ! zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) - ! OBLdepth_0d = min( OBLdepth_0d, -zBottomMinusOffset ) - ! kOBL = CVmix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) - ! endif - - ! ! apply some constraints on OBLdepth - ! if (CS%fixedOBLdepth) OBLdepth_0d = CS%fixedOBLdepth_value - ! OBLdepth_0d = max( OBLdepth_0d, -iFaceHeight(2) ) ! no shallower than top layer - ! OBLdepth_0d = min( OBLdepth_0d, -iFaceHeight(G%ke+1) ) ! no deep than bottom - ! kOBL = CVmix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) - - ! endif ! endif for "correction" step -! smg: remove code above -! ********************************************************************** - ! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports ! Unlike LMD94, we do not match to interior diffusivities. If using the original @@ -963,7 +868,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & if (CS%SW_METHOD == SW_METHOD_ALL_SW) then surfBuoyFlux = buoyFlux(i,j,1) elseif (CS%SW_METHOD == SW_METHOD_MXL_SW) then - surfBuoyFlux = buoyFlux(i,j,1) - buoyFlux(i,j,int(kOBL)+1) ! We know the actual buoyancy flux into the OBL + surfBuoyFlux = buoyFlux(i,j,1) - buoyFlux(i,j,int(CS%kOBL(i,j))+1) ! We know the actual buoyancy flux into the OBL elseif (CS%SW_METHOD == SW_METHOD_LV1_SW) then surfBuoyFlux = buoyFlux(i,j,1) - buoyFlux(i,j,2) endif @@ -986,8 +891,8 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & Kviscosity, & ! (in) Original viscosity (m2/s) Kdiffusivity(:,1), & ! (in) Original heat diffusivity (m2/s) Kdiffusivity(:,2), & ! (in) Original salt diffusivity (m2/s) - OBLdepth_0d, & ! (in) OBL depth (m) - kOBL, & ! (in) level (+fraction) of OBL extent + CS%OBLdepth(i,j), & ! (in) OBL depth (m) + CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent nonLocalTrans(:,1),& ! (out) Non-local heat transport (non-dimensional) nonLocalTrans(:,2),& ! (out) Non-local salt transport (non-dimensional) surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) @@ -1016,7 +921,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & Kdiffusivity(k,2) = Kdiffusivity(k,2) * LangEnhK Kviscosity(k) = Kviscosity(k) * LangEnhK elseif (CS%LT_K_SHAPE == LT_K_SCALED) then - sigma = min(1.0,-iFaceHeight(k)/OBLdepth_0d) + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) SigmaRatio = sigma * (1. - sigma)**2. / 0.148148037 if (CS%id_EnhK > 0) CS%EnhK(i,j,k) = (1.0 + (LangEnhK - 1.)*sigmaRatio) Kdiffusivity(k,1) = Kdiffusivity(k,1) * ( 1. + & @@ -1041,26 +946,26 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & if (surfBuoyFlux < 0.0) then if (CS%NLT_shape == NLT_SHAPE_CUBIC) then do k = 2, G%ke - sigma = min(1.0,-iFaceHeight(k)/OBLdepth_0d) + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) nonLocalTrans(k,1) = (1.0 - sigma)**2 * (1.0 + 2.0*sigma) !* nonLocalTrans(k,2) = nonLocalTrans(k,1) enddo elseif (CS%NLT_shape == NLT_SHAPE_PARABOLIC) then do k = 2, G%ke - sigma = min(1.0,-iFaceHeight(k)/OBLdepth_0d) + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) nonLocalTrans(k,1) = (1.0 - sigma)**2 !*CS%CS2 nonLocalTrans(k,2) = nonLocalTrans(k,1) enddo elseif (CS%NLT_shape == NLT_SHAPE_LINEAR) then do k = 2, G%ke - sigma = min(1.0,-iFaceHeight(k)/OBLdepth_0d) + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) nonLocalTrans(k,1) = (1.0 - sigma)!*CS%CS2 nonLocalTrans(k,2) = nonLocalTrans(k,1) enddo elseif (CS%NLT_shape == NLT_SHAPE_CUBIC_LMD) then ! Sanity check (should agree with CVMix result using simple matching) do k = 2, G%ke - sigma = min(1.0,-iFaceHeight(k)/OBLdepth_0d) + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) nonLocalTrans(k,1) = CS%CS2 * sigma*(1.0 -sigma)**2 nonLocalTrans(k,2) = nonLocalTrans(k,1) enddo @@ -1082,13 +987,13 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & ! recompute wscale for diagnostics, now that we in fact know boundary layer depth !BGR consider if LTEnhancement is wanted for diagnostics if (CS%id_Ws > 0) then - call CVmix_kpp_compute_turbulent_scales( & - -CellHeight/OBLdepth_0d, & ! (in) Normalized boundary layer coordinate - OBLdepth_0d, & ! (in) OBL depth (m) - surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) - surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) - w_s=Ws_1d, & ! (out) Turbulent velocity scale profile (m/s) - CVmix_kpp_params_user=CS%KPP_params) ! KPP parameters + call CVMix_kpp_compute_turbulent_scales( & + -CellHeight/CS%OBLdepth(i,j), & ! (in) Normalized boundary layer coordinate + CS%OBLdepth(i,j), & ! (in) OBL depth (m) + surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) + surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) + w_s=Ws_1d, & ! (out) Turbulent velocity scale profile (m/s) + CVMix_kpp_params_user=CS%KPP_params) ! KPP parameters CS%Ws(i,j,:) = Ws_1d(:) endif @@ -1106,15 +1011,13 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & ! Copy 1d data into 3d diagnostic arrays !/ grabbing obldepth_0d for next time step. - CS%OBLdepthprev(i,j)=OBLdepth_0d - !\ this can replace the other and be allocated independent of diagnostic output - if (CS%id_OBLdepth > 0) CS%OBLdepth(i,j) = OBLdepth_0d + CS%OBLdepthprev(i,j)=CS%OBLdepth(i,j) if (CS%id_BulkDrho > 0) CS%dRho(i,j,:) = deltaRho(:) if (CS%id_BulkUz2 > 0) CS%Uz2(i,j,:) = deltaU2(:) if (CS%id_BulkRi > 0) CS%BulkRi(i,j,:) = BulkRi_1d(:) if (CS%id_sigma > 0) then CS%sigma(i,j,:) = 0. - if (OBLdepth_0d>0.) CS%sigma(i,j,:) = -iFaceHeight/OBLdepth_0d + if (CS%OBLdepth(i,j)>0.) CS%sigma(i,j,:) = -iFaceHeight/CS%OBLdepth(i,j) endif if (CS%id_N > 0) CS%N(i,j,:) = N_1d(:) if (CS%id_N2 > 0) CS%N2(i,j,:) = N2_1d(:) @@ -1185,12 +1088,492 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & end subroutine KPP_calculate + +!> Compute OBL depth +subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Waves) + + ! Arguments + type(KPP_CS), pointer :: CS !< Control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid + type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses (units of H) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Temp !< potential/cons temp (deg C) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Salt !< Salinity (ppt) + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Velocity i-component (m/s) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Velocity j-component (m/s) + type(EOS_type), pointer :: EOS !< Equation of state + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity (m/s) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: buoyFlux !< Surface buoyancy flux (m2/s3) + + ! Local variables + integer :: i, j, k, km1 ! Loop indices + real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface (m) (negative in ocean) + real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface (m) (negative in ocean) + real, dimension( G%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces (1/s2) + real, dimension( G%ke+1 ) :: N_1d ! Brunt-Vaisala frequency at interfaces (1/s) (floored at 0) + real, dimension( G%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars (m/s) + !real, dimension( G%ke ) :: Wm_1d ! Profile of vertical velocity scale for momentum (m/s) + real, dimension( G%ke ) :: Vt2_1d ! Unresolved velocity for bulk Ri calculation/diagnostic (m2/s2) + real, dimension( G%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number + real, dimension( G%ke ) :: deltaU2 ! square of delta U (shear) in denominator of Bulk Ri (m2/s2) + real, dimension( G%ke+1, 2) :: nonLocalTrans ! Non-local transport for heat/salt at interfaces (non-dimensional) + real, dimension( G%ke ) :: surfBuoyFlux2 + real, dimension( G%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer + + ! for EOS calculation + real, dimension( 3*G%ke ) :: rho_1D + real, dimension( 3*G%ke ) :: pres_1D + real, dimension( 3*G%ke ) :: Temp_1D + real, dimension( 3*G%ke ) :: Salt_1D + + real :: surfFricVel, surfBuoyFlux, Coriolis + real :: GoRho, pRef, rho1, rhoK, Uk, Vk, sigma, sigmaRatio + + real :: zBottomMinusOffset ! Height of bottom plus a little bit (m) + real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth. + real :: hTot ! Running sum of thickness used in the surface layer average (m) + real :: delH ! Thickness of a layer (m) + real :: surfHtemp, surfTemp ! Integral and average of temp over the surface layer + real :: surfHsalt, surfSalt ! Integral and average of saln over the surface layer + real :: surfHu, surfU ! Integral and average of u over the surface layer + real :: surfHv, surfV ! Integral and average of v over the surface layer + real :: dh ! The local thickness used for calculating interface positions (m) + real :: hcorr ! A cumulative correction arising from inflation of vanished layers (m) + integer :: kk, ksfc, ktmp + + ! For Langmuir Calculations + real :: LangEnhW ! Langmuir enhancement for turbulent velocity scale + real, dimension(G%ke) :: LangEnhVt2 ! Langmuir enhancement for unresolved shear + real, dimension(G%ke) :: U_H, V_H + real :: MLD_GUESS, LA + real :: LangEnhK ! Langmuir enhancement for mixing coefficient + real :: surfHuS, surfHvS, surfUs, surfVs, wavedir, currentdir + real :: VarUp, VarDn, M, VarLo, VarAvg + real :: H10pct, H20pct,CMNFACT, USx20pct, USy20pct + integer :: B + real :: WST + + ! some constants + GoRho = GV%g_Earth / GV%Rho0 + nonLocalTrans(:,:) = 0.0 + +!$OMP parallel do default(private) shared(G,GV,CS,EOS,uStar,Temp,Salt,u,v,h,GoRho, & +!$OMP Waves,buoyFlux) & + + ! loop over horizontal points on processor + do j = G%jsc, G%jec + do i = G%isc, G%iec + + ! skip calling KPP for land points + if (G%mask2dT(i,j)==0.) cycle + + do k=1,G%ke + U_H(k) = 0.5 * (U(i,j,k)+U(i-1,j,k)) + V_H(k) = 0.5 * (V(i,j,k)+V(i,j-1,k)) + enddo + + ! things independent of position within the column + Coriolis = 0.25*( (G%CoriolisBu(i,j) + G%CoriolisBu(i-1,j-1)) & + +(G%CoriolisBu(i-1,j) + G%CoriolisBu(i,j-1)) ) + surfFricVel = uStar(i,j) + + ! Bullk Richardson number computed for each cell in a column, + ! assuming OBLdepth = grid cell depth. After Rib(k) is + ! known for the column, then CVMix interpolates to find + ! the actual OBLdepth. This approach avoids need to iterate + ! on the OBLdepth calculation. It follows that used in MOM5 + ! and POP. + iFaceHeight(1) = 0.0 ! BBL is all relative to the surface + pRef = 0. + hcorr = 0. + do k=1,G%ke + + ! cell center and cell bottom in meters (negative values in the ocean) + dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) + hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 + dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * dh + iFaceHeight(k+1) = iFaceHeight(k) - dh + + ! find ksfc for cell where "surface layer" sits + SLdepth_0d = CS%surf_layer_ext*max( max(-cellHeight(k),-iFaceHeight(2) ), CS%minOBLdepth ) + ksfc = k + do ktmp = 1,k + if (-1.0*iFaceHeight(ktmp+1) >= SLdepth_0d) then + ksfc = ktmp + exit + endif + enddo + + ! average temp, saln, u, v over surface layer + ! use C-grid average to get u,v on T-points. + surfHtemp=0.0 + surfHsalt=0.0 + surfHu =0.0 + surfHv =0.0 + surfHuS =0.0 + surfHvS =0.0 + hTot =0.0 + do ktmp = 1,ksfc + + ! SLdepth_0d can be between cell interfaces + delH = min( max(0.0, SLdepth_0d - hTot), h(i,j,ktmp)*GV%H_to_m ) + + ! surface layer thickness + hTot = hTot + delH + + ! surface averaged fields + surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH + surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH + surfHu = surfHu + 0.5*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delH + surfHv = surfHv + 0.5*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delH + if (CS%Stokes_Mixing) then + surfHus = surfHus + 0.5*(WAVES%US_x(i,j,ktmp)+WAVES%US_x(i-1,j,ktmp)) * delH + surfHvs = surfHvs + 0.5*(WAVES%US_y(i,j,ktmp)+WAVES%US_y(i,j-1,ktmp)) * delH + endif + + enddo + surfTemp = surfHtemp / hTot + surfSalt = surfHsalt / hTot + surfU = surfHu / hTot + surfV = surfHv / hTot + surfUs = surfHus / hTot + surfVs = surfHvs / hTot + + ! vertical shear between present layer and + ! surface layer averaged surfU,surfV. + ! C-grid average to get Uk and Vk on T-points. + Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU + Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV + + if (CS%Stokes_Mixing) then + ! If momentum is mixed down the Stokes drift gradient, then + ! the Stokes drift must be included in the bulk Richardson number + ! calculation. + Uk = Uk + (0.5*(Waves%Us_x(i,j,k)+Waves%US_x(i-1,j,k)) -surfUs ) + Vk = Vk + (0.5*(Waves%Us_y(i,j,k)+Waves%Us_y(i,j-1,k)) -surfVs ) + endif + + deltaU2(k) = Uk**2 + Vk**2 + + ! pressure, temp, and saln for EOS + ! kk+1 = surface fields + ! kk+2 = k fields + ! kk+3 = km1 fields + km1 = max(1, k-1) + kk = 3*(k-1) + pres_1D(kk+1) = pRef + pres_1D(kk+2) = pRef + pres_1D(kk+3) = pRef + Temp_1D(kk+1) = surfTemp + Temp_1D(kk+2) = Temp(i,j,k) + Temp_1D(kk+3) = Temp(i,j,km1) + Salt_1D(kk+1) = surfSalt + Salt_1D(kk+2) = Salt(i,j,k) + Salt_1D(kk+3) = Salt(i,j,km1) + + ! pRef is pressure at interface between k and km1. + ! iterate pRef for next pass through k-loop. + pRef = pRef + GV%H_to_Pa * h(i,j,k) + + ! this difference accounts for penetrating SW + surfBuoyFlux2(k) = buoyFlux(i,j,1) - buoyFlux(i,j,k+1) + + enddo ! k-loop finishes + + if (CS%LT_K_ENHANCEMENT .or. CS%LT_VT2_ENHANCEMENT) then + if (.not.(present(WAVES).and.associated(WAVES))) then + call MOM_error(FATAL,"Trying to use input WAVES information in KPP\n"//& + "without activating USEWAVES") + endif + !For now get Langmuir number based on prev. MLD (otherwise must compute 3d LA) + MLD_GUESS = max( 1., abs(CS%OBLdepthprev(i,j) ) ) + call get_Langmuir_Number( LA, G, GV, MLD_guess, surfFricVel, I, J, & + H=H(i,j,:), U_H=U_H, V_H=V_H, WAVES=WAVES) + WAVES%LangNum(i,j)=LA + endif + + + ! compute in-situ density + call calculate_density(Temp_1D, Salt_1D, pres_1D, rho_1D, 1, 3*G%ke, EOS) + + ! N2 (can be negative) and N (non-negative) on interfaces. + ! deltaRho is non-local rho difference used for bulk Richardson number. + ! N_1d is local N (with floor) used for unresolved shear calculation. + do k = 1, G%ke + km1 = max(1, k-1) + kk = 3*(k-1) + deltaRho(k) = rho_1D(kk+2) - rho_1D(kk+1) + N2_1d(k) = (GoRho * (rho_1D(kk+2) - rho_1D(kk+3)) ) / & + ((0.5*(h(i,j,km1) + h(i,j,k))+GV%H_subroundoff)*GV%H_to_m) + N_1d(k) = sqrt( max( N2_1d(k), 0.) ) + enddo + N2_1d(G%ke+1 ) = 0.0 + N_1d(G%ke+1 ) = 0.0 + + ! turbulent velocity scales w_s and w_m computed at the cell centers. + ! Note that if sigma > CS%surf_layer_ext, then CVMix_kpp_compute_turbulent_scales + ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass + ! sigma=CS%surf_layer_ext for this calculation. + call CVMix_kpp_compute_turbulent_scales( & + CS%surf_layer_ext, & ! (in) Normalized surface layer depth; sigma = CS%surf_layer_ext + -cellHeight, & ! (in) Assume here that OBL depth (m) = -cellHeight(k) + surfBuoyFlux2, & ! (in) Buoyancy flux at surface (m2/s3) + surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) + w_s=Ws_1d, & ! (out) Turbulent velocity scale profile (m/s) + CVMix_kpp_params_user=CS%KPP_params ) + + !Compute CVMix VT2 + Vt2_1d(:) = CVmix_kpp_compute_unresolved_shear( & + zt_cntr=cellHeight(1:G%ke), & ! Depth of cell center (m) + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers (m/s) + N_iface=N_1d, & ! Buoyancy frequency at interface (1/s) + CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters + + !Modify CVMix VT2 + IF (CS%LT_VT2_ENHANCEMENT) then + IF (CS%LT_VT2_METHOD==LT_VT2_MODE_CONSTANT) then + do k=1,G%ke + LangEnhVT2(k) = CS%KPP_VT2_ENH_FAC + enddo + elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_VR12) then + do k=1,G%ke + LangEnhVT2(k) = min(10.,sqrt(1.+(1.5*WAVES%LangNum(i,j))**(-2) + & + (5.4*WAVES%LangNum(i,j))**(-4))) + enddo + elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_RW16) then + do k=1,G%ke + LangEnhVT2(k) = min(2.25, 1. + 1./WAVES%LangNum(i,j)) + enddo + elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_LF17) then + CS%CS=cvmix_get_kpp_real('c_s',CS%KPP_params) + do k=1,G%ke + WST = (max(0.,-buoyflux(i,j,1))*(-cellHeight(k)))**(1./3.) + LangEnhVT2(k) = sqrt((0.15*WST**3. + 0.17*surfFricVel**3.* & + (1.+0.49*WAVES%LangNum(i,j)**(-2.))) / & + (0.2*ws_1d(k)**3/(CS%cs*CS%surf_layer_ext*CS%vonKarman**4.))) + enddo + else + !This shouldn't be reached. + !call MOM_error(WARNING,"Unexpected behavior in MOM_KPP, see error in Vt2") + LangEnhVT2(:) = 1.0 + endif + else + LangEnhVT2(:) = 1.0 + endif + + do k=1,G%ke + Vt2_1d(k)=Vt2_1d(k)*LangEnhVT2(k) + if (CS%id_EnhVt2 > 0) CS%EnhVt2(i,j,k)=LangEnhVT2(k) + enddo + + ! Calculate Bulk Richardson number from eq (21) of LMD94 + BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( & + zt_cntr = cellHeight(1:G%ke), & ! Depth of cell center (m) + delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) (1/s) + delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference (m2/s2) + Vt_sqr_cntr=Vt2_1d, & + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile (m/s) + N_iface=N_1d) ! Buoyancy frequency (1/s) + + + surfBuoyFlux = buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit + ! h to Monin-Obukov (default is false, ie. not used) + + call CVMix_kpp_compute_OBL_depth( & + BulkRi_1d, & ! (in) Bulk Richardson number + iFaceHeight, & ! (in) Height of interfaces (m) + CS%OBLdepth(i,j), & ! (out) OBL depth (m) + CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent + zt_cntr=cellHeight, & ! (in) Height of cell centers (m) + surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) + surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) + Coriolis=Coriolis, & ! (in) Coriolis parameter (1/s) + CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters + + ! A hack to avoid KPP reaching the bottom. It was needed during development + ! because KPP was unable to handle vanishingly small layers near the bottom. + if (CS%deepOBLoffset>0.) then + zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) + CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset ) + endif + + ! apply some constraints on OBLdepth + if(CS%fixedOBLdepth) CS%OBLdepth(i,j) = CS%fixedOBLdepth_value + CS%OBLdepth(i,j) = max( CS%OBLdepth(i,j), -iFaceHeight(2) ) ! no shallower than top layer + CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(G%ke+1) ) ! no deeper than bottom + CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + +!************************************************************************* +! smg: remove code below + +! Following "correction" step has been found to be unnecessary. +! Code should be removed after further testing. +! BGR: 03/15/2018-> Restructured code (Vt2 changed to compute from call in MOM_KPP now) +! I have not taken this restructuring into account here. +! Do we ever run with correctSurfLayerAvg? +! smg's suggested testing and removal is advised, in the meantime +! I have added warning if correctSurfLayerAvg is attempted. + ! if (CS%correctSurfLayerAvg) then + + ! SLdepth_0d = CS%surf_layer_ext * CS%OBLdepth(i,j) + ! hTot = h(i,j,1) + ! surfTemp = Temp(i,j,1) ; surfHtemp = surfTemp * hTot + ! surfSalt = Salt(i,j,1) ; surfHsalt = surfSalt * hTot + ! surfU = 0.5*(u(i,j,1)+u(i-1,j,1)) ; surfHu = surfU * hTot + ! surfV = 0.5*(v(i,j,1)+v(i,j-1,1)) ; surfHv = surfV * hTot + ! pRef = 0.0 + + ! do k = 2, G%ke + + ! ! Recalculate differences with surface layer + ! Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU + ! Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV + ! deltaU2(k) = Uk**2 + Vk**2 + ! pRef = pRef + GV%H_to_Pa * h(i,j,k) + ! call calculate_density(surfTemp, surfSalt, pRef, rho1, EOS) + ! call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS) + ! deltaRho(k) = rhoK - rho1 + + ! ! Surface layer averaging (needed for next k+1 iteration of this loop) + ! if (hTot < SLdepth_0d) then + ! delH = min( max(0., SLdepth_0d - hTot), h(i,j,k)*GV%H_to_m ) + ! hTot = hTot + delH + ! surfHtemp = surfHtemp + Temp(i,j,k) * delH ; surfTemp = surfHtemp / hTot + ! surfHsalt = surfHsalt + Salt(i,j,k) * delH ; surfSalt = surfHsalt / hTot + ! surfHu = surfHu + 0.5*(u(i,j,k)+u(i-1,j,k)) * delH ; surfU = surfHu / hTot + ! surfHv = surfHv + 0.5*(v(i,j,k)+v(i,j-1,k)) * delH ; surfV = surfHv / hTot + ! endif + + ! enddo + + ! BulkRi_1d = CVMix_kpp_compute_bulk_Richardson( & + ! cellHeight(1:G%ke), & ! Depth of cell center (m) + ! GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) (1/s) + ! deltaU2, & ! Square of resolved velocity difference (m2/s2) + ! ws_cntr=Ws_1d, & ! Turbulent velocity scale profile (m/s) + ! N_iface=N_1d ) ! Buoyancy frequency (1/s) + + ! surfBuoyFlux = buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit + ! ! h to Monin-Obukov (default is false, ie. not used) + + ! call CVMix_kpp_compute_OBL_depth( & + ! BulkRi_1d, & ! (in) Bulk Richardson number + ! iFaceHeight, & ! (in) Height of interfaces (m) + ! CS%OBLdepth(i,j), & ! (out) OBL depth (m) + ! CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent + ! zt_cntr=cellHeight, & ! (in) Height of cell centers (m) + ! surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) + ! surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) + ! Coriolis=Coriolis, & ! (in) Coriolis parameter (1/s) + ! CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters + + ! if (CS%deepOBLoffset>0.) then + ! zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) + ! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset ) + ! CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + ! endif + + ! ! apply some constraints on OBLdepth + ! if(CS%fixedOBLdepth) CS%OBLdepth(i,j) = CS%fixedOBLdepth_value + ! CS%OBLdepth(i,j) = max( CS%OBLdepth(i,j), -iFaceHeight(2) ) ! no shallower than top layer + ! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(G%ke+1) ) ! no deep than bottom + ! CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + + ! endif ! endif for "correction" step + +! smg: remove code above +! ********************************************************************** + + enddo + enddo + + if (CS%smoothBLD) call KPP_smooth_BLD(CS,G,GV,h) + +end subroutine KPP_compute_BLD + + +!> Apply a 1-1-4-1-1 Laplacian filter one time on BLD to reduce any horizontal two-grid-point noise +subroutine KPP_smooth_BLD(CS,G,GV,h) + ! Arguments + type(KPP_CS), pointer :: CS !< Control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses (units of H) + + ! local + real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_original ! Original OBL depths computed by CVMix + real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface (m) (negative in ocean) + real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface (m) (negative in ocean) + real :: wc, ww, we, wn, ws ! averaging weights for smoothing + real :: dh ! The local thickness used for calculating interface positions (m) + real :: hcorr ! A cumulative correction arising from inflation of vanished layers (m) + integer :: i, j, k + + ! Update halos + call pass_var(CS%OBLdepth, G%Domain) + + OBLdepth_original = CS%OBLdepth + + ! apply smoothing on OBL depth + do j = G%jsc, G%jec + do i = G%isc, G%iec + + ! skip land points + if (G%mask2dT(i,j)==0.) cycle + + ! compute weights + ww = 0.125 * G%mask2dT(i-1,j) + we = 0.125 * G%mask2dT(i+1,j) + ws = 0.125 * G%mask2dT(i,j-1) + wn = 0.125 * G%mask2dT(i,j+1) + wc = 1.0 - (ww+we+wn+ws) + + CS%OBLdepth(i,j) = wc * OBLdepth_original(i,j) & + + ww * OBLdepth_original(i-1,j) & + + we * OBLdepth_original(i+1,j) & + + ws * OBLdepth_original(i,j-1) & + + wn * OBLdepth_original(i,j+1) + enddo + enddo + + ! Update kOBL for smoothed OBL depths + do j = G%jsc, G%jec + do i = G%isc, G%iec + + ! skip land points + if (G%mask2dT(i,j)==0.) cycle + + iFaceHeight(1) = 0.0 ! BBL is all relative to the surface + hcorr = 0. + do k=1,G%ke + + ! cell center and cell bottom in meters (negative values in the ocean) + dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) + hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 + dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * dh + iFaceHeight(k+1) = iFaceHeight(k) - dh + enddo + + CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + + enddo + enddo + +end subroutine KPP_smooth_BLD + + + !> Copies KPP surface boundary layer depth into BLD subroutine KPP_get_BLD(CS, BLD, G) type(KPP_CS), pointer :: CS !< Control structure for !! this module type(ocean_grid_type), intent(in) :: G !< Grid structure - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: BLD!< bnd. layer depth (m) + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD!< bnd. layer depth (m) ! Local variables integer :: i,j do j = G%jsc, G%jec ; do i = G%isc, G%iec diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 61c212db8b..bb1e0b11c1 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -408,7 +408,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, kd_lay, kv, j, G, GV, CS) CS%kd_bkgnd(i,j,1) = 0.0; CS%kv_bkgnd(i,j,1) = 0.0 CS%kd_bkgnd(i,j,nz+1) = 0.0; CS%kv_bkgnd(i,j,nz+1) = 0.0 do k=2,nz - CS%kd_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) + 0.5*(kd_lay(i,j,K-1) + kd_lay(i,j,K)) + CS%kd_bkgnd(i,j,k) = 0.5*(kd_lay(i,j,K-1) + kd_lay(i,j,K)) CS%kv_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) * CS%prandtl_bkgnd enddo enddo diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 6eb3b854f4..528dc33135 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -2,53 +2,6 @@ module MOM_diabatic_aux ! This file is part of MOM6. See LICENSE.md for the license. -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, April 1994 - July 2000 * -!* Alistair Adcroft, and Stephen Griffies * -!* * -!* This program contains the subroutine that, along with the * -!* subroutines that it calls, implements diapycnal mass and momentum * -!* fluxes and a bulk mixed layer. The diapycnal diffusion can be * -!* used without the bulk mixed layer. * -!* * -!* diabatic first determines the (diffusive) diapycnal mass fluxes * -!* based on the convergence of the buoyancy fluxes within each layer. * -!* The dual-stream entrainment scheme of MacDougall and Dewar (JPO, * -!* 1997) is used for combined diapycnal advection and diffusion, * -!* calculated implicitly and potentially with the Richardson number * -!* dependent mixing, as described by Hallberg (MWR, 2000). Diapycnal * -!* advection is fundamentally the residual of diapycnal diffusion, * -!* so the fully implicit upwind differencing scheme that is used is * -!* entirely appropriate. The downward buoyancy flux in each layer * -!* is determined from an implicit calculation based on the previously * -!* calculated flux of the layer above and an estimated flux in the * -!* layer below. This flux is subject to the following conditions: * -!* (1) the flux in the top and bottom layers are set by the boundary * -!* conditions, and (2) no layer may be driven below an Angstrom thick-* -!* ness. If there is a bulk mixed layer, the buffer layer is treat- * -!* ed as a fixed density layer with vanishingly small diffusivity. * -!* * -!* diabatic takes 5 arguments: the two velocities (u and v), the * -!* thicknesses (h), a structure containing the forcing fields, and * -!* the length of time over which to act (dt). The velocities and * -!* thickness are taken as inputs and modified within the subroutine. * -!* There is no limit on the time step. * -!* * -!* A small fragment of the grid is shown below: * -!* * -!* j+1 x ^ x ^ x At x: q * -!* j+1 > o > o > At ^: v * -!* j x ^ x ^ x At >: u * -!* j > o > o > At o: h, T, S, buoy, ustar, ea, eb, etc. * -!* j-1 x ^ x ^ x * -!* i-1 i i+1 At x & ^: * -!* i i+1 At > & o: * -!* * -!* The boundaries always run through q grid points (x). * -!* * -!********+*********+*********+*********+*********+*********+*********+** - use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr @@ -239,26 +192,19 @@ subroutine make_frazil(h, tv, G, GV, CS, p_surf) end subroutine make_frazil +!> Applies double diffusion to T & S, assuming no diapycal mass +!! fluxes, using a simple triadiagonal solver. subroutine differential_diffuse_T_S(h, tv, visc, dt, G, GV) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) - type(thermo_var_ptrs), intent(inout) :: tv - type(vertvisc_type), intent(in) :: visc - real, intent(in) :: dt - -! This subroutine applies double diffusion to T & S, assuming no diapycal mass -! fluxes, using a simple triadiagonal solver. - -! Arguments: h - Layer thickness, in m or kg m-2. -! (in) tv - A structure containing pointers to any available -! thermodynamic fields. Absent fields have NULL ptrs. -! (in) visc - A structure containing vertical viscosities, bottom boundary -! layer properies, and related fields. -! (in) dt - Time increment, in s. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. + type(thermo_var_ptrs), intent(inout) :: tv !< pointers to any available modynamic fields. + !! Absent fields have NULL ptrs. + type(vertvisc_type), intent(in) :: visc !< structure containing vertical viscosities, + !! layer properies, and related fields. + real, intent(in) :: dt !< Time increment, in s. + ! local variables real, dimension(SZI_(G)) :: & b1_T, b1_S, & ! Variables used by the tridiagonal solvers of T & S, in H. d1_T, d1_S ! Variables used by the tridiagonal solvers, nondim. @@ -345,30 +291,25 @@ subroutine differential_diffuse_T_S(h, tv, visc, dt, G, GV) S(i,j,k) = S(i,j,k) + c1_S(i,k+1)*S(i,j,k+1) enddo ; enddo enddo - end subroutine differential_diffuse_T_S +!> Keep salinity from falling below a small but positive threshold +!! This occurs when the ice model attempts to extract more salt then +!! is actually available to it from the ocean. subroutine adjust_salt(h, tv, G, GV, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) - type(thermo_var_ptrs), intent(inout) :: tv - type(diabatic_aux_CS), intent(in) :: CS - -! Keep salinity from falling below a small but positive threshold -! This occurs when the ice model attempts to extract more salt then -! is actually available to it from the ocean. - -! Arguments: h - Layer thickness, in m. -! (in/out) tv - A structure containing pointers to any available -! thermodynamic fields. Absent fields have NULL ptrs. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) CS - The control structure returned by a previous call to -! diabatic_driver_init. - real :: salt_add_col(SZI_(G),SZJ_(G)) ! The accumulated salt requirement - real :: S_min ! The minimum salinity - real :: mc ! A layer's mass kg m-2 . + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m + !! or kg m-2) + type(thermo_var_ptrs), intent(inout) :: tv !< structure containing pointers to any + !! available thermodynamic fields. + type(diabatic_aux_CS), intent(in) :: CS !< control structure returned by + !! a previous call to diabatic_driver_init. + + ! local variables + real :: salt_add_col(SZI_(G),SZJ_(G)) !< The accumulated salt requirement + real :: S_min !< The minimum salinity + real :: mc !< A layer's mass kg m-2 . integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -410,33 +351,29 @@ subroutine adjust_salt(h, tv, G, GV, CS) end subroutine adjust_salt +!> Insert salt from brine rejection into the first layer below +!! the mixed layer which both contains mass and in which the +!! change in layer density remains stable after the addition +!! of salt via brine rejection. subroutine insert_brine(h, tv, G, GV, fluxes, nkmb, CS, dt, id_brine_lay) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) - type(thermo_var_ptrs), intent(inout) :: tv - type(forcing), intent(in) :: fluxes - integer, intent(in) :: nkmb - type(diabatic_aux_CS), intent(in) :: CS - real, intent(in) :: dt + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m + !! or kg m-2) + type(thermo_var_ptrs), intent(inout) :: tv !< structure containing pointers to + !! any available hermodynamic fields. + type(forcing), intent(in) :: fluxes !< tructure containing pointers + !! any possible forcing fields + integer, intent(in) :: nkmb !< number of layers in the mixed and + !! buffer layers + type(diabatic_aux_CS), intent(in) :: CS !< control structure returned by a + !! previous call to diabatic_driver_init. + real, intent(in) :: dt !< time step between calls to this + !! function (s) ?? integer, intent(in) :: id_brine_lay -! Insert salt from brine rejection into the first layer below -! the mixed layer which both contains mass and in which the -! change in layer density remains stable after the addition -! of salt via brine rejection. - -! Arguments: h - Layer thickness, in m. -! (in/out) tv - A structure containing pointers to any available -! thermodynamic fields. Absent fields have NULL ptrs. -! (in) fluxes = A structure containing pointers to any possible -! forcing fields; unused fields have NULL ptrs. -! (in) nkmb - The number of layers in the mixed and buffer layers. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) CS - The control structure returned by a previous call to -! diabatic_driver_init. + ! local variables real :: salt(SZI_(G)) ! The amount of salt rejected from ! sea ice. [grams] real :: dzbr(SZI_(G)) ! cumulative depth over which brine is distributed @@ -539,10 +476,9 @@ subroutine insert_brine(h, tv, G, GV, fluxes, nkmb, CS, dt, id_brine_lay) end subroutine insert_brine +!> Simple tri-diagnonal solver for T and S. +!! "Simple" means it only uses arrays hold, ea and eb. subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) -! Simple tri-diagnonal solver for T and S -! "Simple" means it only uses arrays hold, ea and eb - ! Arguments type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure integer, intent(in) :: is, ie, js, je @@ -579,37 +515,22 @@ subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) enddo end subroutine triDiagTS - +!> Calculates u_h and v_h (velocities at thickness points), +!! optionally using the entrainments (in m) passed in as arguments. subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, ea, eb) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< The zonal velocity, in m s-1 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< The meridional velocity, in m s-1 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: u_h - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: v_h - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: ea - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: eb -! This subroutine calculates u_h and v_h (velocities at thickness -! points), optionally using the entrainments (in m) passed in as arguments. - -! Arguments: u - Zonal velocity, in m s-1. -! (in) v - Meridional velocity, in m s-1. -! (in) h - Layer thickness, in m or kg m-2. -! (out) u_h - The zonal velocity at thickness points after -! entrainment, in m s-1. -! (out) v_h - The meridional velocity at thickness points after -! entrainment, in m s-1. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in, opt) ea - The amount of fluid entrained from the layer above within -! this time step, in units of m or kg m-2. Omitting ea is the -! same as setting it to 0. -! (in, opt) eb - The amount of fluid entrained from the layer below within -! this time step, in units of m or kg m-2. Omitting eb is the -! same as setting it to 0. ea and eb must either be both -! present or both absent. - + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: u_h, v_h !< zonal and meridional velocity at thickness + !! points entrainment, in m s-1. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in), optional :: ea, eb !< The amount of fluid entrained + !! from the layer above within this time step + !! , in units of m or kg m-2. Omitting ea is the + !! same as setting it to 0. + + ! local variables real :: b_denom_1 ! The first term in the denominator of b1 in m or kg m-2. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected, in m or kg m-2. @@ -1318,26 +1239,20 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, & end subroutine applyBoundaryFluxesInOut +!> Initializes this module. subroutine diabatic_aux_init(Time, G, GV, param_file, diag, CS, useALEalgorithm, use_ePBL) type(time_type), intent(in) :: Time type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - type(diag_ctrl), target, intent(inout) :: diag - type(diabatic_aux_CS), pointer :: CS - logical, intent(in) :: useALEalgorithm - logical, intent(in) :: use_ePBL - -! Arguments: -! (in) Time = current model time -! (in) G = ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) param_file = structure indicating the open file to parse for parameter values -! (in) diag = structure used to regulate diagnostic output -! (in/out) CS = pointer set to point to the control structure for this module -! (in) use_ePBL = If true, use the implicit energetics planetary boundary -! layer scheme to determine the diffusivity in the -! surface boundary layer. + type(diag_ctrl), target, intent(inout) :: diag !< structure used to regulate diagnostic output + type(diabatic_aux_CS), pointer :: CS !< pointer set to point to the ontrol structure for + !! this module + logical, intent(in) :: useALEalgorithm !< If True, uses ALE. + logical, intent(in) :: use_ePBL !< If true, use the implicit energetics + !! planetary boundary layer scheme to determine the + !! diffusivity in the surface boundary layer. + ! local variables type(vardesc) :: vd ! This "include" declares and sets the variable "version". @@ -1460,4 +1375,48 @@ subroutine diabatic_aux_end(CS) end subroutine diabatic_aux_end +!> \namespace MOM_diabatic_aux +!! +!! This module contains the subroutines that, along with the * +!! subroutines that it calls, implements diapycnal mass and momentum * +!! fluxes and a bulk mixed layer. The diapycnal diffusion can be * +!! used without the bulk mixed layer. * +!! * +!! diabatic first determines the (diffusive) diapycnal mass fluxes * +!! based on the convergence of the buoyancy fluxes within each layer. * +!! The dual-stream entrainment scheme of MacDougall and Dewar (JPO, * +!! 1997) is used for combined diapycnal advection and diffusion, * +!! calculated implicitly and potentially with the Richardson number * +!! dependent mixing, as described by Hallberg (MWR, 2000). Diapycnal * +!! advection is fundamentally the residual of diapycnal diffusion, * +!! so the fully implicit upwind differencing scheme that is used is * +!! entirely appropriate. The downward buoyancy flux in each layer * +!! is determined from an implicit calculation based on the previously * +!! calculated flux of the layer above and an estimated flux in the * +!! layer below. This flux is subject to the following conditions: * +!! (1) the flux in the top and bottom layers are set by the boundary * +!! conditions, and (2) no layer may be driven below an Angstrom thick-* +!! ness. If there is a bulk mixed layer, the buffer layer is treat- * +!! ed as a fixed density layer with vanishingly small diffusivity. * +!! * +!! diabatic takes 5 arguments: the two velocities (u and v), the * +!! thicknesses (h), a structure containing the forcing fields, and * +!! the length of time over which to act (dt). The velocities and * +!! thickness are taken as inputs and modified within the subroutine. * +!! There is no limit on the time step. * +!! * +!! A small fragment of the grid is shown below: * +!! * +!! j+1 x ^ x ^ x At x: q * +!! j+1 > o > o > At ^: v * +!! j x ^ x ^ x At >: u * +!! j > o > o > At o: h, T, S, buoy, ustar, ea, eb, etc. * +!! j-1 x ^ x ^ x * +!! i-1 i i+1 At x & ^: * +!! i i+1 At > & o: * +!! * +!! The boundaries always run through q grid points (x). * +!! * +!!********+*********+*********+*********+*********+*********+*********+** + end module MOM_diabatic_aux diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 0db5fbd5b3..698243a7f6 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -10,6 +10,7 @@ module MOM_diabatic_driver use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_CVMix_shear, only : CVMix_shear_is_used +use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used use MOM_diabatic_aux, only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_CS use MOM_diabatic_aux, only : make_frazil, adjust_salt, insert_brine, differential_diffuse_T_S, triDiagTS use MOM_diabatic_aux, only : find_uv_at_h, diagnoseMLDbyDensityDifference, applyBoundaryFluxesInOut @@ -49,7 +50,8 @@ module MOM_diabatic_driver use MOM_internal_tides, only : propagate_int_tide use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS use MOM_kappa_shear, only : kappa_shear_is_used -use MOM_KPP, only : KPP_CS, KPP_init, KPP_calculate, KPP_end, KPP_get_BLD +use MOM_KPP, only : KPP_CS, KPP_init, KPP_compute_BLD, KPP_calculate +use MOM_KPP, only : KPP_end, KPP_get_BLD use MOM_KPP, only : KPP_NonLocalTransport_temp, KPP_NonLocalTransport_saln use MOM_opacity, only : opacity_init, set_opacity, opacity_end, opacity_CS use MOM_regularize_layers, only : regularize_layers, regularize_layers_init, regularize_layers_CS @@ -83,7 +85,10 @@ module MOM_diabatic_driver public adiabatic_driver_init !> Control structure for this module -type, public:: diabatic_CS ; private +! GMM, I've made the following type public so it work with the legacy version of +! diabatic. This type should be made private once the legacy code is deleted. +!type, public:: diabatic_CS; private +type, public:: diabatic_CS logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with !! nkml sublayers (and additional buffer layers). logical :: use_energetic_PBL !< If true, use the implicit energetics planetary @@ -93,6 +98,7 @@ module MOM_diabatic_driver !! shear-driven diapycnal diffusivity. logical :: use_CVMix_shear !< If true, use the CVMix module to find the !! shear-driven diapycnal diffusivity. + logical :: use_CVMix_ddiff !< If true, use the CVMix double diffusion module. logical :: use_tidal_mixing !< If true, activate tidal mixing diffusivity. logical :: use_CVMix_conv !< If true, use the CVMix module to get enhanced !! mixing due to convection. @@ -245,7 +251,7 @@ module MOM_diabatic_driver integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap -integer :: id_clock_kpp +integer :: id_clock_kpp, id_clock_CVMix_ddiff contains @@ -381,7 +387,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect*h_neglect Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 - if (nz == 1) return showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("diabatic(), MOM_diabatic_driver.F90") @@ -483,13 +488,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif if (CS%ML_mix_first > 0.0) then -! This subroutine -! (1) Cools the mixed layer. -! (2) Performs convective adjustment by mixed layer entrainment. -! (3) Heats the mixed layer and causes it to detrain to -! Monin-Obukhov depth or minimum mixed layer depth. -! (4) Uses any remaining TKE to drive mixed layer entrainment. -! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) + ! This subroutine: + ! (1) Cools the mixed layer. + ! (2) Performs convective adjustment by mixed layer entrainment. + ! (3) Heats the mixed layer and causes it to detrain to + ! Monin-Obukhov depth or minimum mixed layer depth. + ! (4) Uses any remaining TKE to drive mixed layer entrainment. + ! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) call find_uv_at_h(u, v, h, u_h, v_h, G, GV) call cpu_clock_begin(id_clock_mixedlayer) @@ -524,11 +529,12 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)") if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G) endif - endif + endif ! end CS%bulkmixedlayer if (CS%debug) then call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, haloshift=0) endif + if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then call find_uv_at_h(u, v, h_orig, u_h, v_h, G, GV, eaml, ebml) @@ -585,7 +591,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & CS%int_tide_input%tideamp, CS%int_tide_input%Nb, dt, G, GV, CS%int_tide_CSp) endif if (showCallTree) call callTree_waypoint("done with propagate_int_tide (diabatic)") - endif + endif ! end CS%use_int_tides call cpu_clock_begin(id_clock_set_diffusivity) ! Sets: Kd, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S @@ -638,6 +644,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif !$OMP end parallel + call KPP_compute_BLD(CS%KPP_CSp, G, GV, h, tv%T, tv%S, u, v, tv%eqn_of_state, & + fluxes%ustar, CS%KPP_buoy_flux) + call KPP_calculate(CS%KPP_CSp, G, GV, h, tv%T, tv%S, u, v, tv%eqn_of_state, & fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, & CS%KPP_NLTscalar, Waves=Waves) @@ -723,10 +732,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! a diffusivity and happen before KPP. But generally in MOM, we do not match ! KPP boundary layer to interior, so this diffusivity can be computed when convenient. if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T)) then - call cpu_clock_begin(id_clock_differential_diff) + call cpu_clock_begin(id_clock_CVMix_ddiff) call differential_diffuse_T_S(h, tv, visc, dt, G, GV) - call cpu_clock_end(id_clock_differential_diff) + call cpu_clock_end(id_clock_CVMix_ddiff) if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)") if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G) @@ -739,7 +748,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo endif - endif @@ -1374,6 +1382,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! visc%Kv_shear is not in the group pass because it has larger vertical extent. if (associated(visc%Kv_shear)) & call pass_var(visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) + if (associated(visc%Kv_slow)) & + call pass_var(visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) + call cpu_clock_end(id_clock_pass) if (.not. CS%useALEalgorithm) then @@ -1878,7 +1889,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, real :: Kd integer :: num_mode - logical :: use_temperature, differentialDiffusion + logical :: use_temperature type(vardesc) :: vd ! This "include" declares and sets the variable "version". @@ -1913,7 +1924,6 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, ! Set default, read and log parameters call log_version(param_file, mod, version, & "The following parameters are used for diabatic processes.") - call get_param(param_file, mod, "SPONGE", CS%use_sponge, & "If true, sponges may be applied anywhere in the domain. \n"//& "The exact location and properties of those sponges are \n"//& @@ -1930,11 +1940,10 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, the diffusivity from ePBL is added to all\n"//& "other diffusivities. Otherwise, the larger of kappa-\n"//& "shear and ePBL diffusivities are used.", default=.true.) - call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialDiffusion, & - "If true, apply parameterization of double-diffusion.", & - default=.false. ) + CS%use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) CS%use_kappa_shear = kappa_shear_is_used(param_file) CS%use_CVMix_shear = CVMix_shear_is_used(param_file) + if (CS%bulkmixedlayer) then call get_param(param_file, mod, "ML_MIX_FIRST", CS%ML_mix_first, & "The fraction of the mixed layer mixing that is applied \n"//& @@ -2393,8 +2402,8 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=CLOCK_MODULE) id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=CLOCK_ROUTINE) id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=CLOCK_ROUTINE) - id_clock_differential_diff = -1 ; if (differentialDiffusion) & - id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=CLOCK_ROUTINE) + id_clock_CVMix_ddiff = -1 ; if (CS%use_CVMix_ddiff) & + id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion)', grain=CLOCK_ROUTINE) ! initialize the auxiliary diabatic driver module call diabatic_aux_init(Time, G, GV, param_file, diag, CS%diabatic_aux_CSp, & diff --git a/src/parameterizations/vertical/MOM_legacy_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_legacy_diabatic_driver.F90 new file mode 100644 index 0000000000..739c74c80c --- /dev/null +++ b/src/parameterizations/vertical/MOM_legacy_diabatic_driver.F90 @@ -0,0 +1,1660 @@ +!> This routine drives the diabatic/dianeutral physics for MOM. +!! This is a legacy module that will be deleted in the near future. +module MOM_legacy_diabatic_driver + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_bulk_mixed_layer, only : bulkmixedlayer, bulkmixedlayer_init, bulkmixedlayer_CS +use MOM_debugging, only : hchksum +use MOM_checksum_packages, only : MOM_state_chksum, MOM_state_stats +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE +use MOM_CVMix_shear, only : CVMix_shear_is_used +use MOM_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 +use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr +use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_remap_grids +use MOM_diag_mediator, only : diag_ctrl, query_averaging_enabled, enable_averaging, disable_averaging +use MOM_diag_mediator, only : diag_grid_storage, diag_grid_storage_init, diag_grid_storage_end +use MOM_diag_mediator, only : diag_copy_diag_to_storage, diag_copy_storage_to_diag +use MOM_diag_mediator, only : diag_save_grids, diag_restore_grids +use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag, calc_Zint_diags +use MOM_diapyc_energy_req, only : diapyc_energy_req_init, diapyc_energy_req_end +use MOM_diapyc_energy_req, only : diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_CS +use MOM_CVMix_conv, only : CVMix_conv_init, CVMix_conv_cs +use MOM_CVMix_conv, only : CVMix_conv_end, calculate_CVMix_conv +use MOM_domains, only : pass_var, To_West, To_South, To_All, Omit_Corners +use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type +use MOM_tidal_mixing, only : tidal_mixing_init, tidal_mixing_cs +use MOM_tidal_mixing, only : tidal_mixing_end +use MOM_energetic_PBL, only : energetic_PBL, energetic_PBL_init +use MOM_energetic_PBL, only : energetic_PBL_end, energetic_PBL_CS +use MOM_energetic_PBL, only : energetic_PBL_get_MLD +use MOM_entrain_diffusive, only : entrainment_diffusive, entrain_diffusive_init +use MOM_entrain_diffusive, only : entrain_diffusive_end, entrain_diffusive_CS +use MOM_EOS, only : calculate_density, calculate_TFreeze +use MOM_EOS, only : calculate_specific_vol_derivs +use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery,MOM_mesg +use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint +use MOM_file_parser, only : get_param, log_version, param_file_type, read_param +use MOM_forcing_type, only : forcing, MOM_forcing_chksum +use MOM_forcing_type, only : calculateBuoyancyFlux2d, forcing_SinglePointPrint +use MOM_geothermal, only : geothermal, geothermal_init, geothermal_end, geothermal_CS +use MOM_grid, only : ocean_grid_type +use MOM_io, only : vardesc, var_desc +use MOM_int_tide_input, only : set_int_tide_input, int_tide_input_init +use MOM_int_tide_input, only : int_tide_input_end, int_tide_input_CS, int_tide_input_type +use MOM_interface_heights, only : find_eta +use MOM_internal_tides, only : propagate_int_tide +use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS +use MOM_kappa_shear, only : kappa_shear_is_used +use MOM_KPP, only : KPP_CS, KPP_init, KPP_compute_BLD, KPP_calculate +use MOM_KPP, only : KPP_end, KPP_get_BLD +use MOM_KPP, only : KPP_NonLocalTransport_temp, KPP_NonLocalTransport_saln +use MOM_opacity, only : opacity_init, set_opacity, opacity_end, opacity_CS +use MOM_regularize_layers, only : regularize_layers, regularize_layers_init, regularize_layers_CS +use MOM_set_diffusivity, only : set_diffusivity, set_BBL_TKE +use MOM_set_diffusivity, only : set_diffusivity_init, set_diffusivity_end +use MOM_set_diffusivity, only : set_diffusivity_CS +use MOM_shortwave_abs, only : absorbRemainingSW, optics_type +use MOM_sponge, only : apply_sponge, sponge_CS +use MOM_ALE_sponge, only : apply_ALE_sponge, ALE_sponge_CS +use MOM_time_manager, only : operator(-), set_time +use MOM_time_manager, only : operator(<=), time_type ! for testing itides (BDM) +use MOM_tracer_flow_control, only : call_tracer_column_fns, tracer_flow_control_CS +use MOM_tracer_diabatic, only : tracer_vertdiff +use MOM_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs +use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d +use MOM_verticalGrid, only : verticalGrid_type +use MOM_wave_speed, only : wave_speeds +use time_manager_mod, only : increment_time ! for testing itides (BDM) +use MOM_wave_interface, only : wave_parameters_CS +use MOM_diabatic_driver, only : diabatic_CS + +implicit none ; private + +#include + +public legacy_diabatic + +! clock ids +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 + +contains + +!> This subroutine imposes the diapycnal mass fluxes and the +!! accompanying diapycnal advection of momentum and tracers. +subroutine legacy_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & + G, GV, CS, WAVES) + type(ocean_grid_type), intent(inout) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< zonal velocity (m/s) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< meridional velocity (m/s) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness (m for Bouss / kg/m2 for non-Bouss) + type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields + !! unused have NULL ptrs + real, dimension(:,:), pointer :: Hml !< active mixed layer depth + type(forcing), intent(inout) :: fluxes !< points to forcing fields + !! unused fields have NULL ptrs + type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and + type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum + !! equations, to enable the later derived + !! diagnostics, like energy budgets + type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations + real, intent(in) :: dt !< time increment (seconds) + type(time_type), intent(in) :: Time_end !< Time at the end of the interval + type(diabatic_CS), pointer :: CS !< module control structure + type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves + + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & + ea, & ! amount of fluid entrained from the layer above within + ! one time step (m for Bouss, kg/m^2 for non-Bouss) + eb, & ! amount of fluid entrained from the layer below within + ! one time step (m for Bouss, kg/m^2 for non-Bouss) + Kd, & ! diapycnal diffusivity of layers (m^2/sec) + h_orig, & ! initial layer thicknesses (m for Bouss, kg/m^2 for non-Bouss) + h_prebound, & ! initial layer thicknesses (m for Bouss, kg/m^2 for non-Bouss) + hold, & ! layer thickness before diapycnal entrainment, and later + ! the initial layer thicknesses (if a mixed layer is used), + ! (m for Bouss, kg/m^2 for non-Bouss) + dSV_dT, & ! The partial derivatives of specific volume with temperature + dSV_dS, & ! and salinity in m^3/(kg K) and m^3/(kg ppt). + cTKE, & ! convective TKE requirements for each layer in J/m^2. + u_h, & ! zonal and meridional velocities at thickness points after + v_h ! entrainment (m/s) + + real, dimension(SZI_(G),SZJ_(G),CS%nMode) :: & + cn ! baroclinic gravity wave speeds (formerly cg1 - BDM) + + real, dimension(SZI_(G),SZJ_(G)) :: & + Rcv_ml, & ! coordinate density of mixed layer, used for applying sponges + SkinBuoyFlux! 2d surface buoyancy flux (m2/s3), used by ePBL + real, dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag ! diagnostic array for thickness + real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag ! diagnostic array for temp + real, dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag ! diagnostic array for salinity + real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn + real, dimension(SZI_(G),SZJ_(G)) :: TKE_itidal_input_test ! override of energy input for testing (BDM) + + real :: net_ent ! The net of ea-eb at an interface. + + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: & + ! These are targets so that the space can be shared with eaml & ebml. + eatr, & ! The equivalent of ea and eb for tracers, which differ from ea and + ebtr ! eb in that they tend to homogenize tracers in massless layers + ! near the boundaries (m for Bouss and kg/m^2 for non-Bouss) + + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), target :: & + Kd_int, & ! diapycnal diffusivity of interfaces (m^2/s) + Kd_heat, & ! diapycnal diffusivity of heat (m^2/s) + Kd_salt, & ! diapycnal diffusivity of salt and passive tracers (m^2/s) + Kd_ePBL, & ! test array of diapycnal diffusivities at interfaces (m^2/s) + eta, & ! Interface heights before diapycnal mixing, in m. + Tdif_flx, & ! diffusive diapycnal heat flux across interfaces (K m/s) + Tadv_flx, & ! advective diapycnal heat flux across interfaces (K m/s) + Sdif_flx, & ! diffusive diapycnal salt flux across interfaces (ppt m/s) + Sadv_flx ! advective diapycnal salt flux across interfaces (ppt m/s) + + ! The following 5 variables are only used with a bulk mixed layer. + real, pointer, dimension(:,:,:) :: & + eaml, & ! The equivalent of ea and eb due to mixed layer processes, + ebml ! (m for Bouss and kg/m^2 for non-Bouss). These will be + ! pointers to eatr and ebtr so as to reuse the memory as + ! the arrays are not needed at the same time. + + integer :: kb(SZI_(G),SZJ_(G)) ! index of the lightest layer denser + ! than the buffer laye (nondimensional) + + real :: p_ref_cv(SZI_(G)) ! Reference pressure for the potential + ! density which defines the coordinate + ! variable, set to P_Ref, in Pa. + + logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below, + ! where massive is defined as sufficiently thick that + ! the no-flux boundary conditions have not restricted + ! the entrainment - usually sqrt(Kd*dt). + + real :: b_denom_1 ! The first term in the denominator of b1 + ! (m for Bouss, kg/m^2 for non-Bouss) + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected + ! (m for Bouss and kg/m^2 for non-Bouss) + real :: h_neglect2 ! h_neglect^2 (m^2 for Bouss, kg^2/m^4 for non-Bouss) + real :: add_ent ! Entrainment that needs to be added when mixing tracers + ! (m for Bouss and kg/m^2 for non-Bouss) + real :: eaval ! eaval is 2*ea at velocity grid points (m for Bouss, kg/m^2 for non-Bouss) + real :: hval ! hval is 2*h at velocity grid points (m for Bouss, kg/m^2 for non-Bouss) + real :: h_tr ! h_tr is h at tracer points with a tiny thickness + ! added to ensure positive definiteness (m for Bouss, kg/m^2 for non-Bouss) + real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is + ! coupled to the bottom within a timestep (m) + + real :: htot(SZIB_(G)) ! The summed thickness from the bottom, in m. + real :: b1(SZIB_(G)), d1(SZIB_(G)) ! b1, c1, and d1 are variables used by the + real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver. + + real :: Ent_int ! The diffusive entrainment rate at an interface + ! (H units = m for Bouss, kg/m^2 for non-Bouss). + real :: dt_mix ! amount of time over which to apply mixing (seconds) + real :: Idt ! inverse time step (1/s) + + type(p3d) :: z_ptrs(7) ! pointers to diagnostics to be interpolated to depth + integer :: num_z_diags ! number of diagnostics to be interpolated to depth + integer :: z_ids(7) ! id numbers of diagnostics to be interpolated to depth + integer :: dir_flag ! An integer encoding the directions in which to do halo updates. + logical :: showCallTree ! If true, show the call tree + integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m + + integer :: ig, jg ! global indices for testing testing itide point source (BDM) + logical :: avg_enabled ! for testing internal tides (BDM) + real :: Kd_add_here ! An added diffusivity in m2/s + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + nkmb = GV%nk_rho_varies + 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") + + + ! Offer diagnostics of various state varables at the start of diabatic + ! these are mostly for debugging purposes. + if (CS%id_u_predia > 0) call post_data(CS%id_u_predia, u, CS%diag) + if (CS%id_v_predia > 0) call post_data(CS%id_v_predia, v, CS%diag) + if (CS%id_h_predia > 0) call post_data(CS%id_h_predia, h, CS%diag) + if (CS%id_T_predia > 0) call post_data(CS%id_T_predia, tv%T, CS%diag) + if (CS%id_S_predia > 0) call post_data(CS%id_S_predia, tv%S, CS%diag) + if (CS%id_e_predia > 0) then + call find_eta(h, tv, GV%g_Earth, G, GV, eta) + call post_data(CS%id_e_predia, eta, CS%diag) + endif + + + ! set equivalence between the same bits of memory for these arrays + eaml => eatr ; ebml => ebtr + + ! inverse time step + Idt = 1.0 / dt + + if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// & + "Module must be initialized before it is used.") + + if (CS%debug) then + call MOM_state_chksum("Start of diabatic ", u, v, h, G, GV, haloshift=0) + call MOM_forcing_chksum("Start of diabatic", fluxes, G, haloshift=0) + endif + if (CS%debugConservation) call MOM_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, G) + + if (CS%debug_energy_req) & + call diapyc_energy_req_test(h, dt, tv, G, GV, CS%diapyc_en_rec_CSp) + + + call cpu_clock_begin(id_clock_set_diffusivity) + call set_BBL_TKE(u, v, h, fluxes, visc, G, GV, CS%set_diff_CSp) + call cpu_clock_end(id_clock_set_diffusivity) + + ! Frazil formation keeps the temperature above the freezing point. + ! make_frazil is deliberately called at both the beginning and at + ! the end of the diabatic processes. + if (associated(tv%T) .AND. associated(tv%frazil)) then + ! For frazil diagnostic, the first call covers the first half of the time step + call enable_averaging(0.5*dt, Time_end - set_time(int(floor(0.5*dt+0.5))), CS%diag) + if (CS%frazil_tendency_diag) then + do k=1,nz ; do j=js,je ; do i=is,ie + temp_diag(i,j,k) = tv%T(i,j,k) + enddo ; enddo ; enddo + endif + + if (associated(fluxes%p_surf_full)) then + call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp, fluxes%p_surf_full) + else + call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp) + endif + if (showCallTree) call callTree_waypoint("done with 1st make_frazil (diabatic)") + + if (CS%frazil_tendency_diag) then + call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, G, GV, CS) + if (CS%id_frazil_h > 0) call post_data(CS%id_frazil_h, h, CS%diag) + endif + call disable_averaging(CS%diag) + endif + ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep + call enable_averaging(dt, Time_end, CS%diag) + if (CS%debugConservation) call MOM_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, G) + + if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then +!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_orig,h,eaml,ebml) + do k=1,nz ; do j=js,je ; do i=is,ie + h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0 + enddo ; enddo ; enddo + endif + + if (CS%use_geothermal) then + call cpu_clock_begin(id_clock_geothermal) + call geothermal(h, tv, dt, eaml, ebml, G, GV, CS%geothermal_CSp) + call cpu_clock_end(id_clock_geothermal) + if (showCallTree) call callTree_waypoint("geothermal (diabatic)") + if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G) + endif + + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_remap_grids(CS%diag) + + ! Set_opacity estimates the optical properties of the water column. + ! It will need to be modified later to include information about the + ! biological properties and layer thicknesses. + if (associated(CS%optics)) & + call set_opacity(CS%optics, fluxes, G, GV, CS%opacity_CSp) + + if (CS%bulkmixedlayer) then + if (CS%debug) then + call MOM_forcing_chksum("Before mixedlayer", fluxes, G, haloshift=0) + endif + + if (CS%ML_mix_first > 0.0) then +! This subroutine +! (1) Cools the mixed layer. +! (2) Performs convective adjustment by mixed layer entrainment. +! (3) Heats the mixed layer and causes it to detrain to +! Monin-Obukhov depth or minimum mixed layer depth. +! (4) Uses any remaining TKE to drive mixed layer entrainment. +! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) + call find_uv_at_h(u, v, h, u_h, v_h, G, GV) + + call cpu_clock_begin(id_clock_mixedlayer) + if (CS%ML_mix_first < 1.0) then + ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???) + call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*CS%ML_mix_first, & + eaml,ebml, G, GV, CS%bulkmixedlayer_CSp, CS%optics, & + Hml, CS%aggregate_FW_forcing, dt, last_call=.false.) + if (CS%salt_reject_below_ML) & + call insert_brine(h, tv, G, GV, fluxes, nkmb, CS%diabatic_aux_CSp, & + dt*CS%ML_mix_first, CS%id_brine_lay) + else + ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???) + call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, & + G, GV, CS%bulkmixedlayer_CSp, CS%optics, & + Hml, CS%aggregate_FW_forcing, dt, last_call=.true.) + endif + + ! Keep salinity from falling below a small but positive threshold. + ! This constraint is needed for SIS1 ice model, which can extract + ! more salt than is present in the ocean. SIS2 does not suffer + ! from this limitation, in which case we can let salinity=0 and still + ! have salt conserved with SIS2 ice. So for SIS2, we can run with + ! BOUND_SALINITY=False in MOM.F90. + if (associated(tv%S) .and. associated(tv%salt_deficit)) & + call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) + call cpu_clock_end(id_clock_mixedlayer) + if (CS%debug) then + call MOM_state_chksum("After mixedlayer ", u, v, h, G, GV, haloshift=0) + call MOM_forcing_chksum("After mixedlayer", fluxes, G, haloshift=0) + endif + if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)") + if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G) + endif + endif + + if (CS%debug) then + 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) + if (CS%debug) then + call hchksum(eaml, "after find_uv_at_h eaml",G%HI, scale=GV%H_to_m) + call hchksum(ebml, "after find_uv_at_h ebml",G%HI, scale=GV%H_to_m) + endif + else + call find_uv_at_h(u, v, h, u_h, v_h, G, GV) + endif + if (showCallTree) call callTree_waypoint("done with find_uv_at_h (diabatic)") + endif + + if (CS%use_int_tides) then + ! This block provides an interface for the unresolved low-mode internal + ! tide module (BDM). + + ! PROVIDE ENERGY DISTRIBUTION (calculate time-varying energy source) + call set_int_tide_input(u, v, h, tv, fluxes, CS%int_tide_input, dt, G, GV, & + CS%int_tide_input_CSp) + ! CALCULATE MODAL VELOCITIES + cn(:,:,:) = 0.0 + if (CS%uniform_cg) then + ! SET TO CONSTANT VALUE TO TEST PROPAGATE CODE + do m=1,CS%nMode ; cn(:,:,m) = CS%cg_test ; enddo + else + call wave_speeds(h, tv, G, GV, CS%nMode, cn, full_halos=.true.) + ! uncomment the lines below for a hard-coded cn that changes linearly with latitude + !do j=G%jsd,G%jed ; do i=G%isd,G%ied + ! cn(i,j,:) = ((7.-1.)/14000000.)*G%geoLatBu(i,j) + (1.-((7.-1.)/14000000.)*-7000000.) + !enddo ; enddo + endif + + if (CS%int_tide_source_test) then + ! BUILD 2D ARRAY WITH POINT SOURCE FOR TESTING + ! This block of code should be moved into set_int_tide_input. -RWH + TKE_itidal_input_test(:,:) = 0.0 + avg_enabled = query_averaging_enabled(CS%diag,time_end=CS%time_end) + if (CS%time_end <= CS%time_max_source) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + !INPUT ARBITRARY ENERGY POINT SOURCE + if ((G%idg_offset + i == CS%int_tide_source_x) .and. & + (G%jdg_offset + j == CS%int_tide_source_y)) then + TKE_itidal_input_test(i,j) = 1.0 + endif + enddo ; enddo + endif + ! CALL ROUTINE USING PRESCRIBED KE FOR TESTING + call propagate_int_tide(h, tv, cn, TKE_itidal_input_test, & + CS%int_tide_input%tideamp, CS%int_tide_input%Nb, dt, G, GV, CS%int_tide_CSp) + else + ! CALL ROUTINE USING CALCULATED KE INPUT + call propagate_int_tide(h, tv, cn, CS%int_tide_input%TKE_itidal_input, & + 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 + + call cpu_clock_begin(id_clock_set_diffusivity) + ! Sets: Kd, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S + ! Also changes: visc%Kd_shear, visc%TKE_turb (not clear that TKE_turb is used as input ???? + ! And sets visc%Kv_shear + call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, CS%set_diff_CSp, Kd, Kd_int) + call cpu_clock_end(id_clock_set_diffusivity) + if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)") + + if (CS%debug) then + call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, haloshift=0) + call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, haloshift=0) + call MOM_thermovar_chksum("after set_diffusivity ", tv, G) + call hchksum(Kd, "after set_diffusivity Kd",G%HI,haloshift=0) + call hchksum(Kd_Int, "after set_diffusivity Kd_Int",G%HI,haloshift=0) + endif + + + if (CS%useKPP) then + call cpu_clock_begin(id_clock_kpp) + ! KPP needs the surface buoyancy flux but does not update state variables. + ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes. + ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux + ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second) + ! unlike other instances where the fluxes are integrated in time over a time-step. + call calculateBuoyancyFlux2d(G, GV, fluxes, CS%optics, h, tv%T, tv%S, tv, & + CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) + ! The KPP scheme calculates boundary layer diffusivities and non-local transport. + ! MOM6 implementation of KPP matches the boundary layer to zero interior diffusivity, + ! since the matching to nonzero interior diffusivity can be problematic. + ! Changes: Kd_int. Sets: KPP_NLTheat, KPP_NLTscalar + +!$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat) +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_salt(i,j,k) = Kd_int(i,j,k) + Kd_heat(i,j,k) = Kd_int(i,j,k) + enddo ; enddo ; enddo + if (associated(visc%Kd_extra_S)) then +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_salt(i,j,k) = Kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k) + enddo ; enddo ; enddo + endif + if (associated(visc%Kd_extra_T)) then +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_heat(i,j,k) = Kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k) + enddo ; enddo ; enddo + endif +!$OMP end parallel + + call KPP_compute_BLD(CS%KPP_CSp, G, GV, h, tv%T, tv%S, u, v, tv%eqn_of_state, & + fluxes%ustar, CS%KPP_buoy_flux) + + call KPP_calculate(CS%KPP_CSp, G, GV, h, tv%T, tv%S, u, v, tv%eqn_of_state, & + fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, & + CS%KPP_NLTscalar, Waves=Waves) +!$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,G,Kd_heat,Hml) + + if (associated(Hml)) then + call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G) + call pass_var(Hml, G%domain, halo=1) + endif + + if (.not. CS%KPPisPassive) then +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_int(i,j,k) = min( Kd_salt(i,j,k), Kd_heat(i,j,k) ) + enddo ; enddo ; enddo + if (associated(visc%Kd_extra_S)) then +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + visc%Kd_extra_S(i,j,k) = Kd_salt(i,j,k) - Kd_int(i,j,k) + enddo ; enddo ; enddo + endif + if (associated(visc%Kd_extra_T)) then +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + visc%Kd_extra_T(i,j,k) = Kd_heat(i,j,k) - Kd_int(i,j,k) + enddo ; enddo ; enddo + endif + endif ! not passive +!$OMP end parallel + call cpu_clock_end(id_clock_kpp) + if (showCallTree) call callTree_waypoint("done with KPP_calculate (diabatic)") + if (CS%debug) then + call MOM_state_chksum("after KPP", u, v, h, G, GV, haloshift=0) + call MOM_forcing_chksum("after KPP", fluxes, G, haloshift=0) + call MOM_thermovar_chksum("after KPP", tv, G) + call hchksum(Kd, "after KPP Kd",G%HI,haloshift=0) + call hchksum(Kd_Int, "after KPP Kd_Int",G%HI,haloshift=0) + endif + + endif ! endif for KPP + + ! Add vertical diff./visc. due to convection (computed via CVMix) + if (CS%use_CVMix_conv) then + call calculate_CVMix_conv(h, tv, G, GV, CS%CVMix_conv_csp, Hml) + + !!!!!!!! GMM, the following needs to be checked !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + do k=1,nz ; do j=js,je ; do i=is,ie + Kd_int(i,j,k) = Kd_int(i,j,k) + CS%CVMix_conv_csp%kd_conv(i,j,k) + visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + CS%CVMix_conv_csp%kv_conv(i,j,k) + enddo ; enddo ; enddo + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + endif + + if (CS%useKPP) then + + call cpu_clock_begin(id_clock_kpp) + if (CS%debug) then + call hchksum(CS%KPP_temp_flux, "before KPP_applyNLT netHeat",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(CS%KPP_salt_flux, "before KPP_applyNLT netSalt",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(CS%KPP_NLTheat, "before KPP_applyNLT NLTheat",G%HI,haloshift=0) + call hchksum(CS%KPP_NLTscalar, "before KPP_applyNLT NLTscalar",G%HI,haloshift=0) + endif + ! Apply non-local transport of heat and salt + ! Changes: tv%T, tv%S + call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, CS%KPP_NLTheat, CS%KPP_temp_flux, dt, tv%T, tv%C_p) + call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, CS%KPP_NLTscalar, CS%KPP_salt_flux, dt, tv%S) + call cpu_clock_end(id_clock_kpp) + if (showCallTree) call callTree_waypoint("done with KPP_applyNonLocalTransport (diabatic)") + if (CS%debugConservation) call MOM_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, G) + + if (CS%debug) then + call MOM_state_chksum("after KPP_applyNLT ", u, v, h, G, GV, haloshift=0) + call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, haloshift=0) + call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G) + endif + + endif ! endif for KPP + + ! Differential diffusion done here. + ! Changes: tv%T, tv%S + ! If using matching within the KPP scheme, then this step needs to provide + ! a diffusivity and happen before KPP. But generally in MOM, we do not match + ! KPP boundary layer to interior, so this diffusivity can be computed when convenient. + if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T)) then + call cpu_clock_begin(id_clock_differential_diff) + + call differential_diffuse_T_S(h, tv, visc, dt, G, GV) + call cpu_clock_end(id_clock_differential_diff) + if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)") + if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G) + + ! increment heat and salt diffusivity. + ! CS%useKPP==.true. already has extra_T and extra_S included + if (.not. CS%useKPP) then + do K=2,nz ; do j=js,je ; do i=is,ie + Kd_heat(i,j,K) = Kd_heat(i,j,K) + visc%Kd_extra_T(i,j,K) + Kd_salt(i,j,K) = Kd_salt(i,j,K) + visc%Kd_extra_S(i,j,K) + enddo ; enddo ; enddo + endif + + + endif + + + ! This block sets ea, eb from Kd or Kd_int. + ! If using ALE algorithm, set ea=eb=Kd_int on interfaces for + ! use in the tri-diagonal solver. + ! Otherwise, call entrainment_diffusive() which sets ea and eb + ! based on KD and target densities (ie. does remapping as well). + if (CS%useALEalgorithm) then + + do j=js,je ; do i=is,ie + ea(i,j,1) = 0. + enddo ; enddo +!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_neglect,h,ea,GV,dt,Kd_int,eb) & +!$OMP private(hval) + do k=2,nz ; do j=js,je ; do i=is,ie + hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) + ea(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_int(i,j,k) + eb(i,j,k-1) = ea(i,j,k) + enddo ; enddo ; enddo + do j=js,je ; do i=is,ie + eb(i,j,nz) = 0. + enddo ; enddo + if (showCallTree) call callTree_waypoint("done setting ea,eb from Kd_int (diabatic)") + + else ! .not. CS%useALEalgorithm + ! When not using ALE, calculate layer entrainments/detrainments from + ! diffusivities and differences between layer and target densities + call cpu_clock_begin(id_clock_entrain) + ! Calculate appropriately limited diapycnal mass fluxes to account + ! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb + call Entrainment_diffusive(u, v, h, tv, fluxes, dt, G, GV, CS%entrain_diffusive_CSp, & + ea, eb, kb, Kd_Lay=Kd, Kd_int=Kd_int) + call cpu_clock_end(id_clock_entrain) + if (showCallTree) call callTree_waypoint("done with Entrainment_diffusive (diabatic)") + + endif ! endif for (CS%useALEalgorithm) + + if (CS%debug) then + call MOM_forcing_chksum("after calc_entrain ", fluxes, G, haloshift=0) + call MOM_thermovar_chksum("after calc_entrain ", tv, G) + call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, haloshift=0) + call hchksum(ea, "after calc_entrain ea", G%HI, haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "after calc_entrain eb", G%HI, haloshift=0, scale=GV%H_to_m) + endif + + ! Save fields before boundary forcing is applied for tendency diagnostics + if (CS%boundary_forcing_tendency_diag) then + do k=1,nz ; do j=js,je ; do i=is,ie + h_diag(i,j,k) = h(i,j,k) + temp_diag(i,j,k) = tv%T(i,j,k) + saln_diag(i,j,k) = tv%S(i,j,k) + enddo ; enddo ; enddo + endif + + ! Apply forcing when using the ALE algorithm + if (CS%useALEalgorithm) then + call cpu_clock_begin(id_clock_remap) + + ! Changes made to following fields: h, tv%T and tv%S. + + do k=1,nz ; do j=js,je ; do i=is,ie + h_prebound(i,j,k) = h(i,j,k) + enddo ; enddo ; enddo + if (CS%use_energetic_PBL) then + + skinbuoyflux(:,:) = 0.0 + call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & + h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) + + if (CS%debug) then + call hchksum(ea, "after applyBoundaryFluxes ea",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "after applyBoundaryFluxes eb",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(cTKE, "after applyBoundaryFluxes cTKE",G%HI,haloshift=0) + call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT",G%HI,haloshift=0) + call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS",G%HI,haloshift=0) + endif + + call find_uv_at_h(u, v, h, u_h, v_h, G, GV) + call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, & + CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) + + ! If visc%MLD exists, copy the ePBL's MLD into it + if (associated(visc%MLD)) then + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G) + call pass_var(visc%MLD, G%domain, halo=1) + Hml(:,:) = visc%MLD(:,:) + endif + + ! Augment the diffusivities due to those diagnosed in energetic_PBL. + do K=2,nz ; do j=js,je ; do i=is,ie + + if (CS%ePBL_is_additive) then + Kd_add_here = Kd_ePBL(i,j,K) + visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + Kd_ePBL(i,j,K) + else + Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_shear(i,j,K), 0.0) + visc%Kv_shear(i,j,K) = max(visc%Kv_shear(i,j,K), Kd_ePBL(i,j,K)) + endif + Ent_int = Kd_add_here * (GV%m_to_H**2 * dt) / & + (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) + eb(i,j,k-1) = eb(i,j,k-1) + Ent_int + ea(i,j,k) = ea(i,j,k) + Ent_int + Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here + + ! for diagnostics + Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_int(i,j,K) + Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_int(i,j,K) + + enddo ; enddo ; enddo + + if (CS%debug) then + call hchksum(ea, "after ePBL ea",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "after ePBL eb",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(Kd_ePBL, "after ePBL Kd_ePBL",G%HI,haloshift=0) + endif + + else + call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & + h, tv, CS%aggregate_FW_forcing, & + CS%evap_CFL_limit, CS%minimum_forcing_depth) + + endif ! endif for CS%use_energetic_PBL + + ! diagnose the tendencies due to boundary forcing + ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme + ! so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards + if (CS%boundary_forcing_tendency_diag) then + call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, G, GV, CS) + if (CS%id_boundary_forcing_h > 0) call post_data(CS%id_boundary_forcing_h, h, CS%diag, alt_h = h_diag) + endif + ! Boundary fluxes may have changed T, S, and h + call diag_update_remap_grids(CS%diag) + + call cpu_clock_end(id_clock_remap) + if (CS%debug) then + call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, haloshift=0) + call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G) + call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, haloshift=0) + endif + if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)") + if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, G) + + endif ! endif for (CS%useALEalgorithm) + + ! Update h according to divergence of the difference between + ! ea and eb. We keep a record of the original h in hold. + ! In the following, the checks for negative values are to guard + ! against instances where entrainment drives a layer to + ! negative thickness. This situation will never happen if + ! enough iterations are permitted in Calculate_Entrainment. + ! Even if too few iterations are allowed, it is still guarded + ! against. In other words the checks are probably unnecessary. + !$OMP parallel do default(shared) + do j=js,je + do i=is,ie + hold(i,j,1) = h(i,j,1) + h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2)) + hold(i,j,nz) = h(i,j,nz) + h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1)) + if (h(i,j,1) <= 0.0) then + h(i,j,1) = GV%Angstrom + endif + if (h(i,j,nz) <= 0.0) then + h(i,j,nz) = GV%Angstrom + endif + enddo + do k=2,nz-1 ; do i=is,ie + hold(i,j,k) = h(i,j,k) + h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + & + (eb(i,j,k) - ea(i,j,k+1))) + if (h(i,j,k) <= 0.0) then + h(i,j,k) = GV%Angstrom + endif + enddo ; enddo + enddo + ! Checks for negative thickness may have changed layer thicknesses + call diag_update_remap_grids(CS%diag) + + if (CS%debug) then + call MOM_state_chksum("after negative check ", u, v, h, G, GV, haloshift=0) + call MOM_forcing_chksum("after negative check ", fluxes, G, haloshift=0) + call MOM_thermovar_chksum("after negative check ", tv, G) + endif + if (showCallTree) call callTree_waypoint("done with h=ea-eb (diabatic)") + if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G) + + + ! Here, T and S are updated according to ea and eb. + ! If using the bulk mixed layer, T and S are also updated + ! by surface fluxes (in fluxes%*). + ! This is a very long block. + if (CS%bulkmixedlayer) then + + if (associated(tv%T)) then + call cpu_clock_begin(id_clock_tridiag) + ! Temperature and salinity (as state variables) are treated + ! differently from other tracers to insure massless layers that + ! are lighter than the mixed layer have temperatures and salinities + ! that correspond to their prescribed densities. + if (CS%massless_match_targets) then + !$OMP parallel do default (shared) private(h_tr,b1,d1,c1,b_denom_1) + do j=js,je + do i=is,ie + h_tr = hold(i,j,1) + h_neglect + b1(i) = 1.0 / (h_tr + eb(i,j,1)) + d1(i) = h_tr * b1(i) + tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1)) + tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1)) + enddo + do k=2,nkmb ; do i=is,ie + c1(i,k) = eb(i,j,k-1) * b1(i) + h_tr = hold(i,j,k) + h_neglect + b_denom_1 = h_tr + d1(i)*ea(i,j,k) + b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) + if (k kb(i,j)) then + c1(i,k) = eb(i,j,k-1) * b1(i) + h_tr = hold(i,j,k) + h_neglect + b_denom_1 = h_tr + d1(i)*ea(i,j,k) + b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) + d1(i) = b_denom_1 * b1(i) + tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1)) + tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1)) + elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j)) + ! The bottommost buffer layer might entrain all the mass from some + ! of the interior layers that are thin and lighter in the coordinate + ! density than that buffer layer. The T and S of these newly + ! massless interior layers are unchanged. + tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k) + tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k) + endif + enddo ; enddo + + do k=nz-1,nkmb,-1 ; do i=is,ie + if (k >= kb(i,j)) then + tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) + tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) + endif + enddo ; enddo + do i=is,ie ; if (kb(i,j) <= nz) then + tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j)) + tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j)) + endif ; enddo + do k=nkmb-1,1,-1 ; do i=is,ie + tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) + tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) + enddo ; enddo + enddo ! end of j loop + else ! .not. massless_match_targets + ! This simpler form allows T & S to be too dense for the layers + ! between the buffer layers and the interior. + ! Changes: T, S + if (CS%tracer_tridiag) then + call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV) + call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV) + else + call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) + endif + endif ! massless_match_targets + call cpu_clock_end(id_clock_tridiag) + + endif ! endif for associated(T) + if (CS%debugConservation) call MOM_state_stats('BML tridiag', u, v, h, tv%T, tv%S, G) + + if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then + ! The mixed layer code has already been called, but there is some needed + ! bookkeeping. + !$OMP parallel do default(shared) + do k=1,nz ; do j=js,je ; do i=is,ie + hold(i,j,k) = h_orig(i,j,k) + ea(i,j,k) = ea(i,j,k) + eaml(i,j,k) + eb(i,j,k) = eb(i,j,k) + ebml(i,j,k) + enddo ; enddo ; enddo + if (CS%debug) then + call hchksum(ea, "after ea = ea + eaml",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "after eb = eb + ebml",G%HI,haloshift=0, scale=GV%H_to_m) + endif + endif + + if (CS%ML_mix_first < 1.0) then + ! Call the mixed layer code now, perhaps for a second time. + ! This subroutine (1) Cools the mixed layer. + ! (2) Performs convective adjustment by mixed layer entrainment. + ! (3) Heats the mixed layer and causes it to detrain to + ! Monin-Obukhov depth or minimum mixed layer depth. + ! (4) Uses any remaining TKE to drive mixed layer entrainment. + ! (5) Possibly splits the buffer layer into two isopycnal layers. + + call find_uv_at_h(u, v, hold, u_h, v_h, G, GV, ea, eb) + if (CS%debug) call MOM_state_chksum("find_uv_at_h1 ", u, v, h, G, GV, haloshift=0) + + dt_mix = min(dt,dt*(1.0 - CS%ML_mix_first)) + call cpu_clock_begin(id_clock_mixedlayer) + ! Changes: h, tv%T, tv%S, ea and eb (G is also inout???) + call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, & + G, GV, CS%bulkmixedlayer_CSp, CS%optics, & + Hml, CS%aggregate_FW_forcing, dt, last_call=.true.) + + if (CS%salt_reject_below_ML) & + call insert_brine(h, tv, G, GV, fluxes, nkmb, CS%diabatic_aux_CSp, dt_mix, & + CS%id_brine_lay) + + ! Keep salinity from falling below a small but positive threshold. + ! This constraint is needed for SIS1 ice model, which can extract + ! more salt than is present in the ocean. SIS2 does not suffer + ! from this limitation, in which case we can let salinity=0 and still + ! have salt conserved with SIS2 ice. So for SIS2, we can run with + ! BOUND_SALINITY=False in MOM.F90. + if (associated(tv%S) .and. associated(tv%salt_deficit)) & + call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) + + call cpu_clock_end(id_clock_mixedlayer) + if (showCallTree) call callTree_waypoint("done with 2nd bulkmixedlayer (diabatic)") + if (CS%debugConservation) call MOM_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, G) + endif + + else ! following block for when NOT using BULKMIXEDLAYER + + + ! calculate change in temperature & salinity due to dia-coordinate surface diffusion + if (associated(tv%T)) then + + if (CS%debug) then + call hchksum(ea, "before triDiagTS ea ",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "before triDiagTS eb ",G%HI,haloshift=0, scale=GV%H_to_m) + endif + call cpu_clock_begin(id_clock_tridiag) + + ! Keep salinity from falling below a small but positive threshold. + ! This constraint is needed for SIS1 ice model, which can extract + ! more salt than is present in the ocean. SIS2 does not suffer + ! from this limitation, in which case we can let salinity=0 and still + ! have salt conserved with SIS2 ice. So for SIS2, we can run with + ! BOUND_SALINITY=False in MOM.F90. + if (associated(tv%S) .and. associated(tv%salt_deficit)) & + call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) + + if (CS%diabatic_diff_tendency_diag) then + do k=1,nz ; do j=js,je ; do i=is,ie + temp_diag(i,j,k) = tv%T(i,j,k) + saln_diag(i,j,k) = tv%S(i,j,k) + enddo ; enddo ; enddo + endif + + ! Changes T and S via the tridiagonal solver; no change to h + if (CS%tracer_tridiag) then + call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV) + call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV) + else + call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) + endif + + ! diagnose temperature, salinity, heat, and salt tendencies + ! Note: hold here refers to the thicknesses from before the dual-entraintment when using + ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will have changed + ! In either case, tendencies should be posted on hold + if (CS%diabatic_diff_tendency_diag) then + call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, G, GV, CS) + if (CS%id_diabatic_diff_h > 0) call post_data(CS%id_diabatic_diff_h, hold, CS%diag, alt_h = hold) + endif + + call cpu_clock_end(id_clock_tridiag) + if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)") + + endif ! endif corresponding to if (associated(tv%T)) + if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G) + + + endif ! endif for the BULKMIXEDLAYER block + + + if (CS%debug) then + call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, haloshift=0) + call MOM_thermovar_chksum("after mixed layer ", tv, G) + call hchksum(ea, "after mixed layer ea", G%HI, scale=GV%H_to_m) + call hchksum(eb, "after mixed layer eb", G%HI, scale=GV%H_to_m) + endif + + if (.not. CS%useALEalgorithm) then + call cpu_clock_begin(id_clock_remap) + call regularize_layers(h, tv, dt, ea, eb, G, GV, CS%regularize_layers_CSp) + call cpu_clock_end(id_clock_remap) + if (showCallTree) call callTree_waypoint("done with regularize_layers (diabatic)") + if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, tv%T, tv%S, G) + endif + + ! Whenever thickness changes let the diag manager know, as the + ! target grids for vertical remapping may need to be regenerated. + call diag_update_remap_grids(CS%diag) + + ! diagnostics + if ((CS%id_Tdif > 0) .or. (CS%id_Tdif_z > 0) .or. & + (CS%id_Tadv > 0) .or. (CS%id_Tadv_z > 0)) then + do j=js,je ; do i=is,ie + Tdif_flx(i,j,1) = 0.0 ; Tdif_flx(i,j,nz+1) = 0.0 + Tadv_flx(i,j,1) = 0.0 ; Tadv_flx(i,j,nz+1) = 0.0 + enddo ; enddo + !$OMP parallel do default(shared) + do K=2,nz ; do j=js,je ; do i=is,ie + Tdif_flx(i,j,K) = (Idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * & + (tv%T(i,j,k-1) - tv%T(i,j,k)) + Tadv_flx(i,j,K) = (Idt * (ea(i,j,k) - eb(i,j,k-1))) * & + 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k)) + enddo ; enddo ; enddo + endif + if ((CS%id_Sdif > 0) .or. (CS%id_Sdif_z > 0) .or. & + (CS%id_Sadv > 0) .or. (CS%id_Sadv_z > 0)) then + do j=js,je ; do i=is,ie + Sdif_flx(i,j,1) = 0.0 ; Sdif_flx(i,j,nz+1) = 0.0 + Sadv_flx(i,j,1) = 0.0 ; Sadv_flx(i,j,nz+1) = 0.0 + enddo ; enddo + !$OMP parallel do default(shared) + do K=2,nz ; do j=js,je ; do i=is,ie + Sdif_flx(i,j,K) = (Idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * & + (tv%S(i,j,k-1) - tv%S(i,j,k)) + Sadv_flx(i,j,K) = (Idt * (ea(i,j,k) - eb(i,j,k-1))) * & + 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k)) + enddo ; enddo ; enddo + endif + + ! mixing of passive tracers from massless boundary layers to interior + call cpu_clock_begin(id_clock_tracers) + if (CS%mix_boundary_tracers) then + Tr_ea_BBL = sqrt(dt*CS%Kd_BBL_tr) + !$OMP parallel do default(shared) private(htot,in_boundary,add_ent) + do j=js,je + do i=is,ie + ebtr(i,j,nz) = eb(i,j,nz) + htot(i) = 0.0 + in_boundary(i) = (G%mask2dT(i,j) > 0.0) + enddo + do k=nz,2,-1 ; do i=is,ie + if (in_boundary(i)) then + htot(i) = htot(i) + h(i,j,k) + ! If diapycnal mixing has been suppressed because this is a massless + ! layer near the bottom, add some mixing of tracers between these + ! layers. This flux is based on the harmonic mean of the two + ! thicknesses, as this corresponds pretty closely (to within + ! differences in the density jumps between layers) with what is done + ! in the calculation of the fluxes in the first place. Kd_min_tr + ! should be much less than the values that have been set in Kd, + ! perhaps a molecular diffusivity. + add_ent = ((dt * CS%Kd_min_tr) * GV%m_to_H**2) * & + ((h(i,j,k-1)+h(i,j,k)+h_neglect) / & + (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - & + 0.5*(ea(i,j,k) + eb(i,j,k-1)) + if (htot(i) < Tr_ea_BBL) then + add_ent = max(0.0, add_ent, & + (Tr_ea_BBL - htot(i)) - min(ea(i,j,k),eb(i,j,k-1))) + elseif (add_ent < 0.0) then + add_ent = 0.0 ; in_boundary(i) = .false. + endif + + ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent + eatr(i,j,k) = ea(i,j,k) + add_ent + else + ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k) + endif + if (associated(visc%Kd_extra_S)) then ; if (visc%Kd_extra_S(i,j,k) > 0.0) then + add_ent = ((dt * visc%Kd_extra_S(i,j,k)) * GV%m_to_H**2) / & + (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + & + h_neglect) + ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent + eatr(i,j,k) = eatr(i,j,k) + add_ent + endif ; endif + enddo ; enddo + do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo + + enddo + + if (CS%useALEalgorithm) then + ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied + ! so hold should be h_orig + call call_tracer_column_fns(h_prebound, h, ea, eb, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug, & + evap_CFL_limit = CS%evap_CFL_limit, & + minimum_forcing_depth = CS%minimum_forcing_depth) + else + call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug) + endif + + elseif (associated(visc%Kd_extra_S)) then ! extra diffusivity for passive tracers + + do j=js,je ; do i=is,ie + ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1) + enddo ; enddo + !$OMP parallel do default(shared) private(add_ent) + do k=nz,2,-1 ; do j=js,je ; do i=is,ie + if (visc%Kd_extra_S(i,j,k) > 0.0) then + add_ent = ((dt * visc%Kd_extra_S(i,j,k)) * GV%m_to_H**2) / & + (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + & + h_neglect) + else + add_ent = 0.0 + endif + ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent + eatr(i,j,k) = ea(i,j,k) + add_ent + enddo ; enddo ; enddo + + if (CS%useALEalgorithm) then + ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied + call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug,& + evap_CFL_limit = CS%evap_CFL_limit, & + minimum_forcing_depth = CS%minimum_forcing_depth) + else + call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug) + endif + + else + if (CS%useALEalgorithm) then + ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied + call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug, & + evap_CFL_limit = CS%evap_CFL_limit, & + minimum_forcing_depth = CS%minimum_forcing_depth) + else + call call_tracer_column_fns(hold, h, ea, eb, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug) + endif + + endif ! (CS%mix_boundary_tracers) + + + + call cpu_clock_end(id_clock_tracers) + + + ! sponges + if (CS%use_sponge) then + call cpu_clock_begin(id_clock_sponge) + if (associated(CS%ALE_sponge_CSp)) then + ! ALE sponge + call apply_ALE_sponge(h, dt, G, CS%ALE_sponge_CSp, CS%Time) + else + ! Layer mode sponge + if (CS%bulkmixedlayer .and. associated(tv%eqn_of_state)) then + do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo + !$OMP parallel do default(shared) + do j=js,je + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), & + is, ie-is+1, tv%eqn_of_state) + enddo + call apply_sponge(h, dt, G, GV, ea, eb, CS%sponge_CSp, Rcv_ml) + else + call apply_sponge(h, dt, G, GV, ea, eb, CS%sponge_CSp) + endif + endif + call cpu_clock_end(id_clock_sponge) + if (CS%debug) then + call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, haloshift=0) + call MOM_thermovar_chksum("apply_sponge ", tv, G) + endif + endif ! CS%use_sponge + + +! Save the diapycnal mass fluxes as a diagnostic field. + if (associated(CDp%diapyc_vel)) then + !$OMP parallel do default(shared) + do j=js,je + do K=2,nz ; do i=is,ie + CDp%diapyc_vel(i,j,K) = Idt * (GV%H_to_m * (ea(i,j,k) - eb(i,j,k-1))) + enddo ; enddo + do i=is,ie + CDp%diapyc_vel(i,j,1) = 0.0 + CDp%diapyc_vel(i,j,nz+1) = 0.0 + enddo + enddo + endif + +! For momentum, it is only the net flux that homogenizes within +! the mixed layer. Vertical viscosity that is proportional to the +! mixed layer turbulence is applied elsewhere. + if (CS%bulkmixedlayer) then + if (CS%debug) then + call hchksum(ea, "before net flux rearrangement ea",G%HI, scale=GV%H_to_m) + call hchksum(eb, "before net flux rearrangement eb",G%HI, scale=GV%H_to_m) + endif + !$OMP parallel do default(shared) private(net_ent) + do j=js,je + do K=2,GV%nkml ; do i=is,ie + net_ent = ea(i,j,k) - eb(i,j,k-1) + ea(i,j,k) = max(net_ent, 0.0) + eb(i,j,k-1) = max(-net_ent, 0.0) + enddo ; enddo + enddo + if (CS%debug) then + call hchksum(ea, "after net flux rearrangement ea",G%HI, scale=GV%H_to_m) + call hchksum(eb, "after net flux rearrangement eb",G%HI, scale=GV%H_to_m) + endif + endif + +! Initialize halo regions of ea, eb, and hold to default values. + !$OMP parallel do default(shared) + do k=1,nz + do i=is-1,ie+1 + hold(i,js-1,k) = GV%Angstrom ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0 + hold(i,je+1,k) = GV%Angstrom ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0 + enddo + do j=js,je + hold(is-1,j,k) = GV%Angstrom ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0 + hold(ie+1,j,k) = GV%Angstrom ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0 + enddo + enddo + + call cpu_clock_begin(id_clock_pass) + if (G%symmetric) then ; dir_flag = To_All+Omit_Corners + else ; dir_flag = To_West+To_South+Omit_Corners ; endif + call create_group_pass(CS%pass_hold_eb_ea, hold, G%Domain, dir_flag, halo=1) + call create_group_pass(CS%pass_hold_eb_ea, eb, G%Domain, dir_flag, halo=1) + call create_group_pass(CS%pass_hold_eb_ea, ea, G%Domain, dir_flag, halo=1) + call do_group_pass(CS%pass_hold_eb_ea, G%Domain) + ! visc%Kv_shear is not in the group pass because it has larger vertical extent. + if (associated(visc%Kv_shear)) & + call pass_var(visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) + call cpu_clock_end(id_clock_pass) + + if (.not. CS%useALEalgorithm) then + ! Use a tridiagonal solver to determine effect of the diapycnal + ! advection on velocity field. It is assumed that water leaves + ! or enters the ocean with the surface velocity. + if (CS%debug) then + call MOM_state_chksum("before u/v tridiag ", u, v, h, G, GV, haloshift=0) + call hchksum(ea, "before u/v tridiag ea",G%HI, scale=GV%H_to_m) + call hchksum(eb, "before u/v tridiag eb",G%HI, scale=GV%H_to_m) + call hchksum(hold, "before u/v tridiag hold",G%HI, scale=GV%H_to_m) + endif + call cpu_clock_begin(id_clock_tridiag) + !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval) + do j=js,je + do I=Isq,Ieq + if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,1) = u(I,j,1) + hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect + b1(I) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1))) + d1(I) = hval * b1(I) + u(I,j,1) = b1(I) * (hval * u(I,j,1)) + enddo + do k=2,nz ; do I=Isq,Ieq + if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,k) = u(I,j,k) + c1(I,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(I) + eaval = ea(i,j,k) + ea(i+1,j,k) + hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect + b1(I) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(I)*eaval)) + d1(I) = (hval + d1(I)*eaval) * b1(I) + u(I,j,k) = (hval*u(I,j,k) + eaval*u(I,j,k-1))*b1(I) + enddo ; enddo + do k=nz-1,1,-1 ; do I=Isq,Ieq + u(I,j,k) = u(I,j,k) + c1(I,k+1)*u(I,j,k+1) + if (associated(ADp%du_dt_dia)) & + ADp%du_dt_dia(I,j,k) = (u(I,j,k) - ADp%du_dt_dia(I,j,k)) * Idt + enddo ; enddo + if (associated(ADp%du_dt_dia)) then + do I=Isq,Ieq + ADp%du_dt_dia(I,j,nz) = (u(I,j,nz)-ADp%du_dt_dia(I,j,nz)) * Idt + enddo + endif + enddo + if (CS%debug) then + call MOM_state_chksum("aft 1st loop tridiag ", u, v, h, G, GV, haloshift=0) + endif + !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval) + do J=Jsq,Jeq + do i=is,ie + if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,1) = v(i,J,1) + hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect + b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1))) + d1(I) = hval * b1(I) + v(i,J,1) = b1(i) * (hval * v(i,J,1)) + enddo + do k=2,nz ; do i=is,ie + if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,k) = v(i,J,k) + c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i) + eaval = ea(i,j,k) + ea(i,j+1,k) + hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect + b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval)) + d1(i) = (hval + d1(i)*eaval) * b1(i) + v(i,J,k) = (hval*v(i,J,k) + eaval*v(i,J,k-1))*b1(i) + enddo ; enddo + do k=nz-1,1,-1 ; do i=is,ie + v(i,J,k) = v(i,J,k) + c1(i,k+1)*v(i,J,k+1) + if (associated(ADp%dv_dt_dia)) & + ADp%dv_dt_dia(i,J,k) = (v(i,J,k) - ADp%dv_dt_dia(i,J,k)) * Idt + enddo ; enddo + if (associated(ADp%dv_dt_dia)) then + do i=is,ie + ADp%dv_dt_dia(i,J,nz) = (v(i,J,nz)-ADp%dv_dt_dia(i,J,nz)) * Idt + enddo + endif + enddo + call cpu_clock_end(id_clock_tridiag) + if (CS%debug) then + call MOM_state_chksum("after u/v tridiag ", u, v, h, G, GV, haloshift=0) + endif + endif ! useALEalgorithm + + call disable_averaging(CS%diag) + ! Frazil formation keeps temperature above the freezing point. + ! make_frazil is deliberately called at both the beginning and at + ! the end of the diabatic processes. + if (associated(tv%T) .AND. associated(tv%frazil)) then + call enable_averaging(0.5*dt, Time_end, CS%diag) + if (CS%frazil_tendency_diag) then + do k=1,nz ; do j=js,je ; do i=is,ie + temp_diag(i,j,k) = tv%T(i,j,k) + enddo ; enddo ; enddo + endif + + if (associated(fluxes%p_surf_full)) then + call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp, fluxes%p_surf_full) + else + call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp) + endif + + if (CS%frazil_tendency_diag) then + call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, G, GV, CS) + if (CS%id_frazil_h > 0 ) call post_data(CS%id_frazil_h, h, CS%diag) + endif + + if (showCallTree) call callTree_waypoint("done with 2nd make_frazil (diabatic)") + if (CS%debugConservation) call MOM_state_stats('2nd make_frazil', u, v, h, tv%T, tv%S, G) + call disable_averaging(CS%diag) + + endif ! endif for frazil + + ! Diagnose the diapycnal diffusivities and other related quantities. + call enable_averaging(dt, Time_end, CS%diag) + + if (CS%id_Kd_interface > 0) call post_data(CS%id_Kd_interface, Kd_int, CS%diag) + if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) + if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) + if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag) + + if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag) + if (CS%id_eb > 0) call post_data(CS%id_eb, eb, CS%diag) + + if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, ADp%du_dt_dia, CS%diag) + if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, ADp%dv_dt_dia, CS%diag) + if (CS%id_wd > 0) call post_data(CS%id_wd, CDp%diapyc_vel, CS%diag) + + if (CS%id_MLD_003 > 0 .or. CS%id_subMLN2 > 0 .or. CS%id_mlotstsq > 0) then + call diagnoseMLDbyDensityDifference(CS%id_MLD_003, h, tv, 0.03, G, GV, CS%diag, & + id_N2subML=CS%id_subMLN2, id_MLDsq=CS%id_mlotstsq) + endif + if (CS%id_MLD_0125 > 0) then + call diagnoseMLDbyDensityDifference(CS%id_MLD_0125, h, tv, 0.125, G, GV, CS%diag) + endif + if (CS%id_MLD_user > 0) then + call diagnoseMLDbyDensityDifference(CS%id_MLD_user, h, tv, CS%MLDdensityDifference, G, GV, CS%diag) + endif + + if (CS%id_Tdif > 0) call post_data(CS%id_Tdif, Tdif_flx, CS%diag) + if (CS%id_Tadv > 0) call post_data(CS%id_Tadv, Tadv_flx, CS%diag) + if (CS%id_Sdif > 0) call post_data(CS%id_Sdif, Sdif_flx, CS%diag) + if (CS%id_Sadv > 0) call post_data(CS%id_Sadv, Sadv_flx, CS%diag) + if (CS%use_int_tides) then + if (CS%id_cg1 > 0) call post_data(CS%id_cg1, cn(:,:,1),CS%diag) + do m=1,CS%nMode + if (CS%id_cn(m) > 0) call post_data(CS%id_cn(m),cn(:,:,m),CS%diag) + enddo + endif + + call disable_averaging(CS%diag) + + num_z_diags = 0 + if (CS%id_Kd_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_Kd_z ; z_ptrs(num_z_diags)%p => Kd_int + endif + if (CS%id_Tdif_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_Tdif_z ; z_ptrs(num_z_diags)%p => Tdif_flx + endif + if (CS%id_Tadv_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_Tadv_z ; z_ptrs(num_z_diags)%p => Tadv_flx + endif + if (CS%id_Sdif_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_Sdif_z ; z_ptrs(num_z_diags)%p => Sdif_flx + endif + if (CS%id_Sadv_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_Sadv_z ; z_ptrs(num_z_diags)%p => Sadv_flx + endif + + if (num_z_diags > 0) & + call calc_Zint_diags(h, z_ptrs, z_ids, num_z_diags, G, GV, CS%diag_to_Z_CSp) + + if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G) + if (showCallTree) call callTree_leave("diabatic()") + +end subroutine legacy_diabatic + +!> This routine diagnoses tendencies from application of diabatic diffusion +!! using ALE algorithm. Note that layer thickness is not altered by +!! diabatic diffusion. +subroutine diagnose_diabatic_diff_tendency(tv, h, temp_old, saln_old, dt, G, GV, CS) + type(ocean_grid_type), intent(in) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< thickness (m or kg/m2) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: temp_old !< temperature prior to diabatic physics + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: saln_old !< salinity prior to diabatic physics (PPT) + real, intent(in) :: dt !< time step (sec) + type(diabatic_CS), pointer :: CS !< module control structure + + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d + real, dimension(SZI_(G),SZJ_(G)) :: work_2d + real :: Idt + integer :: i, j, k, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + Idt = 1/dt + work_3d(:,:,:) = 0.0 + work_2d(:,:) = 0.0 + + + ! temperature tendency + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*Idt + enddo ; enddo ; enddo + if (CS%id_diabatic_diff_temp_tend > 0) then + call post_data(CS%id_diabatic_diff_temp_tend, work_3d, CS%diag, alt_h = h) + endif + + ! heat tendency + if (CS%id_diabatic_diff_heat_tend > 0 .or. CS%id_diabatic_diff_heat_tend_2d > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = h(i,j,k) * GV%H_to_kg_m2 * tv%C_p * work_3d(i,j,k) + enddo ; enddo ; enddo + if (CS%id_diabatic_diff_heat_tend > 0) then + call post_data(CS%id_diabatic_diff_heat_tend, work_3d, CS%diag, alt_h = h) + endif + if (CS%id_diabatic_diff_heat_tend_2d > 0) then + do j=js,je ; do i=is,ie + work_2d(i,j) = 0.0 + do k=1,nz + work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k) + enddo + enddo ; enddo + call post_data(CS%id_diabatic_diff_heat_tend_2d, work_2d, CS%diag) + endif + endif + + ! salinity tendency + if (CS%id_diabatic_diff_saln_tend > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*Idt + enddo ; enddo ; enddo + call post_data(CS%id_diabatic_diff_saln_tend, work_3d, CS%diag, alt_h = h) + endif + + ! salt tendency + if (CS%id_diabatic_diff_salt_tend > 0 .or. CS%id_diabatic_diff_salt_tend_2d > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = h(i,j,k) * GV%H_to_kg_m2 * CS%ppt2mks * work_3d(i,j,k) + enddo ; enddo ; enddo + if (CS%id_diabatic_diff_salt_tend > 0) then + call post_data(CS%id_diabatic_diff_salt_tend, work_3d, CS%diag, alt_h = h) + endif + if (CS%id_diabatic_diff_salt_tend_2d > 0) then + do j=js,je ; do i=is,ie + work_2d(i,j) = 0.0 + do k=1,nz + work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k) + enddo + enddo ; enddo + call post_data(CS%id_diabatic_diff_salt_tend_2d, work_2d, CS%diag) + endif + endif + +end subroutine diagnose_diabatic_diff_tendency + + +!> This routine diagnoses tendencies from application of boundary fluxes. +!! These impacts are generally 3d, in particular for penetrative shortwave. +!! Other fluxes contribute 3d in cases when the layers vanish or are very thin, +!! in which case we distribute the flux into k > 1 layers. +subroutine diagnose_boundary_forcing_tendency(tv, h, temp_old, saln_old, h_old, & + dt, G, GV, CS) + type(ocean_grid_type), intent(in) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< thickness after boundary flux application (m or kg/m2) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: temp_old !< temperature prior to boundary flux application + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: saln_old !< salinity prior to boundary flux application (PPT) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h_old !< thickness prior to boundary flux application (m or kg/m2) + real, intent(in) :: dt !< time step (sec) + type(diabatic_CS), pointer :: CS !< module control structure + + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d + real, dimension(SZI_(G),SZJ_(G)) :: work_2d + real :: Idt + integer :: i, j, k, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + Idt = 1/dt + work_3d(:,:,:) = 0.0 + work_2d(:,:) = 0.0 + + ! Thickness tendency + if (CS%id_boundary_forcing_h_tendency > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = (h(i,j,k) - h_old(i,j,k))*Idt + enddo ; enddo ; enddo + call post_data(CS%id_boundary_forcing_h_tendency, work_3d, CS%diag, alt_h = h_old) + endif + + ! temperature tendency + if (CS%id_boundary_forcing_temp_tend > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*Idt + enddo ; enddo ; enddo + call post_data(CS%id_boundary_forcing_temp_tend, work_3d, CS%diag, alt_h = h_old) + endif + + ! heat tendency + if (CS%id_boundary_forcing_heat_tend > 0 .or. CS%id_boundary_forcing_heat_tend_2d > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = GV%H_to_kg_m2 * tv%C_p * Idt * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * temp_old(i,j,k)) + enddo ; enddo ; enddo + if (CS%id_boundary_forcing_heat_tend > 0) then + call post_data(CS%id_boundary_forcing_heat_tend, work_3d, CS%diag, alt_h = h_old) + endif + if (CS%id_boundary_forcing_heat_tend_2d > 0) then + do j=js,je ; do i=is,ie + work_2d(i,j) = 0.0 + do k=1,nz + work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k) + enddo + enddo ; enddo + call post_data(CS%id_boundary_forcing_heat_tend_2d, work_2d, CS%diag) + endif + endif + + ! salinity tendency + if (CS%id_boundary_forcing_saln_tend > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*Idt + enddo ; enddo ; enddo + call post_data(CS%id_boundary_forcing_saln_tend, work_3d, CS%diag, alt_h = h_old) + endif + + ! salt tendency + if (CS%id_boundary_forcing_salt_tend > 0 .or. CS%id_boundary_forcing_salt_tend_2d > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = GV%H_to_kg_m2 * CS%ppt2mks * Idt * (h(i,j,k) * tv%S(i,j,k) - h_old(i,j,k) * saln_old(i,j,k)) + enddo ; enddo ; enddo + if (CS%id_boundary_forcing_salt_tend > 0) then + call post_data(CS%id_boundary_forcing_salt_tend, work_3d, CS%diag, alt_h = h_old) + endif + if (CS%id_boundary_forcing_salt_tend_2d > 0) then + do j=js,je ; do i=is,ie + work_2d(i,j) = 0.0 + do k=1,nz + work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k) + enddo + enddo ; enddo + call post_data(CS%id_boundary_forcing_salt_tend_2d, work_2d, CS%diag) + endif + endif + +end subroutine diagnose_boundary_forcing_tendency + + +!> This routine diagnoses tendencies for temperature and heat from frazil formation. +!! This routine is called twice from within subroutine diabatic; at start and at +!! end of the diabatic processes. The impacts from frazil are generally a function +!! of depth. Hence, when checking heat budget, be sure to remove HFSIFRAZIL from HFDS in k=1. +subroutine diagnose_frazil_tendency(tv, h, temp_old, dt, G, GV, CS) + type(ocean_grid_type), intent(in) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(diabatic_CS), pointer :: CS !< module control structure + type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< thickness (m or kg/m2) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: temp_old !< temperature prior to frazil formation + real, intent(in) :: dt !< time step (sec) + + real, dimension(SZI_(G),SZJ_(G)) :: work_2d + real :: Idt + integer :: i, j, k, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + Idt = 1/dt + + ! temperature tendency + if (CS%id_frazil_temp_tend > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + CS%frazil_temp_diag(i,j,k) = Idt * (tv%T(i,j,k)-temp_old(i,j,k)) + enddo ; enddo ; enddo + call post_data(CS%id_frazil_temp_tend, CS%frazil_temp_diag(:,:,:), CS%diag) + endif + + ! heat tendency + if (CS%id_frazil_heat_tend > 0 .or. CS%id_frazil_heat_tend_2d > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + CS%frazil_heat_diag(i,j,k) = GV%H_to_kg_m2 * tv%C_p * h(i,j,k) * Idt * (tv%T(i,j,k)-temp_old(i,j,k)) + enddo ; enddo ; enddo + if (CS%id_frazil_heat_tend > 0) call post_data(CS%id_frazil_heat_tend, CS%frazil_heat_diag(:,:,:), CS%diag) + + ! As a consistency check, we must have + ! FRAZIL_HEAT_TENDENCY_2d = HFSIFRAZIL + if (CS%id_frazil_heat_tend_2d > 0) then + do j=js,je ; do i=is,ie + work_2d(i,j) = 0.0 + do k=1,nz + work_2d(i,j) = work_2d(i,j) + CS%frazil_heat_diag(i,j,k) + enddo + enddo ; enddo + call post_data(CS%id_frazil_heat_tend_2d, work_2d, CS%diag) + endif + endif + +end subroutine diagnose_frazil_tendency + + +!> \namespace mom_diabatic_driver +!! +!! By Robert Hallberg, 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. +!! +!! \section section_diabatic Outline of MOM diabatic +!! +!! * 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 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 treated +!! 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. + +end module MOM_legacy_diabatic_driver diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 9906083597..903868795a 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -23,6 +23,8 @@ module MOM_set_diffusivity use MOM_kappa_shear, only : calculate_kappa_shear, kappa_shear_init, Kappa_shear_CS use MOM_CVMix_shear, only : calculate_CVMix_shear, CVMix_shear_init, CVMix_shear_cs use MOM_CVMix_shear, only : CVMix_shear_end +use MOM_CVMix_ddiff, only : CVMix_ddiff_init, CVMix_ddiff_end, CVMix_ddiff_cs +use MOM_CVMix_ddiff, only : compute_ddiff_coeffs use MOM_bkgnd_mixing, only : calculate_bkgnd_mixing, bkgnd_mixing_init, bkgnd_mixing_cs use MOM_bkgnd_mixing, only : bkgnd_mixing_end, sfc_bkgnd_mixing use MOM_string_functions, only : uppercase @@ -43,104 +45,101 @@ module MOM_set_diffusivity public set_diffusivity_end type, public :: set_diffusivity_CS ; private - logical :: debug ! If true, write verbose checksums for debugging. - - logical :: bulkmixedlayer ! If true, a refined bulk mixed layer is used with - ! GV%nk_rho_varies variable density mixed & buffer - ! layers. - real :: FluxRi_max ! The flux Richardson number where the stratification is - ! large enough that N2 > omega2. The full expression for - ! the Flux Richardson number is usually - ! FLUX_RI_MAX*N2/(N2+OMEGA2). The default is 0.2. - logical :: bottomdraglaw ! If true, the bottom stress is calculated with a - ! drag law c_drag*|u|*u. - logical :: BBL_mixing_as_max ! If true, take the maximum of the diffusivity - ! from the BBL mixing and the other diffusivities. - ! Otherwise, diffusivities from the BBL_mixing is - ! added. - logical :: use_LOTW_BBL_diffusivity ! If true, use simpler/less precise, BBL diffusivity. - logical :: LOTW_BBL_use_omega ! If true, use simpler/less precise, BBL diffusivity. - real :: BBL_effic ! efficiency with which the energy extracted - ! by bottom drag drives BBL diffusion (nondim) - real :: cdrag ! quadratic drag coefficient (nondim) - real :: IMax_decay ! inverse of a maximum decay scale for - ! bottom-drag driven turbulence, (1/m) - - real :: Kd ! interior diapycnal diffusivity (m2/s) - real :: Kd_min ! minimum diapycnal diffusivity (m2/s) - real :: Kd_max ! maximum increment for diapycnal diffusivity (m2/s) - ! Set to a negative value to have no limit. - real :: Kd_add ! uniform diffusivity added everywhere without - ! filtering or scaling (m2/s) - real :: Kv ! interior vertical viscosity (m2/s) - real :: Kdml ! mixed layer diapycnal diffusivity (m2/s) - ! when bulkmixedlayer==.false. - real :: Hmix ! mixed layer thickness (meter) when - ! bulkmixedlayer==.false. + logical :: debug !< If true, write verbose checksums for debugging. + + logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with + !! GV%nk_rho_varies variable density mixed & buffer + !! layers. + real :: FluxRi_max !< The flux Richardson number where the stratification is + !! large enough that N2 > omega2. The full expression for + !! the Flux Richardson number is usually + !! FLUX_RI_MAX*N2/(N2+OMEGA2). The default is 0.2. + logical :: bottomdraglaw !< If true, the bottom stress is calculated with a + !! drag law c_drag*|u|*u. + logical :: BBL_mixing_as_max !< If true, take the maximum of the diffusivity + !! from the BBL mixing and the other diffusivities. + !! Otherwise, diffusivities from the BBL_mixing is + !! added. + logical :: use_LOTW_BBL_diffusivity !< If true, use simpler/less precise, BBL diffusivity. + logical :: LOTW_BBL_use_omega !< If true, use simpler/less precise, BBL diffusivity. + real :: BBL_effic !< efficiency with which the energy extracted + !! by bottom drag drives BBL diffusion (nondim) + real :: cdrag !< quadratic drag coefficient (nondim) + real :: IMax_decay !< inverse of a maximum decay scale for + !! bottom-drag driven turbulence, (1/m) + real :: Kv !< The interior vertical viscosity (m2/s) + real :: Kd !< interior diapycnal diffusivity (m2/s) + real :: Kd_min !< minimum diapycnal diffusivity (m2/s) + real :: Kd_max !< maximum increment for diapycnal diffusivity (m2/s) + !! Set to a negative value to have no limit. + real :: Kd_add !< uniform diffusivity added everywhere without + !! filtering or scaling (m2/s) + real :: Kdml !< mixed layer diapycnal diffusivity (m2/s) + !! when bulkmixedlayer==.false. + real :: Hmix !< mixed layer thickness (meter) when + !! bulkmixedlayer==.false. type(diag_ctrl), pointer :: diag ! structure to regulate diagn output timing - logical :: limit_dissipation ! If enabled, dissipation is limited to be larger - ! than the following: - real :: dissip_min ! Minimum dissipation (W/m3) - real :: dissip_N0 ! Coefficient a in minimum dissipation = a+b*N (W/m3) - real :: dissip_N1 ! Coefficient b in minimum dissipation = a+b*N (J/m3) - real :: dissip_N2 ! Coefficient c in minimum dissipation = c*N2 (W m-3 s2) - real :: dissip_Kd_min ! Minimum Kd (m2/s) with dissipatio Rho0*Kd_min*N^2 - - real :: TKE_itide_max ! maximum internal tide conversion (W m-2) - ! available to mix above the BBL - real :: omega ! Earth's rotation frequency (s-1) - logical :: ML_radiation ! allow a fraction of TKE available from wind work - ! to penetrate below mixed layer base with a vertical - ! decay scale determined by the minimum of - ! (1) The depth of the mixed layer, or - ! (2) An Ekman length scale. - ! Energy availble to drive mixing below the mixed layer is - ! given by E = ML_RAD_COEFF*MSTAR*USTAR**3. Optionally, if - ! ML_rad_TKE_decay is true, this is further reduced by a factor - ! of exp(-h_ML*Idecay_len_TkE), where Idecay_len_TKE is - ! calculated the same way as in the mixed layer code. - ! The diapycnal diffusivity is KD(k) = E/(N2(k)+OMEGA2), - ! where N2 is the squared buoyancy frequency (s-2) and OMEGA2 - ! is the rotation rate of the earth squared. - real :: ML_rad_kd_max ! Maximum diapycnal diffusivity due to turbulence - ! radiated from the base of the mixed layer (m2/s) - real :: ML_rad_efold_coeff ! non-dim coefficient to scale penetration depth - real :: ML_rad_coeff ! coefficient, which scales MSTAR*USTAR^3 to - ! obtain energy available for mixing below - ! mixed layer base (nondimensional) - logical :: ML_rad_TKE_decay ! If true, apply same exponential decay - ! to ML_rad as applied to the other surface - ! sources of TKE in the mixed layer code. - real :: ustar_min ! A minimum value of ustar to avoid numerical - ! problems (m/s). If the value is small enough, - ! this parameter should not affect the solution. - real :: TKE_decay ! ratio of natural Ekman depth to TKE decay scale (nondim) - real :: mstar ! ratio of friction velocity cubed to - ! TKE input to the mixed layer (nondim) - logical :: ML_use_omega ! If true, use absolute rotation rate instead - ! of the vertical component of rotation when - ! setting the decay scale for mixed layer turbulence. - real :: ML_omega_frac ! When setting the decay scale for turbulence, use - ! this fraction of the absolute rotation rate blended - ! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2. - logical :: user_change_diff ! If true, call user-defined code to change diffusivity. - logical :: useKappaShear ! If true, use the kappa_shear module to find the - ! shear-driven diapycnal diffusivity. - logical :: use_CVMix_shear ! If true, use one of the CVMix modules to find - ! shear-driven diapycnal diffusivity. - logical :: double_diffusion ! If true, enable double-diffusive mixing. - logical :: simple_TKE_to_Kd ! If true, uses a simple estimate of Kd/TKE that - ! does not rely on a layer-formulation. - real :: Max_Rrho_salt_fingers ! max density ratio for salt fingering - real :: Max_salt_diff_salt_fingers ! max salt diffusivity for salt fingers (m2/s) - real :: Kv_molecular ! molecular visc for double diff convect (m2/s) + logical :: limit_dissipation !< If enabled, dissipation is limited to be larger + !! than the following: + real :: dissip_min !< Minimum dissipation (W/m3) + real :: dissip_N0 !< Coefficient a in minimum dissipation = a+b*N (W/m3) + real :: dissip_N1 !< Coefficient b in minimum dissipation = a+b*N (J/m3) + real :: dissip_N2 !< Coefficient c in minimum dissipation = c*N2 (W m-3 s2) + real :: dissip_Kd_min !< Minimum Kd (m2/s) with dissipatio Rho0*Kd_min*N^2 + + real :: TKE_itide_max !< maximum internal tide conversion (W m-2) + !! available to mix above the BBL + real :: omega !< Earth's rotation frequency (s-1) + logical :: ML_radiation !< allow a fraction of TKE available from wind work + !! to penetrate below mixed layer base with a vertical + !! decay scale determined by the minimum of + !! (1) The depth of the mixed layer, or + !! (2) An Ekman length scale. + !! Energy availble to drive mixing below the mixed layer is + !! given by E = ML_RAD_COEFF*MSTAR*USTAR**3. Optionally, if + !! ML_rad_TKE_decay is true, this is further reduced by a factor + !! of exp(-h_ML*Idecay_len_TkE), where Idecay_len_TKE is + !! calculated the same way as in the mixed layer code. + !! The diapycnal diffusivity is KD(k) = E/(N2(k)+OMEGA2), + !! where N2 is the squared buoyancy frequency (s-2) and OMEGA2 + !! is the rotation rate of the earth squared. + real :: ML_rad_kd_max !< Maximum diapycnal diffusivity due to turbulence + !! radiated from the base of the mixed layer (m2/s) + real :: ML_rad_efold_coeff !< non-dim coefficient to scale penetration depth + real :: ML_rad_coeff !< coefficient, which scales MSTAR*USTAR^3 to + !! obtain energy available for mixing below + !! mixed layer base (nondimensional) + logical :: ML_rad_TKE_decay !< If true, apply same exponential decay + !! to ML_rad as applied to the other surface + !! sources of TKE in the mixed layer code. + real :: ustar_min !< A minimum value of ustar to avoid numerical + !! problems (m/s). If the value is small enough, + !! this parameter should not affect the solution. + real :: TKE_decay !< ratio of natural Ekman depth to TKE decay scale (nondim) + real :: mstar !! ratio of friction velocity cubed to + !! TKE input to the mixed layer (nondim) + logical :: ML_use_omega !< If true, use absolute rotation rate instead + !! of the vertical component of rotation when + !! setting the decay scale for mixed layer turbulence. + real :: ML_omega_frac !< When setting the decay scale for turbulence, use + !! this fraction of the absolute rotation rate blended + !! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2. + logical :: user_change_diff !< If true, call user-defined code to change diffusivity. + logical :: useKappaShear !< If true, use the kappa_shear module to find the + !! shear-driven diapycnal diffusivity. + logical :: use_CVMix_shear !< If true, use one of the CVMix modules to find + !! shear-driven diapycnal diffusivity. + logical :: use_CVMix_ddiff !< If true, enable double-diffusive mixing via CVMix. + logical :: simple_TKE_to_Kd !< If true, uses a simple estimate of Kd/TKE that + !! does not rely on a layer-formulation. character(len=200) :: inputdir type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() type(Kappa_shear_CS), pointer :: kappaShear_CSp => NULL() type(CVMix_shear_cs), pointer :: CVMix_shear_csp => NULL() + type(CVMix_ddiff_cs), pointer :: CVMix_ddiff_csp => NULL() type(bkgnd_mixing_cs), pointer :: bkgnd_mixing_csp => NULL() type(int_tide_CS), pointer :: int_tide_CSp => NULL() type(tidal_mixing_cs), pointer :: tm_csp => NULL() @@ -158,11 +157,6 @@ module MOM_set_diffusivity integer :: id_N2 = -1 integer :: id_N2_z = -1 - integer :: id_KT_extra = -1 - integer :: id_KS_extra = -1 - integer :: id_KT_extra_z = -1 - integer :: id_KS_extra_z = -1 - end type set_diffusivity_CS type diffusivity_diags @@ -172,12 +166,9 @@ module MOM_set_diffusivity Kd_BBL => NULL(),& ! BBL diffusivity at interfaces (m2/s) Kd_work => NULL(),& ! layer integrated work by diapycnal mixing (W/m2) maxTKE => NULL(),& ! energy required to entrain to h_max (m3/s3) - TKE_to_Kd => NULL(),& ! conversion rate (~1.0 / (G_Earth + dRho_lay)) + TKE_to_Kd => NULL() ! conversion rate (~1.0 / (G_Earth + dRho_lay)) ! between TKE dissipated within a layer and Kd ! in that layer, in m2 s-1 / m3 s-3 = s2 m-1 - KT_extra => NULL(),& ! double diffusion diffusivity for temp (m2/s) - KS_extra => NULL() ! double diffusion diffusivity for saln (m2/s) - end type diffusivity_diags ! Clocks @@ -185,6 +176,17 @@ module MOM_set_diffusivity contains +!> Sets the interior vertical diffusion of scalars due to the following processes: +!! 1) Shear-driven mixing: two options, Jackson et at. and KPP interior; +!! 2) Background mixing via CVMix (Bryan-Lewis profile) or the scheme described by +!! Harrison & Hallberg, JPO 2008; +!! 3) Double-diffusion aplpied via CVMix; +!! 4) Tidal mixing: many options available, see MOM_tidal_mixing.F90; +!! In addition, this subroutine has the option to set the interior vertical +!! viscosity associated with processes 2-4 listed above, which is stored in +!! visc%Kv_slow. Vertical viscosity due to shear-driven mixing is passed via +!! visc%Kv_shear +!! GMM, TODO: add contribution from tidal mixing into visc%Kv_slow subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & G, GV, CS, Kd, Kd_int) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -196,9 +198,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: u_h + intent(in) :: u_h !< zonal thickness transport m^2/s. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: v_h + intent(in) :: v_h !< meridional thickness transport m^2/s. type(thermo_var_ptrs), intent(inout) :: tv !< Structure with pointers to thermodynamic !! fields. Out is for tv%TempxPmE. type(forcing), intent(in) :: fluxes !< Structure of surface fluxes that may be @@ -226,17 +228,15 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! massless layers filled vertically by diffusion. real, dimension(SZI_(G),SZK_(G)) :: & - N2_lay, & ! squared buoyancy frequency associated with layers (1/s2) - maxTKE, & ! energy required to entrain to h_max (m3/s3) - TKE_to_Kd ! conversion rate (~1.0 / (G_Earth + dRho_lay)) between - ! TKE dissipated within a layer and Kd in that layer, in - ! m2 s-1 / m3 s-3 = s2 m-1. + N2_lay, & !< squared buoyancy frequency associated with layers (1/s2) + maxTKE, & !< energy required to entrain to h_max (m3/s3) + TKE_to_Kd !< conversion rate (~1.0 / (G_Earth + dRho_lay)) between + !< TKE dissipated within a layer and Kd in that layer, in + !< m2 s-1 / m3 s-3 = s2 m-1. real, dimension(SZI_(G),SZK_(G)+1) :: & - N2_int, & ! squared buoyancy frequency associated at interfaces (1/s2) - dRho_int, & ! locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? - KT_extra, & ! double difusion diffusivity on temperature (m2/sec) - KS_extra ! double difusion diffusivity on salinity (m2/sec) + N2_int, & !< squared buoyancy frequency associated at interfaces (1/s2) + dRho_int !< locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? real :: I_Rho0 ! inverse of Boussinesq density (m3/kg) real :: dissip ! local variable for dissipation calculations (W/m3) @@ -271,10 +271,16 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & use_EOS = associated(tv%eqn_of_state) - if ((CS%double_diffusion) .and. & + if ((CS%use_CVMix_ddiff) .and. & .not.(associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S)) ) & call MOM_error(FATAL, "set_diffusivity: visc%Kd_extra_T and "//& - "visc%Kd_extra_S must be associated when DOUBLE_DIFFUSION is true.") + "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF is true.") + + ! Set Kd, Kd_int and Kv_slow to constant values. + ! If nothing else is specified, this will be the value used. + Kd(:,:,:) = CS%Kd + Kd_int(:,:,:) = CS%Kd + if (associated(visc%Kv_slow)) visc%Kv_slow(:,:,:) = CS%Kv ! Set up arrays for diagnostics. @@ -293,12 +299,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%id_TKE_to_Kd > 0) then allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz)) ; dd%TKE_to_Kd(:,:,:) = 0.0 endif - if ((CS%id_KT_extra > 0) .or. (CS%id_KT_extra_z > 0)) then - allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1)) ; dd%KT_extra(:,:,:) = 0.0 - endif - if ((CS%id_KS_extra > 0) .or. (CS%id_KS_extra_z > 0)) then - allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1)) ; dd%KS_extra(:,:,:) = 0.0 - endif if ((CS%id_Kd_BBL > 0) .or. (CS%id_Kd_BBL_z > 0)) then allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0 endif @@ -341,6 +341,10 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & elseif (CS%use_CVMix_shear) then !NOTE{BGR}: this needs to be cleaned up. It works in 1D case, but has not been tested outside. call calculate_CVMix_shear(u_h, v_h, h, tv, visc%Kd_shear, visc%Kv_shear,G,GV,CS%CVMix_shear_csp) + if (CS%debug) then + call hchksum(visc%Kd_shear, "after CVMix_shear visc%Kd_shear",G%HI) + call hchksum(visc%Kv_shear, "after CVMix_shear visc%Kv_shear",G%HI) + endif elseif (associated(visc%Kv_shear)) then visc%Kv_shear(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled endif @@ -352,8 +356,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! set surface diffusivities (CS%bkgnd_mixing_csp%Kd_sfc) call sfc_bkgnd_mixing(G, CS%bkgnd_mixing_csp) -! GMM, fix OMP calls below - !$OMP parallel do default(none) shared(is,ie,js,je,nz,G,GV,CS,h,tv,T_f,S_f,fluxes,dd, & !$OMP Kd,visc, & !$OMP Kd_int,dt,u,v,Omega2) & @@ -370,35 +372,13 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo endif - ! add background mixing + ! Add background mixing call calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc%Kv_slow, j, G, GV, CS%bkgnd_mixing_csp) - ! GMM, the following will go into the MOM_CVMix_double_diffusion module - if (CS%double_diffusion) then - call double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, KT_extra, KS_extra) - do K=2,nz ; do i=is,ie - if (KS_extra(i,K) > KT_extra(i,K)) then ! salt fingering - Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KT_extra(i,K) - Kd(i,j,k) = Kd(i,j,k) + 0.5*KT_extra(i,K) - visc%Kd_extra_S(i,j,k) = KS_extra(i,K)-KT_extra(i,K) - visc%Kd_extra_T(i,j,k) = 0.0 - elseif (KT_extra(i,K) > 0.0) then ! double-diffusive convection - Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KS_extra(i,K) - Kd(i,j,k) = Kd(i,j,k) + 0.5*KS_extra(i,K) - visc%Kd_extra_T(i,j,k) = KT_extra(i,K)-KS_extra(i,K) - visc%Kd_extra_S(i,j,k) = 0.0 - else ! There is no double diffusion at this interface. - visc%Kd_extra_T(i,j,k) = 0.0 - visc%Kd_extra_S(i,j,k) = 0.0 - endif - enddo ; enddo - if (associated(dd%KT_extra)) then ; do K=1,nz+1 ; do i=is,ie - dd%KT_extra(i,j,K) = KT_extra(i,K) - enddo ; enddo ; endif - - if (associated(dd%KS_extra)) then ; do K=1,nz+1 ; do i=is,ie - dd%KS_extra(i,j,K) = KS_extra(i,K) - enddo ; enddo ; endif + ! Apply double diffusion + ! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available. + if (CS%use_CVMix_ddiff) then + call compute_ddiff_coeffs(h, tv, G, GV, j, visc%Kd_extra_T, visc%Kd_extra_S, CS%CVMix_ddiff_csp) endif ! Add the input turbulent diffusivity. @@ -496,6 +476,11 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%useKappaShear) call hchksum(visc%Kd_shear,"Turbulent Kd",G%HI,haloshift=0) + if (CS%use_CVMix_ddiff) then + call hchksum(visc%Kd_extra_T, "MOM_set_diffusivity: Kd_extra_T",G%HI,haloshift=0) + call hchksum(visc%Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S",G%HI,haloshift=0) + endif + if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, & G%HI, 0, symmetric=.true.) @@ -512,12 +497,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif - ! send bkgnd_mixing diagnostics to post_data - if (CS%bkgnd_mixing_csp%id_kd_bkgnd > 0) & - call post_data(CS%bkgnd_mixing_csp%id_kd_bkgnd, CS%bkgnd_mixing_csp%kd_bkgnd, CS%bkgnd_mixing_csp%diag) - if (CS%bkgnd_mixing_csp%id_kv_bkgnd > 0) & - call post_data(CS%bkgnd_mixing_csp%id_kv_bkgnd, CS%bkgnd_mixing_csp%kv_bkgnd, CS%bkgnd_mixing_csp%diag) - if (CS%Kd_add > 0.0) then if (present(Kd_int)) then !$OMP parallel do default(none) shared(is,ie,js,je,nz,Kd_int,CS,Kd) @@ -538,13 +517,28 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & T_f, S_f, dd%Kd_user) endif - ! GMM, post diags... - if (CS%id_Kd_layer > 0) call post_data(CS%id_Kd_layer, Kd, CS%diag) + ! post diagnostics - num_z_diags = 0 + ! background mixing + if (CS%bkgnd_mixing_csp%id_kd_bkgnd > 0) & + call post_data(CS%bkgnd_mixing_csp%id_kd_bkgnd, CS%bkgnd_mixing_csp%kd_bkgnd, CS%bkgnd_mixing_csp%diag) + if (CS%bkgnd_mixing_csp%id_kv_bkgnd > 0) & + call post_data(CS%bkgnd_mixing_csp%id_kv_bkgnd, CS%bkgnd_mixing_csp%kv_bkgnd, CS%bkgnd_mixing_csp%diag) + ! double diffusive mixing + if (CS%CVMix_ddiff_csp%id_KT_extra > 0) & + call post_data(CS%CVMix_ddiff_csp%id_KT_extra, visc%Kd_extra_T, CS%CVMix_ddiff_csp%diag) + if (CS%CVMix_ddiff_csp%id_KS_extra > 0) & + call post_data(CS%CVMix_ddiff_csp%id_KS_extra, visc%Kd_extra_S, CS%CVMix_ddiff_csp%diag) + if (CS%CVMix_ddiff_csp%id_R_rho > 0) & + call post_data(CS%CVMix_ddiff_csp%id_R_rho, CS%CVMix_ddiff_csp%R_rho, CS%CVMix_ddiff_csp%diag) + + if (CS%id_Kd_layer > 0) call post_data(CS%id_Kd_layer, Kd, CS%diag) + + ! tidal mixing call post_tidal_diagnostics(G,GV,h,CS%tm_csp) + num_z_diags = 0 if (CS%tm_csp%Int_tide_dissipation .or. CS%tm_csp%Lee_wave_dissipation .or. & CS%tm_csp%Lowmode_itidal_dissipation) then @@ -568,26 +562,11 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif - if (CS%id_KT_extra > 0) call post_data(CS%id_KT_extra, dd%KT_extra, CS%diag) - if (CS%id_KS_extra > 0) call post_data(CS%id_KS_extra, dd%KS_extra, CS%diag) if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, dd%Kd_BBL, CS%diag) - if (CS%id_KT_extra_z > 0) then - num_z_diags = num_z_diags + 1 - z_ids(num_z_diags) = CS%id_KT_extra_z - z_ptrs(num_z_diags)%p => dd%KT_extra - endif - - if (CS%id_KS_extra_z > 0) then - num_z_diags = num_z_diags + 1 - z_ids(num_z_diags) = CS%id_KS_extra_z - z_ptrs(num_z_diags)%p => dd%KS_extra - endif - if (CS%id_Kd_BBL_z > 0) then num_z_diags = num_z_diags + 1 z_ids(num_z_diags) = CS%id_Kd_BBL_z - z_ptrs(num_z_diags)%p => dd%KS_extra endif if (num_z_diags > 0) & @@ -598,8 +577,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (associated(dd%Kd_user)) deallocate(dd%Kd_user) if (associated(dd%maxTKE)) deallocate(dd%maxTKE) if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd) - if (associated(dd%KT_extra)) deallocate(dd%KT_extra) - if (associated(dd%KS_extra)) deallocate(dd%KS_extra) if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL) if (showCallTree) call callTree_leave("set_diffusivity()") @@ -952,119 +929,6 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & end subroutine find_N2 -! GMM, the following will be moved to a new module - -!> This subroutine sets the additional diffusivities of temperature and -!! salinity due to double diffusion, using the same functional form as is -!! used in MOM4.1, and taken from an NCAR technical note (REF?) that updates -!! what was in Large et al. (1994). All the coefficients here should probably -!! be made run-time variables rather than hard-coded constants. -!! -!! \todo Find reference for NCAR tech note above. -subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, Kd_T_dd, Kd_S_dd) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available - !! thermodynamic fields; absent fields have NULL - !! ptrs. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: T_f !< layer temp in C with the values in massless layers - !! filled vertically by diffusion. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: S_f !< Layer salinities in PPT with values in massless - !! layers filled vertically by diffusion. - integer, intent(in) :: j !< Meridional index upon which to work. - type(set_diffusivity_CS), pointer :: CS !< Module control structure. - real, dimension(SZI_(G),SZK_(G)+1), & - intent(out) :: Kd_T_dd !< Interface double diffusion diapycnal - !! diffusivity for temp (m2/sec). - real, dimension(SZI_(G),SZK_(G)+1), & - intent(out) :: Kd_S_dd !< Interface double diffusion diapycnal - !! diffusivity for saln (m2/sec). - -! Arguments: -! (in) tv - structure containing pointers to any available -! thermodynamic fields; absent fields have NULL ptrs -! (in) h - layer thickness (m or kg m-2) -! (in) T_f - layer temp in C with the values in massless layers -! filled vertically by diffusion -! (in) S_f - layer salinities in PPT with values in massless layers -! filled vertically by diffusion -! (in) G - ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) CS - module control structure -! (in) j - meridional index upon which to work -! (out) Kd_T_dd - interface double diffusion diapycnal diffusivity for temp (m2/sec) -! (out) Kd_S_dd - interface double diffusion diapycnal diffusivity for saln (m2/sec) - -! This subroutine sets the additional diffusivities of temperature and -! salinity due to double diffusion, using the same functional form as is -! used in MOM4.1, and taken from an NCAR technical note (###REF?) that updates -! what was in Large et al. (1994). All the coefficients here should probably -! be made run-time variables rather than hard-coded constants. - - real, dimension(SZI_(G)) :: & - dRho_dT, & ! partial derivatives of density wrt temp (kg m-3 degC-1) - dRho_dS, & ! partial derivatives of density wrt saln (kg m-3 PPT-1) - pres, & ! pressure at each interface (Pa) - Temp_int, & ! temp and saln at interfaces - Salin_int - - real :: alpha_dT ! density difference between layers due to temp diffs (kg/m3) - real :: beta_dS ! density difference between layers due to saln diffs (kg/m3) - - real :: Rrho ! vertical density ratio - real :: diff_dd ! factor for double-diffusion - real :: prandtl ! flux ratio for diffusive convection regime - - real, parameter :: Rrho0 = 1.9 ! limit for double-diffusive density ratio - real, parameter :: dsfmax = 1.e-4 ! max diffusivity in case of salt fingering - real, parameter :: Kv_molecular = 1.5e-6 ! molecular viscosity (m2/sec) - - integer :: i, k, is, ie, nz - is = G%isc ; ie = G%iec ; nz = G%ke - - if (associated(tv%eqn_of_state)) then - do i=is,ie - pres(i) = 0.0 ; Kd_T_dd(i,1) = 0.0 ; Kd_S_dd(i,1) = 0.0 - Kd_T_dd(i,nz+1) = 0.0 ; Kd_S_dd(i,nz+1) = 0.0 - enddo - do K=2,nz - do i=is,ie - pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1) - Temp_Int(i) = 0.5 * (T_f(i,j,k-1) + T_f(i,j,k)) - Salin_Int(i) = 0.5 * (S_f(i,j,k-1) + S_f(i,j,k)) - enddo - call calculate_density_derivs(Temp_int, Salin_int, pres, & - dRho_dT(:), dRho_dS(:), is, ie-is+1, tv%eqn_of_state) - - do i=is,ie - alpha_dT = -1.0*dRho_dT(i) * (T_f(i,j,k-1) - T_f(i,j,k)) - beta_dS = dRho_dS(i) * (S_f(i,j,k-1) - S_f(i,j,k)) - - if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then ! salt finger case - Rrho = min(alpha_dT/beta_dS,Rrho0) - diff_dd = 1.0 - ((RRho-1.0)/(RRho0-1.0)) - diff_dd = dsfmax*diff_dd*diff_dd*diff_dd - Kd_T_dd(i,K) = 0.7*diff_dd - Kd_S_dd(i,K) = diff_dd - elseif ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then ! diffusive convection - Rrho = alpha_dT/beta_dS - diff_dd = Kv_molecular*0.909*exp(4.6*exp(-0.54*(1/Rrho-1))) - prandtl = 0.15*Rrho - if (Rrho > 0.5) prandtl = (1.85-0.85/Rrho)*Rrho - Kd_T_dd(i,K) = diff_dd - Kd_S_dd(i,K) = prandtl*diff_dd - else - Kd_T_dd(i,K) = 0.0 ; Kd_S_dd(i,K) = 0.0 - endif - enddo - enddo - endif - -end subroutine double_diffusion !> This routine adds diffusion sustained by flow energy extracted by bottom drag. subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & maxTKE, kb, G, GV, CS, Kd, Kd_int, Kd_BBL) @@ -1974,6 +1838,11 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp ! set params releted to the background mixing call bkgnd_mixing_init(Time, G, GV, param_file, CS%diag, CS%bkgnd_mixing_csp) + call get_param(param_file, mdl, "KV", CS%Kv, & + "The background kinematic viscosity in the interior. \n"//& + "The molecular value, ~1e-6 m2 s-1, may be used.", & + units="m2 s-1", fail_if_missing=.true.) + call get_param(param_file, mdl, "KD", CS%Kd, & "The background diapycnal diffusivity of density in the \n"//& "interior. Zero or the molecular value, ~1e-7 m2 s-1, \n"//& @@ -2076,45 +1945,6 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp endif endif - - ! GMM, the following should be moved to the DD module - call get_param(param_file, mdl, "DOUBLE_DIFFUSION", CS%double_diffusion, & - "If true, increase diffusivitives for temperature or salt \n"//& - "based on double-diffusive paramaterization from MOM4/KPP.", & - default=.false.) - if (CS%double_diffusion) then - call get_param(param_file, mdl, "MAX_RRHO_SALT_FINGERS", CS%Max_Rrho_salt_fingers, & - "Maximum density ratio for salt fingering regime.", & - default=2.55, units="nondim") - call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", CS%Max_salt_diff_salt_fingers, & - "Maximum salt diffusivity for salt fingering regime.", & - default=1.e-4, units="m2 s-1") - call get_param(param_file, mdl, "KV_MOLECULAR", CS%Kv_molecular, & - "Molecular viscosity for calculation of fluxes under \n"//& - "double-diffusive convection.", default=1.5e-6, units="m2 s-1") - ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults. - - CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for temperature', 'm2 s-1') - - CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for salinity', 'm2 s-1') - - if (associated(diag_to_Z_CSp)) then - vd = var_desc("KT_extra", "m2 s-1", & - "Double-Diffusive Temperature Diffusivity, interpolated to z", & - z_grid='z') - CS%id_KT_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("KS_extra", "m2 s-1", & - "Double-Diffusive Salinity Diffusivity, interpolated to z",& - z_grid='z') - CS%id_KS_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("Kd_BBL", "m2 s-1", & - "Bottom Boundary Layer Diffusivity", z_grid='z') - CS%id_Kd_BBL_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - endif - endif - if (CS%user_change_diff) then call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) endif @@ -2130,6 +1960,9 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp ! CVMix shear-driven mixing CS%use_CVMix_shear = CVMix_shear_init(Time, G, GV, param_file, CS%diag, CS%CVMix_shear_csp) + ! CVMix double diffusion mixing + CS%use_CVMix_ddiff = CVMix_ddiff_init(Time, G, GV, param_file, CS%diag, CS%CVMix_ddiff_csp) + end subroutine set_diffusivity_init !> Clear pointers and dealocate memory @@ -2146,6 +1979,9 @@ subroutine set_diffusivity_end(CS) if (CS%use_CVMix_shear) & call CVMix_shear_end(CS%CVMix_shear_csp) + if (CS%use_CVMix_ddiff) & + call CVMix_ddiff_end(CS%CVMix_ddiff_csp) + if (associated(CS)) deallocate(CS) end subroutine set_diffusivity_end diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index ec0b5a80b3..ec1b09a5ad 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -2,38 +2,6 @@ module MOM_set_visc ! This file is part of MOM6. See LICENSE.md for the license. -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, April 1994 - October 2006 * -!* Quadratic Bottom Drag by James Stephens and R. Hallberg. * -!* * -!* This file contains the subroutine that calculates various values * -!* related to the bottom boundary layer, such as the viscosity and * -!* thickness of the BBL (set_viscous_BBL). This would also be the * -!* module in which other viscous quantities that are flow-independent * -!* might be set. This information is transmitted to other modules * -!* via a vertvisc type structure. * -!* * -!* The same code is used for the two velocity components, by * -!* indirectly referencing the velocities and defining a handful of * -!* direction-specific defined variables. * -!* * -!* Macros written all in capital letters are defined in MOM_memory.h. * -!* * -!* A small fragment of the grid is shown below: * -!* * -!* j+1 x ^ x ^ x At x: q * -!* j+1 > o > o > At ^: v, frhatv, tauy * -!* j x ^ x ^ x At >: u, frhatu, taux * -!* j > o > o > At o: h * -!* j-1 x ^ x ^ x * -!* i-1 i i+1 At x & ^: * -!* i i+1 At > & o: * -!* * -!* The boundaries always run through q grid points (x). * -!* * -!********+*********+*********+*********+*********+*********+*********+** - use MOM_debugging, only : uvchksum, hchksum use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr @@ -44,8 +12,9 @@ module MOM_set_visc use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_kappa_shear, only : kappa_shear_is_used -use MOM_CVMix_shear, only : CVMix_shear_is_used -use MOM_CVMix_conv, only : CVMix_conv_is_used +use MOM_cvmix_shear, only : cvmix_shear_is_used +use MOM_cvmix_conv, only : cvmix_conv_is_used +use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used use MOM_io, only : vardesc, var_desc use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_variables, only : thermo_var_ptrs @@ -1791,8 +1760,10 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., & do_not_log=.true.) + use_kappa_shear = .false. ; use_CVMix_shear = .false. useKPP = .false. ; useEPBL = .false. ; use_CVMix_conv = .false. + if (.not.adiabatic) then use_kappa_shear = kappa_shear_is_used(param_file) use_CVMix_shear = CVMix_shear_is_used(param_file) @@ -1811,7 +1782,9 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) allocate(visc%Kd_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kd_shear(:,:,:) = 0.0 allocate(visc%TKE_turb(isd:ied,jsd:jed,nz+1)) ; visc%TKE_turb(:,:,:) = 0.0 allocate(visc%Kv_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kv_shear(:,:,:) = 0.0 - allocate(visc%Kv_slow(isd:ied,jsd:jed,nz+1)) ; visc%Kv_slow(:,:,:) = 0.0 + + ! MOM_bkgnd_mixing is always used, so always allocate visc%Kv_slow. GMM + allocate(visc%Kv_slow(isd:ied,jsd:jed,nz+1)) ; visc%Kv_slow(:,:,:) = 0.0 vd = var_desc("Kd_shear","m2 s-1","Shear-driven turbulent diffusivity at interfaces", & hor_grid='h', z_grid='i') @@ -1854,21 +1827,14 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) type(set_visc_CS), pointer :: CS !< A pointer that is set to point to the control !! structure for this module type(ocean_OBC_type), pointer :: OBC -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (out) visc - A structure containing vertical viscosities and related -! fields. Allocated here. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module + + ! local variables real :: Csmag_chan_dflt, smag_const1, TKE_decay_dflt, bulk_Ri_ML_dflt real :: Kv_background real :: omega_frac_dflt integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, i, j, n - logical :: use_kappa_shear, adiabatic, differential_diffusion, use_omega + logical :: use_kappa_shear, adiabatic, use_omega + logical :: use_CVMix_ddiff type(OBC_segment_type), pointer :: segment ! pointer to OBC segment type ! This include declares and sets the variable "version". #include "version_variable.h" @@ -1891,8 +1857,8 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) ! Set default, read and log parameters call log_version(param_file, mdl, version, "") - CS%RiNo_mix = .false. - use_kappa_shear = .false. ; differential_diffusion = .false. !; adiabatic = .false. ! Needed? -AJA + CS%RiNo_mix = .false. ; use_CVMix_ddiff = .false. + use_kappa_shear = .false. !; adiabatic = .false. ! Needed? -AJA call get_param(param_file, mdl, "BOTTOMDRAGLAW", CS%bottomdraglaw, & "If true, the bottom stress is calculated with a drag \n"//& "law of the form c_drag*|u|*u. The velocity magnitude \n"//& @@ -1919,11 +1885,9 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) if (.not.adiabatic) then use_kappa_shear = kappa_shear_is_used(param_file) CS%RiNo_mix = use_kappa_shear - call get_param(param_file, mdl, "DOUBLE_DIFFUSION", differential_diffusion, & - "If true, increase diffusivitives for temperature or salt \n"//& - "based on double-diffusive paramaterization from MOM4/KPP.", & - default=.false.) + use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) endif + call get_param(param_file, mdl, "PRANDTL_TURB", visc%Prandtl_turb, & "The turbulent Prandtl number applied to shear \n"//& "instability.", units="nondim", default=1.0) @@ -2016,6 +1980,15 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) "The background kinematic viscosity in the interior. \n"//& "The molecular value, ~1e-6 m2 s-1, may be used.", & units="m2 s-1", fail_if_missing=.true.) + + call get_param(param_file, mdl, "ADD_KV_SLOW", visc%add_Kv_slow, & + "If true, the background vertical viscosity in the interior \n"//& + "(i.e., tidal + background + shear + convenction) is addded \n"// & + "when computing the coupling coefficient. The purpose of this \n"// & + "flag is to be able to recover previous answers and it will likely \n"// & + "be removed in the future since this option should always be true.", & + default=.false.) + call get_param(param_file, mdl, "KV_BBL_MIN", CS%KV_BBL_min, & "The minimum viscosities in the bottom boundary layer.", & units="m2 s-1", default=Kv_background) @@ -2065,7 +2038,7 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) Time, 'Rayleigh drag velocity at v points', 'm s-1') endif - if (differential_diffusion) then + if (use_CVMix_ddiff) then allocate(visc%Kd_extra_T(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_T = 0.0 allocate(visc%Kd_extra_S(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_S = 0.0 endif @@ -2113,4 +2086,37 @@ subroutine set_visc_end(visc, CS) deallocate(CS) end subroutine set_visc_end +!> \namespace MOM_set_visc +!!********+*********+*********+*********+*********+*********+*********+** +!!* * +!!* By Robert Hallberg, April 1994 - October 2006 * +!!* Quadratic Bottom Drag by James Stephens and R. Hallberg. * +!!* * +!!* This file contains the subroutine that calculates various values * +!!* related to the bottom boundary layer, such as the viscosity and * +!!* thickness of the BBL (set_viscous_BBL). This would also be the * +!!* module in which other viscous quantities that are flow-independent * +!!* might be set. This information is transmitted to other modules * +!!* via a vertvisc type structure. * +!!* * +!!* The same code is used for the two velocity components, by * +!!* indirectly referencing the velocities and defining a handful of * +!!* direction-specific defined variables. * +!!* * +!!* Macros written all in capital letters are defined in MOM_memory.h. * +!!* * +!!* A small fragment of the grid is shown below: * +!!* * +!!* j+1 x ^ x ^ x At x: q * +!!* j+1 > o > o > At ^: v, frhatv, tauy * +!!* j x ^ x ^ x At >: u, frhatu, taux * +!!* j > o > o > At o: h * +!!* j-1 x ^ x ^ x * +!!* i-1 i i+1 At x & ^: * +!!* i i+1 At > & o: * +!!* * +!!* The boundaries always run through q grid points (x). * +!!* * +!!********+*********+*********+*********+*********+*********+*********+** + end module MOM_set_visc diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 1bb8eb48dd..9c60a05451 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -3,24 +3,27 @@ module MOM_tidal_mixing ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field -use MOM_diag_mediator, only : safe_alloc_ptr, post_data -use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag -use MOM_diag_to_Z, only : calc_Zint_diags -use MOM_EOS, only : calculate_density -use MOM_variables, only : thermo_var_ptrs, p3d -use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE -use MOM_debugging, only : hchksum -use MOM_grid, only : ocean_grid_type -use MOM_verticalGrid, only : verticalGrid_type -use MOM_file_parser, only : openParameterBlock, closeParameterBlock -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_string_functions, only : uppercase, lowercase -use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc -use CVMix_tidal, only : CVMix_init_tidal, CVMix_compute_Simmons_invariant -use CVMix_tidal, only : CVMix_coeffs_tidal, CVMix_tidal_params_type -use CVMix_kinds_and_types, only : CVMix_global_params_type -use CVMix_put_get, only : CVMix_put +use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field +use MOM_diag_mediator, only : safe_alloc_ptr, post_data +use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag +use MOM_diag_to_Z, only : calc_Zint_diags +use MOM_EOS, only : calculate_density +use MOM_variables, only : thermo_var_ptrs, p3d +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_debugging, only : hchksum +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_string_functions, only : uppercase, lowercase +use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc, field_size +use CVMix_tidal, only : CVMix_init_tidal, CVMix_compute_Simmons_invariant +use CVMix_tidal, only : CVMix_coeffs_tidal, CVMix_tidal_params_type +use CVMix_tidal, only : CVMix_compute_Schmittner_invariant, CVMix_compute_SchmittnerCoeff +use CVMix_tidal, only : CVMix_coeffs_tidal_schmittner +use CVMix_kinds_and_types, only : CVMix_global_params_type +use CVMix_put_get, only : CVMix_put implicit none ; private @@ -34,19 +37,22 @@ module MOM_tidal_mixing !> Containers for tidal mixing diagnostics type, public :: tidal_mixing_diags + private real, pointer, dimension(:,:,:) :: & - Kd_itidal => NULL(),& ! internal tide diffusivity at interfaces (m2/s) - Fl_itidal => NULL(),& ! vertical flux of tidal turbulent dissipation (m3/s3) - Kd_lowmode => NULL(),& ! internal tide diffusivity at interfaces - ! due to propagating low modes (m2/s) (BDM) - Fl_lowmode => NULL(),& ! vertical flux of tidal turbulent dissipation - ! due to propagating low modes (m3/s3) (BDM) - Kd_Niku => NULL(),& ! lee-wave diffusivity at interfaces (m2/s) - Kd_Niku_work => NULL(),& ! layer integrated work by lee-wave driven mixing (W/m2) - Kd_Itidal_Work => NULL(),& ! layer integrated work by int tide driven mixing (W/m2) - Kd_Lowmode_Work => NULL(),& ! layer integrated work by low mode driven mixing (W/m2) BDM - N2_int => NULL(),& - vert_dep_3d => NULL() + Kd_itidal => NULL(),& ! internal tide diffusivity at interfaces (m2/s) + Fl_itidal => NULL(),& ! vertical flux of tidal turbulent dissipation (m3/s3) + Kd_lowmode => NULL(),& ! internal tide diffusivity at interfaces + ! due to propagating low modes (m2/s) (BDM) + Fl_lowmode => NULL(),& ! vertical flux of tidal turbulent dissipation + ! due to propagating low modes (m3/s3) (BDM) + Kd_Niku => NULL(),& ! lee-wave diffusivity at interfaces (m2/s) + Kd_Niku_work => NULL(),& ! layer integrated work by lee-wave driven mixing (W/m2) + Kd_Itidal_Work => NULL(),& ! layer integrated work by int tide driven mixing (W/m2) + Kd_Lowmode_Work => NULL(),& ! layer integrated work by low mode driven mixing (W/m2) BDM + N2_int => NULL(),& + vert_dep_3d => NULL(),& + Schmittner_coeff_3d => NULL(),& + tidal_qe_md => NULL() real, pointer, dimension(:,:) :: & TKE_itidal_used => NULL(),& ! internal tide TKE input at ocean bottom (W/m2) @@ -61,6 +67,7 @@ module MOM_tidal_mixing !> Control structure for tidal mixing module. type, public :: tidal_mixing_cs logical :: debug = .true. + ! TODO: private ! Parameters logical :: int_tide_dissipation = .false. ! Internal tide conversion (from barotropic) @@ -128,18 +135,23 @@ module MOM_tidal_mixing real :: min_thickness ! Minimum thickness allowed [m] ! CVMix-specific parameters + integer :: CVMix_tidal_scheme = -1 ! 1 for Simmons, 2 for Schmittner type(CVMix_tidal_params_type) :: CVMix_tidal_params - type(CVMix_global_params_type) :: CVMix_glb_params ! for Prandtl number only - real :: tidal_max_coef ! maximum allowable tidal diffusivity. [m^2/s] + type(CVMix_global_params_type) :: CVMix_glb_params ! for Prandtl number only + real :: tidal_max_coef ! maximum allowable tidal diffusivity. [m^2/s] + real :: tidal_diss_lim_tc ! dissipation limit for tidal-energy-constituent data + type(remapping_CS) :: remap_cs ! Data containers - real, pointer, dimension(:,:) :: TKE_Niku => NULL() - real, pointer, dimension(:,:) :: TKE_itidal => NULL() - real, pointer, dimension(:,:) :: Nb => NULL() - real, pointer, dimension(:,:) :: mask_itidal => NULL() - real, pointer, dimension(:,:) :: h2 => NULL() - real, pointer, dimension(:,:) :: tideamp => NULL() ! RMS tidal amplitude (m/s) - real, allocatable,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) + real, pointer, dimension(:,:) :: TKE_Niku => NULL() + real, pointer, dimension(:,:) :: TKE_itidal => NULL() + real, pointer, dimension(:,:) :: Nb => NULL() + real, pointer, dimension(:,:) :: mask_itidal => NULL() + real, pointer, dimension(:,:) :: h2 => NULL() + real, pointer, dimension(:,:) :: tideamp => NULL() !< RMS tidal amplitude [m/s] + real, allocatable, dimension(:) :: h_src !< tidal constituent input layer thickness [m] + real, allocatable ,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) ! TODO: make this E(x,y) only + real, allocatable ,dimension(:,:,:) :: tidal_qe_3d_in ! q*E(x,y) ! Diagnostics type(diag_ctrl), pointer :: diag => NULL() ! structure to regulate diagn output timing @@ -167,6 +179,8 @@ module MOM_tidal_mixing integer :: id_Polzin_decay_scale_scaled = -1 integer :: id_N2_int = -1 integer :: id_Simmons_coeff = -1 + integer :: id_Schmittner_coeff = -1 + integer :: id_tidal_qe_md = -1 integer :: id_vert_dep = -1 end type tidal_mixing_cs @@ -174,12 +188,12 @@ module MOM_tidal_mixing character(len=40) :: mdl = "MOM_tidal_mixing" !< This module's name. character*(20), parameter :: STLAURENT_PROFILE_STRING = "STLAURENT_02" character*(20), parameter :: POLZIN_PROFILE_STRING = "POLZIN_09" -character*(20), parameter :: SIMMONS_PROFILE_STRING = "SIMMONS" -character*(20), parameter :: SCHMITTNER_PROFILE_STRING = "SCHMITTNER" integer, parameter :: STLAURENT_02 = 1 integer, parameter :: POLZIN_09 = 2 -integer, parameter :: SIMMONS_04 = 3 -integer, parameter :: SCHMITTNER = 4 +character*(20), parameter :: SIMMONS_SCHEME_STRING = "SIMMONS" +character*(20), parameter :: SCHMITTNER_SCHEME_STRING = "SCHMITTNER" +integer, parameter :: SIMMONS = 1 +integer, parameter :: SCHMITTNER = 2 contains @@ -197,7 +211,7 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, ! Local variables logical :: read_tideamp character(len=20) :: tmpstr, int_tide_profile_str - character(len=20) :: default_profile_string, tidal_energy_type + character(len=20) :: CVMix_tidal_scheme_str, tidal_energy_type character(len=200) :: filename, h2_file, Niku_TKE_input_file character(len=200) :: tidal_energy_file, tideamp_file type(vardesc) :: vd @@ -244,40 +258,47 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, if (.not. tidal_mixing_init) return if (CS%int_tide_dissipation) then - default_profile_string = STLAURENT_PROFILE_STRING - if (CS%use_CVMix_tidal) default_profile_string = SIMMONS_PROFILE_STRING - call get_param(param_file, mdl, "INT_TIDE_PROFILE", int_tide_profile_str, & - "INT_TIDE_PROFILE selects the vertical profile of energy \n"//& - "dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//& - "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//& - "\t decay profile.\n"//& - "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//& - "\t decay profile.", & - default=default_profile_string) - ! TODO: list the newly available profile selections - int_tide_profile_str = uppercase(int_tide_profile_str) - select case (int_tide_profile_str) - case (STLAURENT_PROFILE_STRING) ; CS%int_tide_profile = STLAURENT_02 - case (POLZIN_PROFILE_STRING) ; CS%int_tide_profile = POLZIN_09 - case (SIMMONS_PROFILE_STRING) ; CS%int_tide_profile = SIMMONS_04 - case (SCHMITTNER_PROFILE_STRING) ; CS%int_tide_profile = SCHMITTNER - case default - call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & - "#define INT_TIDE_PROFILE "//trim(int_tide_profile_str)//" found in input file.") - end select - ! Check profile consistency - if (CS%use_CVMix_tidal .and. (CS%int_tide_profile == STLAURENT_02 .or. & - CS%int_tide_profile == POLZIN_09)) then - call MOM_error(FATAL, "tidal_mixing_init: Tidal mixing profile"// & - " "//trim(int_tide_profile_str)//" unavailable in CVMix. Available "//& - "profiles in CVMix are "//trim(SIMMONS_PROFILE_STRING)//" and "//& - trim(SCHMITTNER_PROFILE_STRING)//".") - elseif (.not.CS%use_CVMix_tidal .and. (CS%int_tide_profile == SIMMONS_04.or. & - CS%int_tide_profile == SCHMITTNER)) then - call MOM_error(FATAL, "tidal_mixing_init: Tidal mixing profiles "// & - trim(SIMMONS_PROFILE_STRING)//" and "//trim(SCHMITTNER_PROFILE_STRING)//& - " are available only when USE_CVMix_TIDAL is True.") + ! Read in CVMix tidal scheme if CVMix tidal mixing is on + if (CS%use_CVMix_tidal) then + call get_param(param_file, mdl, "CVMIX_TIDAL_SCHEME", CVMix_tidal_scheme_str, & + "CVMIX_TIDAL_SCHEME selects the CVMix tidal mixing\n"//& + "scheme with INT_TIDE_DISSIPATION. Valid values are:\n"//& + "\t SIMMONS - Use the Simmons et al (2004) tidal \n"//& + "\t mixing scheme.\n"//& + "\t SCHMITTNER - Use the Schmittner et al (2014) tidal \n"//& + "\t mixing scheme.", & + default=SIMMONS_SCHEME_STRING) + CVMix_tidal_scheme_str = uppercase(CVMix_tidal_scheme_str) + + select case (CVMix_tidal_scheme_str) + case (SIMMONS_SCHEME_STRING) ; CS%CVMix_tidal_scheme = SIMMONS + case (SCHMITTNER_SCHEME_STRING) ; CS%CVMix_tidal_scheme = SCHMITTNER + case default + call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & + "#define CVMIX_TIDAL_SCHEME "//trim(CVMix_tidal_scheme_str)//" found in input file.") + end select + endif ! CS%use_CVMix_tidal + + ! Read in vertical profile of tidal energy dissipation + if ( CS%CVMix_tidal_scheme.eq.SCHMITTNER .or. .not. CS%use_CVMix_tidal) then + call get_param(param_file, mdl, "INT_TIDE_PROFILE", int_tide_profile_str, & + "INT_TIDE_PROFILE selects the vertical profile of energy \n"//& + "dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//& + "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//& + "\t decay profile.\n"//& + "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//& + "\t decay profile.", & + default=STLAURENT_PROFILE_STRING) + int_tide_profile_str = uppercase(int_tide_profile_str) + + select case (int_tide_profile_str) + case (STLAURENT_PROFILE_STRING) ; CS%int_tide_profile = STLAURENT_02 + case (POLZIN_PROFILE_STRING) ; CS%int_tide_profile = POLZIN_09 + case default + call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & + "#define INT_TIDE_PROFILE "//trim(int_tide_profile_str)//" found in input file.") + end select endif elseif (CS%use_CVMix_tidal) then @@ -480,15 +501,14 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, call get_param(param_file, mdl, "TIDAL_MAX_COEF", CS%tidal_max_coef, & "largest acceptable value for tidal diffusivity", & units="m^2/s", default=50e-4) ! the default is 50e-4 in CVMix, 100e-4 in POP. + call get_param(param_file, mdl, "TIDAL_DISS_LIM_TC", CS%tidal_diss_lim_tc, & + "Min allowable depth for dissipation for tidal-energy-constituent data. \n"//& + "No dissipation contribution is applied above TIDAL_DISS_LIM_TC.", & + units="m", default=0.0) call get_param(param_file, mdl, "TIDAL_ENERGY_FILE",tidal_energy_file, & "The path to the file containing tidal energy \n"//& "dissipation. Used with CVMix tidal mixing schemes.", & fail_if_missing=.true.) - tidal_energy_file = trim(CS%inputdir) // trim(tidal_energy_file) - call get_param(param_file, mdl, "TIDAL_ENERGY_TYPE",tidal_energy_type, & - "The type of input tidal energy flux dataset.",& - fail_if_missing=.true.) - ! TODO: list all available tidal energy types here call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, & do_not_log=.True.) call get_param(param_file, mdl, "PRANDTL_TIDAL", prandtl_tidal, & @@ -498,14 +518,25 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, do_not_log=.true.) call CVMix_put(CS%CVMix_glb_params,'Prandtl',prandtl_tidal) - int_tide_profile_str = lowercase(int_tide_profile_str) - - - ! TODO: check parameter consistency. (see POP::tidal_mixing.F90::tidal_check) + tidal_energy_file = trim(CS%inputdir) // trim(tidal_energy_file) + call get_param(param_file, mdl, "TIDAL_ENERGY_TYPE",tidal_energy_type, & + "The type of input tidal energy flux dataset. Valid values are"//& + "\t Jayne\n"//& + "\t ER03 \n",& + fail_if_missing=.true.) + ! Check whether tidal energy input format and CVMix tidal mixing scheme are consistent + if ( .not. ( & + (uppercase(tidal_energy_type(1:4)).eq.'JAYN' .and. CS%CVMix_tidal_scheme.eq.SIMMONS).or. & + (uppercase(tidal_energy_type(1:4)).eq.'ER03' .and. CS%CVMix_tidal_scheme.eq.SCHMITTNER) ) )then + call MOM_error(FATAL, "tidal_mixing_init: Tidal energy file type ("//& + trim(tidal_energy_type)//") is incompatible with CVMix tidal "//& + " mixing scheme: "//trim(CVMix_tidal_scheme_str) ) + endif + CVMix_tidal_scheme_str = lowercase(CVMix_tidal_scheme_str) ! Set up CVMix - call CVMix_init_tidal(CVMix_tidal_params_user = CS%CVMix_tidal_params, & - mix_scheme = int_tide_profile_str, & + call CVMix_init_tidal(CVmix_tidal_params_user = CS%CVMix_tidal_params, & + mix_scheme = CVMix_tidal_scheme_str, & efficiency = CS%Mu_itides, & vertical_decay_scale = CS%int_tide_decay_scale, & max_coefficient = CS%tidal_max_coef, & @@ -532,6 +563,10 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, ! TODO: add units CS%id_Simmons_coeff = register_diag_field('ocean_model','Simmons_coeff',diag%axesT1,Time, & 'time-invariant portion of the tidal mixing coefficient using the Simmons', '') + CS%id_Schmittner_coeff = register_diag_field('ocean_model','Schmittner_coeff',diag%axesTL,Time, & + 'time-invariant portion of the tidal mixing coefficient using the Schmittner', '') + CS%id_tidal_qe_md = register_diag_field('ocean_model','tidal_qe_md',diag%axesTL,Time, & + 'input tidal energy dissipated locally interpolated to model vertical coordinates', '') CS%id_vert_dep = register_diag_field('ocean_model','vert_dep',diag%axesTi,Time, & 'vertical deposition function needed for Simmons et al tidal mixing', '') @@ -646,20 +681,27 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) ! local real, dimension(SZK_(G)+1) :: Kd_tidal !< tidal diffusivity [m2/s] real, dimension(SZK_(G)+1) :: Kv_tidal !< tidal viscosity [m2/s] - real, dimension(SZK_(G)+1) :: vert_dep !< vertical deposition needed for Simmons tidal mixing. + real, dimension(SZK_(G)+1) :: vert_dep !< vertical deposition real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) + real, dimension(SZK_(G)+1) :: SchmittnerSocn real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) + real, dimension(SZK_(G)) :: tidal_qe_md !< Tidal dissipation energy interpolated from 3d input to model coordinates + real, dimension(SZK_(G)) :: Schmittner_coeff + real, dimension(SZK_(G)) :: h_m !< Cell thickness [m] + real, allocatable, dimension(:,:) :: exp_hab_zetar + integer :: i, k, is, ie real :: dh, hcorr, Simmons_coeff real, parameter :: rho_fw = 1000.0 ! fresh water density [kg/m^3] ! TODO: when coupled, get this from CESM (SHR_CONST_RHOFW) + real :: h_neglect, h_neglect_edge type(tidal_mixing_diags), pointer :: dd is = G%isc ; ie = G%iec dd => CS%dd - select case (CS%int_tide_profile) - case (SIMMONS_04) + select case (CS%CVMix_tidal_scheme) + case (SIMMONS) do i=is,ie if (G%mask2dT(i,j)<1) cycle @@ -723,10 +765,105 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) enddo ! i=is,ie - ! TODO: case (SCHMITTNER) + case (SCHMITTNER) + + ! TODO: correct exp_hab_zetar shapes in CVMix_compute_Schmittner_invariant + ! and CVMix_compute_SchmittnerCoeff low subroutines + + allocate(exp_hab_zetar(G%ke+1,G%ke+1)) + if (GV%Boussinesq) then + h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + else + h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + endif + + + do i=is,ie + + if (G%mask2dT(i,j)<1) cycle + + iFaceHeight = 0.0 ! BBL is all relative to the surface + hcorr = 0.0 + h_m = h(i,j,:)*GV%H_to_m + do k=1,G%ke + ! cell center and cell bottom in meters (negative values in the ocean) + dh = h_m(k) ! 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 + + SchmittnerSocn = 0.0 ! TODO: compute this + + ! form the time-invariant part of Schmittner coefficient term + call CVMix_compute_Schmittner_invariant(nlev = G%ke, & + VertDep = vert_dep, & + efficiency = CS%Mu_itides, & + rho = rho_fw, & + exp_hab_zetar = exp_hab_zetar, & + zw = iFaceHeight, & + CVmix_tidal_params_user = CS%CVMix_tidal_params) + !TODO: in above call, there is no need to pass efficiency, since it gets + ! passed via CVMix_init_tidal and stored in CVMix_tidal_params. Change + ! CVMix API to prevent this redundancy. + + ! remap from input z coordinate to model coordinate: + tidal_qe_md = 0.0 + call remapping_core_h(CS%remap_cs, size(CS%h_src), CS%h_src, CS%tidal_qe_3d_in(i,j,:), & + G%ke, h_m, tidal_qe_md) + + ! form the Schmittner coefficient that is based on 3D q*E, which is formed from + ! summing q_i*TidalConstituent_i over the number of constituents. + call CVMix_compute_SchmittnerCoeff( nlev = G%ke, & + energy_flux = tidal_qe_md(:), & + rho = rho_fw, & + SchmittnerCoeff = Schmittner_coeff, & + exp_hab_zetar = exp_hab_zetar, & + CVmix_tidal_params_user = CS%CVMix_tidal_params) + + + call CVMix_coeffs_tidal_schmittner( Mdiff_out = Kv_tidal, & + Tdiff_out = Kd_tidal, & + Nsqr = N2_int(i,:), & + OceanDepth = -iFaceHeight(G%ke+1), & + vert_dep = vert_dep, & + nlev = G%ke, & + max_nlev = G%ke, & + SchmittnerCoeff = Schmittner_coeff, & + SchmittnerSouthernOcean = SchmittnerSocn, & + CVmix_params = CS%CVMix_glb_params, & + CVmix_tidal_params_user = CS%CVMix_tidal_params) + + do k=1,G%ke + Kd(i,j,k) = Kd(i,j,k) + 0.5*(Kd_tidal(k) + Kd_tidal(k+1) ) + !TODO: Kv(i,j,k) = ???????????? + enddo + + ! diagnostics + if (associated(dd%Kd_itidal)) then + dd%Kd_itidal(i,j,:) = Kd_tidal(:) + endif + if (associated(dd%N2_int)) then + dd%N2_int(i,j,:) = N2_int(i,:) + endif + if (associated(dd%Schmittner_coeff_3d)) then + dd%Schmittner_coeff_3d(i,j,:) = Schmittner_coeff(:) + endif + if (associated(dd%tidal_qe_md)) then + dd%tidal_qe_md(i,j,:) = tidal_qe_md(:) + endif + if (associated(dd%vert_dep_3d)) then + dd%vert_dep_3d(i,j,:) = vert_dep(:) + endif + enddo ! i=is,ie + + deallocate(exp_hab_zetar) + case default - call MOM_error(FATAL, "tidal_mixing_init: The selected"// & - " INT_TIDE_PROFILE is unavailable in CVMix") + call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & + "#define CVMIX_TIDAL_SCHEME found in input file.") end select end subroutine calculate_CVMix_tidal @@ -1193,11 +1330,29 @@ subroutine setup_tidal_diagnostics(G,CS) allocate(dd%N2_int(isd:ied,jsd:jed,nz+1)) ; dd%N2_int(:,:,:) = 0.0 endif if (CS%id_Simmons_coeff > 0) then + if (CS%CVMix_tidal_scheme .ne. SIMMONS) then + call MOM_error(FATAL, "setup_tidal_diagnostics: Simmons_coeff diagnostics is available "//& + "only when CVMix_tidal_scheme is Simmons") + endif allocate(dd%Simmons_coeff_2d(isd:ied,jsd:jed)) ; dd%Simmons_coeff_2d(:,:) = 0.0 endif if (CS%id_vert_dep > 0) then allocate(dd%vert_dep_3d(isd:ied,jsd:jed,nz+1)) ; dd%vert_dep_3d(:,:,:) = 0.0 endif + if (CS%id_Schmittner_coeff > 0) then + if (CS%CVMix_tidal_scheme .ne. SCHMITTNER) then + call MOM_error(FATAL, "setup_tidal_diagnostics: Schmittner_coeff diagnostics is available "//& + "only when CVMix_tidal_scheme is Schmittner.") + endif + allocate(dd%Schmittner_coeff_3d(isd:ied,jsd:jed,nz)) ; dd%Schmittner_coeff_3d(:,:,:) = 0.0 + endif + if (CS%id_tidal_qe_md > 0) then + if (CS%CVMix_tidal_scheme .ne. SCHMITTNER) then + call MOM_error(FATAL, "setup_tidal_diagnostics: tidal_qe_md diagnostics is available "//& + "only when CVMix_tidal_scheme is Schmittner.") + endif + allocate(dd%tidal_qe_md(isd:ied,jsd:jed,nz)) ; dd%tidal_qe_md(:,:,:) = 0.0 + endif end subroutine setup_tidal_diagnostics subroutine post_tidal_diagnostics(G,GV,h,CS) @@ -1232,6 +1387,8 @@ subroutine post_tidal_diagnostics(G,GV,h,CS) if (CS%id_N2_int> 0) call post_data(CS%id_N2_int, dd%N2_int, CS%diag) if (CS%id_vert_dep> 0) call post_data(CS%id_vert_dep, dd%vert_dep_3d, CS%diag) if (CS%id_Simmons_coeff> 0) call post_data(CS%id_Simmons_coeff, dd%Simmons_coeff_2d, CS%diag) + if (CS%id_Schmittner_coeff> 0) call post_data(CS%id_Schmittner_coeff, dd%Schmittner_coeff_3d, CS%diag) + if (CS%id_tidal_qe_md> 0) call post_data(CS%id_tidal_qe_md, dd%tidal_qe_md, CS%diag) if (CS%id_Kd_Itidal_Work > 0) & call post_data(CS%id_Kd_Itidal_Work, dd%Kd_Itidal_Work, CS%diag) @@ -1283,6 +1440,8 @@ subroutine post_tidal_diagnostics(G,GV,h,CS) if (associated(dd%N2_int)) deallocate(dd%N2_int) if (associated(dd%vert_dep_3d)) deallocate(dd%vert_dep_3d) if (associated(dd%Simmons_coeff_2d)) deallocate(dd%Simmons_coeff_2d) + if (associated(dd%Schmittner_coeff_3d)) deallocate(dd%Schmittner_coeff_3d) + if (associated(dd%tidal_qe_md)) deallocate(dd%tidal_qe_md) end subroutine post_tidal_diagnostics @@ -1293,29 +1452,134 @@ subroutine read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) character(len=200), intent(in) :: tidal_energy_file type(tidal_mixing_cs), pointer :: CS ! local - integer :: isd, ied, jsd, jed + integer :: isd, ied, jsd, jed, nz real, allocatable, dimension(:,:) :: tidal_energy_flux_2d ! input tidal energy flux at T-grid points (W/m^2) - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - - if (.not. allocated(CS%tidal_qe_2d)) allocate(CS%tidal_qe_2d(isd:ied,jsd:jed)) - allocate(tidal_energy_flux_2d(isd:ied,jsd:jed)) + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke select case (uppercase(tidal_energy_type(1:4))) - case ('JAYN') ! Jayne 2009 input tidal energy flux + case ('JAYN') ! Jayne 2009 + if (.not. allocated(CS%tidal_qe_2d)) allocate(CS%tidal_qe_2d(isd:ied,jsd:jed)) + allocate(tidal_energy_flux_2d(isd:ied,jsd:jed)) call MOM_read_data(tidal_energy_file,'wave_dissipation',tidal_energy_flux_2d, G%domain) CS%tidal_qe_2d = (CS%Gamma_itides) * tidal_energy_flux_2d + deallocate(tidal_energy_flux_2d) + case ('ER03') ! Egbert & Ray 2003 + call read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) case default call MOM_error(FATAL, "read_tidal_energy: Unknown tidal energy file type.") - ! TODO: add more tidal energy file types, e.g., Arbic, ER03, GN13, LGM0, etc. - ! see POP::tidal_mixing.F90 end select - deallocate(tidal_energy_flux_2d) - end subroutine read_tidal_energy +subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + character(len=20), intent(in) :: tidal_energy_type + character(len=200), intent(in) :: tidal_energy_file + type(tidal_mixing_cs), pointer :: CS + + ! local + integer :: k, isd, ied, jsd, jed, i,j + integer, dimension(4) :: nz_in + real, parameter :: p33 = 1.0/3.0 + real, dimension(SZI_(G),SZJ_(G)) :: & + tidal_qk1, & ! qk1 coefficient used in Schmittner & Egbert + tidal_qo1 ! qo1 coefficient used in Schmittner & Egbert + real, allocatable, dimension(:) :: & + z_t, & ! depth from surface to midpoint of input layer [cm] + z_w ! depth from surface to top of input layer [cm] + real, allocatable, dimension(:,:,:) :: & + tc_m2, & ! input lunar semidiurnal tidal energy flux [W/m^2] + tc_s2, & ! input solar semidiurnal tidal energy flux [W/m^2] + tc_k1, & ! input lunar diurnal tidal energy flux [W/m^2] + tc_o1 ! input lunar diurnal tidal energy flux [W/m^2] + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + + ! get number of input levels: + call field_size(tidal_energy_file, 'z_t',nz_in) + + ! allocate local variables + allocate(z_t(nz_in(1)), z_w(nz_in(1)) ) + allocate(tc_m2(isd:ied,jsd:jed,nz_in(1)), & + tc_s2(isd:ied,jsd:jed,nz_in(1)), & + tc_k1(isd:ied,jsd:jed,nz_in(1)), & + tc_o1(isd:ied,jsd:jed,nz_in(1)) ) + + ! allocate CS variables associated with 3d tidal energy dissipation + if (.not. allocated(CS%tidal_qe_3d_in)) allocate(CS%tidal_qe_3d_in(isd:ied,jsd:jed,nz_in(1))) + if (.not. allocated(CS%h_src)) allocate(CS%h_src(nz_in(1))) + + ! read in tidal constituents + call MOM_read_data(tidal_energy_file, 'M2', tc_m2, G%domain) + call MOM_read_data(tidal_energy_file, 'S2', tc_s2, G%domain) + call MOM_read_data(tidal_energy_file, 'K1', tc_k1, G%domain) + call MOM_read_data(tidal_energy_file, 'O1', tc_o1, G%domain) + call MOM_read_data(tidal_energy_file, 'z_t', z_t) + call MOM_read_data(tidal_energy_file, 'z_w', z_w) + + where (abs(G%geoLatT(:,:)) < 30.0) + tidal_qk1(:,:) = p33 + tidal_qo1(:,:) = p33 + elsewhere + tidal_qk1(:,:) = 1.0 + tidal_qo1(:,:) = 1.0 + endwhere + + CS%tidal_qe_3d_in = 0.0 + do k=1,nz_in(1) + ! input cell thickness + CS%h_src(k) = (z_t(k)-z_w(k))*2.0 *1e-2 + ! form tidal_qe_3d_in from weighted tidal constituents + where ( (z_t(k)*1e-2) <= G%bathyT(:,:) .and. (z_w(k)*1e-2) > CS%tidal_diss_lim_tc) + CS%tidal_qe_3d_in(:,:,k) = p33*tc_m2(:,:,k) + p33*tc_s2(:,:,k) + & + tidal_qk1*tc_k1(:,:,k) + tidal_qo1*tc_o1(:,:,k) + endwhere + enddo + + !open(unit=1905,file="out_1905.txt",access="APPEND") + !do j=G%jsd,G%jed + ! do i=isd,ied + ! if ( i+G%idg_offset .eq. 90 .and. j+G%jdg_offset .eq. 126) then + ! print *, "-------------------------------------------" + ! do k=50,nz_in(1) + ! write(1905,*) i,j,k + ! write(1905,*) CS%tidal_qe_3d_in(i,j,k), tc_m2(i,j,k) + ! write(1905,*) z_t(k), G%bathyT(i,j), z_w(k),CS%tidal_diss_lim_tc + ! end do + ! endif + ! enddo + !enddo + !close(1905) + + ! test if qE is positive + if (any(CS%tidal_qe_3d_in<0.0)) then + call MOM_error(FATAL, "read_tidal_constituents: Negative tidal_qe_3d_in terms.") + endif + + !! collapse 3D q*E to 2D q*E + !CS%tidal_qe_2d = 0.0 + !do k=1,nz_in(1) + ! where (z_t(k) <= G%bathyT(:,:)) + ! CS%tidal_qe_2d(:,:) = CS%tidal_qe_2d(:,:) + CS%tidal_qe_3d_in(:,:,k) + ! endwhere + !enddo + + ! initialize input remapping: + call initialize_remapping(CS%remap_cs, remapping_scheme="PLM", & + boundary_extrapolation=.false., check_remapping=CS%debug) + + deallocate(tc_m2) + deallocate(tc_s2) + deallocate(tc_k1) + deallocate(tc_o1) + deallocate(z_t) + deallocate(z_w) + +end subroutine read_tidal_constituents + + !> Clear pointers and deallocate memory subroutine tidal_mixing_end(CS) type(tidal_mixing_cs), pointer :: CS ! This module's control structure @@ -1323,7 +1587,9 @@ subroutine tidal_mixing_end(CS) if (.not.associated(CS)) return !TODO deallocate all the dynamically allocated members here ... - if (allocated(CS%tidal_qe_2d)) deallocate(CS%tidal_qe_2d) + if (allocated(CS%tidal_qe_2d)) deallocate(CS%tidal_qe_2d) + if (allocated(CS%tidal_qe_3d_in)) deallocate(CS%tidal_qe_3d_in) + if (allocated(CS%h_src)) deallocate(CS%h_src) deallocate(CS%dd) deallocate(CS) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 48a6380ead..bafbe5eb59 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -2,7 +2,7 @@ module MOM_vert_friction ! This file is part of MOM6. See LICENSE.md for the license. - +use MOM_domains, only : pass_var, To_All, Omit_corners use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl use MOM_debugging, only : uvchksum, hchksum @@ -116,6 +116,7 @@ module MOM_vert_friction integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 integer :: id_Ray_u = -1, id_Ray_v = -1, id_taux_bot = -1, id_tauy_bot = -1 + integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 !>@} type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() @@ -614,6 +615,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) ! based on harmonic mean thicknesses, in m or kg m-2. h_ml ! The mixed layer depth, in m or kg m-2. real, allocatable, dimension(:,:) :: hML_u, hML_v + real, allocatable, dimension(:,:,:) :: Kv_v, & !< Total vertical viscosity at u-points + Kv_u !< Total vertical viscosity at v-points real :: zcol(SZI_(G)) ! The height of an interface at h-points, in m or kg m-2. real :: botfn ! A function which goes from 1 at the bottom to 0 much more ! than Hbbl into the interior. @@ -646,6 +649,14 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) I_Hbbl(:) = 1.0 / (CS%Hbbl * GV%m_to_H + h_neglect) I_valBL = 0.0 ; if (CS%harm_BL_val > 0.0) I_valBL = 1.0 / CS%harm_BL_val + if (CS%id_Kv_u > 0) then + allocate(Kv_u(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) ; Kv_u(:,:,:) = 0.0 + endif + + if (CS%id_Kv_v > 0) then + allocate(Kv_v(G%isd:G%ied,G%JsdB:G%JedB,G%ke)) ; Kv_v(:,:,:) = 0.0 + endif + if (CS%debug .or. (CS%id_hML_u > 0)) then allocate(hML_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; hML_u(:,:) = 0.0 endif @@ -821,6 +832,13 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) CS%h_u(I,j,k) = hvel(I,k) ; enddo ; enddo endif + ! Diagnose total Kv at u-points + if (CS%id_Kv_u > 0) then + do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) Kv_u(I,j,k) = 0.5 * (CS%a_u(I,j,K)+CS%a_u(I,j,K+1)) * CS%h_u(I,j,k) + enddo ; enddo + endif + enddo @@ -984,6 +1002,14 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) do K=1,nz+1 ; do i=is,ie ; if (do_i(i)) CS%a_v(i,J,K) = a(i,K) ; enddo ; enddo do k=1,nz ; do i=is,ie ; if (do_i(i)) CS%h_v(i,J,k) = hvel(i,k) ; enddo ; enddo endif + + ! Diagnose total Kv at v-points + if (CS%id_Kv_v > 0) then + do k=1,nz ; do i=is,ie + if (do_i(I)) Kv_v(i,J,k) = 0.5 * (CS%a_v(i,J,K)+CS%a_v(i,J,K+1)) * CS%h_v(i,J,k) + enddo ; enddo + endif + enddo ! end of v-point j loop if (CS%debug) then @@ -997,6 +1023,9 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) endif ! Offer diagnostic fields for averaging. + if (CS%id_Kv_slow > 0) call post_data(CS%id_Kv_slow, visc%Kv_slow, CS%diag) + if (CS%id_Kv_u > 0) call post_data(CS%id_Kv_u, Kv_u, CS%diag) + if (CS%id_Kv_v > 0) call post_data(CS%id_Kv_v, Kv_v, CS%diag) if (CS%id_au_vv > 0) call post_data(CS%id_au_vv, CS%a_u, CS%diag) if (CS%id_av_vv > 0) call post_data(CS%id_av_vv, CS%a_v, CS%diag) if (CS%id_h_u > 0) call post_data(CS%id_h_u, CS%h_u, CS%diag) @@ -1165,6 +1194,44 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m endif endif + ! add "slow" varying vertical viscosity (e.g., from background, tidal etc) + if (associated(visc%Kv_slow) .and. (visc%add_Kv_slow)) then + ! GMM/ A factor of 2 is also needed here, see comment above from BGR. + if (work_on_u) then + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + Kv_add(i,K) = Kv_add(i,K) + 1.0 * (visc%Kv_slow(i,j,k) + visc%Kv_slow(i+1,j,k)) + endif ; enddo ; enddo + if (do_OBCs) then + do I=is,ie ; if (do_i(I) .and. (OBC%segnum_u(I,j) /= OBC_NONE)) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + 2. * visc%Kv_slow(i,j,k) ; enddo + elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + 2. * visc%Kv_slow(i+1,j,k) ; enddo + endif + endif ; enddo + endif + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + a(i,K) = a(i,K) + Kv_add(i,K) + endif ; enddo ; enddo + else + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + Kv_add(i,K) = Kv_add(i,K) + 1.0*(visc%Kv_slow(i,j,k) + visc%Kv_slow(i,j+1,k)) + endif ; enddo ; enddo + if (do_OBCs) then + do i=is,ie ; if (do_i(i) .and. (OBC%segnum_v(i,J) /= OBC_NONE)) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + 2. * visc%Kv_slow(i,j,k) ; enddo + elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + 2. * visc%Kv_slow(i,j+1,k) ; enddo + endif + endif ; enddo + endif + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + a(i,K) = a(i,K) + Kv_add(i,K) + endif ; enddo ; enddo + endif + endif + do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then ! botfn determines when a point is within the influence of the bottom ! boundary layer, going from 1 at the bottom to 0 in the interior. @@ -1671,17 +1738,30 @@ subroutine vertvisc_init(MIS, Time, G, GV, param_file, diag, ADp, dirs, & ALLOC_(CS%a_v(isd:ied,JsdB:JedB,nz+1)) ; CS%a_v(:,:,:) = 0.0 ALLOC_(CS%h_v(isd:ied,JsdB:JedB,nz)) ; CS%h_v(:,:,:) = 0.0 + CS%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, Time, & + 'Slow varying vertical viscosity', 'm2 s-1') + + CS%id_Kv_u = register_diag_field('ocean_model', 'Kv_u', diag%axesCuL, Time, & + 'Total vertical viscosity at u-points', 'm2 s-1') + + CS%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, Time, & + 'Total vertical viscosity at v-points', 'm2 s-1') + CS%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, Time, & 'Zonal Viscous Vertical Coupling Coefficient', 'm s-1') + CS%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, Time, & 'Meridional Viscous Vertical Coupling Coefficient', 'm s-1') CS%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, Time, & 'Thickness at Zonal Velocity Points for Viscosity', thickness_units) + CS%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, Time, & 'Thickness at Meridional Velocity Points for Viscosity', thickness_units) + CS%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, Time, & 'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', thickness_units) + CS%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, Time, & 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', thickness_units)