MOM6
MOM_CVMix_shear.F90
Go to the documentation of this file.
1 !> Interface to CVMix interior shear schemes
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 !> \author Brandon Reichl
7 
8 use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
9 use mom_diag_mediator, only : diag_ctrl, time_type
10 use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
12 use mom_grid, only : ocean_grid_type
17 use cvmix_shear, only : cvmix_init_shear, cvmix_coeffs_shear
19 implicit none ; private
20 
21 #include <MOM_memory.h>
22 
24 
25 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
26 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
27 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
28 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
29 
30 !> Control structure including parameters for CVMix interior shear schemes.
31 type, public :: cvmix_shear_cs ! TODO: private
32  logical :: use_lmd94 !< Flags to use the LMD94 scheme
33  logical :: use_pp81 !< Flags to use Pacanowski and Philander (JPO 1981)
34  logical :: smooth_ri !< If true, smooth Ri using a 1-2-1 filter
35  real :: ri_zero !< LMD94 critical Richardson number
36  real :: nu_zero !< LMD94 maximum interior diffusivity
37  real :: kpp_exp !< Exponent of unitless factor of diff.
38  !! for KPP internal shear mixing scheme.
39  real, allocatable, dimension(:,:,:) :: n2 !< Squared Brunt-Vaisala frequency [s-2]
40  real, allocatable, dimension(:,:,:) :: s2 !< Squared shear frequency [s-2]
41  real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number
42  real, allocatable, dimension(:,:,:) :: ri_grad_smooth !< Gradient Richardson number
43  !! after smoothing
44  character(10) :: mix_scheme !< Mixing scheme name (string)
45 
46  type(diag_ctrl), pointer :: diag => null() !< Pointer to the diagnostics control structure
47  !>@{ Diagnostic handles
48  integer :: id_n2 = -1, id_s2 = -1, id_ri_grad = -1, id_kv = -1, id_kd = -1
49  integer :: id_ri_grad_smooth = -1
50  !!@}
51 
52 end type cvmix_shear_cs
53 
54 character(len=40) :: mdl = "MOM_CVMix_shear" !< This module's name.
55 
56 contains
57 
58 !> Subroutine for calculating (internal) vertical diffusivities/viscosities
59 subroutine calculate_cvmix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS )
60  type(ocean_grid_type), intent(in) :: g !< Grid structure.
61  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
62  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
63  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: u_h !< Initial zonal velocity on T points [L T-1 ~> m s-1]
64  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: v_h !< Initial meridional velocity on T
65  !! points [L T-1 ~> m s-1]
66  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
67  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
68  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd !< The vertical diffusivity at each interface
69  !! (not layer!) [Z2 T-1 ~> m2 s-1].
70  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv !< The vertical viscosity at each interface
71  !! (not layer!) [Z2 T-1 ~> m2 s-1].
72  type(cvmix_shear_cs), pointer :: cs !< The control structure returned by a previous
73  !! call to CVMix_shear_init.
74  ! Local variables
75  integer :: i, j, k, kk, km1
76  real :: gorho ! Gravitational acceleration divided by density in MKS units [m4 s-2]
77  real :: pref, du, dv, drho, dz, n2, s2, dummy
78  real, dimension(2*(G%ke)) :: pres_1d, temp_1d, salt_1d, rho_1d
79  real, dimension(G%ke+1) :: ri_grad !< Gradient Richardson number
80  real, dimension(G%ke+1) :: kvisc !< Vertical viscosity at interfaces [m2 s-1]
81  real, dimension(G%ke+1) :: kdiff !< Diapycnal diffusivity at interfaces [m2 s-1]
82  real, parameter :: epsln = 1.e-10 !< Threshold to identify vanished layers
83 
84  ! some constants
85  gorho = gv%mks_g_Earth / (us%R_to_kg_m3*gv%Rho0)
86 
87  do j = g%jsc, g%jec
88  do i = g%isc, g%iec
89 
90  ! skip calling for land points
91  if (g%mask2dT(i,j)==0.) cycle
92 
93  ! Richardson number computed for each cell in a column.
94  pref = 0.
95  ri_grad(:)=1.e8 !Initialize w/ large Richardson value
96  do k=1,g%ke
97  ! pressure, temp, and saln for EOS
98  ! kk+1 = k fields
99  ! kk+2 = km1 fields
100  km1 = max(1, k-1)
101  kk = 2*(k-1)
102  pres_1d(kk+1) = pref
103  pres_1d(kk+2) = pref
104  temp_1d(kk+1) = tv%T(i,j,k)
105  temp_1d(kk+2) = tv%T(i,j,km1)
106  salt_1d(kk+1) = tv%S(i,j,k)
107  salt_1d(kk+2) = tv%S(i,j,km1)
108 
109  ! pRef is pressure at interface between k and km1.
110  ! iterate pRef for next pass through k-loop.
111  pref = pref + gv%H_to_Pa * h(i,j,k)
112 
113  enddo ! k-loop finishes
114 
115  ! compute in-situ density
116  call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, 1, 2*g%ke, tv%EQN_OF_STATE)
117 
118  ! N2 (can be negative) on interface
119  do k = 1, g%ke
120  km1 = max(1, k-1)
121  kk = 2*(k-1)
122  du = us%L_T_to_m_s*(u_h(i,j,k) - u_h(i,j,km1))
123  dv = us%L_T_to_m_s*(v_h(i,j,k) - v_h(i,j,km1))
124  drho = (gorho * (rho_1d(kk+1) - rho_1d(kk+2)) )
125  dz = ((0.5*(h(i,j,km1) + h(i,j,k))+gv%H_subroundoff)*gv%H_to_m)
126  n2 = drho/dz
127  s2 = (du*du+dv*dv)/(dz*dz)
128  ri_grad(k) = max(0.,n2)/max(s2,1.e-10)
129 
130  ! fill 3d arrays, if user asks for diagsnostics
131  if (cs%id_N2 > 0) cs%N2(i,j,k) = n2
132  if (cs%id_S2 > 0) cs%S2(i,j,k) = s2
133 
134  enddo
135 
136  ri_grad(g%ke+1) = ri_grad(g%ke)
137 
138  if (cs%id_ri_grad > 0) cs%ri_grad(i,j,:) = ri_grad(:)
139 
140  if (cs%smooth_ri) then
141  ! 1) fill Ri_grad in vanished layers with adjacent value
142  do k = 2, g%ke
143  if (h(i,j,k) .le. epsln) ri_grad(k) = ri_grad(k-1)
144  enddo
145 
146  ri_grad(g%ke+1) = ri_grad(g%ke)
147 
148  ! 2) vertically smooth Ri with 1-2-1 filter
149  dummy = 0.25 * ri_grad(2)
150  ri_grad(g%ke+1) = ri_grad(g%ke)
151  do k = 3, g%ke
152  ri_grad(k) = dummy + 0.5 * ri_grad(k) + 0.25 * ri_grad(k+1)
153  dummy = 0.25 * ri_grad(k)
154  enddo
155 
156  if (cs%id_ri_grad_smooth > 0) cs%ri_grad_smooth(i,j,:) = ri_grad(:)
157  endif
158 
159  do k=1,g%ke+1
160  kvisc(k) = us%Z2_T_to_m2_s * kv(i,j,k)
161  kdiff(k) = us%Z2_T_to_m2_s * kd(i,j,k)
162  enddo
163 
164  ! Call to CVMix wrapper for computing interior mixing coefficients.
165  call cvmix_coeffs_shear(mdiff_out=kvisc(:), &
166  tdiff_out=kdiff(:), &
167  rich=ri_grad(:), &
168  nlev=g%ke, &
169  max_nlev=g%ke)
170  do k=1,g%ke+1
171  kv(i,j,k) = us%m2_s_to_Z2_T * kvisc(k)
172  kd(i,j,k) = us%m2_s_to_Z2_T * kdiff(k)
173  enddo
174  enddo
175  enddo
176 
177  ! write diagnostics
178  if (cs%id_kd > 0) call post_data(cs%id_kd, kd, cs%diag)
179  if (cs%id_kv > 0) call post_data(cs%id_kv, kv, cs%diag)
180  if (cs%id_N2 > 0) call post_data(cs%id_N2, cs%N2, cs%diag)
181  if (cs%id_S2 > 0) call post_data(cs%id_S2, cs%S2, cs%diag)
182  if (cs%id_ri_grad > 0) call post_data(cs%id_ri_grad, cs%ri_grad, cs%diag)
183  if (cs%id_ri_grad_smooth > 0) call post_data(cs%id_ri_grad_smooth ,cs%ri_grad_smooth, cs%diag)
184 
185 end subroutine calculate_cvmix_shear
186 
187 
188 !> Initialized the CVMix internal shear mixing routine.
189 !! \note *This is where we test to make sure multiple internal shear
190 !! mixing routines (including JHL) are not enabled at the same time.
191 !! (returns) CVMix_shear_init - True if module is to be used, False otherwise
192 logical function cvmix_shear_init(Time, G, GV, US, param_file, diag, CS)
193  type(time_type), intent(in) :: time !< The current time.
194  type(ocean_grid_type), intent(in) :: g !< Grid structure.
195  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
196  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
197  type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle
198  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure.
199  type(cvmix_shear_cs), pointer :: cs !< This module's control structure.
200  ! Local variables
201  integer :: numbertrue=0
202  logical :: use_jhl
203 ! This include declares and sets the variable "version".
204 #include "version_variable.h"
205 
206  if (associated(cs)) then
207  call mom_error(warning, "CVMix_shear_init called with an associated "// &
208  "control structure.")
209  return
210  endif
211  allocate(cs)
212 
213 ! Set default, read and log parameters
214  call log_version(param_file, mdl, version, &
215  "Parameterization of shear-driven turbulence via CVMix (various options)")
216  call get_param(param_file, mdl, "USE_LMD94", cs%use_LMD94, &
217  "If true, use the Large-McWilliams-Doney (JGR 1994) "//&
218  "shear mixing parameterization.", default=.false.)
219  if (cs%use_LMD94) then
220  numbertrue=numbertrue + 1
221  cs%Mix_Scheme='KPP'
222  endif
223  call get_param(param_file, mdl, "USE_PP81", cs%use_PP81, &
224  "If true, use the Pacanowski and Philander (JPO 1981) "//&
225  "shear mixing parameterization.", default=.false.)
226  if (cs%use_PP81) then
227  numbertrue = numbertrue + 1
228  cs%Mix_Scheme='PP'
229  endif
230  use_jhl=kappa_shear_is_used(param_file)
231  if (use_jhl) numbertrue = numbertrue + 1
232  ! After testing for interior schemes, make sure only 0 or 1 are enabled.
233  ! Otherwise, warn user and kill job.
234  if ((numbertrue) > 1) then
235  call mom_error(fatal, 'MOM_CVMix_shear_init: '// &
236  'Multiple shear driven internal mixing schemes selected,'//&
237  ' please disable all but one scheme to proceed.')
238  endif
239  cvmix_shear_init=(cs%use_PP81.or.cs%use_LMD94)
240 
241 ! Forego remainder of initialization if not using this scheme
242  if (.not. cvmix_shear_init) return
243  call get_param(param_file, mdl, "NU_ZERO", cs%Nu_Zero, &
244  "Leading coefficient in KPP shear mixing.", &
245  units="nondim", default=5.e-3)
246  call get_param(param_file, mdl, "RI_ZERO", cs%Ri_Zero, &
247  "Critical Richardson for KPP shear mixing, "// &
248  "NOTE this the internal mixing and this is "// &
249  "not for setting the boundary layer depth." &
250  ,units="nondim", default=0.8)
251  call get_param(param_file, mdl, "KPP_EXP", cs%KPP_exp, &
252  "Exponent of unitless factor of diffusivities, "// &
253  "for KPP internal shear mixing scheme." &
254  ,units="nondim", default=3.0)
255  call get_param(param_file, mdl, "SMOOTH_RI", cs%smooth_ri, &
256  "If true, vertically smooth the Richardson "// &
257  "number by applying a 1-2-1 filter once.", &
258  default = .false.)
259  call cvmix_init_shear(mix_scheme=cs%Mix_Scheme, &
260  kpp_nu_zero=cs%Nu_Zero, &
261  kpp_ri_zero=cs%Ri_zero, &
262  kpp_exp=cs%KPP_exp)
263 
264  ! Register diagnostics; allocation and initialization
265  cs%diag => diag
266 
267  cs%id_N2 = register_diag_field('ocean_model', 'N2_shear', diag%axesTi, time, &
268  'Square of Brunt-Vaisala frequency used by MOM_CVMix_shear module', '1/s2')
269  if (cs%id_N2 > 0) then
270  allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) ) ; cs%N2(:,:,:) = 0.
271  endif
272 
273  cs%id_S2 = register_diag_field('ocean_model', 'S2_shear', diag%axesTi, time, &
274  'Square of vertical shear used by MOM_CVMix_shear module','1/s2')
275  if (cs%id_S2 > 0) then
276  allocate( cs%S2( szi_(g), szj_(g), szk_(g)+1 ) ) ; cs%S2(:,:,:) = 0.
277  endif
278 
279  cs%id_ri_grad = register_diag_field('ocean_model', 'ri_grad_shear', diag%axesTi, time, &
280  'Gradient Richarson number used by MOM_CVMix_shear module','nondim')
281  if (cs%id_ri_grad > 0) then !Initialize w/ large Richardson value
282  allocate( cs%ri_grad( szi_(g), szj_(g), szk_(g)+1 )) ; cs%ri_grad(:,:,:) = 1.e8
283  endif
284 
285  cs%id_ri_grad_smooth = register_diag_field('ocean_model', 'ri_grad_shear_smooth', &
286  diag%axesTi, time, &
287  'Smoothed gradient Richarson number used by MOM_CVMix_shear module','nondim')
288  if (cs%id_ri_grad_smooth > 0) then !Initialize w/ large Richardson value
289  allocate( cs%ri_grad_smooth( szi_(g), szj_(g), szk_(g)+1 )) ; cs%ri_grad_smooth(:,:,:) = 1.e8
290  endif
291 
292  cs%id_kd = register_diag_field('ocean_model', 'kd_shear_CVMix', diag%axesTi, time, &
293  'Vertical diffusivity added by MOM_CVMix_shear module', 'm2/s', conversion=us%Z2_T_to_m2_s)
294  cs%id_kv = register_diag_field('ocean_model', 'kv_shear_CVMix', diag%axesTi, time, &
295  'Vertical viscosity added by MOM_CVMix_shear module', 'm2/s', conversion=us%Z2_T_to_m2_s)
296 
297 end function cvmix_shear_init
298 
299 !> Reads the parameters "LMD94" and "PP81" and returns state.
300 !! This function allows other modules to know whether this parameterization will
301 !! be used without needing to duplicate the log entry.
302 logical function cvmix_shear_is_used(param_file)
303  type(param_file_type), intent(in) :: param_file !< Run-time parameter files handle.
304  ! Local variables
305  logical :: lmd94, pp81
306  call get_param(param_file, mdl, "USE_LMD94", lmd94, &
307  default=.false., do_not_log = .true.)
308  call get_param(param_file, mdl, "Use_PP81", pp81, &
309  default=.false., do_not_log = .true.)
310  cvmix_shear_is_used = (lmd94 .or. pp81)
311 end function cvmix_shear_is_used
312 
313 !> Clear pointers and dealocate memory
314 subroutine cvmix_shear_end(CS)
315  type(cvmix_shear_cs), pointer :: cs !< Control structure for this module that
316  !! will be deallocated in this subroutine
317 
318  if (.not. associated(cs)) return
319 
320  if (cs%id_N2 > 0) deallocate(cs%N2)
321  if (cs%id_S2 > 0) deallocate(cs%S2)
322  if (cs%id_ri_grad > 0) deallocate(cs%ri_grad)
323  deallocate(cs)
324 
325 end subroutine cvmix_shear_end
326 
327 end module mom_cvmix_shear
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_cvmix_shear::cvmix_shear_cs
Control structure including parameters for CVMix interior shear schemes.
Definition: MOM_CVMix_shear.F90:31
mom_cvmix_shear::cvmix_shear_init
logical function, public cvmix_shear_init(Time, G, GV, US, param_file, diag, CS)
Initialized the CVMix internal shear mixing routine.
Definition: MOM_CVMix_shear.F90:193
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_cvmix_shear::cvmix_shear_is_used
logical function, public cvmix_shear_is_used(param_file)
Reads the parameters "LMD94" and "PP81" and returns state. This function allows other modules to know...
Definition: MOM_CVMix_shear.F90:303
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_kappa_shear::kappa_shear_is_used
logical function, public kappa_shear_is_used(param_file)
This function indicates to other modules whether the Jackson et al shear mixing parameterization will...
Definition: MOM_kappa_shear.F90:2109
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_cvmix_shear::mdl
character(len=40) mdl
This module's name.
Definition: MOM_CVMix_shear.F90:54
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_cvmix_shear
Interface to CVMix interior shear schemes.
Definition: MOM_CVMix_shear.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_cvmix_shear::cvmix_shear_end
subroutine, public cvmix_shear_end(CS)
Clear pointers and dealocate memory.
Definition: MOM_CVMix_shear.F90:315
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_kappa_shear
Shear-dependent mixing following Jackson et al. 2008.
Definition: MOM_kappa_shear.F90:2
mom_cvmix_shear::calculate_cvmix_shear
subroutine, public calculate_cvmix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS)
Subroutine for calculating (internal) vertical diffusivities/viscosities.
Definition: MOM_CVMix_shear.F90:60
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60