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>
23 public cvmix_ddiff_init, cvmix_ddiff_end, cvmix_ddiff_is_used, compute_ddiff_coeffs
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"
61 logical function cvmix_ddiff_init(Time, G, GV, US, param_file, diag, CS)
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")
84 call get_param(param_file, mdl,
"USE_CVMIX_DDIFF", cvmix_ddiff_init, &
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.", &
90 if (.not. cvmix_ddiff_init)
return
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')
148 if (cs%id_R_rho > 0) &
149 allocate(cs%R_rho( szi_(g), szj_(g), szk_(g)+1)); cs%R_rho(:,:,:) = 0.0
151 call cvmix_init_ddiff(strat_param_max=cs%strat_param_max, &
152 kappa_ddiff_s=cs%kappa_ddiff_s, &
153 ddiff_exp1=cs%ddiff_exp1, &
154 ddiff_exp2=cs%ddiff_exp2, &
155 mol_diff=cs%mol_diff, &
156 kappa_ddiff_param1=cs%kappa_ddiff_param1, &
157 kappa_ddiff_param2=cs%kappa_ddiff_param2, &
158 kappa_ddiff_param3=cs%kappa_ddiff_param3, &
159 diff_conv_type=cs%diff_conv_type)
161 end function cvmix_ddiff_init
165 subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS)
170 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
172 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(out) :: kd_t
174 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(out) :: kd_s
178 integer,
intent(in) :: j
180 real,
dimension(SZK_(G)) :: &
181 cellheight, & !< Height of cell centers [m]
182 drho_dt, & !< partial derivatives of density wrt temp [kg m-3 degC-1]
183 drho_ds, & !< partial derivatives of density wrt saln [kg m-3 ppt-1]
184 pres_int, & !< pressure at each interface [Pa]
185 temp_int, & !< temp and at interfaces [degC]
186 salt_int, & !< salt at at interfaces [ppt]
187 alpha_dt, & !< alpha*dT across interfaces
188 beta_ds, & !< beta*dS across interfaces
189 dt, & !< temp. difference between adjacent layers [degC]
191 real,
dimension(SZK_(G)+1) :: &
192 kd1_t, & !< Diapycanal diffusivity of temperature [m2 s-1].
195 real,
dimension(SZK_(G)+1) :: ifaceheight
197 real :: pref, g_o_rho0, rhok, rhokm1, dz, dh, hcorr
201 pres_int(:) = 0.0; temp_int(:) = 0.0; salt_int(:) = 0.0
202 alpha_dt(:) = 0.0; beta_ds(:) = 0.0; drho_dt(:) = 0.0
203 drho_ds(:) = 0.0; dt(:) = 0.0; ds(:) = 0.0
217 if (g%mask2dT(i,j) == 0.) cycle
222 temp_int(1) = tv%T(i,j,1); salt_int(1) = tv%S(i,j,1)
225 pres_int(k) = pref + gv%H_to_Pa * h(i,j,k-1)
228 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))
229 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))
231 dt(k) = (tv%T(i,j,k-1)-tv%T(i,j,k))
232 ds(k) = (tv%S(i,j,k-1)-tv%S(i,j,k))
233 pref = pref + gv%H_to_Pa * h(i,j,k-1)
237 g%ke, tv%EQN_OF_STATE)
243 alpha_dt(k) = -1.0*drho_dt(k) * dt(k)
244 beta_ds(k) = drho_ds(k) * ds(k)
247 if (cs%id_R_rho > 0.0)
then
249 cs%R_rho(i,j,k) = alpha_dt(k)/beta_ds(k)
251 if(cs%R_rho(i,j,k) /= cs%R_rho(i,j,k)) cs%R_rho(i,j,k) = 0.0
259 dh = h(i,j,k) * gv%H_to_m
261 hcorr = min( dh - cs%min_thickness, 0. )
262 dh = max( dh, cs%min_thickness )
263 cellheight(k) = ifaceheight(k) - 0.5 * dh
264 ifaceheight(k+1) = ifaceheight(k) - dh
270 kd1_t(:) = 0.0 ; kd1_s(:) = 0.0
271 call cvmix_coeffs_ddiff(tdiff_out=kd1_t(:), &
272 sdiff_out=kd1_s(:), &
273 strat_param_num=alpha_dt(:), &
274 strat_param_denom=beta_ds(:), &
278 kd_t(i,j,k) = us%m2_s_to_Z2_T * kd1_t(k)
279 kd_s(i,j,k) = us%m2_s_to_Z2_T * kd1_s(k)
290 end subroutine compute_ddiff_coeffs
295 logical function cvmix_ddiff_is_used(param_file)
297 call get_param(param_file, mdl,
"USE_CVMIX_DDIFF", cvmix_ddiff_is_used, &
298 default=.false., do_not_log = .true.)
300 end function cvmix_ddiff_is_used
303 subroutine cvmix_ddiff_end(CS)
309 end subroutine cvmix_ddiff_end