MOM_CVMix_ddiff.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 double diffusion scheme.
7
9use mom_diag_mediator, only : post_data
10use mom_eos, only : calculate_density_derivs
11use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
13use mom_file_parser, only : get_param, log_version, param_file_type
14use mom_debugging, only : hchksum
15use mom_grid, only : ocean_grid_type
19use cvmix_ddiff, only : cvmix_init_ddiff, cvmix_coeffs_ddiff
20use cvmix_kpp, only : cvmix_kpp_compute_kobl_depth
21implicit none ; private
22
23#include <MOM_memory.h>
24
26
27!> Control structure including parameters for CVMix double diffusion.
28type, public :: cvmix_ddiff_cs ; private
29
30 ! Parameters
31 real :: strat_param_max !< maximum value for the stratification parameter [nondim]
32 real :: kappa_ddiff_s !< leading coefficient in formula for salt-fingering regime
33 !! for salinity diffusion [Z2 T-1 ~> m2 s-1]
34 real :: ddiff_exp1 !< interior exponent in salt-fingering regime formula [nondim]
35 real :: ddiff_exp2 !< exterior exponent in salt-fingering regime formula [nondim]
36 real :: mol_diff !< molecular diffusivity [Z2 T-1 ~> m2 s-1]
37 real :: kappa_ddiff_param1 !< exterior coefficient in diffusive convection regime [nondim]
38 real :: kappa_ddiff_param2 !< middle coefficient in diffusive convection regime [nondim]
39 real :: kappa_ddiff_param3 !< interior coefficient in diffusive convection regime [nondim]
40 real :: min_thickness !< Minimum thickness allowed [H ~> m or kg m-2]
41 character(len=4) :: diff_conv_type !< type of diffusive convection to use. Options are Marmorino &
42 !! Caldwell 1976 ("MC76"; default) and Kelley 1988, 1990 ("K90")
43 logical :: debug !< If true, turn on debugging
44
45end type cvmix_ddiff_cs
46
47character(len=40) :: mdl = "MOM_CVMix_ddiff" !< This module's name.
48
49contains
50
51!> Initialized the CVMix double diffusion module.
52logical function cvmix_ddiff_init(Time, G, GV, US, param_file, diag, CS)
53
54 type(time_type), intent(in) :: time !< The current time.
55 type(ocean_grid_type), intent(in) :: g !< Grid structure.
56 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
57 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
58 type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle
59 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure.
60 type(cvmix_ddiff_cs), pointer :: cs !< This module's control structure.
61
62 ! This include declares and sets the variable "version".
63# include "version_variable.h"
64
65 if (associated(cs)) then
66 call mom_error(warning, "CVMix_ddiff_init called with an associated "// &
67 "control structure.")
68 return
69 endif
70
71 ! Read parameters
72 call get_param(param_file, mdl, "USE_CVMIX_DDIFF", cvmix_ddiff_init, default=.false., do_not_log=.true.)
73 call log_version(param_file, mdl, version, &
74 "Parameterization of mixing due to double diffusion processes via CVMix", &
75 all_default=.not.cvmix_ddiff_init)
76 call get_param(param_file, mdl, "USE_CVMIX_DDIFF", cvmix_ddiff_init, &
77 "If true, turns on double diffusive processes via CVMix. "//&
78 "Note that double diffusive processes on viscosity are ignored "//&
79 "in CVMix, see http://cvmix.github.io/ for justification.", &
80 default=.false.)
81
82 if (.not. cvmix_ddiff_init) return
83 allocate(cs)
84
85 call get_param(param_file, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)
86
87 call get_param(param_file, mdl, 'MIN_THICKNESS', cs%min_thickness, &
88 units="m", scale=gv%m_to_H, default=0.001, do_not_log=.true.)
89
90 call openparameterblock(param_file,'CVMIX_DDIFF')
91
92 call get_param(param_file, mdl, "STRAT_PARAM_MAX", cs%strat_param_max, &
93 "The maximum value for the double dissusion stratification parameter", &
94 units="nondim", default=2.55)
95
96 call get_param(param_file, mdl, "KAPPA_DDIFF_S", cs%kappa_ddiff_s, &
97 "Leading coefficient in formula for salt-fingering regime for salinity diffusion.", &
98 units="m2 s-1", default=1.0e-4, scale=us%m2_s_to_Z2_T)
99
100 call get_param(param_file, mdl, "DDIFF_EXP1", cs%ddiff_exp1, &
101 "Interior exponent in salt-fingering regime formula.", &
102 units="nondim", default=1.0)
103
104 call get_param(param_file, mdl, "DDIFF_EXP2", cs%ddiff_exp2, &
105 "Exterior exponent in salt-fingering regime formula.", &
106 units="nondim", default=3.0)
107
108 call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM1", cs%kappa_ddiff_param1, &
109 "Exterior coefficient in diffusive convection regime.", &
110 units="nondim", default=0.909)
111
112 call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM2", cs%kappa_ddiff_param2, &
113 "Middle coefficient in diffusive convection regime.", &
114 units="nondim", default=4.6)
115
116 call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM3", cs%kappa_ddiff_param3, &
117 "Interior coefficient in diffusive convection regime.", &
118 units="nondim", default=-0.54)
119
120 call get_param(param_file, mdl, "MOL_DIFF", cs%mol_diff, &
121 "Molecular diffusivity used in CVMix double diffusion.", &
122 units="m2 s-1", default=1.5e-6, scale=us%m2_s_to_Z2_T)
123
124 call get_param(param_file, mdl, "DIFF_CONV_TYPE", cs%diff_conv_type, &
125 "type of diffusive convection to use. Options are Marmorino \n" //&
126 "and Caldwell 1976 (MC76) and Kelley 1988, 1990 (K90).", &
127 default="MC76")
128
129 call closeparameterblock(param_file)
130
131 call cvmix_init_ddiff(strat_param_max=cs%strat_param_max, &
132 kappa_ddiff_s=us%Z2_T_to_m2_s*cs%kappa_ddiff_s, &
133 ddiff_exp1=cs%ddiff_exp1, &
134 ddiff_exp2=cs%ddiff_exp2, &
135 mol_diff=us%Z2_T_to_m2_s*cs%mol_diff, &
136 kappa_ddiff_param1=cs%kappa_ddiff_param1, &
137 kappa_ddiff_param2=cs%kappa_ddiff_param2, &
138 kappa_ddiff_param3=cs%kappa_ddiff_param3, &
139 diff_conv_type=cs%diff_conv_type)
140
141end function cvmix_ddiff_init
142
143!> Subroutine for computing vertical diffusion coefficients for the
144!! double diffusion mixing parameterization.
145subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS, R_rho)
146
147 type(ocean_grid_type), intent(in) :: g !< Grid structure.
148 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
149 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
150 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
151 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
152 integer, intent(in) :: j !< Meridional grid index to work on.
153 ! Kd_T and Kd_S are intent inout because only one j-row is set here, but they are essentially outputs.
154 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: kd_t !< Interface double diffusion diapycnal
155 !! diffusivity for temperature
156 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
157 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: kd_s !< Interface double diffusion diapycnal
158 !! diffusivity for salinity
159 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
160 type(cvmix_ddiff_cs), pointer :: cs !< The control structure returned
161 !! by a previous call to CVMix_ddiff_init.
162 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
163 optional, intent(inout) :: r_rho !< The density ratios at interfaces [nondim].
164
165 ! Local variables
166 real, dimension(SZK_(GV)) :: &
167 cellheight, & !< Height of cell centers relative to the sea surface [H ~> m or kg m-2]
168 drho_dt, & !< partial derivatives of density with temperature [R C-1 ~> kg m-3 degC-1]
169 drho_ds, & !< partial derivatives of density with salinity [R S-1 ~> kg m-3 ppt-1]
170 pres_int, & !< pressure at each interface [R L2 T-2 ~> Pa]
171 temp_int, & !< temp and at interfaces [C ~> degC]
172 salt_int, & !< salt at at interfaces [S ~> ppt]
173 alpha_dt, & !< alpha*dT across interfaces [kg m-3]
174 beta_ds, & !< beta*dS across interfaces [kg m-3]
175 dt, & !< temperature difference between adjacent layers [C ~> degC]
176 ds !< salinity difference between adjacent layers [S ~> ppt]
177 real, dimension(SZK_(GV)+1) :: &
178 kd1_t, & !< Diapycanal diffusivity of temperature [m2 s-1].
179 kd1_s !< Diapycanal diffusivity of salinity [m2 s-1].
180
181 real, dimension(SZK_(GV)+1) :: ifaceheight !< Height of interfaces relative to the sea surface [H ~> m or kg m-2]
182 real :: dh, hcorr ! Limited thicknesses and a cumulative correction [H ~> m or kg m-2]
183 integer :: i, k
184
185 ! initialize dummy variables
186 pres_int(:) = 0.0 ; temp_int(:) = 0.0 ; salt_int(:) = 0.0
187 alpha_dt(:) = 0.0 ; beta_ds(:) = 0.0 ; drho_dt(:) = 0.0
188 drho_ds(:) = 0.0 ; dt(:) = 0.0 ; ds(:) = 0.0
189
190
191 ! GMM, I am leaving some code commented below. We need to pass BLD to
192 ! this subroutine to avoid adding diffusivity above that. This needs
193 ! to be done once we re-structure the order of the calls.
194 !if (.not. associated(hbl)) then
195 ! allocate(hbl(SZI_(G), SZJ_(G)))
196 ! hbl(:,:) = 0.0
197 !endif
198
199 do i = g%isc, g%iec
200
201 ! skip calling at land points
202 if (g%mask2dT(i,j) == 0.) cycle
203
204 pres_int(1) = 0. ; if (associated(tv%p_surf)) pres_int(1) = tv%p_surf(i,j)
205 ! we don't have SST and SSS, so let's use values at top-most layer
206 temp_int(1) = tv%T(i,j,1) ; salt_int(1) = tv%S(i,j,1)
207 do k=2,gv%ke
208 ! pressure at interface
209 pres_int(k) = pres_int(k-1) + (gv%g_Earth * gv%H_to_RZ) * h(i,j,k-1)
210 ! temp and salt at interface
211 ! for temp: (t1*h1 + t2*h2)/(h1+h2)
212 temp_int(k) = (tv%T(i,j,k-1)*h(i,j,k-1) + tv%T(i,j,k)*h(i,j,k)) / (h(i,j,k-1)+h(i,j,k))
213 salt_int(k) = (tv%S(i,j,k-1)*h(i,j,k-1) + tv%S(i,j,k)*h(i,j,k)) / (h(i,j,k-1)+h(i,j,k))
214 ! dT and dS
215 dt(k) = (tv%T(i,j,k-1)-tv%T(i,j,k))
216 ds(k) = (tv%S(i,j,k-1)-tv%S(i,j,k))
217 enddo ! k-loop finishes
218
219 call calculate_density_derivs(temp_int, salt_int, pres_int, drho_dt, drho_ds, tv%eqn_of_state)
220
221 ! The "-1.0" below is needed so that the following criteria is satisfied:
222 ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then "salt finger"
223 ! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then "diffusive convection"
224 do k=1,gv%ke
225 alpha_dt(k) = -1.0*us%R_to_kg_m3*drho_dt(k) * dt(k)
226 beta_ds(k) = us%R_to_kg_m3*drho_ds(k) * ds(k)
227 enddo
228
229 if (present(r_rho)) then
230 do k=1,gv%ke
231 ! Set R_rho using Adcroft's rule of reciprocals.
232 r_rho(i,j,k) = 0.0 ; if (abs(beta_ds(k)) > 0.0) r_rho(i,j,k) = alpha_dt(k) / beta_ds(k)
233 ! avoid NaN's again for safety, perhaps unnecessarily.
234 if (r_rho(i,j,k) /= r_rho(i,j,k)) r_rho(i,j,k) = 0.0
235 enddo
236 endif
237
238 ifaceheight(1) = 0.0 ! BBL is all relative to the surface
239 hcorr = 0.0
240 ! compute heights at cell center and interfaces
241 do k=1,gv%ke
242 dh = h(i,j,k) ! Nominal thickness to use for increment, in height units
243 dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
244 hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
245 dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
246 cellheight(k) = ifaceheight(k) - 0.5 * dh
247 ifaceheight(k+1) = ifaceheight(k) - dh
248 enddo
249
250 ! gets index of the level and interface above hbl in [H ~> m or kg m-2]
251 !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight, hbl(i,j))
252
253 kd1_t(:) = 0.0 ; kd1_s(:) = 0.0
254 call cvmix_coeffs_ddiff(tdiff_out=kd1_t(:), &
255 sdiff_out=kd1_s(:), &
256 strat_param_num=alpha_dt(:), &
257 strat_param_denom=beta_ds(:), &
258 nlev=gv%ke, &
259 max_nlev=gv%ke)
260 do k=1,gv%ke+1
261 kd_t(i,j,k) = gv%m2_s_to_HZ_T * kd1_t(k)
262 kd_s(i,j,k) = gv%m2_s_to_HZ_T * kd1_s(k)
263 enddo
264
265 ! Do not apply mixing due to convection within the boundary layer
266 !do k=1,kOBL
267 ! Kd_T(i,j,k) = 0.0
268 ! Kd_S(i,j,k) = 0.0
269 !enddo
270
271 enddo ! i-loop
272
273end subroutine compute_ddiff_coeffs
274
275!> Reads the parameter "USE_CVMIX_DDIFF" and returns state.
276!! This function allows other modules to know whether this parameterization will
277!! be used without needing to duplicate the log entry.
278logical function cvmix_ddiff_is_used(param_file)
279 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
280 call get_param(param_file, mdl, "USE_CVMIX_DDIFF", cvmix_ddiff_is_used, &
281 default=.false., do_not_log=.true.)
282
283end function cvmix_ddiff_is_used
284
285!> Clear pointers and deallocate memory
286! NOTE: Placeholder destructor
287subroutine cvmix_ddiff_end(CS)
288 type(cvmix_ddiff_cs), pointer :: cs !< Control structure for this module that
289 !! will be deallocated in this subroutine
290end subroutine cvmix_ddiff_end
291
292end module mom_cvmix_ddiff