MOM_CVMix_conv.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Interface to CVMix convection scheme.
7
8use mom_debugging, only : hchksum
10use mom_diag_mediator, only : post_data
11use mom_eos, only : calculate_density
12use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
14use mom_file_parser, only : get_param, log_version, param_file_type
15use mom_grid, only : ocean_grid_type
16use mom_interface_heights, only : thickness_to_dz
20use cvmix_convection, only : cvmix_init_conv, cvmix_coeffs_conv
21use cvmix_kpp, only : cvmix_kpp_compute_kobl_depth
22
23implicit none ; private
24
25#include <MOM_memory.h>
26
28
29!> Control structure including parameters for CVMix convection.
30type, public :: cvmix_conv_cs ; private
31
32 ! Parameters
33 real :: kd_conv_const !< diffusivity constant used in convective regime [Z2 T-1 ~> m2 s-1]
34 real :: kv_conv_const !< viscosity constant used in convective regime [Z2 T-1 ~> m2 s-1]
35 real :: bv_sqr_conv !< Threshold for squared buoyancy frequency
36 !! needed to trigger Brunt-Vaisala parameterization [T-2 ~> s-2]
37 real :: min_thickness !< Minimum thickness allowed [Z ~> m]
38 logical :: debug !< If true, turn on debugging
39
40 ! Diagnostic handles and pointers
41 type(diag_ctrl), pointer :: diag => null() !< Pointer to diagnostics control structure
42 !>@{ Diagnostics handles
43 integer :: id_n2 = -1, id_kd_conv = -1, id_kv_conv = -1
44 !>@}
45
46end type cvmix_conv_cs
47
48character(len=40) :: mdl = "MOM_CVMix_conv" !< This module's name.
49
50contains
51
52!> Initialized the CVMix convection mixing routine.
53logical function cvmix_conv_init(Time, G, GV, US, param_file, diag, CS)
54
55 type(time_type), intent(in) :: time !< The current time.
56 type(ocean_grid_type), intent(in) :: g !< Grid structure.
57 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
58 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
59 type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle
60 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure.
61 type(cvmix_conv_cs), intent(inout) :: cs !< CVMix convection control structure
62
63 real :: prandtl_conv !< Turbulent Prandtl number used in convective instabilities [nondim]
64 logical :: useepbl !< If True, use the ePBL boundary layer scheme.
65
66 ! This include declares and sets the variable "version".
67# include "version_variable.h"
68
69 ! Read parameters
70 call get_param(param_file, mdl, "USE_CVMix_CONVECTION", cvmix_conv_init, default=.false., do_not_log=.true.)
71 call log_version(param_file, mdl, version, &
72 "Parameterization of enhanced mixing due to convection via CVMix", &
73 all_default=.not.cvmix_conv_init)
74 call get_param(param_file, mdl, "USE_CVMix_CONVECTION", cvmix_conv_init, &
75 "If true, turns on the enhanced mixing due to convection "//&
76 "via CVMix. This scheme increases diapycnal diffs./viscs. "//&
77 "at statically unstable interfaces. Relevant parameters are "//&
78 "contained in the CVMix_CONVECTION% parameter block.", &
79 default=.false.)
80
81 if (.not. cvmix_conv_init) return
82
83 call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useepbl, default=.false., &
84 do_not_log=.true.)
85
86 ! Warn user if EPBL is being used, since in this case mixing due to convection will
87 ! be aplied in the boundary layer
88 if (useepbl) then
89 call mom_error(warning, 'MOM_CVMix_conv_init: '// &
90 'CVMix convection may not be properly applied when ENERGETICS_SFC_PBL = True '//&
91 'as convective mixing might occur in the boundary layer.')
92 endif
93
94 call get_param(param_file, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)
95
96 call get_param(param_file, mdl, 'MIN_THICKNESS', cs%min_thickness, &
97 units="m", scale=us%m_to_Z, default=0.001, do_not_log=.true.)
98
99 call openparameterblock(param_file,'CVMix_CONVECTION')
100
101 call get_param(param_file, mdl, "PRANDTL_CONV", prandtl_conv, &
102 "The turbulent Prandtl number applied to convective "//&
103 "instabilities (i.e., used to convert KD_CONV into KV_CONV)", &
104 units="nondim", default=1.0)
105
106 call get_param(param_file, mdl, 'KD_CONV', cs%kd_conv_const, &
107 "Diffusivity used in convective regime. Corresponding viscosity "//&
108 "(KV_CONV) will be set to KD_CONV * PRANDTL_CONV.", &
109 units='m2/s', default=1.00, scale=us%m2_s_to_Z2_T)
110
111 call get_param(param_file, mdl, 'BV_SQR_CONV', cs%bv_sqr_conv, &
112 "Threshold for squared buoyancy frequency needed to trigger "//&
113 "Brunt-Vaisala parameterization.", &
114 units='1/s^2', default=0.0, scale=us%T_to_s**2)
115
116 call closeparameterblock(param_file)
117
118 ! set kv_conv_const based on kd_conv_const and prandtl_conv
119 cs%kv_conv_const = cs%kd_conv_const * prandtl_conv
120
121 ! Register diagnostics
122 cs%diag => diag
123 cs%id_N2 = register_diag_field('ocean_model', 'N2_conv', diag%axesTi, time, &
124 'Square of Brunt-Vaisala frequency used by MOM_CVMix_conv module', '1/s2', conversion=us%s_to_T**2)
125 cs%id_kd_conv = register_diag_field('ocean_model', 'kd_conv', diag%axesTi, time, &
126 'Additional diffusivity added by MOM_CVMix_conv module', 'm2/s', conversion=us%Z2_T_to_m2_s)
127 cs%id_kv_conv = register_diag_field('ocean_model', 'kv_conv', diag%axesTi, time, &
128 'Additional viscosity added by MOM_CVMix_conv module', 'm2/s', conversion=us%Z2_T_to_m2_s)
129
130 call cvmix_init_conv(convect_diff=us%Z2_T_to_m2_s*cs%kd_conv_const, &
131 convect_visc=us%Z2_T_to_m2_s*cs%kv_conv_const, &
132 lbruntvaisala=.true., &
133 bvsqr_convect=us%s_to_T**2*cs%bv_sqr_conv)
134
135end function cvmix_conv_init
136
137!> Subroutine for calculating enhanced diffusivity/viscosity
138!! due to convection via CVMix
139subroutine calculate_cvmix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux)
140
141 type(ocean_grid_type), intent(in) :: g !< Grid structure.
142 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
143 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
144 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
145 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
146 type(cvmix_conv_cs), intent(in) :: cs !< CVMix convection control structure
147 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hbl !< Depth of ocean boundary layer [Z ~> m]
148 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
149 intent(inout) :: kd !< Diapycnal diffusivity at each interface
150 !! that will be incremented here
151 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
152 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
153 intent(inout) :: kv !< Viscosity at each interface that will be
154 !! incremented here [H Z T-1 ~> m2 s-1 or Pa s]
155 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
156 optional, intent(inout) :: kd_aux !< A second diapycnal diffusivity at each
157 !! interface that will also be incremented
158 !! here [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
159
160 ! local variables
161 real, dimension(SZK_(GV)) :: rho_lwr !< Adiabatic Water Density [kg m-3], this is a dummy
162 !! variable since here convection is always
163 !! computed based on Brunt Vaisala.
164 real, dimension(SZK_(GV)) :: rho_1d !< water density in a column [kg m-3], this is also
165 !! a dummy variable, same reason as above.
166 real, dimension(SZK_(GV)+1) :: n2 !< Squared buoyancy frequency [s-2]
167 real, dimension(SZK_(GV)+1) :: kv_col !< Viscosities at interfaces in the column [m2 s-1]
168 real, dimension(SZK_(GV)+1) :: kd_col !< Diffusivities at interfaces in the column [m2 s-1]
169 real, dimension(SZK_(GV)+1) :: ifaceheight !< Height of interfaces [Z ~> m]
170 real, dimension(SZK_(GV)) :: cellheight !< Height of cell centers [Z ~> m]
171 real, dimension(SZI_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m]
172 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
173 kd_conv, & !< Diffusivity added by convection for diagnostics [Z2 T-1 ~> m2 s-1]
174 kv_conv, & !< Viscosity added by convection for diagnostics [Z2 T-1 ~> m2 s-1]
175 n2_3d !< Squared buoyancy frequency for diagnostics [T-2 ~> s-2]
176 integer :: kobl !< level of ocean boundary layer extent
177 real :: g_o_rho0 ! Gravitational acceleration, perhaps divided by density, times unit conversion factors
178 ! [H s-2 R-1 ~> m4 s-2 kg-1 or m s-2]
179 real :: pref ! Interface pressures [R L2 T-2 ~> Pa]
180 real :: rhok, rhokm1 ! In situ densities of the layers above and below at the interface pressure [R ~> kg m-3]
181 real :: dh_int ! The distance between layer centers [H ~> m or kg m-2]
182 real :: dh, hcorr ! Limited thicknesses and a cumulative correction [Z ~> m]
183 integer :: i, j, k
184
185 if (gv%Boussinesq) then
186 g_o_rho0 = (us%s_to_T**2*gv%Z_to_H) * gv%g_Earth_Z_T2 / gv%Rho0
187 else
188 g_o_rho0 = (us%s_to_T**2*gv%RZ_to_H) * gv%g_Earth_Z_T2
189 endif
190
191 ! initialize dummy variables
192 rho_lwr(:) = 0.0 ; rho_1d(:) = 0.0
193
194 ! set N2 to zero at the top- and bottom-most interfaces
195 n2(1) = 0.0 ; n2(gv%ke+1) = 0.0
196
197 if (cs%id_N2 > 0) n2_3d(:,:,:) = 0.0
198 if (cs%id_kv_conv > 0) kv_conv(:,:,:) = 0.0
199 if (cs%id_kd_conv > 0) kd_conv(:,:,:) = 0.0
200
201 do j = g%jsc, g%jec
202
203 ! Find the vertical distances across layers.
204 call thickness_to_dz(h, tv, dz, j, g, gv)
205
206 do i = g%isc, g%iec
207
208 ! skip calling at land points
209 !if (G%mask2dT(i,j) == 0.) cycle
210
211 pref = 0. ; if (associated(tv%p_surf)) pref = tv%p_surf(i,j)
212 ! Compute Brunt-Vaisala frequency (static stability) on interfaces
213 do k=2,gv%ke
214
215 ! pRef is pressure at interface between k and km1 [R L2 T-2 ~> Pa].
216 pref = pref + (gv%H_to_RZ*gv%g_Earth) * h(i,j,k)
217 call calculate_density(tv%t(i,j,k), tv%s(i,j,k), pref, rhok, tv%eqn_of_state)
218 call calculate_density(tv%t(i,j,k-1), tv%s(i,j,k-1), pref, rhokm1, tv%eqn_of_state)
219
220 dh_int = 0.5*(h(i,j,k-1) + h(i,j,k)) + gv%H_subroundoff
221 n2(k) = g_o_rho0 * (rhok - rhokm1) / dh_int ! Can be negative
222
223 enddo
224
225 ifaceheight(1) = 0.0 ! BBL is all relative to the surface
226 hcorr = 0.0
227 ! compute heights at cell center and interfaces
228 do k=1,gv%ke
229 dh = dz(i,k) ! Nominal thickness to use for increment, in the units of heights
230 dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
231 hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
232 dh = max(dh, cs%min_thickness) ! Limited increment dh>=min_thickness
233 cellheight(k) = ifaceheight(k) - 0.5 * dh
234 ifaceheight(k+1) = ifaceheight(k) - dh
235 enddo
236
237 ! gets index of the level and interface above hbl
238 kobl = cvmix_kpp_compute_kobl_depth(ifaceheight, cellheight, hbl(i,j))
239
240 kv_col(:) = 0.0 ; kd_col(:) = 0.0
241 call cvmix_coeffs_conv(mdiff_out=kv_col(:), &
242 tdiff_out=kd_col(:), &
243 nsqr=n2(:), &
244 dens=rho_1d(:), &
245 dens_lwr=rho_lwr(:), &
246 nlev=gv%ke, &
247 max_nlev=gv%ke, &
248 obl_ind=kobl)
249
250 ! Increment the diffusivity outside of the boundary layer.
251 do k=max(1,kobl+1),gv%ke+1
252 kd(i,j,k) = kd(i,j,k) + gv%m2_s_to_HZ_T * kd_col(k)
253 enddo
254 if (present(kd_aux)) then
255 ! Increment the other diffusivity outside of the boundary layer.
256 do k=max(1,kobl+1),gv%ke+1
257 kd_aux(i,j,k) = kd_aux(i,j,k) + gv%m2_s_to_HZ_T * kd_col(k)
258 enddo
259 endif
260
261 ! Increment the viscosity outside of the boundary layer.
262 do k=max(1,kobl+1),gv%ke+1
263 kv(i,j,k) = kv(i,j,k) + gv%m2_s_to_HZ_T * kv_col(k)
264 enddo
265
266 ! Store 3-d arrays for diagnostics.
267 if (cs%id_kv_conv > 0) then
268 ! Do not apply mixing due to convection within the boundary layer
269 do k=max(1,kobl+1),gv%ke+1
270 kv_conv(i,j,k) = us%m2_s_to_Z2_T * kv_col(k)
271 enddo
272 endif
273 if (cs%id_kd_conv > 0) then
274 ! Do not apply mixing due to convection within the boundary layer
275 do k=max(1,kobl+1),gv%ke+1
276 kd_conv(i,j,k) = us%m2_s_to_Z2_T * kd_col(k)
277 enddo
278 endif
279
280 if (cs%id_N2 > 0) then ; do k=2,gv%ke ; n2_3d(i,j,k) = us%T_to_s**2*n2(k) ; enddo ; endif
281
282 enddo
283 enddo
284
285 if (cs%debug) then
286 ! if (CS%id_N2 > 0) call hchksum(N2_3d, "MOM_CVMix_conv: N2", G%HI, haloshift=0, unscale=US%s_to_T**2)
287 ! if (CS%id_kd_conv > 0) &
288 ! call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, unscale=US%Z2_T_to_m2_s)
289 ! if (CS%id_kv_conv > 0) &
290 ! call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, unscale=US%Z2_T_to_m2_s)
291 call hchksum(kd, "MOM_CVMix_conv: Kd", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)
292 call hchksum(kv, "MOM_CVMix_conv: Kv", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)
293 endif
294
295 ! send diagnostics to post_data
296 if (cs%id_N2 > 0) call post_data(cs%id_N2, n2_3d, cs%diag)
297 if (cs%id_kd_conv > 0) call post_data(cs%id_kd_conv, kd_conv, cs%diag)
298 if (cs%id_kv_conv > 0) call post_data(cs%id_kv_conv, kv_conv, cs%diag)
299
300end subroutine calculate_cvmix_conv
301
302!> Reads the parameter "USE_CVMix_CONVECTION" and returns state.
303!! This function allows other modules to know whether this parameterization will
304!! be used without needing to duplicate the log entry.
305logical function cvmix_conv_is_used(param_file)
306 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
307 call get_param(param_file, mdl, "USE_CVMix_CONVECTION", cvmix_conv_is_used, &
308 default=.false., do_not_log=.true.)
309
310end function cvmix_conv_is_used
311
312end module mom_cvmix_conv