src/parameterizations/vertical/MOM_diabatic_driver.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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		!> This routine drives the diabatic/dianeutral physics for MOM
6		module MOM_diabatic_driver
7
8		use MOM_bulk_mixed_layer, only : bulkmixedlayer, bulkmixedlayer_init, bulkmixedlayer_CS
9		use MOM_debugging, only : hchksum
10		use MOM_checksum_packages, only : MOM_state_chksum, MOM_state_stats
11		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
12		use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
13		use MOM_CVMix_shear, only : CVMix_shear_is_used
14		use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used
15		use MOM_diabatic_aux, only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_CS
16		use MOM_diabatic_aux, only : make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS
17		use MOM_diabatic_aux, only : triDiagTS_Eulerian, find_uv_at_h
18		use MOM_diabatic_aux, only : applyBoundaryFluxesInOut, set_pen_shortwave
19		use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
20		use MOM_diag_mediator, only : post_product_sum_u, post_product_sum_v
21		use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_remap_grids
22		use MOM_diag_mediator, only : diag_ctrl, query_averaging_enabled, enable_averages, disable_averaging
23		use MOM_diag_mediator, only : diag_grid_storage, diag_grid_storage_init, diag_grid_storage_end
24		use MOM_diag_mediator, only : diag_copy_diag_to_storage, diag_copy_storage_to_diag
25		use MOM_diag_mediator, only : diag_save_grids, diag_restore_grids
26		use MOM_diagnose_mld, only : diagnoseMLDbyDensityDifference, diagnoseMLDbyEnergy
27		use MOM_diagnose_kdwork, only : vbf_CS, KdWork_init, KdWork_end, KdWork_diagnostics
28		use MOM_diagnose_kdwork, only : Allocate_VBF_CS, Deallocate_VBF_CS
29		use MOM_diapyc_energy_req, only : diapyc_energy_req_init, diapyc_energy_req_end
30		use MOM_diapyc_energy_req, only : diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_CS
31		use MOM_CVMix_conv, only : CVMix_conv_init, CVMix_conv_cs
32		use MOM_CVMix_conv, only : calculate_CVMix_conv
33		use MOM_domains, only : pass_var, To_West, To_South, To_All, Omit_Corners
34		use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
35		use MOM_energetic_PBL, only : energetic_PBL, energetic_PBL_init
36		use MOM_energetic_PBL, only : energetic_PBL_end, energetic_PBL_CS
37		use MOM_energetic_PBL, only : energetic_PBL_get_MLD
38		use MOM_entrain_diffusive, only : entrainment_diffusive, entrain_diffusive_init
39		use MOM_entrain_diffusive, only : entrain_diffusive_CS
40		use MOM_EOS, only : calculate_density, calculate_density_derivs, calculate_TFreeze
41		use MOM_EOS, only : calculate_specific_vol_derivs, EOS_domain
42		use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery,MOM_mesg
43		use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
44		use MOM_file_parser, only : get_param, log_version, param_file_type, read_param
45		use MOM_forcing_type, only : forcing, MOM_forcing_chksum, find_ustar
46		use MOM_forcing_type, only : calculateBuoyancyFlux2d, forcing_SinglePointPrint
47		use MOM_geothermal, only : geothermal_entraining, geothermal_in_place
48		use MOM_geothermal, only : geothermal_init, geothermal_end, geothermal_CS
49		use MOM_grid, only : ocean_grid_type
50		use MOM_int_tide_input, only : set_int_tide_input, int_tide_input_init
51		use MOM_int_tide_input, only : int_tide_input_end, int_tide_input_CS, int_tide_input_type
52		use MOM_interface_heights, only : find_eta, calc_derived_thermo, thickness_to_dz
53		use MOM_interface_heights, only : convert_MLD_to_ML_thickness
54		use MOM_internal_tides, only : propagate_int_tide, register_int_tide_restarts
55		use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS
56		use MOM_kappa_shear, only : kappa_shear_is_used
57		use MOM_CVMix_KPP, only : KPP_CS, KPP_init, KPP_compute_BLD, KPP_calculate
58		use MOM_CVMix_KPP, only : KPP_end, KPP_get_BLD, register_KPP_restarts
59		use MOM_CVMix_KPP, only : KPP_NonLocalTransport_temp, KPP_NonLocalTransport_saln
60		use MOM_oda_incupd, only : apply_oda_incupd, oda_incupd_CS
61		use MOM_opacity, only : opacity_init, opacity_end, opacity_CS
62		use MOM_opacity, only : absorbRemainingSW, optics_type, optics_nbands
63		use MOM_open_boundary, only : ocean_OBC_type
64		use MOM_regularize_layers, only : regularize_layers, regularize_layers_init, regularize_layers_CS
65		use MOM_restart, only : MOM_restart_CS
66		use MOM_set_diffusivity, only : set_diffusivity, set_BBL_TKE
67		use MOM_set_diffusivity, only : set_diffusivity_init, set_diffusivity_end
68		use MOM_set_diffusivity, only : set_diffusivity_CS
69		use MOM_sponge, only : apply_sponge, sponge_CS
70		use MOM_ALE_sponge, only : apply_ALE_sponge, ALE_sponge_CS
71		use MOM_time_manager, only : time_type, real_to_time, operator(-), operator(<=)
72		use MOM_tracer_flow_control, only : call_tracer_column_fns, tracer_flow_control_CS
73		use MOM_tracer_diabatic, only : tracer_vertdiff, tracer_vertdiff_Eulerian
74		use MOM_unit_scaling, only : unit_scale_type
75		use MOM_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs
76		use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d
77		use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units
78		use MOM_wave_interface, only : wave_parameters_CS
79		use MOM_stochastics, only : stochastic_CS
80
81		implicit none ; private
82
83		#include <MOM_memory.h>
84
85		public diabatic
86		public diabatic_driver_init
87		public diabatic_driver_end
88		public extract_diabatic_member
89		public adiabatic
90		public adiabatic_driver_init
91		public register_diabatic_restarts
92
93		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
94		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
95		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
96		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
97
98		!> Control structure for this module
99		type, public :: diabatic_CS ; private
100		logical :: initialized = .false. !< True if this control structure has been initialized.
101
102		logical :: use_legacy_diabatic !< If true (default), use a legacy version of the diabatic
103		!! algorithm. This is temporary and is needed to avoid change
104		!! in answers.
105		logical :: use_bulkmixedlayer !< If true, a refined bulk mixed layer is used with
106		!! nkml sublayers (and additional buffer layers).
107		logical :: use_energetic_PBL !< If true, use the implicit energetics planetary
108		!! boundary layer scheme to determine the diffusivity
109		!! in the surface boundary layer.
110		logical :: use_KPP !< If true, use CVMix/KPP boundary layer scheme to determine the
111		!! OBLD and the diffusivities within this layer.
112		logical :: use_kappa_shear !< If true, use the kappa_shear module to find the
113		!! shear-driven diapycnal diffusivity.
114		logical :: use_CVMix_shear !< If true, use the CVMix module to find the
115		!! shear-driven diapycnal diffusivity.
116		logical :: use_CVMix_ddiff !< If true, use the CVMix double diffusion module.
117		logical :: use_CVMix_conv !< If true, use the CVMix module to get enhanced
118		!! mixing due to convection.
119		logical :: double_diffuse !< If true, some form of double-diffusive mixing is used.
120		logical :: use_sponge !< If true, sponges may be applied anywhere in the
121		!! domain. The exact location and properties of
122		!! those sponges are set by calls to
123		!! initialize_sponge and set_up_sponge_field.
124		logical :: use_oda_incupd !< If True, DA incremental update is
125		!! applied everywhere
126		logical :: use_geothermal !< If true, apply geothermal heating.
127		logical :: use_int_tides !< If true, use the code that advances a separate set
128		!! of equations for the internal tide energy density.
129		logical :: ePBL_is_additive !< If true, the diffusivity from ePBL is added to all
130		!! other diffusivities. Otherwise, the larger of kappa-
131		!! shear and ePBL diffusivities are used.
132		real :: ePBL_Prandtl !< The Prandtl number used by ePBL to convert vertical
133		!! diffusivities into viscosities [nondim].
134		logical :: useALEalgorithm !< If true, use the ALE algorithm rather than layered
135		!! isopycnal/stacked shallow water mode. This logical
136		!! passed by argument to diabatic_driver_init.
137		logical :: aggregate_FW_forcing !< Determines whether net incoming/outgoing surface
138		!! FW fluxes are applied separately or combined before
139		!! being applied.
140		real :: ML_mix_first !< The nondimensional fraction of the mixed layer
141		!! algorithm that is applied before diffusive mixing [nondim].
142		!! The default is 0, while 0.5 gives Strang splitting
143		!! and 1 is a sensible value too. Note that if there
144		!! are convective instabilities in the initial state,
145		!! the first call may do much more than the second.
146		integer :: NKBL !< The number of buffer layers (if bulk_mixed_layer)
147		logical :: massless_match_targets !< If true (the default), keep the T & S
148		!! consistent with the target values.
149		logical :: mix_boundary_tracers !< If true, mix the passive tracers in massless layers at the
150		!! bottom into the interior as though a diffusivity of
151		!! Kd_min_tr (see below) were operating.
152		logical :: mix_boundary_tracer_ALE !< If true, in ALE mode mix the passive tracers in massless
153		!! layers at the bottom into the interior as though a
154		!! diffusivity of Kd_min_tr (see below) were operating.
155		real :: Kd_BBL_tr !< A bottom boundary layer tracer diffusivity that
156		!! will allow for explicitly specified bottom fluxes
157		!! [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1]. The entrainment at the
158		!! bottom is at least sqrt(Kd_BBL_tr*dt) over the same distance.
159		real :: Kd_min_tr !< A minimal diffusivity that should always be
160		!! applied to tracers, especially in massless layers
161		!! near the bottom [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
162		real :: minimum_forcing_depth !< The smallest depth over which heat and freshwater
163		!! fluxes are applied [H ~> m or kg m-2].
164		real :: evap_CFL_limit = 0.8 !< The largest fraction of a layer that can be
165		!! evaporated in one time-step [nondim].
166		integer :: halo_TS_diff = 0 !< The temperature, salinity and thickness halo size that
167		!! must be valid for the diffusivity calculations.
168		integer :: halo_diabatic = 0 !< The temperature, salinity, specific volume and thickness
169		!! halo size that must be valid for the diabatic calculations,
170		!! including vertical mixing and internal tide propagation.
171		logical :: useKPP = .false. !< use CVMix/KPP diffusivities and non-local transport
172		logical :: KPPisPassive !< If true, KPP is in passive mode, not changing answers.
173		logical :: debug !< If true, write verbose checksums for debugging purposes.
174		logical :: debugConservation !< If true, monitor conservation and extrema.
175		logical :: tracer_tridiag !< If true, use tracer_vertdiff instead of tridiagTS for
176		!< vertical diffusion of T and S
177		logical :: debug_energy_req !< If true, test the mixing energy requirement code.
178		type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output
179		real :: MLDdensityDifference !< Density difference used to determine MLD_user [R ~> kg m-3]
180		real :: dz_subML_N2 !< The distance over which to calculate a diagnostic of the
181		!! average stratification at the base of the mixed layer [Z ~> m].
182		real :: MLD_En_vals(3) !< Energy values for energy mixed layer diagnostics [R Z3 T-2 ~> J m-2]
183		real :: BMLD_En_vals(3) !< Energy values for energy bottom mixed layer diagnostics [R Z3 T-2 ~> J m-2]
184		logical :: use_OM4_MLD_En_iter !< If true, uses an older iteration in the energetics MLD calculation to bitwise
185		!! reproduce OM4 era models
186		real :: ref_h_mld = 0.0 !< The depth of the "surface" density used in a difference mixed based
187		!! MLD calculation [Z ~> m].
188		logical :: Use_KdWork_diag = .false. !< Logical flag to indicate if any Kd_work diagnostics are on.
189		logical :: Use_N2_diag = .false. !< Logical flag to indicate if any N2 diagnostics are on.
190
191		!>@{ Diagnostic IDs
192		integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic
193		integer :: id_ea_t = -1, id_eb_t = -1, id_ea_s = -1, id_eb_s = -1
194		integer :: id_Kd_heat = -1, id_Kd_salt = -1, id_Kd_int = -1, id_Kd_ePBL = -1
195		integer :: id_Tdif = -1, id_Sdif = -1, id_Tadv = -1, id_Sadv = -1
196		integer :: id_N2_dd = -1, id_N2_salt_dd = -1, id_N2_temp_dd
197		! These are handles to diagnostics related to the mixed layer properties.
198		integer :: id_MLD_003 = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_mlotstsq = -1
199		integer :: id_MLD_003_zr = -1, id_MLD_003_rr = -1
200		integer :: id_MLD_EN1 = -1, id_MLD_EN2 = -1, id_MLD_EN3 = -1, id_subMLN2 = -1
201		integer :: id_BMLD_EN1 = -1, id_BMLD_EN2 = -1, id_BMLD_EN3 = -1
202
203		! These are handles to diagnostics that are only available in non-ALE layered mode.
204		integer :: id_wd = -1
205		integer :: id_dudt_dia = -1, id_dvdt_dia = -1
206		integer :: id_hf_dudt_dia_2d = -1, id_hf_dvdt_dia_2d = -1
207
208		! diagnostic for fields prior to applying diapycnal physics
209		integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1
210		integer :: id_T_predia = -1, id_S_predia = -1, id_e_predia = -1
211
212		integer :: id_diabatic_diff_temp_tend = -1
213		integer :: id_diabatic_diff_saln_tend = -1
214		integer :: id_diabatic_diff_heat_tend = -1
215		integer :: id_diabatic_diff_salt_tend = -1
216		integer :: id_diabatic_diff_heat_tend_2d = -1
217		integer :: id_diabatic_diff_salt_tend_2d = -1
218		integer :: id_diabatic_diff_h = -1
219
220		integer :: id_boundary_forcing_h = -1
221		integer :: id_boundary_forcing_h_tendency = -1
222		integer :: id_boundary_forcing_temp_tend = -1
223		integer :: id_boundary_forcing_saln_tend = -1
224		integer :: id_boundary_forcing_heat_tend = -1
225		integer :: id_boundary_forcing_salt_tend = -1
226		integer :: id_boundary_forcing_heat_tend_2d = -1
227		integer :: id_boundary_forcing_salt_tend_2d = -1
228
229		integer :: id_frazil_h = -1
230		integer :: id_frazil_temp_tend = -1
231		integer :: id_frazil_heat_tend = -1
232		integer :: id_frazil_heat_tend_2d = -1
233		!>@}
234
235		logical :: diabatic_diff_tendency_diag = .false. !< If true calculate diffusive tendency diagnostics
236		logical :: boundary_forcing_tendency_diag = .false. !< If true calculate frazil diagnostics
237		logical :: frazil_tendency_diag = .false. !< If true calculate frazil tendency diagnostics
238
239		type(diabatic_aux_CS), pointer :: diabatic_aux_CSp => NULL() !< Control structure for a child module
240		type(int_tide_input_CS), pointer :: int_tide_input_CSp => NULL() !< Control structure for a child module
241		type(int_tide_input_type), pointer :: int_tide_input => NULL() !< Control structure for a child module
242		type(set_diffusivity_CS), pointer :: set_diff_CSp => NULL() !< Control structure for a child module
243		type(sponge_CS), pointer :: sponge_CSp => NULL() !< Control structure for a child module
244		type(ALE_sponge_CS), pointer :: ALE_sponge_CSp => NULL() !< Control structure for a child module
245		type(tracer_flow_control_CS), pointer :: tracer_flow_CSp => NULL() !< Control structure for a child module
246		type(optics_type), pointer :: optics => NULL() !< Control structure for a child module
247		type(KPP_CS), pointer :: KPP_CSp => NULL() !< Control structure for a child module
248		type(diapyc_energy_req_CS), pointer :: diapyc_en_rec_CSp => NULL() !< Control structure for a child module
249		type(oda_incupd_CS), pointer :: oda_incupd_CSp => NULL() !< Control structure for a child module
250		type(int_tide_CS), pointer :: int_tide_CSp => NULL() !< Control structure for a child module
251		type(vbf_CS), pointer :: VBF => NULL() !< Control structure for a child module
252
253
254		type(bulkmixedlayer_CS) :: bulkmixedlayer !< Bulk mixed layer control structure
255		type(CVMix_conv_CS) :: CVMix_conv !< CVMix convection control structure
256		type(energetic_PBL_CS) :: ePBL !< Energetic PBL control structure
257		type(entrain_diffusive_CS) :: entrain_diffusive !< Diffusive entrainment control structure
258		type(geothermal_CS) :: geothermal !< Geothermal control structure
259		type(opacity_CS) :: opacity !< Opacity control structure
260		type(regularize_layers_CS) :: regularize_layers !< Regularize layer control structure
261
262		type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass
263		type(diag_grid_storage) :: diag_grids_prev !< Stores diagnostic grids at some previous point in the algorithm
264
265		type(time_type), pointer :: Time !< Pointer to model time (needed for sponges)
266		end type diabatic_CS
267
268		!>@{ clock ids
269		integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity
270		integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge
271		integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap
272		integer :: id_clock_kpp, id_clock_oda_incupd
273		!>@}
274
275		contains
276
277		!> This subroutine imposes the diapycnal mass fluxes and the
278		!! accompanying diapycnal advection of momentum and tracers.
279	12	subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
280		G, GV, US, CS, stoch_CS, OBC, Waves)
281		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
282		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
283		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
284		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
285		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2]
286		type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
287		!! unused have NULL ptrs
288		real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Active mixed layer depth [Z ~> m]
289		type(forcing), intent(inout) :: fluxes !< points to forcing fields
290		!! unused fields have NULL ptrs
291		type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities,
292		!! BBL properties and related fields
293		type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum
294		!! equations, to enable the later derived
295		!! diagnostics, like energy budgets
296		type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations
297		real, intent(in) :: dt !< time increment [T ~> s]
298		type(time_type), intent(in) :: Time_end !< Time at the end of the interval
299		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
300		type(diabatic_CS), pointer :: CS !< module control structure
301		type(stochastic_CS), pointer :: stoch_CS !< stochastic control structure
302		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
303		type(Wave_parameters_CS), pointer :: Waves !< Surface gravity waves
304
305		! local variables
306		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
307	24	eta ! Interface heights before diapycnal mixing [Z ~> m]
308	24	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: temp_diag ! Previous temperature for diagnostics [C ~> degC]
309	24	real, dimension(SZI_(G)) :: T_freeze, & ! The freezing potential temperature at the current salinity [C ~> degC].
310	24	ps ! Surface pressure [R L2 T-2 ~> Pa]
311		real, dimension(SZI_(G),SZK_(GV)) :: &
312	24	pressure ! The pressure at the middle of each layer [R L2 T-2 ~> Pa].
313		real :: H_to_RL2_T2 ! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
314		integer :: i, j, k, is, ie, js, je, nz
315		logical :: showCallTree ! If true, show the call tree
316
317	12	real, allocatable, dimension(:,:,:) :: h_in ! thickness before thermodynamics [H ~> m or kg m-2]
318	12	real, allocatable, dimension(:,:,:) :: t_in ! temperature before thermodynamics [C ~> degC]
319	12	real, allocatable, dimension(:,:,:) :: s_in ! salinity before thermodynamics [S ~> ppt]
320
321	12	if (GV%ke == 1) return
322
323	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
324
325	12	if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
326	0	"Module must be initialized before it is used.")
327
328	12	if (.not. CS%initialized) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
329	0	"Module must be initialized before it is used.")
330
331	12	if (dt == 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
332	0	"diabatic was called with a zero length timestep.")
333	12	if (dt < 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
334	0	"diabatic was called with a negative timestep.")
335
336	12	showCallTree = callTree_showQuery()
337
338		! Offer diagnostics of various state variables at the start of diabatic
339		! these are mostly for debugging purposes.
340	12	if (CS%id_u_predia > 0) call post_data(CS%id_u_predia, u, CS%diag)
341	12	if (CS%id_v_predia > 0) call post_data(CS%id_v_predia, v, CS%diag)
342	12	if (CS%id_h_predia > 0) call post_data(CS%id_h_predia, h, CS%diag)
343	12	if (CS%id_T_predia > 0) call post_data(CS%id_T_predia, tv%T, CS%diag)
344	12	if (CS%id_S_predia > 0) call post_data(CS%id_S_predia, tv%S, CS%diag)
345	12	if (CS%id_e_predia > 0) then
346		!$omp target update to(h)
347		!$omp target enter data map(alloc: eta)
348	0	call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref)
349		!$omp target exit data map(from: eta)
350	0	call post_data(CS%id_e_predia, eta, CS%diag)
351		endif
352
353		! Save a copy of the initial state if stochastic perturbations are active.
354	12	if (stoch_CS%do_sppt) then
355	0	allocate(h_in(G%isd:G%ied, G%jsd:G%jed, GV%ke)) ; h_in(:,:,:) = h(:,:,:)
356	0	allocate(t_in(G%isd:G%ied, G%jsd:G%jed, GV%ke)) ; t_in(:,:,:) = tv%T(:,:,:)
357	0	allocate(s_in(G%isd:G%ied, G%jsd:G%jed, GV%ke)) ; s_in(:,:,:) = tv%S(:,:,:)
358		endif
359
360	12	if (CS%debug) then
361	0	call MOM_state_chksum("Start of diabatic ", u, v, h, G, GV, US, haloshift=0)
362	0	call MOM_forcing_chksum("Start of diabatic", fluxes, G, US, haloshift=0)
363		endif
364	12	if (CS%debugConservation) call MOM_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, G, GV, US)
365
366	12	if (CS%debug_energy_req) &
367	0	call diapyc_energy_req_test(h, dt, tv, G, GV, US, CS%diapyc_en_rec_CSp)
368
369	12	call cpu_clock_begin(id_clock_set_diffusivity)
370	12	call set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS%set_diff_CSp, OBC=OBC)
371	12	call cpu_clock_end(id_clock_set_diffusivity)
372
373		! Frazil formation keeps the temperature above the freezing point.
374		! make_frazil is deliberately called at both the beginning and at
375		! the end of the diabatic processes.
376	12	if (associated(tv%T) .AND. associated(tv%frazil)) then
377		! For frazil diagnostic, the first call covers the first half of the time step
378	12	call enable_averages(0.5dt, Time_end - real_to_time(0.5dt, unscale=US%T_to_s), CS%diag)
379	12	if (CS%frazil_tendency_diag) then
380	0	do k=1,nz ; do j=js,je ; do i=is,ie
381	0	temp_diag(i,j,k) = tv%T(i,j,k)
382		enddo ; enddo ; enddo
383		endif
384
385	12	if (associated(fluxes%p_surf_full)) then
386	0	call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff)
387		else
388	12	call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, halo=CS%halo_TS_diff)
389		endif
390	12	if (showCallTree) call callTree_waypoint("done with 1st make_frazil (diabatic)")
391
392	12	if (CS%frazil_tendency_diag) then
393	0	call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, G, GV, US, CS)
394	0	if (CS%id_frazil_h > 0) call post_data(CS%id_frazil_h, h, CS%diag)
395		endif
396	12	call disable_averaging(CS%diag)
397		endif ! associated(tv%T) .AND. associated(tv%frazil)
398	12	if (CS%debugConservation) call MOM_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, G, GV, US)
399
400	12	if (CS%use_int_tides) then
401		! This block provides an interface for the unresolved low-mode internal tide module.
402		call set_int_tide_input(u, v, h, tv, fluxes, CS%int_tide_input, dt, G, GV, US, &
403	0	CS%int_tide_input_CSp)
404
405		call propagate_int_tide(h, tv, CS%int_tide_input%Nb, CS%int_tide_input%Rho_bot, dt, &
406	0	G, GV, US, CS%int_tide_input_CSp, CS%int_tide_CSp)
407
408	0	if (showCallTree) call callTree_waypoint("done with propagate_int_tide (diabatic)")
409		endif ! end CS%use_int_tides
410
411	12	if (CS%useALEalgorithm .and. CS%use_legacy_diabatic) then
412		call diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
413	12	G, GV, US, CS, stoch_CS, Waves)
414	0	elseif (CS%useALEalgorithm) then
415		call diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
416	0	G, GV, US, CS, stoch_CS, Waves)
417		else
418		call layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
419	0	G, GV, US, CS, Waves)
420		endif
421
422	12	call cpu_clock_begin(id_clock_pass)
423	12	if (associated(visc%sfc_buoy_flx)) &
424	12	call pass_var(visc%sfc_buoy_flx, G%domain, halo=1, complete=.not.associated(visc%MLD))
425	12	if (associated(visc%h_ML)) &
426	12	call pass_var(visc%h_ML, G%Domain, halo=1, complete=.not.associated(visc%MLD))
427	12	if (associated(visc%MLD)) &
428	12	call pass_var(visc%MLD, G%Domain, halo=1, complete=.true.)
429	12	if (associated(visc%Kv_shear)) &
430	12	call pass_var(visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1)
431	12	call cpu_clock_end(id_clock_pass)
432
433	12	call disable_averaging(CS%diag)
434		! Frazil formation keeps temperature above the freezing point.
435		! make_frazil is deliberately called at both the beginning and at
436		! the end of the diabatic processes.
437	12	if (associated(tv%T) .AND. associated(tv%frazil)) then
438	12	call enable_averages(0.5*dt, Time_end, CS%diag)
439	12	if (CS%frazil_tendency_diag) then
440	0	do k=1,nz ; do j=js,je ; do i=is,ie
441	0	temp_diag(i,j,k) = tv%T(i,j,k)
442		enddo ; enddo ; enddo
443		endif
444
445	12	if (associated(fluxes%p_surf_full)) then
446	0	call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full)
447		else
448	12	call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp)
449		endif
450
451	12	if (CS%frazil_tendency_diag) then
452	0	call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, G, GV, US, CS)
453	0	if (CS%id_frazil_h > 0 ) call post_data(CS%id_frazil_h, h, CS%diag)
454		endif
455
456	12	if (showCallTree) call callTree_waypoint("done with 2nd make_frazil (diabatic)")
457	12	if (CS%debugConservation) call MOM_state_stats('2nd make_frazil', u, v, h, tv%T, tv%S, G, GV, US)
458	12	call disable_averaging(CS%diag)
459
460		endif ! endif for frazil
461
462	12	if (stoch_CS%do_sppt) then
463		! perturb diabatic tendencies.
464		! These stochastic perturbations do not conserve heat, salt or mass.
465	0	do k=1,nz ; do j=js,je ; do i=is,ie
466	0	h(i,j,k) = max(h_in(i,j,k) + (h(i,j,k)-h_in(i,j,k)) * stoch_CS%sppt_wts(i,j), GV%Angstrom_H)
467	0	tv%S(i,j,k) = max(s_in(i,j,k) + (tv%S(i,j,k)-s_in(i,j,k)) * stoch_CS%sppt_wts(i,j), 0.0)
468		enddo ; enddo ; enddo
469		! now that we have updated thickness and salinity, calculate freeing point
470	0	H_to_RL2_T2 = GV%H_to_RZ * GV%g_Earth
471	0	do j=js,je
472	0	ps(:) = 0.0
473	0	if (associated(fluxes%p_surf)) then
474	0	do i=is,ie
475	0	ps(i) = fluxes%p_surf(i,j)
476		enddo
477		endif
478
479	0	do i=is,ie
480	0	pressure(i,1) = ps(i) + (0.5H_to_RL2_T2)h(i,j,1)
481		enddo
482	0	do k=2,nz ; do i=is,ie
483		pressure(i,k) = pressure(i,k-1) + &
484	0	(0.5H_to_RL2_T2) (h(i,j,k) + h(i,j,k-1))
485		enddo ; enddo
486	0	do k=1,nz
487		call calculate_TFreeze(tv%S(is:ie,j,k), pressure(is:ie,k), T_freeze(is:ie), &
488	0	tv%eqn_of_state)
489	0	do i=is,ie
490	0	tv%T(i,j,k) = max(t_in(i,j,k) + (tv%T(i,j,k)-t_in(i,j,k)) * stoch_CS%sppt_wts(i,j), T_freeze(i))
491		enddo
492		enddo
493		enddo
494
495	0	deallocate(h_in, t_in, s_in)
496		endif
497
498		! Diagnose mixed layer depths.
499	12	call enable_averages(dt, Time_end, CS%diag)
500	12	if (CS%id_MLD_003 > 0 .or. CS%id_subMLN2 > 0 .or. CS%id_mlotstsq > 0) then
501		call diagnoseMLDbyDensityDifference(CS%id_MLD_003, h, tv, 0.03*US%kg_m3_to_R, G, GV, US, CS%diag, &
502		CS%ref_h_mld, CS%id_MLD_003_zr, CS%id_MLD_003_rr, &
503	0	id_N2subML=CS%id_subMLN2, id_MLDsq=CS%id_mlotstsq, dz_subML=CS%dz_subML_N2)
504		endif
505	12	if (CS%id_MLD_0125 > 0) then
506		call diagnoseMLDbyDensityDifference(CS%id_MLD_0125, h, tv, 0.125*US%kg_m3_to_R, G, GV, US, CS%diag, &
507	0	ref_H_MLD=0.0, id_ref_z=-1, id_ref_rho=-1)
508		endif
509	12	if (CS%id_MLD_user > 0) then
510		call diagnoseMLDbyDensityDifference(CS%id_MLD_user, h, tv, CS%MLDdensityDifference, G, GV, US, CS%diag, &
511	0	ref_H_MLD=0.0, id_ref_z=-1, id_ref_rho=-1)
512		endif
513	12	if ((CS%id_MLD_EN1 > 0) .or. (CS%id_MLD_EN2 > 0) .or. (CS%id_MLD_EN3 > 0)) then
514		! Surface Mixed Layer diagnostic
515		call diagnoseMLDbyEnergy((/CS%id_MLD_EN1, CS%id_MLD_EN2, CS%id_MLD_EN3/), h, tv, G, GV, US, CS%MLD_En_vals, &
516	0	(/1,nz/), CS%diag, OM4_iteration=CS%use_OM4_MLD_En_iter)
517		endif
518	12	if ((CS%id_BMLD_EN1 > 0) .or. (CS%id_BMLD_EN2 > 0) .or. (CS%id_BMLD_EN3 > 0)) then
519		! Bottom Mixed Layer diagnostic
520		call diagnoseMLDbyEnergy((/CS%id_BMLD_EN1, CS%id_BMLD_EN2, CS%id_BMLD_EN3/), h, tv, G, GV, US, CS%BMLD_En_vals, &
521	0	(/nz,1/), CS%diag, OM4_iteration=.false.)
522		endif
523	12	if (stoch_CS%do_sppt .and. stoch_CS%id_sppt_wts > 0) &
524	0	call post_data(stoch_CS%id_sppt_wts, stoch_CS%sppt_wts, CS%diag)
525
526	12	call disable_averaging(CS%diag)
527
528	12	if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G, GV, US)
529
530	12	end subroutine diabatic
531
532
533		!> Applies diabatic forcing and diapycnal mixing of temperature, salinity and other tracers for use
534		!! with an ALE algorithm. This version uses an older set of algorithms compared with diabatic_ALE.
535	12	subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
536		G, GV, US, CS, stoch_CS, Waves)
537		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
538		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
539		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
540		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
541		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
542		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2]
543		type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
544		!! unused have NULL ptrs
545		real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Active mixed layer depth [Z ~> m]
546		type(forcing), intent(inout) :: fluxes !< points to forcing fields
547		!! unused fields have NULL ptrs
548		type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities,
549		!! BBL properties and related fields
550		type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum
551		!! equations, to enable the later derived
552		!! diagnostics, like energy budgets
553		type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations
554		real, intent(in) :: dt !< time increment [T ~> s]
555		type(time_type), intent(in) :: Time_end !< Time at the end of the interval
556		type(diabatic_CS), pointer :: CS !< module control structure
557		type(stochastic_CS), pointer :: stoch_CS !< stochastic control structure
558		type(Wave_parameters_CS), pointer :: Waves !< Surface gravity waves
559
560		! local variables
561		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
562	24	h_orig, & ! Initial layer thicknesses [H ~> m or kg m-2]
563	24	dz, & ! The vertical distance between interfaces around a layer [Z ~> m]
564	24	dSV_dT, & ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
565	24	dSV_dS, & ! The partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
566	24	cTKE, & ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
567	24	u_h, & ! Zonal velocities interpolated to thickness points [L T-1 ~> m s-1]
568	24	v_h, & ! Meridional velocities interpolated to thickness points [L T-1 ~> m s-1]
569	24	temp_diag, & ! Diagnostic array of previous temperatures [C ~> degC]
570	24	saln_diag ! Diagnostic array of previous salinity [S ~> ppt]
571
572		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
573	24	ent_s, & ! The diffusive coupling across interfaces within one time step for
574		! salinity and passive tracers [H ~> m or kg m-2]
575	24	ent_t, & ! The diffusive coupling across interfaces within one time step for
576		! temperature [H ~> m or kg m-2]
577	24	Kd_int, & ! diapycnal diffusivity of interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
578	24	Kd_heat, & ! diapycnal diffusivity of heat [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
579	24	Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
580	24	Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to
581		! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
582	24	Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to
583		! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
584	24	Kd_ePBL, & ! test array of diapycnal diffusivities at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
585	24	KPP_NLTheat, & ! KPP non-local transport for heat [nondim]
586	24	KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim]
587	24	KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]
588	24	Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
589	24	Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
590	24	N2_salt, & !< Salinity contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]
591	24	N2_temp !< Temperature contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]
592
593		real, dimension(SZI_(G),SZJ_(G)) :: &
594	24	U_star, & ! The friction velocity [Z T-1 ~> m s-1].
595	24	KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
596	24	KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
597	24	SkinBuoyFlux, & ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
598	24	BBL_BuoyFlux ! 2d bottom buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
599
600		real, dimension(SZI_(G)) :: &
601	24	p_i ,& ! Pressure at the interface [R L2 T-2 ~> Pa]
602	24	T_i, & ! Temperature at the interface [C ~> degC]
603	24	S_i, & ! Salinity at the interface [S ~> ppt]
604	24	drhodS, & ! Local change in density w.r.t. salinity using model EOS & state [R C-1 ~> kg m-3 ppt-1]
605	24	drhodT, & ! Local change in density w.r.t. temperature using model EOS & state [R C-1 ~> kg m-3 degC-1]
606	24	dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
607	24	dSpV_dS ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
608
609		logical, dimension(SZI_(G)) :: &
610	24	in_boundary ! True if there are no massive layers below, where massive is defined as
611		! sufficiently thick that the no-flux boundary conditions have not restricted
612		! the entrainment - usually sqrt(Kd*dt).
613
614		real :: h_neglect ! A thickness that is so small it is usually lost
615		! in roundoff and can be neglected [H ~> m or kg m-2]
616		real :: dz_neglect ! A vertical distance that is so small it is usually lost
617		! in roundoff and can be neglected [Z ~> m]
618		real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2]
619		real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2]
620		real :: I_dzval ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1]
621		real :: I_h ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1]
622		real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is
623		! coupled to the bottom within a timestep [H ~> m or kg m-2]
624		real :: Kd_add_here ! An added diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
625	24	real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].
626
627		real :: Ent_int ! The diffusive entrainment rate at an interface [H ~> m or kg m-2]
628		real :: Idt ! The inverse time step [T-1 ~> s-1]
629		real :: g_Rho0 ! G_Earth/Rho0 [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2]
630		real :: H_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
631		real :: alt_H_to_pres! A conversion factor from thicknesses to pressure w/ alternative scaling [R Z T-2 ~> Pa m-1]
632		logical :: nonBous ! True if not using the Boussinesq approximation
633
634		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
635
636		logical :: showCallTree ! If true, show the call tree
637		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
638
639	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
640	12	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
641	12	dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect
642	12	h_neglect = GV%H_subroundoff
643
644	12	nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
645	12	g_Rho0 = GV%g_Earth_Z_T2 / GV%H_to_RZ
646	12	H_to_pres = GV%H_to_RZ * GV%g_Earth
647	12	alt_H_to_pres = H_to_pres * US%L_to_Z*2 GV%Z_to_H
648
649	16001964	Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0
650
651	12	showCallTree = callTree_showQuery()
652	12	if (showCallTree) call callTree_enter("diabatic_ALE_legacy(), MOM_diabatic_driver.F90")
653
654		! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
655	12	call enable_averages(dt, Time_end, CS%diag)
656
657	12	if (CS%Use_KdWork_diag) call Allocate_VBF_CS(G, GV, CS%VBF)
658
659	12	if (CS%use_geothermal) then
660	12	call cpu_clock_begin(id_clock_geothermal)
661	12	call geothermal_in_place(h, tv, dt, G, GV, US, CS%geothermal, BBL_BuoyFlux, halo=CS%halo_TS_diff)
662	12	call cpu_clock_end(id_clock_geothermal)
663	12	if (showCallTree) call callTree_waypoint("geothermal (diabatic)")
664	12	if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G, GV, US)
665		endif
666
667		! Whenever thickness changes let the diag manager know, target grids
668		! for vertical remapping may need to be regenerated.
669	12	call diag_update_remap_grids(CS%diag)
670
671		! Set_pen_shortwave estimates the optical properties of the water column.
672		! It will need to be modified later to include information about the
673		! biological properties and layer thicknesses.
674	12	if (associated(CS%optics)) &
675	12	call set_pen_shortwave(CS%optics, fluxes, G, GV, US, CS%diabatic_aux_CSp, CS%opacity, CS%tracer_flow_CSp)
676
677	12	if (CS%debug) call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, US, haloshift=0)
678
679	12	if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then
680	12	if (CS%use_geothermal) then
681	12	call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, zero_mix=.true.)
682		else
683	0	call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
684		endif
685	12	if (showCallTree) call callTree_waypoint("done with find_uv_at_h (diabatic)")
686		endif
687
688	12	call cpu_clock_begin(id_clock_set_diffusivity)
689		! Sets: Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb
690		! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow
691	12	if (CS%debug) &
692	0	call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff)
693	12	if (CS%double_diffuse) then
694		call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, &
695	0	CS%set_diff_CSp, CS%VBF, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S)
696		else
697		call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, &
698	12	CS%set_diff_CSp, CS%VBF)
699		endif
700	12	call cpu_clock_end(id_clock_set_diffusivity)
701	12	if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)")
702
703
704	12	if (CS%debug) then
705	0	call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, US, haloshift=0)
706	0	call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, US, haloshift=0)
707	0	call MOM_thermovar_chksum("after set_diffusivity ", tv, G, US)
708	0	call hchksum(Kd_Int, "after set_diffusivity Kd_Int", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
709		endif
710
711		! Set diffusivities for heat and salt separately
712	12	if (CS%useKPP) then
713		! Add contribution from double diffusion
714	0	if (CS%double_diffuse) then
715		!$OMP parallel do default(shared)
716	0	do K=1,nz+1 ; do j=js,je ; do i=is,ie
717	0	Kd_salt(i,j,K) = Kd_int(i,j,K) + Kd_extra_S(i,j,K)
718	0	Kd_heat(i,j,K) = Kd_int(i,j,K) + Kd_extra_T(i,j,K)
719		enddo ; enddo ; enddo
720		else
721		!$OMP parallel do default(shared)
722	0	do K=1,nz+1 ; do j=js,je ; do i=is,ie
723	0	Kd_salt(i,j,K) = Kd_int(i,j,K)
724	0	Kd_heat(i,j,K) = Kd_int(i,j,K)
725		enddo ; enddo ; enddo
726		endif
727
728	0	if (CS%debug) then
729	0	call hchksum(Kd_heat, "after set_diffusivity Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
730	0	call hchksum(Kd_salt, "after set_diffusivity Kd_salt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
731		endif
732
733	0	call cpu_clock_begin(id_clock_kpp)
734		! total vertical viscosity in the interior is represented via visc%Kv_shear
735
736		! NOTE: The following do not require initialization, but their checksums do
737		! require initialization, and past versions were initialized to zero.
738	0	KPP_NLTheat(:,:,:) = 0.
739	0	KPP_NLTscalar(:,:,:) = 0.
740	0	KPP_buoy_flux(:,:,:) = 0.
741	0	KPP_temp_flux(:,:) = 0.
742	0	KPP_salt_flux(:,:) = 0.
743
744		! KPP needs the surface buoyancy flux but does not update state variables.
745		! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
746		! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux
747		! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second)
748		! unlike other instances where the fluxes are integrated in time over a time-step.
749		call calculateBuoyancyFlux2d(G, GV, US, fluxes, CS%optics, h, tv%T, tv%S, tv, &
750	0	KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux)
751
752		! Determine the friction velocity, perhaps using the evovling surface density.
753	0	call find_ustar(fluxes, tv, U_star, G, GV, US)
754
755		! The KPP scheme calculates boundary layer diffusivities and non-local transport.
756	0	if ( associated(fluxes%lamult) ) then
757		call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, &
758	0	U_star, KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult)
759
760		call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, &
761	0	Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult)
762		else
763		call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, &
764	0	U_star, KPP_buoy_flux, Waves=Waves)
765
766		call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, &
767	0	Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves)
768		endif
769
770	0	call KPP_get_BLD(CS%KPP_CSp, BLD(:,:), G, US)
771		! If visc%MLD or visc%h_ML exist, copy KPP's BLD into them with appropriate conversions.
772	0	if (associated(visc%h_ML)) call convert_MLD_to_ML_thickness(BLD, h, visc%h_ML, tv, G, GV)
773	0	if (associated(visc%MLD)) visc%MLD(:,:) = BLD(:,:)
774	0	if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = KPP_buoy_flux(:,:,1)
775
776	0	if (.not.CS%KPPisPassive) then
777		!$OMP parallel do default(shared)
778	0	do K=1,nz+1 ; do j=js,je ; do i=is,ie
779	0	Kd_int(i,j,K) = min( Kd_salt(i,j,K), Kd_heat(i,j,K) )
780		enddo ; enddo ; enddo
781	0	if (CS%double_diffuse) then
782		!$OMP parallel do default(shared)
783	0	do K=1,nz+1 ; do j=js,je ; do i=is,ie
784	0	Kd_extra_S(i,j,K) = (Kd_salt(i,j,K) - Kd_int(i,j,K))
785	0	Kd_extra_T(i,j,K) = (Kd_heat(i,j,K) - Kd_int(i,j,K))
786		enddo ; enddo ; enddo
787		endif
788		endif ! not passive
789
790	0	if (showCallTree) call callTree_waypoint("done with KPP_calculate (diabatic)")
791	0	if (CS%debug) then
792	0	call MOM_state_chksum("after KPP", u, v, h, G, GV, US, haloshift=0)
793	0	call MOM_forcing_chksum("after KPP", fluxes, G, US, haloshift=0)
794	0	call MOM_thermovar_chksum("after KPP", tv, G, US)
795	0	call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
796	0	call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
797		call hchksum(KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, &
798	0	unscale=US%C_to_degCGV%H_to_mUS%s_to_T)
799		call hchksum(KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, &
800	0	unscale=US%S_to_pptGV%H_to_mUS%s_to_T)
801	0	call hchksum(KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0)
802	0	call hchksum(KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0)
803		endif
804		! Apply non-local transport of heat and salt
805		! Changes: tv%T, tv%S
806		call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, KPP_NLTheat, KPP_temp_flux, &
807	0	dt, tv%tr_T, tv%T, tv%C_p)
808		call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, KPP_NLTscalar, KPP_salt_flux, &
809	0	dt, tv%tr_S, tv%S)
810	0	call cpu_clock_end(id_clock_kpp)
811	0	if (showCallTree) call callTree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
812	0	if (CS%debugConservation) call MOM_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, G, GV, US)
813
814	0	if (CS%debug) then
815	0	call MOM_state_chksum("after KPP_applyNLT ", u, v, h, G, GV, US, haloshift=0)
816	0	call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, US, haloshift=0)
817	0	call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G, US)
818		endif
819		endif ! endif for KPP
820
821		! This is the "old" method for applying differential diffusion.
822		! Changes: tv%T, tv%S
823	12	if (CS%double_diffuse .and. associated(tv%T)) then
824
825	0	call cpu_clock_begin(id_clock_differential_diff)
826	0	call differential_diffuse_T_S(h, tv%T, tv%S, Kd_extra_T, Kd_extra_S, tv, dt, G, GV)
827	0	call cpu_clock_end(id_clock_differential_diff)
828
829	0	if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)")
830	0	if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G, GV, US)
831
832		! increment heat and salt diffusivity.
833		! CS%useKPP==.true. already has extra_T and extra_S included
834	0	if (.not. CS%useKPP) then
835		!$OMP parallel do default(shared)
836	0	do K=2,nz ; do j=js,je ; do i=is,ie
837	0	Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_extra_T(i,j,K)
838	0	Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_extra_S(i,j,K)
839		enddo ; enddo ; enddo
840		endif
841
842		endif
843
844		! Calculate vertical mixing due to convection (computed via CVMix)
845	12	if (CS%use_CVMix_conv) then
846		! Increment vertical diffusion and viscosity due to convection
847	0	call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv, BLD, Kd_int, visc%Kv_shear)
848		endif
849
850		! Find the vertical distances across layers.
851	12	call thickness_to_dz(h, tv, dz, G, GV, US)
852
853		! This block sets ent_t and ent_s from h and Kd_int.
854	87132	do j=js,je ; do i=is,ie
855	86400	ent_s(i,j,1) = 0.0 ; ent_s(i,j,nz+1) = 0.0
856	87120	ent_t(i,j,1) = 0.0 ; ent_t(i,j,nz+1) = 0.0
857		enddo ; enddo
858		!$OMP parallel do default(shared) private(I_dzval)
859	6447780	do K=2,nz ; do j=js,je ; do i=is,ie
860	6393600	I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))
861	6393600	ent_s(i,j,K) = dt * I_dzval * Kd_int(i,j,K)
862	6446880	ent_t(i,j,K) = ent_s(i,j,K)
863		enddo ; enddo ; enddo
864	12	if (showCallTree) call callTree_waypoint("done setting ent_s and ent_t from Kd_int (diabatic)")
865
866	12	if (CS%debug) then
867	0	call MOM_forcing_chksum("after calc_entrain ", fluxes, G, US, haloshift=0)
868	0	call MOM_thermovar_chksum("after calc_entrain ", tv, G, US)
869	0	call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, US, haloshift=0)
870	0	call hchksum(ent_s, "after calc_entrain ent_s", G%HI, haloshift=0, unscale=GV%H_to_MKS)
871		endif
872
873		! Save fields before boundary forcing is applied for tendency diagnostics
874	6534912	do k=1,nz ; do j=js,je ; do i=is,ie
875	6534000	h_orig(i,j,k) = h(i,j,k)
876		enddo ; enddo ; enddo
877	12	if (CS%boundary_forcing_tendency_diag) then
878	0	do k=1,nz ; do j=js,je ; do i=is,ie
879	0	temp_diag(i,j,k) = tv%T(i,j,k)
880	0	saln_diag(i,j,k) = tv%S(i,j,k)
881		enddo ; enddo ; enddo
882		endif
883
884		! Apply forcing
885		! Changes made to following fields: h, tv%T and tv%S.
886	12	call cpu_clock_begin(id_clock_remap)
887
888	12	if (CS%use_energetic_PBL) then
889
890	105276	skinbuoyflux(:,:) = 0.0
891		call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, &
892		optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, &
893	12	CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD_h=visc%h_ML)
894
895	12	if (CS%debug) then
896	0	call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, unscale=GV%H_to_mks)
897	0	call hchksum(ent_s, "after applyBoundaryFluxes ent_s", G%HI, haloshift=0, unscale=GV%H_to_mks)
898		call hchksum(cTKE, "after applyBoundaryFluxes cTKE", G%HI, haloshift=0, &
899	0	unscale=US%RZ3_T3_to_W_m2*US%T_to_s)
900		call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT", G%HI, haloshift=0, &
901	0	unscale=US%kg_m3_to_R*US%degC_to_C)
902		call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS", G%HI, haloshift=0, &
903	0	unscale=US%kg_m3_to_R*US%ppt_to_S)
904	0	call hchksum(h, "after applyBoundaryFluxes h", G%HI, haloshift=0, unscale=GV%H_to_mks)
905	0	call hchksum(tv%T, "after applyBoundaryFluxes tv%T", G%HI, haloshift=0, unscale=US%C_to_degC)
906	0	call hchksum(tv%S, "after applyBoundaryFluxes tv%S", G%HI, haloshift=0, unscale=US%S_to_ppt)
907		call hchksum(SkinBuoyFlux, "after applyBdryFlux SkinBuoyFlux", G%HI, haloshift=0, &
908	0	unscale=US%Z_to_m*2US%s_to_T**3)
909		endif
910
911	12	call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
912		call energetic_PBL(h, u_h, v_h, tv, fluxes, visc, dt, Kd_ePBL, G, GV, US, &
913	12	CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, BBL_BuoyFlux, waves=waves)
914
915	12	call energetic_PBL_get_MLD(CS%ePBL, BLD(:,:), G, US)
916		! If visc%MLD or visc%h_ML exist, copy ePBL's BLD into them with appropriate conversions.
917	12	if (associated(visc%h_ML)) call convert_MLD_to_ML_thickness(BLD, h, visc%h_ML, tv, G, GV)
918	105276	if (associated(visc%MLD)) visc%MLD(:,:) = BLD(:,:)
919	105276	if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = SkinBuoyFlux(:,:)
920
921		! Find the vertical distances across layers, which may have been modified by the net surface flux
922	12	call thickness_to_dz(h, tv, dz, G, GV, US)
923
924		! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
925	6447780	do K=2,nz ; do j=js,je ; do i=is,ie
926	6393600	if (CS%ePBL_is_additive) then
927	0	Kd_add_here = Kd_ePBL(i,j,K)
928	0	visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + CS%ePBL_Prandtl*Kd_ePBL(i,j,K)
929		else
930	6393600	Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_shear(i,j,K), 0.0)
931	6393600	visc%Kv_shear(i,j,K) = max(visc%Kv_shear(i,j,K), CS%ePBL_Prandtl*Kd_ePBL(i,j,K))
932		endif
933
934	6393600	Ent_int = Kd_add_here * dt / (0.5*(dz(i,j,k-1) + dz(i,j,k)) + dz_neglect)
935	6393600	ent_s(i,j,K) = ent_s(i,j,K) + Ent_int
936	6393600	Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here
937
938		! for diagnostics
939	6393600	Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_int(i,j,K)
940	6446880	Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_int(i,j,K)
941
942		enddo ; enddo ; enddo
943
944	12	if (CS%debug) then
945	0	call hchksum(ent_t, "after ePBL ent_t", G%HI, haloshift=0, unscale=GV%H_to_MKS)
946	0	call hchksum(ent_s, "after ePBL ent_s", G%HI, haloshift=0, unscale=GV%H_to_MKS)
947	0	call hchksum(Kd_ePBL, "after ePBL Kd_ePBL", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
948		endif
949
950		else
951		call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, &
952		optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, &
953	0	CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD_h=visc%h_ML)
954
955		! Find the vertical distances across layers, which may have been modified by the net surface flux
956	0	call thickness_to_dz(h, tv, dz, G, GV, US)
957
958		endif ! endif for CS%use_energetic_PBL
959
960		! diagnose the tendencies due to boundary forcing
961		! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme
962		! so all tendency diagnostics need to be posted on h_orig, and grids rebuilt afterwards
963	12	if (CS%boundary_forcing_tendency_diag) then
964	0	call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_orig, dt, G, GV, US, CS)
965	0	if (CS%id_boundary_forcing_h > 0) call post_data(CS%id_boundary_forcing_h, h, CS%diag, alt_h=h_orig)
966		endif
967		! Boundary fluxes may have changed T, S, and h
968	12	call diag_update_remap_grids(CS%diag)
969	12	call cpu_clock_end(id_clock_remap)
970	12	if (CS%debug) then
971	0	call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, US, haloshift=0)
972	0	call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G, US)
973	0	call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, US, haloshift=0)
974		endif
975	12	if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)")
976	12	if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, G, GV, US)
977
978	12	if (CS%debug) then
979	0	call MOM_state_chksum("after negative check ", u, v, h, G, GV, US, haloshift=0)
980	0	call MOM_forcing_chksum("after negative check ", fluxes, G, US, haloshift=0)
981	0	call MOM_thermovar_chksum("after negative check ", tv, G, US)
982		endif
983	12	if (showCallTree) call callTree_waypoint("done with h=ea-eb (diabatic)")
984	12	if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G, GV, US)
985
986		! calculate change in temperature & salinity due to dia-coordinate surface diffusion
987	12	if (associated(tv%T)) then
988
989	12	if (CS%debug) then
990	0	call hchksum(ent_t, "before triDiagTS ent_t ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
991	0	call hchksum(ent_s, "before triDiagTS ent_s ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
992		endif
993
994	12	call cpu_clock_begin(id_clock_tridiag)
995		! Keep salinity from falling below a small but positive threshold.
996		! This constraint is needed for SIS1 ice model, which can extract
997		! more salt than is present in the ocean. SIS2 does not suffer
998		! from this limitation, in which case we can let salinity=0 and still
999		! have salt conserved with SIS2 ice. So for SIS2, we can run with
1000		! BOUND_SALINITY=False in MOM.F90.
1001	12	if (associated(tv%S) .and. associated(tv%salt_deficit)) &
1002	0	call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp)
1003
1004	12	if (CS%diabatic_diff_tendency_diag) then
1005	0	do k=1,nz ; do j=js,je ; do i=is,ie
1006	0	temp_diag(i,j,k) = tv%T(i,j,k)
1007	0	saln_diag(i,j,k) = tv%S(i,j,k)
1008		enddo ; enddo ; enddo
1009		endif
1010
1011		! Changes T and S via the tridiagonal solver; no change to h
1012	6534912	do k=1,nz ; do j=js,je ; do i=is,ie
1013	6534000	ent_t(i,j,K) = ent_s(i,j,K) ; ent_t(i,j,K+1) = ent_s(i,j,K+1)
1014		enddo ; enddo ; enddo
1015	12	if (CS%tracer_tridiag) then
1016	0	call tracer_vertdiff_Eulerian(h, ent_t, dt, tv%T, G, GV)
1017	0	call tracer_vertdiff_Eulerian(h, ent_s, dt, tv%S, G, GV)
1018		else
1019	12	call triDiagTS_Eulerian(G, GV, is, ie, js, je, h, ent_s, tv%T, tv%S)
1020		endif
1021
1022		! diagnose temperature, salinity, heat, and salt tendencies
1023	12	if (CS%diabatic_diff_tendency_diag) then
1024	0	call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, G, GV, US, CS)
1025	0	if (CS%id_diabatic_diff_h > 0) call post_data(CS%id_diabatic_diff_h, h, CS%diag, alt_h=h)
1026		endif
1027
1028	12	call cpu_clock_end(id_clock_tridiag)
1029
1030	12	if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)")
1031
1032		endif ! endif corresponding to if (associated(tv%T))
1033
1034	12	if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G, GV, US)
1035
1036	12	if (CS%debug) then
1037	0	call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, US, haloshift=0)
1038	0	call MOM_thermovar_chksum("after mixed layer ", tv, G, US)
1039		endif
1040
1041		! Whenever thickness changes let the diag manager know, as the
1042		! target grids for vertical remapping may need to be regenerated.
1043	12	call diag_update_remap_grids(CS%diag)
1044
1045		! Set diffusivities for VBF diagnostics if enabled
1046	12	if (CS%use_energetic_PBL .and. associated(CS%VBF%Kd_ePBL)) CS%VBF%Kd_ePBL(:,:,:) = Kd_ePBL(:,:,:)
1047	12	if (associated(CS%VBF%Kd_salt)) CS%VBF%Kd_temp(:,:,:) = Kd_heat(:,:,:)
1048	12	if (associated(CS%VBF%Kd_temp)) CS%VBF%Kd_salt(:,:,:) = Kd_salt(:,:,:)
1049
1050
1051		! Diagnose the diapycnal diffusivities and other related quantities.
1052	12	if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_int, CS%diag)
1053	12	if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag)
1054	12	if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag)
1055	12	if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag)
1056
1057	12	if (CS%id_ea > 0) call post_data(CS%id_ea, ent_s(:,:,1:nz), CS%diag)
1058	12	if (CS%id_eb > 0) call post_data(CS%id_eb, ent_s(:,:,2:nz+1), CS%diag)
1059	12	if (CS%id_ea_t > 0) call post_data(CS%id_ea_t, ent_t(:,:,1:nz), CS%diag)
1060	12	if (CS%id_eb_t > 0) call post_data(CS%id_eb_t, ent_t(:,:,2:nz+1), CS%diag)
1061	12	if (CS%id_ea_s > 0) call post_data(CS%id_ea_s, ent_s(:,:,1:nz), CS%diag)
1062	12	if (CS%id_eb_s > 0) call post_data(CS%id_eb_s, ent_s(:,:,2:nz+1), CS%diag)
1063	12	Idt = 1.0 / dt
1064	12	if (CS%id_Tdif > 0) then
1065	0	do j=js,je ; do i=is,ie
1066	0	Tdif_flx(i,j,1) = 0.0 ; Tdif_flx(i,j,nz+1) = 0.0
1067		enddo ; enddo
1068		!$OMP parallel do default(shared)
1069	0	do K=2,nz ; do j=js,je ; do i=is,ie
1070	0	Tdif_flx(i,j,K) = (Idt * ent_t(i,j,K)) * (tv%T(i,j,k-1) - tv%T(i,j,k))
1071		enddo ; enddo ; enddo
1072	0	if (CS%id_Tdif > 0) call post_data(CS%id_Tdif, Tdif_flx, CS%diag)
1073		endif
1074	12	if (CS%id_Sdif > 0) then
1075	0	do j=js,je ; do i=is,ie
1076	0	Sdif_flx(i,j,1) = 0.0 ; Sdif_flx(i,j,nz+1) = 0.0
1077		enddo ; enddo
1078		!$OMP parallel do default(shared)
1079	0	do K=2,nz ; do j=js,je ; do i=is,ie
1080	0	Sdif_flx(i,j,K) = (Idt * ent_s(i,j,K)) * (tv%S(i,j,k-1) - tv%S(i,j,k))
1081		enddo ; enddo ; enddo
1082	0	if (CS%id_Sdif > 0) call post_data(CS%id_Sdif, Sdif_flx, CS%diag)
1083		endif
1084
1085	12	if (CS%Use_KdWork_diag .or. CS%Use_N2_diag) then
1086	0	N2_salt(:,:,:) = 0.0
1087	0	N2_temp(:,:,:) = 0.0
1088		!Compute N2 and don't mask negatives here
1089	0	EOSdom(:) = EOS_domain(G%HI)
1090	0	if (nonBous) then
1091		!$OMP parallel do default(shared)
1092	0	do j=js,je
1093	0	if (associated(tv%p_surf)) then
1094	0	do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo
1095		else
1096	0	do i=is,ie ; p_i(i) = 0.0 ; enddo
1097		endif
1098	0	do K=2,nz
1099	0	do i=is,ie
1100	0	p_i(i) = p_i(i) + H_to_pres * h(i,j,k-1)
1101		enddo
1102	0	T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))
1103	0	S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))
1104	0	call calculate_specific_vol_derivs(T_i, S_i, p_i, dSpV_dT, dSpV_dS, tv%eqn_of_state, EOSdom)
1105	0	do i=is,ie
1106	0	I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))
1107	0	N2_salt(i,j,K) = (tv%S(i,j,k-1) - tv%S(i,j,k)) * (dSpv_dS(i) * (alt_H_to_pres * I_dzval))
1108	0	N2_temp(i,j,K) = (tv%T(i,j,k-1) - tv%T(i,j,k)) * (dSpV_dT(i) * (alt_H_to_pres * I_dzval))
1109		enddo
1110		enddo
1111		enddo
1112		else
1113		!$OMP parallel do default(shared)
1114	0	do j=js,je
1115	0	if (associated(tv%p_surf)) then
1116	0	do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo
1117		else
1118	0	do i=is,ie ; p_i(i) = 0.0 ; enddo
1119		endif
1120	0	do K=2,nz
1121	0	do i=is,ie
1122	0	p_i(i) = p_i(i) + H_to_pres* h(i,j,k-1)
1123		enddo
1124	0	T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))
1125	0	S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))
1126	0	call calculate_density_derivs(T_i, S_i, p_i, dRhodT, dRhodS, tv%eqn_of_state, EOSdom)
1127	0	do i=is,ie
1128	0	I_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
1129	0	N2_salt(i,j,K) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (dRhodS(i) * (g_rho0 * I_h))
1130	0	N2_temp(i,j,K) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (dRhodT(i) * (g_rho0 * I_h))
1131		enddo
1132		enddo
1133		enddo
1134		endif
1135	0	if (CS%id_N2_dd>0) call post_data(CS%id_N2_dd, N2_salt(:,:,:)+N2_temp(:,:,:), CS%diag)
1136	0	if (CS%id_N2_salt_dd>0) call post_data(CS%id_N2_salt_dd, N2_salt, CS%diag)
1137	0	if (CS%id_N2_temp_dd>0) call post_data(CS%id_N2_temp_dd, N2_temp, CS%diag)
1138
1139	0	if (CS%Use_KdWork_diag) then
1140	0	call KdWork_diagnostics(G,GV,US,CS%diag,CS%VBF,N2_salt,N2_temp,dz)
1141		endif
1142
1143	0	call deallocate_VBF_CS(CS%VBF)
1144
1145		endif
1146
1147		! mixing of passive tracers from massless boundary layers to interior
1148	12	call cpu_clock_begin(id_clock_tracers)
1149
1150	12	if (CS%mix_boundary_tracer_ALE) then
1151	0	Tr_ea_BBL = sqrt(dt * CS%Kd_BBL_tr)
1152
1153		!$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
1154	0	do j=js,je
1155	0	do i=is,ie
1156	0	htot(i) = 0.0
1157	0	in_boundary(i) = (G%mask2dT(i,j) > 0.0)
1158		enddo
1159	0	do k=nz,2,-1 ; do i=is,ie
1160	0	if (in_boundary(i)) then
1161	0	htot(i) = htot(i) + h(i,j,k)
1162		! If diapycnal mixing has been suppressed because this is a massless
1163		! layer near the bottom, add some mixing of tracers between these
1164		! layers. This flux is based on the harmonic mean of the two
1165		! thicknesses, as this corresponds pretty closely (to within
1166		! differences in the density jumps between layers) with what is done
1167		! in the calculation of the fluxes in the first place. Kd_min_tr
1168		! should be much less than the values that have been set in Kd_int,
1169		! perhaps a molecular diffusivity.
1170		add_ent = ((dt * CS%Kd_min_tr)) * &
1171		((dz(i,j,k-1)+dz(i,j,k)+dz_neglect) / (dz(i,j,k-1)*dz(i,j,k)+dz_neglect2)) - &
1172	0	0.5*(ent_s(i,j,K) + ent_s(i,j,K))
1173	0	if (htot(i) < Tr_ea_BBL) then
1174	0	add_ent = max(0.0, add_ent, (Tr_ea_BBL - htot(i)) - ent_s(i,j,K))
1175	0	elseif (add_ent < 0.0) then
1176	0	add_ent = 0.0 ; in_boundary(i) = .false.
1177		endif
1178
1179	0	ent_s(i,j,K) = ent_s(i,j,K) + add_ent
1180		endif
1181
1182	0	if (CS%double_diffuse) then ; if (Kd_extra_S(i,j,k) > 0.0) then
1183		add_ent = (dt * Kd_extra_S(i,j,k)) / &
1184	0	(0.5 * (dz(i,j,k-1) + dz(i,j,k)) + dz_neglect)
1185	0	ent_s(i,j,K) = ent_s(i,j,K) + add_ent
1186		endif ; endif
1187		enddo ; enddo
1188
1189		enddo
1190	12	elseif (CS%double_diffuse .and. .not.CS%mix_boundary_tracers) then ! extra diffusivity for passive tracers
1191		!$OMP parallel do default(shared) private(add_ent)
1192	0	do k=nz,2,-1 ; do j=js,je ; do i=is,ie
1193	0	if (Kd_extra_S(i,j,k) > 0.0) then
1194		add_ent = (dt * Kd_extra_S(i,j,k)) / &
1195	0	(0.5 * (dz(i,j,k-1) + dz(i,j,k)) + dz_neglect)
1196		else
1197	0	add_ent = 0.0
1198		endif
1199	0	ent_s(i,j,K) = ent_s(i,j,K) + add_ent
1200		enddo ; enddo ; enddo
1201		endif ! (CS%mix_boundary_tracers)
1202
1203		! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
1204		call call_tracer_column_fns(h_orig, h, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), fluxes, BLD, dt, &
1205		G, GV, US, tv, CS%optics, CS%tracer_flow_CSp, CS%debug, &
1206		KPP_CSp=CS%KPP_CSp, &
1207		nonLocalTrans=KPP_NLTscalar, &
1208		evap_CFL_limit=CS%evap_CFL_limit, &
1209	12	minimum_forcing_depth=CS%minimum_forcing_depth, h_BL=visc%h_ML)
1210
1211	12	call cpu_clock_end(id_clock_tracers)
1212
1213		! Apply ALE sponge
1214	12	if (CS%use_sponge .and. associated(CS%ALE_sponge_CSp)) then
1215	0	call cpu_clock_begin(id_clock_sponge)
1216	0	call apply_ALE_sponge(h, tv, dt, G, GV, US, CS%ALE_sponge_CSp, CS%Time)
1217	0	call cpu_clock_end(id_clock_sponge)
1218	0	if (CS%debug) then
1219	0	call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, US, haloshift=0)
1220	0	call MOM_thermovar_chksum("apply_sponge ", tv, G, US)
1221		endif
1222		endif ! CS%use_sponge
1223
1224		! Apply data assimilation incremental update -oda_incupd-
1225	12	if (CS%use_oda_incupd .and. associated(CS%oda_incupd_CSp)) then
1226	0	call MOM_mesg("Starting ODA_INCUPD legacy ", 5)
1227	0	call cpu_clock_begin(id_clock_oda_incupd)
1228	0	call apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS%oda_incupd_CSp)
1229	0	call cpu_clock_end(id_clock_oda_incupd)
1230	0	if (CS%debug) then
1231	0	call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0)
1232	0	call MOM_thermovar_chksum("apply_oda_incupd ", tv, G, US)
1233		endif
1234		endif ! CS%use_oda_incupd
1235
1236
1237
1238	12	call disable_averaging(CS%diag)
1239
1240	12	if (showCallTree) call callTree_leave("diabatic_ALE_legacy()")
1241
1242	12	end subroutine diabatic_ALE_legacy
1243
1244
1245		!> This subroutine imposes the diapycnal mass fluxes and the
1246		!! accompanying diapycnal advection of momentum and tracers.
1247	0	subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
1248		G, GV, US, CS, stoch_CS, Waves)
1249		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1250		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1251		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1252		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1253		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
1254		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2]
1255		type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
1256		!! unused have NULL ptrs
1257		real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Active mixed layer depth [Z ~> m]
1258		type(forcing), intent(inout) :: fluxes !< points to forcing fields
1259		!! unused fields have NULL ptrs
1260		type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities,
1261		!! BBL properties and related fields
1262		type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum
1263		!! equations, to enable the later derived
1264		!! diagnostics, like energy budgets
1265		type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations
1266		real, intent(in) :: dt !< time increment [T ~> s]
1267		type(time_type), intent(in) :: Time_end !< Time at the end of the interval
1268		type(diabatic_CS), pointer :: CS !< module control structure
1269		type(stochastic_CS), pointer :: stoch_CS !< stochastic control structure
1270		type(Wave_parameters_CS), pointer :: Waves !< Surface gravity waves
1271
1272		! local variables
1273		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
1274	0	h_orig, & ! Initial layer thicknesses [H ~> m or kg m-2]
1275	0	dz, & ! The vertical distance between interfaces around a layer [Z ~> m]
1276	0	dSV_dT, & ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
1277	0	dSV_dS, & ! The partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
1278	0	cTKE, & ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
1279	0	u_h, & ! Zonal velocities interpolated to thickness points [L T-1 ~> m s-1]
1280	0	v_h, & ! Meridional velocities interpolated to thickness points [L T-1 ~> m s-1]
1281	0	temp_diag, & ! Diagnostic array of previous temperatures [C ~> degC]
1282	0	saln_diag ! Diagnostic array of previous salinity [S ~> ppt]
1283
1284		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
1285	0	ent_s, & ! The diffusive coupling across interfaces within one time step for
1286		! salinity and passive tracers [H ~> m or kg m-2]
1287	0	ent_t, & ! The diffusive coupling across interfaces within one time step for
1288		! temperature [H ~> m or kg m-2]
1289	0	Kd_heat, & ! diapycnal diffusivity of heat or the smaller of the diapycnal diffusivities of
1290		! heat and salt [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1291	0	Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1292	0	Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to
1293		! Kd_int returned from set_diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1294	0	Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to
1295		! Kd_int returned from set_diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1296	0	Kd_ePBL, & ! boundary layer or convective diapycnal diffusivities at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1297	0	KPP_NLTheat, & ! KPP non-local transport for heat [nondim]
1298	0	KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim]
1299	0	KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]
1300	0	Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1301	0	Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1302	0	N2_salt, & !< Salinity contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]
1303	0	N2_temp !< Temperature contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]
1304
1305		real, dimension(SZI_(G),SZJ_(G)) :: &
1306	0	U_star, & ! The friction velocity [Z T-1 ~> m s-1].
1307	0	KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1308	0	KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1309	0	SkinBuoyFlux, & ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
1310	0	BBL_BuoyFlux ! 2d bottom buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
1311
1312		logical, dimension(SZI_(G)) :: &
1313	0	in_boundary ! True if there are no massive layers below, where massive is defined as
1314		! sufficiently thick that the no-flux boundary conditions have not restricted
1315		! the entrainment - usually sqrt(Kd*dt).
1316
1317		real, dimension(SZI_(G)) :: &
1318	0	p_i ,& ! Pressure at the interface [R L2 T-2 ~> Pa]
1319	0	T_i, & ! Temperature at the interface [C ~> degC]
1320	0	S_i, & ! Salinity at the interface [S ~> ppt]
1321	0	drhodS, & ! Local change in density w.r.t. salinity using model EOS & state [R C-1 ~> kg m-3 ppt-1]
1322	0	drhodT, & ! Local change in density w.r.t. temperature using model EOS & state [R C-1 ~> kg m-3 degC-1]
1323	0	dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
1324	0	dSpV_dS ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
1325
1326		real :: h_neglect ! A thickness that is so small it is usually lost
1327		! in roundoff and can be neglected [H ~> m or kg m-2]
1328		real :: dz_neglect ! A vertical distance that is so small it is usually lost
1329		! in roundoff and can be neglected [Z ~> m]
1330		real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2]
1331		real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2]
1332		real :: I_dzval ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1]
1333		real :: I_h ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1]
1334		real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is
1335		! coupled to the bottom within a timestep [H ~> m or kg m-2]
1336	0	real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].
1337		real :: Kd_add_here ! An added diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1338		real :: Idt ! The inverse time step [T-1 ~> s-1]
1339		real :: g_Rho0 ! G_Earth/Rho0 [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2]
1340		real :: H_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
1341		real :: alt_H_to_pres! A conversion factor from thicknesses to pressure w/ alternative scaling [R Z T-2 ~> Pa m-1]
1342		logical :: nonBous ! True if not using the Boussinesq approximation
1343
1344		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1345
1346		logical :: showCallTree ! If true, show the call tree
1347		integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, Isq, Ieq, Jsq, Jeq, nz
1348
1349	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1350	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1351	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1352	0	dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect
1353	0	h_neglect = GV%H_subroundoff
1354
1355	0	nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
1356	0	g_Rho0 = GV%g_Earth_Z_T2 / GV%H_to_RZ
1357	0	H_to_pres = GV%H_to_RZ * GV%g_Earth
1358	0	alt_H_to_pres = H_to_pres * US%L_to_Z*2 GV%Z_to_H
1359
1360	0	Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0
1361	0	ent_s(:,:,:) = 0.0 ; ent_t(:,:,:) = 0.0
1362
1363	0	showCallTree = callTree_showQuery()
1364	0	if (showCallTree) call callTree_enter("diabatic_ALE(), MOM_diabatic_driver.F90")
1365
1366	0	if (.not. (CS%useALEalgorithm)) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
1367	0	"The ALE algorithm must be enabled when using MOM_diabatic_driver.")
1368
1369		! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
1370	0	call enable_averages(dt, Time_end, CS%diag)
1371
1372	0	if (CS%Use_KdWork_diag) call Allocate_VBF_CS(G, GV, CS%VBF)
1373
1374	0	if (CS%use_geothermal) then
1375	0	call cpu_clock_begin(id_clock_geothermal)
1376	0	call geothermal_in_place(h, tv, dt, G, GV, US, CS%geothermal, BBL_buoyflux, halo=CS%halo_TS_diff)
1377	0	call cpu_clock_end(id_clock_geothermal)
1378	0	if (showCallTree) call callTree_waypoint("geothermal (diabatic)")
1379	0	if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G, GV, US)
1380		endif
1381
1382		! Whenever thickness changes let the diag manager know, target grids
1383		! for vertical remapping may need to be regenerated.
1384	0	call diag_update_remap_grids(CS%diag)
1385
1386		! Set_pen_shortwave estimates the optical properties of the water column.
1387		! It will need to be modified later to include information about the
1388		! biological properties and layer thicknesses.
1389	0	if (associated(CS%optics)) &
1390	0	call set_pen_shortwave(CS%optics, fluxes, G, GV, US, CS%diabatic_aux_CSp, CS%opacity, CS%tracer_flow_CSp)
1391
1392	0	if (CS%debug) call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, US, haloshift=0)
1393
1394	0	if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then
1395	0	if (CS%use_geothermal) then
1396	0	call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, zero_mix=.true.)
1397		else
1398	0	call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
1399		endif
1400	0	if (showCallTree) call callTree_waypoint("done with find_uv_at_h (diabatic)")
1401		endif
1402
1403	0	call cpu_clock_begin(id_clock_set_diffusivity)
1404		! Sets: Kd_heat, Kd_extra_T, Kd_extra_S and visc%TKE_turb
1405		! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow
1406	0	if (CS%debug) &
1407	0	call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff)
1408	0	if (CS%double_diffuse) then
1409		call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_heat, G, GV, US, &
1410	0	CS%set_diff_CSp, CS%VBF, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S)
1411		else
1412		call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_heat, G, GV, US, &
1413	0	CS%set_diff_CSp, CS%VBF)
1414		endif
1415	0	call cpu_clock_end(id_clock_set_diffusivity)
1416	0	if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)")
1417
1418	0	if (CS%debug) then
1419	0	call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, US, haloshift=0)
1420	0	call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, US, haloshift=0)
1421	0	call MOM_thermovar_chksum("after set_diffusivity ", tv, G, US)
1422	0	call hchksum(Kd_heat, "after set_diffusivity Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
1423		endif
1424
1425		! Set diffusivities for heat and salt separately, and possibly change the meaning of Kd_heat.
1426	0	if (CS%double_diffuse) then
1427		! Add contributions from double diffusion
1428		!$OMP parallel do default(shared)
1429	0	do K=1,nz+1 ; do j=js,je ; do i=is,ie
1430	0	Kd_salt(i,j,K) = Kd_heat(i,j,K) + Kd_extra_S(i,j,K)
1431	0	Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_extra_T(i,j,K)
1432		enddo ; enddo ; enddo
1433		else
1434		!$OMP parallel do default(shared)
1435	0	do K=1,nz+1 ; do j=js,je ; do i=is,ie
1436	0	Kd_salt(i,j,K) = Kd_heat(i,j,K)
1437		enddo ; enddo ; enddo
1438		endif
1439
1440	0	if (CS%debug) then
1441	0	call hchksum(Kd_heat, "after double diffuse Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
1442	0	call hchksum(Kd_salt, "after double diffuse Kd_salt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
1443		endif
1444
1445	0	if (CS%useKPP) then
1446	0	call cpu_clock_begin(id_clock_kpp)
1447		! total vertical viscosity in the interior is represented via visc%Kv_shear
1448	0	do k=1,nz+1 ; do j=js,je ; do i=is,ie
1449	0	visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k)
1450		enddo ; enddo ; enddo
1451
1452		! NOTE: The following do not require initialization, but their checksums do
1453		! require initialization, and past versions were initialized to zero.
1454	0	KPP_NLTheat(:,:,:) = 0.
1455	0	KPP_NLTscalar(:,:,:) = 0.
1456	0	KPP_buoy_flux(:,:,:) = 0.
1457	0	KPP_temp_flux(:,:) = 0.
1458	0	KPP_salt_flux(:,:) = 0.
1459
1460		! KPP needs the surface buoyancy flux but does not update state variables.
1461		! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
1462		! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux
1463		! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second)
1464		! unlike other instances where the fluxes are integrated in time over a time-step.
1465		call calculateBuoyancyFlux2d(G, GV, US, fluxes, CS%optics, h, tv%T, tv%S, tv, &
1466	0	KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux)
1467
1468		! Determine the friction velocity, perhaps using the evovling surface density.
1469	0	call find_ustar(fluxes, tv, U_star, G, GV, US)
1470
1471		! The KPP scheme calculates boundary layer diffusivities and non-local transport.
1472	0	if ( associated(fluxes%lamult) ) then
1473		call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, &
1474	0	U_star, KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult)
1475
1476		call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, &
1477	0	Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult)
1478		else
1479		call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, &
1480	0	U_star, KPP_buoy_flux, Waves=Waves)
1481
1482		call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, &
1483	0	Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves)
1484		endif
1485
1486	0	call KPP_get_BLD(CS%KPP_CSp, BLD(:,:), G, US)
1487		! If visc%MLD or visc%h_ML exist, copy KPP's BLD into them with appropriate conversions.
1488	0	if (associated(visc%h_ML)) call convert_MLD_to_ML_thickness(BLD, h, visc%h_ML, tv, G, GV)
1489	0	if (associated(visc%MLD)) visc%MLD(:,:) = BLD(:,:)
1490	0	if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = KPP_buoy_flux(:,:,1)
1491
1492	0	if (showCallTree) call callTree_waypoint("done with KPP_calculate (diabatic)")
1493	0	if (CS%debug) then
1494	0	call MOM_state_chksum("after KPP", u, v, h, G, GV, US, haloshift=0)
1495	0	call MOM_forcing_chksum("after KPP", fluxes, G, US, haloshift=0)
1496	0	call MOM_thermovar_chksum("after KPP", tv, G, US)
1497	0	call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
1498	0	call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
1499		call hchksum(KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, &
1500	0	unscale=US%C_to_degCGV%H_to_mUS%s_to_T)
1501		call hchksum(KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, &
1502	0	unscale=US%S_to_pptGV%H_to_mUS%s_to_T)
1503	0	call hchksum(KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0)
1504	0	call hchksum(KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0)
1505		endif
1506		! Apply non-local transport of heat and salt
1507		! Changes: tv%T, tv%S
1508		call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, KPP_NLTheat, KPP_temp_flux, &
1509	0	dt, tv%tr_T, tv%T, tv%C_p)
1510		call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, KPP_NLTscalar, KPP_salt_flux, &
1511	0	dt, tv%tr_S, tv%S)
1512	0	call cpu_clock_end(id_clock_kpp)
1513	0	if (showCallTree) call callTree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
1514	0	if (CS%debugConservation) call MOM_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, G, GV, US)
1515
1516	0	if (CS%debug) then
1517	0	call MOM_state_chksum("after KPP_applyNLT ", u, v, h, G, GV, US, haloshift=0)
1518	0	call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, US, haloshift=0)
1519	0	call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G, US)
1520		endif
1521		endif ! endif for KPP
1522
1523		! Calculate vertical mixing due to convection (computed via CVMix)
1524	0	if (CS%use_CVMix_conv) then
1525		! Increment vertical diffusion and viscosity due to convection
1526	0	call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv, BLD, Kd_heat, visc%Kv_shear, Kd_aux=Kd_salt)
1527		endif
1528
1529		! Save fields before boundary forcing is applied for tendency diagnostics
1530	0	do k=1,nz ; do j=js,je ; do i=is,ie
1531	0	h_orig(i,j,k) = h(i,j,k)
1532		enddo ; enddo ; enddo
1533	0	if (CS%boundary_forcing_tendency_diag) then
1534	0	do k=1,nz ; do j=js,je ; do i=is,ie
1535	0	temp_diag(i,j,k) = tv%T(i,j,k)
1536	0	saln_diag(i,j,k) = tv%S(i,j,k)
1537		enddo ; enddo ; enddo
1538		endif
1539
1540		! Apply forcing
1541		! Changes made to following fields: h, tv%T and tv%S.
1542	0	call cpu_clock_begin(id_clock_remap)
1543
1544	0	if (CS%use_energetic_PBL) then
1545
1546	0	skinbuoyflux(:,:) = 0.0
1547		call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, &
1548		optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, &
1549	0	CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD_h=visc%h_ML)
1550
1551	0	if (CS%debug) then
1552	0	call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, unscale=GV%H_to_MKS)
1553	0	call hchksum(ent_s, "after applyBoundaryFluxes ent_s", G%HI, haloshift=0, unscale=GV%H_to_MKS)
1554		call hchksum(cTKE, "after applyBoundaryFluxes cTKE", G%HI, haloshift=0, &
1555	0	unscale=US%RZ3_T3_to_W_m2*US%T_to_s)
1556		call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT", G%HI, haloshift=0, &
1557	0	unscale=US%kg_m3_to_R*US%degC_to_C)
1558		call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS", G%HI, haloshift=0, &
1559	0	unscale=US%kg_m3_to_R*US%ppt_to_S)
1560		endif
1561
1562	0	call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
1563		call energetic_PBL(h, u_h, v_h, tv, fluxes, visc, dt, Kd_ePBL, G, GV, US, &
1564	0	CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, BBL_BuoyFlux, waves=waves)
1565
1566	0	call energetic_PBL_get_MLD(CS%ePBL, BLD(:,:), G, US)
1567		! If visc%MLD or visc%h_ML exist, copy ePBL's BLD into them with appropriate conversions.
1568	0	if (associated(visc%h_ML)) call convert_MLD_to_ML_thickness(BLD, h, visc%h_ML, tv, G, GV)
1569	0	if (associated(visc%MLD)) visc%MLD(:,:) = BLD(:,:)
1570	0	if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = SkinBuoyFlux(:,:)
1571
1572		! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
1573	0	do K=2,nz ; do j=js,je ; do i=is,ie
1574	0	if (CS%ePBL_is_additive) then
1575	0	Kd_add_here = Kd_ePBL(i,j,K)
1576	0	visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + CS%ePBL_Prandtl*Kd_ePBL(i,j,K)
1577		else
1578	0	Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_shear(i,j,K), 0.0)
1579	0	visc%Kv_shear(i,j,K) = max(visc%Kv_shear(i,j,K), CS%ePBL_Prandtl*Kd_ePBL(i,j,K))
1580		endif
1581
1582	0	Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_add_here
1583	0	Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_add_here
1584		enddo ; enddo ; enddo
1585
1586	0	if (CS%debug) then
1587	0	call hchksum(ent_t, "after ePBL ent_t", G%HI, haloshift=0, unscale=GV%H_to_MKS)
1588	0	call hchksum(ent_s, "after ePBL ent_s", G%HI, haloshift=0, unscale=GV%H_to_MKS)
1589	0	call hchksum(Kd_ePBL, "after ePBL Kd_ePBL", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
1590		endif
1591
1592		else
1593		call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, &
1594		optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, &
1595	0	CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD_h=visc%h_ML)
1596
1597		endif ! endif for CS%use_energetic_PBL
1598
1599		! diagnose the tendencies due to boundary forcing
1600		! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme
1601		! so all tendency diagnostics need to be posted on h_orig, and grids rebuilt afterwards
1602	0	if (CS%boundary_forcing_tendency_diag) then
1603	0	call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_orig, dt, G, GV, US, CS)
1604	0	if (CS%id_boundary_forcing_h > 0) call post_data(CS%id_boundary_forcing_h, h, CS%diag, alt_h=h_orig)
1605		endif
1606		! Boundary fluxes may have changed T, S, and h
1607	0	call diag_update_remap_grids(CS%diag)
1608	0	call cpu_clock_end(id_clock_remap)
1609	0	if (CS%debug) then
1610	0	call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, US, haloshift=0)
1611	0	call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G, US)
1612	0	call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, US, haloshift=0)
1613		endif
1614	0	if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)")
1615	0	if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, G, GV, US)
1616
1617	0	if (showCallTree) call callTree_waypoint("done with h=ea-eb (diabatic)")
1618	0	if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G, GV, US)
1619
1620		! calculate change in temperature & salinity due to dia-coordinate surface diffusion
1621	0	if (associated(tv%T)) then
1622
1623	0	if (CS%debug) then
1624	0	call hchksum(ent_t, "before triDiagTS ent_t ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
1625	0	call hchksum(ent_s, "before triDiagTS ent_s ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
1626		endif
1627
1628	0	call cpu_clock_begin(id_clock_tridiag)
1629		! Keep salinity from falling below a small but positive threshold.
1630		! This constraint is needed for SIS1 ice model, which can extract
1631		! more salt than is present in the ocean. SIS2 does not suffer
1632		! from this limitation, in which case we can let salinity=0 and still
1633		! have salt conserved with SIS2 ice. So for SIS2, we can run with
1634		! BOUND_SALINITY=False in MOM.F90.
1635	0	if (associated(tv%S) .and. associated(tv%salt_deficit)) &
1636	0	call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp)
1637
1638	0	if (CS%diabatic_diff_tendency_diag) then
1639	0	do k=1,nz ; do j=js,je ; do i=is,ie
1640	0	temp_diag(i,j,k) = tv%T(i,j,k)
1641	0	saln_diag(i,j,k) = tv%S(i,j,k)
1642		enddo ; enddo ; enddo
1643		endif
1644
1645		! Find the vertical distances across layers, which may have been modified by the net surface flux
1646	0	call thickness_to_dz(h, tv, dz, G, GV, US)
1647
1648		! set ent_t=dtKd_heat/h_int and est_s=dtKd_salt/h_int on interfaces for use in the tridiagonal solver.
1649	0	do j=js,je ; do i=is,ie
1650	0	ent_t(i,j,1) = 0. ; ent_t(i,j,nz+1) = 0.
1651	0	ent_s(i,j,1) = 0. ; ent_s(i,j,nz+1) = 0.
1652		enddo ; enddo
1653
1654		!$OMP parallel do default(shared) private(I_dzval)
1655	0	do K=2,nz ; do j=js,je ; do i=is,ie
1656	0	I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))
1657	0	ent_t(i,j,K) = dt * I_dzval * Kd_heat(i,j,k)
1658	0	ent_s(i,j,K) = dt * I_dzval * Kd_salt(i,j,k)
1659		enddo ; enddo ; enddo
1660	0	if (showCallTree) call callTree_waypoint("done setting ent_t and ent_t from Kd_heat and " //&
1661	0	"Kd_salt (diabatic_ALE)")
1662
1663		! Changes T and S via the tridiagonal solver; no change to h
1664	0	call tracer_vertdiff_Eulerian(h, ent_t, dt, tv%T, G, GV)
1665	0	call tracer_vertdiff_Eulerian(h, ent_s, dt, tv%S, G, GV)
1666
1667		! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below
1668	0	if (CS%diabatic_diff_tendency_diag) then
1669	0	call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, G, GV, US, CS)
1670		endif
1671	0	call cpu_clock_end(id_clock_tridiag)
1672
1673	0	if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)")
1674
1675		endif ! endif corresponding to if (associated(tv%T))
1676
1677	0	if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G, GV, US)
1678
1679	0	if (CS%debug) then
1680	0	call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, US, haloshift=0)
1681	0	call MOM_thermovar_chksum("after mixed layer ", tv, G, US)
1682		endif
1683
1684		! Whenever thickness changes let the diag manager know, as the
1685		! target grids for vertical remapping may need to be regenerated.
1686	0	call diag_update_remap_grids(CS%diag)
1687
1688		! Set diffusivities for VBF diagnostics if enabled
1689	0	if (CS%use_energetic_PBL .and. associated(CS%VBF%Kd_ePBL)) CS%VBF%Kd_ePBL(:,:,:) = Kd_ePBL(:,:,:)
1690	0	if (associated(CS%VBF%Kd_salt)) CS%VBF%Kd_temp(:,:,:) = Kd_heat(:,:,:)
1691	0	if (associated(CS%VBF%Kd_temp)) CS%VBF%Kd_salt(:,:,:) = Kd_salt(:,:,:)
1692
1693		! Diagnose the diapycnal diffusivities and other related quantities.
1694	0	if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag)
1695	0	if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag)
1696	0	if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag)
1697	0	if (CS%id_Kd_int > 0) then
1698	0	if (CS%double_diffuse .or. CS%useKPP) then
1699		! Using this as a work array might cause confusion.
1700	0	do K=1,nz ; do j=js,je ; do i=is,ie
1701	0	Kd_heat(i,j,k) = min(Kd_heat(i,j,k), Kd_salt(i,j,k))
1702		enddo ; enddo ; enddo
1703		endif
1704	0	call post_data(CS%id_Kd_int, Kd_heat, CS%diag)
1705		endif
1706
1707	0	if (CS%id_ea_t > 0) call post_data(CS%id_ea_t, ent_t(:,:,1:nz), CS%diag)
1708	0	if (CS%id_eb_t > 0) call post_data(CS%id_eb_t, ent_t(:,:,2:nz+1), CS%diag)
1709	0	if (CS%id_ea_s > 0) call post_data(CS%id_ea_s, ent_s(:,:,1:nz), CS%diag)
1710	0	if (CS%id_eb_s > 0) call post_data(CS%id_eb_s, ent_s(:,:,2:nz+1), CS%diag)
1711
1712	0	Idt = 1.0 / dt
1713	0	if (CS%id_Tdif > 0) then
1714	0	do j=js,je ; do i=is,ie
1715	0	Tdif_flx(i,j,1) = 0.0 ; Tdif_flx(i,j,nz+1) = 0.0
1716		enddo ; enddo
1717		!$OMP parallel do default(shared)
1718	0	do K=2,nz ; do j=js,je ; do i=is,ie
1719	0	Tdif_flx(i,j,K) = (Idt * ent_t(i,j,K)) * (tv%T(i,j,k-1) - tv%T(i,j,k))
1720		enddo ; enddo ; enddo
1721	0	if (CS%id_Tdif > 0) call post_data(CS%id_Tdif, Tdif_flx, CS%diag)
1722		endif
1723	0	if (CS%id_Sdif > 0) then
1724	0	do j=js,je ; do i=is,ie
1725	0	Sdif_flx(i,j,1) = 0.0 ; Sdif_flx(i,j,nz+1) = 0.0
1726		enddo ; enddo
1727		!$OMP parallel do default(shared)
1728	0	do K=2,nz ; do j=js,je ; do i=is,ie
1729	0	Sdif_flx(i,j,K) = (Idt * ent_s(i,j,K)) * (tv%S(i,j,k-1) - tv%S(i,j,k))
1730		enddo ; enddo ; enddo
1731	0	if (CS%id_Sdif > 0) call post_data(CS%id_Sdif, Sdif_flx, CS%diag)
1732		endif
1733
1734	0	if (CS%Use_KdWork_diag .or. CS%Use_N2_diag) then
1735	0	N2_salt(:,:,:) = 0.0
1736	0	N2_temp(:,:,:) = 0.0
1737		!Compute N2 and don't mask negatives here
1738	0	EOSdom(:) = EOS_domain(G%HI)
1739	0	if (nonBous) then
1740		!$OMP parallel do default(shared)
1741	0	do j=js,je
1742	0	if (associated(tv%p_surf)) then
1743	0	do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo
1744		else
1745	0	do i=is,ie ; p_i(i) = 0.0 ; enddo
1746		endif
1747	0	do K=2,nz
1748	0	do i=is,ie
1749	0	p_i(i) = p_i(i) + H_to_pres * h(i,j,k-1)
1750		enddo
1751	0	T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))
1752	0	S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))
1753	0	call calculate_specific_vol_derivs(T_i, S_i, p_i, dSpV_dT, dSpV_dS, tv%eqn_of_state, EOSdom)
1754	0	do i=is,ie
1755	0	I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))
1756	0	N2_salt(i,j,K) = (tv%S(i,j,k-1) - tv%S(i,j,k)) * (dSpv_dS(i) * (alt_H_to_pres * I_dzval))
1757	0	N2_temp(i,j,K) = (tv%T(i,j,k-1) - tv%T(i,j,k)) * (dSpV_dT(i) * (alt_H_to_pres * I_dzval))
1758		enddo
1759		enddo
1760		enddo
1761		else
1762		!$OMP parallel do default(shared)
1763	0	do j=js,je
1764	0	if (associated(tv%p_surf)) then
1765	0	do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo
1766		else
1767	0	do i=is,ie ; p_i(i) = 0.0 ; enddo
1768		endif
1769	0	do K=2,nz
1770	0	do i=is,ie
1771	0	p_i(i) = p_i(i) + H_to_pres* h(i,j,k-1)
1772		enddo
1773	0	T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))
1774	0	S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))
1775	0	call calculate_density_derivs(T_i, S_i, p_i, dRhodT, dRhodS, tv%eqn_of_state, EOSdom)
1776	0	do i=is,ie
1777	0	I_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
1778	0	N2_salt(i,j,K) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (dRhodS(i) * (g_rho0 * I_h))
1779	0	N2_temp(i,j,K) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (dRhodT(i) * (g_rho0 * I_h))
1780		enddo
1781		enddo
1782		enddo
1783		endif
1784	0	if (CS%id_N2_dd>0) call post_data(CS%id_N2_dd, N2_salt(:,:,:)+N2_temp(:,:,:), CS%diag)
1785	0	if (CS%id_N2_salt_dd>0) call post_data(CS%id_N2_salt_dd, N2_salt, CS%diag)
1786	0	if (CS%id_N2_temp_dd>0) call post_data(CS%id_N2_temp_dd, N2_temp, CS%diag)
1787
1788	0	if (CS%Use_KdWork_diag) then
1789	0	call KdWork_diagnostics(G,GV,US,CS%diag,CS%VBF,N2_salt,N2_temp,dz)
1790		endif
1791
1792	0	call deallocate_VBF_CS(CS%VBF)
1793
1794		endif
1795
1796		! mixing of passive tracers from massless boundary layers to interior
1797	0	call cpu_clock_begin(id_clock_tracers)
1798
1799	0	if (CS%mix_boundary_tracer_ALE) then
1800	0	Tr_ea_BBL = sqrt(dt * CS%Kd_BBL_tr)
1801		!$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
1802	0	do j=js,je
1803	0	do i=is,ie
1804	0	htot(i) = 0.0
1805	0	in_boundary(i) = (G%mask2dT(i,j) > 0.0)
1806		enddo
1807	0	do k=nz,2,-1 ; do i=is,ie
1808	0	if (in_boundary(i)) then
1809	0	htot(i) = htot(i) + h(i,j,k)
1810		! If diapycnal mixing has been suppressed because this is a massless layer near the
1811		! bottom, add some mixing of tracers between these layers. This flux is based on the
1812		! harmonic mean of the two thicknesses, following what is done in layered mode. Kd_min_tr
1813		! should be much less than the values in Kd_salt, perhaps a molecular diffusivity.
1814		add_ent = (dt * CS%Kd_min_tr) * &
1815		((dz(i,j,k-1)+dz(i,j,k) + dz_neglect) / (dz(i,j,k-1)*dz(i,j,k) + dz_neglect2)) - &
1816	0	ent_s(i,j,K)
1817	0	if (htot(i) < Tr_ea_BBL) then
1818	0	add_ent = max(0.0, add_ent, (Tr_ea_BBL - htot(i)) - ent_s(i,j,K))
1819	0	elseif (add_ent < 0.0) then
1820	0	add_ent = 0.0 ; in_boundary(i) = .false.
1821		endif
1822
1823	0	ent_s(i,j,K) = ent_s(i,j,K) + add_ent
1824		endif
1825		enddo ; enddo
1826		enddo
1827		endif ! (CS%mix_boundary_tracer_ALE)
1828
1829		! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
1830		call call_tracer_column_fns(h_orig, h, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), fluxes, BLD, dt, &
1831		G, GV, US, tv, CS%optics, CS%tracer_flow_CSp, CS%debug, &
1832		KPP_CSp=CS%KPP_CSp, &
1833		nonLocalTrans=KPP_NLTscalar, &
1834		evap_CFL_limit=CS%evap_CFL_limit, &
1835	0	minimum_forcing_depth=CS%minimum_forcing_depth, h_BL=visc%h_ML)
1836
1837	0	call cpu_clock_end(id_clock_tracers)
1838
1839		! Apply ALE sponge
1840	0	if (CS%use_sponge .and. associated(CS%ALE_sponge_CSp)) then
1841	0	call cpu_clock_begin(id_clock_sponge)
1842	0	call apply_ALE_sponge(h, tv, dt, G, GV, US, CS%ALE_sponge_CSp, CS%Time)
1843	0	call cpu_clock_end(id_clock_sponge)
1844	0	if (CS%debug) then
1845	0	call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, US, haloshift=0)
1846	0	call MOM_thermovar_chksum("apply_sponge ", tv, G, US)
1847		endif
1848		endif ! CS%use_sponge
1849
1850		! Apply data assimilation incremental update -oda_incupd-
1851	0	if (CS%use_oda_incupd .and. associated(CS%oda_incupd_CSp)) then
1852	0	call MOM_mesg("Starting ODA_INCUPD ", 5)
1853	0	call cpu_clock_begin(id_clock_oda_incupd)
1854	0	call apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS%oda_incupd_CSp)
1855	0	call cpu_clock_end(id_clock_oda_incupd)
1856	0	if (CS%debug) then
1857	0	call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0)
1858	0	call MOM_thermovar_chksum("apply_oda_incupd ", tv, G, US)
1859		endif
1860		endif ! CS%use_oda_incupd
1861
1862
1863	0	call cpu_clock_begin(id_clock_pass)
1864		! visc%Kv_slow is not in the group pass because it has larger vertical extent.
1865	0	if (associated(visc%Kv_slow)) &
1866	0	call pass_var(visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1)
1867	0	call cpu_clock_end(id_clock_pass)
1868
1869	0	call disable_averaging(CS%diag)
1870
1871	0	if (showCallTree) call callTree_leave("diabatic_ALE()")
1872
1873	0	end subroutine diabatic_ALE
1874
1875		!> Imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers
1876		!! using the original MOM6 algorithms.
1877	0	subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
1878		G, GV, US, CS, Waves)
1879		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1880		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1881		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1882		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1883		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
1884		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2]
1885		type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
1886		!! unused have NULL ptrs
1887		real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Active mixed layer depth [Z ~> m]
1888		type(forcing), intent(inout) :: fluxes !< points to forcing fields
1889		!! unused fields have NULL ptrs
1890		type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities,
1891		!! BBL properties and related fields
1892		type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum
1893		!! equations, to enable the later derived
1894		!! diagnostics, like energy budgets
1895		type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations
1896		real, intent(in) :: dt !< time increment [T ~> s]
1897		type(time_type), intent(in) :: Time_end !< Time at the end of the interval
1898		type(diabatic_CS), pointer :: CS !< module control structure
1899		type(Wave_parameters_CS), pointer :: Waves !< Surface gravity waves
1900
1901		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
1902	0	ea, & ! amount of fluid entrained from the layer above within
1903		! one time step [H ~> m or kg m-2]
1904	0	eb, & ! amount of fluid entrained from the layer below within
1905		! one time step [H ~> m or kg m-2]
1906	0	Kd_lay, & ! diapycnal diffusivity of layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1907	0	h_orig, & ! initial layer thicknesses [H ~> m or kg m-2]
1908	0	dz, & ! The vertical distance between interfaces around a layer [Z ~> m]
1909	0	hold, & ! layer thickness before diapycnal entrainment, and later the initial
1910		! layer thicknesses (if a mixed layer is used) [H ~> m or kg m-2]
1911	0	dz_old, & ! The initial vertical distance between interfaces around a layer
1912		! or the distance before entrainment [Z ~> m]
1913	0	u_h, & ! Zonal velocities at thickness points after entrainment [L T-1 ~> m s-1]
1914	0	v_h, & ! Meridional velocities at thickness points after entrainment [L T-1 ~> m s-1]
1915	0	temp_diag, & ! Diagnostic array of previous temperatures [C ~> degC]
1916	0	saln_diag ! Diagnostic array of previous salinity [S ~> ppt]
1917		real, dimension(SZI_(G),SZJ_(G)) :: &
1918	0	h_MLD, & ! Active mixed layer thickness [H ~> m or kg m-2].
1919	0	U_star, & ! The friction velocity [Z T-1 ~> m s-1].
1920	0	KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1921	0	KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1922	0	Rcv_ml ! Coordinate density of mixed layer [R ~> kg m-3], used for applying sponges
1923
1924		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &
1925		! These are targets so that the space can be shared with eaml & ebml.
1926	0	eatr, & ! The equivalent of ea for tracers, which differs from ea in that it tends to
1927		! homogenize tracers in massless layers near the boundaries [H ~> m or kg m-2]
1928	0	ebtr ! The equivalent of eb for tracers, which differs from eb in that it tends to
1929		! homogenize tracers in massless layers near the boundaries [H ~> m or kg m-2]
1930
1931		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
1932	0	Kd_int, & ! diapycnal diffusivity of interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1933	0	Kd_heat, & ! diapycnal diffusivity of heat [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1934	0	Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1935	0	Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to
1936		! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1937	0	Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to
1938		! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1939	0	KPP_NLTheat, & ! KPP non-local transport for heat [nondim]
1940	0	KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim]
1941	0	KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]
1942	0	Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1943	0	Tadv_flx, & ! advective diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1944	0	Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1945	0	Sadv_flx ! advective diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1946
1947		! The following 3 variables are only used with a bulk mixed layer.
1948		real, pointer, dimension(:,:,:) :: &
1949	0	eaml, & ! The equivalent of ea due to mixed layer processes [H ~> m or kg m-2].
1950	0	ebml ! The equivalent of eb due to mixed layer processes [H ~> m or kg m-2].
1951		! eaml and ebml are pointers to eatr and ebtr so as to reuse the memory as
1952		! the arrays are not needed at the same time.
1953
1954	0	integer :: kb(SZI_(G),SZJ_(G)) ! index of the lightest layer denser
1955		! than the buffer layer [nondim]
1956
1957	0	real :: p_ref_cv(SZI_(G)) ! Reference pressure for the potential density that defines the
1958		! coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].
1959
1960	0	logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below,
1961		! where massive is defined as sufficiently thick that
1962		! the no-flux boundary conditions have not restricted
1963		! the entrainment - usually sqrt(Kd*dt).
1964
1965		real :: h_neglect ! A thickness that is so small it is usually lost
1966		! in roundoff and can be neglected [H ~> m or kg m-2]
1967		real :: dz_neglect ! A vertical distance that is so small it is usually lost
1968		! in roundoff and can be neglected [Z ~> m]
1969		real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2]
1970		real :: net_ent ! The net of ea-eb at an interface [H ~> m or kg m-2]
1971		real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2]
1972		real :: eaval ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
1973		real :: hval ! hval is 2*h at velocity grid points [H ~> m or kg m-2]
1974		real :: h_tr ! h_tr is h at tracer points with a tiny thickness
1975		! added to ensure positive definiteness [H ~> m or kg m-2]
1976		real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is
1977		! coupled to the bottom within a timestep [H ~> m or kg m-2]
1978
1979	0	real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].
1980	0	real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]
1981		real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2]
1982	0	real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim]
1983	0	real :: c1(SZIB_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim]
1984
1985		real :: dt_mix ! The amount of time over which to apply mixing [T ~> s]
1986		real :: Idt ! The inverse time step [T-1 ~> s-1]
1987
1988		integer :: dir_flag ! An integer encoding the directions in which to do halo updates.
1989		logical :: showCallTree ! If true, show the call tree
1990		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1991		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, halo
1992
1993	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1994	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1995	0	nkmb = GV%nk_rho_varies
1996	0	h_neglect = GV%H_subroundoff
1997	0	dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect
1998	0	Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0
1999
2000	0	showCallTree = callTree_showQuery()
2001	0	if (showCallTree) call callTree_enter("layered_diabatic(), MOM_diabatic_driver.F90")
2002
2003		! set equivalence between the same bits of memory for these arrays
2004	0	eaml => eatr ; ebml => ebtr
2005
2006		! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
2007	0	call enable_averages(dt, Time_end, CS%diag)
2008
2009	0	if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then
2010	0	halo = CS%halo_TS_diff
2011		!$OMP parallel do default(shared)
2012	0	do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo
2013	0	h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0
2014		enddo ; enddo ; enddo
2015		endif
2016
2017	0	if (CS%use_geothermal) then
2018	0	call cpu_clock_begin(id_clock_geothermal)
2019	0	call geothermal_entraining(h, tv, dt, eaml, ebml, G, GV, US, CS%geothermal, halo=CS%halo_TS_diff)
2020	0	call cpu_clock_end(id_clock_geothermal)
2021	0	if (showCallTree) call callTree_waypoint("geothermal (diabatic)")
2022	0	if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G, GV, US)
2023		endif
2024
2025		! Whenever thickness changes let the diag manager know, target grids
2026		! for vertical remapping may need to be regenerated.
2027	0	call diag_update_remap_grids(CS%diag)
2028
2029		! Set_pen_shortwave estimates the optical properties of the water column.
2030		! It will need to be modified later to include information about the
2031		! biological properties and layer thicknesses.
2032	0	if (associated(CS%optics)) &
2033	0	call set_pen_shortwave(CS%optics, fluxes, G, GV, US, CS%diabatic_aux_CSp, CS%opacity, CS%tracer_flow_CSp)
2034
2035	0	if (CS%use_bulkmixedlayer) then
2036	0	if (CS%debug) call MOM_forcing_chksum("Before mixedlayer", fluxes, G, US, haloshift=0)
2037
2038	0	if (CS%ML_mix_first > 0.0) then
2039		! This subroutine
2040		! (1) Cools the mixed layer.
2041		! (2) Performs convective adjustment by mixed layer entrainment.
2042		! (3) Heats the mixed layer and causes it to detrain to
2043		! Monin-Obukhov depth or minimum mixed layer depth.
2044		! (4) Uses any remaining TKE to drive mixed layer entrainment.
2045		! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate)
2046	0	call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
2047
2048	0	call cpu_clock_begin(id_clock_mixedlayer)
2049	0	if (CS%ML_mix_first < 1.0) then
2050		! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???)
2051		call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*CS%ML_mix_first, &
2052		eaml, ebml, G, GV, US, CS%bulkmixedlayer, CS%optics, &
2053	0	BLD, h_MLD, CS%aggregate_FW_forcing, dt, last_call=.false.)
2054		else
2055		! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???)
2056		call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, &
2057		G, GV, US, CS%bulkmixedlayer, CS%optics, &
2058	0	BLD, h_MLD, CS%aggregate_FW_forcing, dt, last_call=.true.)
2059	0	if (associated(visc%h_ML)) visc%h_ML(:,:) = h_MLD(:,:)
2060	0	if (associated(visc%MLD)) visc%MLD(:,:) = BLD(:,:)
2061		endif
2062
2063		! Keep salinity from falling below a small but positive threshold.
2064		! This constraint is needed for SIS1 ice model, which can extract
2065		! more salt than is present in the ocean. SIS2 does not suffer
2066		! from this limitation, in which case we can let salinity=0 and still
2067		! have salt conserved with SIS2 ice. So for SIS2, we can run with
2068		! BOUND_SALINITY=False in MOM.F90.
2069	0	if (associated(tv%S) .and. associated(tv%salt_deficit)) &
2070	0	call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp)
2071	0	call cpu_clock_end(id_clock_mixedlayer)
2072	0	if (CS%debug) then
2073	0	call MOM_state_chksum("After mixedlayer ", u, v, h, G, GV, US, haloshift=0)
2074	0	call MOM_forcing_chksum("After mixedlayer", fluxes, G, US, haloshift=0)
2075		endif
2076	0	if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)")
2077	0	if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G, GV, US)
2078		endif
2079		endif
2080
2081	0	if (CS%debug) &
2082	0	call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, US, haloshift=0)
2083	0	if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then
2084	0	if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then
2085	0	call find_uv_at_h(u, v, h_orig, u_h, v_h, G, GV, US, eaml, ebml)
2086	0	if (CS%debug) then
2087	0	call hchksum(eaml, "after find_uv_at_h eaml", G%HI, unscale=GV%H_to_MKS)
2088	0	call hchksum(ebml, "after find_uv_at_h ebml", G%HI, unscale=GV%H_to_MKS)
2089		endif
2090		else
2091	0	call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
2092		endif
2093	0	if (showCallTree) call callTree_waypoint("done with find_uv_at_h (diabatic)")
2094		endif
2095
2096	0	call cpu_clock_begin(id_clock_set_diffusivity)
2097		! Sets: Kd_lay, Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb
2098		! Also changes: visc%Kd_shear and visc%Kv_shear
2099	0	if ((CS%halo_TS_diff > 0) .and. (CS%ML_mix_first > 0.0)) then
2100	0	if (associated(tv%T)) call pass_var(tv%T, G%Domain, halo=CS%halo_TS_diff, complete=.false.)
2101	0	if (associated(tv%S)) call pass_var(tv%S, G%Domain, halo=CS%halo_TS_diff, complete=.false.)
2102	0	call pass_var(h, G%domain, halo=CS%halo_TS_diff, complete=.true.)
2103		endif
2104
2105		! Update derived thermodynamic quantities.
2106	0	if ((CS%ML_mix_first > 0.0) .and. allocated(tv%SpV_avg)) then
2107	0	call calc_derived_thermo(tv, h, G, GV, US, halo=CS%halo_TS_diff)
2108		endif
2109
2110	0	if (CS%debug) &
2111	0	call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff)
2112	0	if (CS%double_diffuse) then
2113		call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, &
2114	0	CS%set_diff_CSp, CS%VBF, Kd_lay=Kd_lay, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S)
2115		else
2116		call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, &
2117	0	CS%set_diff_CSp, CS%VBF, Kd_lay=Kd_lay)
2118		endif
2119	0	call cpu_clock_end(id_clock_set_diffusivity)
2120	0	if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)")
2121
2122	0	if (CS%debug) then
2123	0	call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, US, haloshift=0)
2124	0	call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, US, haloshift=0)
2125	0	call MOM_thermovar_chksum("after set_diffusivity ", tv, G, US)
2126	0	call hchksum(Kd_lay, "after set_diffusivity Kd_lay", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
2127	0	call hchksum(Kd_Int, "after set_diffusivity Kd_Int", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
2128		endif
2129
2130
2131	0	if (CS%useKPP) then
2132	0	call cpu_clock_begin(id_clock_kpp)
2133
2134		! NOTE: The following do not require initialization, but their checksums do
2135		! require initialization, and past versions were initialized to zero.
2136	0	KPP_NLTheat(:,:,:) = 0.
2137	0	KPP_NLTscalar(:,:,:) = 0.
2138	0	KPP_buoy_flux(:,:,:) = 0.
2139	0	KPP_temp_flux(:,:) = 0.
2140	0	KPP_salt_flux(:,:) = 0.
2141
2142		! KPP needs the surface buoyancy flux but does not update state variables.
2143		! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
2144		! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux
2145		! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second)
2146		! unlike other instances where the fluxes are integrated in time over a time-step.
2147		call calculateBuoyancyFlux2d(G, GV, US, fluxes, CS%optics, h, tv%T, tv%S, tv, &
2148	0	KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux)
2149		! The KPP scheme calculates boundary layer diffusivities and non-local transport.
2150
2151		! Set diffusivities for heat and salt separately
2152
2153	0	if (CS%double_diffuse) then
2154		! Add contribution from double diffusion
2155		!$OMP parallel do default(shared)
2156	0	do K=1,nz+1 ; do j=js,je ; do i=is,ie
2157	0	Kd_salt(i,j,K) = Kd_int(i,j,K) + Kd_extra_S(i,j,K)
2158	0	Kd_heat(i,j,K) = Kd_int(i,j,K) + Kd_extra_T(i,j,K)
2159		enddo ; enddo ; enddo
2160		else
2161		!$OMP parallel do default(shared)
2162	0	do K=1,nz+1 ; do j=js,je ; do i=is,ie
2163	0	Kd_salt(i,j,K) = Kd_int(i,j,K)
2164	0	Kd_heat(i,j,K) = Kd_int(i,j,K)
2165		enddo ; enddo ; enddo
2166		endif
2167
2168		! Determine the friction velocity, perhaps using the evovling surface density.
2169	0	call find_ustar(fluxes, tv, U_star, G, GV, US)
2170
2171	0	if ( associated(fluxes%lamult) ) then
2172		call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, &
2173	0	U_star, KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult)
2174
2175		call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, &
2176	0	Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult)
2177		else
2178		call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, &
2179	0	U_star, KPP_buoy_flux, Waves=Waves)
2180
2181		call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, &
2182	0	Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves)
2183		endif
2184
2185	0	call KPP_get_BLD(CS%KPP_CSp, BLD(:,:), G, US)
2186		! If visc%MLD or visc%h_ML exist, copy KPP's BLD into them with appropriate conversions.
2187	0	if (associated(visc%h_ML)) call convert_MLD_to_ML_thickness(BLD, h, visc%h_ML, tv, G, GV)
2188	0	if (associated(visc%MLD)) visc%MLD(:,:) = BLD(:,:)
2189	0	if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = KPP_buoy_flux(:,:,1)
2190
2191	0	if (.not. CS%KPPisPassive) then
2192		!$OMP parallel do default(shared)
2193	0	do k=1,nz+1 ; do j=js,je ; do i=is,ie
2194	0	Kd_int(i,j,K) = min( Kd_salt(i,j,k), Kd_heat(i,j,k) )
2195		enddo ; enddo ; enddo
2196	0	if (CS%double_diffuse) then
2197		!$OMP parallel do default(shared)
2198	0	do k=1,nz+1 ; do j=js,je ; do i=is,ie
2199	0	Kd_extra_S(i,j,k) = (Kd_salt(i,j,k) - Kd_int(i,j,K))
2200	0	Kd_extra_T(i,j,k) = (Kd_heat(i,j,k) - Kd_int(i,j,K))
2201		enddo ; enddo ; enddo
2202		endif
2203		endif ! not passive
2204
2205	0	call cpu_clock_end(id_clock_kpp)
2206	0	if (showCallTree) call callTree_waypoint("done with KPP_calculate (diabatic)")
2207	0	if (CS%debug) then
2208	0	call MOM_state_chksum("after KPP", u, v, h, G, GV, US, haloshift=0)
2209	0	call MOM_forcing_chksum("after KPP", fluxes, G, US, haloshift=0)
2210	0	call MOM_thermovar_chksum("after KPP", tv, G, US)
2211	0	call hchksum(Kd_lay, "after KPP Kd_lay", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
2212	0	call hchksum(Kd_Int, "after KPP Kd_Int", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
2213		endif
2214		endif ! endif for KPP
2215
2216		! Add vertical diff./visc. due to convection (computed via CVMix)
2217	0	if (CS%use_CVMix_conv) then
2218	0	call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv, BLD, Kd_int, visc%Kv_shear)
2219		endif
2220
2221	0	if (CS%useKPP) then
2222	0	call cpu_clock_begin(id_clock_kpp)
2223	0	if (CS%debug) then
2224		call hchksum(KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, &
2225	0	unscale=US%C_to_degCGV%H_to_mUS%s_to_T)
2226		call hchksum(KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, &
2227	0	unscale=US%S_to_pptGV%H_to_mUS%s_to_T)
2228	0	call hchksum(KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0)
2229	0	call hchksum(KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0)
2230		endif
2231		! Apply non-local transport of heat and salt
2232		! Changes: tv%T, tv%S
2233		call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, KPP_NLTheat, KPP_temp_flux, &
2234	0	dt, tv%tr_T, tv%T, tv%C_p)
2235		call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, KPP_NLTscalar, KPP_salt_flux, &
2236	0	dt, tv%tr_S, tv%S)
2237	0	call cpu_clock_end(id_clock_kpp)
2238	0	if (showCallTree) call callTree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
2239	0	if (CS%debugConservation) call MOM_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, G, GV, US)
2240
2241	0	if (CS%debug) then
2242	0	call MOM_state_chksum("after KPP_applyNLT ", u, v, h, G, GV, US, haloshift=0)
2243	0	call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, US, haloshift=0)
2244	0	call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G, US)
2245		endif
2246		endif ! endif for KPP
2247
2248		! Differential diffusion done here.
2249		! Changes: tv%T, tv%S
2250	0	if (CS%double_diffuse .and. associated(tv%T)) then
2251
2252	0	call cpu_clock_begin(id_clock_differential_diff)
2253	0	call differential_diffuse_T_S(h, tv%T, tv%S, Kd_extra_T, Kd_extra_S, tv, dt, G, GV)
2254	0	call cpu_clock_end(id_clock_differential_diff)
2255	0	if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)")
2256	0	if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G, GV, US)
2257
2258		! increment heat and salt diffusivity.
2259		! CS%useKPP==.true. already has extra_T and extra_S included
2260	0	if (.not. CS%useKPP) then
2261		!$OMP parallel do default(shared)
2262	0	do K=2,nz ; do j=js,je ; do i=is,ie
2263	0	Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_extra_T(i,j,K)
2264	0	Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_extra_S(i,j,K)
2265		enddo ; enddo ; enddo
2266		endif
2267
2268		endif
2269
2270		! Calculate layer entrainments and detrainments from diffusivities and differences between
2271		! layer and target densities (i.e. do remapping as well as diffusion).
2272	0	call cpu_clock_begin(id_clock_entrain)
2273		! Calculate appropriately limited diapycnal mass fluxes to account
2274		! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb
2275		call Entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS%entrain_diffusive, &
2276	0	ea, eb, kb, Kd_lay=Kd_lay, Kd_int=Kd_int)
2277	0	call cpu_clock_end(id_clock_entrain)
2278	0	if (showCallTree) call callTree_waypoint("done with Entrainment_diffusive (diabatic)")
2279
2280	0	if (CS%debug) then
2281	0	call MOM_forcing_chksum("after calc_entrain ", fluxes, G, US, haloshift=0)
2282	0	call MOM_thermovar_chksum("after calc_entrain ", tv, G, US)
2283	0	call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, US, haloshift=0)
2284	0	call hchksum(ea, "after calc_entrain ea", G%HI, haloshift=0, unscale=GV%H_to_MKS)
2285	0	call hchksum(eb, "after calc_entrain eb", G%HI, haloshift=0, unscale=GV%H_to_MKS)
2286		endif
2287
2288		! Save fields before boundary forcing is applied for tendency diagnostics
2289	0	if (CS%boundary_forcing_tendency_diag) then
2290	0	do k=1,nz ; do j=js,je ; do i=is,ie
2291	0	temp_diag(i,j,k) = tv%T(i,j,k)
2292	0	saln_diag(i,j,k) = tv%S(i,j,k)
2293		enddo ; enddo ; enddo
2294		endif
2295
2296		! Update h according to divergence of the difference between
2297		! ea and eb. We keep a record of the original h in hold.
2298		! In the following, the checks for negative values are to guard
2299		! against instances where entrainment drives a layer to
2300		! negative thickness. This situation will never happen if
2301		! enough iterations are permitted in Calculate_Entrainment.
2302		! Even if too few iterations are allowed, it is still guarded
2303		! against. In other words the checks are probably unnecessary.
2304		!$OMP parallel do default(shared)
2305	0	do j=js,je
2306	0	do i=is,ie
2307	0	hold(i,j,1) = h(i,j,1)
2308	0	h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2))
2309	0	hold(i,j,nz) = h(i,j,nz)
2310	0	h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1))
2311	0	if (h(i,j,1) <= 0.0) then
2312	0	h(i,j,1) = GV%Angstrom_H
2313		endif
2314	0	if (h(i,j,nz) <= 0.0) then
2315	0	h(i,j,nz) = GV%Angstrom_H
2316		endif
2317		enddo
2318	0	do k=2,nz-1 ; do i=is,ie
2319	0	hold(i,j,k) = h(i,j,k)
2320		h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
2321	0	(eb(i,j,k) - ea(i,j,k+1)))
2322	0	if (h(i,j,k) <= 0.0) then
2323	0	h(i,j,k) = GV%Angstrom_H
2324		endif
2325		enddo ; enddo
2326		enddo
2327		! Checks for negative thickness may have changed layer thicknesses
2328	0	call diag_update_remap_grids(CS%diag)
2329
2330	0	if (CS%debug) then
2331	0	call MOM_state_chksum("after negative check ", u, v, h, G, GV, US, haloshift=0)
2332	0	call MOM_forcing_chksum("after negative check ", fluxes, G, US, haloshift=0)
2333	0	call MOM_thermovar_chksum("after negative check ", tv, G, US)
2334		endif
2335	0	if (showCallTree) call callTree_waypoint("done with h=ea-eb (diabatic)")
2336	0	if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G, GV, US)
2337
2338		! Here, T and S are updated according to ea and eb.
2339		! If using the bulk mixed layer, T and S are also updated
2340		! by surface fluxes (in fluxes%*).
2341		! This is a very long block.
2342	0	if (CS%use_bulkmixedlayer) then
2343
2344	0	if (associated(tv%T)) then
2345	0	call cpu_clock_begin(id_clock_tridiag)
2346		! Temperature and salinity (as state variables) are treated
2347		! differently from other tracers to insure massless layers that
2348		! are lighter than the mixed layer have temperatures and salinities
2349		! that correspond to their prescribed densities.
2350	0	if (CS%massless_match_targets) then
2351		!$OMP parallel do default (shared) private(h_tr,b1,d1,c1,b_denom_1)
2352	0	do j=js,je
2353	0	do i=is,ie
2354	0	h_tr = hold(i,j,1) + h_neglect
2355	0	b1(i) = 1.0 / (h_tr + eb(i,j,1))
2356	0	d1(i) = h_tr * b1(i)
2357	0	tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1))
2358	0	tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1))
2359		enddo
2360	0	do k=2,nkmb ; do i=is,ie
2361	0	c1(i,k) = eb(i,j,k-1) * b1(i)
2362	0	h_tr = hold(i,j,k) + h_neglect
2363	0	b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2364	0	b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2365	0	if (k<nkmb) d1(i) = b_denom_1 * b1(i)
2366	0	tv%T(i,j,k) = b1(i) * (h_trtv%T(i,j,k) + ea(i,j,k)tv%T(i,j,k-1))
2367	0	tv%S(i,j,k) = b1(i) * (h_trtv%S(i,j,k) + ea(i,j,k)tv%S(i,j,k-1))
2368		enddo ; enddo
2369
2370	0	do k=nkmb+1,nz ; do i=is,ie
2371	0	if (k == kb(i,j)) then
2372	0	c1(i,k) = eb(i,j,k-1) * b1(i)
2373		d1(i) = (((eb(i,j,nkmb)-eb(i,j,k-1)) + hold(i,j,nkmb) + h_neglect) + &
2374	0	d1(i)ea(i,j,nkmb)) b1(i)
2375	0	h_tr = hold(i,j,k) + h_neglect
2376	0	b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2377	0	b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2378	0	d1(i) = b_denom_1 * b1(i)
2379	0	tv%T(i,j,k) = b1(i) * (h_trtv%T(i,j,k) + ea(i,j,k)tv%T(i,j,nkmb))
2380	0	tv%S(i,j,k) = b1(i) * (h_trtv%S(i,j,k) + ea(i,j,k)tv%S(i,j,nkmb))
2381	0	elseif (k > kb(i,j)) then
2382	0	c1(i,k) = eb(i,j,k-1) * b1(i)
2383	0	h_tr = hold(i,j,k) + h_neglect
2384	0	b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2385	0	b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2386	0	d1(i) = b_denom_1 * b1(i)
2387	0	tv%T(i,j,k) = b1(i) * (h_trtv%T(i,j,k) + ea(i,j,k)tv%T(i,j,k-1))
2388	0	tv%S(i,j,k) = b1(i) * (h_trtv%S(i,j,k) + ea(i,j,k)tv%S(i,j,k-1))
2389	0	elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j))
2390		! The bottommost buffer layer might entrain all the mass from some
2391		! of the interior layers that are thin and lighter in the coordinate
2392		! density than that buffer layer. The T and S of these newly
2393		! massless interior layers are unchanged.
2394	0	tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k)
2395	0	tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k)
2396		endif
2397		enddo ; enddo
2398
2399	0	do k=nz-1,nkmb,-1 ; do i=is,ie
2400	0	if (k >= kb(i,j)) then
2401	0	tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
2402	0	tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
2403		endif
2404		enddo ; enddo
2405	0	do i=is,ie ; if (kb(i,j) <= nz) then
2406	0	tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j))
2407	0	tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j))
2408		endif ; enddo
2409	0	do k=nkmb-1,1,-1 ; do i=is,ie
2410	0	tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
2411	0	tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
2412		enddo ; enddo
2413		enddo ! end of j loop
2414		else ! .not. massless_match_targets
2415		! This simpler form allows T & S to be too dense for the layers
2416		! between the buffer layers and the interior.
2417		! Changes: T, S
2418	0	if (CS%tracer_tridiag) then
2419	0	call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV)
2420	0	call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV)
2421		else
2422	0	call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
2423		endif
2424		endif ! massless_match_targets
2425	0	call cpu_clock_end(id_clock_tridiag)
2426
2427		endif ! endif for associated(T)
2428	0	if (CS%debugConservation) call MOM_state_stats('BML tridiag', u, v, h, tv%T, tv%S, G, GV, US)
2429
2430	0	if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then
2431		! The mixed layer code has already been called, but there is some needed
2432		! bookkeeping.
2433		!$OMP parallel do default(shared)
2434	0	do k=1,nz ; do j=js,je ; do i=is,ie
2435	0	hold(i,j,k) = h_orig(i,j,k)
2436	0	ea(i,j,k) = ea(i,j,k) + eaml(i,j,k)
2437	0	eb(i,j,k) = eb(i,j,k) + ebml(i,j,k)
2438		enddo ; enddo ; enddo
2439	0	if (CS%debug) then
2440	0	call hchksum(ea, "after ea = ea + eaml", G%HI, haloshift=0, unscale=GV%H_to_MKS)
2441	0	call hchksum(eb, "after eb = eb + ebml", G%HI, haloshift=0, unscale=GV%H_to_MKS)
2442		endif
2443		endif
2444
2445	0	if (CS%ML_mix_first < 1.0) then
2446		! Call the mixed layer code now, perhaps for a second time.
2447		! This subroutine (1) Cools the mixed layer.
2448		! (2) Performs convective adjustment by mixed layer entrainment.
2449		! (3) Heats the mixed layer and causes it to detrain to
2450		! Monin-Obukhov depth or minimum mixed layer depth.
2451		! (4) Uses any remaining TKE to drive mixed layer entrainment.
2452		! (5) Possibly splits the buffer layer into two isopycnal layers.
2453
2454	0	call find_uv_at_h(u, v, hold, u_h, v_h, G, GV, US, ea, eb)
2455	0	if (CS%debug) call MOM_state_chksum("find_uv_at_h1 ", u, v, h, G, GV, US, haloshift=0)
2456
2457	0	dt_mix = min(dt, dt*(1.0 - CS%ML_mix_first))
2458	0	call cpu_clock_begin(id_clock_mixedlayer)
2459		! Changes: h, tv%T, tv%S, ea and eb (G is also inout???)
2460		call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, &
2461		G, GV, US, CS%bulkmixedlayer, CS%optics, &
2462	0	BLD, h_MLD, CS%aggregate_FW_forcing, dt, last_call=.true.)
2463	0	if (associated(visc%h_ML)) visc%h_ML(:,:) = h_MLD(:,:)
2464	0	if (associated(visc%MLD)) visc%MLD(:,:) = BLD(:,:)
2465
2466		! Keep salinity from falling below a small but positive threshold.
2467		! This constraint is needed for SIS1 ice model, which can extract
2468		! more salt than is present in the ocean. SIS2 does not suffer
2469		! from this limitation, in which case we can let salinity=0 and still
2470		! have salt conserved with SIS2 ice. So for SIS2, we can run with
2471		! BOUND_SALINITY=False in MOM.F90.
2472	0	if (associated(tv%S) .and. associated(tv%salt_deficit)) &
2473	0	call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp)
2474
2475	0	call cpu_clock_end(id_clock_mixedlayer)
2476	0	if (showCallTree) call callTree_waypoint("done with 2nd bulkmixedlayer (diabatic)")
2477	0	if (CS%debugConservation) call MOM_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, G, GV, US)
2478		endif
2479
2480		else ! following block for when NOT using BULKMIXEDLAYER
2481
2482		! calculate change in temperature & salinity due to dia-coordinate surface diffusion
2483	0	if (associated(tv%T)) then
2484
2485	0	if (CS%debug) then
2486	0	call hchksum(ea, "before triDiagTS ea ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
2487	0	call hchksum(eb, "before triDiagTS eb ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
2488		endif
2489	0	call cpu_clock_begin(id_clock_tridiag)
2490
2491		! Keep salinity from falling below a small but positive threshold.
2492		! This constraint is needed for SIS1 ice model, which can extract
2493		! more salt than is present in the ocean. SIS2 does not suffer
2494		! from this limitation, in which case we can let salinity=0 and still
2495		! have salt conserved with SIS2 ice. So for SIS2, we can run with
2496		! BOUND_SALINITY=False in MOM.F90.
2497	0	if (associated(tv%S) .and. associated(tv%salt_deficit)) &
2498	0	call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp)
2499
2500	0	if (CS%diabatic_diff_tendency_diag) then
2501	0	do k=1,nz ; do j=js,je ; do i=is,ie
2502	0	temp_diag(i,j,k) = tv%T(i,j,k)
2503	0	saln_diag(i,j,k) = tv%S(i,j,k)
2504		enddo ; enddo ; enddo
2505		endif
2506
2507		! Changes T and S via the tridiagonal solver; no change to h
2508	0	if (CS%tracer_tridiag) then
2509	0	call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV)
2510	0	call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV)
2511		else
2512	0	call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
2513		endif
2514
2515		! diagnose temperature, salinity, heat, and salt tendencies
2516		! Note: hold here refers to the thicknesses from before the dual-entrainment when using
2517		! the bulk mixed layer scheme, so tendencies should be posted on hold.
2518	0	if (CS%diabatic_diff_tendency_diag) then
2519	0	call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, G, GV, US, CS)
2520	0	if (CS%id_diabatic_diff_h > 0) call post_data(CS%id_diabatic_diff_h, hold, CS%diag, alt_h=hold)
2521		endif
2522
2523	0	call cpu_clock_end(id_clock_tridiag)
2524	0	if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)")
2525
2526		endif ! endif corresponding to if (associated(tv%T))
2527	0	if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G, GV, US)
2528
2529		endif ! endif for the BULKMIXEDLAYER block
2530
2531	0	if (CS%debug) then
2532	0	call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, US, haloshift=0)
2533	0	call MOM_thermovar_chksum("after mixed layer ", tv, G, US)
2534	0	call hchksum(ea, "after mixed layer ea", G%HI, unscale=GV%H_to_MKS)
2535	0	call hchksum(eb, "after mixed layer eb", G%HI, unscale=GV%H_to_MKS)
2536		endif
2537
2538	0	call cpu_clock_begin(id_clock_remap)
2539	0	call regularize_layers(h, tv, dt, ea, eb, G, GV, US, CS%regularize_layers)
2540	0	call cpu_clock_end(id_clock_remap)
2541	0	if (showCallTree) call callTree_waypoint("done with regularize_layers (diabatic)")
2542	0	if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, tv%T, tv%S, G, GV, US)
2543
2544		! Whenever thickness changes let the diag manager know, as the
2545		! target grids for vertical remapping may need to be regenerated.
2546	0	if (associated(ADp%du_dt_dia) .or. associated(ADp%dv_dt_dia)) &
2547		! Remapped d[uv]dt_dia require east/north halo updates of h
2548	0	call pass_var(h, G%domain, To_West+To_South+Omit_Corners, halo=1)
2549	0	call diag_update_remap_grids(CS%diag)
2550
2551		! diagnostics
2552	0	Idt = 1.0 / dt
2553	0	if ((CS%id_Tdif > 0) .or. (CS%id_Tadv > 0)) then
2554	0	do j=js,je ; do i=is,ie
2555	0	Tdif_flx(i,j,1) = 0.0 ; Tdif_flx(i,j,nz+1) = 0.0
2556	0	Tadv_flx(i,j,1) = 0.0 ; Tadv_flx(i,j,nz+1) = 0.0
2557		enddo ; enddo
2558		!$OMP parallel do default(shared)
2559	0	do K=2,nz ; do j=js,je ; do i=is,ie
2560		Tdif_flx(i,j,K) = (Idt * 0.5(ea(i,j,k) + eb(i,j,k-1))) &
2561	0	(tv%T(i,j,k-1) - tv%T(i,j,k))
2562		Tadv_flx(i,j,K) = (Idt * (ea(i,j,k) - eb(i,j,k-1))) * &
2563	0	0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
2564		enddo ; enddo ; enddo
2565		endif
2566	0	if ((CS%id_Sdif > 0) .or. (CS%id_Sadv > 0)) then
2567	0	do j=js,je ; do i=is,ie
2568	0	Sdif_flx(i,j,1) = 0.0 ; Sdif_flx(i,j,nz+1) = 0.0
2569	0	Sadv_flx(i,j,1) = 0.0 ; Sadv_flx(i,j,nz+1) = 0.0
2570		enddo ; enddo
2571		!$OMP parallel do default(shared)
2572	0	do K=2,nz ; do j=js,je ; do i=is,ie
2573		Sdif_flx(i,j,K) = (Idt * 0.5(ea(i,j,k) + eb(i,j,k-1))) &
2574	0	(tv%S(i,j,k-1) - tv%S(i,j,k))
2575		Sadv_flx(i,j,K) = (Idt * (ea(i,j,k) - eb(i,j,k-1))) * &
2576	0	0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
2577		enddo ; enddo ; enddo
2578		endif
2579
2580		! mixing of passive tracers from massless boundary layers to interior
2581	0	call cpu_clock_begin(id_clock_tracers)
2582
2583		! Find the vertical distances across layers.
2584	0	if (CS%mix_boundary_tracers .or. CS%double_diffuse) &
2585	0	call thickness_to_dz(h, tv, dz, G, GV, US)
2586	0	if (CS%double_diffuse) &
2587	0	call thickness_to_dz(hold, tv, dz_old, G, GV, US)
2588
2589	0	if (CS%mix_boundary_tracers) then
2590	0	Tr_ea_BBL = sqrt(dt * CS%Kd_BBL_tr)
2591		!$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
2592	0	do j=js,je
2593	0	do i=is,ie
2594	0	ebtr(i,j,nz) = eb(i,j,nz)
2595	0	htot(i) = 0.0
2596	0	in_boundary(i) = (G%mask2dT(i,j) > 0.0)
2597		enddo
2598	0	do k=nz,2,-1 ; do i=is,ie
2599	0	if (in_boundary(i)) then
2600	0	htot(i) = htot(i) + h(i,j,k)
2601		! If diapycnal mixing has been suppressed because this is a massless
2602		! layer near the bottom, add some mixing of tracers between these
2603		! layers. This flux is based on the harmonic mean of the two
2604		! thicknesses, as this corresponds pretty closely (to within
2605		! differences in the density jumps between layers) with what is done
2606		! in the calculation of the fluxes in the first place. Kd_min_tr
2607		! should be much less than the values that have been set in Kd_lay,
2608		! perhaps a molecular diffusivity.
2609		add_ent = (dt * CS%Kd_min_tr) * &
2610		((dz(i,j,k-1) + dz(i,j,k) + dz_neglect) / &
2611		(dz(i,j,k-1)*dz(i,j,k) + dz_neglect2)) - &
2612	0	0.5*(ea(i,j,k) + eb(i,j,k-1))
2613	0	if (htot(i) < Tr_ea_BBL) then
2614		add_ent = max(0.0, add_ent, &
2615	0	(Tr_ea_BBL - htot(i)) - min(ea(i,j,k), eb(i,j,k-1)))
2616	0	elseif (add_ent < 0.0) then
2617	0	add_ent = 0.0 ; in_boundary(i) = .false.
2618		endif
2619
2620	0	ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
2621	0	eatr(i,j,k) = ea(i,j,k) + add_ent
2622		else
2623	0	ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k)
2624		endif
2625	0	if (CS%double_diffuse) then ; if (Kd_extra_S(i,j,K) > 0.0) then
2626		add_ent = (dt * Kd_extra_S(i,j,K)) / &
2627	0	(0.25 * ((dz(i,j,k-1) + dz(i,j,k)) + (dz_old(i,j,k-1) + dz_old(i,j,k))) + dz_neglect)
2628	0	ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
2629	0	eatr(i,j,k) = eatr(i,j,k) + add_ent
2630		endif ; endif
2631		enddo ; enddo
2632	0	do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo
2633
2634		enddo
2635
2636		call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, BLD, dt, G, GV, US, tv, &
2637		CS%optics, CS%tracer_flow_CSp, CS%debug, &
2638	0	KPP_CSp=CS%KPP_CSp, nonLocalTrans=KPP_NLTscalar, h_BL=visc%h_ML)
2639
2640	0	elseif (CS%double_diffuse) then ! extra diffusivity for passive tracers
2641
2642	0	do j=js,je ; do i=is,ie
2643	0	ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1)
2644		enddo ; enddo
2645		!$OMP parallel do default(shared) private(add_ent)
2646	0	do k=nz,2,-1 ; do j=js,je ; do i=is,ie
2647	0	if (Kd_extra_S(i,j,K) > 0.0) then
2648		add_ent = (dt * Kd_extra_S(i,j,K)) / &
2649	0	(0.25 * ((dz(i,j,k-1) + dz(i,j,k)) + (dz_old(i,j,k-1) + dz_old(i,j,k))) + dz_neglect)
2650		else
2651	0	add_ent = 0.0
2652		endif
2653	0	ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
2654	0	eatr(i,j,k) = ea(i,j,k) + add_ent
2655		enddo ; enddo ; enddo
2656
2657		call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, BLD, dt, G, GV, US, tv, &
2658		CS%optics, CS%tracer_flow_CSp, CS%debug, &
2659	0	KPP_CSp=CS%KPP_CSp, nonLocalTrans=KPP_NLTscalar, h_BL=visc%h_ML)
2660
2661		else
2662		call call_tracer_column_fns(hold, h, ea, eb, fluxes, BLD, dt, G, GV, US, tv, &
2663		CS%optics, CS%tracer_flow_CSp, CS%debug, &
2664	0	KPP_CSp=CS%KPP_CSp, nonLocalTrans=KPP_NLTscalar, h_BL=visc%h_ML)
2665
2666		endif ! (CS%mix_boundary_tracers)
2667
2668	0	call cpu_clock_end(id_clock_tracers)
2669
2670		! sponges
2671	0	if (CS%use_sponge) then
2672	0	call cpu_clock_begin(id_clock_sponge)
2673		! Layer mode sponge
2674	0	if (CS%use_bulkmixedlayer .and. associated(tv%eqn_of_state)) then
2675	0	do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo
2676	0	EOSdom(:) = EOS_domain(G%HI)
2677		!$OMP parallel do default(shared)
2678	0	do j=js,je
2679		call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), &
2680	0	tv%eqn_of_state, EOSdom)
2681		enddo
2682	0	call apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS%sponge_CSp, Rcv_ml)
2683		else
2684	0	call apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS%sponge_CSp)
2685		endif
2686	0	call cpu_clock_end(id_clock_sponge)
2687	0	if (CS%debug) then
2688	0	call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, US, haloshift=0)
2689	0	call MOM_thermovar_chksum("apply_sponge ", tv, G, US)
2690		endif
2691		endif ! CS%use_sponge
2692
2693		! Apply data assimilation incremental update -oda_incupd-
2694	0	if (CS%use_oda_incupd .and. associated(CS%oda_incupd_CSp)) then
2695	0	call cpu_clock_begin(id_clock_oda_incupd)
2696	0	call apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS%oda_incupd_CSp)
2697	0	call cpu_clock_end(id_clock_oda_incupd)
2698	0	if (CS%debug) then
2699	0	call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0)
2700	0	call MOM_thermovar_chksum("apply_oda_incupd ", tv, G, US)
2701		endif
2702		endif ! CS%use_oda_incupd
2703
2704
2705
2706		! Save the diapycnal mass fluxes as a diagnostic field.
2707	0	if (associated(CDp%diapyc_vel)) then
2708		!$OMP parallel do default(shared)
2709	0	do j=js,je
2710	0	do K=2,nz ; do i=is,ie
2711	0	CDp%diapyc_vel(i,j,K) = Idt * (ea(i,j,k) - eb(i,j,k-1))
2712		enddo ; enddo
2713	0	do i=is,ie
2714	0	CDp%diapyc_vel(i,j,1) = 0.0
2715	0	CDp%diapyc_vel(i,j,nz+1) = 0.0
2716		enddo
2717		enddo
2718		endif
2719
2720		! For momentum, it is only the net flux that homogenizes within
2721		! the mixed layer. Vertical viscosity that is proportional to the
2722		! mixed layer turbulence is applied elsewhere.
2723	0	if (CS%use_bulkmixedlayer) then
2724	0	if (CS%debug) then
2725	0	call hchksum(ea, "before net flux rearrangement ea", G%HI, unscale=GV%H_to_MKS)
2726	0	call hchksum(eb, "before net flux rearrangement eb", G%HI, unscale=GV%H_to_MKS)
2727		endif
2728		!$OMP parallel do default(shared) private(net_ent)
2729	0	do j=js,je
2730	0	do K=2,GV%nkml ; do i=is,ie
2731	0	net_ent = ea(i,j,k) - eb(i,j,k-1)
2732	0	ea(i,j,k) = max(net_ent, 0.0)
2733	0	eb(i,j,k-1) = max(-net_ent, 0.0)
2734		enddo ; enddo
2735		enddo
2736	0	if (CS%debug) then
2737	0	call hchksum(ea, "after net flux rearrangement ea", G%HI, unscale=GV%H_to_MKS)
2738	0	call hchksum(eb, "after net flux rearrangement eb", G%HI, unscale=GV%H_to_MKS)
2739		endif
2740		endif
2741
2742		! Initialize halo regions of ea, eb, and hold to default values.
2743		!$OMP parallel do default(shared)
2744	0	do k=1,nz
2745	0	do i=is-1,ie+1
2746	0	hold(i,js-1,k) = GV%Angstrom_H ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0
2747	0	hold(i,je+1,k) = GV%Angstrom_H ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0
2748		enddo
2749	0	do j=js,je
2750	0	hold(is-1,j,k) = GV%Angstrom_H ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0
2751	0	hold(ie+1,j,k) = GV%Angstrom_H ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0
2752		enddo
2753		enddo
2754
2755	0	call cpu_clock_begin(id_clock_pass)
2756	0	if (G%symmetric) then ; dir_flag = To_All+Omit_Corners
2757	0	else ; dir_flag = To_West+To_South+Omit_Corners ; endif
2758	0	call create_group_pass(CS%pass_hold_eb_ea, hold, G%Domain, dir_flag, halo=1)
2759	0	call create_group_pass(CS%pass_hold_eb_ea, eb, G%Domain, dir_flag, halo=1)
2760	0	call create_group_pass(CS%pass_hold_eb_ea, ea, G%Domain, dir_flag, halo=1)
2761	0	call do_group_pass(CS%pass_hold_eb_ea, G%Domain)
2762	0	call cpu_clock_end(id_clock_pass)
2763
2764		! Use a tridiagonal solver to determine effect of the diapycnal
2765		! advection on velocity field. It is assumed that water leaves
2766		! or enters the ocean with the surface velocity.
2767	0	if (CS%debug) then
2768	0	call MOM_state_chksum("before u/v tridiag ", u, v, h, G, GV, US, haloshift=0)
2769	0	call hchksum(ea, "before u/v tridiag ea", G%HI, unscale=GV%H_to_MKS)
2770	0	call hchksum(eb, "before u/v tridiag eb", G%HI, unscale=GV%H_to_MKS)
2771	0	call hchksum(hold, "before u/v tridiag hold", G%HI, unscale=GV%H_to_MKS)
2772		endif
2773	0	call cpu_clock_begin(id_clock_tridiag)
2774
2775		!$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)
2776	0	do j=js,je
2777	0	do I=Isq,Ieq
2778	0	if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,1) = u(I,j,1)
2779	0	hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect
2780	0	b1(I) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))
2781	0	d1(I) = hval * b1(I)
2782	0	u(I,j,1) = b1(I) * (hval * u(I,j,1))
2783		enddo
2784	0	do k=2,nz ; do I=Isq,Ieq
2785	0	if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,k) = u(I,j,k)
2786	0	c1(I,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(I)
2787	0	eaval = ea(i,j,k) + ea(i+1,j,k)
2788	0	hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect
2789	0	b1(I) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(I)*eaval))
2790	0	d1(I) = (hval + d1(I)eaval) b1(I)
2791	0	u(I,j,k) = (hvalu(I,j,k) + eavalu(I,j,k-1))*b1(I)
2792		enddo ; enddo
2793	0	do k=nz-1,1,-1 ; do I=Isq,Ieq
2794	0	u(I,j,k) = u(I,j,k) + c1(I,k+1)*u(I,j,k+1)
2795	0	if (associated(ADp%du_dt_dia)) &
2796	0	ADp%du_dt_dia(I,j,k) = (u(I,j,k) - ADp%du_dt_dia(I,j,k)) * Idt
2797		enddo ; enddo
2798	0	if (associated(ADp%du_dt_dia)) then
2799	0	do I=Isq,Ieq
2800	0	ADp%du_dt_dia(I,j,nz) = (u(I,j,nz)-ADp%du_dt_dia(I,j,nz)) * Idt
2801		enddo
2802		endif
2803		enddo
2804	0	if (CS%debug) then
2805	0	call MOM_state_chksum("aft 1st loop tridiag ", u, v, h, G, GV, US, haloshift=0)
2806		endif
2807		!$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)
2808	0	do J=Jsq,Jeq
2809	0	do i=is,ie
2810	0	if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,1) = v(i,J,1)
2811	0	hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect
2812	0	b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))
2813	0	d1(I) = hval * b1(I)
2814	0	v(i,J,1) = b1(i) * (hval * v(i,J,1))
2815		enddo
2816	0	do k=2,nz ; do i=is,ie
2817	0	if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,k) = v(i,J,k)
2818	0	c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)
2819	0	eaval = ea(i,j,k) + ea(i,j+1,k)
2820	0	hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect
2821	0	b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))
2822	0	d1(i) = (hval + d1(i)eaval) b1(i)
2823	0	v(i,J,k) = (hvalv(i,J,k) + eavalv(i,J,k-1))*b1(i)
2824		enddo ; enddo
2825	0	do k=nz-1,1,-1 ; do i=is,ie
2826	0	v(i,J,k) = v(i,J,k) + c1(i,k+1)*v(i,J,k+1)
2827	0	if (associated(ADp%dv_dt_dia)) &
2828	0	ADp%dv_dt_dia(i,J,k) = (v(i,J,k) - ADp%dv_dt_dia(i,J,k)) * Idt
2829		enddo ; enddo
2830	0	if (associated(ADp%dv_dt_dia)) then
2831	0	do i=is,ie
2832	0	ADp%dv_dt_dia(i,J,nz) = (v(i,J,nz)-ADp%dv_dt_dia(i,J,nz)) * Idt
2833		enddo
2834		endif
2835		enddo
2836	0	call cpu_clock_end(id_clock_tridiag)
2837	0	if (CS%debug) then
2838	0	call MOM_state_chksum("after u/v tridiag ", u, v, h, G, GV, US, haloshift=0)
2839		endif
2840
2841		! Diagnose the diapycnal diffusivities and other related quantities.
2842	0	if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_int, CS%diag)
2843	0	if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag)
2844	0	if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag)
2845
2846	0	if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag)
2847	0	if (CS%id_eb > 0) call post_data(CS%id_eb, eb, CS%diag)
2848
2849	0	if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, ADp%du_dt_dia, CS%diag)
2850	0	if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, ADp%dv_dt_dia, CS%diag)
2851	0	if (CS%id_wd > 0) call post_data(CS%id_wd, CDp%diapyc_vel, CS%diag)
2852
2853	0	if (CS%id_Tdif > 0) call post_data(CS%id_Tdif, Tdif_flx, CS%diag)
2854	0	if (CS%id_Tadv > 0) call post_data(CS%id_Tadv, Tadv_flx, CS%diag)
2855	0	if (CS%id_Sdif > 0) call post_data(CS%id_Sdif, Sdif_flx, CS%diag)
2856	0	if (CS%id_Sadv > 0) call post_data(CS%id_Sadv, Sadv_flx, CS%diag)
2857
2858		! Diagnostics for thickness-weighted vertically averaged diapycnal accelerations
2859	0	if (CS%id_hf_dudt_dia_2d > 0) &
2860	0	call post_product_sum_u(CS%id_hf_dudt_dia_2d, ADp%du_dt_dia, ADp%diag_hfrac_u, G, nz, CS%diag)
2861	0	if (CS%id_hf_dvdt_dia_2d > 0) &
2862	0	call post_product_sum_v(CS%id_hf_dvdt_dia_2d, ADp%dv_dt_dia, ADp%diag_hfrac_v, G, nz, CS%diag)
2863
2864	0	call disable_averaging(CS%diag)
2865
2866	0	if (showCallTree) call callTree_leave("layered_diabatic()")
2867
2868	0	end subroutine layered_diabatic
2869
2870		!> Returns pointers or values of members within the diabatic_CS type. For extensibility,
2871		!! each returned argument is an optional argument
2872	25	subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit, minimum_forcing_depth, &
2873		KPP_CSp, energetic_PBL_CSp, diabatic_aux_CSp, diabatic_halo, use_KPP)
2874		type(diabatic_CS), target, intent(in) :: CS !< module control structure
2875		! All output arguments are optional
2876		type(opacity_CS), optional, pointer :: opacity_CSp !< A pointer to be set to the opacity control structure
2877		type(optics_type), optional, pointer :: optics_CSp !< A pointer to be set to the optics control structure
2878		type(KPP_CS), optional, pointer :: KPP_CSp !< A pointer to be set to the KPP CS
2879		type(energetic_PBL_CS), optional, pointer :: energetic_PBL_CSp !< A pointer to be set to the ePBL CS
2880		real, optional, intent( out) :: evap_CFL_limit !<The largest fraction of a layer that can be
2881		!! evaporated in one time-step [nondim].
2882		real, optional, intent( out) :: minimum_forcing_depth !< The smallest depth over which heat
2883		!! and freshwater fluxes are applied [H ~> m or kg m-2].
2884		type(diabatic_aux_CS), optional, pointer :: diabatic_aux_CSp !< A pointer to be set to the diabatic_aux
2885		!! control structure
2886		integer, optional, intent( out) :: diabatic_halo !< The halo size where the diabatic algorithms
2887		!! assume thermodynamics properties are valid.
2888		logical, optional, intent( out) :: use_KPP !< If true, diabatic is using KPP vertical mixing
2889
2890		! Pointers to control structures
2891	25	if (present(opacity_CSp)) opacity_CSp => CS%opacity
2892	25	if (present(optics_CSp)) optics_CSp => CS%optics
2893	25	if (present(KPP_CSp)) KPP_CSp => CS%KPP_CSp
2894	25	if (present(energetic_PBL_CSp) .and. CS%use_energetic_PBL) energetic_PBL_CSp => CS%ePBL
2895	25	if (present(diabatic_aux_CSp)) diabatic_aux_CSp => CS%diabatic_aux_CSp
2896
2897		! Constants within diabatic_CS
2898	25	if (present(evap_CFL_limit)) evap_CFL_limit = CS%evap_CFL_limit
2899	25	if (present(minimum_forcing_depth)) minimum_forcing_depth = CS%minimum_forcing_depth
2900	25	if (present(diabatic_halo)) diabatic_halo = CS%halo_diabatic
2901	25	if (present(use_KPP)) use_KPP = CS%use_KPP
2902	25	end subroutine extract_diabatic_member
2903
2904		!> Routine called for adiabatic physics
2905	0	subroutine adiabatic(h, tv, fluxes, dt, G, GV, US, CS)
2906		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
2907		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2908		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2909		intent(inout) :: h !< thickness [H ~> m or kg m-2]
2910		type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
2911		type(forcing), intent(inout) :: fluxes !< boundary fluxes
2912		real, intent(in) :: dt !< time step [T ~> s]
2913		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2914		type(diabatic_CS), pointer :: CS !< module control structure
2915
2916	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: zeros ! An array of zeros with units of [H ~> m or kg m-2]
2917
2918	0	zeros(:,:,:) = 0.0
2919
2920		call call_tracer_column_fns(h, h, zeros, zeros, fluxes, zeros(:,:,1), dt, G, GV, US, tv, &
2921	0	CS%optics, CS%tracer_flow_CSp, CS%debug)
2922
2923	0	end subroutine adiabatic
2924
2925
2926		!> This routine diagnoses tendencies from application of diabatic diffusion
2927		!! using ALE algorithm. Note that layer thickness is not altered by
2928		!! diabatic diffusion.
2929	0	subroutine diagnose_diabatic_diff_tendency(tv, h, temp_old, saln_old, dt, G, GV, US, CS)
2930		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2931		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2932		type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields
2933		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< thickness [H ~> m or kg m-2]
2934		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: temp_old !< temperature prior to diabatic
2935		!! physics [C ~> degC]
2936		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: saln_old !< salinity prior to diabatic physics [S ~> ppt]
2937		real, intent(in) :: dt !< time step [T ~> s]
2938		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2939		type(diabatic_CS), pointer :: CS !< module control structure
2940
2941		! Local variables
2942	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_3d ! A 3-d work array for diagnostics [various]
2943	0	real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array for diagnostics [various]
2944		real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
2945		integer :: i, j, k, is, ie, js, je, nz
2946		logical :: do_saln_tend ! Calculate salinity-based tendency diagnostics
2947
2948	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2949	0	Idt = 0.0 ; if (dt > 0.0) Idt = 1. / dt
2950	0	work_3d(:,:,:) = 0.0
2951	0	work_2d(:,:) = 0.0
2952
2953
2954		! temperature tendency
2955	0	do k=1,nz ; do j=js,je ; do i=is,ie
2956	0	work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*Idt
2957		enddo ; enddo ; enddo
2958	0	if (CS%id_diabatic_diff_temp_tend > 0) then
2959	0	call post_data(CS%id_diabatic_diff_temp_tend, work_3d, CS%diag, alt_h=h)
2960		endif
2961
2962		! heat tendency
2963	0	if (CS%id_diabatic_diff_heat_tend > 0 .or. CS%id_diabatic_diff_heat_tend_2d > 0) then
2964	0	do k=1,nz ; do j=js,je ; do i=is,ie
2965	0	work_3d(i,j,k) = h(i,j,k)GV%H_to_RZ tv%C_p * work_3d(i,j,k)
2966		enddo ; enddo ; enddo
2967	0	if (CS%id_diabatic_diff_heat_tend > 0) then
2968	0	call post_data(CS%id_diabatic_diff_heat_tend, work_3d, CS%diag, alt_h=h)
2969		endif
2970	0	if (CS%id_diabatic_diff_heat_tend_2d > 0) then
2971	0	work_2d(:,:) = 0.0
2972	0	do k=1,nz ; do j=js,je ; do i=is,ie
2973	0	work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
2974		enddo ; enddo ; enddo
2975	0	call post_data(CS%id_diabatic_diff_heat_tend_2d, work_2d, CS%diag)
2976		endif
2977		endif
2978
2979		! salinity tendency
2980		do_saln_tend = CS%id_diabatic_diff_saln_tend > 0 &
2981		.or. CS%id_diabatic_diff_salt_tend > 0 &
2982	0	.or. CS%id_diabatic_diff_salt_tend_2d > 0
2983
2984	0	if (do_saln_tend) then
2985	0	do k=1,nz ; do j=js,je ; do i=is,ie
2986	0	work_3d(i,j,k) = (tv%S(i,j,k) - saln_old(i,j,k)) * Idt
2987		enddo ; enddo ; enddo
2988
2989	0	if (CS%id_diabatic_diff_saln_tend > 0) &
2990	0	call post_data(CS%id_diabatic_diff_saln_tend, work_3d, CS%diag, alt_h=h)
2991
2992		! salt tendency
2993	0	if (CS%id_diabatic_diff_salt_tend > 0 .or. CS%id_diabatic_diff_salt_tend_2d > 0) then
2994	0	do k=1,nz ; do j=js,je ; do i=is,ie
2995	0	work_3d(i,j,k) = h(i,j,k) * work_3d(i,j,k)
2996		enddo ; enddo ; enddo
2997	0	if (CS%id_diabatic_diff_salt_tend > 0) then
2998	0	call post_data(CS%id_diabatic_diff_salt_tend, work_3d, CS%diag, alt_h=h)
2999		endif
3000	0	if (CS%id_diabatic_diff_salt_tend_2d > 0) then
3001	0	work_2d(:,:) = 0.0
3002	0	do k=1,nz ; do j=js,je ; do i=is,ie
3003	0	work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
3004		enddo ; enddo ; enddo
3005	0	call post_data(CS%id_diabatic_diff_salt_tend_2d, work_2d, CS%diag)
3006		endif
3007		endif
3008		endif
3009
3010	0	end subroutine diagnose_diabatic_diff_tendency
3011
3012
3013		!> This routine diagnoses tendencies from application of boundary fluxes.
3014		!! These impacts are generally 3d, in particular for penetrative shortwave.
3015		!! Other fluxes contribute 3d in cases when the layers vanish or are very thin,
3016		!! in which case we distribute the flux into k > 1 layers.
3017	0	subroutine diagnose_boundary_forcing_tendency(tv, h, temp_old, saln_old, h_old, &
3018		dt, G, GV, US, CS)
3019		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
3020		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
3021		type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields
3022		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
3023		intent(in) :: h !< thickness after boundary flux application [H ~> m or kg m-2]
3024		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
3025		intent(in) :: temp_old !< temperature prior to boundary flux application [C ~> degC]
3026		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
3027		intent(in) :: saln_old !< salinity prior to boundary flux application [S ~> ppt]
3028		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
3029		intent(in) :: h_old !< thickness prior to boundary flux application [H ~> m or kg m-2]
3030		real, intent(in) :: dt !< time step [T ~> s]
3031		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3032		type(diabatic_CS), pointer :: CS !< module control structure
3033
3034		! Local variables
3035	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_3d ! A 3-d work array for diagnostics [various]
3036	0	real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array for diagnostics [various]
3037		real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
3038		integer :: i, j, k, is, ie, js, je, nz
3039
3040	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
3041	0	Idt = 0.0 ; if (dt > 0.0) Idt = 1. / dt
3042	0	work_3d(:,:,:) = 0.0
3043	0	work_2d(:,:) = 0.0
3044
3045		! Thickness tendency
3046	0	if (CS%id_boundary_forcing_h_tendency > 0) then
3047	0	do k=1,nz ; do j=js,je ; do i=is,ie
3048	0	work_3d(i,j,k) = (h(i,j,k) - h_old(i,j,k))*Idt
3049		enddo ; enddo ; enddo
3050	0	call post_data(CS%id_boundary_forcing_h_tendency, work_3d, CS%diag, alt_h=h_old)
3051		endif
3052
3053		! temperature tendency
3054	0	if (CS%id_boundary_forcing_temp_tend > 0) then
3055	0	do k=1,nz ; do j=js,je ; do i=is,ie
3056	0	work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*Idt
3057		enddo ; enddo ; enddo
3058	0	call post_data(CS%id_boundary_forcing_temp_tend, work_3d, CS%diag, alt_h=h_old)
3059		endif
3060
3061		! heat tendency
3062	0	if (CS%id_boundary_forcing_heat_tend > 0 .or. CS%id_boundary_forcing_heat_tend_2d > 0) then
3063	0	do k=1,nz ; do j=js,je ; do i=is,ie
3064	0	work_3d(i,j,k) = GV%H_to_RZ * tv%C_p * Idt * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * temp_old(i,j,k))
3065		enddo ; enddo ; enddo
3066	0	if (CS%id_boundary_forcing_heat_tend > 0) then
3067	0	call post_data(CS%id_boundary_forcing_heat_tend, work_3d, CS%diag, alt_h=h_old)
3068		endif
3069	0	if (CS%id_boundary_forcing_heat_tend_2d > 0) then
3070	0	work_2d(:,:) = 0.0
3071	0	do k=1,nz ; do j=js,je ; do i=is,ie
3072	0	work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
3073		enddo ; enddo ; enddo
3074	0	call post_data(CS%id_boundary_forcing_heat_tend_2d, work_2d, CS%diag)
3075		endif
3076		endif
3077
3078		! salinity tendency
3079	0	if (CS%id_boundary_forcing_saln_tend > 0) then
3080	0	do k=1,nz ; do j=js,je ; do i=is,ie
3081	0	work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*Idt
3082		enddo ; enddo ; enddo
3083	0	call post_data(CS%id_boundary_forcing_saln_tend, work_3d, CS%diag, alt_h=h_old)
3084		endif
3085
3086		! salt tendency
3087	0	if (CS%id_boundary_forcing_salt_tend > 0 .or. CS%id_boundary_forcing_salt_tend_2d > 0) then
3088	0	do k=1,nz ; do j=js,je ; do i=is,ie
3089	0	work_3d(i,j,k) = Idt * (h(i,j,k) * tv%S(i,j,k) - h_old(i,j,k) * saln_old(i,j,k))
3090		enddo ; enddo ; enddo
3091	0	if (CS%id_boundary_forcing_salt_tend > 0) then
3092	0	call post_data(CS%id_boundary_forcing_salt_tend, work_3d, CS%diag, alt_h=h_old)
3093		endif
3094	0	if (CS%id_boundary_forcing_salt_tend_2d > 0) then
3095	0	work_2d(:,:) = 0.0
3096	0	do k=1,nz ; do j=js,je ; do i=is,ie
3097	0	work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
3098		enddo ; enddo ; enddo
3099	0	call post_data(CS%id_boundary_forcing_salt_tend_2d, work_2d, CS%diag)
3100		endif
3101		endif
3102
3103	0	end subroutine diagnose_boundary_forcing_tendency
3104
3105
3106		!> This routine diagnoses tendencies for temperature and heat from frazil formation.
3107		!! This routine is called twice from within subroutine diabatic; at start and at
3108		!! end of the diabatic processes. The impacts from frazil are generally a function
3109		!! of depth. Hence, when checking heat budget, be sure to remove HFSIFRAZIL from HFDS in k=1.
3110	0	subroutine diagnose_frazil_tendency(tv, h, temp_old, dt, G, GV, US, CS)
3111		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
3112		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
3113		type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields
3114		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< thickness [H ~> m or kg m-2]
3115		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: temp_old !< temperature prior to frazil formation [C ~> degC]
3116		real, intent(in) :: dt !< time step [T ~> s]
3117		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3118		type(diabatic_CS), pointer :: CS !< module control structure
3119
3120	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_3d ! A 3-d work array for diagnostics [various]
3121	0	real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array for diagnostics [various]
3122		real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
3123		integer :: i, j, k, is, ie, js, je, nz
3124
3125	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
3126	0	Idt = 0.0 ; if (dt > 0.0) Idt = 1. / dt
3127
3128		! temperature tendency
3129	0	if (CS%id_frazil_temp_tend > 0) then
3130	0	do k=1,nz ; do j=js,je ; do i=is,ie
3131	0	work_3d(i,j,k) = Idt * (tv%T(i,j,k)-temp_old(i,j,k))
3132		enddo ; enddo ; enddo
3133	0	call post_data(CS%id_frazil_temp_tend, work_3d, CS%diag)
3134		endif
3135
3136		! heat tendency
3137	0	if (CS%id_frazil_heat_tend > 0 .or. CS%id_frazil_heat_tend_2d > 0) then
3138	0	do k=1,nz ; do j=js,je ; do i=is,ie
3139	0	work_3d(i,j,k) = GV%H_to_RZ * tv%C_p * h(i,j,k) * Idt * (tv%T(i,j,k)-temp_old(i,j,k))
3140		enddo ; enddo ; enddo
3141	0	if (CS%id_frazil_heat_tend > 0) call post_data(CS%id_frazil_heat_tend, work_3d, CS%diag)
3142
3143		! As a consistency check, we must have
3144		! FRAZIL_HEAT_TENDENCY_2d = HFSIFRAZIL
3145	0	if (CS%id_frazil_heat_tend_2d > 0) then
3146	0	work_2d(:,:) = 0.0
3147	0	do k=1,nz ; do j=js,je ; do i=is,ie
3148	0	work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
3149		enddo ; enddo ; enddo
3150	0	call post_data(CS%id_frazil_heat_tend_2d, work_2d, CS%diag)
3151		endif
3152		endif
3153
3154	0	end subroutine diagnose_frazil_tendency
3155
3156
3157		!> A simplified version of diabatic_driver_init that will allow
3158		!! tracer column functions to be called without allowing any
3159		!! of the diabatic processes to be used.
3160	0	subroutine adiabatic_driver_init(Time, G, param_file, diag, CS, &
3161		tracer_flow_CSp)
3162		type(time_type), intent(in) :: Time !< current model time
3163		type(ocean_grid_type), intent(in) :: G !< model grid structure
3164		type(param_file_type), intent(in) :: param_file !< the file to parse for parameter values
3165		type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
3166		type(diabatic_CS), pointer :: CS !< module control structure
3167		type(tracer_flow_control_CS), pointer :: tracer_flow_CSp !< pointer to control structure of the
3168		!! tracer flow control module
3169
3170		! This "include" declares and sets the variable "version".
3171		# include "version_variable.h"
3172		character(len=40) :: mdl = "MOM_diabatic_driver" ! This module's name.
3173
3174	0	if (associated(CS)) then
3175		call MOM_error(WARNING, "adiabatic_driver_init called with an "// &
3176	0	"associated control structure.")
3177	0	return
3178	0	else ; allocate(CS) ; endif
3179
3180	0	CS%diag => diag
3181	0	if (associated(tracer_flow_CSp)) CS%tracer_flow_CSp => tracer_flow_CSp
3182
3183		! Set default, read and log parameters
3184		call log_version(param_file, mdl, version, &
3185	0	"The following parameters are used for diabatic processes.")
3186
3187		! Check for any subsidiary parameters that are inconsistent with the adiabatic mode.
3188		call get_param(param_file, mdl, "SPONGE", CS%use_sponge, &
3189		"If true, sponges may be applied anywhere in the domain. "//&
3190		"The exact location and properties of those sponges are "//&
3191		"specified via calls to initialize_sponge and possibly "//&
3192	0	"set_up_sponge_field.", default=.false., do_not_log=.true.)
3193		call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", CS%use_energetic_PBL, &
3194		"If true, use an implied energetics planetary boundary "//&
3195		"layer scheme to determine the diffusivity and viscosity "//&
3196	0	"in the surface boundary layer.", default=.false., do_not_log=.true.)
3197		call get_param(param_file, mdl, "USE_KPP", CS%use_KPP, &
3198		"If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
3199		"to calculate diffusivities and non-local transport in the OBL.", &
3200	0	default=.false., do_not_log=.true.)
3201
3202	0	if (CS%use_sponge) call MOM_error(WARNING, &
3203	0	"When ADIABATIC = True, it is inconsistent to set SPONGE = True.")
3204	0	if (CS%use_energetic_PBL) call MOM_error(WARNING, &
3205	0	"When ADIABATIC = True, it is inconsistent to set ENERGETICS_SFC_PBL = True.")
3206	0	if (CS%use_KPP) call MOM_error(WARNING, &
3207	0	"When ADIABATIC = True, it is inconsistent to set USE_KPP = True.")
3208
3209	0	if (CS%use_sponge .or. CS%use_energetic_PBL .or. CS%use_KPP) &
3210	0	call MOM_error(FATAL, "adiabatic_driver_init is aborting due to inconsistent parameter settings.")
3211
3212		end subroutine adiabatic_driver_init
3213
3214
3215		!> This routine initializes the diabatic driver module.
3216	1	subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, diag, &
3217		ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, &
3218		ALE_sponge_CSp, oda_incupd_CSp, int_tide_CSp)
3219		type(time_type), target :: Time !< model time
3220		type(ocean_grid_type), intent(inout) :: G !< model grid structure
3221		type(verticalGrid_type), intent(in) :: GV !< model vertical grid structure
3222		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3223		type(param_file_type), intent(in) :: param_file !< file to parse for parameter values
3224		logical, intent(in) :: useALEalgorithm !< logical for whether to use ALE remapping
3225		type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output
3226		type(accel_diag_ptrs), intent(inout) :: ADp !< pointers to accelerations in momentum equations,
3227		!! to enable diagnostics, like energy budgets
3228		type(cont_diag_ptrs), intent(inout) :: CDp !< pointers to terms in continuity equations
3229		type(diabatic_CS), pointer :: CS !< module control structure
3230		type(tracer_flow_control_CS), pointer :: tracer_flow_CSp !< pointer to control structure of the
3231		!! tracer flow control module
3232		type(sponge_CS), pointer :: sponge_CSp !< pointer to the sponge module control structure
3233		type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< pointer to the ALE sponge module control structure
3234		type(oda_incupd_CS), pointer :: oda_incupd_CSp !< pointer to the ocean data assimilation incremental
3235		!! update module control structure
3236		type(int_tide_CS), pointer :: int_tide_CSp !< pointer to the internal tide structure
3237
3238		! Local variables
3239		real :: Kd ! A diffusivity used in the default for other tracer diffusivities [Z2 T-1 ~> m2 s-1]
3240		logical :: use_temperature
3241		character(len=20) :: EN1, EN2, EN3
3242
3243		! This "include" declares and sets the variable "version".
3244		# include "version_variable.h"
3245		character(len=40) :: mdl = "MOM_diabatic_driver" ! This module's name.
3246		character(len=48) :: thickness_units
3247		logical :: physical_OBL_scheme
3248		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nbands
3249	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
3250	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
3251
3252	1	CS%initialized = .true.
3253
3254	1	CS%diag => diag
3255	1	CS%Time => Time
3256
3257	1	if (associated(tracer_flow_CSp)) CS%tracer_flow_CSp => tracer_flow_CSp
3258	1	if (associated(sponge_CSp)) CS%sponge_CSp => sponge_CSp
3259	1	if (associated(ALE_sponge_CSp)) CS%ALE_sponge_CSp => ALE_sponge_CSp
3260	1	if (associated(oda_incupd_CSp)) CS%oda_incupd_CSp => oda_incupd_CSp
3261	1	if (associated(int_tide_CSp)) CS%int_tide_CSp => int_tide_CSp
3262
3263	1	CS%useALEalgorithm = useALEalgorithm
3264	1	CS%use_bulkmixedlayer = (GV%nkml > 0)
3265
3266		! Set default, read and log parameters
3267		call log_version(param_file, mdl, version, &
3268		"The following parameters are used for diabatic processes.", &
3269	1	log_to_all=.true., debugging=.true.)
3270		call get_param(param_file, mdl, "USE_LEGACY_DIABATIC_DRIVER", CS%use_legacy_diabatic, &
3271		"If true, use a legacy version of the diabatic subroutine. "//&
3272		"This is temporary and is needed to avoid change in answers.", &
3273	1	default=.true.)
3274		call get_param(param_file, mdl, "SPONGE", CS%use_sponge, &
3275		"If true, sponges may be applied anywhere in the domain. "//&
3276		"The exact location and properties of those sponges are "//&
3277		"specified via calls to initialize_sponge and possibly "//&
3278	1	"set_up_sponge_field.", default=.false.)
3279		call get_param(param_file, mdl, "ODA_INCUPD", CS%use_oda_incupd, &
3280		"If true, oda incremental updates will be applied "//&
3281	1	"everywhere in the domain.", default=.false.)
3282		call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
3283		"If true, temperature and salinity are used as state "//&
3284	1	"variables.", default=.true.)
3285		call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", CS%use_energetic_PBL, &
3286		"If true, use an implied energetics planetary boundary "//&
3287		"layer scheme to determine the diffusivity and viscosity "//&
3288	1	"in the surface boundary layer.", default=.false.)
3289		call get_param(param_file, mdl, "EPBL_IS_ADDITIVE", CS%ePBL_is_additive, &
3290		"If true, the diffusivity from ePBL is added to all "//&
3291		"other diffusivities. Otherwise, the larger of kappa-shear "//&
3292	1	"and ePBL diffusivities are used.", default=.true.)
3293		call get_param(param_file, mdl, "PRANDTL_EPBL", CS%ePBL_Prandtl, &
3294		"The Prandtl number used by ePBL to convert vertical diffusivities into "//&
3295	1	"viscosities.", default=1.0, units="nondim", do_not_log=.not.CS%use_energetic_PBL)
3296		call get_param(param_file, mdl, "USE_KPP", CS%use_KPP, &
3297		"If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
3298		"to calculate diffusivities and non-local transport in the OBL.", &
3299	1	default=.false., do_not_log=.true.)
3300	1	CS%use_CVMix_ddiff = CVMix_ddiff_is_used(param_file)
3301
3302	1	CS%use_kappa_shear = kappa_shear_is_used(param_file)
3303	1	CS%use_CVMix_shear = CVMix_shear_is_used(param_file)
3304
3305	1	if (CS%use_bulkmixedlayer) then
3306		call get_param(param_file, mdl, "ML_MIX_FIRST", CS%ML_mix_first, &
3307		"The fraction of the mixed layer mixing that is applied "//&
3308		"before interior diapycnal mixing. 0 by default.", &
3309	0	units="nondim", default=0.0)
3310	0	call get_param(param_file, mdl, "NKBL", CS%nkbl, default=2, do_not_log=.true.)
3311		else
3312	1	CS%ML_mix_first = 0.0
3313		endif
3314	1	if (use_temperature) then
3315		call get_param(param_file, mdl, "DO_GEOTHERMAL", CS%use_geothermal, &
3316	1	"If true, apply geothermal heating.", default=.false.)
3317		else
3318	0	CS%use_geothermal = .false.
3319		endif
3320		call get_param(param_file, mdl, "INTERNAL_TIDES", CS%use_int_tides, &
3321		"If true, use the code that advances a separate set of "//&
3322	1	"equations for the internal tide energy density.", default=.false.)
3323
3324		call get_param(param_file, mdl, "MASSLESS_MATCH_TARGETS", CS%massless_match_targets, &
3325		"If true, the temperature and salinity of massless layers "//&
3326		"are kept consistent with their target densities. "//&
3327		"Otherwise the properties of massless layers evolve "//&
3328		"diffusively to match massive neighboring layers.", &
3329	1	default=.true.)
3330
3331		call get_param(param_file, mdl, "AGGREGATE_FW_FORCING", CS%aggregate_FW_forcing, &
3332		"If true, the net incoming and outgoing fresh water fluxes are combined "//&
3333		"and applied as either incoming or outgoing depending on the sign of the net. "//&
3334		"If false, the net incoming fresh water flux is added to the model and "//&
3335		"thereafter the net outgoing is removed from the topmost non-vanished "//&
3336	1	"layers of the updated state.", default=.true.)
3337
3338		call get_param(param_file, mdl, "DEBUG", CS%debug, &
3339		"If true, write out verbose debugging data.", &
3340	1	default=.false., debuggingParam=.true.)
3341		call get_param(param_file, mdl, "DEBUG_CONSERVATION", CS%debugConservation, &
3342		"If true, monitor conservation and extrema.", &
3343	1	default=.false., debuggingParam=.true.)
3344
3345		call get_param(param_file, mdl, "DEBUG_ENERGY_REQ", CS%debug_energy_req, &
3346	1	"If true, debug the energy requirements.", default=.false., do_not_log=.true.)
3347		call get_param(param_file, mdl, "MIX_BOUNDARY_TRACERS", CS%mix_boundary_tracers, &
3348		"If true, mix the passive tracers in massless layers at "//&
3349		"the bottom into the interior as though a diffusivity of "//&
3350	1	"KD_MIN_TR were operating.", default=.true.)
3351		call get_param(param_file, mdl, "MIX_BOUNDARY_TRACER_ALE", CS%mix_boundary_tracer_ALE, &
3352		"If true and in ALE mode, mix the passive tracers in massless layers at "//&
3353		"the bottom into the interior as though a diffusivity of "//&
3354	1	"KD_MIN_TR were operating.", default=.false., do_not_log=.not.CS%useALEalgorithm)
3355
3356	1	if (CS%mix_boundary_tracers .or. CS%mix_boundary_tracer_ALE) then
3357	1	call get_param(param_file, mdl, "KD", Kd, units="m2 s-1", default=0.0, scale=US%m2_s_to_Z2_T)
3358		call get_param(param_file, mdl, "KD_MIN_TR", CS%Kd_min_tr, &
3359		"A minimal diffusivity that should always be applied to "//&
3360		"tracers, especially in massless layers near the bottom. "//&
3361		"The default is 0.1*KD.", &
3362	1	units="m2 s-1", default=0.1KdUS%Z2_T_to_m2_s, scale=GV%m2_s_to_HZ_T)
3363		call get_param(param_file, mdl, "KD_BBL_TR", CS%Kd_BBL_tr, &
3364		"A bottom boundary layer tracer diffusivity that will "//&
3365		"allow for explicitly specified bottom fluxes. The "//&
3366		"entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) "//&
3367		"over the same distance.", &
3368	1	units="m2 s-1", default=0., scale=GV%m2_s_to_HZ_T(US%Z_to_mGV%m_to_H))
3369		! The scaling factor here is usually equivalent to GV%m2_s_to_HZ_T*GV%Z_to_H.
3370		endif
3371
3372		call get_param(param_file, mdl, "TRACER_TRIDIAG", CS%tracer_tridiag, &
3373		"If true, use the passive tracer tridiagonal solver for T and S", &
3374	1	default=.false.)
3375
3376		call get_param(param_file, mdl, "MINIMUM_FORCING_DEPTH", CS%minimum_forcing_depth, &
3377		"The smallest depth over which forcing can be applied. This "//&
3378		"only takes effect when near-surface layers become thin "//&
3379		"relative to this scale, in which case the forcing tendencies "//&
3380		"scaled down by distributing the forcing over this depth scale.", &
3381	1	units="m", default=0.001, scale=GV%m_to_H)
3382		call get_param(param_file, mdl, "EVAP_CFL_LIMIT", CS%evap_CFL_limit, &
3383		"The largest fraction of a layer than can be lost to forcing "//&
3384		"(e.g. evaporation, sea-ice formation) in one time-step. The unused "//&
3385		"mass loss is passed down through the column.", &
3386	1	units="nondim", default=0.8)
3387
3388	1	if (CS%use_energetic_PBL .and. .not.CS%useALEalgorithm) &
3389		call MOM_error(FATAL, "diabatic_driver_init: "//&
3390	0	"ENERGETICS_SFC_PBL = True is only coded to work when USE_REGRIDDING = True.")
3391
3392		! Register all available diagnostics for this module.
3393	1	thickness_units = get_thickness_units(GV)
3394
3395		CS%id_ea_t = register_diag_field('ocean_model', 'ea_t', diag%axesTL, Time, &
3396	1	'Layer (heat) entrainment from above per timestep', 'm', conversion=GV%H_to_m)
3397		CS%id_eb_t = register_diag_field('ocean_model', 'eb_t', diag%axesTL, Time, &
3398	1	'Layer (heat) entrainment from below per timestep', 'm', conversion=GV%H_to_m)
3399		CS%id_ea_s = register_diag_field('ocean_model', 'ea_s', diag%axesTL, Time, &
3400	1	'Layer (salt) entrainment from above per timestep', 'm', conversion=GV%H_to_m)
3401		CS%id_eb_s = register_diag_field('ocean_model', 'eb_s', diag%axesTL, Time, &
3402	1	'Layer (salt) entrainment from below per timestep', 'm', conversion=GV%H_to_m)
3403		! used by layer diabatic
3404		CS%id_ea = register_diag_field('ocean_model', 'ea', diag%axesTL, Time, &
3405	1	'Layer entrainment from above per timestep', 'm', conversion=GV%H_to_m)
3406		CS%id_eb = register_diag_field('ocean_model', 'eb', diag%axesTL, Time, &
3407	1	'Layer entrainment from below per timestep', 'm', conversion=GV%H_to_m)
3408	1	if (.not.CS%useALEalgorithm) then
3409		CS%id_wd = register_diag_field('ocean_model', 'wd', diag%axesTi, Time, &
3410	0	'Diapycnal velocity', 'm s-1', conversion=GV%H_to_m*US%s_to_T)
3411	0	if (CS%id_wd > 0) call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1)
3412
3413		CS%id_dudt_dia = register_diag_field('ocean_model', 'dudt_dia', diag%axesCuL, Time, &
3414	0	'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2)
3415		CS%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, Time, &
3416	0	'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2)
3417
3418		CS%id_hf_dudt_dia_2d = register_diag_field('ocean_model', 'hf_dudt_dia_2d', diag%axesCu1, Time, &
3419		'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Diapycnal Mixing', &
3420	0	'm s-2', conversion=US%L_T2_to_m_s2)
3421	0	if (CS%id_hf_dudt_dia_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
3422
3423		CS%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, Time, &
3424		'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', &
3425	0	'm s-2', conversion=US%L_T2_to_m_s2)
3426	0	if (CS%id_hf_dvdt_dia_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
3427
3428	0	if ((CS%id_dudt_dia > 0) .or. (CS%id_hf_dudt_dia_2d > 0)) &
3429	0	call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz)
3430	0	if ((CS%id_dvdt_dia > 0) .or. (CS%id_hf_dvdt_dia_2d > 0)) &
3431	0	call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz)
3432		endif
3433
3434	1	if (use_temperature) then
3435		CS%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff", diag%axesTi, &
3436		Time, "Diffusive diapycnal temperature flux across interfaces", &
3437	1	units="degC m s-1", conversion=US%C_to_degCGV%H_to_mUS%s_to_T)
3438	1	if (.not.CS%useALEalgorithm) then
3439		CS%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv", diag%axesTi, &
3440		Time, "Advective diapycnal temperature flux across interfaces", &
3441	0	units="degC m s-1", conversion=US%C_to_degCGV%H_to_mUS%s_to_T)
3442		endif
3443		CS%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff", diag%axesTi, &
3444		Time, "Diffusive diapycnal salnity flux across interfaces", &
3445	1	units="psu m s-1", conversion=US%S_to_pptGV%H_to_mUS%s_to_T)
3446	1	if (.not.CS%useALEalgorithm) then
3447		CS%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv", diag%axesTi, &
3448		Time, "Advective diapycnal salnity flux across interfaces", &
3449	0	units="psu m s-1", conversion=US%S_to_pptGV%H_to_mUS%s_to_T)
3450		endif
3451
3452		CS%id_N2_dd = register_diag_field('ocean_model',"N2_diabatic", diag%axesTi, &
3453		Time, "Squared buoyancy frequency diagnosed after diffusion applied in thermodynamic timestep.", &
3454	1	"s-2", conversion=US%s_to_T**2)
3455		CS%id_N2_temp_dd = register_diag_field('ocean_model',"N2_temp_diabatic", diag%axesTi, &
3456		Time, "Squared buoyancy frequency due to temperature stratification diagnosed after diffusion applied "//&
3457	1	"in thermodynamic timestep.", "s-2", conversion=US%s_to_T**2)
3458		CS%id_N2_salt_dd = register_diag_field('ocean_model',"N2_salt_diabatic", diag%axesTi, &
3459		Time, "Squared buoyancy frequency due to salinity stratification diagnosed after diffusion applied "//&
3460	1	"in thermodynamic timestep.", "s-2", conversion=US%s_to_T**2)
3461	1	if (CS%id_N2_dd>0 .or. CS%id_N2_temp_dd>0 .or. CS%id_N2_salt_dd>0) CS%Use_N2_diag = .true.
3462
3463		CS%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, Time, &
3464		'Mixed layer depth (delta rho = 0.03)', units='m', conversion=US%Z_to_m, &
3465		cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', &
3466	1	cmor_standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t')
3467		CS%id_mlotstsq = register_diag_field('ocean_model', 'mlotstsq', diag%axesT1, Time, &
3468		long_name='Square of Ocean Mixed Layer Thickness Defined by Sigma T', &
3469		standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t', &
3470	1	units='m2', conversion=US%Z_to_m**2)
3471		CS%id_MLD_0125 = register_diag_field('ocean_model', 'MLD_0125', diag%axesT1, Time, &
3472	1	'Mixed layer depth (delta rho = 0.125)', 'm', conversion=US%Z_to_m)
3473		call get_param(param_file, mdl, "MLD_EN_VALS", CS%MLD_En_vals, &
3474		"The energy values used to compute MLDs. If not set (or all set to 0.), the "//&
3475		"default will overwrite to 25., 2500., 250000.", units='J/m2', &
3476	1	defaults=(/25., 2500., 250000./), scale=US%W_m2_to_RZ3_T3*US%s_to_T)
3477	1	write(EN1,'(F10.2)') CS%MLD_En_vals(1)US%RZ3_T3_to_W_m2US%T_to_s
3478	1	write(EN2,'(F10.2)') CS%MLD_En_vals(2)US%RZ3_T3_to_W_m2US%T_to_s
3479	1	write(EN3,'(F10.2)') CS%MLD_En_vals(3)US%RZ3_T3_to_W_m2US%T_to_s
3480		CS%id_MLD_EN1 = register_diag_field('ocean_model', 'MLD_EN1', diag%axesT1, Time, &
3481		'Mixed layer depth for energy value set to '//trim(EN1)//' J/m2 (Energy set by 1st MLD_EN_VALS)', &
3482	1	units='m', conversion=US%Z_to_m)
3483		CS%id_MLD_EN2 = register_diag_field('ocean_model', 'MLD_EN2', diag%axesT1, Time, &
3484		'Mixed layer depth for energy value set to '//trim(EN2)//' J/m2 (Energy set by 2nd MLD_EN_VALS)', &
3485	1	units='m', conversion=US%Z_to_m)
3486		CS%id_MLD_EN3 = register_diag_field('ocean_model', 'MLD_EN3', diag%axesT1, Time, &
3487		'Mixed layer depth for energy value set to '//trim(EN3)//' J/m2 (Energy set by 3rd MLD_EN_VALS)', &
3488	1	units='m', conversion=US%Z_to_m)
3489		call get_param(param_file, mdl, "BMLD_EN_VALS", CS%BMLD_En_vals, &
3490		"The energy values used to compute Bottom MLDs. If not set (or all set to 0.), the "//&
3491		"default will overwrite to 25., 2500., 250000.", units='J/m2', &
3492	1	defaults=(/25., 2500., 250000./), scale=US%W_m2_to_RZ3_T3*US%s_to_T)
3493	1	write(EN1,'(F10.2)') CS%BMLD_En_vals(1)US%RZ3_T3_to_W_m2US%T_to_s
3494	1	write(EN2,'(F10.2)') CS%BMLD_En_vals(2)US%RZ3_T3_to_W_m2US%T_to_s
3495	1	write(EN3,'(F10.2)') CS%BMLD_En_vals(3)US%RZ3_T3_to_W_m2US%T_to_s
3496		CS%id_BMLD_EN1 = register_diag_field('ocean_model', 'BMLD_EN1', diag%axesT1, Time, &
3497		'Bottom mixed layer depth for energy value set to '//trim(EN1)//' J/m2 (Energy set by 1st MLD_EN_VALS)', &
3498	1	units='m', conversion=US%Z_to_m)
3499		CS%id_BMLD_EN2 = register_diag_field('ocean_model', 'BMLD_EN2', diag%axesT1, Time, &
3500		'Bottom mixed layer depth for energy value set to '//trim(EN2)//' J/m2 (Energy set by 2nd MLD_EN_VALS)', &
3501	1	units='m', conversion=US%Z_to_m)
3502		CS%id_BMLD_EN3 = register_diag_field('ocean_model', 'BMLD_EN3', diag%axesT1, Time, &
3503		'Bottom mixed layer depth for energy value set to '//trim(EN3)//' J/m2 (Energy set by 3rd MLD_EN_VALS)', &
3504	1	units='m', conversion=US%Z_to_m)
3505		if ((CS%id_MLD_EN1>0).or. (CS%id_MLD_EN2>0).or. (CS%id_MLD_EN3>0).or. &
3506	1	(CS%id_BMLD_EN1>0).or.(CS%id_BMLD_EN2>0).or.(CS%id_BMLD_EN3>0)) then
3507		call get_param(param_file, mdl, "USE_OM4_MLD_EN_ITER", CS%use_OM4_MLD_En_iter, &
3508		"If true, uses an older set of iteration coefficients in computing the PE based "//&
3509		"surface MLD to reproduce OM4 era models. False uses an updated (general) method.",&
3510	0	default=.true.)
3511		endif
3512		CS%id_subMLN2 = register_diag_field('ocean_model', 'subML_N2', diag%axesT1, Time, &
3513	1	'Squared buoyancy frequency below mixed layer', units='s-2', conversion=US%s_to_T**2)
3514		CS%id_MLD_user = register_diag_field('ocean_model', 'MLD_user', diag%axesT1, Time, &
3515	1	'Mixed layer depth (used defined)', units='m', conversion=US%Z_to_m)
3516	1	if (CS%id_MLD_003 > 0) then
3517		call get_param(param_file, mdl, "HREF_FOR_MLD", CS%ref_h_mld, &
3518		"Reference depth used to calculate the potential density used to find the mixed layer depth "//&
3519	0	"based on a delta rho = 0.03 kg/m3.", units='m', default=0.0, scale=US%m_to_Z)
3520		CS%id_MLD_003_zr = register_diag_field('ocean_model', 'MLD_003_refZ', diag%axesT1, Time, &
3521	0	'Depth of reference density for MLD (delta rho = 0.03)', units='m', conversion=US%Z_to_m)
3522		CS%id_MLD_003_rr = register_diag_field('ocean_model', 'MLD_003_refRho', diag%axesT1, Time, &
3523	0	'Reference density for MLD (delta rho = 0.03)', units='kg/m3', conversion=US%R_to_kg_m3)
3524		endif
3525		endif
3526
3527	1	call KdWork_init(Time, G,GV,US,diag,CS%VBF,CS%Use_KdWork_diag)
3528	1	if (CS%Use_KdWork_diag.and.(.not.useALEalgorithm)) &
3529	0	call MOM_error(WARNING,"The KdWork diagnostics are not fully implemented for use in layer mode.")
3530	1	if (CS%Use_KdWork_diag.and.(CS%use_legacy_diabatic)) &
3531	0	call MOM_error(WARNING,"The KdWork diagnostics are only approximate with the legacy diabatic driver.")
3532
3533		call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", CS%MLDdensityDifference, &
3534		"The density difference used to determine a diagnostic mixed "//&
3535		"layer depth, MLD_user, following the definition of Levitus 1982. "//&
3536		"The MLD is the depth at which the density is larger than the "//&
3537		"surface density by the specified amount.", &
3538	1	units='kg/m3', default=0.1, scale=US%kg_m3_to_R)
3539		call get_param(param_file, mdl, "DIAG_DEPTH_SUBML_N2", CS%dz_subML_N2, &
3540		"The distance over which to calculate a diagnostic of the "//&
3541		"stratification at the base of the mixed layer.", &
3542	1	units='m', default=50.0, scale=US%m_to_Z)
3543
3544		! diagnostics for values prior to diabatic and prior to ALE
3545		CS%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, Time, &
3546	1	'Zonal velocity before diabatic forcing', 'm s-1', conversion=US%L_T_to_m_s)
3547		CS%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, Time, &
3548	1	'Meridional velocity before diabatic forcing', 'm s-1', conversion=US%L_T_to_m_s)
3549		CS%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, Time, &
3550		'Layer Thickness before diabatic forcing', &
3551	1	trim(thickness_units), conversion=GV%H_to_MKS, v_extensive=.true.)
3552		CS%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, Time, &
3553	1	'Interface Heights before diabatic forcing', 'm', conversion=US%Z_to_m)
3554	1	if (use_temperature) then
3555		CS%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, Time, &
3556	1	'Potential Temperature', 'degC', conversion=US%C_to_degC)
3557		CS%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, Time, &
3558	1	'Salinity', 'PSU', conversion=US%S_to_ppt)
3559		endif
3560
3561		CS%id_Kd_int = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, Time, &
3562	1	'Total diapycnal diffusivity at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
3563	1	if (CS%use_energetic_PBL) then
3564		CS%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, Time, &
3565	1	'ePBL diapycnal diffusivity at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
3566		endif
3567
3568		CS%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, Time, &
3569		'Total diapycnal diffusivity for heat at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s, &
3570		cmor_field_name='difvho', &
3571		cmor_standard_name='ocean_vertical_heat_diffusivity', &
3572	1	cmor_long_name='Ocean vertical heat diffusivity')
3573		CS%id_Kd_salt = register_diag_field('ocean_model', 'Kd_salt', diag%axesTi, Time, &
3574		'Total diapycnal diffusivity for salt at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s, &
3575		cmor_field_name='difvso', &
3576		cmor_standard_name='ocean_vertical_salt_diffusivity', &
3577	1	cmor_long_name='Ocean vertical salt diffusivity')
3578
3579		! CS%useKPP is set to True if KPP-scheme is to be used, False otherwise.
3580		! KPP_init() allocated CS%KPP_Csp and also sets CS%KPPisPassive
3581	1	CS%useKPP = KPP_init(param_file, G, GV, US, diag, Time, CS%KPP_CSp, passive=CS%KPPisPassive)
3582
3583		! Diagnostics for tendencies of temperature and salinity due to diabatic processes,
3584		! available only for ALE algorithm.
3585		! Diagnostics for tendencies of temperature and heat due to frazil
3586		CS%id_diabatic_diff_h = register_diag_field('ocean_model', 'diabatic_diff_h', diag%axesTL, Time, &
3587		'Cell thickness used during diabatic diffusion', &
3588	1	thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
3589
3590	1	if (CS%useALEalgorithm) then
3591		CS%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', &
3592		'diabatic_diff_temp_tendency', diag%axesTL, Time, &
3593	1	'Diabatic diffusion temperature tendency', 'degC s-1', conversion=US%C_to_degC*US%s_to_T)
3594	1	if (CS%id_diabatic_diff_temp_tend > 0) then
3595	0	CS%diabatic_diff_tendency_diag = .true.
3596		endif
3597
3598		CS%id_diabatic_diff_saln_tend = register_diag_field('ocean_model',&
3599		'diabatic_diff_saln_tendency', diag%axesTL, Time, &
3600	1	'Diabatic diffusion salinity tendency', 'psu s-1', conversion=US%S_to_ppt*US%s_to_T)
3601	1	if (CS%id_diabatic_diff_saln_tend > 0) then
3602	0	CS%diabatic_diff_tendency_diag = .true.
3603		endif
3604
3605		CS%id_diabatic_diff_heat_tend = register_diag_field('ocean_model', &
3606		'diabatic_heat_tendency', diag%axesTL, Time, &
3607		'Diabatic diffusion heat tendency', &
3608		'W m-2', conversion=US%QRZ_T_to_W_m2, cmor_field_name='opottempdiff', &
3609		cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'// &
3610		'due_to_parameterized_dianeutral_mixing', &
3611		cmor_long_name='Tendency of sea water potential temperature expressed as heat content '// &
3612		'due to parameterized dianeutral mixing', &
3613	1	v_extensive=.true.)
3614	1	if (CS%id_diabatic_diff_heat_tend > 0) then
3615	0	CS%diabatic_diff_tendency_diag = .true.
3616		endif
3617
3618		CS%id_diabatic_diff_salt_tend = register_diag_field('ocean_model', &
3619		'diabatic_salt_tendency', diag%axesTL, Time, &
3620		'Diabatic diffusion of salt tendency', &
3621		'kg m-2 s-1', conversion=US%S_to_ppt0.001GV%H_to_RZ*US%RZ_T_to_kg_m2s, &
3622		cmor_field_name='osaltdiff', &
3623		cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &
3624		'due_to_parameterized_dianeutral_mixing', &
3625		cmor_long_name='Tendency of sea water salinity expressed as salt content '// &
3626		'due to parameterized dianeutral mixing', &
3627	1	v_extensive=.true.)
3628	1	if (CS%id_diabatic_diff_salt_tend > 0) then
3629	0	CS%diabatic_diff_tendency_diag = .true.
3630		endif
3631
3632		! This diagnostic should equal to roundoff if all is working well.
3633		CS%id_diabatic_diff_heat_tend_2d = register_diag_field('ocean_model', &
3634		'diabatic_heat_tendency_2d', diag%axesT1, Time, &
3635		'Depth integrated diabatic diffusion heat tendency', &
3636		'W m-2', conversion=US%QRZ_T_to_W_m2, cmor_field_name='opottempdiff_2d', &
3637		cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'//&
3638		'due_to_parameterized_dianeutral_mixing_depth_integrated', &
3639		cmor_long_name='Tendency of sea water potential temperature expressed as heat content '//&
3640	1	'due to parameterized dianeutral mixing depth integrated')
3641	1	if (CS%id_diabatic_diff_heat_tend_2d > 0) then
3642	0	CS%diabatic_diff_tendency_diag = .true.
3643		endif
3644
3645		! This diagnostic should equal to roundoff if all is working well.
3646		CS%id_diabatic_diff_salt_tend_2d = register_diag_field('ocean_model', &
3647		'diabatic_salt_tendency_2d', diag%axesT1, Time, &
3648		'Depth integrated diabatic diffusion salt tendency', &
3649		'kg m-2 s-1', conversion=US%S_to_ppt0.001GV%H_to_RZ*US%RZ_T_to_kg_m2s, &
3650		cmor_field_name='osaltdiff_2d', &
3651		cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &
3652		'due_to_parameterized_dianeutral_mixing_depth_integrated', &
3653		cmor_long_name='Tendency of sea water salinity expressed as salt content '// &
3654	1	'due to parameterized dianeutral mixing depth integrated')
3655	1	if (CS%id_diabatic_diff_salt_tend_2d > 0) then
3656	0	CS%diabatic_diff_tendency_diag = .true.
3657		endif
3658
3659		! Diagnostics for tendencies of thickness temperature and salinity due to boundary forcing,
3660		! available only for ALE algorithm.
3661		! Diagnostics for tendencies of temperature and heat due to frazil
3662		CS%id_boundary_forcing_h = register_diag_field('ocean_model', 'boundary_forcing_h', diag%axesTL, Time, &
3663		'Cell thickness after applying boundary forcing', &
3664	1	thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
3665		CS%id_boundary_forcing_h_tendency = register_diag_field('ocean_model', &
3666		'boundary_forcing_h_tendency', diag%axesTL, Time, &
3667		'Cell thickness tendency due to boundary forcing', &
3668	1	trim(thickness_units)//" s-1", conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.)
3669	1	if (CS%id_boundary_forcing_h_tendency > 0) then
3670	0	CS%boundary_forcing_tendency_diag = .true.
3671		endif
3672
3673		CS%id_boundary_forcing_temp_tend = register_diag_field('ocean_model',&
3674		'boundary_forcing_temp_tendency', diag%axesTL, Time, &
3675	1	'Boundary forcing temperature tendency', 'degC s-1', conversion=US%C_to_degC*US%s_to_T)
3676	1	if (CS%id_boundary_forcing_temp_tend > 0) then
3677	0	CS%boundary_forcing_tendency_diag = .true.
3678		endif
3679
3680		CS%id_boundary_forcing_saln_tend = register_diag_field('ocean_model',&
3681		'boundary_forcing_saln_tendency', diag%axesTL, Time, &
3682	1	'Boundary forcing saln tendency', 'psu s-1', conversion=US%S_to_ppt*US%s_to_T)
3683	1	if (CS%id_boundary_forcing_saln_tend > 0) then
3684	0	CS%boundary_forcing_tendency_diag = .true.
3685		endif
3686
3687		CS%id_boundary_forcing_heat_tend = register_diag_field('ocean_model',&
3688		'boundary_forcing_heat_tendency', diag%axesTL, Time, &
3689		'Boundary forcing heat tendency', &
3690	1	'W m-2', conversion=US%QRZ_T_to_W_m2, v_extensive=.true.)
3691	1	if (CS%id_boundary_forcing_heat_tend > 0) then
3692	0	CS%boundary_forcing_tendency_diag = .true.
3693		endif
3694
3695		CS%id_boundary_forcing_salt_tend = register_diag_field('ocean_model',&
3696		'boundary_forcing_salt_tendency', diag%axesTL, Time, &
3697		'Boundary forcing salt tendency', &
3698		'kg m-2 s-1', conversion=US%S_to_ppt0.001GV%H_to_RZ*US%RZ_T_to_kg_m2s, &
3699	1	v_extensive = .true.)
3700	1	if (CS%id_boundary_forcing_salt_tend > 0) then
3701	0	CS%boundary_forcing_tendency_diag = .true.
3702		endif
3703
3704		! This diagnostic should equal to surface heat flux if all is working well.
3705		CS%id_boundary_forcing_heat_tend_2d = register_diag_field('ocean_model',&
3706		'boundary_forcing_heat_tendency_2d', diag%axesT1, Time, &
3707		'Depth integrated boundary forcing of ocean heat', &
3708	1	'W m-2', conversion=US%QRZ_T_to_W_m2)
3709	1	if (CS%id_boundary_forcing_heat_tend_2d > 0) then
3710	0	CS%boundary_forcing_tendency_diag = .true.
3711		endif
3712
3713		! This diagnostic should equal to surface salt flux if all is working well.
3714		CS%id_boundary_forcing_salt_tend_2d = register_diag_field('ocean_model',&
3715		'boundary_forcing_salt_tendency_2d', diag%axesT1, Time, &
3716		'Depth integrated boundary forcing of ocean salt', &
3717	1	'kg m-2 s-1', conversion=US%S_to_ppt0.001GV%H_to_RZ*US%RZ_T_to_kg_m2s)
3718	1	if (CS%id_boundary_forcing_salt_tend_2d > 0) then
3719	0	CS%boundary_forcing_tendency_diag = .true.
3720		endif
3721		endif
3722
3723		! diagnostics for tendencies of temp and heat due to frazil
3724		CS%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, Time, &
3725		long_name='Cell Thickness', standard_name='cell_thickness', &
3726	1	units=thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
3727
3728		! diagnostic for tendency of temp due to frazil
3729		CS%id_frazil_temp_tend = register_diag_field('ocean_model',&
3730		'frazil_temp_tendency', diag%axesTL, Time, &
3731	1	'Temperature tendency due to frazil formation', 'degC s-1', conversion=US%C_to_degC*US%s_to_T)
3732	1	if (CS%id_frazil_temp_tend > 0) then
3733	0	CS%frazil_tendency_diag = .true.
3734		endif
3735
3736		! diagnostic for tendency of heat due to frazil
3737		CS%id_frazil_heat_tend = register_diag_field('ocean_model',&
3738		'frazil_heat_tendency', diag%axesTL, Time, &
3739		'Heat tendency due to frazil formation', &
3740	1	'W m-2', conversion=US%QRZ_T_to_W_m2, v_extensive=.true.)
3741	1	if (CS%id_frazil_heat_tend > 0) then
3742	0	CS%frazil_tendency_diag = .true.
3743		endif
3744
3745		! If all is working properly, this diagnostic should equal to hfsifrazil.
3746		CS%id_frazil_heat_tend_2d = register_diag_field('ocean_model',&
3747		'frazil_heat_tendency_2d', diag%axesT1, Time, &
3748		'Depth integrated heat tendency due to frazil formation', &
3749	1	'W m-2', conversion=US%QRZ_T_to_W_m2)
3750	1	if (CS%id_frazil_heat_tend_2d > 0) then
3751	0	CS%frazil_tendency_diag = .true.
3752		endif
3753
3754		! CS%use_CVMix_conv is set to True if CVMix convection will be used, otherwise it is False.
3755	1	CS%use_CVMix_conv = CVMix_conv_init(Time, G, GV, US, param_file, diag, CS%CVMix_conv)
3756
3757		call entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS%entrain_diffusive, &
3758	1	just_read_params=CS%useALEalgorithm)
3759
3760		! initialize the geothermal heating module
3761	1	if (CS%use_geothermal) &
3762	1	call geothermal_init(Time, G, GV, US, param_file, diag, CS%geothermal, useALEalgorithm)
3763
3764		! initialize module for internal tide induced mixing
3765	1	if (CS%use_int_tides) then
3766		call int_tide_input_init(Time, G, GV, US, param_file, diag, CS%int_tide_input_CSp, &
3767	0	CS%int_tide_input)
3768	0	call internal_tides_init(Time, G, GV, US, param_file, diag, int_tide_CSp)
3769		endif
3770
3771		!if (associated(int_tide_CSp)) CS%int_tide_CSp => int_tide_CSp
3772
3773	1	physical_OBL_scheme = (CS%use_bulkmixedlayer .or. CS%use_KPP .or. CS%use_energetic_PBL)
3774		! initialize module for setting diffusivities
3775		call set_diffusivity_init(Time, G, GV, US, param_file, diag, CS%set_diff_CSp, CS%int_tide_CSp, &
3776		halo_TS=CS%halo_TS_diff, double_diffuse=CS%double_diffuse, &
3777	1	physical_OBL_scheme=physical_OBL_scheme)
3778
3779	1	CS%halo_diabatic = CS%halo_TS_diff
3780	1	if (CS%use_int_tides) CS%halo_diabatic = max(CS%halo_TS_diff, 2)
3781
3782	1	if (CS%useKPP .and. (CS%double_diffuse .and. .not.CS%use_CVMix_ddiff)) &
3783		call MOM_error(FATAL, 'diabatic_driver_init: DOUBLE_DIFFUSION (old method) does not work '//&
3784	0	'with KPP. Please set DOUBLE_DIFFUSION=False and USE_CVMIX_DDIFF=True.')
3785
3786		! set up the clocks for this module
3787	1	id_clock_entrain = cpu_clock_id('(Ocean diabatic entrain)', grain=CLOCK_MODULE)
3788	1	if (CS%use_bulkmixedlayer) &
3789	0	id_clock_mixedlayer = cpu_clock_id('(Ocean mixed layer)', grain=CLOCK_MODULE)
3790	1	id_clock_remap = cpu_clock_id('(Ocean vert remap)', grain=CLOCK_MODULE)
3791	1	if (CS%use_geothermal) &
3792	1	id_clock_geothermal = cpu_clock_id('(Ocean geothermal)', grain=CLOCK_ROUTINE)
3793	1	id_clock_set_diffusivity = cpu_clock_id('(Ocean set_diffusivity)', grain=CLOCK_MODULE)
3794	1	id_clock_kpp = cpu_clock_id('(Ocean KPP)', grain=CLOCK_MODULE)
3795	1	id_clock_tracers = cpu_clock_id('(Ocean tracer_columns)', grain=CLOCK_MODULE_DRIVER+5)
3796	1	if (CS%use_sponge) &
3797	0	id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=CLOCK_MODULE)
3798	1	if (CS%use_oda_incupd) &
3799	0	id_clock_oda_incupd = cpu_clock_id('(Ocean inc. update data assimilation)', grain=CLOCK_MODULE)
3800	1	id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=CLOCK_ROUTINE)
3801	1	id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=CLOCK_ROUTINE)
3802	1	id_clock_differential_diff = -1 ; if (CS%double_diffuse .and. .not.CS%use_CVMix_ddiff) &
3803	0	id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=CLOCK_ROUTINE)
3804
3805		! initialize the auxiliary diabatic driver module
3806		call diabatic_aux_init(Time, G, GV, US, param_file, diag, CS%diabatic_aux_CSp, &
3807	1	CS%useALEalgorithm, CS%use_energetic_PBL)
3808
3809		! initialize the boundary layer modules
3810	1	if (CS%use_bulkmixedlayer) &
3811	0	call bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS%bulkmixedlayer)
3812	1	if (CS%use_energetic_PBL) &
3813	1	call energetic_PBL_init(Time, G, GV, US, param_file, diag, CS%ePBL)
3814
3815	1	call regularize_layers_init(Time, G, GV, param_file, diag, CS%regularize_layers)
3816
3817	1	if (CS%debug_energy_req) &
3818	0	call diapyc_energy_req_init(Time, G, GV, US, param_file, diag, CS%diapyc_en_rec_CSp)
3819
3820		! obtain information about the number of bands for penetrative shortwave
3821	1	if (use_temperature) then
3822	1	call get_param(param_file, mdl, "PEN_SW_NBANDS", nbands, default=1)
3823	1	if (nbands > 0) then
3824	1	allocate(CS%optics)
3825	1	call opacity_init(Time, G, GV, US, param_file, diag, CS%opacity, CS%optics)
3826		endif
3827		endif
3828
3829		! Initialize the diagnostic grid storage
3830	1	call diag_grid_storage_init(CS%diag_grids_prev, G, GV, diag)
3831
3832	1	end subroutine diabatic_driver_init
3833
3834		!> Routine to register restarts, pass-through to children modules
3835	1	subroutine register_diabatic_restarts(G, GV, US, param_file, int_tide_CSp, restart_CSp, CS)
3836		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
3837		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
3838		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3839		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
3840		type(int_tide_CS), pointer :: int_tide_CSp !< Internal tide control structure
3841		type(MOM_restart_CS), pointer :: restart_CSp !< MOM restart control structure
3842		type(diabatic_CS), pointer :: CS !< module control structure
3843
3844		logical :: use_int_tides
3845
3846	1	if (associated(CS)) then
3847		call MOM_error(WARNING, "diabatic_driver_init called with an "// &
3848	0	"associated control structure.")
3849	0	return
3850		else
3851	1	allocate(CS)
3852		endif
3853
3854	1	use_int_tides=.false.
3855
3856	1	call read_param(param_file, "INTERNAL_TIDES", use_int_tides)
3857
3858	1	if (use_int_tides) then
3859	0	call register_int_tide_restarts(G, GV, US, param_file, int_tide_CSp, restart_CSp)
3860		endif
3861
3862	1	call register_KPP_restarts(G, param_file, restart_CSp, CS%KPP_CSp)
3863
3864		end subroutine register_diabatic_restarts
3865
3866		!> Routine to close the diabatic driver module
3867	1	subroutine diabatic_driver_end(CS)
3868		type(diabatic_CS), intent(inout) :: CS !< module control structure
3869
3870	1	if (associated(CS%VBF)) then
3871	1	call KdWork_end(CS%VBF)
3872		endif
3873
3874	1	if (associated(CS%optics)) then
3875	1	call opacity_end(CS%opacity, CS%optics)
3876	1	deallocate(CS%optics)
3877		endif
3878
3879	1	if (CS%debug_energy_req) &
3880	0	call diapyc_energy_req_end(CS%diapyc_en_rec_CSp)
3881
3882	1	if (CS%use_energetic_PBL) &
3883	1	call energetic_PBL_end(CS%ePBL)
3884
3885	1	call diabatic_aux_end(CS%diabatic_aux_CSp)
3886
3887	1	call set_diffusivity_end(CS%set_diff_CSp)
3888
3889	1	deallocate(CS%set_diff_CSp)
3890
3891	1	if (CS%use_geothermal) &
3892	1	call geothermal_end(CS%geothermal)
3893
3894	1	if (CS%useKPP) &
3895	0	call KPP_end(CS%KPP_CSp)
3896
3897		! GMM, the following is commented out because arrays in
3898		! CS%diag_grids_prev are neither pointers or allocatables
3899		! and, therefore, cannot be deallocated.
3900
3901		!call diag_grid_storage_end(CS%diag_grids_prev)
3902	1	end subroutine diabatic_driver_end
3903
3904
3905		!> \namespace mom_diabatic_driver
3906		!!
3907		!! By Robert Hallberg, Alistair Adcroft, and Stephen Griffies
3908		!!
3909		!! This program contains the subroutine that, along with the
3910		!! subroutines that it calls, implements diapycnal mass and momentum
3911		!! fluxes and a bulk mixed layer. The diapycnal diffusion can be
3912		!! used without the bulk mixed layer.
3913		!!
3914		!! \section section_diabatic Outline of MOM diabatic
3915		!!
3916		!! * diabatic first determines the (diffusive) diapycnal mass fluxes
3917		!! based on the convergence of the buoyancy fluxes within each layer.
3918		!!
3919		!! * The dual-stream entrainment scheme of MacDougall and Dewar (JPO,
3920		!! 1997) is used for combined diapycnal advection and diffusion,
3921		!! calculated implicitly and potentially with the Richardson number
3922		!! dependent mixing, as described by Hallberg (MWR, 2000).
3923		!!
3924		!! * Diapycnal advection is the residual of diapycnal diffusion,
3925		!! so the fully implicit upwind differencing scheme that is used is
3926		!! entirely appropriate.
3927		!!
3928		!! * The downward buoyancy flux in each layer is determined from
3929		!! an implicit calculation based on the previously
3930		!! calculated flux of the layer above and an estimated flux in the
3931		!! layer below. This flux is subject to the following conditions:
3932		!! (1) the flux in the top and bottom layers are set by the boundary
3933		!! conditions, and (2) no layer may be driven below a minimal thickness.
3934		!! If there is a bulk mixed layer, the buffer layer is treated
3935		!! as a fixed density layer with vanishingly small diffusivity.
3936		!!
3937		!! diabatic takes 5 arguments: the two velocities (u and v), the
3938		!! thicknesses (h), a structure containing the forcing fields, and
3939		!! the length of time over which to act (dt). The velocities and
3940		!! thickness are taken as inputs and modified within the subroutine.
3941		!! There is no limit on the time step.
3942
3943	0	end module MOM_diabatic_driver