17 use cvmix_ddiff,
only : cvmix_init_ddiff, cvmix_coeffs_ddiff
18 use cvmix_kpp,
only : cvmix_kpp_compute_kobl_depth
19 implicit none ;
private
21 #include <MOM_memory.h>
29 real :: strat_param_max
35 real :: kappa_ddiff_param1
36 real :: kappa_ddiff_param2
37 real :: kappa_ddiff_param3
39 character(len=4) :: diff_conv_type
46 integer :: id_kt_extra = -1, id_ks_extra = -1, id_r_rho = -1
52 real,
allocatable,
dimension(:,:,:) :: r_rho
56 character(len=40) ::
mdl =
"MOM_CVMix_ddiff"
63 type(time_type),
intent(in) :: time
68 type(
diag_ctrl),
target,
intent(inout) :: diag
72 #include "version_variable.h"
74 if (
associated(cs))
then
75 call mom_error(warning,
"CVMix_ddiff_init called with an associated "// &
83 "Parameterization of mixing due to double diffusion processes via CVMix")
85 "If true, turns on double diffusive processes via CVMix. "//&
86 "Note that double diffusive processes on viscosity are ignored "//&
87 "in CVMix, see http://cvmix.github.io/ for justification.", &
92 call get_param(param_file,
mdl,
'DEBUG', cs%debug, default=.false., do_not_log=.true.)
94 call get_param(param_file,
mdl,
'MIN_THICKNESS', cs%min_thickness, default=0.001, do_not_log=.true.)
96 call openparameterblock(param_file,
'CVMIX_DDIFF')
98 call get_param(param_file,
mdl,
"STRAT_PARAM_MAX", cs%strat_param_max, &
99 "The maximum value for the double dissusion stratification parameter", &
100 units=
"nondim", default=2.55)
102 call get_param(param_file,
mdl,
"KAPPA_DDIFF_S", cs%kappa_ddiff_s, &
103 "Leading coefficient in formula for salt-fingering regime "//&
104 "for salinity diffusion.", units=
"m2 s-1", default=1.0e-4)
106 call get_param(param_file,
mdl,
"DDIFF_EXP1", cs%ddiff_exp1, &
107 "Interior exponent in salt-fingering regime formula.", &
108 units=
"nondim", default=1.0)
110 call get_param(param_file,
mdl,
"DDIFF_EXP2", cs%ddiff_exp2, &
111 "Exterior exponent in salt-fingering regime formula.", &
112 units=
"nondim", default=3.0)
114 call get_param(param_file,
mdl,
"KAPPA_DDIFF_PARAM1", cs%kappa_ddiff_param1, &
115 "Exterior coefficient in diffusive convection regime.", &
116 units=
"nondim", default=0.909)
118 call get_param(param_file,
mdl,
"KAPPA_DDIFF_PARAM2", cs%kappa_ddiff_param2, &
119 "Middle coefficient in diffusive convection regime.", &
120 units=
"nondim", default=4.6)
122 call get_param(param_file,
mdl,
"KAPPA_DDIFF_PARAM3", cs%kappa_ddiff_param3, &
123 "Interior coefficient in diffusive convection regime.", &
124 units=
"nondim", default=-0.54)
126 call get_param(param_file,
mdl,
"MOL_DIFF", cs%mol_diff, &
127 "Molecular diffusivity used in CVMix double diffusion.", &
128 units=
"m2 s-1", default=1.5e-6)
130 call get_param(param_file,
mdl,
"DIFF_CONV_TYPE", cs%diff_conv_type, &
131 "type of diffusive convection to use. Options are Marmorino \n" //&
132 "and Caldwell 1976 (MC76) and Kelley 1988, 1990 (K90).", &
135 call closeparameterblock(param_file)
140 cs%id_KT_extra = register_diag_field(
'ocean_model',
'KT_extra',diag%axesTi,time, &
141 'Double-diffusive diffusivity for temperature',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
143 cs%id_KS_extra = register_diag_field(
'ocean_model',
'KS_extra',diag%axesTi,time, &
144 'Double-diffusive diffusivity for salinity',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
146 cs%id_R_rho = register_diag_field(
'ocean_model',
'R_rho',diag%axesTi,time, &
147 'Double-diffusion density ratio',
'nondim')
149 if (cs%id_R_rho > 0)
then
150 allocate(cs%R_rho( szi_(g), szj_(g), szk_(g)+1))
151 cs%R_rho(:,:,:) = 0.0
154 call cvmix_init_ddiff(strat_param_max=cs%strat_param_max, &
155 kappa_ddiff_s=cs%kappa_ddiff_s, &
156 ddiff_exp1=cs%ddiff_exp1, &
157 ddiff_exp2=cs%ddiff_exp2, &
158 mol_diff=cs%mol_diff, &
159 kappa_ddiff_param1=cs%kappa_ddiff_param1, &
160 kappa_ddiff_param2=cs%kappa_ddiff_param2, &
161 kappa_ddiff_param3=cs%kappa_ddiff_param3, &
162 diff_conv_type=cs%diff_conv_type)
173 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
175 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(out) :: kd_t
177 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(out) :: kd_s
181 integer,
intent(in) :: j
183 real,
dimension(SZK_(G)) :: &
184 cellheight, & !< Height of cell centers [m]
185 drho_dt, & !< partial derivatives of density wrt temp [kg m-3 degC-1]
186 drho_ds, & !< partial derivatives of density wrt saln [kg m-3 ppt-1]
187 pres_int, & !< pressure at each interface [Pa]
188 temp_int, & !< temp and at interfaces [degC]
189 salt_int, & !< salt at at interfaces [ppt]
190 alpha_dt, & !< alpha*dT across interfaces
191 beta_ds, & !< beta*dS across interfaces
192 dt, & !< temp. difference between adjacent layers [degC]
194 real,
dimension(SZK_(G)+1) :: &
195 kd1_t, & !< Diapycanal diffusivity of temperature [m2 s-1].
198 real,
dimension(SZK_(G)+1) :: ifaceheight
200 real :: pref, g_o_rho0, rhok, rhokm1, dz, dh, hcorr
204 pres_int(:) = 0.0; temp_int(:) = 0.0; salt_int(:) = 0.0
205 alpha_dt(:) = 0.0; beta_ds(:) = 0.0; drho_dt(:) = 0.0
206 drho_ds(:) = 0.0; dt(:) = 0.0; ds(:) = 0.0
220 if (g%mask2dT(i,j) == 0.) cycle
225 temp_int(1) = tv%T(i,j,1); salt_int(1) = tv%S(i,j,1)
228 pres_int(k) = pref + gv%H_to_Pa * h(i,j,k-1)
231 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))
232 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))
234 dt(k) = (tv%T(i,j,k-1)-tv%T(i,j,k))
235 ds(k) = (tv%S(i,j,k-1)-tv%S(i,j,k))
236 pref = pref + gv%H_to_Pa * h(i,j,k-1)
240 g%ke, tv%EQN_OF_STATE)
246 alpha_dt(k) = -1.0*drho_dt(k) * dt(k)
247 beta_ds(k) = drho_ds(k) * ds(k)
250 if (cs%id_R_rho > 0.0)
then
252 cs%R_rho(i,j,k) = alpha_dt(k)/beta_ds(k)
254 if(cs%R_rho(i,j,k) /= cs%R_rho(i,j,k)) cs%R_rho(i,j,k) = 0.0
262 dh = h(i,j,k) * gv%H_to_m
264 hcorr = min( dh - cs%min_thickness, 0. )
265 dh = max( dh, cs%min_thickness )
266 cellheight(k) = ifaceheight(k) - 0.5 * dh
267 ifaceheight(k+1) = ifaceheight(k) - dh
273 kd1_t(:) = 0.0 ; kd1_s(:) = 0.0
274 call cvmix_coeffs_ddiff(tdiff_out=kd1_t(:), &
275 sdiff_out=kd1_s(:), &
276 strat_param_num=alpha_dt(:), &
277 strat_param_denom=beta_ds(:), &
281 kd_t(i,j,k) = us%m2_s_to_Z2_T * kd1_t(k)
282 kd_s(i,j,k) = us%m2_s_to_Z2_T * kd1_s(k)
301 default=.false., do_not_log = .true.)