MOM6
MOM_CVMix_ddiff.F90
Go to the documentation of this file.
1 !> Interface to CVMix double diffusion scheme.
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_diag_mediator, only : diag_ctrl, time_type, register_diag_field
9 use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
12 use mom_debugging, only : hchksum
13 use mom_grid, only : ocean_grid_type
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
20 
21 #include <MOM_memory.h>
22 
24 
25 !> Control structure including parameters for CVMix double diffusion.
26 type, public :: cvmix_ddiff_cs
27 
28  ! Parameters
29  real :: strat_param_max !< maximum value for the stratification parameter [nondim]
30  real :: kappa_ddiff_s !< leading coefficient in formula for salt-fingering regime
31  !! for salinity diffusion [m2 s-1]
32  real :: ddiff_exp1 !< interior exponent in salt-fingering regime formula [nondim]
33  real :: ddiff_exp2 !< exterior exponent in salt-fingering regime formula [nondim]
34  real :: mol_diff !< molecular diffusivity [m2 s-1]
35  real :: kappa_ddiff_param1 !< exterior coefficient in diffusive convection regime [nondim]
36  real :: kappa_ddiff_param2 !< middle coefficient in diffusive convection regime [nondim]
37  real :: kappa_ddiff_param3 !< interior coefficient in diffusive convection regime [nondim]
38  real :: min_thickness !< Minimum thickness allowed [m]
39  character(len=4) :: diff_conv_type !< type of diffusive convection to use. Options are Marmorino &
40  !! Caldwell 1976 ("MC76"; default) and Kelley 1988, 1990 ("K90")
41  logical :: debug !< If true, turn on debugging
42 
43  ! Daignostic handles and pointers
44  type(diag_ctrl), pointer :: diag => null() !< Pointer to diagnostics control structure
45  !>@{ Diagnostics handles
46  integer :: id_kt_extra = -1, id_ks_extra = -1, id_r_rho = -1
47  !!@}
48 
49  ! Diagnostics arrays
50 ! real, allocatable, dimension(:,:,:) :: KT_extra !< Double diffusion diffusivity for temp [Z2 s-1 ~> m2 s-1]
51 ! real, allocatable, dimension(:,:,:) :: KS_extra !< Double diffusion diffusivity for salt [Z2 s-1 ~> m2 s-1]
52  real, allocatable, dimension(:,:,:) :: r_rho !< Double-diffusion density ratio [nondim]
53 
54 end type cvmix_ddiff_cs
55 
56 character(len=40) :: mdl = "MOM_CVMix_ddiff" !< This module's name.
57 
58 contains
59 
60 !> Initialized the CVMix double diffusion module.
61 logical function cvmix_ddiff_init(Time, G, GV, US, param_file, diag, CS)
62 
63  type(time_type), intent(in) :: time !< The current time.
64  type(ocean_grid_type), intent(in) :: g !< Grid structure.
65  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
66  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
67  type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle
68  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure.
69  type(cvmix_ddiff_cs), pointer :: cs !< This module's control structure.
70 
71 ! This include declares and sets the variable "version".
72 #include "version_variable.h"
73 
74  if (associated(cs)) then
75  call mom_error(warning, "CVMix_ddiff_init called with an associated "// &
76  "control structure.")
77  return
78  endif
79  allocate(cs)
80 
81  ! Read parameters
82  call log_version(param_file, mdl, version, &
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.", &
88  default=.false.)
89 
90  if (.not. cvmix_ddiff_init) return
91 
92  call get_param(param_file, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)
93 
94  call get_param(param_file, mdl, 'MIN_THICKNESS', cs%min_thickness, default=0.001, do_not_log=.true.)
95 
96  call openparameterblock(param_file,'CVMIX_DDIFF')
97 
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)
101 
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)
105 
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)
109 
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)
113 
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)
117 
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)
121 
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)
125 
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)
129 
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).", &
133  default="MC76")
134 
135  call closeparameterblock(param_file)
136 
137  ! Register diagnostics
138  cs%diag => diag
139 
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)
142 
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)
145 
146  cs%id_R_rho = register_diag_field('ocean_model','R_rho',diag%axesTi,time, &
147  'Double-diffusion density ratio', 'nondim')
148 
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
152  endif
153 
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)
163 
164 end function cvmix_ddiff_init
165 
166 !> Subroutine for computing vertical diffusion coefficients for the
167 !! double diffusion mixing parameterization.
168 subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS)
169 
170  type(ocean_grid_type), intent(in) :: g !< Grid structure.
171  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
172  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
173  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
174  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
175  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd_t !< Interface double diffusion diapycnal
176  !! diffusivity for temp [Z2 T-1 ~> m2 s-1].
177  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd_s !< Interface double diffusion diapycnal
178  !! diffusivity for salt [Z2 T-1 ~> m2 s-1].
179  type(cvmix_ddiff_cs), pointer :: cs !< The control structure returned
180  !! by a previous call to CVMix_ddiff_init.
181  integer, intent(in) :: j !< Meridional grid indice.
182  ! Local variables
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]
193  ds !< salt difference between adjacent layers [ppt]
194  real, dimension(SZK_(G)+1) :: &
195  kd1_t, & !< Diapycanal diffusivity of temperature [m2 s-1].
196  kd1_s !< Diapycanal diffusivity of salinity [m2 s-1].
197 
198  real, dimension(SZK_(G)+1) :: ifaceheight !< Height of interfaces [m]
199  integer :: kobl !< level of OBL extent
200  real :: pref, g_o_rho0, rhok, rhokm1, dz, dh, hcorr
201  integer :: i, k
202 
203  ! initialize dummy variables
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
207 
208 
209  ! GMM, I am leaving some code commented below. We need to pass BLD to
210  ! this soubroutine to avoid adding diffusivity above that. This needs
211  ! to be done once we re-structure the order of the calls.
212  !if (.not. associated(hbl)) then
213  ! allocate(hbl(SZI_(G), SZJ_(G)));
214  ! hbl(:,:) = 0.0
215  !endif
216 
217  do i = g%isc, g%iec
218 
219  ! skip calling at land points
220  if (g%mask2dT(i,j) == 0.) cycle
221 
222  pref = 0.
223  pres_int(1) = pref
224  ! we don't have SST and SSS, so let's use values at top-most layer
225  temp_int(1) = tv%T(i,j,1); salt_int(1) = tv%S(i,j,1)
226  do k=2,g%ke
227  ! pressure at interface
228  pres_int(k) = pref + gv%H_to_Pa * h(i,j,k-1)
229  ! temp and salt at interface
230  ! for temp: (t1*h1 + t2*h2)/(h1+h2)
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))
233  ! dT and dS
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)
237  enddo ! k-loop finishes
238 
239  call calculate_density_derivs(temp_int(:), salt_int(:), pres_int(:), drho_dt(:), drho_ds(:), 1, &
240  g%ke, tv%EQN_OF_STATE)
241 
242  ! The "-1.0" below is needed so that the following criteria is satisfied:
243  ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then "salt finger"
244  ! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then "diffusive convection"
245  do k=1,g%ke
246  alpha_dt(k) = -1.0*drho_dt(k) * dt(k)
247  beta_ds(k) = drho_ds(k) * ds(k)
248  enddo
249 
250  if (cs%id_R_rho > 0.0) then
251  do k=1,g%ke
252  cs%R_rho(i,j,k) = alpha_dt(k)/beta_ds(k)
253  ! avoid NaN's
254  if(cs%R_rho(i,j,k) /= cs%R_rho(i,j,k)) cs%R_rho(i,j,k) = 0.0
255  enddo
256  endif
257 
258  ifaceheight(1) = 0.0 ! BBL is all relative to the surface
259  hcorr = 0.0
260  ! compute heights at cell center and interfaces
261  do k=1,g%ke
262  dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
263  dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
264  hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
265  dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
266  cellheight(k) = ifaceheight(k) - 0.5 * dh
267  ifaceheight(k+1) = ifaceheight(k) - dh
268  enddo
269 
270  ! gets index of the level and interface above hbl
271  !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j))
272 
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(:), &
278  nlev=g%ke, &
279  max_nlev=g%ke)
280  do k=1,g%ke+1
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)
283  enddo
284 
285  ! Do not apply mixing due to convection within the boundary layer
286  !do k=1,kOBL
287  ! Kd_T(i,j,k) = 0.0
288  ! Kd_S(i,j,k) = 0.0
289  !enddo
290 
291  enddo ! i-loop
292 
293 end subroutine compute_ddiff_coeffs
294 
295 !> Reads the parameter "USE_CVMIX_DDIFF" and returns state.
296 !! This function allows other modules to know whether this parameterization will
297 !! be used without needing to duplicate the log entry.
298 logical function cvmix_ddiff_is_used(param_file)
299  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
300  call get_param(param_file, mdl, "USE_CVMIX_DDIFF", cvmix_ddiff_is_used, &
301  default=.false., do_not_log = .true.)
302 
303 end function cvmix_ddiff_is_used
304 
305 !> Clear pointers and dealocate memory
306 subroutine cvmix_ddiff_end(CS)
307  type(cvmix_ddiff_cs), pointer :: cs !< Control structure for this module that
308  !! will be deallocated in this subroutine
309 
310  deallocate(cs)
311 
312 end subroutine cvmix_ddiff_end
313 
314 end module mom_cvmix_ddiff
mom_cvmix_ddiff::cvmix_ddiff_is_used
logical function, public cvmix_ddiff_is_used(param_file)
Reads the parameter "USE_CVMIX_DDIFF" and returns state. This function allows other modules to know w...
Definition: MOM_CVMix_ddiff.F90:299
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_cvmix_ddiff::cvmix_ddiff_cs
Control structure including parameters for CVMix double diffusion.
Definition: MOM_CVMix_ddiff.F90:26
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_cvmix_ddiff::cvmix_ddiff_init
logical function, public cvmix_ddiff_init(Time, G, GV, US, param_file, diag, CS)
Initialized the CVMix double diffusion module.
Definition: MOM_CVMix_ddiff.F90:62
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_cvmix_ddiff::compute_ddiff_coeffs
subroutine, public compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS)
Subroutine for computing vertical diffusion coefficients for the double diffusion mixing parameteriza...
Definition: MOM_CVMix_ddiff.F90:169
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
mom_cvmix_ddiff
Interface to CVMix double diffusion scheme.
Definition: MOM_CVMix_ddiff.F90:2
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.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_file_parser::openparameterblock
subroutine, public openparameterblock(CS, blockName, desc)
Tags blockName onto the end of the active parameter block name.
Definition: MOM_file_parser.F90:2015
mom_cvmix_ddiff::cvmix_ddiff_end
subroutine, public cvmix_ddiff_end(CS)
Clear pointers and dealocate memory.
Definition: MOM_CVMix_ddiff.F90:307
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_cvmix_ddiff::mdl
character(len=40) mdl
This module's name.
Definition: MOM_CVMix_ddiff.F90:56
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_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_file_parser::closeparameterblock
subroutine, public closeparameterblock(CS)
Remove the lowest level of recursion from the active block name.
Definition: MOM_file_parser.F90:2033
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