From 71ff26ba73f886348240027fb0843854fa887a8e Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 9 Mar 2018 10:18:56 -0700 Subject: [PATCH 01/40] Changed multi-word names from camelCase to snake_case --- src/parameterizations/vertical/MOM_cvmix_shear.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_shear.F90 b/src/parameterizations/vertical/MOM_cvmix_shear.F90 index 460dde7c47..7704069d78 100644 --- a/src/parameterizations/vertical/MOM_cvmix_shear.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_shear.F90 @@ -25,10 +25,10 @@ module MOM_cvmix_shear #include -public Calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_is_used +public calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_is_used !> Control structure including parameters for CVMix interior shear schemes. -type, public :: CVMix_shear_CS +type, public :: cvmix_shear_cs logical :: use_LMD94, use_PP81 !< Flags for various schemes real :: Ri_zero !< LMD94 critical Richardson number real :: Nu_zero !< LMD94 maximum interior diffusivity @@ -36,14 +36,14 @@ module MOM_cvmix_shear real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) real, allocatable, dimension(:,:,:) :: S2 !< Squared shear frequency (1/s2) character(10) :: Mix_Scheme !< Mixing scheme name (string) -end type CVMix_shear_CS +end type cvmix_shear_cs character(len=40) :: mdl = "MOM_CVMix_shear" !< This module's name. contains !> Subroutine for calculating (internal) diffusivity -subroutine Calculate_cvmix_shear(u_H, v_H, h, tv, KH, & +subroutine calculate_cvmix_shear(u_H, v_H, h, tv, KH, & KM, G, GV, CS ) type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. @@ -55,7 +55,7 @@ subroutine Calculate_cvmix_shear(u_H, v_H, h, tv, KH, & !! (not layer!) in m2 s-1. real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: KM !< The vertical viscosity at each interface !! (not layer!) in m2 s-1. - type(CVMix_shear_CS), pointer :: CS !< The control structure returned by a previous call to + type(cvmix_shear_cs), pointer :: CS !< The control structure returned by a previous call to !! CVMix_shear_init. ! Local variables integer :: i, j, k, kk, km1 @@ -120,7 +120,7 @@ subroutine Calculate_cvmix_shear(u_H, v_H, h, tv, KH, & enddo enddo -end subroutine Calculate_cvmix_shear +end subroutine calculate_cvmix_shear !> Initialized the cvmix internal shear mixing routine. @@ -133,7 +133,7 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. - type(CVMix_shear_CS), pointer :: CS !< This module's control structure. + type(cvmix_shear_cs), pointer :: CS !< This module's control structure. ! Local variables integer :: NumberTrue=0 logical :: use_JHL From c6327d69cb7001de76eb21f970023eca9a9f87b3 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 12 Mar 2018 14:06:25 -0600 Subject: [PATCH 02/40] Copy OBLdepth from KPP into visc%MLD and Hml --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 91b3c343e0..0c1b3d1256 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -635,6 +635,14 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_turb, CS%KPP_NLTheat, CS%KPP_NLTscalar) !$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat) + ! If visc%MLD exists, copy the KPP BLD into it + if (associated(visc%MLD)) then + call pass_var(CS%KPP_CSp%OBLdepth, G%domain, halo=1) + visc%MLD(:,:) = CS%KPP_CSp%OBLdepth(:,:) + Hml(:,:) = CS%KPP_CSp%OBLdepth(:,:) + endif + + if (.not. CS%KPPisPassive) then !$OMP do do k=1,nz+1 ; do j=js,je ; do i=is,ie From dca5736441c39253a0052903779a1a841126565a Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 12 Mar 2018 15:25:34 -0600 Subject: [PATCH 03/40] Add comments for adding additonal CVMix components and clean module --- .../vertical/MOM_set_diffusivity.F90 | 83 +++++-------------- 1 file changed, 22 insertions(+), 61 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 15cbd6bb1c..2af9b8ef79 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -2,28 +2,6 @@ module MOM_set_diffusivity ! This file is part of MOM6. See LICENSE.md for the license. -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, September 1997 - June 2007 * -!* * -!* This file contains the subroutines that sets the diapycnal * -!* diffusivity, perhaps adding up pieces that are calculated in other * -!* files and passed in via the vertvisc type argument. * -!* * -!* A small fragment of the grid is shown below: * -!* * -!* j+1 x ^ x ^ x At x: q * -!* j+1 > o > o > At ^: v * -!* j x ^ x ^ x At >: u * -!* j > o > o > At o: h, buoy, ustar, T, S, Kd, ea, eb, etc. * -!* j-1 x ^ x ^ x * -!* i-1 i i+1 At x & ^: * -!* i i+1 At > & o: * -!* * -!* The boundaries always run through q grid points (x). * -!* * -!********+*********+*********+*********+*********+*********+*********+** - use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, time_type @@ -41,12 +19,11 @@ module MOM_set_diffusivity use MOM_intrinsic_functions, only : invcosh use MOM_io, only : slasher, vardesc, var_desc, MOM_read_data use MOM_kappa_shear, only : calculate_kappa_shear, kappa_shear_init, Kappa_shear_CS -use MOM_cvmix_shear, only : calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_CS +use MOM_cvmix_shear, only : calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_cs use MOM_string_functions, only : uppercase use MOM_thickness_diffuse, only : vert_fill_TS use MOM_variables, only : thermo_var_ptrs, vertvisc_type, p3d -use MOM_verticalGrid, only : verticalGrid_type - +use MOM_verticalGrid, only : verticalGrid_type use user_change_diffusivity, only : user_change_diff, user_change_diff_init use user_change_diffusivity, only : user_change_diff_end, user_change_diff_CS @@ -244,11 +221,9 @@ module MOM_set_diffusivity logical :: user_change_diff ! If true, call user-defined code to change diffusivity. logical :: useKappaShear ! If true, use the kappa_shear module to find the ! shear-driven diapycnal diffusivity. - - logical :: useCVmix ! If true, use one of the CVMix modules to find + logical :: use_cvmix_shear ! If true, use one of the CVMix modules to find ! shear-driven diapycnal diffusivity. - - logical :: double_diffusion ! If true, enable double-diffusive mixing. + logical :: double_diffusion ! If true, enable double-diffusive mixing. logical :: simple_TKE_to_Kd ! If true, uses a simple estimate of Kd/TKE that ! does not rely on a layer-formulation. real :: Max_Rrho_salt_fingers ! max density ratio for salt fingering @@ -266,7 +241,7 @@ module MOM_set_diffusivity type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() type(Kappa_shear_CS), pointer :: kappaShear_CSp => NULL() - type(CVMix_shear_CS), pointer :: CVMix_Shear_CSp => NULL() + type(cvmix_shear_cs), pointer :: cvmix_shear_csp => NULL() type(int_tide_CS), pointer :: int_tide_CSp => NULL() integer :: id_TKE_itidal = -1 @@ -380,22 +355,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & optional, intent(out) :: Kd_int !< Diapycnal diffusivity at each interface !! (m2/sec). -! Arguments: -! (in) u - zonal velocity (m/s) -! (in) v - meridional velocity (m/s) -! (in) h - Layer thickness (m or kg/m2) -! (in) tv - structure with pointers to thermodynamic fields -! (in) fluxes - structure of surface fluxes that may be used -! (in) visc - structure containing vertical viscosities, bottom boundary -! layer properies, and related fields -! (in) dt - time increment (sec) -! (in) G - ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) CS - module control structure -! (in) j - meridional index upon which to work -! (out) Kd - diapycnal diffusivity of each layer (m2/sec) -! (out,opt) Kd_int - diapycnal diffusivity at each interface (m2/sec) - + ! local variables real, dimension(SZI_(G)) :: & depth, & ! distance from surface of an interface (meter) N2_bot ! bottom squared buoyancy frequency (1/s2) @@ -580,9 +540,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & call hchksum(visc%TKE_turb, "after calc_KS visc%TKE_turb",G%HI) endif if (showCallTree) call callTree_waypoint("done with calculate_kappa_shear (set_diffusivity)") - elseif (CS%useCVMix) then + elseif (CS%use_cvmix_shear) then !NOTE{BGR}: this needs to be cleaned up. It works in 1D case, but has not been tested outside. - call calculate_cvmix_shear(u_h, v_h, h, tv, visc%Kd_turb, visc%Kv_turb,G,GV,CS%CVMix_shear_CSp) + call calculate_cvmix_shear(u_h, v_h, h, tv, visc%Kd_turb, visc%Kv_turb,G,GV,CS%cvmix_shear_csp) elseif (associated(visc%Kv_turb)) then visc%Kv_turb(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled endif @@ -674,6 +634,10 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & Kd(i,j,k) = max(CS%Kd_min, Kd_sfc(i,j) * & ((abs_sin * invcosh(N_2Omega/abs_sin)) * I_x30)*N02_N2) enddo ; enddo + ! GMM, CVMix "internal" bg mixing can go here + !elseif (CS%use_cvmix_internal??) then + + else do k=1,nz ; do i=is,ie Kd(i,j,k) = Kd_sfc(i,j) @@ -708,7 +672,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif ! Add the input turbulent diffusivity. - if (CS%useKappaShear .or. CS%useCVMix) then + if (CS%useKappaShear .or. CS%use_cvmix_shear) then if (present(Kd_int)) then do K=2,nz ; do i=is,ie Kd_int(i,j,K) = visc%Kd_turb(i,j,K) + 0.5*(Kd(i,j,k-1) + Kd(i,j,k)) @@ -834,6 +798,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & T_f, S_f, dd%Kd_user) endif + ! GMM, post diags... if (CS%id_Kd_layer > 0) call post_data(CS%id_Kd_layer, Kd, CS%diag) num_z_diags = 0 @@ -2533,6 +2498,7 @@ subroutine set_density_ratios(h, tv, kb, G, GV, CS, j, ds_dsp1, rho_0) end subroutine set_density_ratios +!> Initialized the set_diffusivity module subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp, int_tide_CSp) type(time_type), intent(in) :: Time type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. @@ -2547,21 +2513,14 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp type(int_tide_CS), pointer :: int_tide_CSp !< pointer to the internal tides control !! structure (BDM) -! Arguments: -! (in) Time - current model time -! (in) G - ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) param_file - structure indicating open file to parse for params -! (in) diag - structure used to regulate diagnostic output -! (in/out) CS - pointer set to point to the module control structure -! (in) diag_to_Z_CSp - pointer to the Z-diagnostics control structure -! (in) int_tide_CSp - pointer to the internal tides control structure (BDM) - + ! local variables real :: decay_length, utide, zbot, hamp type(vardesc) :: vd logical :: read_tideamp, ML_use_omega + ! This include declares and sets the variable "version". #include "version_variable.h" + character(len=40) :: mdl = "MOM_set_diffusivity" ! This module's name. character(len=20) :: tmpstr character(len=200) :: filename, tideamp_file, h2_file, Niku_TKE_input_file @@ -3131,6 +3090,8 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp endif endif + + ! GMM, the following should be moved to the DD module call get_param(param_file, mdl, "DOUBLE_DIFFUSION", CS%double_diffusion, & "If true, increase diffusivitives for temperature or salt \n"//& "based on double-diffusive paramaterization from MOM4/KPP.", & @@ -3181,9 +3142,9 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp CS%useKappaShear = kappa_shear_init(Time, G, GV, param_file, CS%diag, CS%kappaShear_CSp) if (CS%useKappaShear) & id_clock_kappaShear = cpu_clock_id('(Ocean kappa_shear)', grain=CLOCK_MODULE) - CS%useCVMix = CVMix_shear_init(Time, G, GV, param_file, CS%diag, CS%CVMix_shear_CSp) - + ! CVMix shear-driven mixing + CS%use_cvmix_shear = cvmix_shear_init(Time, G, GV, param_file, CS%diag, CS%cvmix_shear_csp) end subroutine set_diffusivity_init From 97152b1c4bd12f85a067377f75c68677558101bc Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 12 Mar 2018 15:26:05 -0600 Subject: [PATCH 04/40] Adding first version of convection calls via CVMix --- .../vertical/MOM_cvmix_conv.F90 | 245 ++++++++++++++++++ .../vertical/MOM_diabatic_driver.F90 | 26 +- 2 files changed, 259 insertions(+), 12 deletions(-) create mode 100644 src/parameterizations/vertical/MOM_cvmix_conv.F90 diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 new file mode 100644 index 0000000000..7da68faa49 --- /dev/null +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -0,0 +1,245 @@ +!> Interface to CVMix convection scheme. +module MOM_cvmix_conv + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field +use MOM_diag_mediator, only : post_data +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_debugging, only : hchksum +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_file_parser, only : get_param, log_version, param_file_type +use cvmix_convection, only : cvmix_init_conv, cvmix_coeffs_conv +use cvmix_kpp, only : CVmix_kpp_compute_kOBL_depth +implicit none ; private + +#include + +public cvmix_conv_init, calculate_cvmix_conv, cvmix_conv_end + +!> Control structure including parameters for CVMix convection. +type, public :: cvmix_conv_cs + + ! Parameters + real :: kd_conv !< diffusivity constant used in convective regime (m2/s) + real :: kv_conv !< viscosity constant used in convective regime (m2/s) + real :: bv_sqr_conv !< Threshold for squared buoyancy frequency + !! needed to trigger Brunt-Vaisala parameterization (1/s^2) + real :: min_thickness !< Minimum thickness allowed (m) + logical :: debug !< If true, turn on debugging + + ! Daignostic handles and pointers + type(diag_ctrl), pointer :: diag => NULL() + integer :: id_N2 = -1, id_kd_conv = -1, id_kv_conv = -1 + + ! Diagnostics arrays + real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) + real, allocatable, dimension(:,:,:) :: kd_conv_3d !< Diffusivity added by convection (m2/s) + real, allocatable, dimension(:,:,:) :: kv_conv_3d !< Viscosity added by convection (m2/s) + +end type cvmix_conv_cs + +character(len=40) :: mdl = "MOM_cvmix_conv" !< This module's name. + +contains + +!> Initialized the cvmix convection mixing routine. +logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) + + type(time_type), intent(in) :: Time !< The current time. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. + type(cvmix_conv_cs), pointer :: CS !< This module's control structure. + + ! Local variables + real :: prandtl_turb + +! This include declares and sets the variable "version". +#include "version_variable.h" + + if (associated(CS)) then + call MOM_error(WARNING, "cvmix_conv_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + ! Read parameters + call log_version(param_file, mdl, version, & + "Parameterization of enhanced mixing due to convection via CVMix") + call get_param(param_file, mdl, "USE_CVMIX_CONVECTION", cvmix_conv_init, & + "If true, turns on the enhanced mixing due to convection \n"// & + "via CVMix. This scheme increases diapycnal diffs./viscs. \n"// & + " at statically unstable interfaces. Relevant parameters are \n"// & + "contained in the CVMIX_CONVECTION% parameter block.", & + default=.false.) + + if (.not. cvmix_conv_init) return + + call get_param(param_file, mdl, "PRANDTL_TURB", Prandtl_turb, & + "The turbulent Prandtl number applied to shear/conv. \n"//& + "instabilities.", units="nondim", default=1.0, do_not_log=.true.) + + call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) + + call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, do_not_log=.True.) + + call openParameterBlock(param_file,'CVMIX_CONVECTION') + + call get_param(param_file, mdl, 'KD_CONV', CS%kd_conv, & + "Diffusivity used in convective regime. Corresponding viscosity \n" // & + "(KV_CONV) will be set to KD_CONV * PRANDTL_TURB.", & + units='m2/s', default=1.00) + + call get_param(param_file, mdl, 'BV_SQR_CONV', CS%bv_sqr_conv, & + "Threshold for squared buoyancy frequency needed to trigger \n" // & + "Brunt-Vaisala parameterization.", & + units='1/s^2', default=0.0) + + call closeParameterBlock(param_file) + + ! set kv_conv based on kd_conv and Prandtl_turb + CS%kv_conv = CS%kd_conv * Prandtl_turb + + ! Register diagnostics + CS%diag => diag + CS%id_N2 = register_diag_field('ocean_model', 'conv_N2', diag%axesTi, Time, & + 'Square of Brunt-Vaisala frequency used by MOM_cvmix_conv module', '1/s2') + if (CS%id_N2 > 0) allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) + CS%id_kd_conv = register_diag_field('ocean_model', 'conv_kd', diag%axesTi, Time, & + 'Additional diffusivity added by MOM_cvmix_conv module', 'm2/s') + if (CS%id_kd_conv > 0) allocate( CS%kd_conv_3d( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) + CS%id_kv_conv = register_diag_field('ocean_model', 'conv_kv', diag%axesTi, Time, & + 'Additional viscosity added by MOM_cvmix_conv module', 'm2/s') + if (CS%id_kv_conv > 0) allocate( CS%kv_conv_3d( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) + + if (CS%id_N2 > 0) CS%N2(:,:,:) = 0. + if (CS%id_kd_conv > 0) CS%kd_conv_3d(:,:,:) = 0. + if (CS%id_kv_conv > 0) CS%kv_conv_3d(:,:,:) = 0. + + call cvmix_init_conv(convect_diff=CS%kd_conv, & + convect_visc=CS%kv_conv, & + lBruntVaisala=.true., & + BVsqr_convect=CS%bv_sqr_conv) + +end function cvmix_conv_init + +!> Subroutine for calculating enhanced diffusivity/viscosity +!! due to convection via CVMix +subroutine calculate_cvmix_conv(h, tv, G, GV, hbl, 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)) intent(in) :: hbl!< Depth of ocean boundary layer (m) + type(cvmix_conv_cs), pointer(inout) :: CS !< The control structure returned by a previous call to + !! CVMix_conv_init. + + ! local variables + real, dimension(SZK_(G)) :: rho_lwr !< Adiabatic Water Density, this is a dummy + !! variable since here convection is always + !! computed based on Brunt Vaisala. + real, dimension(SZK_(G)) :: rho_1d !< water density in a column, this is also + !! a dummy variable, same reason as above. + real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) + real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) + real :: kOBL !< level (+fraction) of OBL extent + real :: pref, g_o_rho0, rhok, , rhokm1, dz, dh, hcorr + integer :: i, j, k + + g_o_rho0 = GV%g_Earth / GV%Rho0 + + ! initialize dummy variables + rho_lwr(:) = 0.0; rho_1d(:) = 0.0 + + do j = G%jsc, G%jec + do i = G%isc, G%iec + + ! set N2 to zero at the top- and bottom-most interfaces + CS%N2(i,j,1) = 0. + CS%N2(i,j,G%ke+1) =0. + + ! skip calling at land points + if (G%mask2dT(i,j)==0.) cycle + + pRef = 0. + ! Compute Brunt-Vaisala frequency (static stability) on interfaces + do k=2,G%ke + + ! pRef is pressure at interface between k and km1. + pRef = pRef + GV%H_to_Pa * h(i,j,k) + call calculate_density(tv%t(i,j,k), tv%s(i,j,k), pref, rhok, tv%eqn_of_state) + call calculate_density(tv%t(i,j,k-1), tv%s(i,j,k-1), pref, rhokm1, tv%eqn_of_state) + + dz = ((0.5*(h(i,j,k-1) + h(i,j,k))+GV%H_subroundoff)*GV%H_to_m) + CS%N2(i,j,k) = g_o_rho0 * (rhok - rhokm1) / dz ! Can be negative + + enddo + + iFaceHeight(1) = 0.0 ! BBL is all relative to the surface + hcorr = 0.0 + ! compute heights at cell center and interfaces + do k=1,G%ke + dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) + hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 + dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * dh + iFaceHeight(k+1) = iFaceHeight(k) - dh + enddo + + kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl) + + call cvmix_coeffs_conv(Mdiff_out = CS%kv_conv_3d(i,j,:), & + Tdiff_out = CS%kd_conv_3d(i,j,:), & + Nsqr = CS%N2(i,j,:), & + dens = rho_1d(:), & + dens_lwr = rho_lwr(:), & + nlev = G%ke, & + max_nlev = G%ke, & + OBL_ind = kOBL) + + enddo + enddo + + if (CS%debug) then + call hchksum(CS%N2, "CVMix convection: N2",G%HI,haloshift=0) + call hchksum(CS%kd_conv_3d, "CVMix convection: kd_conv_3d",G%HI,haloshift=0) + call hchksum(CS%kv_conv_3d, "CVMix convection: kv_conv_3d",G%HI,haloshift=0) + endif + + ! send diagnostics to post_data + if (CS%id_N2 > 0) call post_data(CS%id_N2, CS%N2, CS%diag) + if (CS%id_kd_conv > 0) call post_data(CS%id_kd_conv, CS%kd_conv_3d, CS%diag) + if (CS%id_kv_conv > 0) call post_data(CS%id_kv_conv, CS%kv_conv_3d, CS%diag) + +end subroutine calculate_cvmix_conv + +! GMM, not sure if we need the code below - DELETE???? +!!logical function cvmix_conv_is_used(param_file) +! Reads the parameter "USE_CVMIX_CONVECTION" and returns state. +! This function allows other modules to know whether this parameterization will +! be used without needing to duplicate the log entry. +!! type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters +!! call get_param(param_file, mdl, "USE_CVMIX_CONVECTION", kappa_shear_is_used, & +!! default=.false., do_not_log = .true.) +!!end function cvmix_conv_is_used + +!> Clear pointers and dealocate memory +subroutine cvmix_conv_end(CS) + type(cvmix_conv_cs), pointer :: CS ! Control structure + + if (CS%id_N2 > 0) deallocate(CS%N2, CS%diag) + if (CS%id_kd_conv > 0) deallocate(CS%kd_conv_3d, CS%diag) + if (CS%id_kv_conv > 0) deallocate(CS%kv_conv_3d, CS%diag) + deallocate(CS) + +end subroutine cvmix_conv_end + + +end module MOM_cvmix_conv diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 0c1b3d1256..56824059a3 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -22,8 +22,8 @@ module MOM_diabatic_driver use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag, calc_Zint_diags use MOM_diapyc_energy_req, only : diapyc_energy_req_init, diapyc_energy_req_end use MOM_diapyc_energy_req, only : diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_CS -use MOM_diffConvection, only : diffConvection_CS, diffConvection_init -use MOM_diffConvection, only : diffConvection_calculate, diffConvection_end +use MOM_cvmix_conv, only : cvmix_conv_init, cvmix_conv_cs +use MOM_cvmix_conv, only : cvmix_conv_end, calculate_cvmix_conv use MOM_domains, only : pass_var, To_West, To_South, To_All, Omit_Corners use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_energetic_PBL, only : energetic_PBL, energetic_PBL_init @@ -90,6 +90,8 @@ module MOM_diabatic_driver !! shear-driven diapycnal diffusivity. logical :: use_cvmix_shear !< If true, use the CVMix module to find the !! shear-driven diapycnal diffusivity. + logical :: use_cvmix_conv !< If true, use the CVMix module to get enhanced + !! mixing due to convection. logical :: use_sponge !< If true, sponges may be applied anywhere in the !! domain. The exact location and properties of !! those sponges are set by calls to @@ -149,8 +151,6 @@ module MOM_diabatic_driver logical :: useKPP !< use CVmix/KPP diffusivities and non-local transport logical :: salt_reject_below_ML !< If true, add salt below mixed layer (layer mode only) logical :: KPPisPassive !< If true, KPP is in passive mode, not changing answers. - logical :: useConvection !< If true, calculate large diffusivities when column - !! is statically unstable. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: debugConservation !< If true, monitor conservation and extrema. logical :: tracer_tridiag !< If true, use tracer_vertdiff instead of tridiagTS for @@ -220,7 +220,7 @@ module MOM_diabatic_driver type(optics_type), pointer :: optics => NULL() type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() type(KPP_CS), pointer :: KPP_CSp => NULL() - type(diffConvection_CS), pointer :: Conv_CSp => NULL() + type(cvmix_conv_cs), pointer :: cvmix_conv_csp => NULL() type(diapyc_energy_req_CS), pointer :: diapyc_en_rec_CSp => NULL() type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass @@ -674,9 +674,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G endif ! endif for KPP - ! Check for static instabilities and increase Kd_int where unstable - if (CS%useConvection) call diffConvection_calculate(CS%Conv_CSp, & - G, GV, h, tv%T, tv%S, tv%eqn_of_state, Kd_int) + ! Add diffusivity due to convection (computed via CVMix) + if (CS%use_cvmix_conv) & + call calculate_cvmix_conv(h, tv, G, GV, Hml, CS%cvmix_conv_csp) if (CS%useKPP) then @@ -2332,9 +2332,9 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, endif - ! CS%useConvection is set to True IF convection will be used, otherwise False. - ! CS%Conv_CSp is allocated by diffConvection_init() - CS%useConvection = diffConvection_init(param_file, G, diag, Time, CS%Conv_CSp) + ! CS%use_cvmix_conv is set to True if CVMix convection will be used, otherwise + ! False. + CS%use_cvmix_conv = cvmix_conv_init(Time, G, GV, param_file, diag, CS%cvmix_conv_csp) call entrain_diffusive_init(Time, G, GV, param_file, diag, CS%entrain_diffusive_CSp) @@ -2422,7 +2422,9 @@ subroutine diabatic_driver_end(CS) deallocate( CS%KPP_NLTscalar ) call KPP_end(CS%KPP_CSp) endif - if (CS%useConvection) call diffConvection_end(CS%Conv_CSp) + + if (CS%use_cvmix_conv) call cvmix_conv_end(CS%cvmix_conv_csp) + if (CS%use_energetic_PBL) & call energetic_PBL_end(CS%energetic_PBL_CSp) if (CS%debug_energy_req) & From c99e94ba2bfee186b4f9abe056882345239a71fa Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 12 Mar 2018 15:26:55 -0600 Subject: [PATCH 05/40] Updates visc%Kd_turb and visc%Kv_turb after convection is applied --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 56824059a3..bfe9c3a662 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -675,9 +675,16 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G endif ! endif for KPP ! Add diffusivity due to convection (computed via CVMix) - if (CS%use_cvmix_conv) & + if (CS%use_cvmix_conv) then call calculate_cvmix_conv(h, tv, G, GV, Hml, CS%cvmix_conv_csp) + do k=1,nz ; do j=js,je ; do i=is,ie + visc%Kd_turb(i,j,k) = visc%Kd_turb(i,j,k) + CS%cvmix_conv_csp%kd_conv_3d(i,j,k) + visc%Kv_turb(i,j,k) = visc%Kv_turb(i,j,k) + CS%cvmix_conv_csp%kv_conv_3d(i,j,k) + enddo ; enddo ; enddo + + endif + if (CS%useKPP) then call cpu_clock_begin(id_clock_kpp) From fc00e77aa83471c017a24825b0f909a42cf2a55c Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 13 Mar 2018 09:11:32 -0600 Subject: [PATCH 06/40] Added function to copy KPP surface boundary layer depth into BLD --- src/parameterizations/vertical/MOM_KPP.F90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index df64bb7b21..87ce532a28 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -32,6 +32,7 @@ module MOM_KPP public :: KPP_end public :: KPP_NonLocalTransport_temp public :: KPP_NonLocalTransport_saln +public :: KPP_get_BLD ! Enumerated constants integer, private, parameter :: NLT_SHAPE_CVMIX = 0 !< Use the CVmix profile @@ -936,7 +937,18 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & end subroutine KPP_calculate - +!> Copies KPP surface boundary layer depth into BLD +subroutine KPP_get_BLD(CS, BLD, G) + type(KPP_CS), pointer :: CS !< Control structure for + !! this module + type(ocean_grid_type), intent(in) :: G !< Grid structure + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: BLD!< bnd. layer depth (m) + ! Local variables + integer :: i,j + do j = G%jsc, G%jec ; do i = G%isc, G%iec + BLD(i,j) = CS%OBLdepth(i,j) + enddo ; enddo +end subroutine KPP_get_BLD !> Apply KPP non-local transport of surface fluxes for temperature. subroutine KPP_NonLocalTransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, & From 7f6d8f13bcc8305efe5ac1c39ddd6117dc4495f1 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 13 Mar 2018 09:12:20 -0600 Subject: [PATCH 07/40] Added call to KPP_get_BLD --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index bfe9c3a662..bba8f657f2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -47,7 +47,7 @@ module MOM_diabatic_driver use MOM_internal_tides, only : propagate_int_tide use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS use MOM_kappa_shear, only : kappa_shear_is_used -use MOM_KPP, only : KPP_CS, KPP_init, KPP_calculate, KPP_end +use MOM_KPP, only : KPP_CS, KPP_init, KPP_calculate, KPP_end, KPP_get_BLD use MOM_KPP, only : KPP_NonLocalTransport_temp, KPP_NonLocalTransport_saln use MOM_opacity, only : opacity_init, set_opacity, opacity_end, opacity_CS use MOM_regularize_layers, only : regularize_layers, regularize_layers_init, regularize_layers_CS @@ -637,9 +637,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G ! If visc%MLD exists, copy the KPP BLD into it if (associated(visc%MLD)) then - call pass_var(CS%KPP_CSp%OBLdepth, G%domain, halo=1) - visc%MLD(:,:) = CS%KPP_CSp%OBLdepth(:,:) - Hml(:,:) = CS%KPP_CSp%OBLdepth(:,:) + call KPP_get_BLD(CS%KPP_CSp, visc%MLD(:,:), G) + call pass_var(visc%MLD, G%domain, halo=1) + Hml(:,:) = visc%MLD(:,:) endif From e84f706de2f32036b2786a8cebc43a700d699771 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 13 Mar 2018 09:12:56 -0600 Subject: [PATCH 08/40] Fix array allocattion/dealocation * TODO: I am not sure if visc/diff due to convection is being added properly into the total visc./diff. This needs to be checked! --- .../vertical/MOM_cvmix_conv.F90 | 58 ++++++++++--------- 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index 7da68faa49..5fab297910 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -5,6 +5,8 @@ module MOM_cvmix_conv use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field use MOM_diag_mediator, only : post_data +use MOM_EOS, only : calculate_density +use MOM_variables, only : thermo_var_ptrs use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_debugging, only : hchksum @@ -13,6 +15,7 @@ module MOM_cvmix_conv use MOM_file_parser, only : get_param, log_version, param_file_type use cvmix_convection, only : cvmix_init_conv, cvmix_coeffs_conv use cvmix_kpp, only : CVmix_kpp_compute_kOBL_depth + implicit none ; private #include @@ -105,21 +108,19 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) ! set kv_conv based on kd_conv and Prandtl_turb CS%kv_conv = CS%kd_conv * Prandtl_turb + ! allocate arrays and set them to zero + allocate(CS%N2(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%N2(:,:,:) = 0. + allocate(CS%kd_conv_3d(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kd_conv_3d(:,:,:) = 0. + allocate(CS%kv_conv_3d(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kv_conv_3d(:,:,:) = 0. + ! Register diagnostics CS%diag => diag CS%id_N2 = register_diag_field('ocean_model', 'conv_N2', diag%axesTi, Time, & 'Square of Brunt-Vaisala frequency used by MOM_cvmix_conv module', '1/s2') - if (CS%id_N2 > 0) allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) CS%id_kd_conv = register_diag_field('ocean_model', 'conv_kd', diag%axesTi, Time, & 'Additional diffusivity added by MOM_cvmix_conv module', 'm2/s') - if (CS%id_kd_conv > 0) allocate( CS%kd_conv_3d( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) CS%id_kv_conv = register_diag_field('ocean_model', 'conv_kv', diag%axesTi, Time, & 'Additional viscosity added by MOM_cvmix_conv module', 'm2/s') - if (CS%id_kv_conv > 0) allocate( CS%kv_conv_3d( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) - - if (CS%id_N2 > 0) CS%N2(:,:,:) = 0. - if (CS%id_kd_conv > 0) CS%kd_conv_3d(:,:,:) = 0. - if (CS%id_kv_conv > 0) CS%kv_conv_3d(:,:,:) = 0. call cvmix_init_conv(convect_diff=CS%kd_conv, & convect_visc=CS%kv_conv, & @@ -136,8 +137,9 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, hbl, CS) 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)) intent(in) :: hbl!< Depth of ocean boundary layer (m) - type(cvmix_conv_cs), pointer(inout) :: CS !< The control structure returned by a previous call to + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hbl!< Depth of ocean boundary layer (m) + !type(cvmix_conv_cs), intent(inout) :: CS !< The control structure returned by a previous call to + type(cvmix_conv_cs), pointer :: CS !< The control structure returned by a previous call to !! CVMix_conv_init. ! local variables @@ -148,8 +150,8 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, hbl, CS) !! a dummy variable, same reason as above. real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) - real :: kOBL !< level (+fraction) of OBL extent - real :: pref, g_o_rho0, rhok, , rhokm1, dz, dh, hcorr + integer :: kOBL !< level of OBL extent + real :: pref, g_o_rho0, rhok, rhokm1, dz, dh, hcorr integer :: i, j, k g_o_rho0 = GV%g_Earth / GV%Rho0 @@ -162,10 +164,10 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, hbl, CS) ! set N2 to zero at the top- and bottom-most interfaces CS%N2(i,j,1) = 0. - CS%N2(i,j,G%ke+1) =0. + CS%N2(i,j,G%ke+1) = 0. ! skip calling at land points - if (G%mask2dT(i,j)==0.) cycle + !if (G%mask2dT(i,j) == 0.) cycle pRef = 0. ! Compute Brunt-Vaisala frequency (static stability) on interfaces @@ -193,24 +195,24 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, hbl, CS) iFaceHeight(k+1) = iFaceHeight(k) - dh enddo - kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl) + kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) - call cvmix_coeffs_conv(Mdiff_out = CS%kv_conv_3d(i,j,:), & - Tdiff_out = CS%kd_conv_3d(i,j,:), & - Nsqr = CS%N2(i,j,:), & - dens = rho_1d(:), & - dens_lwr = rho_lwr(:), & - nlev = G%ke, & - max_nlev = G%ke, & - OBL_ind = kOBL) + call cvmix_coeffs_conv(Mdiff_out=CS%kv_conv_3d(i,j,:), & + Tdiff_out=CS%kd_conv_3d(i,j,:), & + Nsqr=CS%N2(i,j,:), & + dens=rho_1d(:), & + dens_lwr=rho_lwr(:), & + nlev=G%ke, & + max_nlev=G%ke, & + OBL_ind=kOBL) enddo enddo if (CS%debug) then - call hchksum(CS%N2, "CVMix convection: N2",G%HI,haloshift=0) - call hchksum(CS%kd_conv_3d, "CVMix convection: kd_conv_3d",G%HI,haloshift=0) - call hchksum(CS%kv_conv_3d, "CVMix convection: kv_conv_3d",G%HI,haloshift=0) + call hchksum(CS%N2, "MOM_cvmix_conv: N2",G%HI,haloshift=0) + call hchksum(CS%kd_conv_3d, "MOM_cvmix_conv: kd_conv_3d",G%HI,haloshift=0) + call hchksum(CS%kv_conv_3d, "MOM_cvmix_conv: kv_conv_3d",G%HI,haloshift=0) endif ! send diagnostics to post_data @@ -234,9 +236,9 @@ end subroutine calculate_cvmix_conv subroutine cvmix_conv_end(CS) type(cvmix_conv_cs), pointer :: CS ! Control structure - if (CS%id_N2 > 0) deallocate(CS%N2, CS%diag) - if (CS%id_kd_conv > 0) deallocate(CS%kd_conv_3d, CS%diag) - if (CS%id_kv_conv > 0) deallocate(CS%kv_conv_3d, CS%diag) + deallocate(CS%N2, CS%diag) + deallocate(CS%kd_conv_3d, CS%diag) + deallocate(CS%kv_conv_3d, CS%diag) deallocate(CS) end subroutine cvmix_conv_end From a209b7a8a02a61201f3c505ae48e904c2f26d2a5 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 13 Mar 2018 10:36:59 -0600 Subject: [PATCH 09/40] Added a call to diabatic_driver_end --- src/core/MOM.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 60ce1d9c99..bcf6516505 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2900,6 +2900,8 @@ subroutine MOM_end(CS) call tracer_registry_end(CS%tracer_Reg) call tracer_flow_control_end(CS%tracer_flow_CSp) + if (associated(CS%diabatic_CSp)) call diabatic_driver_end(CS%diabatic_CSp) + if (CS%offline_tracer_mode) call offline_transport_end(CS%offline_CSp) DEALLOC_(CS%uhtr) ; DEALLOC_(CS%vhtr) From f12702c58a54393473dbb18890b4a03d128efe94 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 13 Mar 2018 10:40:51 -0600 Subject: [PATCH 10/40] Dealocate memory via cvmix_shear_end and changed some var names to snake_case convention --- src/parameterizations/vertical/MOM_cvmix_shear.F90 | 14 ++++++++++++-- .../vertical/MOM_diabatic_driver.F90 | 6 +++--- .../vertical/MOM_set_diffusivity.F90 | 7 ++++++- 3 files changed, 21 insertions(+), 6 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_shear.F90 b/src/parameterizations/vertical/MOM_cvmix_shear.F90 index 7704069d78..487b1cbf6f 100644 --- a/src/parameterizations/vertical/MOM_cvmix_shear.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_shear.F90 @@ -25,7 +25,7 @@ module MOM_cvmix_shear #include -public calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_is_used +public calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_is_used, cvmix_shear_end !> Control structure including parameters for CVMix interior shear schemes. type, public :: cvmix_shear_cs @@ -38,7 +38,7 @@ module MOM_cvmix_shear character(10) :: Mix_Scheme !< Mixing scheme name (string) end type cvmix_shear_cs -character(len=40) :: mdl = "MOM_CVMix_shear" !< This module's name. +character(len=40) :: mdl = "MOM_cvmix_shear" !< This module's name. contains @@ -213,4 +213,14 @@ logical function cvmix_shear_is_used(param_file) cvmix_shear_is_used = (LMD94 .or. PP81) end function cvmix_shear_is_used +!> Clear pointers and dealocate memory +subroutine cvmix_shear_end(CS) + type(cvmix_shear_cs), pointer :: CS ! Control structure + + deallocate(CS%N2, CS%diag) + deallocate(CS%S2, CS%diag) + deallocate(CS) + +end subroutine cvmix_shear_end + end module MOM_cvmix_shear diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index bba8f657f2..aa4f9f072e 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -9,7 +9,7 @@ module MOM_diabatic_driver use MOM_checksum_packages, only : MOM_state_chksum, MOM_state_stats use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE -use MOM_CVMix_shear, only : CVMix_shear_is_used +use MOM_CVMix_shear, only : cvmix_shear_is_used use MOM_diabatic_aux, only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_CS use MOM_diabatic_aux, only : make_frazil, adjust_salt, insert_brine, differential_diffuse_T_S, triDiagTS use MOM_diabatic_aux, only : find_uv_at_h, diagnoseMLDbyDensityDifference, applyBoundaryFluxesInOut @@ -522,7 +522,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G if (CS%debug) then call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, haloshift=0) endif - if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then + if (CS%use_kappa_shear .or. CS%use_cvmix_shear) then if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then call find_uv_at_h(u, v, h_orig, u_h, v_h, G, GV, eaml, ebml) if (CS%debug) then @@ -1920,7 +1920,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, apply parameterization of double-diffusion.", & default=.false. ) CS%use_kappa_shear = kappa_shear_is_used(param_file) - CS%use_CVMix_shear = cvmix_shear_is_used(param_file) + CS%use_cvmix_shear = cvmix_shear_is_used(param_file) if (CS%bulkmixedlayer) then call get_param(param_file, mod, "ML_MIX_FIRST", CS%ML_mix_first, & "The fraction of the mixed layer mixing that is applied \n"//& diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 2af9b8ef79..345f602902 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -20,6 +20,7 @@ module MOM_set_diffusivity use MOM_io, only : slasher, vardesc, var_desc, MOM_read_data use MOM_kappa_shear, only : calculate_kappa_shear, kappa_shear_init, Kappa_shear_CS use MOM_cvmix_shear, only : calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_cs +use MOM_cvmix_shear, only : cvmix_shear_end use MOM_string_functions, only : uppercase use MOM_thickness_diffuse, only : vert_fill_TS use MOM_variables, only : thermo_var_ptrs, vertvisc_type, p3d @@ -3148,12 +3149,16 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp end subroutine set_diffusivity_init +!> Clear pointers and dealocate memory subroutine set_diffusivity_end(CS) - type(set_diffusivity_CS), pointer :: CS + type(set_diffusivity_CS), pointer :: CS !< Control structure for this module if (CS%user_change_diff) & call user_change_diff_end(CS%user_change_diff_CSp) + if (CS%use_cvmix_shear) & + call cvmix_shear_end(CS%cvmix_shear_csp) + if (associated(CS)) deallocate(CS) end subroutine set_diffusivity_end From a77d1143b6e848f541dbdcb129519ed70b80f3f5 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 13 Mar 2018 14:12:54 -0600 Subject: [PATCH 11/40] Fixed a bug in the deallocate call --- src/parameterizations/vertical/MOM_cvmix_conv.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index 5fab297910..792f0deb1b 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -236,9 +236,9 @@ end subroutine calculate_cvmix_conv subroutine cvmix_conv_end(CS) type(cvmix_conv_cs), pointer :: CS ! Control structure - deallocate(CS%N2, CS%diag) - deallocate(CS%kd_conv_3d, CS%diag) - deallocate(CS%kv_conv_3d, CS%diag) + deallocate(CS%N2) + deallocate(CS%kd_conv_3d) + deallocate(CS%kv_conv_3d) deallocate(CS) end subroutine cvmix_conv_end From afd8f591c562667c19e67b70022b988f6c3659ca Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 13 Mar 2018 14:35:14 -0600 Subject: [PATCH 12/40] Fixed a bug in the deallocate call --- .../vertical/MOM_cvmix_shear.F90 | 35 ++++++++++++++++--- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_shear.F90 b/src/parameterizations/vertical/MOM_cvmix_shear.F90 index 487b1cbf6f..4a91524c76 100644 --- a/src/parameterizations/vertical/MOM_cvmix_shear.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_shear.F90 @@ -35,14 +35,21 @@ module MOM_cvmix_shear real :: KPP_exp !< real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) real, allocatable, dimension(:,:,:) :: S2 !< Squared shear frequency (1/s2) + real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number + real, allocatable, dimension(:,:,:) :: km !< vertical viscosity at interface (m2/s) + real, allocatable, dimension(:,:,:) :: kh !< vertical diffusivity at interface (m2/s) character(10) :: Mix_Scheme !< Mixing scheme name (string) + ! Daignostic handles and pointers + type(diag_ctrl), pointer :: diag => NULL() + integer :: id_N2 = -1, id_S2 = -1, id_ri_grad = -1, id_km = -1, id_kh = -1 + end type cvmix_shear_cs character(len=40) :: mdl = "MOM_cvmix_shear" !< This module's name. contains -!> Subroutine for calculating (internal) diffusivity +!> Subroutine for calculating (internal) vertical diffusivities/viscosities subroutine calculate_cvmix_shear(u_H, v_H, h, tv, KH, & KM, G, GV, CS ) type(ocean_grid_type), intent(in) :: G !< Grid structure. @@ -51,7 +58,7 @@ subroutine calculate_cvmix_shear(u_H, v_H, h, tv, KH, & real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: v_H !< Initial meridional velocity on T points, in m s-1. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: KH !< The vertical viscosity at each interface + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: KH !< The vertical diffusivity at each interface !! (not layer!) in m2 s-1. real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: KM !< The vertical viscosity at each interface !! (not layer!) in m2 s-1. @@ -196,6 +203,23 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) ! Allocation and initialization allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%N2(:,:,:) = 0. allocate( CS%S2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%S2(:,:,:) = 0. + !Initialize w/ large Richardson value + allocate( CS%ri_grad( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%ri_grad(:,:,:) = 1.e8 + + ! Register diagnostics + CS%diag => diag + CS%id_N2 = register_diag_field('ocean_model', 'shear_N2', diag%axesTi, Time, & + 'Square of Brunt-Vaisala frequency used by MOM_cvmix_shear module', '1/s2') + CS%id_S2 = register_diag_field('ocean_model', 'shear_S2', diag%axesTi, Time, & + 'Square of vertical shear used by MOM_cvmix_shear module','1/s2') + CS%id_ri_grad = register_diag_field('ocean_model', 'shear_ri_grad', diag%axesTi, Time, & + 'Gradient Richarson number used by MOM_cvmix_shear module','nondim') + CS%id_kh = register_diag_field('ocean_model', 'shear_KH', diag%axesTi, Time, & + 'Vertical diffusivity added by MOM_cvmix_shear module', 'm2/s') + if (CS%id_kh > 0) allocate(CS%kh(SZI_(G), SZJ_(G), SZK_(G)+1)) + CS%id_km = register_diag_field('ocean_model', 'shear_KM', diag%axesTi, Time, & + 'Vertical viscosity added by MOM_cvmix_shear module', 'm2/s') + if (CS%id_km > 0) allocate(CS%km(SZI_(G), SZJ_(G), SZK_(G)+1)) end function cvmix_shear_init @@ -217,8 +241,11 @@ end function cvmix_shear_is_used subroutine cvmix_shear_end(CS) type(cvmix_shear_cs), pointer :: CS ! Control structure - deallocate(CS%N2, CS%diag) - deallocate(CS%S2, CS%diag) + deallocate(CS%N2) + deallocate(CS%S2) + deallocate(CS%ri_grad) + if (CS%id_kh > 0) deallocate(CS%kh, CS%diag) + if (CS%id_km > 0) deallocate(CS%km, CS%diag) deallocate(CS) end subroutine cvmix_shear_end From 1d850a60845c582e78a4899a3179f7e54bcb6fe4 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 14 Mar 2018 16:56:05 -0600 Subject: [PATCH 13/40] Changed hbl to an optional pointer; do not apply mixing due to convection within the boundary layer --- .../vertical/MOM_cvmix_conv.F90 | 28 +++++++++++++++++-- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index 792f0deb1b..49ec101b44 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -60,6 +60,7 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) ! Local variables real :: prandtl_turb + logical :: useEPBL ! This include declares and sets the variable "version". #include "version_variable.h" @@ -83,6 +84,17 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) if (.not. cvmix_conv_init) return + call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useEPBL, default=.false., & + do_not_log=.true.) + + ! Warn user if EPBL is being used, since in this case mixing due to convection will + ! be aplied in the boundary layer + if (useEPBL) then + call MOM_error(WARNING, 'MOM_cvmix_conv_init: '// & + 'CVMix convection may not be properly applied when ENERGETICS_SFC_PBL = True'//& + 'as convective mixing might occur in the boundary layer.') + endif + call get_param(param_file, mdl, "PRANDTL_TURB", Prandtl_turb, & "The turbulent Prandtl number applied to shear/conv. \n"//& "instabilities.", units="nondim", default=1.0, do_not_log=.true.) @@ -131,16 +143,15 @@ end function cvmix_conv_init !> Subroutine for calculating enhanced diffusivity/viscosity !! due to convection via CVMix -subroutine calculate_cvmix_conv(h, tv, G, GV, hbl, CS) +subroutine calculate_cvmix_conv(h, tv, G, GV, CS, hbl) type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hbl!< Depth of ocean boundary layer (m) - !type(cvmix_conv_cs), intent(inout) :: CS !< The control structure returned by a previous call to type(cvmix_conv_cs), pointer :: CS !< The control structure returned by a previous call to !! CVMix_conv_init. + real, dimension(:,:), optional, pointer :: hbl!< Depth of ocean boundary layer (m) ! local variables real, dimension(SZK_(G)) :: rho_lwr !< Adiabatic Water Density, this is a dummy @@ -159,6 +170,11 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, hbl, CS) ! initialize dummy variables rho_lwr(:) = 0.0; rho_1d(:) = 0.0 + if (.not. associated(hbl)) then + allocate(hbl(SZI_(G), SZJ_(G))); + hbl(:,:) = 0.0 + endif + do j = G%jsc, G%jec do i = G%isc, G%iec @@ -206,6 +222,12 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, hbl, CS) max_nlev=G%ke, & OBL_ind=kOBL) + ! Do not apply mixing due to convection within the boundary layer + do k=1,NINT(hbl(i,j)) + CS%kv_conv_3d(i,j,k) = 0.0 + CS%kd_conv_3d(i,j,k) = 0.0 + enddo + enddo enddo From ba2a0a4a547e490f6f3d01f5937811b8ed4c088b Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 15 Mar 2018 11:17:45 -0600 Subject: [PATCH 14/40] Added function cvmix_conv_is_used and register restart fields. This function allows other modules to know whether this parameterization will be used without needing to duplicate the log entry. This commit also registers restart fields when only USE_CVMIX_CONVECTION is used (i.e., all other parame- terizations associated with vertical mixing are not enabled). --- .../vertical/MOM_cvmix_conv.F90 | 20 +++++++++---------- .../vertical/MOM_set_viscosity.F90 | 14 +++++++------ 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index 49ec101b44..88dece60fb 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -20,7 +20,7 @@ module MOM_cvmix_conv #include -public cvmix_conv_init, calculate_cvmix_conv, cvmix_conv_end +public cvmix_conv_init, calculate_cvmix_conv, cvmix_conv_end, cvmix_conv_is_used !> Control structure including parameters for CVMix convection. type, public :: cvmix_conv_cs @@ -244,15 +244,15 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, CS, hbl) end subroutine calculate_cvmix_conv -! GMM, not sure if we need the code below - DELETE???? -!!logical function cvmix_conv_is_used(param_file) -! Reads the parameter "USE_CVMIX_CONVECTION" and returns state. -! This function allows other modules to know whether this parameterization will -! be used without needing to duplicate the log entry. -!! type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters -!! call get_param(param_file, mdl, "USE_CVMIX_CONVECTION", kappa_shear_is_used, & -!! default=.false., do_not_log = .true.) -!!end function cvmix_conv_is_used +!> Reads the parameter "USE_CVMIX_CONVECTION" and returns state. +!! This function allows other modules to know whether this parameterization will +!! be used without needing to duplicate the log entry. +logical function cvmix_conv_is_used(param_file) + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + call get_param(param_file, mdl, "USE_CVMIX_CONVECTION", cvmix_conv_is_used, & + default=.false., do_not_log = .true.) + +end function cvmix_conv_is_used !> Clear pointers and dealocate memory subroutine cvmix_conv_end(CS) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 394b17dbd2..60c96fbb37 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -44,7 +44,8 @@ module MOM_set_visc use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_kappa_shear, only : kappa_shear_is_used -use MOM_CVMix_shear, only : CVMix_shear_is_used +use MOM_cvmix_shear, only : cvmix_shear_is_used +use MOM_cvmix_conv, only : cvmix_conv_is_used use MOM_io, only : vardesc, var_desc use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_variables, only : thermo_var_ptrs @@ -1783,18 +1784,19 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) ! (in) restart_CS - A pointer to the restart control structure. type(vardesc) :: vd logical :: use_kappa_shear, adiabatic, useKPP, useEPBL - logical :: use_CVMix, MLE_use_PBL_MLD + logical :: use_cvmix_shear, MLE_use_PBL_MLD, use_cvmix_conv integer :: isd, ied, jsd, jed, nz character(len=40) :: mdl = "MOM_set_visc" ! This module's name. isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., & do_not_log=.true.) - use_kappa_shear = .false. ; use_CVMix = .false. ; - useKPP = .false. ; useEPBL = .false. + use_kappa_shear = .false. ; use_cvmix_shear = .false. ; + useKPP = .false. ; useEPBL = .false. ; use_cvmix_conv = .false. ; if (.not.adiabatic) then use_kappa_shear = kappa_shear_is_used(param_file) - use_CVMix = CVMix_shear_is_used(param_file) + use_cvmix_shear = cvmix_shear_is_used(param_file) + use_cvmix_conv = cvmix_conv_is_used(param_file) call get_param(param_file, mdl, "USE_KPP", useKPP, & "If true, turns on the [CVmix] KPP scheme of Large et al., 1984,\n"// & "to calculate diffusivities and non-local transport in the OBL.", & @@ -1805,7 +1807,7 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) "in the surface boundary layer.", default=.false., do_not_log=.true.) endif - if (use_kappa_shear .or. useKPP .or. useEPBL .or. use_CVMix) then + if (use_kappa_shear .or. useKPP .or. useEPBL .or. use_cvmix_shear .or. use_cvmix_conv) then allocate(visc%Kd_turb(isd:ied,jsd:jed,nz+1)) ; visc%Kd_turb(:,:,:) = 0.0 allocate(visc%TKE_turb(isd:ied,jsd:jed,nz+1)) ; visc%TKE_turb(:,:,:) = 0.0 allocate(visc%Kv_turb(isd:ied,jsd:jed,nz+1)) ; visc%Kv_turb(:,:,:) = 0.0 From 941bfb6e29f3525e81c9ba081286bd138bf2ced2 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 15 Mar 2018 11:22:47 -0600 Subject: [PATCH 15/40] Updated mixing coeff. due to convection Vertical diff. is updated via Kd_int, but I could not find an equivalent for viscosity. I am also updating visc%Kv_turb and visc%Kv_turb, but I am not sure if this is correct. Further checking is needed! --- .../vertical/MOM_diabatic_driver.F90 | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index aa4f9f072e..cc0e9e7501 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -635,14 +635,11 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_turb, CS%KPP_NLTheat, CS%KPP_NLTscalar) !$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat) - ! If visc%MLD exists, copy the KPP BLD into it - if (associated(visc%MLD)) then - call KPP_get_BLD(CS%KPP_CSp, visc%MLD(:,:), G) - call pass_var(visc%MLD, G%domain, halo=1) - Hml(:,:) = visc%MLD(:,:) + 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 @@ -676,12 +673,16 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G ! Add diffusivity due to convection (computed via CVMix) if (CS%use_cvmix_conv) then - call calculate_cvmix_conv(h, tv, G, GV, Hml, CS%cvmix_conv_csp) + 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 - visc%Kd_turb(i,j,k) = visc%Kd_turb(i,j,k) + CS%cvmix_conv_csp%kd_conv_3d(i,j,k) + Kd_int(i,j,k) = Kd_int(i,j,k) + CS%cvmix_conv_csp%kd_conv_3d(i,j,k) + ! GMM, I am not sure if Kv_turb is the right place to add kv_conv_3d visc%Kv_turb(i,j,k) = visc%Kv_turb(i,j,k) + CS%cvmix_conv_csp%kv_conv_3d(i,j,k) + visc%Kd_turb(i,j,k) = visc%Kd_turb(i,j,k) + CS%cvmix_conv_csp%kd_conv_3d(i,j,k) enddo ; enddo ; enddo + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! endif From 53004ec665089d2fab707ced506d7501da62ffd0 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 19 Mar 2018 16:35:05 -0600 Subject: [PATCH 16/40] Deleted module MOM_diffConvection This module was not being used and convection is now applied via MOM_cvmix_conv. --- .../vertical/MOM_diffConvection.F90 | 163 ------------------ 1 file changed, 163 deletions(-) delete mode 100644 src/parameterizations/vertical/MOM_diffConvection.F90 diff --git a/src/parameterizations/vertical/MOM_diffConvection.F90 b/src/parameterizations/vertical/MOM_diffConvection.F90 deleted file mode 100644 index b63dbe4472..0000000000 --- a/src/parameterizations/vertical/MOM_diffConvection.F90 +++ /dev/null @@ -1,163 +0,0 @@ -module MOM_diffConvection - -! This file is part of MOM6. See LICENSE.md for the license. - -use MOM_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data -use MOM_diag_mediator, only : query_averaging_enabled, register_diag_field -use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_PE -use MOM_EOS, only : EOS_type, calculate_density -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_file_parser, only : openParameterBlock, closeParameterBlock -use MOM_grid, only : ocean_grid_type, isPointInCell -use MOM_verticalGrid, only : verticalGrid_type - -implicit none ; private - -#include "MOM_memory.h" - -public :: diffConvection_init, diffConvection_calculate, diffConvection_end - -! Control structure for containing KPP parameters/data -type, public :: diffConvection_CS ; private - - ! Parameters - real :: Kd_convection ! The value of diffusivity to add at statically unstable interfaces (m2/s) - logical :: debug ! If true, turn on debugging - logical :: passiveMode ! If true, make the motions but go nowhere - - ! Daignostic handles and pointers - type(diag_ctrl), pointer :: diag => NULL() - integer :: id_N2 = -1, id_Kd_conv = -1 - - ! Diagnostics arrays - real, allocatable, dimension(:,:,:) :: N2 ! Brunt-Vaisala frequency (1/s2) - real, allocatable, dimension(:,:,:) :: Kd_conv ! Diffusivity added by convection (m2/s) - -end type diffConvection_CS - -! Module data used for debugging only -logical, parameter :: verbose = .False. - -contains - -logical function diffConvection_init(paramFile, G, diag, Time, CS) -!< Initialize the CVmix KPP module and set up diagnostics -!! Returns True if module is to be used, otherwise returns False. - -! Arguments - type(param_file_type), intent(in) :: paramFile !< File parser - type(ocean_grid_type), intent(in) :: G !< Ocean grid - type(diag_ctrl), target, intent(in) :: diag !< Diagnostics - type(time_type), intent(in) :: Time !< Time - type(diffConvection_CS), pointer :: CS !< Control structure -! Local variables -#include "version_variable.h" - character(len=40) :: mdl = 'MOM_diffConvection' ! This module's name. - - if (associated(CS)) call MOM_error(FATAL, 'MOM_diffConvection, diffConvection_init: '// & - 'Control structure has already been initialized') - allocate(CS) - -! Read parameters - call log_version(paramFile, mdl, version, & - 'This module implements enhanced diffusivity as a\n' // & - 'function of static stability, N^2.') - call get_param(paramFile, mdl, "USE_CONVECTION", diffConvection_init, & - "If true, turns on the diffusive convection scheme that\n"// & - "increases diapycnal diffusivities at statically unstable\n"// & - "interfaces. Relevant parameters are contained in the\n"// & - "CONVECTION% parameter block.", & - default=.false.) - - call openParameterBlock(paramFile,'CONVECTION') - call get_param(paramFile, mdl, 'PASSIVE', CS%passiveMode, & - 'If True, puts KPP into a passive-diagnostic mode.', & - default=.False.) - call get_param(paramFile, mdl, 'KD_CONV', CS%Kd_convection, & - 'DIffusivity used in statically unstable regions of column.', & - units='m2/s', default=1.00) - call closeParameterBlock(paramFile) - call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) - -! Forego remainder of initialization if not using this scheme - if (.not. diffConvection_init) return - -! Register diagnostics - CS%diag => diag - CS%id_N2 = register_diag_field('ocean_model', 'Conv_N2', diag%axesTi, Time, & - 'Square of Brunt-Vaisala frequency used by diffConvection module', '1/s2') - if (CS%id_N2 > 0) allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) - CS%id_Kd_conv = register_diag_field('ocean_model', 'Conv_Kd', diag%axesTi, Time, & - 'Additional diffusivity added by diffConvection module', 'm2/s') - if (CS%id_Kd_conv > 0) allocate( CS%Kd_conv( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) - - if (CS%id_N2 > 0) CS%N2(:,:,:) = 0. - if (CS%id_Kd_conv > 0) CS%Kd_conv(:,:,:) = 0. - -end function diffConvection_init - - -subroutine diffConvection_calculate(CS, G, GV, h, Temp, Salt, EOS, Kd_int) -!< Calculates diffusivity and non-local transport for KPP parameterization - -! Arguments - type(diffConvection_CS), pointer :: CS !< Control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses (units of H) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Temp !< Pot. temperature (degrees C) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Salt !< Salinity (ppt) - type(EOS_type), pointer :: EOS !< Equation of state - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kd_int !< (in) Vertical diffusivity on interfaces (m2/s) - !! (out) Modified vertical diffusivity (m2/s) -! Local variables - integer :: i, j, k - real, dimension( G%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces (1/s2) - real, dimension( G%ke+1 ) :: Kd_1d ! Vertical diffusivity at interfaces (m2/s) - real :: GoRho, pRef, rhoK, rhoKm1 - - GoRho = GV%g_Earth / GV%Rho0 - - N2_1d( 1 ) = 0. - N2_1d( G%ke+1 ) = 0. - Kd_1d( 1 ) = 0. - Kd_1d( G%ke+1 ) = 0. - do j=G%jsc,G%jec ; do i=G%isc,G%iec - ! This k-loop calculates external quantities independent of any iterations - ! Start at bottom of top level - pRef = 0. ! Ignore atmospheric pressure - do K = 2, G%ke - ! Pressure at interface K is incremented by mass of level above - pRef = pRef + GV%g_Earth * GV%Rho0 * h(i,j,k-1) * GV%H_to_m ! Boussinesq approximation!!!! ????? - ! Compute Brunt-Vaisala frequency (static stability) on interfaces - call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS) - call calculate_density(Temp(i,j,k-1), Salt(i,j,k-1), pRef, rhoKm1, EOS) - N2_1d(K) = GoRho * (rhoK - rhoKm1) / & - (0.5*(h(i,j,k-1) + h(i,j,k)) + GV%H_subroundoff) ! Can be negative - Kd_1d(K) = 0. - if (N2_1d(K) < 0.) Kd_1d(K) = CS%Kd_convection - enddo ! k - - if (.not. CS%passiveMode) Kd_int(i,j,:) = Kd_int(i,j,:) + Kd_1d(:) - - if (CS%id_N2 > 0) CS%N2(i,j,:) = N2_1d(:) - if (CS%id_Kd_conv > 0) CS%Kd_conv(i,j,:) = Kd_1d(:) - - enddo ; enddo ! j - - if (CS%id_N2 > 0) call post_data(CS%id_N2, CS%N2, CS%diag) - if (CS%id_Kd_conv > 0) call post_data(CS%id_Kd_conv, CS%Kd_conv, CS%diag) - -end subroutine diffConvection_calculate - - -subroutine diffConvection_end(CS) -! Clear pointers, dealocate memory - type(diffConvection_CS), pointer :: CS ! Control structure - - if (CS%id_N2 > 0) deallocate(CS%N2, CS%diag) - if (CS%id_Kd_conv > 0) deallocate(CS%Kd_conv, CS%diag) - deallocate(CS) -end subroutine diffConvection_end - -end module MOM_diffConvection From a101e49ae3a0b319577fc4aa159832726103559f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 20 Mar 2018 16:02:21 -0600 Subject: [PATCH 17/40] Reduced the numbers of characters in a line --- src/parameterizations/vertical/MOM_cvmix_conv.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index 88dece60fb..721c35faa2 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -149,8 +149,8 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, CS, hbl) type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. - type(cvmix_conv_cs), pointer :: CS !< The control structure returned by a previous call to - !! CVMix_conv_init. + type(cvmix_conv_cs), pointer :: CS !< The control structure returned + !! by a previous call to CVMix_conv_init. real, dimension(:,:), optional, pointer :: hbl!< Depth of ocean boundary layer (m) ! local variables From 83089ff18e4d2de1240fdf48fb4926c9cff49c63 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 21 Mar 2018 10:44:30 -0600 Subject: [PATCH 18/40] Updated CVMix --- pkg/CVMix-src | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pkg/CVMix-src b/pkg/CVMix-src index d83f582714..66857c94be 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit d83f582714e7f0f98d20efd8fac8fab01fa3bfe6 +Subproject commit 66857c94bed214c32ccb5791010c6611ac5ae270 From 1116f49ce2192620e3b303e27e4bc8f08090cc42 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 21 Mar 2018 18:10:35 -0600 Subject: [PATCH 19/40] Fixed a problem in cvmix_init_bkgnd_BryanLewis_low --- pkg/CVMix-src | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pkg/CVMix-src b/pkg/CVMix-src index 66857c94be..c3a711e4e4 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit 66857c94bed214c32ccb5791010c6611ac5ae270 +Subproject commit c3a711e4e45f5ebcdc528f8ac690ad9a843375e9 From 1be0248846dfaa38c1152fb870ba618f6ef48232 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 21 Mar 2018 18:29:31 -0600 Subject: [PATCH 20/40] Added first version of MOM_bkgnd_mixing * everything related to setting background mixing that was previously in MOM_set_diffusivity, has been moved to this new module * Bryan and Lewis background mixing is now applied via CVMix --- .../vertical/MOM_bkgnd_mixing.F90 | 603 ++++++++++++++++++ 1 file changed, 603 insertions(+) create mode 100644 src/parameterizations/vertical/MOM_bkgnd_mixing.F90 diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 new file mode 100644 index 0000000000..3ad9dfd638 --- /dev/null +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -0,0 +1,603 @@ +!> Interface to background mixing schemes, including the Bryan and Lewis (1979) +!! which is applied via CVMix. + +module MOM_bkgnd_mixing + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field +use MOM_diag_mediator, only : post_data +use MOM_EOS, only : calculate_density, calculate_density_derivs +use MOM_variables, only : thermo_var_ptrs +use MOM_forcing_type, only : forcing +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_debugging, only : hchksum +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_file_parser, only : get_param, log_version, param_file_type +use cvmix_background, only : cvmix_init_bkgnd, cvmix_coeffs_bkgnd +use MOM_intrinsic_functions, only : invcosh + +implicit none ; private + +#include + +public bkgnd_mixing_init, bkgnd_mixing_end, calculate_bkgnd_mixing + +!> Control structure including parameters for this module. +type, public :: bkgnd_mixing_cs + + ! Parameters + real :: Kd_Bryan_Lewis_deep !< The abyssal value of a Bryan-Lewis diffusivity profile + !! (m2/s) + real :: Kd_Bryan_Lewis_surface !< "The surface value of a Bryan-Lewis diffusivity profile + !! (m2/s) + real :: Bryan_Lewis_depth_cent !< The depth about which the transition in the Bryan-Lewis + !! is centered (m) + real :: Bryan_Lewis_width_trans!< The width of the transition in the Bryan-Lewis profile (m) + real :: Kd_min !< minimum diapycnal diffusivity (m2/s) + real :: Kd !< interior diapycnal diffusivity (m2/s) + real :: N0_2Omega !< ratio of the typical Buoyancy frequency to + !! twice the Earth's rotation period, used with the + !! Henyey scaling from the mixing + real :: prandtl_turb !< Turbulent Prandtl number + real :: Kd_tanh_lat_scale !< A nondimensional scaling for the range of + !! diffusivities with Kd_tanh_lat_fn. Valid values + !! are in the range of -2 to 2; 0.4 reproduces CM2M. + real :: Kdml !< mixed layer diapycnal diffusivity (m2/s) + !! when bulkmixedlayer==.false. + real :: Hmix !< mixed layer thickness (meter) when + !! bulkmixedlayer==.false. + logical :: Kd_tanh_lat_fn !< If true, use the tanh dependence of Kd_sfc on + !! latitude, like GFDL CM2.1/CM2M. There is no + !! physical justification for this form, and it can + !! not be used with Henyey_IGW_background. + logical :: Bryan_Lewis_diffusivity!< If true, background vertical diffusivity + !! uses Bryan-Lewis (1979) like tanh profile. + logical :: Henyey_IGW_background !< If true, use a simplified variant of the + !! Henyey et al, JGR (1986) latitudinal scaling for the background diapycnal diffusivity, + !! which gives a marked decrease in the diffusivity near the equator. The simplification + !! here is to assume that the in-situ stratification is the same as the reference stratificaiton. + logical :: Henyey_IGW_background_new !< same as Henyey_IGW_background + !! but incorporate the effect of stratification on TKE dissipation, + !! e = f/f_0 * acosh(N/f) / acosh(N_0/f_0) * e_0 + !! where e is the TKE dissipation, and N_0 and f_0 + !! are the reference buoyancy frequency and inertial frequencies respectively. + !! e_0 is the reference dissipation at (N_0,f_0). In the previous version, N=N_0. + !! Additionally, the squared inverse relationship between diapycnal diffusivities + !! and stratification is included: + !! + !! kd = e/N^2 + !! + !! where kd is the diapycnal diffusivity. This approach assumes that work done + !! against gravity is uniformly distributed throughout the column. Whereas, kd=kd_0*e, + !! as in the original version, concentrates buoyancy work in regions of strong stratification. + logical :: bulkmixedlayer !< If true, a refined bulk mixed layer scheme is used + logical :: debug !< If true, turn on debugging in this module + ! Daignostic handles and pointers + type(diag_ctrl), pointer :: diag => NULL() + integer :: id_kd_bkgnd = -1, id_kv_bkgnd = -1 + + ! Diagnostics arrays + real, allocatable, dimension(:,:,:) :: kd_bkgnd !< Background diffusivity (m2/s) + real, allocatable, dimension(:,:,:) :: kv_bkgnd !< Background viscosity (m2/s) + +end type bkgnd_mixing_cs + +character(len=40) :: mdl = "MOM_bkgnd_mixing" !< This module's name. + +contains + +!> Initialize the background mixing routine. +subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) + + type(time_type), intent(in) :: Time !< The current time. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. + type(bkgnd_mixing_cs), pointer :: CS !< This module's control structure. + + ! Local variables + +! This include declares and sets the variable "version". +#include "version_variable.h" + + if (associated(CS)) then + call MOM_error(WARNING, "bkgnd_mixing_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + ! Read parameters + call log_version(param_file, mdl, version, & + "Adding static vertical background mixing coefficients") + + call get_param(param_file, mdl, "KD", CS%Kd, & + "The background diapycnal diffusivity of density in the \n"//& + "interior. Zero or the molecular value, ~1e-7 m2 s-1, \n"//& + "may be used.", units="m2 s-1", fail_if_missing=.true.) + + call get_param(param_file, mdl, "KD_MIN", CS%Kd_min, & + "The minimum diapycnal diffusivity.", & + units="m2 s-1", default=0.01*CS%Kd) + + ! The following is needed to set one of the choices of vertical background mixing + call get_param(param_file, mdl, "BULKMIXEDLAYER", CS%bulkmixedlayer, & + do_not_log=.true.) + if (CS%bulkmixedlayer) then + ! Check that Kdml is not set when using bulk mixed layer + call get_param(param_file, mdl, "KDML", CS%Kdml, default=-1.) + if (CS%Kdml>0.) call MOM_error(FATAL, & + "bkgnd_mixing_init: KDML cannot be set when using"// & + "bulk mixed layer.") + CS%Kdml = CS%Kd ! This is not used with a bulk mixed layer, but also + ! cannot be a NaN. + else + call get_param(param_file, mdl, "KDML", CS%Kdml, & + "If BULKMIXEDLAYER is false, KDML is the elevated \n"//& + "diapycnal diffusivity in the topmost HMIX of fluid. \n"//& + "KDML is only used if BULKMIXEDLAYER is false.", & + units="m2 s-1", default=CS%Kd) + call get_param(param_file, mdl, "HMIX_FIXED", CS%Hmix, & + "The prescribed depth over which the near-surface \n"//& + "viscosity and diffusivity are elevated when the bulk \n"//& + "mixed layer is not used.", units="m", fail_if_missing=.true.) + endif + + + call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) + + call get_param(param_file, mdl, "PRANDTL_TURB", CS%prandtl_turb, & + units="nondim", default=1.0, do_not_log=.true.) + + call openParameterBlock(param_file,'MOM_BACKGROUND_MIXING') + + call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", & + CS%Bryan_Lewis_diffusivity, & + "If true, use a Bryan & Lewis (JGR 1979) like tanh \n"//& + "profile of background diapycnal diffusivity with depth. \n"//& + "This is done via CVMix.", default=.false.) + + if (CS%Bryan_Lewis_diffusivity) then + + call get_param(param_file, mdl, "KD_BRYAN_LEWIS_DEEP", & + CS%Kd_Bryan_Lewis_deep, & + "The abyssal value of a Bryan-Lewis diffusivity profile.", & + units="m2 s-1", fail_if_missing=.true.) + + call get_param(param_file, mdl, "KD_BRYAN_LEWIS_SURFACE", & + CS%Kd_Bryan_Lewis_surface, & + "The surface value of a Bryan-Lewis diffusivity profile.", & + units="m2 s-1", fail_if_missing=.true.) + + call get_param(param_file, mdl, "BRYAN_LEWIS_DEPTH_CENT", & + CS%Bryan_Lewis_depth_cent, & + "The depth about which the transition in the Bryan-Lewis.", & + units="m", fail_if_missing=.true.) + + call get_param(param_file, mdl, "BRYAN_LEWIS_WIDTH_TRANS", & + CS%Bryan_Lewis_width_trans, & + "The width of the transition in the Bryan-Lewis.",& + units="m", fail_if_missing=.true.) + + endif ! CS%Bryan_Lewis_diffusivity + + call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND", & + CS%Henyey_IGW_background, & + "If true, use a latitude-dependent scaling for the near \n"//& + "surface background diffusivity, as described in \n"//& + "Harrison & Hallberg, JPO 2008.", default=.false.) + + call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND_NEW", & + CS%Henyey_IGW_background_new, & + "If true, use a better latitude-dependent scaling for the\n"//& + "background diffusivity, as described in \n"//& + "Harrison & Hallberg, JPO 2008.", default=.false.) + + if (CS%Henyey_IGW_background .and. CS%Henyey_IGW_background_new) & + call MOM_error(FATAL, "set_diffusivity_init: HENYEY_IGW_BACKGROUND and \n"//& + "HENYEY_IGW_BACKGROUND_NEW are mutually exclusive. Set only one or none.") + + if (CS%Henyey_IGW_background) & + call get_param(param_file, mdl, "HENYEY_N0_2OMEGA", CS%N0_2Omega, & + "The ratio of the typical Buoyancy frequency to twice \n"//& + "the Earth's rotation period, used with the Henyey \n"//& + "scaling from the mixing.", units="nondim", default=20.0) + + call get_param(param_file, mdl, "KD_TANH_LAT_FN", & + CS%Kd_tanh_lat_fn, & + "If true, use a tanh dependence of Kd_sfc on latitude, \n"//& + "like CM2.1/CM2M. There is no physical justification \n"//& + "for this form, and it can not be used with \n"//& + "HENYEY_IGW_BACKGROUND.", default=.false.) + + if (CS%Kd_tanh_lat_fn) & + call get_param(param_file, mdl, "KD_TANH_LAT_SCALE", & + CS%Kd_tanh_lat_scale, & + "A nondimensional scaling for the range ofdiffusivities \n"//& + "with KD_TANH_LAT_FN. Valid values are in the range of \n"//& + "-2 to 2; 0.4 reproduces CM2M.", units="nondim", default=0.0) + + if (CS%Henyey_IGW_background .and. CS%Kd_tanh_lat_fn) call MOM_error(FATAL, & + "MOM_bkgnd_mixing: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.") + + call closeParameterBlock(param_file) + + ! allocate arrays and set them to zero + allocate(CS%kd_bkgnd(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kd_bkgnd(:,:,:) = 0. + allocate(CS%kv_bkgnd(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kv_bkgnd(:,:,:) = 0. + + ! Register diagnostics + CS%diag => diag + CS%id_kd_bkgnd = register_diag_field('ocean_model', 'bkgnd_kd', diag%axesTi, Time, & + 'Background diffusivity added by MOM_bkgnd_mixing module', 'm2/s') + CS%id_kv_bkgnd = register_diag_field('ocean_model', 'bkgnd_kv', diag%axesTi, Time, & + 'Background viscosity added by MOM_bkgnd_mixing module', 'm2/s') + +end subroutine bkgnd_mixing_init + +!> Subroutine for calculating vertical background diffusivities/viscosities +subroutine calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS) + + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T_f!< temperature (deg C), after massless + !! layers filled vertically by diffusion. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: S_f!< salinity, after massless + !! layers filled vertically by diffusion. + type(forcing), intent(in) :: fluxes !< Structure of surface fluxes that may be + !! used. + type(bkgnd_mixing_cs), pointer :: CS !< The control structure returned by + !! a previous call to bkgnd_mixing_init. + + ! local variables + real, dimension(SZI_(G), SZJ_(G)) :: Kd_sfc !< surface value of the diffusivity (m2/s) + real, dimension(SZI_(G)) :: & + depth_i, & !< distance from surface of an interface (meter) + N2_bot !< bottom squared buoyancy frequency (1/s2) + 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? + real, dimension(SZI_(G),SZK_(G)) :: & + N2_lay !< squared buoyancy frequency associated with layers (1/s2) + + real, dimension(SZK_(G)) :: depth_k !< distance from surface of an interface (meter) + real :: I_x30 !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg))) + real :: deg_to_rad !< factor converting degrees to radians, pi/180. + real :: abs_sin !< absolute value of sine of latitude (nondim) + real :: depth_c !< depth of the center of a layer (meter) + real :: I_Hmix !< inverse of fixed mixed layer thickness (1/m) + real :: epsilon + real :: I_2Omega !< 1/(2 Omega) (sec) + real :: N_2Omega + real :: N02_N2 + integer :: i, j, k, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + + ! set some parameters + deg_to_rad = atan(1.0)/45.0 ! = PI/180 + epsilon = 1.e-10 + + if (CS%Bryan_Lewis_diffusivity) then +!$OMP parallel do default(none) shared(is,ie,js,je,nz, CS) +!$OMP private(cvmix_init_bkgnd,cvmix_coeffs_bkgnd) + + ! Bryan & Lewis is computed via CVMix + do j=js,je; do i=is,ie + + depth_k(:) = 0.0 + do k=1,nz + depth_k(k) = depth_k(k) + GV%H_to_m*h(i,j,k) + enddo + + call cvmix_init_bkgnd(max_nlev=nz, & + zw = depth_k(:), & !< interface depth, must be positive. + bl1 = CS%Kd_Bryan_Lewis_deep, & + bl2 = CS%Kd_Bryan_Lewis_surface, & + bl3 = 1.0/CS%Bryan_Lewis_depth_cent , & + bl4 = CS%Bryan_Lewis_width_trans, & + prandtl = CS%prandtl_turb) + + call cvmix_coeffs_bkgnd(Mdiff_out=CS%kv_bkgnd(i,j,:), & + Tdiff_out=CS%kd_bkgnd(i,j,:), & + nlev=nz, & + max_nlev=nz) + enddo; enddo + else +!$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc) + do j=js,je ; do i=is,ie + Kd_sfc(i,j) = CS%Kd + enddo ; enddo + endif + + if (CS%Henyey_IGW_background) then + I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. +!$OMP parallel do default(none) +!shared(is,ie,js,je,Kd_sfc,CS,G,deg_to_rad,epsilon,I_x30) & +!$OMP private(abs_sin) + do j=js,je ; do i=is,ie + abs_sin = abs(sin(G%geoLatT(i,j)*deg_to_rad)) + Kd_sfc(i,j) = max(CS%Kd_min, Kd_sfc(i,j) * & + ((abs_sin * invcosh(CS%N0_2Omega/max(epsilon,abs_sin))) * I_x30) ) + enddo ; enddo + elseif (CS%Kd_tanh_lat_fn) then +!$OMP parallel do default(none) shared(is,ie,js,je,Kd_sfc,CS,G) + do j=js,je ; do i=is,ie + ! The transition latitude and latitude range are hard-scaled here, since + ! this is not really intended for wide-spread use, but rather for + ! comparison with CM2M / CM2.1 settings. + Kd_sfc(i,j) = max(CS%Kd_min, Kd_sfc(i,j) * (1.0 + & + CS%Kd_tanh_lat_scale * 0.5*tanh((abs(G%geoLatT(i,j)) - 35.0)/5.0) )) + enddo ; enddo + endif + +!$OMP parallel do default(none) +!shared(is,ie,js,je,nz,G,GV,CS,h,tv,T_f,S_f,fluxes,dd, & +!$OMP +!Kd,Kd_sfc,epsilon,deg_to_rad,I_2Omega,visc, & +!$OMP Kd_int,dt,u,v,Omega2) & +!$OMP +!private(dRho_int,I_trans,atan_fn_sfc,I_atan_fn,atan_fn_lay, & +!$OMP I_Hmix,depth_c,depth,N2_lay, N2_int, +!N2_bot, & +!$OMP I_x30,abs_sin,N_2Omega,N02_N2,KT_extra, +!KS_extra, & +!$OMP TKE_to_Kd,maxTKE,dissip,kb) + do j=js,je + ! Set up variables related to the stratification. + call find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, dRho_int, N2_lay, N2_int, N2_bot) + !if (associated(dd%N2_3d)) then + ! do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo + !endif + + ! Set up the background diffusivity. + if ((.not. CS%Bryan_Lewis_diffusivity) .and. (.not.CS%bulkmixedlayer) .and. & + (CS%Kd/= CS%Kdml)) then + + I_Hmix = 1.0 / CS%Hmix + do i=is,ie ; depth_i(i) = 0.0 ; enddo + do k=1,nz ; do i=is,ie + depth_c = depth_i(i) + 0.5*GV%H_to_m*h(i,j,k) + + if (depth_c <= CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = CS%Kdml + elseif (depth_c >= 2.0*CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = Kd_sfc(i,j) + else + CS%kd_bkgnd(i,j,k) = ((Kd_sfc(i,j) - CS%Kdml) * I_Hmix) * depth_c + & + (2.0*CS%Kdml - Kd_sfc(i,j)) + endif + + depth_i(i) = depth_i(i) + GV%H_to_m*h(i,j,k) + enddo ; enddo + elseif (CS%Henyey_IGW_background_new) then + I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. + do k=1,nz ; do i=is,ie + abs_sin = max(epsilon,abs(sin(G%geoLatT(i,j)*deg_to_rad))) + N_2Omega = max(abs_sin,sqrt(N2_lay(i,k))*I_2Omega) + N02_N2 = (CS%N0_2Omega/N_2Omega)**2 + CS%kd_bkgnd(i,j,k) = max(CS%Kd_min, Kd_sfc(i,j) * & + ((abs_sin * invcosh(N_2Omega/abs_sin)) * I_x30)*N02_N2) + enddo ; enddo + + else + do k=1,nz ; do i=is,ie + CS%kd_bkgnd(i,j,k) = Kd_sfc(i,j) + enddo ; enddo + endif + enddo ! j-loop + + if (CS%debug) then + call hchksum(Kd_sfc,"Kd_sfc",G%HI,haloshift=0) + call hchksum(CS%kd_bkgnd, "MOM_bkgnd_mixing: kd_bkgnd",G%HI,haloshift=0) + call hchksum(CS%kv_bkgnd, "MOM_bkgnd_mixing: kv_bkgnd",G%HI,haloshift=0) + endif + + ! send diagnostics to post_data + if (CS%id_kd_bkgnd > 0) call post_data(CS%id_kd_bkgnd, CS%kd_bkgnd, CS%diag) + if (CS%id_kv_bkgnd > 0) call post_data(CS%id_kv_bkgnd, CS%kv_bkgnd, CS%diag) + +end subroutine calculate_bkgnd_mixing + + +!> Computes N2 +subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, dRho_int, & + N2_lay, N2_int, N2_bot) + 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(in) :: tv + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T_f, S_f + type(forcing), intent(in) :: fluxes + integer, intent(in) :: j + real, dimension(SZI_(G),SZK_(G)+1), intent(out) :: dRho_int, N2_int + real, dimension(SZI_(G),SZK_(G)), intent(out) :: N2_lay + real, dimension(SZI_(G)), intent(out) :: N2_bot + + real, dimension(SZI_(G),SZK_(G)+1) :: & + dRho_int_unfilt, & ! unfiltered density differences across interfaces + dRho_dT, & ! partial derivative of density wrt temp (kg m-3 degC-1) + dRho_dS ! partial derivative of density wrt saln (kg m-3 PPT-1) + + real, dimension(SZI_(G)) :: & + pres, & ! pressure at each interface (Pa) + Temp_int, & ! temperature at each interface (degC) + Salin_int, & ! salinity at each interface (PPT) + drho_bot, & + h_amp, & + hb, & + z_from_bot + + real :: Rml_base ! density of the deepest variable density layer + real :: dz_int ! thickness associated with an interface (meter) + real :: G_Rho0 ! gravitation acceleration divided by Bouss reference density (m4 s-2 kg-1) + real :: H_neglect ! negligibly small thickness, in the same units as h. + + logical :: do_i(SZI_(G)), do_any + integer :: i, k, is, ie, nz + + is = G%isc ; ie = G%iec ; nz = G%ke + G_Rho0 = GV%g_Earth / GV%Rho0 + H_neglect = GV%H_subroundoff + + ! Find the (limited) density jump across each interface. + do i=is,ie + dRho_int(i,1) = 0.0 ; dRho_int(i,nz+1) = 0.0 + dRho_int_unfilt(i,1) = 0.0 ; dRho_int_unfilt(i,nz+1) = 0.0 + enddo + if (associated(tv%eqn_of_state)) then + if (associated(fluxes%p_surf)) then + do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo + else + do i=is,ie ; pres(i) = 0.0 ; enddo + endif + 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) + T_f(i,j,k-1)) + Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1)) + enddo + call calculate_density_derivs(Temp_int, Salin_int, pres, & + dRho_dT(:,K), dRho_dS(:,K), is, ie-is+1, tv%eqn_of_state) + do i=is,ie + dRho_int(i,K) = max(dRho_dT(i,K)*(T_f(i,j,k) - T_f(i,j,k-1)) + & + dRho_dS(i,K)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0) + dRho_int_unfilt(i,K) = max(dRho_dT(i,K)*(tv%T(i,j,k) - tv%T(i,j,k-1)) + & + dRho_dS(i,K)*(tv%S(i,j,k) - tv%S(i,j,k-1)), 0.0) + enddo + enddo + else + do K=2,nz ; do i=is,ie + dRho_int(i,K) = GV%Rlay(k) - GV%Rlay(k-1) + enddo ; enddo + endif + + ! Set the buoyancy frequencies. + do k=1,nz ; do i=is,ie + N2_lay(i,k) = G_Rho0 * 0.5*(dRho_int(i,K) + dRho_int(i,K+1)) / & + (GV%H_to_m*(h(i,j,k) + H_neglect)) + enddo ; enddo + do i=is,ie ; N2_int(i,1) = 0.0 ; N2_int(i,nz+1) = 0.0 ; enddo + do K=2,nz ; do i=is,ie + N2_int(i,K) = G_Rho0 * dRho_int(i,K) / & + (0.5*GV%H_to_m*(h(i,j,k-1) + h(i,j,k) + H_neglect)) + enddo ; enddo + + ! Find the bottom boundary layer stratification, and use this in the deepest layers. + do i=is,ie + hb(i) = 0.0 ; dRho_bot(i) = 0.0 + z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) + do_i(i) = (G%mask2dT(i,j) > 0.5) + h_amp(i) = 0.0 + enddo + + do k=nz,2,-1 + do_any = .false. + do i=is,ie ; if (do_i(i)) then + dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) + z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above + + hb(i) = hb(i) + dz_int + dRho_bot(i) = dRho_bot(i) + dRho_int(i,K) + + if (z_from_bot(i) > h_amp(i)) then + if (k>2) then + ! Always include at least one full layer. + hb(i) = hb(i) + 0.5*GV%H_to_m*(h(i,j,k-1) + h(i,j,k-2)) + dRho_bot(i) = dRho_bot(i) + dRho_int(i,K-1) + endif + do_i(i) = .false. + else + do_any = .true. + endif + endif ; enddo + if (.not.do_any) exit + enddo + + do i=is,ie + if (hb(i) > 0.0) then + N2_bot(i) = (G_Rho0 * dRho_bot(i)) / hb(i) + else ; N2_bot(i) = 0.0 ; endif + z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) + do_i(i) = (G%mask2dT(i,j) > 0.5) + enddo + + do k=nz,2,-1 + do_any = .false. + do i=is,ie ; if (do_i(i)) then + dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) + z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above + + N2_int(i,K) = N2_bot(i) + if (k>2) N2_lay(i,k-1) = N2_bot(i) + + if (z_from_bot(i) > h_amp(i)) then + if (k>2) N2_int(i,K-1) = N2_bot(i) + do_i(i) = .false. + else + do_any = .true. + endif + endif ; enddo + if (.not.do_any) exit + enddo + do i=is,ie + if (hb(i) > 0.0) then + N2_bot(i) = (G_Rho0 * dRho_bot(i)) / hb(i) + else ; N2_bot(i) = 0.0 ; endif + z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) + do_i(i) = (G%mask2dT(i,j) > 0.5) + enddo + + do k=nz,2,-1 + do_any = .false. + do i=is,ie ; if (do_i(i)) then + dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) + z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above + + N2_int(i,K) = N2_bot(i) + if (k>2) N2_lay(i,k-1) = N2_bot(i) + + if (z_from_bot(i) > h_amp(i)) then + if (k>2) N2_int(i,K-1) = N2_bot(i) + do_i(i) = .false. + else + do_any = .true. + endif + endif ; enddo + if (.not.do_any) exit + enddo + + if (associated(tv%eqn_of_state)) then + do K=1,nz+1 ; do i=is,ie + dRho_int(i,K) = dRho_int_unfilt(i,K) + enddo ; enddo + endif + +end subroutine find_N2 + +!> Reads the parameter "USE_CVMIX_BACKGROUND" and returns state. +!! This function allows other modules to know whether this parameterization will +!! be used without needing to duplicate the log entry. +logical function cvmix_bkgnd_is_used(param_file) + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + call get_param(param_file, mdl, "USE_CVMIX_BACKGROUND", cvmix_bkgnd_is_used, & + default=.false., do_not_log = .true.) + +end function cvmix_bkgnd_is_used + +!> Clear pointers and dealocate memory +subroutine bkgnd_mixing_end(CS) + type(bkgnd_mixing_cs), pointer :: CS ! Control structure + + deallocate(CS%kd_bkgnd) + deallocate(CS%kv_bkgnd) + deallocate(CS) + +end subroutine bkgnd_mixing_end + + +end module MOM_bkgnd_mixing From 08db85c6c925cf4605b0a6aeb03701f553bf0bae Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 21 Mar 2018 18:32:50 -0600 Subject: [PATCH 21/40] Deleted code that has been moved to MOM_bkgnd_mixing.F90 --- .../vertical/MOM_set_diffusivity.F90 | 414 +----------------- 1 file changed, 24 insertions(+), 390 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 345f602902..fe07edea8a 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -21,6 +21,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_bkgnd_mixing, only : calculate_bkgnd_mixing, bkgnd_mixing_init, bkgnd_mixing_cs +use MOM_bkgnd_mixing, only : bkgnd_mixing_end use MOM_string_functions, only : uppercase use MOM_thickness_diffuse, only : vert_fill_TS use MOM_variables, only : thermo_var_ptrs, vertvisc_type, p3d @@ -48,44 +50,6 @@ module MOM_set_diffusivity ! large enough that N2 > omega2. The full expression for ! the Flux Richardson number is usually ! FLUX_RI_MAX*N2/(N2+OMEGA2). The default is 0.2. - logical :: Henyey_IGW_background ! If true, use a simplified variant of the - ! Henyey et al, JGR (1986) latitudinal scaling for - ! the background diapycnal diffusivity, which gives - ! a marked decrease in the diffusivity near the - ! equator. The simplification here is to assume - ! that the in-situ stratification is the same as - ! the reference stratificaiton. - logical :: Henyey_IGW_background_new ! same as Henyey_IGW_background - ! but incorporate the effect of - ! stratification on TKE dissipation, - ! - ! e = f/f_0 * acosh(N/f) / acosh(N_0/f_0) * e_0 - ! - ! where e is the TKE dissipation, and N_0 and f_0 are the - ! reference buoyancy frequency and inertial frequencies respectively. - ! e_0 is the reference dissipation at (N_0,f_0). In the - ! previous version, N=N_0. - ! - ! Additionally, the squared inverse relationship between - ! diapycnal diffusivities and stratification is included - ! - ! kd = e/N^2 - ! - ! where kd is the diapycnal diffusivity. - ! This approach assumes that work done - ! against gravity is uniformly distributed - ! throughout the column. Whereas, kd=kd_0*e, - ! as in the original version, concentrates buoyancy - ! work in regions of strong stratification. - - logical :: Kd_tanh_lat_fn ! If true, use the tanh dependence of Kd_sfc on - ! latitude, like GFDL CM2.1/CM2M. There is no physical - ! justification for this form, and it can not be - ! used with Henyey_IGW_background. - real :: Kd_tanh_lat_scale ! A nondimensional scaling for the range of - ! diffusivities with Kd_tanh_lat_fn. Valid values - ! are in the range of -2 to 2; 0.4 reproduces CM2M. - logical :: bottomdraglaw ! If true, the bottom stress is calculated with a ! drag law c_drag*|u|*u. logical :: BBL_mixing_as_max ! If true, take the maximum of the diffusivity @@ -111,21 +75,6 @@ module MOM_set_diffusivity ! when bulkmixedlayer==.false. real :: Hmix ! mixed layer thickness (meter) when ! bulkmixedlayer==.false. - - logical :: Bryan_Lewis_diffusivity ! If true, background vertical diffusivity - ! uses Bryan-Lewis (1979) like tanh profile. - real :: Kd_Bryan_Lewis_deep ! abyssal value of Bryan-Lewis profile (m2/s) - real :: Kd_Bryan_Lewis_surface ! surface value of Bryan-Lewis profile (m2/s) - real :: Bryan_Lewis_depth_cent ! center of transition depth in Bryan-Lewis (meter) - real :: Bryan_Lewis_width_trans ! width of transition for Bryan-Lewis (meter) - - real :: N0_2Omega ! ratio of the typical Buoyancy frequency to - ! twice the Earth's rotation period, used with the - ! Henyey scaling from the mixing - real :: N2_FLOOR_IOMEGA2 ! floor applied to N2(k) scaled by Omega^2 - ! If =0., N2(k) is positive definite - ! If =1., N2(k) > Omega^2 everywhere - type(diag_ctrl), pointer :: diag ! structure to regulate diagn output timing real :: Int_tide_decay_scale ! decay scale for internal wave TKE (meter) @@ -243,6 +192,7 @@ module MOM_set_diffusivity 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(bkgnd_mixing_cs), pointer :: bkgnd_mixing_csp => NULL() type(int_tide_CS), pointer :: int_tide_CSp => NULL() integer :: id_TKE_itidal = -1 @@ -358,10 +308,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! local variables real, dimension(SZI_(G)) :: & - depth, & ! distance from surface of an interface (meter) N2_bot ! bottom squared buoyancy frequency (1/s2) - real, dimension(SZI_(G), SZJ_(G)) :: & - Kd_sfc ! surface value of the diffusivity (m2/s) type(diffusivity_diags) :: dd ! structure w/ arrays of pointers to avail diags @@ -382,22 +329,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & KT_extra, & ! double difusion diffusivity on temperature (m2/sec) KS_extra ! double difusion diffusivity on salinity (m2/sec) - real :: I_trans ! inverse of the transitional for Bryan-Lewis (1/m) - real :: depth_c ! depth of the center of a layer (meter) - real :: I_Hmix ! inverse of fixed mixed layer thickness (1/m) real :: I_Rho0 ! inverse of Boussinesq density (m3/kg) - real :: I_x30 ! 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg))) - real :: abs_sin ! absolute value of sine of latitude (nondim) - real :: atan_fn_sfc ! surface value of Bryan-Lewis profile (nondim) - real :: atan_fn_lay ! value of Bryan-Lewis profile in layer middle (nondim) - real :: I_atan_fn ! inverse of change in Bryan-Lewis profile from surface to infinite depth (nondim) - real :: deg_to_rad ! factor converting degrees to radians, pi/180. real :: dissip ! local variable for dissipation calculations (W/m3) real :: Omega2 ! squared absolute rotation rate (1/s2) - real :: I_2Omega ! 1/(2 Omega) (sec) - real :: N_2Omega - real :: N02_N2 - real :: epsilon logical :: use_EOS ! If true, compute density from T/S using equation of state. type(p3d) :: z_ptrs(6) ! pointers to diagns to be interpolated into depth space @@ -424,10 +358,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & I_Rho0 = 1.0/GV%Rho0 kappa_fill = 1.e-3 ! m2 s-1 dt_fill = 7200. - deg_to_rad = atan(1.0)/45.0 ! = PI/180 Omega2 = CS%Omega*CS%Omega - I_2Omega = 0.5/CS%Omega - epsilon = 1.e-10 use_EOS = associated(tv%eqn_of_state) @@ -548,42 +479,15 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & visc%Kv_turb(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled endif -! Calculate the diffusivity, Kd, for each layer. This would be -! the appropriate place to add a depth-dependent parameterization or -! another explicit parameterization of Kd. + ! Calculate the diffusivity, Kd, for each layer. This would be + ! the appropriate place to add a depth-dependent parameterization or + ! another explicit parameterization of Kd. - if (CS%Bryan_Lewis_diffusivity) then -!$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc) - do j=js,je ; do i=is,ie - Kd_sfc(i,j) = CS%Kd_Bryan_Lewis_surface - enddo ; enddo - else -!$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc) - do j=js,je ; do i=is,ie - Kd_sfc(i,j) = CS%Kd - enddo ; enddo - endif - if (CS%Henyey_IGW_background) then - I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. -!$OMP parallel do default(none) shared(is,ie,js,je,Kd_sfc,CS,G,deg_to_rad,epsilon,I_x30) & -!$OMP private(abs_sin) - do j=js,je ; do i=is,ie - abs_sin = abs(sin(G%geoLatT(i,j)*deg_to_rad)) - Kd_sfc(i,j) = max(CS%Kd_min, Kd_sfc(i,j) * & - ((abs_sin * invcosh(CS%N0_2Omega/max(epsilon,abs_sin))) * I_x30) ) - enddo ; enddo - elseif (CS%Kd_tanh_lat_fn) then -!$OMP parallel do default(none) shared(is,ie,js,je,Kd_sfc,CS,G) - do j=js,je ; do i=is,ie - ! The transition latitude and latitude range are hard-scaled here, since - ! this is not really intended for wide-spread use, but rather for - ! comparison with CM2M / CM2.1 settings. - Kd_sfc(i,j) = max(CS%Kd_min, Kd_sfc(i,j) * (1.0 + & - CS%Kd_tanh_lat_scale * 0.5*tanh((abs(G%geoLatT(i,j)) - 35.0)/5.0) )) - enddo ; enddo - endif + ! GMM, call MOM_bkgnd_mixing_calc here + call calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS%bkgnd_mixing_csp) + +! GMM, fix OMP calls below - if (CS%debug) call hchksum(Kd_sfc,"Kd_sfc",G%HI,haloshift=0) !$OMP parallel do default(none) shared(is,ie,js,je,nz,G,GV,CS,h,tv,T_f,S_f,fluxes,dd, & !$OMP Kd,Kd_sfc,epsilon,deg_to_rad,I_2Omega,visc, & !$OMP Kd_int,dt,u,v,Omega2) & @@ -592,59 +496,15 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & !$OMP I_x30,abs_sin,N_2Omega,N02_N2,KT_extra, KS_extra, & !$OMP TKE_to_Kd,maxTKE,dissip,kb) do j=js,je - ! Set up variables related to the stratification. - call find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, N2_lay, N2_int, N2_bot) - if (associated(dd%N2_3d)) then - do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo - endif - - ! Set up the background diffusivity. - if (CS%Bryan_Lewis_diffusivity) then - I_trans = 1.0 / CS%Bryan_Lewis_width_trans - atan_fn_sfc = atan(CS%Bryan_Lewis_depth_cent*I_trans) - I_atan_fn = 1.0 / (2.0*atan(1.0) + atan_fn_sfc) - do i=is,ie ; depth(i) = 0.0 ; enddo - do k=1,nz ; do i=is,ie - atan_fn_lay = atan((CS%Bryan_Lewis_depth_cent - & - (depth(i)+0.5*GV%H_to_m*h(i,j,k)))*I_trans) - Kd(i,j,k) = Kd_sfc(i,j) + (CS%Kd_Bryan_Lewis_deep - Kd_sfc(i,j)) * & - (atan_fn_sfc - atan_fn_lay) * I_atan_fn - depth(i) = depth(i) + GV%H_to_m*h(i,j,k) - enddo ; enddo - elseif ((.not.CS%bulkmixedlayer) .and. (CS%Kd /= CS%Kdml)) then - I_Hmix = 1.0 / CS%Hmix - do i=is,ie ; depth(i) = 0.0 ; enddo - do k=1,nz ; do i=is,ie - depth_c = depth(i) + 0.5*GV%H_to_m*h(i,j,k) - - if (depth_c <= CS%Hmix) then ; Kd(i,j,k) = CS%Kdml - elseif (depth_c >= 2.0*CS%Hmix) then ; Kd(i,j,k) = Kd_sfc(i,j) - else - Kd(i,j,k) = ((Kd_sfc(i,j) - CS%Kdml) * I_Hmix) * depth_c + & - (2.0*CS%Kdml - Kd_sfc(i,j)) - endif - - depth(i) = depth(i) + GV%H_to_m*h(i,j,k) - enddo ; enddo - elseif (CS%Henyey_IGW_background_new) then - I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. - do k=1,nz ; do i=is,ie - abs_sin = max(epsilon,abs(sin(G%geoLatT(i,j)*deg_to_rad))) - N_2Omega = max(abs_sin,sqrt(N2_lay(i,k))*I_2Omega) - N02_N2 = (CS%N0_2Omega/N_2Omega)**2 - Kd(i,j,k) = max(CS%Kd_min, Kd_sfc(i,j) * & - ((abs_sin * invcosh(N_2Omega/abs_sin)) * I_x30)*N02_N2) - enddo ; enddo - ! GMM, CVMix "internal" bg mixing can go here - !elseif (CS%use_cvmix_internal??) then - - - else - do k=1,nz ; do i=is,ie - Kd(i,j,k) = Kd_sfc(i,j) - enddo ; enddo - endif + ! Add vertical background diffusivities and viscosities + do k=1,nz ; do i=is,ie + ! diffusivities + Kd(i,j,k) = Kd(i,j,k) + 0.5*(CS%bkgnd_mixing_csp%kd_bkgnd(i,j,K) + CS%bkgnd_mixing_csp%kd_bkgnd(i,j,K+1)) + ! viscosities + ! GMM, will this be done here? + enddo ; enddo + ! 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 @@ -1105,159 +965,7 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, CS, & end subroutine find_TKE_to_Kd -subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & - N2_lay, N2_int, N2_bot) - 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(in) :: tv - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T_f, S_f - type(forcing), intent(in) :: fluxes - integer, intent(in) :: j - type(set_diffusivity_CS), pointer :: CS - real, dimension(SZI_(G),SZK_(G)+1), intent(out) :: dRho_int, N2_int - real, dimension(SZI_(G),SZK_(G)), intent(out) :: N2_lay - real, dimension(SZI_(G)), intent(out) :: N2_bot - - real, dimension(SZI_(G),SZK_(G)+1) :: & - dRho_int_unfilt, & ! unfiltered density differences across interfaces - dRho_dT, & ! partial derivative of density wrt temp (kg m-3 degC-1) - dRho_dS ! partial derivative of density wrt saln (kg m-3 PPT-1) - - real, dimension(SZI_(G)) :: & - pres, & ! pressure at each interface (Pa) - Temp_int, & ! temperature at each interface (degC) - Salin_int, & ! salinity at each interface (PPT) - drho_bot, & - h_amp, & - hb, & - z_from_bot - - real :: Rml_base ! density of the deepest variable density layer - real :: dz_int ! thickness associated with an interface (meter) - real :: G_Rho0 ! gravitation acceleration divided by Bouss reference density (m4 s-2 kg-1) - real :: H_neglect ! negligibly small thickness, in the same units as h. - - logical :: do_i(SZI_(G)), do_any - integer :: i, k, is, ie, nz - - is = G%isc ; ie = G%iec ; nz = G%ke - G_Rho0 = GV%g_Earth / GV%Rho0 - H_neglect = GV%H_subroundoff - - ! Find the (limited) density jump across each interface. - do i=is,ie - dRho_int(i,1) = 0.0 ; dRho_int(i,nz+1) = 0.0 - dRho_int_unfilt(i,1) = 0.0 ; dRho_int_unfilt(i,nz+1) = 0.0 - enddo - if (associated(tv%eqn_of_state)) then - if (associated(fluxes%p_surf)) then - do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo - else - do i=is,ie ; pres(i) = 0.0 ; enddo - endif - 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) + T_f(i,j,k-1)) - Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1)) - enddo - call calculate_density_derivs(Temp_int, Salin_int, pres, & - dRho_dT(:,K), dRho_dS(:,K), is, ie-is+1, tv%eqn_of_state) - do i=is,ie - dRho_int(i,K) = max(dRho_dT(i,K)*(T_f(i,j,k) - T_f(i,j,k-1)) + & - dRho_dS(i,K)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0) - dRho_int_unfilt(i,K) = max(dRho_dT(i,K)*(tv%T(i,j,k) - tv%T(i,j,k-1)) + & - dRho_dS(i,K)*(tv%S(i,j,k) - tv%S(i,j,k-1)), 0.0) - enddo - enddo - else - do K=2,nz ; do i=is,ie - dRho_int(i,K) = GV%Rlay(k) - GV%Rlay(k-1) - enddo ; enddo - endif - - ! Set the buoyancy frequencies. - do k=1,nz ; do i=is,ie - N2_lay(i,k) = G_Rho0 * 0.5*(dRho_int(i,K) + dRho_int(i,K+1)) / & - (GV%H_to_m*(h(i,j,k) + H_neglect)) - enddo ; enddo - do i=is,ie ; N2_int(i,1) = 0.0 ; N2_int(i,nz+1) = 0.0 ; enddo - do K=2,nz ; do i=is,ie - N2_int(i,K) = G_Rho0 * dRho_int(i,K) / & - (0.5*GV%H_to_m*(h(i,j,k-1) + h(i,j,k) + H_neglect)) - enddo ; enddo - - ! Find the bottom boundary layer stratification, and use this in the deepest layers. - do i=is,ie - hb(i) = 0.0 ; dRho_bot(i) = 0.0 - z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) - do_i(i) = (G%mask2dT(i,j) > 0.5) - - if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation) then - h_amp(i) = sqrt(CS%h2(i,j)) ! for computing Nb - else - h_amp(i) = 0.0 - endif - enddo - - do k=nz,2,-1 - do_any = .false. - do i=is,ie ; if (do_i(i)) then - dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) - z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above - - hb(i) = hb(i) + dz_int - dRho_bot(i) = dRho_bot(i) + dRho_int(i,K) - - if (z_from_bot(i) > h_amp(i)) then - if (k>2) then - ! Always include at least one full layer. - hb(i) = hb(i) + 0.5*GV%H_to_m*(h(i,j,k-1) + h(i,j,k-2)) - dRho_bot(i) = dRho_bot(i) + dRho_int(i,K-1) - endif - do_i(i) = .false. - else - do_any = .true. - endif - endif ; enddo - if (.not.do_any) exit - enddo - - do i=is,ie - if (hb(i) > 0.0) then - N2_bot(i) = (G_Rho0 * dRho_bot(i)) / hb(i) - else ; N2_bot(i) = 0.0 ; endif - z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) - do_i(i) = (G%mask2dT(i,j) > 0.5) - enddo - - do k=nz,2,-1 - do_any = .false. - do i=is,ie ; if (do_i(i)) then - dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) - z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above - - N2_int(i,K) = N2_bot(i) - if (k>2) N2_lay(i,k-1) = N2_bot(i) - - if (z_from_bot(i) > h_amp(i)) then - if (k>2) N2_int(i,K-1) = N2_bot(i) - do_i(i) = .false. - else - do_any = .true. - endif - endif ; enddo - if (.not.do_any) exit - enddo - - if (associated(tv%eqn_of_state)) then - do K=1,nz+1 ; do i=is,ie - dRho_int(i,K) = dRho_int_unfilt(i,K) - enddo ; enddo - endif - -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 @@ -2248,6 +1956,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, endif ! Polzin end subroutine add_int_tide_diffusivity + !> This subroutine calculates several properties related to bottom !! boundary layer turbulence. subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, CS) @@ -2666,79 +2375,8 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp "for an isopycnal layer-formulation.", & default=.false.) - call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", & - CS%Bryan_Lewis_diffusivity, & - "If true, use a Bryan & Lewis (JGR 1979) like tanh \n"//& - "profile of background diapycnal diffusivity with depth.", & - default=.false.) - if (CS%Bryan_Lewis_diffusivity) then - call get_param(param_file, mdl, "KD_BRYAN_LEWIS_DEEP", & - CS%Kd_Bryan_Lewis_deep, & - "The abyssal value of a Bryan-Lewis diffusivity profile. \n"//& - "KD_BRYAN_LEWIS_DEEP is only used if \n"//& - "BRYAN_LEWIS_DIFFUSIVITY is true.", units="m2 s-1", & - fail_if_missing=.true.) - call get_param(param_file, mdl, "KD_BRYAN_LEWIS_SURFACE", & - CS%Kd_Bryan_Lewis_surface, & - "The surface value of a Bryan-Lewis diffusivity profile. \n"//& - "KD_BRYAN_LEWIS_SURFACE is only used if \n"//& - "BRYAN_LEWIS_DIFFUSIVITY is true.", units="m2 s-1", & - fail_if_missing=.true.) - call get_param(param_file, mdl, "BRYAN_LEWIS_DEPTH_CENT", & - CS%Bryan_Lewis_depth_cent, & - "The depth about which the transition in the Bryan-Lewis \n"//& - "profile is centered. BRYAN_LEWIS_DEPTH_CENT is only \n"//& - "used if BRYAN_LEWIS_DIFFUSIVITY is true.", units="m", & - fail_if_missing=.true.) - call get_param(param_file, mdl, "BRYAN_LEWIS_WIDTH_TRANS", & - CS%Bryan_Lewis_width_trans, & - "The width of the transition in the Bryan-Lewis \n"//& - "profile. BRYAN_LEWIS_WIDTH_TRANS is only \n"//& - "used if BRYAN_LEWIS_DIFFUSIVITY is true.", units="m", & - fail_if_missing=.true.) - endif - - call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND", & - CS%Henyey_IGW_background, & - "If true, use a latitude-dependent scaling for the near \n"//& - "surface background diffusivity, as described in \n"//& - "Harrison & Hallberg, JPO 2008.", default=.false.) - call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND_NEW", & - CS%Henyey_IGW_background_new, & - "If true, use a better latitude-dependent scaling for the\n"//& - "background diffusivity, as described in \n"//& - "Harrison & Hallberg, JPO 2008.", default=.false.) - if (CS%Henyey_IGW_background .and. CS%Henyey_IGW_background_new) call MOM_error(FATAL, & - "set_diffusivity_init: HENYEY_IGW_BACKGROUND and HENYEY_IGW_BACKGROUND_NEW "// & - "are mutually exclusive. Set only one or none.") - if (CS%Henyey_IGW_background) & - call get_param(param_file, mdl, "HENYEY_N0_2OMEGA", CS%N0_2Omega, & - "The ratio of the typical Buoyancy frequency to twice \n"//& - "the Earth's rotation period, used with the Henyey \n"//& - "scaling from the mixing.", units="nondim", default=20.0) - call get_param(param_file, mdl, "N2_FLOOR_IOMEGA2", CS%N2_FLOOR_IOMEGA2, & - "The floor applied to N2(k) scaled by Omega^2:\n"//& - "\tIf =0., N2(k) is simply positive definite.\n"//& - "\tIf =1., N2(k) > Omega^2 everywhere.", units="nondim", & - default=1.0) - - call get_param(param_file, mdl, "KD_TANH_LAT_FN", & - CS%Kd_tanh_lat_fn, & - "If true, use a tanh dependence of Kd_sfc on latitude, \n"//& - "like CM2.1/CM2M. There is no physical justification \n"//& - "for this form, and it can not be used with \n"//& - "HENYEY_IGW_BACKGROUND.", default=.false.) - if (CS%Kd_tanh_lat_fn) & - call get_param(param_file, mdl, "KD_TANH_LAT_SCALE", & - CS%Kd_tanh_lat_scale, & - "A nondimensional scaling for the range ofdiffusivities \n"//& - "with KD_TANH_LAT_FN. Valid values are in the range of \n"//& - "-2 to 2; 0.4 reproduces CM2M.", units="nondim", default=0.0) - - call get_param(param_file, mdl, "KV", CS%Kv, & - "The background kinematic viscosity in the interior. \n"//& - "The molecular value, ~1e-6 m2 s-1, may be used.", & - units="m2 s-1", fail_if_missing=.true.) + ! set params releted to the background mixing + call bkgnd_mixing_init(Time, G, GV, param_file, CS%diag, CS%bkgnd_mixing_csp) call get_param(param_file, mdl, "KD", CS%Kd, & "The background diapycnal diffusivity of density in the \n"//& @@ -3130,12 +2768,6 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp endif endif - if (CS%Int_tide_dissipation .and. CS%Bryan_Lewis_diffusivity) & - call MOM_error(FATAL,"MOM_Set_Diffusivity: "// & - "Bryan-Lewis and internal tidal dissipation are both enabled. Choose one.") - if (CS%Henyey_IGW_background .and. CS%Kd_tanh_lat_fn) call MOM_error(FATAL, & - "Set_diffusivity: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.") - if (CS%user_change_diff) then call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) endif @@ -3153,6 +2785,8 @@ end subroutine set_diffusivity_init subroutine set_diffusivity_end(CS) type(set_diffusivity_CS), pointer :: CS !< Control structure for this module + call bkgnd_mixing_end(CS%bkgnd_mixing_csp) + if (CS%user_change_diff) & call user_change_diff_end(CS%user_change_diff_CSp) From 8ef5a1b93d03957b895221ac4db08bbb0025e7d9 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 23 Mar 2018 15:06:41 -0600 Subject: [PATCH 22/40] Renamed Bryan&Lewis coeffs and added a 3D array for depth --- .../vertical/MOM_bkgnd_mixing.F90 | 104 +++++++++--------- 1 file changed, 53 insertions(+), 51 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 3ad9dfd638..e561b75a3e 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -11,6 +11,7 @@ module MOM_bkgnd_mixing use MOM_variables, only : thermo_var_ptrs use MOM_forcing_type, only : forcing use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_error_handler, only : is_root_pe use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_debugging, only : hchksum use MOM_grid, only : ocean_grid_type @@ -29,13 +30,14 @@ module MOM_bkgnd_mixing type, public :: bkgnd_mixing_cs ! Parameters - real :: Kd_Bryan_Lewis_deep !< The abyssal value of a Bryan-Lewis diffusivity profile - !! (m2/s) - real :: Kd_Bryan_Lewis_surface !< "The surface value of a Bryan-Lewis diffusivity profile - !! (m2/s) - real :: Bryan_Lewis_depth_cent !< The depth about which the transition in the Bryan-Lewis - !! is centered (m) - real :: Bryan_Lewis_width_trans!< The width of the transition in the Bryan-Lewis profile (m) + real :: Bryan_Lewis_c1 !< The vertical diffusivity values for Bryan-Lewis profile + !! at |z|=D (m2/s) + real :: Bryan_Lewis_c2 !< The amplitude of variation in diffusivity for the + !! Bryan-Lewis diffusivity profile (m2/s) + real :: Bryan_Lewis_c3 !< The inverse length scale for transition region in the + !! Bryan-Lewis diffusivity profile (1/m) + real :: Bryan_Lewis_c4 !< The depth where diffusivity is Bryan_Lewis_bl1 in the + !! Bryan-Lewis profile (m) real :: Kd_min !< minimum diapycnal diffusivity (m2/s) real :: Kd !< interior diapycnal diffusivity (m2/s) real :: N0_2Omega !< ratio of the typical Buoyancy frequency to @@ -153,7 +155,7 @@ subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mdl, "PRANDTL_TURB", CS%prandtl_turb, & units="nondim", default=1.0, do_not_log=.true.) - call openParameterBlock(param_file,'MOM_BACKGROUND_MIXING') +! call openParameterBlock(param_file,'MOM_BACKGROUND_MIXING') call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", & CS%Bryan_Lewis_diffusivity, & @@ -163,24 +165,24 @@ subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) if (CS%Bryan_Lewis_diffusivity) then - call get_param(param_file, mdl, "KD_BRYAN_LEWIS_DEEP", & - CS%Kd_Bryan_Lewis_deep, & - "The abyssal value of a Bryan-Lewis diffusivity profile.", & + call get_param(param_file, mdl, "BRYAN_LEWIS_C1", & + CS%Bryan_Lewis_c1, & + "The vertical diffusivity values for Bryan-Lewis profile at |z|=D.", & units="m2 s-1", fail_if_missing=.true.) - call get_param(param_file, mdl, "KD_BRYAN_LEWIS_SURFACE", & - CS%Kd_Bryan_Lewis_surface, & - "The surface value of a Bryan-Lewis diffusivity profile.", & + call get_param(param_file, mdl, "BRYAN_LEWIS_C2", & + CS%Bryan_Lewis_c2, & + "The amplitude of variation in diffusivity for the Bryan-Lewis profile", & units="m2 s-1", fail_if_missing=.true.) - call get_param(param_file, mdl, "BRYAN_LEWIS_DEPTH_CENT", & - CS%Bryan_Lewis_depth_cent, & - "The depth about which the transition in the Bryan-Lewis.", & - units="m", fail_if_missing=.true.) + call get_param(param_file, mdl, "BRYAN_LEWIS_C3", & + CS%Bryan_Lewis_c3, & + "The inverse length scale for transition region in the Bryan-Lewis profile", & + units="m-1", fail_if_missing=.true.) - call get_param(param_file, mdl, "BRYAN_LEWIS_WIDTH_TRANS", & - CS%Bryan_Lewis_width_trans, & - "The width of the transition in the Bryan-Lewis.",& + call get_param(param_file, mdl, "BRYAN_LEWIS_C4", & + CS%Bryan_Lewis_c4, & + "The depth where diffusivity is BRYAN_LEWIS_C1 in the Bryan-Lewis profile",& units="m", fail_if_missing=.true.) endif ! CS%Bryan_Lewis_diffusivity @@ -224,7 +226,7 @@ subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) if (CS%Henyey_IGW_background .and. CS%Kd_tanh_lat_fn) call MOM_error(FATAL, & "MOM_bkgnd_mixing: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.") - call closeParameterBlock(param_file) +! call closeParameterBlock(param_file) ! allocate arrays and set them to zero allocate(CS%kd_bkgnd(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kd_bkgnd(:,:,:) = 0. @@ -232,9 +234,9 @@ subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) ! Register diagnostics CS%diag => diag - CS%id_kd_bkgnd = register_diag_field('ocean_model', 'bkgnd_kd', diag%axesTi, Time, & + CS%id_kd_bkgnd = register_diag_field('ocean_model', 'Kd_bkgnd', diag%axesTi, Time, & 'Background diffusivity added by MOM_bkgnd_mixing module', 'm2/s') - CS%id_kv_bkgnd = register_diag_field('ocean_model', 'bkgnd_kv', diag%axesTi, Time, & + CS%id_kv_bkgnd = register_diag_field('ocean_model', 'Kv_bkgnd', diag%axesTi, Time, & 'Background viscosity added by MOM_bkgnd_mixing module', 'm2/s') end subroutine bkgnd_mixing_init @@ -256,6 +258,7 @@ subroutine calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS) !! a previous call to bkgnd_mixing_init. ! local variables + real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: depth_3d !< distance from surface of an interface (m) real, dimension(SZI_(G), SZJ_(G)) :: Kd_sfc !< surface value of the diffusivity (m2/s) real, dimension(SZI_(G)) :: & depth_i, & !< distance from surface of an interface (meter) @@ -284,32 +287,7 @@ subroutine calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS) deg_to_rad = atan(1.0)/45.0 ! = PI/180 epsilon = 1.e-10 - if (CS%Bryan_Lewis_diffusivity) then -!$OMP parallel do default(none) shared(is,ie,js,je,nz, CS) -!$OMP private(cvmix_init_bkgnd,cvmix_coeffs_bkgnd) - - ! Bryan & Lewis is computed via CVMix - do j=js,je; do i=is,ie - - depth_k(:) = 0.0 - do k=1,nz - depth_k(k) = depth_k(k) + GV%H_to_m*h(i,j,k) - enddo - - call cvmix_init_bkgnd(max_nlev=nz, & - zw = depth_k(:), & !< interface depth, must be positive. - bl1 = CS%Kd_Bryan_Lewis_deep, & - bl2 = CS%Kd_Bryan_Lewis_surface, & - bl3 = 1.0/CS%Bryan_Lewis_depth_cent , & - bl4 = CS%Bryan_Lewis_width_trans, & - prandtl = CS%prandtl_turb) - - call cvmix_coeffs_bkgnd(Mdiff_out=CS%kv_bkgnd(i,j,:), & - Tdiff_out=CS%kd_bkgnd(i,j,:), & - nlev=nz, & - max_nlev=nz) - enddo; enddo - else + if (.not. CS%Bryan_Lewis_diffusivity) then !$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc) do j=js,je ; do i=is,ie Kd_sfc(i,j) = CS%Kd @@ -349,6 +327,8 @@ subroutine calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS) !$OMP I_x30,abs_sin,N_2Omega,N02_N2,KT_extra, !KS_extra, & !$OMP TKE_to_Kd,maxTKE,dissip,kb) + + depth_3d(:,:,:) = 0.0 do j=js,je ! Set up variables related to the stratification. call find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, dRho_int, N2_lay, N2_int, N2_bot) @@ -357,7 +337,29 @@ subroutine calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS) !endif ! Set up the background diffusivity. - if ((.not. CS%Bryan_Lewis_diffusivity) .and. (.not.CS%bulkmixedlayer) .and. & + if (CS%Bryan_Lewis_diffusivity) then + + do i=is,ie + !depth_k(:) = 0.0 + do k=2,nz+1 + depth_3d(i,j,k) = depth_3d(i,j,k-1) + GV%H_to_m*h(i,j,k-1) + enddo + ! if (is_root_pe()) write(*,*)'depth_3d(i,j,:)',depth_3d(i,j,:) + + call cvmix_init_bkgnd(max_nlev=nz, & + zw = depth_3d(i,j,:), & !< interface depth, must bepositive. + bl1 = CS%Bryan_Lewis_c1, & + bl2 = CS%Bryan_Lewis_c2, & + bl3 = CS%Bryan_Lewis_c3, & + bl4 = CS%Bryan_Lewis_c4, & + prandtl = CS%prandtl_turb) + + call cvmix_coeffs_bkgnd(Mdiff_out=CS%kv_bkgnd(i,j,:), & + Tdiff_out=CS%kd_bkgnd(i,j,:), & + nlev=nz, & + max_nlev=nz) + enddo + elseif ((.not. CS%Bryan_Lewis_diffusivity) .and. (.not.CS%bulkmixedlayer) .and. & (CS%Kd/= CS%Kdml)) then I_Hmix = 1.0 / CS%Hmix From 4e8dd5a9bacbc47988d2fd99b50defc00580948e Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 23 Mar 2018 15:07:41 -0600 Subject: [PATCH 23/40] Added a if statement to check if Int_tide_dissipation and Bryan_Lewis are used at the same time. --- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index fe07edea8a..c127ad9a72 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -2772,6 +2772,10 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) endif + if (CS%Int_tide_dissipation .and. CS%bkgnd_mixing_csp%Bryan_Lewis_diffusivity) & + call MOM_error(FATAL,"MOM_Set_Diffusivity: "// & + "Bryan-Lewis and internal tidal dissipation are both enabled. Choose one.") + CS%useKappaShear = kappa_shear_init(Time, G, GV, param_file, CS%diag, CS%kappaShear_CSp) if (CS%useKappaShear) & id_clock_kappaShear = cpu_clock_id('(Ocean kappa_shear)', grain=CLOCK_MODULE) From 7a81871673fa0458bf367dbbf991de84e0d52aaf Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 26 Mar 2018 13:47:12 -0600 Subject: [PATCH 24/40] Re-structured MOM_bkgnd_mixing * Moved find_N2 back to MOM_set_diffusivity * Added a new func. in MOM_bkgnd_mixing (sfc_bkgnd_mixing) * Modified calculate_bkgnd_mixing These changes *do not* change answers for ocean_only/global_ALE/z --- .../vertical/MOM_bkgnd_mixing.F90 | 404 ++++++------------ .../vertical/MOM_set_diffusivity.F90 | 196 ++++++++- 2 files changed, 304 insertions(+), 296 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index e561b75a3e..ba7711586c 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -24,7 +24,10 @@ module MOM_bkgnd_mixing #include -public bkgnd_mixing_init, bkgnd_mixing_end, calculate_bkgnd_mixing +public bkgnd_mixing_init +public bkgnd_mixing_end +public calculate_bkgnd_mixing +public sfc_bkgnd_mixing !> Control structure including parameters for this module. type, public :: bkgnd_mixing_cs @@ -81,6 +84,7 @@ module MOM_bkgnd_mixing type(diag_ctrl), pointer :: diag => NULL() integer :: id_kd_bkgnd = -1, id_kv_bkgnd = -1 + real, allocatable, dimension(:,:) :: Kd_sfc !< surface value of the diffusivity (m2/s) ! Diagnostics arrays real, allocatable, dimension(:,:,:) :: kd_bkgnd !< Background diffusivity (m2/s) real, allocatable, dimension(:,:,:) :: kv_bkgnd !< Background viscosity (m2/s) @@ -149,7 +153,6 @@ subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) "mixed layer is not used.", units="m", fail_if_missing=.true.) endif - call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) call get_param(param_file, mdl, "PRANDTL_TURB", CS%prandtl_turb, & @@ -231,6 +234,7 @@ subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) ! allocate arrays and set them to zero allocate(CS%kd_bkgnd(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kd_bkgnd(:,:,:) = 0. allocate(CS%kv_bkgnd(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kv_bkgnd(:,:,:) = 0. + allocate(CS%Kd_sfc(SZI_(G), SZJ_(G))); CS%Kd_sfc(:,:) = 0. ! Register diagnostics CS%diag => diag @@ -241,56 +245,30 @@ subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) end subroutine bkgnd_mixing_init -!> Subroutine for calculating vertical background diffusivities/viscosities -subroutine calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS) +!> Get surface vertical background diffusivities/viscosities. +subroutine sfc_bkgnd_mixing(G, 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)), intent(in) :: T_f!< temperature (deg C), after massless - !! layers filled vertically by diffusion. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: S_f!< salinity, after massless - !! layers filled vertically by diffusion. - type(forcing), intent(in) :: fluxes !< Structure of surface fluxes that may be - !! used. - type(bkgnd_mixing_cs), pointer :: CS !< The control structure returned by + type(bkgnd_mixing_cs), pointer, intent(inout) :: CS !< The control structure returned by !! a previous call to bkgnd_mixing_init. - ! local variables - real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: depth_3d !< distance from surface of an interface (m) - real, dimension(SZI_(G), SZJ_(G)) :: Kd_sfc !< surface value of the diffusivity (m2/s) - real, dimension(SZI_(G)) :: & - depth_i, & !< distance from surface of an interface (meter) - N2_bot !< bottom squared buoyancy frequency (1/s2) - 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? - real, dimension(SZI_(G),SZK_(G)) :: & - N2_lay !< squared buoyancy frequency associated with layers (1/s2) - - real, dimension(SZK_(G)) :: depth_k !< distance from surface of an interface (meter) real :: I_x30 !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg))) real :: deg_to_rad !< factor converting degrees to radians, pi/180. real :: abs_sin !< absolute value of sine of latitude (nondim) - real :: depth_c !< depth of the center of a layer (meter) - real :: I_Hmix !< inverse of fixed mixed layer thickness (1/m) real :: epsilon - real :: I_2Omega !< 1/(2 Omega) (sec) - real :: N_2Omega - real :: N02_N2 - integer :: i, j, k, is, ie, js, je, nz + integer :: i, j, k, is, ie, js, je - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ! set some parameters deg_to_rad = atan(1.0)/45.0 ! = PI/180 epsilon = 1.e-10 + if (.not. CS%Bryan_Lewis_diffusivity) then !$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc) do j=js,je ; do i=is,ie - Kd_sfc(i,j) = CS%Kd + CS%Kd_sfc(i,j) = CS%Kd enddo ; enddo endif @@ -301,7 +279,7 @@ subroutine calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS) !$OMP private(abs_sin) do j=js,je ; do i=is,ie abs_sin = abs(sin(G%geoLatT(i,j)*deg_to_rad)) - Kd_sfc(i,j) = max(CS%Kd_min, Kd_sfc(i,j) * & + CS%Kd_sfc(i,j) = max(CS%Kd_min, CS%Kd_sfc(i,j) * & ((abs_sin * invcosh(CS%N0_2Omega/max(epsilon,abs_sin))) * I_x30) ) enddo ; enddo elseif (CS%Kd_tanh_lat_fn) then @@ -310,276 +288,132 @@ subroutine calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS) ! The transition latitude and latitude range are hard-scaled here, since ! this is not really intended for wide-spread use, but rather for ! comparison with CM2M / CM2.1 settings. - Kd_sfc(i,j) = max(CS%Kd_min, Kd_sfc(i,j) * (1.0 + & + CS%Kd_sfc(i,j) = max(CS%Kd_min, CS%Kd_sfc(i,j) * (1.0 + & CS%Kd_tanh_lat_scale * 0.5*tanh((abs(G%geoLatT(i,j)) - 35.0)/5.0) )) enddo ; enddo endif -!$OMP parallel do default(none) -!shared(is,ie,js,je,nz,G,GV,CS,h,tv,T_f,S_f,fluxes,dd, & -!$OMP -!Kd,Kd_sfc,epsilon,deg_to_rad,I_2Omega,visc, & -!$OMP Kd_int,dt,u,v,Omega2) & -!$OMP -!private(dRho_int,I_trans,atan_fn_sfc,I_atan_fn,atan_fn_lay, & -!$OMP I_Hmix,depth_c,depth,N2_lay, N2_int, -!N2_bot, & -!$OMP I_x30,abs_sin,N_2Omega,N02_N2,KT_extra, -!KS_extra, & -!$OMP TKE_to_Kd,maxTKE,dissip,kb) - - depth_3d(:,:,:) = 0.0 - do j=js,je - ! Set up variables related to the stratification. - call find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, dRho_int, N2_lay, N2_int, N2_bot) - !if (associated(dd%N2_3d)) then - ! do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo - !endif - - ! Set up the background diffusivity. - if (CS%Bryan_Lewis_diffusivity) then - - do i=is,ie - !depth_k(:) = 0.0 - do k=2,nz+1 - depth_3d(i,j,k) = depth_3d(i,j,k-1) + GV%H_to_m*h(i,j,k-1) - enddo - ! if (is_root_pe()) write(*,*)'depth_3d(i,j,:)',depth_3d(i,j,:) - - call cvmix_init_bkgnd(max_nlev=nz, & - zw = depth_3d(i,j,:), & !< interface depth, must bepositive. - bl1 = CS%Bryan_Lewis_c1, & - bl2 = CS%Bryan_Lewis_c2, & - bl3 = CS%Bryan_Lewis_c3, & - bl4 = CS%Bryan_Lewis_c4, & - prandtl = CS%prandtl_turb) - - call cvmix_coeffs_bkgnd(Mdiff_out=CS%kv_bkgnd(i,j,:), & - Tdiff_out=CS%kd_bkgnd(i,j,:), & - nlev=nz, & - max_nlev=nz) - enddo - elseif ((.not. CS%Bryan_Lewis_diffusivity) .and. (.not.CS%bulkmixedlayer) .and. & - (CS%Kd/= CS%Kdml)) then - - I_Hmix = 1.0 / CS%Hmix - do i=is,ie ; depth_i(i) = 0.0 ; enddo - do k=1,nz ; do i=is,ie - depth_c = depth_i(i) + 0.5*GV%H_to_m*h(i,j,k) - - if (depth_c <= CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = CS%Kdml - elseif (depth_c >= 2.0*CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = Kd_sfc(i,j) - else - CS%kd_bkgnd(i,j,k) = ((Kd_sfc(i,j) - CS%Kdml) * I_Hmix) * depth_c + & - (2.0*CS%Kdml - Kd_sfc(i,j)) - endif - - depth_i(i) = depth_i(i) + GV%H_to_m*h(i,j,k) - enddo ; enddo - elseif (CS%Henyey_IGW_background_new) then - I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. - do k=1,nz ; do i=is,ie - abs_sin = max(epsilon,abs(sin(G%geoLatT(i,j)*deg_to_rad))) - N_2Omega = max(abs_sin,sqrt(N2_lay(i,k))*I_2Omega) - N02_N2 = (CS%N0_2Omega/N_2Omega)**2 - CS%kd_bkgnd(i,j,k) = max(CS%Kd_min, Kd_sfc(i,j) * & - ((abs_sin * invcosh(N_2Omega/abs_sin)) * I_x30)*N02_N2) - enddo ; enddo - - else - do k=1,nz ; do i=is,ie - CS%kd_bkgnd(i,j,k) = Kd_sfc(i,j) - enddo ; enddo - endif - enddo ! j-loop - - if (CS%debug) then - call hchksum(Kd_sfc,"Kd_sfc",G%HI,haloshift=0) - call hchksum(CS%kd_bkgnd, "MOM_bkgnd_mixing: kd_bkgnd",G%HI,haloshift=0) - call hchksum(CS%kv_bkgnd, "MOM_bkgnd_mixing: kv_bkgnd",G%HI,haloshift=0) - endif + if (CS%debug) call hchksum(CS%Kd_sfc,"After sfc_bkgnd_mixing: Kd_sfc",G%HI,haloshift=0) - ! send diagnostics to post_data - if (CS%id_kd_bkgnd > 0) call post_data(CS%id_kd_bkgnd, CS%kd_bkgnd, CS%diag) - if (CS%id_kv_bkgnd > 0) call post_data(CS%id_kv_bkgnd, CS%kv_bkgnd, CS%diag) +end subroutine sfc_bkgnd_mixing -end subroutine calculate_bkgnd_mixing +!> Calculates the vertical background diffusivities/viscosities +subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, j, G, GV, CS) -!> Computes N2 -subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, dRho_int, & - N2_lay, N2_int, N2_bot) - 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(in) :: tv - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T_f, S_f - type(forcing), intent(in) :: fluxes + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. + real, dimension(SZI_(G),SZK_(G)), intent(in) :: N2_lay!< squared buoyancy frequency associated + !! with layers (1/s2) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd !< Diapycnal diffusivity of each layer (m2/sec). integer, intent(in) :: j - real, dimension(SZI_(G),SZK_(G)+1), intent(out) :: dRho_int, N2_int - real, dimension(SZI_(G),SZK_(G)), intent(out) :: N2_lay - real, dimension(SZI_(G)), intent(out) :: N2_bot - - real, dimension(SZI_(G),SZK_(G)+1) :: & - dRho_int_unfilt, & ! unfiltered density differences across interfaces - dRho_dT, & ! partial derivative of density wrt temp (kg m-3 degC-1) - dRho_dS ! partial derivative of density wrt saln (kg m-3 PPT-1) + type(bkgnd_mixing_cs), pointer :: CS !< The control structure returned by + !! a previous call to bkgnd_mixing_init. + ! local variables + real, dimension(SZI_(G), SZK_(G)+1) :: depth_2d !< distance from surface of an interface (m) real, dimension(SZI_(G)) :: & - pres, & ! pressure at each interface (Pa) - Temp_int, & ! temperature at each interface (degC) - Salin_int, & ! salinity at each interface (PPT) - drho_bot, & - h_amp, & - hb, & - z_from_bot - - real :: Rml_base ! density of the deepest variable density layer - real :: dz_int ! thickness associated with an interface (meter) - real :: G_Rho0 ! gravitation acceleration divided by Bouss reference density (m4 s-2 kg-1) - real :: H_neglect ! negligibly small thickness, in the same units as h. - - logical :: do_i(SZI_(G)), do_any - integer :: i, k, is, ie, nz - - is = G%isc ; ie = G%iec ; nz = G%ke - G_Rho0 = GV%g_Earth / GV%Rho0 - H_neglect = GV%H_subroundoff - - ! Find the (limited) density jump across each interface. - do i=is,ie - dRho_int(i,1) = 0.0 ; dRho_int(i,nz+1) = 0.0 - dRho_int_unfilt(i,1) = 0.0 ; dRho_int_unfilt(i,nz+1) = 0.0 - enddo - if (associated(tv%eqn_of_state)) then - if (associated(fluxes%p_surf)) then - do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo - else - do i=is,ie ; pres(i) = 0.0 ; enddo - endif - 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) + T_f(i,j,k-1)) - Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1)) + depth !< distance from surface of an interface (meter) + real :: depth_c !< depth of the center of a layer (meter) + real :: I_Hmix !< inverse of fixed mixed layer thickness (1/m) + real :: I_2Omega !< 1/(2 Omega) (sec) + real :: N_2Omega + real :: N02_N2 + real :: I_x30 !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg))) + real :: deg_to_rad !< factor converting degrees to radians, pi/180. + real :: abs_sin !< absolute value of sine of latitude (nondim) + real :: epsilon + integer :: i, k, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + + ! set some parameters + deg_to_rad = atan(1.0)/45.0 ! = PI/180 + epsilon = 1.e-10 + + depth_2d(:,:) = 0.0 + ! Set up the background diffusivity. + if (CS%Bryan_Lewis_diffusivity) then + + do i=is,ie + do k=2,nz+1 + depth_2d(i,k) = depth_2d(i,k-1) + GV%H_to_m*h(i,j,k-1) enddo - call calculate_density_derivs(Temp_int, Salin_int, pres, & - dRho_dT(:,K), dRho_dS(:,K), is, ie-is+1, tv%eqn_of_state) - do i=is,ie - dRho_int(i,K) = max(dRho_dT(i,K)*(T_f(i,j,k) - T_f(i,j,k-1)) + & - dRho_dS(i,K)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0) - dRho_int_unfilt(i,K) = max(dRho_dT(i,K)*(tv%T(i,j,k) - tv%T(i,j,k-1)) + & - dRho_dS(i,K)*(tv%S(i,j,k) - tv%S(i,j,k-1)), 0.0) + ! if (is_root_pe()) write(*,*)'depth_3d(i,j,:)',depth_3d(i,j,:) + + call cvmix_init_bkgnd(max_nlev=nz, & + zw = depth_2d(i,:), & !< interface depth, must bepositive. + bl1 = CS%Bryan_Lewis_c1, & + bl2 = CS%Bryan_Lewis_c2, & + bl3 = CS%Bryan_Lewis_c3, & + bl4 = CS%Bryan_Lewis_c4, & + prandtl = CS%prandtl_turb) + + call cvmix_coeffs_bkgnd(Mdiff_out=CS%kv_bkgnd(i,j,:), & + Tdiff_out=CS%kd_bkgnd(i,j,:), & + nlev=nz, & + max_nlev=nz) + + do k=1,nz + ! Update Kd + Kd(i,j,k) = Kd(i,j,k) + 0.5*(CS%kd_bkgnd(i,j,K) + CS%kd_bkgnd(i,j,K+1)) + ! ######## CHECK ############### + ! GMM, we could update Kv here????? + ! Kv(i,j,k) = Kv(i,j,k) + 0.5*(CS%bkgnd_mixing_csp%kv_bkgnd(i,j,K) + & + ! CS%bkgnd_mixing_csp%kv_bkgnd(i,j,K+1)) enddo enddo - else - do K=2,nz ; do i=is,ie - dRho_int(i,K) = GV%Rlay(k) - GV%Rlay(k-1) - enddo ; enddo - endif - ! Set the buoyancy frequencies. - do k=1,nz ; do i=is,ie - N2_lay(i,k) = G_Rho0 * 0.5*(dRho_int(i,K) + dRho_int(i,K+1)) / & - (GV%H_to_m*(h(i,j,k) + H_neglect)) - enddo ; enddo - do i=is,ie ; N2_int(i,1) = 0.0 ; N2_int(i,nz+1) = 0.0 ; enddo - do K=2,nz ; do i=is,ie - N2_int(i,K) = G_Rho0 * dRho_int(i,K) / & - (0.5*GV%H_to_m*(h(i,j,k-1) + h(i,j,k) + H_neglect)) - enddo ; enddo - - ! Find the bottom boundary layer stratification, and use this in the deepest layers. - do i=is,ie - hb(i) = 0.0 ; dRho_bot(i) = 0.0 - z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) - do_i(i) = (G%mask2dT(i,j) > 0.5) - h_amp(i) = 0.0 - enddo - - do k=nz,2,-1 - do_any = .false. - do i=is,ie ; if (do_i(i)) then - dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) - z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above - - hb(i) = hb(i) + dz_int - dRho_bot(i) = dRho_bot(i) + dRho_int(i,K) - - if (z_from_bot(i) > h_amp(i)) then - if (k>2) then - ! Always include at least one full layer. - hb(i) = hb(i) + 0.5*GV%H_to_m*(h(i,j,k-1) + h(i,j,k-2)) - dRho_bot(i) = dRho_bot(i) + dRho_int(i,K-1) - endif - do_i(i) = .false. + elseif ((.not. CS%Bryan_Lewis_diffusivity) .and. (.not.CS%bulkmixedlayer) .and. & + (CS%Kd/= CS%Kdml)) then + I_Hmix = 1.0 / CS%Hmix + do i=is,ie ; depth(i) = 0.0 ; enddo + do k=1,nz ; do i=is,ie + depth_c = depth(i) + 0.5*GV%H_to_m*h(i,j,k) + if (depth_c <= CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = CS%Kdml + elseif (depth_c >= 2.0*CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = CS%Kd_sfc(i,j) else - do_any = .true. + Kd(i,j,k) = ((CS%Kd_sfc(i,j) - CS%Kdml) * I_Hmix) * depth_c + & + (2.0*CS%Kdml - CS%Kd_sfc(i,j)) endif - endif ; enddo - if (.not.do_any) exit - enddo - - do i=is,ie - if (hb(i) > 0.0) then - N2_bot(i) = (G_Rho0 * dRho_bot(i)) / hb(i) - else ; N2_bot(i) = 0.0 ; endif - z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) - do_i(i) = (G%mask2dT(i,j) > 0.5) - enddo - - do k=nz,2,-1 - do_any = .false. - do i=is,ie ; if (do_i(i)) then - dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) - z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above - - N2_int(i,K) = N2_bot(i) - if (k>2) N2_lay(i,k-1) = N2_bot(i) - - if (z_from_bot(i) > h_amp(i)) then - if (k>2) N2_int(i,K-1) = N2_bot(i) - do_i(i) = .false. - else - do_any = .true. - endif - endif ; enddo - if (.not.do_any) exit - enddo - do i=is,ie - if (hb(i) > 0.0) then - N2_bot(i) = (G_Rho0 * dRho_bot(i)) / hb(i) - else ; N2_bot(i) = 0.0 ; endif - z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) - do_i(i) = (G%mask2dT(i,j) > 0.5) - enddo - - do k=nz,2,-1 - do_any = .false. - do i=is,ie ; if (do_i(i)) then - dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) - z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above - - N2_int(i,K) = N2_bot(i) - if (k>2) N2_lay(i,k-1) = N2_bot(i) - - if (z_from_bot(i) > h_amp(i)) then - if (k>2) N2_int(i,K-1) = N2_bot(i) - do_i(i) = .false. - else - do_any = .true. - endif - endif ; enddo - if (.not.do_any) exit - enddo - if (associated(tv%eqn_of_state)) then - do K=1,nz+1 ; do i=is,ie - dRho_int(i,K) = dRho_int_unfilt(i,K) + depth(i) = depth(i) + GV%H_to_m*h(i,j,k) + enddo ; enddo + + elseif (CS%Henyey_IGW_background_new) then + I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. + do k=1,nz ; do i=is,ie + abs_sin = max(epsilon,abs(sin(G%geoLatT(i,j)*deg_to_rad))) + N_2Omega = max(abs_sin,sqrt(N2_lay(i,k))*I_2Omega) + N02_N2 = (CS%N0_2Omega/N_2Omega)**2 + Kd(i,j,k) = max(CS%Kd_min, CS%Kd_sfc(i,j) * & + ((abs_sin * invcosh(N_2Omega/abs_sin)) * I_x30)*N02_N2) + enddo ; enddo + + else + do k=1,nz ; do i=is,ie + Kd(i,j,k) = CS%Kd_sfc(i,j) enddo ; enddo endif -end subroutine find_N2 + ! Update CS%kd_bkgnd + ! GMM, we could update CS%kv_bkgnd here????? + if (.not. CS%Bryan_Lewis_diffusivity) then + do i=is,ie + CS%kd_bkgnd(i,j,1) = 0.0 + CS%kd_bkgnd(i,j,nz+1) = 0.0 + do k=2,nz + ! Update CS%kd_bkgnd + CS%kd_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) + 0.5*(Kd(i,j,K-1) + Kd(i,j,K)) + ! ######## CHECK ############### + ! GMM, we could update CS%kv_bkgnd here????? + enddo + enddo + endif + +end subroutine calculate_bkgnd_mixing !> Reads the parameter "USE_CVMIX_BACKGROUND" and returns state. !! This function allows other modules to know whether this parameterization will diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index c127ad9a72..81babcd2bd 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -22,7 +22,7 @@ module MOM_set_diffusivity use MOM_cvmix_shear, only : calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_cs use MOM_cvmix_shear, only : cvmix_shear_end use MOM_bkgnd_mixing, only : calculate_bkgnd_mixing, bkgnd_mixing_init, bkgnd_mixing_cs -use MOM_bkgnd_mixing, only : bkgnd_mixing_end +use MOM_bkgnd_mixing, only : bkgnd_mixing_end, sfc_bkgnd_mixing use MOM_string_functions, only : uppercase use MOM_thickness_diffuse, only : vert_fill_TS use MOM_variables, only : thermo_var_ptrs, vertvisc_type, p3d @@ -483,8 +483,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! the appropriate place to add a depth-dependent parameterization or ! another explicit parameterization of Kd. - ! GMM, call MOM_bkgnd_mixing_calc here - call calculate_bkgnd_mixing(h, tv, T_f, S_f, fluxes, G, GV, CS%bkgnd_mixing_csp) + ! set surface diffusivities (CS%bkgnd_mixing_csp%Kd_sfc) + call sfc_bkgnd_mixing(G, CS%bkgnd_mixing_csp) ! GMM, fix OMP calls below @@ -496,13 +496,16 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & !$OMP I_x30,abs_sin,N_2Omega,N02_N2,KT_extra, KS_extra, & !$OMP TKE_to_Kd,maxTKE,dissip,kb) do j=js,je - ! Add vertical background diffusivities and viscosities - do k=1,nz ; do i=is,ie - ! diffusivities - Kd(i,j,k) = Kd(i,j,k) + 0.5*(CS%bkgnd_mixing_csp%kd_bkgnd(i,j,K) + CS%bkgnd_mixing_csp%kd_bkgnd(i,j,K+1)) - ! viscosities - ! GMM, will this be done here? - enddo ; enddo + + ! Set up variables related to the stratification. + call find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, N2_lay, N2_int, N2_bot) + + if (associated(dd%N2_3d)) then + do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo + endif + + ! add background mixing + call calculate_bkgnd_mixing(h, tv, N2_lay, Kd, j, G, GV, CS%bkgnd_mixing_csp) ! GMM, the following will go into the MOM_cvmix_double_diffusion module if (CS%double_diffusion) then @@ -624,21 +627,38 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & enddo ! j-loop if (CS%debug) then - call hchksum(Kd,"BBL Kd",G%HI,haloshift=0) + call hchksum(Kd ,"Kd",G%HI,haloshift=0) + + !!if (associated(CS%bkgnd_mixing_csp%kd_bkgnd)) & + !1 call hchksum(CS%bkgnd_mixing_csp%kd_bkgnd, "kd_bkgnd",G%HI,haloshift=0) + + !!if (associated(CS%bkgnd_mixing_csp%kv_bkgnd)) & + !! call hchksum(CS%bkgnd_mixing_csp%kv_bkgnd, "kv_bkgnd",G%HI,haloshift=0) + if (CS%useKappaShear) call hchksum(visc%Kd_turb,"Turbulent Kd",G%HI,haloshift=0) + if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, & G%HI, 0, symmetric=.true.) endif + if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) then call uvchksum("BBL bbl_thick_[uv]", visc%bbl_thick_u, & visc%bbl_thick_v, G%HI, 0, symmetric=.true.) endif + if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) then call uvchksum("Ray_[uv]", visc%Ray_u, visc%Ray_v, G%HI, 0, symmetric=.true.) endif + endif + ! send bkgnd_mixing diagnostics to post_data + if (CS%bkgnd_mixing_csp%id_kd_bkgnd > 0) & + call post_data(CS%bkgnd_mixing_csp%id_kd_bkgnd, CS%bkgnd_mixing_csp%kd_bkgnd, CS%bkgnd_mixing_csp%diag) + if (CS%bkgnd_mixing_csp%id_kv_bkgnd > 0) & + call post_data(CS%bkgnd_mixing_csp%id_kv_bkgnd, CS%bkgnd_mixing_csp%kv_bkgnd, CS%bkgnd_mixing_csp%diag) + if (CS%Kd_add > 0.0) then if (present(Kd_int)) then !$OMP parallel do default(none) shared(is,ie,js,je,nz,Kd_int,CS,Kd) @@ -965,6 +985,160 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, CS, & end subroutine find_TKE_to_Kd +subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & + N2_lay, N2_int, N2_bot) + 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(in) :: tv + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T_f, S_f + type(forcing), intent(in) :: fluxes + integer, intent(in) :: j + type(set_diffusivity_CS), pointer :: CS + real, dimension(SZI_(G),SZK_(G)+1), intent(out) :: dRho_int, N2_int + real, dimension(SZI_(G),SZK_(G)), intent(out) :: N2_lay + real, dimension(SZI_(G)), intent(out) :: N2_bot + + real, dimension(SZI_(G),SZK_(G)+1) :: & + dRho_int_unfilt, & ! unfiltered density differences across interfaces + dRho_dT, & ! partial derivative of density wrt temp (kg m-3 degC-1) + dRho_dS ! partial derivative of density wrt saln (kg m-3 PPT-1) + + real, dimension(SZI_(G)) :: & + pres, & ! pressure at each interface (Pa) + Temp_int, & ! temperature at each interface (degC) + Salin_int, & ! salinity at each interface (PPT) + drho_bot, & + h_amp, & + hb, & + z_from_bot + + real :: Rml_base ! density of the deepest variable density layer + real :: dz_int ! thickness associated with an interface (meter) + real :: G_Rho0 ! gravitation acceleration divided by Bouss reference density (m4 s-2 kg-1) + real :: H_neglect ! negligibly small thickness, in the same units as h. + + logical :: do_i(SZI_(G)), do_any + integer :: i, k, is, ie, nz + + is = G%isc ; ie = G%iec ; nz = G%ke + G_Rho0 = GV%g_Earth / GV%Rho0 + H_neglect = GV%H_subroundoff + + ! Find the (limited) density jump across each interface. + do i=is,ie + dRho_int(i,1) = 0.0 ; dRho_int(i,nz+1) = 0.0 + dRho_int_unfilt(i,1) = 0.0 ; dRho_int_unfilt(i,nz+1) = 0.0 + enddo + if (associated(tv%eqn_of_state)) then + if (associated(fluxes%p_surf)) then + do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo + else + do i=is,ie ; pres(i) = 0.0 ; enddo + endif + 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) + T_f(i,j,k-1)) + Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1)) + enddo + call calculate_density_derivs(Temp_int, Salin_int, pres, & + dRho_dT(:,K), dRho_dS(:,K), is, ie-is+1, tv%eqn_of_state) + do i=is,ie + dRho_int(i,K) = max(dRho_dT(i,K)*(T_f(i,j,k) - T_f(i,j,k-1)) + & + dRho_dS(i,K)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0) + dRho_int_unfilt(i,K) = max(dRho_dT(i,K)*(tv%T(i,j,k) - tv%T(i,j,k-1)) + & + dRho_dS(i,K)*(tv%S(i,j,k) - tv%S(i,j,k-1)), 0.0) + enddo + enddo + else + do K=2,nz ; do i=is,ie + dRho_int(i,K) = GV%Rlay(k) - GV%Rlay(k-1) + enddo ; enddo + endif + + ! Set the buoyancy frequencies. + do k=1,nz ; do i=is,ie + N2_lay(i,k) = G_Rho0 * 0.5*(dRho_int(i,K) + dRho_int(i,K+1)) / & + (GV%H_to_m*(h(i,j,k) + H_neglect)) + enddo ; enddo + do i=is,ie ; N2_int(i,1) = 0.0 ; N2_int(i,nz+1) = 0.0 ; enddo + do K=2,nz ; do i=is,ie + N2_int(i,K) = G_Rho0 * dRho_int(i,K) / & + (0.5*GV%H_to_m*(h(i,j,k-1) + h(i,j,k) + H_neglect)) + enddo ; enddo + + ! Find the bottom boundary layer stratification, and use this in the deepest layers. + do i=is,ie + hb(i) = 0.0 ; dRho_bot(i) = 0.0 + z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) + do_i(i) = (G%mask2dT(i,j) > 0.5) + + if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation) then + h_amp(i) = sqrt(CS%h2(i,j)) ! for computing Nb + else + h_amp(i) = 0.0 + endif + enddo + + do k=nz,2,-1 + do_any = .false. + do i=is,ie ; if (do_i(i)) then + dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) + z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above + + hb(i) = hb(i) + dz_int + dRho_bot(i) = dRho_bot(i) + dRho_int(i,K) + + if (z_from_bot(i) > h_amp(i)) then + if (k>2) then + ! Always include at least one full layer. + hb(i) = hb(i) + 0.5*GV%H_to_m*(h(i,j,k-1) + h(i,j,k-2)) + dRho_bot(i) = dRho_bot(i) + dRho_int(i,K-1) + endif + do_i(i) = .false. + else + do_any = .true. + endif + endif ; enddo + if (.not.do_any) exit + enddo + + do i=is,ie + if (hb(i) > 0.0) then + N2_bot(i) = (G_Rho0 * dRho_bot(i)) / hb(i) + else ; N2_bot(i) = 0.0 ; endif + z_from_bot(i) = 0.5*GV%H_to_m*h(i,j,nz) + do_i(i) = (G%mask2dT(i,j) > 0.5) + enddo + + do k=nz,2,-1 + do_any = .false. + do i=is,ie ; if (do_i(i)) then + dz_int = 0.5*GV%H_to_m*(h(i,j,k) + h(i,j,k-1)) + z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above + + N2_int(i,K) = N2_bot(i) + if (k>2) N2_lay(i,k-1) = N2_bot(i) + + if (z_from_bot(i) > h_amp(i)) then + if (k>2) N2_int(i,K-1) = N2_bot(i) + do_i(i) = .false. + else + do_any = .true. + endif + endif ; enddo + if (.not.do_any) exit + enddo + + if (associated(tv%eqn_of_state)) then + do K=1,nz+1 ; do i=is,ie + dRho_int(i,K) = dRho_int_unfilt(i,K) + enddo ; enddo + endif + +end subroutine find_N2 + ! GMM, the following will be moved to a new module !> This subroutine sets the additional diffusivities of temperature and From 4cac1e17d9024b42d475264d1526ceeefeac6cca Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 30 Mar 2018 09:49:11 -0600 Subject: [PATCH 25/40] Renamed variables (KH to kd; KM to kv) and added diagnostics --- .../vertical/MOM_cvmix_shear.F90 | 68 +++++++++++-------- 1 file changed, 41 insertions(+), 27 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_shear.F90 b/src/parameterizations/vertical/MOM_cvmix_shear.F90 index 4a91524c76..b99382974d 100644 --- a/src/parameterizations/vertical/MOM_cvmix_shear.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_shear.F90 @@ -36,12 +36,12 @@ module MOM_cvmix_shear 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(:,:,:) :: km !< vertical viscosity at interface (m2/s) - real, allocatable, dimension(:,:,:) :: kh !< vertical diffusivity at interface (m2/s) +! 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() - integer :: id_N2 = -1, id_S2 = -1, id_ri_grad = -1, id_km = -1, id_kh = -1 + integer :: id_N2 = -1, id_S2 = -1, id_ri_grad = -1, id_kv = -1, id_kd = -1 end type cvmix_shear_cs @@ -50,17 +50,17 @@ module MOM_cvmix_shear contains !> Subroutine for calculating (internal) vertical diffusivities/viscosities -subroutine calculate_cvmix_shear(u_H, v_H, h, tv, KH, & - KM, G, GV, CS ) +subroutine calculate_cvmix_shear(u_H, v_H, h, tv, kd, & + kv, G, GV, CS ) type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: u_H !< Initial zonal velocity on T points, in m s-1. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: v_H !< Initial meridional velocity on T points, in m s-1. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: KH !< The vertical diffusivity at each interface + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd !< The vertical diffusivity at each interface !! (not layer!) in m2 s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: KM !< The vertical viscosity at each interface + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv !< The vertical viscosity at each interface !! (not layer!) in m2 s-1. type(cvmix_shear_cs), pointer :: CS !< The control structure returned by a previous call to !! CVMix_shear_init. @@ -116,17 +116,30 @@ subroutine calculate_cvmix_shear(u_H, v_H, h, tv, KH, & N2 = DRHO/DZ S2 = (DU*DU+DV*DV)/(DZ*DZ) Ri_Grad(k) = max(0.,N2)/max(S2,1.e-16) + + ! fill 3d arrays, if user asks for diagsnostics + if (CS%id_N2 > 0) CS%N2(i,j,k) = N2 + if (CS%id_S2 > 0) CS%S2(i,j,k) = S2 + if (CS%id_ri_grad > 0) CS%ri_grad(i,j,k) = Ri_Grad(k) + enddo ! Call to CVMix wrapper for computing interior mixing coefficients. - call cvmix_coeffs_shear(Mdiff_out=KM(i,j,:), & - Tdiff_out=KH(i,j,:), & + call cvmix_coeffs_shear(Mdiff_out=kv(i,j,:), & + Tdiff_out=kd(i,j,:), & RICH=Ri_Grad, & nlev=G%ke, & max_nlev=G%ke) enddo enddo + ! write diagnostics + if (CS%id_kd > 0) call post_data(CS%id_kd,kd, CS%diag) + if (CS%id_kv > 0) call post_data(CS%id_kv,kv, CS%diag) + if (CS%id_N2 > 0) call post_data(CS%id_N2,CS%N2, CS%diag) + if (CS%id_S2 > 0) call post_data(CS%id_S2,CS%S2, CS%diag) + if (CS%id_ri_grad > 0) call post_data(CS%id_ri_grad,CS%ri_grad, CS%diag) + end subroutine calculate_cvmix_shear @@ -200,26 +213,29 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) KPP_nu_zero=CS%Nu_Zero, & KPP_Ri_zero=CS%Ri_zero, & KPP_exp=CS%KPP_exp) - ! Allocation and initialization - allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%N2(:,:,:) = 0. - allocate( CS%S2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%S2(:,:,:) = 0. - !Initialize w/ large Richardson value - allocate( CS%ri_grad( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%ri_grad(:,:,:) = 1.e8 - ! Register diagnostics + ! Register diagnostics; allocation and initialization CS%diag => diag - CS%id_N2 = register_diag_field('ocean_model', 'shear_N2', diag%axesTi, Time, & + + CS%id_N2 = register_diag_field('ocean_model', 'N2_shear', diag%axesTi, Time, & 'Square of Brunt-Vaisala frequency used by MOM_cvmix_shear module', '1/s2') - CS%id_S2 = register_diag_field('ocean_model', 'shear_S2', diag%axesTi, Time, & + if (CS%id_N2 > 0) & + allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%N2(:,:,:) = 0. + + CS%id_S2 = register_diag_field('ocean_model', 'S2_shear', diag%axesTi, Time, & 'Square of vertical shear used by MOM_cvmix_shear module','1/s2') - CS%id_ri_grad = register_diag_field('ocean_model', 'shear_ri_grad', diag%axesTi, Time, & + if (CS%id_S2 > 0) & + allocate( CS%S2( SZI_(G), SZJ_(G), SZK_(G)+1 ) );CS%S2(:,:,:) = 0. + + CS%id_ri_grad = register_diag_field('ocean_model', 'ri_grad_shear', diag%axesTi, Time, & 'Gradient Richarson number used by MOM_cvmix_shear module','nondim') - CS%id_kh = register_diag_field('ocean_model', 'shear_KH', diag%axesTi, Time, & + if (CS%id_ri_grad > 0) & !Initialize w/ large Richardson value + allocate( CS%ri_grad( SZI_(G), SZJ_(G), SZK_(G)+1 ));CS%ri_grad(:,:,:) = 1.e8 + + CS%id_kd = register_diag_field('ocean_model', 'kd_shear', diag%axesTi, Time, & 'Vertical diffusivity added by MOM_cvmix_shear module', 'm2/s') - if (CS%id_kh > 0) allocate(CS%kh(SZI_(G), SZJ_(G), SZK_(G)+1)) - CS%id_km = register_diag_field('ocean_model', 'shear_KM', diag%axesTi, Time, & + CS%id_kv = register_diag_field('ocean_model', 'kv_shear', diag%axesTi, Time, & 'Vertical viscosity added by MOM_cvmix_shear module', 'm2/s') - if (CS%id_km > 0) allocate(CS%km(SZI_(G), SZJ_(G), SZK_(G)+1)) end function cvmix_shear_init @@ -241,11 +257,9 @@ end function cvmix_shear_is_used subroutine cvmix_shear_end(CS) type(cvmix_shear_cs), pointer :: CS ! Control structure - deallocate(CS%N2) - deallocate(CS%S2) - deallocate(CS%ri_grad) - if (CS%id_kh > 0) deallocate(CS%kh, CS%diag) - if (CS%id_km > 0) deallocate(CS%km, CS%diag) + if (CS%id_N2 > 0) deallocate(CS%N2) + if (CS%id_S2 > 0) deallocate(CS%S2) + if (CS%id_ri_grad > 0) deallocate(CS%ri_grad) deallocate(CS) end subroutine cvmix_shear_end From 69ab99895da933072fd3b9d9d0d7f2df751c8481 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 30 Mar 2018 10:18:16 -0600 Subject: [PATCH 26/40] Added a unique Prandtl # for convection and background --- .../vertical/MOM_bkgnd_mixing.F90 | 11 +++++++---- .../vertical/MOM_cvmix_conv.F90 | 17 +++++++++-------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index ba7711586c..d249cec5a8 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -46,7 +46,8 @@ module MOM_bkgnd_mixing real :: N0_2Omega !< ratio of the typical Buoyancy frequency to !! twice the Earth's rotation period, used with the !! Henyey scaling from the mixing - real :: prandtl_turb !< Turbulent Prandtl number + real :: prandtl_bkgnd !< Turbulent Prandtl number used to convert + !! vertical background diffusivity into viscosity real :: Kd_tanh_lat_scale !< A nondimensional scaling for the range of !! diffusivities with Kd_tanh_lat_fn. Valid values !! are in the range of -2 to 2; 0.4 reproduces CM2M. @@ -155,8 +156,10 @@ subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) - call get_param(param_file, mdl, "PRANDTL_TURB", CS%prandtl_turb, & - units="nondim", default=1.0, do_not_log=.true.) + call get_param(param_file, mdl, "PRANDTL_BKGND", CS%prandtl_bkgnd, & + "Turbulent Prandtl number used to convert vertical \n"//& + "background diffusivities into viscosities.", & + units="nondim", default=1.0) ! call openParameterBlock(param_file,'MOM_BACKGROUND_MIXING') @@ -349,7 +352,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, j, G, GV, CS) bl2 = CS%Bryan_Lewis_c2, & bl3 = CS%Bryan_Lewis_c3, & bl4 = CS%Bryan_Lewis_c4, & - prandtl = CS%prandtl_turb) + prandtl = CS%prandtl_bkgnd) call cvmix_coeffs_bkgnd(Mdiff_out=CS%kv_bkgnd(i,j,:), & Tdiff_out=CS%kd_bkgnd(i,j,:), & diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index 721c35faa2..38ea8e6851 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -59,8 +59,8 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) type(cvmix_conv_cs), pointer :: CS !< This module's control structure. ! Local variables - real :: prandtl_turb - logical :: useEPBL + real :: prandtl_conv !< Turbulent Prandtl number used in convective instabilities. + logical :: useEPBL !< If True, use the ePBL boundary layer scheme. ! This include declares and sets the variable "version". #include "version_variable.h" @@ -95,16 +95,17 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) 'as convective mixing might occur in the boundary layer.') endif - call get_param(param_file, mdl, "PRANDTL_TURB", Prandtl_turb, & - "The turbulent Prandtl number applied to shear/conv. \n"//& - "instabilities.", units="nondim", default=1.0, do_not_log=.true.) - call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, do_not_log=.True.) call openParameterBlock(param_file,'CVMIX_CONVECTION') + call get_param(param_file, mdl, "PRANDTL_CONV", prandtl_conv, & + "The turbulent Prandtl number applied to convective \n"//& + "instabilities (i.e., used to convert KD_CONV into KV_CONV)", & + units="nondim", default=1.0) + call get_param(param_file, mdl, 'KD_CONV', CS%kd_conv, & "Diffusivity used in convective regime. Corresponding viscosity \n" // & "(KV_CONV) will be set to KD_CONV * PRANDTL_TURB.", & @@ -117,8 +118,8 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) call closeParameterBlock(param_file) - ! set kv_conv based on kd_conv and Prandtl_turb - CS%kv_conv = CS%kd_conv * Prandtl_turb + ! set kv_conv based on kd_conv and prandtl_conv + CS%kv_conv = CS%kd_conv * prandtl_conv ! allocate arrays and set them to zero allocate(CS%N2(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%N2(:,:,:) = 0. From 6a2e864f3732e7eafcb2f51b25ac7696ad5fffcd Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 30 Mar 2018 10:49:33 -0600 Subject: [PATCH 27/40] Replaced Kd_turb -> Kd_shear and Kv_turb -> Kv_shear --- src/core/MOM_variables.F90 | 10 +++---- .../vertical/MOM_diabatic_driver.F90 | 24 ++++++++--------- .../vertical/MOM_kappa_shear.F90 | 6 ++--- .../vertical/MOM_set_diffusivity.F90 | 26 +++++++++---------- .../vertical/MOM_set_viscosity.F90 | 16 ++++++------ 5 files changed, 41 insertions(+), 41 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index f0f3437f4b..4b85111a21 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -173,7 +173,7 @@ module MOM_variables !! coefficients, and related fields. type, public :: vertvisc_type real :: Prandtl_turb !< The Prandtl number for the turbulent diffusion - !! that is captured in Kd_turb. + !! that is captured in Kd_shear. real, pointer, dimension(:,:) :: & bbl_thick_u => NULL(), & !< The bottom boundary layer thickness at the !! u-points, in m. @@ -224,10 +224,10 @@ module MOM_variables ! Kd_extra_S is positive for salt fingering; Kd_extra_T ! is positive for double diffusive convection. These ! are only allocated if DOUBLE_DIFFUSION is true. - Kd_turb => NULL(), &!< The turbulent diapycnal diffusivity at the interfaces - !! between each layer, in m2 s-1. - Kv_turb => NULL(), &!< The turbulent vertical viscosity at the interfaces - !! between each layer, in m2 s-1. + Kd_shear => NULL(), &!< The shear-driven turbulent diapycnal diffusivity + !! at the interfaces between each layer, in m2 s-1. + Kv_shear => NULL(), &!< The shear-driven turbulent vertical viscosity + !! at the interfaces between each layer, in m2 s-1. TKE_turb => NULL() !< The turbulent kinetic energy per unit mass defined !! at the interfaces between each layer, in m2 s-2. end type vertvisc_type diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index cc0e9e7501..de8c01ac7f 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -582,8 +582,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G call cpu_clock_begin(id_clock_set_diffusivity) ! Sets: Kd, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S - ! Also changes: visc%Kd_turb, visc%TKE_turb (not clear that TKE_turb is used as input ????) - ! And sets visc%Kv_turb + ! Also changes: visc%Kd_shear, visc%TKE_turb (not clear that TKE_turb is used as input ???? + ! And sets visc%Kv_shear call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, CS%set_diff_CSp, Kd, Kd_int) call cpu_clock_end(id_clock_set_diffusivity) if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)") @@ -632,7 +632,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G !$OMP end parallel call KPP_calculate(CS%KPP_CSp, G, GV, h, tv%T, tv%S, u, v, tv%eqn_of_state, & - fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_turb, CS%KPP_NLTheat, CS%KPP_NLTscalar) + fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar) !$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat) if (associated(Hml)) then @@ -678,9 +678,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G !!!!!!!! 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_3d(i,j,k) - ! GMM, I am not sure if Kv_turb is the right place to add kv_conv_3d - visc%Kv_turb(i,j,k) = visc%Kv_turb(i,j,k) + CS%cvmix_conv_csp%kv_conv_3d(i,j,k) - visc%Kd_turb(i,j,k) = visc%Kd_turb(i,j,k) + CS%cvmix_conv_csp%kd_conv_3d(i,j,k) + ! GMM, I am not sure if Kv_shear is the right place to add kv_conv_3d + visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + CS%cvmix_conv_csp%kv_conv_3d(i,j,k) + visc%Kd_shear(i,j,k) = visc%Kd_shear(i,j,k) + CS%cvmix_conv_csp%kd_conv_3d(i,j,k) enddo ; enddo ; enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -829,10 +829,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G if (CS%ePBL_is_additive) then Kd_add_here = Kd_ePBL(i,j,K) - visc%Kv_turb(i,j,K) = visc%Kv_turb(i,j,K) + Kd_ePBL(i,j,K) + visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + Kd_ePBL(i,j,K) else - Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_turb(i,j,K), 0.0) - visc%Kv_turb(i,j,K) = max(visc%Kv_turb(i,j,K), Kd_ePBL(i,j,K)) + Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_shear(i,j,K), 0.0) + visc%Kv_shear(i,j,K) = max(visc%Kv_shear(i,j,K), Kd_ePBL(i,j,K)) endif Ent_int = Kd_add_here * (GV%m_to_H**2 * dt) / & (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) @@ -1365,9 +1365,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G call create_group_pass(CS%pass_hold_eb_ea, eb, G%Domain, dir_flag, halo=1) call create_group_pass(CS%pass_hold_eb_ea, ea, G%Domain, dir_flag, halo=1) call do_group_pass(CS%pass_hold_eb_ea, G%Domain) - ! visc%Kv_turb is not in the group pass because it has larger vertical extent. - if (associated(visc%Kv_turb)) & - call pass_var(visc%Kv_turb, G%Domain, To_All+Omit_Corners, halo=1) + ! visc%Kv_shear is not in the group pass because it has larger vertical extent. + if (associated(visc%Kv_shear)) & + call pass_var(visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) call cpu_clock_end(id_clock_pass) if (.not. CS%useALEalgorithm) then diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index ce32e25b49..6794d7b45b 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -71,7 +71,7 @@ module MOM_kappa_shear real :: TKE_bg ! The background level of TKE, in m2 s-2. real :: kappa_0 ! The background diapycnal diffusivity, in m2 s-1. real :: kappa_tol_err ! The fractional error in kappa that is tolerated. - real :: Prandtl_turb ! Prandtl number used to convert Kd_turb into viscosity. + real :: Prandtl_turb ! Prandtl number used to convert Kd_shear into viscosity. integer :: nkml ! The number of layers in the mixed layer, as ! treated in this routine. If the pieces of the ! mixed layer are not to be treated collectively, @@ -130,7 +130,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & intent(inout) :: kv_io !< The vertical viscosity at each interface !! (not layer!) in m2 s-1. This discards any !! previous value i.e. intent(out) and simply - !! sets Kv = Prandtl * Kd_turb + !! sets Kv = Prandtl * Kd_shear real, intent(in) :: dt !< Time increment, in s. type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous !! call to kappa_shear_init. @@ -156,7 +156,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & ! the iteration toward convergence. ! (in/out) kv_io - The vertical viscosity at each interface ! (not layer!) in m2 s-1. This discards any previous value -! i.e. intent(out) and simply sets Kv = Prandtl * Kd_turb +! i.e. intent(out) and simply sets Kv = Prandtl * Kd_shear ! (in) dt - Time increment, in s. ! (in) G - The ocean's grid structure. ! (in) GV - The ocean's vertical grid structure. diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 81babcd2bd..a672faecc1 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -461,22 +461,22 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & call hchksum(v_h, "before calc_KS v_h",G%HI) endif call cpu_clock_begin(id_clock_kappaShear) - ! Changes: visc%Kd_turb, visc%TKE_turb (not clear that TKE_turb is used as input ????) - ! Sets visc%Kv_turb - call calculate_kappa_shear(u_h, v_h, h, tv, fluxes%p_surf, visc%Kd_turb, visc%TKE_turb, & - visc%Kv_turb, dt, G, GV, CS%kappaShear_CSp) + ! Changes: visc%Kd_shear, visc%TKE_turb (not clear that TKE_turb is used as input ????) + ! Sets visc%Kv_shear + call calculate_kappa_shear(u_h, v_h, h, tv, fluxes%p_surf, visc%Kd_shear, visc%TKE_turb, & + visc%Kv_shear, dt, G, GV, CS%kappaShear_CSp) call cpu_clock_end(id_clock_kappaShear) if (CS%debug) then - call hchksum(visc%Kd_turb, "after calc_KS visc%Kd_turb",G%HI) - call hchksum(visc%Kv_turb, "after calc_KS visc%Kv_turb",G%HI) + call hchksum(visc%Kd_shear, "after calc_KS visc%Kd_shear",G%HI) + call hchksum(visc%Kv_shear, "after calc_KS visc%Kv_shear",G%HI) call hchksum(visc%TKE_turb, "after calc_KS visc%TKE_turb",G%HI) endif if (showCallTree) call callTree_waypoint("done with calculate_kappa_shear (set_diffusivity)") elseif (CS%use_cvmix_shear) then !NOTE{BGR}: this needs to be cleaned up. It works in 1D case, but has not been tested outside. - call calculate_cvmix_shear(u_h, v_h, h, tv, visc%Kd_turb, visc%Kv_turb,G,GV,CS%cvmix_shear_csp) - elseif (associated(visc%Kv_turb)) then - visc%Kv_turb(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled + call calculate_cvmix_shear(u_h, v_h, h, tv, visc%Kd_shear, visc%Kv_shear,G,GV,CS%cvmix_shear_csp) + elseif (associated(visc%Kv_shear)) then + visc%Kv_shear(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled endif ! Calculate the diffusivity, Kd, for each layer. This would be @@ -539,15 +539,15 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%useKappaShear .or. CS%use_cvmix_shear) then if (present(Kd_int)) then do K=2,nz ; do i=is,ie - Kd_int(i,j,K) = visc%Kd_turb(i,j,K) + 0.5*(Kd(i,j,k-1) + Kd(i,j,k)) + Kd_int(i,j,K) = visc%Kd_shear(i,j,K) + 0.5*(Kd(i,j,k-1) + Kd(i,j,k)) enddo ; enddo do i=is,ie - Kd_int(i,j,1) = visc%Kd_turb(i,j,1) ! This isn't actually used. It could be 0. + Kd_int(i,j,1) = visc%Kd_shear(i,j,1) ! This isn't actually used. It could be 0. Kd_int(i,j,nz+1) = 0.0 enddo endif do k=1,nz ; do i=is,ie - Kd(i,j,k) = Kd(i,j,k) + 0.5*(visc%Kd_turb(i,j,K) + visc%Kd_turb(i,j,K+1)) + Kd(i,j,k) = Kd(i,j,k) + 0.5*(visc%Kd_shear(i,j,K) + visc%Kd_shear(i,j,K+1)) enddo ; enddo else if (present(Kd_int)) then @@ -635,7 +635,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & !!if (associated(CS%bkgnd_mixing_csp%kv_bkgnd)) & !! call hchksum(CS%bkgnd_mixing_csp%kv_bkgnd, "kv_bkgnd",G%HI,haloshift=0) - if (CS%useKappaShear) call hchksum(visc%Kd_turb,"Turbulent Kd",G%HI,haloshift=0) + if (CS%useKappaShear) call hchksum(visc%Kd_shear,"Turbulent Kd",G%HI,haloshift=0) if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, & diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 60c96fbb37..8ab4bafb75 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1808,20 +1808,20 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) endif if (use_kappa_shear .or. useKPP .or. useEPBL .or. use_cvmix_shear .or. use_cvmix_conv) then - allocate(visc%Kd_turb(isd:ied,jsd:jed,nz+1)) ; visc%Kd_turb(:,:,:) = 0.0 + allocate(visc%Kd_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kd_shear(:,:,:) = 0.0 allocate(visc%TKE_turb(isd:ied,jsd:jed,nz+1)) ; visc%TKE_turb(:,:,:) = 0.0 - allocate(visc%Kv_turb(isd:ied,jsd:jed,nz+1)) ; visc%Kv_turb(:,:,:) = 0.0 + allocate(visc%Kv_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kv_shear(:,:,:) = 0.0 - vd = var_desc("Kd_turb","m2 s-1","Turbulent diffusivity at interfaces", & + vd = var_desc("Kd_shear","m2 s-1","Shear-driven turbulent diffusivity at interfaces", & hor_grid='h', z_grid='i') - call register_restart_field(visc%Kd_turb, vd, .false., restart_CS) + call register_restart_field(visc%Kd_shear, vd, .false., restart_CS) vd = var_desc("TKE_turb","m2 s-2","Turbulent kinetic energy per unit mass at interfaces", & hor_grid='h', z_grid='i') call register_restart_field(visc%TKE_turb, vd, .false., restart_CS) - vd = var_desc("Kv_turb","m2 s-1","Turbulent viscosity at interfaces", & + vd = var_desc("Kv_shear","m2 s-1","Shear-driven turbulent viscosity at interfaces", & hor_grid='h', z_grid='i') - call register_restart_field(visc%Kv_turb, vd, .false., restart_CS) + call register_restart_field(visc%Kv_shear, vd, .false., restart_CS) endif ! visc%MLD is used to communicate the state of the (e)PBL to the rest of the model @@ -2092,9 +2092,9 @@ subroutine set_visc_end(visc, CS) if (CS%dynamic_viscous_ML) then deallocate(visc%nkml_visc_u) ; deallocate(visc%nkml_visc_v) endif - if (associated(visc%Kd_turb)) deallocate(visc%Kd_turb) + if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear) if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb) - if (associated(visc%Kv_turb)) deallocate(visc%Kv_turb) + if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear) if (associated(visc%ustar_bbl)) deallocate(visc%ustar_bbl) if (associated(visc%TKE_bbl)) deallocate(visc%TKE_bbl) if (associated(visc%taux_shelf)) deallocate(visc%taux_shelf) From 644f200d9ac682c0348a60f800acec0e65ce0b2f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 30 Mar 2018 11:05:41 -0600 Subject: [PATCH 28/40] Replaced Kv_turb -> Kv_shear --- .../vertical/MOM_vert_friction.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index fee1fb456a..ff14a698ed 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1094,21 +1094,21 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m endif endif ; enddo - if (associated(visc%Kv_turb)) then - ! BGR/ Add factor of 2. * the averaged Kv_turb. + if (associated(visc%Kv_shear)) then + ! BGR/ Add factor of 2. * the averaged Kv_shear. ! this is needed to reproduce the analytical solution to ! a simple diffusion problem, likely due to h_shear being ! equal to 2 x \delta z if (work_on_u) then do K=2,nz ; do i=is,ie ; if (do_i(i)) then - Kv_add(i,K) = (2.*0.5)*(visc%Kv_turb(i,j,k) + visc%Kv_turb(i+1,j,k)) + Kv_add(i,K) = (2.*0.5)*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k)) endif ; enddo ; enddo if (do_OBCs) then do I=is,ie ; if (do_i(I) .and. (OBC%segnum_u(I,j) /= OBC_NONE)) then if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then - do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_turb(i,j,k) ; enddo + do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_shear(i,j,k) ; enddo elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then - do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_turb(i+1,j,k) ; enddo + do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_shear(i+1,j,k) ; enddo endif endif ; enddo endif @@ -1117,14 +1117,14 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m endif ; enddo ; enddo else do K=2,nz ; do i=is,ie ; if (do_i(i)) then - Kv_add(i,K) = (2.*0.5)*(visc%Kv_turb(i,j,k) + visc%Kv_turb(i,j+1,k)) + Kv_add(i,K) = (2.*0.5)*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i,j+1,k)) endif ; enddo ; enddo if (do_OBCs) then do i=is,ie ; if (do_i(i) .and. (OBC%segnum_v(i,J) /= OBC_NONE)) then if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then - do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_turb(i,j,k) ; enddo + do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_shear(i,j,k) ; enddo elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then - do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_turb(i,j+1,k) ; enddo + do K=2,nz ; Kv_add(i,K) = 2.*visc%Kv_shear(i,j+1,k) ; enddo endif endif ; enddo endif From 80c33a469bfb5111e3f2a8182ab504d10fa8bce0 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 30 Mar 2018 11:10:49 -0600 Subject: [PATCH 29/40] Replaced Kv_turb -> Kv_shear --- src/core/MOM.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index bcf6516505..22dbb86b15 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1448,7 +1448,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & real, dimension(:,:), pointer :: shelf_area type(MOM_restart_CS), pointer :: restart_CSp_tmp => NULL() type(group_pass_type) :: tmp_pass_uv_T_S_h, pass_uv_T_S_h - type(group_pass_type) :: tmp_pass_Kv_turb + ! GMM, the following *is not* used. Should we delete it? + type(group_pass_type) :: tmp_pass_Kv_shear real :: default_val ! default value for a parameter logical :: write_geom_files ! If true, write out the grid geometry files. @@ -2288,8 +2289,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call do_group_pass(pass_uv_T_S_h, G%Domain) - if (associated(CS%visc%Kv_turb)) & - call pass_var(CS%visc%Kv_turb, G%Domain, To_All+Omit_Corners, halo=1) + if (associated(CS%visc%Kv_shear)) & + call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) call cpu_clock_end(id_clock_pass_init) From 801a86a06a978af6ddde612f5d1b3aab6d8b2cf5 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 2 Apr 2018 16:38:41 -0600 Subject: [PATCH 30/40] Added variable Kv_slow in visc type --- src/core/MOM_variables.F90 | 3 +++ src/parameterizations/vertical/MOM_set_viscosity.F90 | 1 + 2 files changed, 4 insertions(+) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 4b85111a21..f7fa45f12c 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -228,6 +228,9 @@ module MOM_variables !! at the interfaces between each layer, in m2 s-1. Kv_shear => NULL(), &!< The shear-driven turbulent vertical viscosity !! at the interfaces between each layer, in m2 s-1. + Kv_slow => NULL(), &!< The turbulent vertical viscosity component due to + !! "slow" processes (e.g., tidal, background, + !! convection etc). TKE_turb => NULL() !< The turbulent kinetic energy per unit mass defined !! at the interfaces between each layer, in m2 s-2. end type vertvisc_type diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 8ab4bafb75..19f956f1b1 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1811,6 +1811,7 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) allocate(visc%Kd_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kd_shear(:,:,:) = 0.0 allocate(visc%TKE_turb(isd:ied,jsd:jed,nz+1)) ; visc%TKE_turb(:,:,:) = 0.0 allocate(visc%Kv_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kv_shear(:,:,:) = 0.0 + ! GMM, allocate visc%Kv_slow here? vd = var_desc("Kd_shear","m2 s-1","Shear-driven turbulent diffusivity at interfaces", & hor_grid='h', z_grid='i') From e54429028de1f94833c81c1f27abb1e27cdf5747 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 3 Apr 2018 09:13:38 -0600 Subject: [PATCH 31/40] Allocated/deallocate visc%Kv_slow --- src/parameterizations/vertical/MOM_set_viscosity.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 19f956f1b1..e38d484714 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1811,7 +1811,7 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) allocate(visc%Kd_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kd_shear(:,:,:) = 0.0 allocate(visc%TKE_turb(isd:ied,jsd:jed,nz+1)) ; visc%TKE_turb(:,:,:) = 0.0 allocate(visc%Kv_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kv_shear(:,:,:) = 0.0 - ! GMM, allocate visc%Kv_slow here? + 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') @@ -1823,6 +1823,10 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) vd = var_desc("Kv_shear","m2 s-1","Shear-driven turbulent viscosity at interfaces", & hor_grid='h', z_grid='i') call register_restart_field(visc%Kv_shear, vd, .false., restart_CS) + vd = var_desc("Kv_slow","m2 s-1","Vertical turbulent viscosity at interfaces due to slow" & + " processes", hor_grid='h', z_grid='i') + call register_restart_field(visc%Kv_slow, vd, .false., restart_CS) + endif ! visc%MLD is used to communicate the state of the (e)PBL to the rest of the model @@ -2094,6 +2098,7 @@ subroutine set_visc_end(visc, CS) deallocate(visc%nkml_visc_u) ; deallocate(visc%nkml_visc_v) endif if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear) + if (associated(visc%Kv_slow)) deallocate(visc%Kv_slow) if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb) if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear) if (associated(visc%ustar_bbl)) deallocate(visc%ustar_bbl) From 105094f2affebdaeffbd914486c53c00c1da5d20 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 3 Apr 2018 09:15:12 -0600 Subject: [PATCH 32/40] Add viscosities due to convection into visc%Kv_slow --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index de8c01ac7f..c1ae3fe7e7 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -671,16 +671,14 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G endif ! endif for KPP - ! Add diffusivity due to convection (computed via CVMix) + ! 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_3d(i,j,k) - ! GMM, I am not sure if Kv_shear is the right place to add kv_conv_3d - visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + CS%cvmix_conv_csp%kv_conv_3d(i,j,k) - visc%Kd_shear(i,j,k) = visc%Kd_shear(i,j,k) + CS%cvmix_conv_csp%kd_conv_3d(i,j,k) + visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + CS%cvmix_conv_csp%kv_conv_3d(i,j,k) enddo ; enddo ; enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From 37ca38fdcdfb845484631700c7355237c5568167 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 3 Apr 2018 09:55:16 -0600 Subject: [PATCH 33/40] Renamed kd_bkgnd_3d -> kd_bkgnd; kv_bkgnd_3d -> kv_bkgnd --- .../vertical/MOM_cvmix_conv.F90 | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index 38ea8e6851..ea50bd1a16 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -39,8 +39,8 @@ module MOM_cvmix_conv ! Diagnostics arrays real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) - real, allocatable, dimension(:,:,:) :: kd_conv_3d !< Diffusivity added by convection (m2/s) - real, allocatable, dimension(:,:,:) :: kv_conv_3d !< Viscosity added by convection (m2/s) + real, allocatable, dimension(:,:,:) :: kd_conv !< Diffusivity added by convection (m2/s) + real, allocatable, dimension(:,:,:) :: kv_conv !< Viscosity added by convection (m2/s) end type cvmix_conv_cs @@ -123,8 +123,8 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) ! allocate arrays and set them to zero allocate(CS%N2(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%N2(:,:,:) = 0. - allocate(CS%kd_conv_3d(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kd_conv_3d(:,:,:) = 0. - allocate(CS%kv_conv_3d(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kv_conv_3d(:,:,:) = 0. + allocate(CS%kd_conv(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kd_conv(:,:,:) = 0. + allocate(CS%kv_conv(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%kv_conv(:,:,:) = 0. ! Register diagnostics CS%diag => diag @@ -214,8 +214,8 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, CS, hbl) kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) - call cvmix_coeffs_conv(Mdiff_out=CS%kv_conv_3d(i,j,:), & - Tdiff_out=CS%kd_conv_3d(i,j,:), & + call cvmix_coeffs_conv(Mdiff_out=CS%kv_conv(i,j,:), & + Tdiff_out=CS%kd_conv(i,j,:), & Nsqr=CS%N2(i,j,:), & dens=rho_1d(:), & dens_lwr=rho_lwr(:), & @@ -225,8 +225,8 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, CS, hbl) ! Do not apply mixing due to convection within the boundary layer do k=1,NINT(hbl(i,j)) - CS%kv_conv_3d(i,j,k) = 0.0 - CS%kd_conv_3d(i,j,k) = 0.0 + CS%kv_conv(i,j,k) = 0.0 + CS%kd_conv(i,j,k) = 0.0 enddo enddo @@ -234,14 +234,14 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, CS, hbl) if (CS%debug) then call hchksum(CS%N2, "MOM_cvmix_conv: N2",G%HI,haloshift=0) - call hchksum(CS%kd_conv_3d, "MOM_cvmix_conv: kd_conv_3d",G%HI,haloshift=0) - call hchksum(CS%kv_conv_3d, "MOM_cvmix_conv: kv_conv_3d",G%HI,haloshift=0) + call hchksum(CS%kd_conv, "MOM_cvmix_conv: kd_conv",G%HI,haloshift=0) + call hchksum(CS%kv_conv, "MOM_cvmix_conv: kv_conv",G%HI,haloshift=0) endif ! send diagnostics to post_data if (CS%id_N2 > 0) call post_data(CS%id_N2, CS%N2, CS%diag) - if (CS%id_kd_conv > 0) call post_data(CS%id_kd_conv, CS%kd_conv_3d, CS%diag) - if (CS%id_kv_conv > 0) call post_data(CS%id_kv_conv, CS%kv_conv_3d, CS%diag) + if (CS%id_kd_conv > 0) call post_data(CS%id_kd_conv, CS%kd_conv, CS%diag) + if (CS%id_kv_conv > 0) call post_data(CS%id_kv_conv, CS%kv_conv, CS%diag) end subroutine calculate_cvmix_conv @@ -260,8 +260,8 @@ subroutine cvmix_conv_end(CS) type(cvmix_conv_cs), pointer :: CS ! Control structure deallocate(CS%N2) - deallocate(CS%kd_conv_3d) - deallocate(CS%kv_conv_3d) + deallocate(CS%kd_conv) + deallocate(CS%kv_conv) deallocate(CS) end subroutine cvmix_conv_end From a2498a115a4596387b378d1edcd2efe54a59f510 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 3 Apr 2018 10:26:10 -0600 Subject: [PATCH 34/40] Changed name in register_diag_field kd_shear -> kd_shear_cvmix; kv_shear -> kv_shear_cvmix --- src/parameterizations/vertical/MOM_cvmix_shear.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_shear.F90 b/src/parameterizations/vertical/MOM_cvmix_shear.F90 index b99382974d..345522126b 100644 --- a/src/parameterizations/vertical/MOM_cvmix_shear.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_shear.F90 @@ -232,9 +232,9 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) if (CS%id_ri_grad > 0) & !Initialize w/ large Richardson value allocate( CS%ri_grad( SZI_(G), SZJ_(G), SZK_(G)+1 ));CS%ri_grad(:,:,:) = 1.e8 - CS%id_kd = register_diag_field('ocean_model', 'kd_shear', diag%axesTi, Time, & + CS%id_kd = register_diag_field('ocean_model', 'kd_shear_cvmix', diag%axesTi, Time, & 'Vertical diffusivity added by MOM_cvmix_shear module', 'm2/s') - CS%id_kv = register_diag_field('ocean_model', 'kv_shear', diag%axesTi, Time, & + CS%id_kv = register_diag_field('ocean_model', 'kv_shear_cvmix', diag%axesTi, Time, & 'Vertical viscosity added by MOM_cvmix_shear module', 'm2/s') end function cvmix_shear_init From c1acec8033095e17a16487f50f6fef036333aedd Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 3 Apr 2018 10:28:12 -0600 Subject: [PATCH 35/40] Add vertical background viscosity into visc%Kv_slow --- .../vertical/MOM_bkgnd_mixing.F90 | 56 ++++++++++--------- .../vertical/MOM_diabatic_driver.F90 | 4 +- .../vertical/MOM_set_diffusivity.F90 | 2 +- .../vertical/MOM_set_viscosity.F90 | 4 +- 4 files changed, 36 insertions(+), 30 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index d249cec5a8..2e4bba1110 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -18,6 +18,7 @@ module MOM_bkgnd_mixing use MOM_verticalGrid, only : verticalGrid_type use MOM_file_parser, only : get_param, log_version, param_file_type use cvmix_background, only : cvmix_init_bkgnd, cvmix_coeffs_bkgnd +use MOM_variables, only : vertvisc_type use MOM_intrinsic_functions, only : invcosh implicit none ; private @@ -302,23 +303,26 @@ end subroutine sfc_bkgnd_mixing !> Calculates the vertical background diffusivities/viscosities -subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, j, G, GV, CS) +subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc, j, G, GV, CS) - type(ocean_grid_type), intent(in) :: G !< Grid structure. - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. - type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. real, dimension(SZI_(G),SZK_(G)), intent(in) :: N2_lay!< squared buoyancy frequency associated - !! with layers (1/s2) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd !< Diapycnal diffusivity of each layer (m2/sec). - integer, intent(in) :: j - type(bkgnd_mixing_cs), pointer :: CS !< The control structure returned by - !! a previous call to bkgnd_mixing_init. + !! with layers (1/s2) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd !< Diapycnal diffusivity of each layer (m2/sec). + type(vertvisc_type), intent(inout) :: visc!< Structure containing vertical viscosities, + !! bottom boundary layer properies, and related + !! fields. + integer, intent(in) :: j !< Meridional grid indice. + type(bkgnd_mixing_cs), pointer :: CS !< The control structure returned by + !! a previous call to bkgnd_mixing_init. ! local variables - real, dimension(SZI_(G), SZK_(G)+1) :: depth_2d !< distance from surface of an interface (m) + real, dimension(SZI_(G), SZK_(G)+1) :: depth_2d !< distance from surface of an interface (m) real, dimension(SZI_(G)) :: & - depth !< distance from surface of an interface (meter) + depth !< distance from surface of an interface (meter) real :: depth_c !< depth of the center of a layer (meter) real :: I_Hmix !< inverse of fixed mixed layer thickness (1/m) real :: I_2Omega !< 1/(2 Omega) (sec) @@ -359,15 +363,11 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, j, G, GV, CS) nlev=nz, & max_nlev=nz) - do k=1,nz ! Update Kd + do k=1,nz Kd(i,j,k) = Kd(i,j,k) + 0.5*(CS%kd_bkgnd(i,j,K) + CS%kd_bkgnd(i,j,K+1)) - ! ######## CHECK ############### - ! GMM, we could update Kv here????? - ! Kv(i,j,k) = Kv(i,j,k) + 0.5*(CS%bkgnd_mixing_csp%kv_bkgnd(i,j,K) + & - ! CS%bkgnd_mixing_csp%kv_bkgnd(i,j,K+1)) enddo - enddo + enddo ! i loop elseif ((.not. CS%Bryan_Lewis_diffusivity) .and. (.not.CS%bulkmixedlayer) .and. & (CS%Kd/= CS%Kdml)) then @@ -401,17 +401,23 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, j, G, GV, CS) enddo ; enddo endif - ! Update CS%kd_bkgnd - ! GMM, we could update CS%kv_bkgnd here????? + ! Update CS%kd_bkgnd and CS%kv_bkgnd for diagnostic purposes if (.not. CS%Bryan_Lewis_diffusivity) then do i=is,ie - CS%kd_bkgnd(i,j,1) = 0.0 - CS%kd_bkgnd(i,j,nz+1) = 0.0 + 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 - ! Update CS%kd_bkgnd CS%kd_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) + 0.5*(Kd(i,j,K-1) + Kd(i,j,K)) - ! ######## CHECK ############### - ! GMM, we could update CS%kv_bkgnd here????? + CS%kv_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) * CS%prandtl_bkgnd + enddo + enddo + endif + + ! Update visc%Kv_slow, if associated + if (associated(visc%Kv_slow)) then + do i=is,ie + do k=1,nz+1 + visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + CS%kv_bkgnd(i,j,k) enddo enddo endif diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index c1ae3fe7e7..25c2464b56 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -677,8 +677,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G !!!!!!!! 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_3d(i,j,k) - visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + CS%cvmix_conv_csp%kv_conv_3d(i,j,k) + 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index a672faecc1..a66a94334b 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -505,7 +505,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif ! add background mixing - call calculate_bkgnd_mixing(h, tv, N2_lay, Kd, j, G, GV, CS%bkgnd_mixing_csp) + call calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc, j, G, GV, CS%bkgnd_mixing_csp) ! GMM, the following will go into the MOM_cvmix_double_diffusion module if (CS%double_diffusion) then diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index e38d484714..3df3e7b780 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1823,8 +1823,8 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) vd = var_desc("Kv_shear","m2 s-1","Shear-driven turbulent viscosity at interfaces", & hor_grid='h', z_grid='i') call register_restart_field(visc%Kv_shear, vd, .false., restart_CS) - vd = var_desc("Kv_slow","m2 s-1","Vertical turbulent viscosity at interfaces due to slow" & - " processes", hor_grid='h', z_grid='i') + vd = var_desc("Kv_slow","m2 s-1","Vertical turbulent viscosity at interfaces due \n" // & + " to slow processes", hor_grid='h', z_grid='i') call register_restart_field(visc%Kv_slow, vd, .false., restart_CS) endif From 648018e4afa58b13fa22b736daeb9c163de6e0bb Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 3 Apr 2018 10:29:04 -0600 Subject: [PATCH 36/40] Rename kv_conv -> kv_conv_const; kd_conv -> kd_conv_const --- src/parameterizations/vertical/MOM_cvmix_conv.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index ea50bd1a16..55e7d55d6e 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -26,8 +26,8 @@ module MOM_cvmix_conv type, public :: cvmix_conv_cs ! Parameters - real :: kd_conv !< diffusivity constant used in convective regime (m2/s) - real :: kv_conv !< viscosity constant used in convective regime (m2/s) + real :: kd_conv_const !< diffusivity constant used in convective regime (m2/s) + real :: kv_conv_const !< viscosity constant used in convective regime (m2/s) real :: bv_sqr_conv !< Threshold for squared buoyancy frequency !! needed to trigger Brunt-Vaisala parameterization (1/s^2) real :: min_thickness !< Minimum thickness allowed (m) @@ -106,7 +106,7 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) "instabilities (i.e., used to convert KD_CONV into KV_CONV)", & units="nondim", default=1.0) - call get_param(param_file, mdl, 'KD_CONV', CS%kd_conv, & + call get_param(param_file, mdl, 'KD_CONV', CS%kd_conv_const, & "Diffusivity used in convective regime. Corresponding viscosity \n" // & "(KV_CONV) will be set to KD_CONV * PRANDTL_TURB.", & units='m2/s', default=1.00) @@ -118,8 +118,8 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) call closeParameterBlock(param_file) - ! set kv_conv based on kd_conv and prandtl_conv - CS%kv_conv = CS%kd_conv * prandtl_conv + ! set kv_conv_const based on kd_conv_const and prandtl_conv + CS%kv_conv_const = CS%kd_conv_const * prandtl_conv ! allocate arrays and set them to zero allocate(CS%N2(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%N2(:,:,:) = 0. @@ -135,8 +135,8 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) CS%id_kv_conv = register_diag_field('ocean_model', 'conv_kv', diag%axesTi, Time, & 'Additional viscosity added by MOM_cvmix_conv module', 'm2/s') - call cvmix_init_conv(convect_diff=CS%kd_conv, & - convect_visc=CS%kv_conv, & + call cvmix_init_conv(convect_diff=CS%kd_conv_const, & + convect_visc=CS%kv_conv_const, & lBruntVaisala=.true., & BVsqr_convect=CS%bv_sqr_conv) From 1218973de6e7e167c1e0545a1d6685acdf195653 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 3 Apr 2018 13:46:17 -0600 Subject: [PATCH 37/40] Raname variables Kd -> kd_lay; visc%Kv_slow - > kv --- .../vertical/MOM_bkgnd_mixing.F90 | 31 +++++++++---------- .../vertical/MOM_set_diffusivity.F90 | 8 +---- 2 files changed, 15 insertions(+), 24 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 2e4bba1110..763b3665d4 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -303,7 +303,7 @@ end subroutine sfc_bkgnd_mixing !> Calculates the vertical background diffusivities/viscosities -subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc, j, G, GV, CS) +subroutine calculate_bkgnd_mixing(h, tv, N2_lay, kd_lay, kv, j, G, GV, CS) type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. @@ -311,10 +311,9 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc, j, G, GV, CS) type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. real, dimension(SZI_(G),SZK_(G)), intent(in) :: N2_lay!< squared buoyancy frequency associated !! with layers (1/s2) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd !< Diapycnal diffusivity of each layer (m2/sec). - type(vertvisc_type), intent(inout) :: visc!< Structure containing vertical viscosities, - !! bottom boundary layer properies, and related - !! fields. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: kd_lay!< Diapycnal diffusivity of each layer m2 s-1. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv !< The "slow" vertical viscosity at each interface + !! (not layer!) in m2 s-1. integer, intent(in) :: j !< Meridional grid indice. type(bkgnd_mixing_cs), pointer :: CS !< The control structure returned by !! a previous call to bkgnd_mixing_init. @@ -365,7 +364,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc, j, G, GV, CS) ! Update Kd do k=1,nz - Kd(i,j,k) = Kd(i,j,k) + 0.5*(CS%kd_bkgnd(i,j,K) + CS%kd_bkgnd(i,j,K+1)) + kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5*(CS%kd_bkgnd(i,j,K) + CS%kd_bkgnd(i,j,K+1)) enddo enddo ! i loop @@ -378,7 +377,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc, j, G, GV, CS) if (depth_c <= CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = CS%Kdml elseif (depth_c >= 2.0*CS%Hmix) then ; CS%kd_bkgnd(i,j,k) = CS%Kd_sfc(i,j) else - Kd(i,j,k) = ((CS%Kd_sfc(i,j) - CS%Kdml) * I_Hmix) * depth_c + & + kd_lay(i,j,k) = ((CS%Kd_sfc(i,j) - CS%Kdml) * I_Hmix) * depth_c + & (2.0*CS%Kdml - CS%Kd_sfc(i,j)) endif @@ -391,13 +390,13 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc, j, G, GV, CS) abs_sin = max(epsilon,abs(sin(G%geoLatT(i,j)*deg_to_rad))) N_2Omega = max(abs_sin,sqrt(N2_lay(i,k))*I_2Omega) N02_N2 = (CS%N0_2Omega/N_2Omega)**2 - Kd(i,j,k) = max(CS%Kd_min, CS%Kd_sfc(i,j) * & + kd_lay(i,j,k) = max(CS%Kd_min, CS%Kd_sfc(i,j) * & ((abs_sin * invcosh(N_2Omega/abs_sin)) * I_x30)*N02_N2) enddo ; enddo else do k=1,nz ; do i=is,ie - Kd(i,j,k) = CS%Kd_sfc(i,j) + kd_lay(i,j,k) = CS%Kd_sfc(i,j) enddo ; enddo endif @@ -407,20 +406,18 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc, 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(i,j,K-1) + Kd(i,j,K)) + CS%kd_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) + 0.5*(kd_lay(i,j,K-1) + kd_lay(i,j,K)) CS%kv_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) * CS%prandtl_bkgnd enddo enddo endif - ! Update visc%Kv_slow, if associated - if (associated(visc%Kv_slow)) then - do i=is,ie - do k=1,nz+1 - visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + CS%kv_bkgnd(i,j,k) - enddo + ! Update kv + do i=is,ie + do k=1,nz+1 + kv(i,j,k) = kv(i,j,k) + CS%kv_bkgnd(i,j,k) enddo - endif + enddo end subroutine calculate_bkgnd_mixing diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index a66a94334b..08c438bf24 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -505,7 +505,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif ! add background mixing - call calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc, j, G, GV, CS%bkgnd_mixing_csp) + 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 @@ -629,12 +629,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%debug) then call hchksum(Kd ,"Kd",G%HI,haloshift=0) - !!if (associated(CS%bkgnd_mixing_csp%kd_bkgnd)) & - !1 call hchksum(CS%bkgnd_mixing_csp%kd_bkgnd, "kd_bkgnd",G%HI,haloshift=0) - - !!if (associated(CS%bkgnd_mixing_csp%kv_bkgnd)) & - !! call hchksum(CS%bkgnd_mixing_csp%kv_bkgnd, "kv_bkgnd",G%HI,haloshift=0) - if (CS%useKappaShear) call hchksum(visc%Kd_shear,"Turbulent Kd",G%HI,haloshift=0) if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then From 0ce9980a4c69aba51a1d08c63458a6d2cb77fd96 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 4 Apr 2018 14:06:37 -0600 Subject: [PATCH 38/40] By pass reading BULKMIXEDLAYER via get_param BULKMIXEDLAYER is not always defined in MOM_input (e.g., in CM2G63L), so adding new code where get_param is used to determine if BULKMIXEDLAYER is used won't always work. To by pass that the following is now used: CS%bulkmixedlayer = (GV%nkml > 0) --- src/parameterizations/vertical/MOM_bkgnd_mixing.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 763b3665d4..91f20e4ab2 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -133,8 +133,10 @@ subroutine bkgnd_mixing_init(Time, G, GV, param_file, diag, CS) units="m2 s-1", default=0.01*CS%Kd) ! The following is needed to set one of the choices of vertical background mixing - call get_param(param_file, mdl, "BULKMIXEDLAYER", CS%bulkmixedlayer, & - do_not_log=.true.) + + ! BULKMIXEDLAYER is not always defined (e.g., CM2G63L), so the following by pass + ! the need to include BULKMIXEDLAYER in MOM_input + CS%bulkmixedlayer = (GV%nkml > 0) if (CS%bulkmixedlayer) then ! Check that Kdml is not set when using bulk mixed layer call get_param(param_file, mdl, "KDML", CS%Kdml, default=-1.) From 36f7a0652120552c523a0dde6353343e9dc509bf Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 4 Apr 2018 15:17:52 -0600 Subject: [PATCH 39/40] Changed kv to a pointer since it might not always be associated --- src/parameterizations/vertical/MOM_bkgnd_mixing.F90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 91f20e4ab2..14c6c3412e 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -314,7 +314,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, kd_lay, kv, j, G, GV, CS) real, dimension(SZI_(G),SZK_(G)), intent(in) :: N2_lay!< squared buoyancy frequency associated !! with layers (1/s2) real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: kd_lay!< Diapycnal diffusivity of each layer m2 s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv !< The "slow" vertical viscosity at each interface + real, dimension(:,:,:), pointer :: kv !< The "slow" vertical viscosity at each interface !! (not layer!) in m2 s-1. integer, intent(in) :: j !< Meridional grid indice. type(bkgnd_mixing_cs), pointer :: CS !< The control structure returned by @@ -415,11 +415,13 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, kd_lay, kv, j, G, GV, CS) endif ! Update kv - do i=is,ie - do k=1,nz+1 - kv(i,j,k) = kv(i,j,k) + CS%kv_bkgnd(i,j,k) + if (associated(kv)) then + do i=is,ie + do k=1,nz+1 + kv(i,j,k) = kv(i,j,k) + CS%kv_bkgnd(i,j,k) + enddo enddo - enddo + endif end subroutine calculate_bkgnd_mixing From 4dde6037bb65e7301948d9008283c9d36b3602b0 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 4 Apr 2018 15:19:53 -0600 Subject: [PATCH 40/40] Avoided a possible seg. fault in set_diffusivity_end --- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 08c438bf24..eb12eb23d6 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -2957,6 +2957,8 @@ end subroutine set_diffusivity_init subroutine set_diffusivity_end(CS) type(set_diffusivity_CS), pointer :: CS !< Control structure for this module + if (.not.associated(CS)) return + call bkgnd_mixing_end(CS%bkgnd_mixing_csp) if (CS%user_change_diff) &