MOM_diabatic_driver.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.

2! See the LICENSE file for licensing information.

3! SPDX-License-Identifier: Apache-2.0

4

5!> This routine drives the diabatic/dianeutral physics for MOM

6module mom_diabatic_driver

7

8use mom_bulk_mixed_layer, only : bulkmixedlayer, bulkmixedlayer_init, bulkmixedlayer_cs

9use mom_debugging, only : hchksum

10use mom_checksum_packages, only : mom_state_chksum, mom_state_stats

11use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end

12use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine

13use mom_cvmix_shear, only : cvmix_shear_is_used

14use mom_cvmix_ddiff, only : cvmix_ddiff_is_used

15use mom_diabatic_aux, only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_cs

16use mom_diabatic_aux, only : make_frazil, adjust_salt, differential_diffuse_t_s, tridiagts

17use mom_diabatic_aux, only : tridiagts_eulerian, find_uv_at_h

18use mom_diabatic_aux, only : applyboundaryfluxesinout, set_pen_shortwave

19use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr

20use mom_diag_mediator, only : post_product_sum_u, post_product_sum_v

21use mom_diag_mediator, only : diag_ctrl, time_type, diag_update_remap_grids

22use mom_diag_mediator, only : diag_ctrl, query_averaging_enabled, enable_averages, disable_averaging

23use mom_diag_mediator, only : diag_grid_storage, diag_grid_storage_init, diag_grid_storage_end

24use mom_diag_mediator, only : diag_copy_diag_to_storage, diag_copy_storage_to_diag

25use mom_diag_mediator, only : diag_save_grids, diag_restore_grids

26use mom_diagnose_mld, only : diagnosemldbydensitydifference, diagnosemldbyenergy

27use mom_diagnose_kdwork, only : vbf_cs, kdwork_init, kdwork_end, kdwork_diagnostics

28use mom_diagnose_kdwork, only : allocate_vbf_cs, deallocate_vbf_cs

29use mom_diapyc_energy_req, only : diapyc_energy_req_init, diapyc_energy_req_end

30use mom_diapyc_energy_req, only : diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_cs

31use mom_cvmix_conv, only : cvmix_conv_init, cvmix_conv_cs

32use mom_cvmix_conv, only : calculate_cvmix_conv

33use mom_domains, only : pass_var, to_west, to_south, to_all, omit_corners

34use mom_domains, only : create_group_pass, do_group_pass, group_pass_type

35use mom_energetic_pbl, only : energetic_pbl, energetic_pbl_init

36use mom_energetic_pbl, only : energetic_pbl_end, energetic_pbl_cs

37use mom_energetic_pbl, only : energetic_pbl_get_mld

38use mom_entrain_diffusive, only : entrainment_diffusive, entrain_diffusive_init

39use mom_entrain_diffusive, only : entrain_diffusive_cs

40use mom_eos, only : calculate_density, calculate_density_derivs, calculate_tfreeze

41use mom_eos, only : calculate_specific_vol_derivs, eos_domain

42use mom_error_handler, only : mom_error, fatal, warning, calltree_showquery,mom_mesg

43use mom_error_handler, only : calltree_enter, calltree_leave, calltree_waypoint

44use mom_file_parser, only : get_param, log_version, param_file_type, read_param

45use mom_forcing_type, only : forcing, mom_forcing_chksum, find_ustar

46use mom_forcing_type, only : calculatebuoyancyflux2d, forcing_singlepointprint

47use mom_geothermal, only : geothermal_entraining, geothermal_in_place

48use mom_geothermal, only : geothermal_init, geothermal_end, geothermal_cs

49use mom_grid, only : ocean_grid_type

50use mom_int_tide_input, only : set_int_tide_input, int_tide_input_init

51use mom_int_tide_input, only : int_tide_input_end, int_tide_input_cs, int_tide_input_type

52use mom_interface_heights, only : find_eta, calc_derived_thermo, thickness_to_dz

53use mom_interface_heights, only : convert_mld_to_ml_thickness

54use mom_internal_tides, only : propagate_int_tide, register_int_tide_restarts

55use mom_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_cs

56use mom_kappa_shear, only : kappa_shear_is_used

57use mom_cvmix_kpp, only : kpp_cs, kpp_init, kpp_compute_bld, kpp_calculate

58use mom_cvmix_kpp, only : kpp_end, kpp_get_bld, register_kpp_restarts

59use mom_cvmix_kpp, only : kpp_nonlocaltransport_temp, kpp_nonlocaltransport_saln

60use mom_oda_incupd, only : apply_oda_incupd, oda_incupd_cs

61use mom_opacity, only : opacity_init, opacity_end, opacity_cs

62use mom_opacity, only : absorbremainingsw, optics_type, optics_nbands

63use mom_open_boundary, only : ocean_obc_type

64use mom_regularize_layers, only : regularize_layers, regularize_layers_init, regularize_layers_cs

65use mom_restart, only : mom_restart_cs

66use mom_set_diffusivity, only : set_diffusivity, set_bbl_tke

67use mom_set_diffusivity, only : set_diffusivity_init, set_diffusivity_end

68use mom_set_diffusivity, only : set_diffusivity_cs

69use mom_sponge, only : apply_sponge, sponge_cs

70use mom_ale_sponge, only : apply_ale_sponge, ale_sponge_cs

71use mom_time_manager, only : time_type, real_to_time, operator(-), operator(<=)

72use mom_tracer_flow_control, only : call_tracer_column_fns, tracer_flow_control_cs

73use mom_tracer_diabatic, only : tracer_vertdiff, tracer_vertdiff_eulerian

74use mom_unit_scaling, only : unit_scale_type

75use mom_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs

76use mom_variables, only : cont_diag_ptrs, mom_thermovar_chksum, p3d

77use mom_verticalgrid, only : verticalgrid_type, get_thickness_units

78use mom_wave_interface, only : wave_parameters_cs

79use mom_stochastics, only : stochastic_cs

80

81implicit none ; private

82

83#include <MOM_memory.h>

84

85public diabatic

86public diabatic_driver_init

87public diabatic_driver_end

88public extract_diabatic_member

89public adiabatic

90public adiabatic_driver_init

91public 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

99type, 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 logical :: mld_param_003 = .false. !< Logical flag if MLD in brine plume should use the 0.03 mixed layer depth

191 logical :: mld_param_en1 = .false. !< Logical flag if MLD in brine plume should use the EN1 mixed layer depth

192 logical :: mld_param_epbl = .false.!< Logical flag if MLD in brine plume should use the ePBL boundary layer depth

193

194 !>@{ Diagnostic IDs

195 integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic

196 integer :: id_ea_t = -1, id_eb_t = -1, id_ea_s = -1, id_eb_s = -1

197 integer :: id_kd_heat = -1, id_kd_salt = -1, id_kd_int = -1, id_kd_epbl = -1

198 integer :: id_tdif = -1, id_sdif = -1, id_tadv = -1, id_sadv = -1

199 integer :: id_n2_dd = -1, id_n2_salt_dd = -1, id_n2_temp_dd

200 ! These are handles to diagnostics related to the mixed layer properties.

201 integer :: id_mld_003 = -1, id_mld_0125 = -1, id_mld_user = -1, id_mlotstsq = -1

202 integer :: id_mld_003_zr = -1, id_mld_003_rr = -1

203 integer :: id_mld_en1 = -1, id_mld_en2 = -1, id_mld_en3 = -1, id_submln2 = -1

204 integer :: id_bmld_en1 = -1, id_bmld_en2 = -1, id_bmld_en3 = -1

205

206 ! These are handles to diagnostics that are only available in non-ALE layered mode.

207 integer :: id_wd = -1

208 integer :: id_dudt_dia = -1, id_dvdt_dia = -1

209 integer :: id_hf_dudt_dia_2d = -1, id_hf_dvdt_dia_2d = -1

210

211 ! diagnostic for fields prior to applying diapycnal physics

212 integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1

213 integer :: id_t_predia = -1, id_s_predia = -1, id_e_predia = -1

214

215 integer :: id_diabatic_diff_temp_tend = -1

216 integer :: id_diabatic_diff_saln_tend = -1

217 integer :: id_diabatic_diff_heat_tend = -1

218 integer :: id_diabatic_diff_salt_tend = -1

219 integer :: id_diabatic_diff_heat_tend_2d = -1

220 integer :: id_diabatic_diff_salt_tend_2d = -1

221 integer :: id_diabatic_diff_h = -1

222

223 integer :: id_boundary_forcing_h = -1

224 integer :: id_boundary_forcing_h_tendency = -1

225 integer :: id_boundary_forcing_temp_tend = -1

226 integer :: id_boundary_forcing_saln_tend = -1

227 integer :: id_boundary_forcing_heat_tend = -1

228 integer :: id_boundary_forcing_salt_tend = -1

229 integer :: id_boundary_forcing_heat_tend_2d = -1

230 integer :: id_boundary_forcing_salt_tend_2d = -1

231

232 integer :: id_frazil_h = -1

233 integer :: id_frazil_temp_tend = -1

234 integer :: id_frazil_heat_tend = -1

235 integer :: id_frazil_heat_tend_2d = -1

236 !>@}

237

238 logical :: diabatic_diff_tendency_diag = .false. !< If true calculate diffusive tendency diagnostics

239 logical :: boundary_forcing_tendency_diag = .false. !< If true calculate frazil diagnostics

240 logical :: frazil_tendency_diag = .false. !< If true calculate frazil tendency diagnostics

241

242 type(diabatic_aux_cs), pointer :: diabatic_aux_csp => null() !< Control structure for a child module

243 type(int_tide_input_cs), pointer :: int_tide_input_csp => null() !< Control structure for a child module

244 type(int_tide_input_type), pointer :: int_tide_input => null() !< Control structure for a child module

245 type(set_diffusivity_cs), pointer :: set_diff_csp => null() !< Control structure for a child module

246 type(sponge_cs), pointer :: sponge_csp => null() !< Control structure for a child module

247 type(ale_sponge_cs), pointer :: ale_sponge_csp => null() !< Control structure for a child module

248 type(tracer_flow_control_cs), pointer :: tracer_flow_csp => null() !< Control structure for a child module

249 type(optics_type), pointer :: optics => null() !< Control structure for a child module

250 type(kpp_cs), pointer :: kpp_csp => null() !< Control structure for a child module

251 type(diapyc_energy_req_cs), pointer :: diapyc_en_rec_csp => null() !< Control structure for a child module

252 type(oda_incupd_cs), pointer :: oda_incupd_csp => null() !< Control structure for a child module

253 type(int_tide_cs), pointer :: int_tide_csp => null() !< Control structure for a child module

254 type(vbf_cs), pointer :: vbf => null() !< Control structure for a child module

255

256

257 type(bulkmixedlayer_cs) :: bulkmixedlayer !< Bulk mixed layer control structure

258 type(cvmix_conv_cs) :: cvmix_conv !< CVMix convection control structure

259 type(energetic_pbl_cs) :: epbl !< Energetic PBL control structure

260 type(entrain_diffusive_cs) :: entrain_diffusive !< Diffusive entrainment control structure

261 type(geothermal_cs) :: geothermal !< Geothermal control structure

262 type(opacity_cs) :: opacity !< Opacity control structure

263 type(regularize_layers_cs) :: regularize_layers !< Regularize layer control structure

264

265 type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass

266 type(diag_grid_storage) :: diag_grids_prev !< Stores diagnostic grids at some previous point in the algorithm

267

268 type(time_type), pointer :: time !< Pointer to model time (needed for sponges)

269end type diabatic_cs

270

271!>@{ clock ids

272integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity

273integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge

274integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap

275integer :: id_clock_kpp, id_clock_oda_incupd

276!>@}

277

278contains

279

280!> This subroutine imposes the diapycnal mass fluxes and the

281!! accompanying diapycnal advection of momentum and tracers.

282subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &

283 G, GV, US, CS, stoch_CS, OBC, Waves)

284 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure

285 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure

286 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]

287 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]

288 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2]

289 type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields

290 !! unused have NULL ptrs

291 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: bld !< Active mixed layer depth [Z ~> m]

292 type(forcing), intent(inout) :: fluxes !< points to forcing fields

293 !! unused fields have NULL ptrs

294 type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities,

295 !! BBL properties and related fields

296 type(accel_diag_ptrs), intent(inout) :: adp !< Points to accelerations in momentum

297 !! equations, to enable the later derived

298 !! diagnostics, like energy budgets

299 type(cont_diag_ptrs), intent(inout) :: cdp !< points to terms in continuity equations

300 real, intent(in) :: dt !< time increment [T ~> s]

301 type(time_type), intent(in) :: time_end !< Time at the end of the interval

302 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

303 type(diabatic_cs), pointer :: cs !< module control structure

304 type(stochastic_cs), pointer :: stoch_cs !< stochastic control structure

305 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure.

306 type(wave_parameters_cs), pointer :: waves !< Surface gravity waves

307

308 ! local variables

309 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &

310 eta ! Interface heights before diapycnal mixing [Z ~> m]

311 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: temp_diag ! Previous temperature for diagnostics [C ~> degC]

312 real, dimension(SZI_(G)) :: t_freeze, & ! The freezing potential temperature at the current salinity [C ~> degC].

313 ps ! Surface pressure [R L2 T-2 ~> Pa]

314 real, dimension(SZI_(G),SZK_(GV)) :: &

315 pressure ! The pressure at the middle of each layer [R L2 T-2 ~> Pa].

316 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]

317 integer :: i, j, k, is, ie, js, je, nz

318 logical :: showcalltree ! If true, show the call tree

319

320 real, allocatable, dimension(:,:,:) :: h_in ! thickness before thermodynamics [H ~> m or kg m-2]

321 real, allocatable, dimension(:,:,:) :: t_in ! temperature before thermodynamics [C ~> degC]

322 real, allocatable, dimension(:,:,:) :: s_in ! salinity before thermodynamics [S ~> ppt]

323

324 if (gv%ke == 1) return

325

326 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

327

328 if (.not. associated(cs)) call mom_error(fatal, "MOM_diabatic_driver: "// &

329 "Module must be initialized before it is used.")

330

331 if (.not. cs%initialized) call mom_error(fatal, "MOM_diabatic_driver: "// &

332 "Module must be initialized before it is used.")

333

334 if (dt == 0.0) call mom_error(fatal, "MOM_diabatic_driver: "// &

335 "diabatic was called with a zero length timestep.")

336 if (dt < 0.0) call mom_error(fatal, "MOM_diabatic_driver: "// &

337 "diabatic was called with a negative timestep.")

338

339 showcalltree = calltree_showquery()

340

341 ! Offer diagnostics of various state variables at the start of diabatic

342 ! these are mostly for debugging purposes.

343 if (cs%id_u_predia > 0) call post_data(cs%id_u_predia, u, cs%diag)

344 if (cs%id_v_predia > 0) call post_data(cs%id_v_predia, v, cs%diag)

345 if (cs%id_h_predia > 0) call post_data(cs%id_h_predia, h, cs%diag)

346 if (cs%id_T_predia > 0) call post_data(cs%id_T_predia, tv%T, cs%diag)

347 if (cs%id_S_predia > 0) call post_data(cs%id_S_predia, tv%S, cs%diag)

348 if (cs%id_e_predia > 0) then

349 call find_eta(h, tv, g, gv, us, eta, dzref=g%Z_ref)

350 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 if (stoch_cs%do_sppt) then

355 allocate(h_in(g%isd:g%ied, g%jsd:g%jed, gv%ke)) ; h_in(:,:,:) = h(:,:,:)

356 allocate(t_in(g%isd:g%ied, g%jsd:g%jed, gv%ke)) ; t_in(:,:,:) = tv%T(:,:,:)

357 allocate(s_in(g%isd:g%ied, g%jsd:g%jed, gv%ke)) ; s_in(:,:,:) = tv%S(:,:,:)

358 endif

359

360 if (cs%debug) then

361 call mom_state_chksum("Start of diabatic ", u, v, h, g, gv, us, haloshift=0)

362 call mom_forcing_chksum("Start of diabatic", fluxes, g, us, haloshift=0)

363 endif

364 if (cs%debugConservation) call mom_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, g, gv, us)

365

366 if (cs%debug_energy_req) &

367 call diapyc_energy_req_test(h, dt, tv, g, gv, us, cs%diapyc_en_rec_CSp)

368

369 call cpu_clock_begin(id_clock_set_diffusivity)

370 call set_bbl_tke(u, v, h, tv, fluxes, visc, g, gv, us, cs%set_diff_CSp, obc=obc)

371 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 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 call enable_averages(0.5*dt, time_end - real_to_time(0.5*dt, unscale=us%T_to_s), cs%diag)

379 if (cs%frazil_tendency_diag) then

380 do k=1,nz ; do j=js,je ; do i=is,ie

381 temp_diag(i,j,k) = tv%T(i,j,k)

382 enddo ; enddo ; enddo

383 endif

384

385 if (associated(fluxes%p_surf_full)) then

386 call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp, fluxes%p_surf_full, halo=cs%halo_TS_diff)

387 else

388 call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp, halo=cs%halo_TS_diff)

389 endif

390 if (showcalltree) call calltree_waypoint("done with 1st make_frazil (diabatic)")

391

392 if (cs%frazil_tendency_diag) then

393 call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, g, gv, us, cs)

394 if (cs%id_frazil_h > 0) call post_data(cs%id_frazil_h, h, cs%diag)

395 endif

396 call disable_averaging(cs%diag)

397 endif ! associated(tv%T) .AND. associated(tv%frazil)

398 if (cs%debugConservation) call mom_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, g, gv, us)

399

400 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 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 g, gv, us, cs%int_tide_input_CSp, cs%int_tide_CSp)

407

408 if (showcalltree) call calltree_waypoint("done with propagate_int_tide (diabatic)")

409 endif ! end CS%use_int_tides

410

411 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 g, gv, us, cs, stoch_cs, waves)

414 elseif (cs%useALEalgorithm) then

415 call diabatic_ale(u, v, h, tv, bld, fluxes, visc, adp, cdp, dt, time_end, &

416 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 g, gv, us, cs, waves)

420 endif

421

422 call cpu_clock_begin(id_clock_pass)

423 if (associated(visc%sfc_buoy_flx)) &

424 call pass_var(visc%sfc_buoy_flx, g%domain, halo=1, complete=.not.associated(visc%MLD))

425 if (associated(visc%h_ML)) &

426 call pass_var(visc%h_ML, g%Domain, halo=1, complete=.not.associated(visc%MLD))

427 if (associated(visc%MLD)) &

428 call pass_var(visc%MLD, g%Domain, halo=1, complete=.true.)

429 if (associated(visc%Kv_shear)) &

430 call pass_var(visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)

431 call cpu_clock_end(id_clock_pass)

432

433 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 if (associated(tv%T) .AND. associated(tv%frazil)) then

438 call enable_averages(0.5*dt, time_end, cs%diag)

439 if (cs%frazil_tendency_diag) then

440 do k=1,nz ; do j=js,je ; do i=is,ie

441 temp_diag(i,j,k) = tv%T(i,j,k)

442 enddo ; enddo ; enddo

443 endif

444

445 if (associated(fluxes%p_surf_full)) then

446 call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp, fluxes%p_surf_full)

447 else

448 call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp)

449 endif

450

451 if (cs%frazil_tendency_diag) then

452 call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, g, gv, us, cs)

453 if (cs%id_frazil_h > 0 ) call post_data(cs%id_frazil_h, h, cs%diag)

454 endif

455

456 if (showcalltree) call calltree_waypoint("done with 2nd make_frazil (diabatic)")

457 if (cs%debugConservation) call mom_state_stats('2nd make_frazil', u, v, h, tv%T, tv%S, g, gv, us)

458 call disable_averaging(cs%diag)

459

460 endif ! endif for frazil

461

462 if (stoch_cs%do_sppt) then

463 ! perturb diabatic tendencies.

464 ! These stochastic perturbations do not conserve heat, salt or mass.

465 do k=1,nz ; do j=js,je ; do i=is,ie

466 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 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 h_to_rl2_t2 = gv%H_to_RZ * gv%g_Earth

471 do j=js,je

472 ps(:) = 0.0

473 if (associated(fluxes%p_surf)) then

474 do i=is,ie

475 ps(i) = fluxes%p_surf(i,j)

476 enddo

477 endif

478

479 do i=is,ie

480 pressure(i,1) = ps(i) + (0.5*h_to_rl2_t2)*h(i,j,1)

481 enddo

482 do k=2,nz ; do i=is,ie

483 pressure(i,k) = pressure(i,k-1) + &

484 (0.5*h_to_rl2_t2) * (h(i,j,k) + h(i,j,k-1))

485 enddo ; enddo

486 do k=1,nz

487 call calculate_tfreeze(tv%S(is:ie,j,k), pressure(is:ie,k), t_freeze(is:ie), &

488 tv%eqn_of_state)

489 do i=is,ie

490 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 deallocate(h_in, t_in, s_in)

496 endif

497

498 ! Diagnose mixed layer depths.

499 call enable_averages(dt, time_end, cs%diag)

500 if (cs%id_MLD_003 > 0 .or. cs%id_subMLN2 > 0 .or. cs%id_mlotstsq > 0 .or. cs%MLD_param_003 ) then

501 if (cs%MLD_param_003) then

502 call diagnosemldbydensitydifference(cs%id_MLD_003, h, tv, 0.03*us%kg_m3_to_R, g, gv, us, cs%diag, &

503 cs%ref_h_mld, cs%id_MLD_003_zr, cs%id_MLD_003_rr, &

504 id_n2subml=cs%id_subMLN2, id_mldsq=cs%id_mlotstsq, dz_subml=cs%dz_subML_N2, &

505 mld_out=visc%MLD_param)

506 call convert_mld_to_ml_thickness(visc%MLD_param, h, visc%h_ML_param, tv, g, gv)

507 else

508 call diagnosemldbydensitydifference(cs%id_MLD_003, h, tv, 0.03*us%kg_m3_to_R, g, gv, us, cs%diag, &

509 cs%ref_h_mld, cs%id_MLD_003_zr, cs%id_MLD_003_rr, &

510 id_n2subml=cs%id_subMLN2, id_mldsq=cs%id_mlotstsq, dz_subml=cs%dz_subML_N2)

511 endif

512 endif

513 if (cs%id_MLD_0125 > 0) then

514 call diagnosemldbydensitydifference(cs%id_MLD_0125, h, tv, 0.125*us%kg_m3_to_R, g, gv, us, cs%diag, &

515 ref_h_mld=0.0, id_ref_z=-1, id_ref_rho=-1)

516 endif

517 if (cs%id_MLD_user > 0) then

518 call diagnosemldbydensitydifference(cs%id_MLD_user, h, tv, cs%MLDdensityDifference, g, gv, us, cs%diag, &

519 ref_h_mld=0.0, id_ref_z=-1, id_ref_rho=-1)

520 endif

521 if ((cs%id_MLD_EN1 > 0) .or. (cs%id_MLD_EN2 > 0) .or. (cs%id_MLD_EN3 > 0) .or. (cs%MLD_param_EN1)) then

522 ! Surface Mixed Layer diagnostic

523 if (cs%MLD_param_EN1) then

524 call diagnosemldbyenergy((/cs%id_MLD_EN1, cs%id_MLD_EN2, cs%id_MLD_EN3/), h, tv, g, gv, us, cs%MLD_En_vals, &

525 (/1,nz/), cs%diag, om4_iteration=cs%use_OM4_MLD_En_iter,mld_out=visc%MLD_param)

526 call convert_mld_to_ml_thickness(visc%MLD_param, h, visc%h_ML_param, tv, g, gv)

527 else

528 call diagnosemldbyenergy((/cs%id_MLD_EN1, cs%id_MLD_EN2, cs%id_MLD_EN3/), h, tv, g, gv, us, cs%MLD_En_vals, &

529 (/1,nz/), cs%diag, om4_iteration=cs%use_OM4_MLD_En_iter)

530 endif

531 endif

532 if ((cs%id_BMLD_EN1 > 0) .or. (cs%id_BMLD_EN2 > 0) .or. (cs%id_BMLD_EN3 > 0)) then

533 ! Bottom Mixed Layer diagnostic

534 call diagnosemldbyenergy((/cs%id_BMLD_EN1, cs%id_BMLD_EN2, cs%id_BMLD_EN3/), h, tv, g, gv, us, cs%BMLD_En_vals, &

535 (/nz,1/), cs%diag, om4_iteration=.false.)

536 endif

537 if (stoch_cs%do_sppt .and. stoch_cs%id_sppt_wts > 0) &

538 call post_data(stoch_cs%id_sppt_wts, stoch_cs%sppt_wts, cs%diag)

539

540 call disable_averaging(cs%diag)

541

542 if (cs%debugConservation) call mom_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, g, gv, us)

543

544end subroutine diabatic

545

546

547!> Applies diabatic forcing and diapycnal mixing of temperature, salinity and other tracers for use

548!! with an ALE algorithm. This version uses an older set of algorithms compared with diabatic_ALE.

549subroutine diabatic_ale_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &

550 G, GV, US, CS, stoch_CS, Waves)

551 type(ocean_grid_type), intent(inout) :: G !< ocean grid structure

552 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

553 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

554 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]

555 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]

556 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2]

557 type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields

558 !! unused have NULL ptrs

559 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Active mixed layer depth [Z ~> m]

560 type(forcing), intent(inout) :: fluxes !< points to forcing fields

561 !! unused fields have NULL ptrs

562 type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities,

563 !! BBL properties and related fields

564 type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum

565 !! equations, to enable the later derived

566 !! diagnostics, like energy budgets

567 type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations

568 real, intent(in) :: dt !< time increment [T ~> s]

569 type(time_type), intent(in) :: Time_end !< Time at the end of the interval

570 type(diabatic_cs), pointer :: CS !< module control structure

571 type(stochastic_cs), pointer :: stoch_CS !< stochastic control structure

572 type(wave_parameters_cs), pointer :: Waves !< Surface gravity waves

573

574 ! local variables

575 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &

576 h_orig, & ! Initial layer thicknesses [H ~> m or kg m-2]

577 dz, & ! The vertical distance between interfaces around a layer [Z ~> m]

578 dSV_dT, & ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]

579 dSV_dS, & ! The partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1].

580 cTKE, & ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].

581 u_h, & ! Zonal velocities interpolated to thickness points [L T-1 ~> m s-1]

582 v_h, & ! Meridional velocities interpolated to thickness points [L T-1 ~> m s-1]

583 temp_diag, & ! Diagnostic array of previous temperatures [C ~> degC]

584 saln_diag ! Diagnostic array of previous salinity [S ~> ppt]

585

586 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &

587 ent_s, & ! The diffusive coupling across interfaces within one time step for

588 ! salinity and passive tracers [H ~> m or kg m-2]

589 ent_t, & ! The diffusive coupling across interfaces within one time step for

590 ! temperature [H ~> m or kg m-2]

591 kd_int, & ! diapycnal diffusivity of interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

592 kd_heat, & ! diapycnal diffusivity of heat [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

593 kd_salt, & ! diapycnal diffusivity of salt and passive tracers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

594 kd_extra_t , & ! The extra diffusivity of temperature due to double diffusion relative to

595 ! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

596 kd_extra_s , & ! The extra diffusivity of salinity due to double diffusion relative to

597 ! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

598 kd_epbl, & ! test array of diapycnal diffusivities at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

599 kpp_nltheat, & ! KPP non-local transport for heat [nondim]

600 kpp_nltscalar, & ! KPP non-local transport for scalars [nondim]

601 kpp_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]

602 tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

603 sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

604 n2_salt, & !< Salinity contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]

605 n2_temp !< Temperature contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]

606

607 real, dimension(SZI_(G),SZJ_(G)) :: &

608 U_star, & ! The friction velocity [Z T-1 ~> m s-1].

609 KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

610 KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

611 SkinBuoyFlux, & ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL

612 BBL_BuoyFlux ! 2d bottom buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL

613

614 real, dimension(SZI_(G)) :: &

615 p_i ,& ! Pressure at the interface [R L2 T-2 ~> Pa]

616 T_i, & ! Temperature at the interface [C ~> degC]

617 S_i, & ! Salinity at the interface [S ~> ppt]

618 drhodS, & ! Local change in density w.r.t. salinity using model EOS & state [R C-1 ~> kg m-3 ppt-1]

619 drhodT, & ! Local change in density w.r.t. temperature using model EOS & state [R C-1 ~> kg m-3 degC-1]

620 dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]

621 dSpV_dS ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]

622

623 logical, dimension(SZI_(G)) :: &

624 in_boundary ! True if there are no massive layers below, where massive is defined as

625 ! sufficiently thick that the no-flux boundary conditions have not restricted

626 ! the entrainment - usually sqrt(Kd*dt).

627

628 real :: h_neglect ! A thickness that is so small it is usually lost

629 ! in roundoff and can be neglected [H ~> m or kg m-2]

630 real :: dz_neglect ! A vertical distance that is so small it is usually lost

631 ! in roundoff and can be neglected [Z ~> m]

632 real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2]

633 real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2]

634 real :: I_dzval ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1]

635 real :: I_h ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1]

636 real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is

637 ! coupled to the bottom within a timestep [H ~> m or kg m-2]

638 real :: Kd_add_here ! An added diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1].

639 real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].

640

641 real :: Ent_int ! The diffusive entrainment rate at an interface [H ~> m or kg m-2]

642 real :: Idt ! The inverse time step [T-1 ~> s-1]

643 real :: g_Rho0 ! G_Earth/Rho0 [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2]

644 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]

645 real :: alt_H_to_pres! A conversion factor from thicknesses to pressure w/ alternative scaling [R Z T-2 ~> Pa m-1]

646 logical :: nonBous ! True if not using the Boussinesq approximation

647

648 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state

649

650 logical :: showCallTree ! If true, show the call tree

651 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz

652

653 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

654 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

655 dz_neglect = gv%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect

656 h_neglect = gv%H_subroundoff

657

658 nonbous = .not.(gv%Boussinesq .or. gv%semi_Boussinesq)

659 g_rho0 = gv%g_Earth_Z_T2 / gv%H_to_RZ

660 h_to_pres = gv%H_to_RZ * gv%g_Earth

661 alt_h_to_pres = h_to_pres * us%L_to_Z**2 * gv%Z_to_H

662

663 kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0

664

665 showcalltree = calltree_showquery()

666 if (showcalltree) call calltree_enter("diabatic_ALE_legacy(), MOM_diabatic_driver.F90")

667

668 ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep

669 call enable_averages(dt, time_end, cs%diag)

670

671 if (cs%Use_KdWork_diag) call allocate_vbf_cs(g, gv, cs%VBF)

672

673 if (cs%use_geothermal) then

674 call cpu_clock_begin(id_clock_geothermal)

675 call geothermal_in_place(h, tv, dt, g, gv, us, cs%geothermal, bbl_buoyflux, halo=cs%halo_TS_diff)

676 call cpu_clock_end(id_clock_geothermal)

677 if (showcalltree) call calltree_waypoint("geothermal (diabatic)")

678 if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)

679 endif

680

681 ! Whenever thickness changes let the diag manager know, target grids

682 ! for vertical remapping may need to be regenerated.

683 call diag_update_remap_grids(cs%diag)

684

685 ! Set_pen_shortwave estimates the optical properties of the water column.

686 ! It will need to be modified later to include information about the

687 ! biological properties and layer thicknesses.

688 if (associated(cs%optics)) &

689 call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity, cs%tracer_flow_CSp)

690

691 if (cs%debug) call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)

692

693 if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then

694 if (cs%use_geothermal) then

695 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us, zero_mix=.true.)

696 else

697 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)

698 endif

699 if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")

700 endif

701

702 call cpu_clock_begin(id_clock_set_diffusivity)

703 ! Sets: Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb

704 ! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow

705 if (cs%debug) &

706 call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)

707 if (cs%double_diffuse) then

708 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, kd_int, g, gv, us, &

709 cs%set_diff_CSp, cs%VBF, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)

710 else

711 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, kd_int, g, gv, us, &

712 cs%set_diff_CSp, cs%VBF)

713 endif

714 call cpu_clock_end(id_clock_set_diffusivity)

715 if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")

716

717

718 if (cs%debug) then

719 call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)

720 call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)

721 call mom_thermovar_chksum("after set_diffusivity ", tv, g, us)

722 call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

723 endif

724

725 ! Set diffusivities for heat and salt separately

726 if (cs%useKPP) then

727 ! Add contribution from double diffusion

728 if (cs%double_diffuse) then

729 !$OMP parallel do default(shared)

730 do k=1,nz+1 ; do j=js,je ; do i=is,ie

731 kd_salt(i,j,k) = kd_int(i,j,k) + kd_extra_s(i,j,k)

732 kd_heat(i,j,k) = kd_int(i,j,k) + kd_extra_t(i,j,k)

733 enddo ; enddo ; enddo

734 else

735 !$OMP parallel do default(shared)

736 do k=1,nz+1 ; do j=js,je ; do i=is,ie

737 kd_salt(i,j,k) = kd_int(i,j,k)

738 kd_heat(i,j,k) = kd_int(i,j,k)

739 enddo ; enddo ; enddo

740 endif

741

742 if (cs%debug) then

743 call hchksum(kd_heat, "after set_diffusivity Kd_heat", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

744 call hchksum(kd_salt, "after set_diffusivity Kd_salt", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

745 endif

746

747 call cpu_clock_begin(id_clock_kpp)

748 ! total vertical viscosity in the interior is represented via visc%Kv_shear

749

750 ! NOTE: The following do not require initialization, but their checksums do

751 ! require initialization, and past versions were initialized to zero.

752 kpp_nltheat(:,:,:) = 0.

753 kpp_nltscalar(:,:,:) = 0.

754 kpp_buoy_flux(:,:,:) = 0.

755 kpp_temp_flux(:,:) = 0.

756 kpp_salt_flux(:,:) = 0.

757

758 ! KPP needs the surface buoyancy flux but does not update state variables.

759 ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.

760 ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux

761 ! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second)

762 ! unlike other instances where the fluxes are integrated in time over a time-step.

763 call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &

764 kpp_buoy_flux, kpp_temp_flux, kpp_salt_flux)

765

766 ! Determine the friction velocity, perhaps using the evovling surface density.

767 call find_ustar(fluxes, tv, u_star, g, gv, us)

768

769 ! The KPP scheme calculates boundary layer diffusivities and non-local transport.

770 if ( associated(fluxes%lamult) ) then

771 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &

772 u_star, kpp_buoy_flux, waves=waves, lamult=fluxes%lamult)

773

774 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, tv, u_star, kpp_buoy_flux, kd_heat, &

775 kd_salt, visc%Kv_shear, kpp_nltheat, kpp_nltscalar, waves=waves, lamult=fluxes%lamult)

776 else

777 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &

778 u_star, kpp_buoy_flux, waves=waves)

779

780 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, tv, u_star, kpp_buoy_flux, kd_heat, &

781 kd_salt, visc%Kv_shear, kpp_nltheat, kpp_nltscalar, waves=waves)

782 endif

783

784 call kpp_get_bld(cs%KPP_CSp, bld(:,:), g, us)

785 ! If visc%MLD or visc%h_ML exist, copy KPP's BLD into them with appropriate conversions.

786 if (associated(visc%h_ML)) call convert_mld_to_ml_thickness(bld, h, visc%h_ML, tv, g, gv)

787 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

788 if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = kpp_buoy_flux(:,:,1)

789

790 if (.not.cs%KPPisPassive) then

791 !$OMP parallel do default(shared)

792 do k=1,nz+1 ; do j=js,je ; do i=is,ie

793 kd_int(i,j,k) = min( kd_salt(i,j,k), kd_heat(i,j,k) )

794 enddo ; enddo ; enddo

795 if (cs%double_diffuse) then

796 !$OMP parallel do default(shared)

797 do k=1,nz+1 ; do j=js,je ; do i=is,ie

798 kd_extra_s(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))

799 kd_extra_t(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))

800 enddo ; enddo ; enddo

801 endif

802 endif ! not passive

803

804 if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")

805 if (cs%debug) then

806 call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)

807 call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)

808 call mom_thermovar_chksum("after KPP", tv, g, us)

809 call hchksum(kd_heat, "after KPP Kd_heat", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

810 call hchksum(kd_salt, "after KPP Kd_salt", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

811 call hchksum(kpp_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, &

812 unscale=us%C_to_degC*gv%H_to_m*us%s_to_T)

813 call hchksum(kpp_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, &

814 unscale=us%S_to_ppt*gv%H_to_m*us%s_to_T)

815 call hchksum(kpp_nltheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)

816 call hchksum(kpp_nltscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)

817 endif

818 ! Apply non-local transport of heat and salt

819 ! Changes: tv%T, tv%S

820 call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, kpp_nltheat, kpp_temp_flux, &

821 dt, tv%tr_T, tv%T, tv%C_p)

822 call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, kpp_nltscalar, kpp_salt_flux, &

823 dt, tv%tr_S, tv%S)

824 call cpu_clock_end(id_clock_kpp)

825 if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")

826 if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)

827

828 if (cs%debug) then

829 call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)

830 call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)

831 call mom_thermovar_chksum("after KPP_applyNLT ", tv, g, us)

832 endif

833 endif ! endif for KPP

834

835 ! This is the "old" method for applying differential diffusion.

836 ! Changes: tv%T, tv%S

837 if (cs%double_diffuse .and. associated(tv%T)) then

838

839 call cpu_clock_begin(id_clock_differential_diff)

840 call differential_diffuse_t_s(h, tv%T, tv%S, kd_extra_t, kd_extra_s, tv, dt, g, gv)

841 call cpu_clock_end(id_clock_differential_diff)

842

843 if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")

844 if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g, gv, us)

845

846 ! increment heat and salt diffusivity.

847 ! CS%useKPP==.true. already has extra_T and extra_S included

848 if (.not. cs%useKPP) then

849 !$OMP parallel do default(shared)

850 do k=2,nz ; do j=js,je ; do i=is,ie

851 kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)

852 kd_salt(i,j,k) = kd_salt(i,j,k) + kd_extra_s(i,j,k)

853 enddo ; enddo ; enddo

854 endif

855

856 endif

857

858 ! Calculate vertical mixing due to convection (computed via CVMix)

859 if (cs%use_CVMix_conv) then

860 ! Increment vertical diffusion and viscosity due to convection

861 call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv, bld, kd_int, visc%Kv_shear)

862 endif

863

864 ! Find the vertical distances across layers.

865 call thickness_to_dz(h, tv, dz, g, gv, us)

866

867 ! This block sets ent_t and ent_s from h and Kd_int.

868 do j=js,je ; do i=is,ie

869 ent_s(i,j,1) = 0.0 ; ent_s(i,j,nz+1) = 0.0

870 ent_t(i,j,1) = 0.0 ; ent_t(i,j,nz+1) = 0.0

871 enddo ; enddo

872 !$OMP parallel do default(shared) private(I_dzval)

873 do k=2,nz ; do j=js,je ; do i=is,ie

874 i_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))

875 ent_s(i,j,k) = dt * i_dzval * kd_int(i,j,k)

876 ent_t(i,j,k) = ent_s(i,j,k)

877 enddo ; enddo ; enddo

878 if (showcalltree) call calltree_waypoint("done setting ent_s and ent_t from Kd_int (diabatic)")

879

880 if (cs%debug) then

881 call mom_forcing_chksum("after calc_entrain ", fluxes, g, us, haloshift=0)

882 call mom_thermovar_chksum("after calc_entrain ", tv, g, us)

883 call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, us, haloshift=0)

884 call hchksum(ent_s, "after calc_entrain ent_s", g%HI, haloshift=0, unscale=gv%H_to_MKS)

885 endif

886

887 ! Save fields before boundary forcing is applied for tendency diagnostics

888 do k=1,nz ; do j=js,je ; do i=is,ie

889 h_orig(i,j,k) = h(i,j,k)

890 enddo ; enddo ; enddo

891 if (cs%boundary_forcing_tendency_diag) then

892 do k=1,nz ; do j=js,je ; do i=is,ie

893 temp_diag(i,j,k) = tv%T(i,j,k)

894 saln_diag(i,j,k) = tv%S(i,j,k)

895 enddo ; enddo ; enddo

896 endif

897

898 ! Apply forcing

899 ! Changes made to following fields: h, tv%T and tv%S.

900 call cpu_clock_begin(id_clock_remap)

901

902 if (cs%use_energetic_PBL) then

903

904 skinbuoyflux(:,:) = 0.0

905 call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &

906 optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, cs%evap_CFL_limit, &

907 cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux, mld_h=visc%h_ML_param)

908

909 if (cs%debug) then

910 call hchksum(ent_t, "after applyBoundaryFluxes ent_t", g%HI, haloshift=0, unscale=gv%H_to_mks)

911 call hchksum(ent_s, "after applyBoundaryFluxes ent_s", g%HI, haloshift=0, unscale=gv%H_to_mks)

912 call hchksum(ctke, "after applyBoundaryFluxes cTKE", g%HI, haloshift=0, &

913 unscale=us%RZ3_T3_to_W_m2*us%T_to_s)

914 call hchksum(dsv_dt, "after applyBoundaryFluxes dSV_dT", g%HI, haloshift=0, &

915 unscale=us%kg_m3_to_R*us%degC_to_C)

916 call hchksum(dsv_ds, "after applyBoundaryFluxes dSV_dS", g%HI, haloshift=0, &

917 unscale=us%kg_m3_to_R*us%ppt_to_S)

918 call hchksum(h, "after applyBoundaryFluxes h", g%HI, haloshift=0, unscale=gv%H_to_mks)

919 call hchksum(tv%T, "after applyBoundaryFluxes tv%T", g%HI, haloshift=0, unscale=us%C_to_degC)

920 call hchksum(tv%S, "after applyBoundaryFluxes tv%S", g%HI, haloshift=0, unscale=us%S_to_ppt)

921 call hchksum(skinbuoyflux, "after applyBdryFlux SkinBuoyFlux", g%HI, haloshift=0, &

922 unscale=us%Z_to_m**2*us%s_to_T**3)

923 endif

924

925 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)

926 call energetic_pbl(h, u_h, v_h, tv, fluxes, visc, dt, kd_epbl, g, gv, us, &

927 cs%ePBL, stoch_cs, dsv_dt, dsv_ds, ctke, skinbuoyflux, bbl_buoyflux, waves=waves)

928

929 call energetic_pbl_get_mld(cs%ePBL, bld(:,:), g, us)

930 ! If visc%MLD or visc%h_ML exist, copy ePBL's BLD into them with appropriate conversions.

931 if (associated(visc%h_ML)) call convert_mld_to_ml_thickness(bld, h, visc%h_ML, tv, g, gv)

932 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

933 if (cs%MLD_param_ePBL) visc%h_ML_param = visc%h_ML

934 if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = skinbuoyflux(:,:)

935

936 ! Find the vertical distances across layers, which may have been modified by the net surface flux

937 call thickness_to_dz(h, tv, dz, g, gv, us)

938

939 ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.

940 do k=2,nz ; do j=js,je ; do i=is,ie

941 if (cs%ePBL_is_additive) then

942 kd_add_here = kd_epbl(i,j,k)

943 visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%ePBL_Prandtl*kd_epbl(i,j,k)

944 else

945 kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)

946 visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), cs%ePBL_Prandtl*kd_epbl(i,j,k))

947 endif

948

949 ent_int = kd_add_here * dt / (0.5*(dz(i,j,k-1) + dz(i,j,k)) + dz_neglect)

950 ent_s(i,j,k) = ent_s(i,j,k) + ent_int

951 kd_int(i,j,k) = kd_int(i,j,k) + kd_add_here

952

953 ! for diagnostics

954 kd_heat(i,j,k) = kd_heat(i,j,k) + kd_int(i,j,k)

955 kd_salt(i,j,k) = kd_salt(i,j,k) + kd_int(i,j,k)

956

957 enddo ; enddo ; enddo

958

959 if (cs%debug) then

960 call hchksum(ent_t, "after ePBL ent_t", g%HI, haloshift=0, unscale=gv%H_to_MKS)

961 call hchksum(ent_s, "after ePBL ent_s", g%HI, haloshift=0, unscale=gv%H_to_MKS)

962 call hchksum(kd_epbl, "after ePBL Kd_ePBL", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

963 endif

964

965 else

966 call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &

967 optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, &

968 cs%evap_CFL_limit, cs%minimum_forcing_depth, mld_h=visc%h_ML_param)

969

970 ! Find the vertical distances across layers, which may have been modified by the net surface flux

971 call thickness_to_dz(h, tv, dz, g, gv, us)

972

973 endif ! endif for CS%use_energetic_PBL

974

975 ! diagnose the tendencies due to boundary forcing

976 ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme

977 ! so all tendency diagnostics need to be posted on h_orig, and grids rebuilt afterwards

978 if (cs%boundary_forcing_tendency_diag) then

979 call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_orig, dt, g, gv, us, cs)

980 if (cs%id_boundary_forcing_h > 0) call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h=h_orig)

981 endif

982 ! Boundary fluxes may have changed T, S, and h

983 call diag_update_remap_grids(cs%diag)

984 call cpu_clock_end(id_clock_remap)

985 if (cs%debug) then

986 call mom_forcing_chksum("after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)

987 call mom_thermovar_chksum("after applyBoundaryFluxes ", tv, g, us)

988 call mom_state_chksum("after applyBoundaryFluxes ", u, v, h, g, gv, us, haloshift=0)

989 endif

990 if (showcalltree) call calltree_waypoint("done with applyBoundaryFluxes (diabatic)")

991 if (cs%debugConservation) call mom_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, g, gv, us)

992

993 if (cs%debug) then

994 call mom_state_chksum("after negative check ", u, v, h, g, gv, us, haloshift=0)

995 call mom_forcing_chksum("after negative check ", fluxes, g, us, haloshift=0)

996 call mom_thermovar_chksum("after negative check ", tv, g, us)

997 endif

998 if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")

999 if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)

1000

1001 ! calculate change in temperature & salinity due to dia-coordinate surface diffusion

1002 if (associated(tv%T)) then

1003

1004 if (cs%debug) then

1005 call hchksum(ent_t, "before triDiagTS ent_t ", g%HI, haloshift=0, unscale=gv%H_to_MKS)

1006 call hchksum(ent_s, "before triDiagTS ent_s ", g%HI, haloshift=0, unscale=gv%H_to_MKS)

1007 endif

1008

1009 call cpu_clock_begin(id_clock_tridiag)

1010 ! Keep salinity from falling below a small but positive threshold.

1011 ! This constraint is needed for SIS1 ice model, which can extract

1012 ! more salt than is present in the ocean. SIS2 does not suffer

1013 ! from this limitation, in which case we can let salinity=0 and still

1014 ! have salt conserved with SIS2 ice. So for SIS2, we can run with

1015 ! BOUND_SALINITY=False in MOM.F90.

1016 if (associated(tv%S) .and. associated(tv%salt_deficit)) &

1017 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

1018

1019 if (cs%diabatic_diff_tendency_diag) then

1020 do k=1,nz ; do j=js,je ; do i=is,ie

1021 temp_diag(i,j,k) = tv%T(i,j,k)

1022 saln_diag(i,j,k) = tv%S(i,j,k)

1023 enddo ; enddo ; enddo

1024 endif

1025

1026 ! Changes T and S via the tridiagonal solver; no change to h

1027 do k=1,nz ; do j=js,je ; do i=is,ie

1028 ent_t(i,j,k) = ent_s(i,j,k) ; ent_t(i,j,k+1) = ent_s(i,j,k+1)

1029 enddo ; enddo ; enddo

1030 if (cs%tracer_tridiag) then

1031 call tracer_vertdiff_eulerian(h, ent_t, dt, tv%T, g, gv)

1032 call tracer_vertdiff_eulerian(h, ent_s, dt, tv%S, g, gv)

1033 else

1034 call tridiagts_eulerian(g, gv, is, ie, js, je, h, ent_s, tv%T, tv%S)

1035 endif

1036

1037 ! diagnose temperature, salinity, heat, and salt tendencies

1038 if (cs%diabatic_diff_tendency_diag) then

1039 call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, g, gv, us, cs)

1040 if (cs%id_diabatic_diff_h > 0) call post_data(cs%id_diabatic_diff_h, h, cs%diag, alt_h=h)

1041 endif

1042

1043 call cpu_clock_end(id_clock_tridiag)

1044

1045 if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")

1046

1047 endif ! endif corresponding to if (associated(tv%T))

1048

1049 if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)

1050

1051 if (cs%debug) then

1052 call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)

1053 call mom_thermovar_chksum("after mixed layer ", tv, g, us)

1054 endif

1055

1056 ! Whenever thickness changes let the diag manager know, as the

1057 ! target grids for vertical remapping may need to be regenerated.

1058 call diag_update_remap_grids(cs%diag)

1059

1060 ! Set diffusivities for VBF diagnostics if enabled

1061 if (cs%use_energetic_PBL .and. associated(cs%VBF%Kd_ePBL)) cs%VBF%Kd_ePBL(:,:,:) = kd_epbl(:,:,:)

1062 if (associated(cs%VBF%Kd_salt)) cs%VBF%Kd_temp(:,:,:) = kd_heat(:,:,:)

1063 if (associated(cs%VBF%Kd_temp)) cs%VBF%Kd_salt(:,:,:) = kd_salt(:,:,:)

1064

1065

1066 ! Diagnose the diapycnal diffusivities and other related quantities.

1067 if (cs%id_Kd_int > 0) call post_data(cs%id_Kd_int, kd_int, cs%diag)

1068 if (cs%id_Kd_heat > 0) call post_data(cs%id_Kd_heat, kd_heat, cs%diag)

1069 if (cs%id_Kd_salt > 0) call post_data(cs%id_Kd_salt, kd_salt, cs%diag)

1070 if (cs%id_Kd_ePBL > 0) call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)

1071

1072 if (cs%id_ea > 0) call post_data(cs%id_ea, ent_s(:,:,1:nz), cs%diag)

1073 if (cs%id_eb > 0) call post_data(cs%id_eb, ent_s(:,:,2:nz+1), cs%diag)

1074 if (cs%id_ea_t > 0) call post_data(cs%id_ea_t, ent_t(:,:,1:nz), cs%diag)

1075 if (cs%id_eb_t > 0) call post_data(cs%id_eb_t, ent_t(:,:,2:nz+1), cs%diag)

1076 if (cs%id_ea_s > 0) call post_data(cs%id_ea_s, ent_s(:,:,1:nz), cs%diag)

1077 if (cs%id_eb_s > 0) call post_data(cs%id_eb_s, ent_s(:,:,2:nz+1), cs%diag)

1078 idt = 1.0 / dt

1079 if (cs%id_Tdif > 0) then

1080 do j=js,je ; do i=is,ie

1081 tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0

1082 enddo ; enddo

1083 !$OMP parallel do default(shared)

1084 do k=2,nz ; do j=js,je ; do i=is,ie

1085 tdif_flx(i,j,k) = (idt * ent_t(i,j,k)) * (tv%T(i,j,k-1) - tv%T(i,j,k))

1086 enddo ; enddo ; enddo

1087 if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)

1088 endif

1089 if (cs%id_Sdif > 0) then

1090 do j=js,je ; do i=is,ie

1091 sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0

1092 enddo ; enddo

1093 !$OMP parallel do default(shared)

1094 do k=2,nz ; do j=js,je ; do i=is,ie

1095 sdif_flx(i,j,k) = (idt * ent_s(i,j,k)) * (tv%S(i,j,k-1) - tv%S(i,j,k))

1096 enddo ; enddo ; enddo

1097 if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)

1098 endif

1099

1100 if (cs%Use_KdWork_diag .or. cs%Use_N2_diag) then

1101 n2_salt(:,:,:) = 0.0

1102 n2_temp(:,:,:) = 0.0

1103 !Compute N2 and don't mask negatives here

1104 eosdom(:) = eos_domain(g%HI)

1105 if (nonbous) then

1106 !$OMP parallel do default(shared)

1107 do j=js,je

1108 if (associated(tv%p_surf)) then

1109 do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo

1110 else

1111 do i=is,ie ; p_i(i) = 0.0 ; enddo

1112 endif

1113 do k=2,nz

1114 do i=is,ie

1115 p_i(i) = p_i(i) + h_to_pres * h(i,j,k-1)

1116 enddo

1117 t_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))

1118 s_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))

1119 call calculate_specific_vol_derivs(t_i, s_i, p_i, dspv_dt, dspv_ds, tv%eqn_of_state, eosdom)

1120 do i=is,ie

1121 i_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))

1122 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))

1123 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))

1124 enddo

1125 enddo

1126 enddo

1127 else

1128 !$OMP parallel do default(shared)

1129 do j=js,je

1130 if (associated(tv%p_surf)) then

1131 do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo

1132 else

1133 do i=is,ie ; p_i(i) = 0.0 ; enddo

1134 endif

1135 do k=2,nz

1136 do i=is,ie

1137 p_i(i) = p_i(i) + h_to_pres* h(i,j,k-1)

1138 enddo

1139 t_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))

1140 s_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))

1141 call calculate_density_derivs(t_i, s_i, p_i, drhodt, drhods, tv%eqn_of_state, eosdom)

1142 do i=is,ie

1143 i_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))

1144 n2_salt(i,j,k) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (drhods(i) * (g_rho0 * i_h))

1145 n2_temp(i,j,k) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (drhodt(i) * (g_rho0 * i_h))

1146 enddo

1147 enddo

1148 enddo

1149 endif

1150 if (cs%id_N2_dd>0) call post_data(cs%id_N2_dd, n2_salt(:,:,:)+n2_temp(:,:,:), cs%diag)

1151 if (cs%id_N2_salt_dd>0) call post_data(cs%id_N2_salt_dd, n2_salt, cs%diag)

1152 if (cs%id_N2_temp_dd>0) call post_data(cs%id_N2_temp_dd, n2_temp, cs%diag)

1153

1154 if (cs%Use_KdWork_diag) then

1155 call kdwork_diagnostics(g,gv,us,cs%diag,cs%VBF,n2_salt,n2_temp,dz)

1156 endif

1157

1158 call deallocate_vbf_cs(cs%VBF)

1159

1160 endif

1161

1162 ! mixing of passive tracers from massless boundary layers to interior

1163 call cpu_clock_begin(id_clock_tracers)

1164

1165 if (cs%mix_boundary_tracer_ALE) then

1166 tr_ea_bbl = sqrt(dt * cs%Kd_BBL_tr)

1167

1168 !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)

1169 do j=js,je

1170 do i=is,ie

1171 htot(i) = 0.0

1172 in_boundary(i) = (g%mask2dT(i,j) > 0.0)

1173 enddo

1174 do k=nz,2,-1 ; do i=is,ie

1175 if (in_boundary(i)) then

1176 htot(i) = htot(i) + h(i,j,k)

1177 ! If diapycnal mixing has been suppressed because this is a massless

1178 ! layer near the bottom, add some mixing of tracers between these

1179 ! layers. This flux is based on the harmonic mean of the two

1180 ! thicknesses, as this corresponds pretty closely (to within

1181 ! differences in the density jumps between layers) with what is done

1182 ! in the calculation of the fluxes in the first place. Kd_min_tr

1183 ! should be much less than the values that have been set in Kd_int,

1184 ! perhaps a molecular diffusivity.

1185 add_ent = ((dt * cs%Kd_min_tr)) * &

1186 ((dz(i,j,k-1)+dz(i,j,k)+dz_neglect) / (dz(i,j,k-1)*dz(i,j,k)+dz_neglect2)) - &

1187 0.5*(ent_s(i,j,k) + ent_s(i,j,k))

1188 if (htot(i) < tr_ea_bbl) then

1189 add_ent = max(0.0, add_ent, (tr_ea_bbl - htot(i)) - ent_s(i,j,k))

1190 elseif (add_ent < 0.0) then

1191 add_ent = 0.0 ; in_boundary(i) = .false.

1192 endif

1193

1194 ent_s(i,j,k) = ent_s(i,j,k) + add_ent

1195 endif

1196

1197 if (cs%double_diffuse) then ; if (kd_extra_s(i,j,k) > 0.0) then

1198 add_ent = (dt * kd_extra_s(i,j,k)) / &

1199 (0.5 * (dz(i,j,k-1) + dz(i,j,k)) + dz_neglect)

1200 ent_s(i,j,k) = ent_s(i,j,k) + add_ent

1201 endif ; endif

1202 enddo ; enddo

1203

1204 enddo

1205 elseif (cs%double_diffuse .and. .not.cs%mix_boundary_tracers) then ! extra diffusivity for passive tracers

1206 !$OMP parallel do default(shared) private(add_ent)

1207 do k=nz,2,-1 ; do j=js,je ; do i=is,ie

1208 if (kd_extra_s(i,j,k) > 0.0) then

1209 add_ent = (dt * kd_extra_s(i,j,k)) / &

1210 (0.5 * (dz(i,j,k-1) + dz(i,j,k)) + dz_neglect)

1211 else

1212 add_ent = 0.0

1213 endif

1214 ent_s(i,j,k) = ent_s(i,j,k) + add_ent

1215 enddo ; enddo ; enddo

1216 endif ! (CS%mix_boundary_tracers)

1217

1218 ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied

1219 call call_tracer_column_fns(h_orig, h, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), fluxes, bld, dt, &

1220 g, gv, us, tv, cs%optics, cs%tracer_flow_CSp, cs%debug, &

1221 kpp_csp=cs%KPP_CSp, &

1222 nonlocaltrans=kpp_nltscalar, &

1223 evap_cfl_limit=cs%evap_CFL_limit, &

1224 minimum_forcing_depth=cs%minimum_forcing_depth, h_bl=visc%h_ML)

1225

1226 call cpu_clock_end(id_clock_tracers)

1227

1228 ! Apply ALE sponge

1229 if (cs%use_sponge .and. associated(cs%ALE_sponge_CSp)) then

1230 call cpu_clock_begin(id_clock_sponge)

1231 call apply_ale_sponge(h, tv, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)

1232 call cpu_clock_end(id_clock_sponge)

1233 if (cs%debug) then

1234 call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)

1235 call mom_thermovar_chksum("apply_sponge ", tv, g, us)

1236 endif

1237 endif ! CS%use_sponge

1238

1239 ! Apply data assimilation incremental update -oda_incupd-

1240 if (cs%use_oda_incupd .and. associated(cs%oda_incupd_CSp)) then

1241 call mom_mesg("Starting ODA_INCUPD legacy ", 5)

1242 call cpu_clock_begin(id_clock_oda_incupd)

1243 call apply_oda_incupd(h, tv, u, v, dt, g, gv, us, cs%oda_incupd_CSp)

1244 call cpu_clock_end(id_clock_oda_incupd)

1245 if (cs%debug) then

1246 call mom_state_chksum("apply_oda_incupd ", u, v, h, g, gv, us, haloshift=0)

1247 call mom_thermovar_chksum("apply_oda_incupd ", tv, g, us)

1248 endif

1249 endif ! CS%use_oda_incupd

1250

1251

1252

1253 call disable_averaging(cs%diag)

1254

1255 if (showcalltree) call calltree_leave("diabatic_ALE_legacy()")

1256

1257end subroutine diabatic_ale_legacy

1258

1259

1260!> This subroutine imposes the diapycnal mass fluxes and the

1261!! accompanying diapycnal advection of momentum and tracers.

1262subroutine diabatic_ale(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &

1263 G, GV, US, CS, stoch_CS, Waves)

1264 type(ocean_grid_type), intent(inout) :: G !< ocean grid structure

1265 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

1266 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

1267 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]

1268 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]

1269 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2]

1270 type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields

1271 !! unused have NULL ptrs

1272 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Active mixed layer depth [Z ~> m]

1273 type(forcing), intent(inout) :: fluxes !< points to forcing fields

1274 !! unused fields have NULL ptrs

1275 type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities,

1276 !! BBL properties and related fields

1277 type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum

1278 !! equations, to enable the later derived

1279 !! diagnostics, like energy budgets

1280 type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations

1281 real, intent(in) :: dt !< time increment [T ~> s]

1282 type(time_type), intent(in) :: Time_end !< Time at the end of the interval

1283 type(diabatic_cs), pointer :: CS !< module control structure

1284 type(stochastic_cs), pointer :: stoch_CS !< stochastic control structure

1285 type(wave_parameters_cs), pointer :: Waves !< Surface gravity waves

1286

1287 ! local variables

1288 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &

1289 h_orig, & ! Initial layer thicknesses [H ~> m or kg m-2]

1290 dz, & ! The vertical distance between interfaces around a layer [Z ~> m]

1291 dSV_dT, & ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]

1292 dSV_dS, & ! The partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1].

1293 cTKE, & ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].

1294 u_h, & ! Zonal velocities interpolated to thickness points [L T-1 ~> m s-1]

1295 v_h, & ! Meridional velocities interpolated to thickness points [L T-1 ~> m s-1]

1296 temp_diag, & ! Diagnostic array of previous temperatures [C ~> degC]

1297 saln_diag ! Diagnostic array of previous salinity [S ~> ppt]

1298

1299 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &

1300 ent_s, & ! The diffusive coupling across interfaces within one time step for

1301 ! salinity and passive tracers [H ~> m or kg m-2]

1302 ent_t, & ! The diffusive coupling across interfaces within one time step for

1303 ! temperature [H ~> m or kg m-2]

1304 kd_heat, & ! diapycnal diffusivity of heat or the smaller of the diapycnal diffusivities of

1305 ! heat and salt [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1306 kd_salt, & ! diapycnal diffusivity of salt and passive tracers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1307 kd_extra_t , & ! The extra diffusivity of temperature due to double diffusion relative to

1308 ! Kd_int returned from set_diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1309 kd_extra_s , & ! The extra diffusivity of salinity due to double diffusion relative to

1310 ! Kd_int returned from set_diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1311 kd_epbl, & ! boundary layer or convective diapycnal diffusivities at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1312 kpp_nltheat, & ! KPP non-local transport for heat [nondim]

1313 kpp_nltscalar, & ! KPP non-local transport for scalars [nondim]

1314 kpp_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]

1315 tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

1316 sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

1317 n2_salt, & !< Salinity contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]

1318 n2_temp !< Temperature contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]

1319

1320 real, dimension(SZI_(G),SZJ_(G)) :: &

1321 U_star, & ! The friction velocity [Z T-1 ~> m s-1].

1322 KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

1323 KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

1324 SkinBuoyFlux, & ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL

1325 BBL_BuoyFlux ! 2d bottom buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL

1326

1327 logical, dimension(SZI_(G)) :: &

1328 in_boundary ! True if there are no massive layers below, where massive is defined as

1329 ! sufficiently thick that the no-flux boundary conditions have not restricted

1330 ! the entrainment - usually sqrt(Kd*dt).

1331

1332 real, dimension(SZI_(G)) :: &

1333 p_i ,& ! Pressure at the interface [R L2 T-2 ~> Pa]

1334 T_i, & ! Temperature at the interface [C ~> degC]

1335 S_i, & ! Salinity at the interface [S ~> ppt]

1336 drhodS, & ! Local change in density w.r.t. salinity using model EOS & state [R C-1 ~> kg m-3 ppt-1]

1337 drhodT, & ! Local change in density w.r.t. temperature using model EOS & state [R C-1 ~> kg m-3 degC-1]

1338 dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]

1339 dSpV_dS ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]

1340

1341 real :: h_neglect ! A thickness that is so small it is usually lost

1342 ! in roundoff and can be neglected [H ~> m or kg m-2]

1343 real :: dz_neglect ! A vertical distance that is so small it is usually lost

1344 ! in roundoff and can be neglected [Z ~> m]

1345 real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2]

1346 real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2]

1347 real :: I_dzval ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1]

1348 real :: I_h ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1]

1349 real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is

1350 ! coupled to the bottom within a timestep [H ~> m or kg m-2]

1351 real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].

1352 real :: Kd_add_here ! An added diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1].

1353 real :: Idt ! The inverse time step [T-1 ~> s-1]

1354 real :: g_Rho0 ! G_Earth/Rho0 [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2]

1355 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]

1356 real :: alt_H_to_pres! A conversion factor from thicknesses to pressure w/ alternative scaling [R Z T-2 ~> Pa m-1]

1357 logical :: nonBous ! True if not using the Boussinesq approximation

1358

1359 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state

1360

1361 logical :: showCallTree ! If true, show the call tree

1362 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, Isq, Ieq, Jsq, Jeq, nz

1363

1364 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

1365 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

1366 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

1367 dz_neglect = gv%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect

1368 h_neglect = gv%H_subroundoff

1369

1370 nonbous = .not.(gv%Boussinesq .or. gv%semi_Boussinesq)

1371 g_rho0 = gv%g_Earth_Z_T2 / gv%H_to_RZ

1372 h_to_pres = gv%H_to_RZ * gv%g_Earth

1373 alt_h_to_pres = h_to_pres * us%L_to_Z**2 * gv%Z_to_H

1374

1375 kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0

1376 ent_s(:,:,:) = 0.0 ; ent_t(:,:,:) = 0.0

1377

1378 showcalltree = calltree_showquery()

1379 if (showcalltree) call calltree_enter("diabatic_ALE(), MOM_diabatic_driver.F90")

1380

1381 if (.not. (cs%useALEalgorithm)) call mom_error(fatal, "MOM_diabatic_driver: "// &

1382 "The ALE algorithm must be enabled when using MOM_diabatic_driver.")

1383

1384 ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep

1385 call enable_averages(dt, time_end, cs%diag)

1386

1387 if (cs%Use_KdWork_diag) call allocate_vbf_cs(g, gv, cs%VBF)

1388

1389 if (cs%use_geothermal) then

1390 call cpu_clock_begin(id_clock_geothermal)

1391 call geothermal_in_place(h, tv, dt, g, gv, us, cs%geothermal, bbl_buoyflux, halo=cs%halo_TS_diff)

1392 call cpu_clock_end(id_clock_geothermal)

1393 if (showcalltree) call calltree_waypoint("geothermal (diabatic)")

1394 if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)

1395 endif

1396

1397 ! Whenever thickness changes let the diag manager know, target grids

1398 ! for vertical remapping may need to be regenerated.

1399 call diag_update_remap_grids(cs%diag)

1400

1401 ! Set_pen_shortwave estimates the optical properties of the water column.

1402 ! It will need to be modified later to include information about the

1403 ! biological properties and layer thicknesses.

1404 if (associated(cs%optics)) &

1405 call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity, cs%tracer_flow_CSp)

1406

1407 if (cs%debug) call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)

1408

1409 if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then

1410 if (cs%use_geothermal) then

1411 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us, zero_mix=.true.)

1412 else

1413 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)

1414 endif

1415 if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")

1416 endif

1417

1418 call cpu_clock_begin(id_clock_set_diffusivity)

1419 ! Sets: Kd_heat, Kd_extra_T, Kd_extra_S and visc%TKE_turb

1420 ! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow

1421 if (cs%debug) &

1422 call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)

1423 if (cs%double_diffuse) then

1424 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, kd_heat, g, gv, us, &

1425 cs%set_diff_CSp, cs%VBF, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)

1426 else

1427 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, kd_heat, g, gv, us, &

1428 cs%set_diff_CSp, cs%VBF)

1429 endif

1430 call cpu_clock_end(id_clock_set_diffusivity)

1431 if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")

1432

1433 if (cs%debug) then

1434 call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)

1435 call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)

1436 call mom_thermovar_chksum("after set_diffusivity ", tv, g, us)

1437 call hchksum(kd_heat, "after set_diffusivity Kd_heat", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

1438 endif

1439

1440 ! Set diffusivities for heat and salt separately, and possibly change the meaning of Kd_heat.

1441 if (cs%double_diffuse) then

1442 ! Add contributions from double diffusion

1443 !$OMP parallel do default(shared)

1444 do k=1,nz+1 ; do j=js,je ; do i=is,ie

1445 kd_salt(i,j,k) = kd_heat(i,j,k) + kd_extra_s(i,j,k)

1446 kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)

1447 enddo ; enddo ; enddo

1448 else

1449 !$OMP parallel do default(shared)

1450 do k=1,nz+1 ; do j=js,je ; do i=is,ie

1451 kd_salt(i,j,k) = kd_heat(i,j,k)

1452 enddo ; enddo ; enddo

1453 endif

1454

1455 if (cs%debug) then

1456 call hchksum(kd_heat, "after double diffuse Kd_heat", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

1457 call hchksum(kd_salt, "after double diffuse Kd_salt", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

1458 endif

1459

1460 if (cs%useKPP) then

1461 call cpu_clock_begin(id_clock_kpp)

1462 ! total vertical viscosity in the interior is represented via visc%Kv_shear

1463 do k=1,nz+1 ; do j=js,je ; do i=is,ie

1464 visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k)

1465 enddo ; enddo ; enddo

1466

1467 ! NOTE: The following do not require initialization, but their checksums do

1468 ! require initialization, and past versions were initialized to zero.

1469 kpp_nltheat(:,:,:) = 0.

1470 kpp_nltscalar(:,:,:) = 0.

1471 kpp_buoy_flux(:,:,:) = 0.

1472 kpp_temp_flux(:,:) = 0.

1473 kpp_salt_flux(:,:) = 0.

1474

1475 ! KPP needs the surface buoyancy flux but does not update state variables.

1476 ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.

1477 ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux

1478 ! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second)

1479 ! unlike other instances where the fluxes are integrated in time over a time-step.

1480 call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &

1481 kpp_buoy_flux, kpp_temp_flux, kpp_salt_flux)

1482

1483 ! Determine the friction velocity, perhaps using the evovling surface density.

1484 call find_ustar(fluxes, tv, u_star, g, gv, us)

1485

1486 ! The KPP scheme calculates boundary layer diffusivities and non-local transport.

1487 if ( associated(fluxes%lamult) ) then

1488 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &

1489 u_star, kpp_buoy_flux, waves=waves, lamult=fluxes%lamult)

1490

1491 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, tv, u_star, kpp_buoy_flux, kd_heat, &

1492 kd_salt, visc%Kv_shear, kpp_nltheat, kpp_nltscalar, waves=waves, lamult=fluxes%lamult)

1493 else

1494 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &

1495 u_star, kpp_buoy_flux, waves=waves)

1496

1497 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, tv, u_star, kpp_buoy_flux, kd_heat, &

1498 kd_salt, visc%Kv_shear, kpp_nltheat, kpp_nltscalar, waves=waves)

1499 endif

1500

1501 call kpp_get_bld(cs%KPP_CSp, bld(:,:), g, us)

1502 ! If visc%MLD or visc%h_ML exist, copy KPP's BLD into them with appropriate conversions.

1503 if (associated(visc%h_ML)) call convert_mld_to_ml_thickness(bld, h, visc%h_ML, tv, g, gv)

1504 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

1505 if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = kpp_buoy_flux(:,:,1)

1506

1507 if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")

1508 if (cs%debug) then

1509 call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)

1510 call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)

1511 call mom_thermovar_chksum("after KPP", tv, g, us)

1512 call hchksum(kd_heat, "after KPP Kd_heat", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

1513 call hchksum(kd_salt, "after KPP Kd_salt", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

1514 call hchksum(kpp_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, &

1515 unscale=us%C_to_degC*gv%H_to_m*us%s_to_T)

1516 call hchksum(kpp_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, &

1517 unscale=us%S_to_ppt*gv%H_to_m*us%s_to_T)

1518 call hchksum(kpp_nltheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)

1519 call hchksum(kpp_nltscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)

1520 endif

1521 ! Apply non-local transport of heat and salt

1522 ! Changes: tv%T, tv%S

1523 call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, kpp_nltheat, kpp_temp_flux, &

1524 dt, tv%tr_T, tv%T, tv%C_p)

1525 call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, kpp_nltscalar, kpp_salt_flux, &

1526 dt, tv%tr_S, tv%S)

1527 call cpu_clock_end(id_clock_kpp)

1528 if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")

1529 if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)

1530

1531 if (cs%debug) then

1532 call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)

1533 call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)

1534 call mom_thermovar_chksum("after KPP_applyNLT ", tv, g, us)

1535 endif

1536 endif ! endif for KPP

1537

1538 ! Calculate vertical mixing due to convection (computed via CVMix)

1539 if (cs%use_CVMix_conv) then

1540 ! Increment vertical diffusion and viscosity due to convection

1541 call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv, bld, kd_heat, visc%Kv_shear, kd_aux=kd_salt)

1542 endif

1543

1544 ! Save fields before boundary forcing is applied for tendency diagnostics

1545 do k=1,nz ; do j=js,je ; do i=is,ie

1546 h_orig(i,j,k) = h(i,j,k)

1547 enddo ; enddo ; enddo

1548 if (cs%boundary_forcing_tendency_diag) then

1549 do k=1,nz ; do j=js,je ; do i=is,ie

1550 temp_diag(i,j,k) = tv%T(i,j,k)

1551 saln_diag(i,j,k) = tv%S(i,j,k)

1552 enddo ; enddo ; enddo

1553 endif

1554

1555 ! Apply forcing

1556 ! Changes made to following fields: h, tv%T and tv%S.

1557 call cpu_clock_begin(id_clock_remap)

1558

1559 if (cs%use_energetic_PBL) then

1560

1561 skinbuoyflux(:,:) = 0.0

1562 call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &

1563 optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, cs%evap_CFL_limit, &

1564 cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux, mld_h=visc%h_ML_param)

1565

1566 if (cs%debug) then

1567 call hchksum(ent_t, "after applyBoundaryFluxes ent_t", g%HI, haloshift=0, unscale=gv%H_to_MKS)

1568 call hchksum(ent_s, "after applyBoundaryFluxes ent_s", g%HI, haloshift=0, unscale=gv%H_to_MKS)

1569 call hchksum(ctke, "after applyBoundaryFluxes cTKE", g%HI, haloshift=0, &

1570 unscale=us%RZ3_T3_to_W_m2*us%T_to_s)

1571 call hchksum(dsv_dt, "after applyBoundaryFluxes dSV_dT", g%HI, haloshift=0, &

1572 unscale=us%kg_m3_to_R*us%degC_to_C)

1573 call hchksum(dsv_ds, "after applyBoundaryFluxes dSV_dS", g%HI, haloshift=0, &

1574 unscale=us%kg_m3_to_R*us%ppt_to_S)

1575 endif

1576

1577 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)

1578 call energetic_pbl(h, u_h, v_h, tv, fluxes, visc, dt, kd_epbl, g, gv, us, &

1579 cs%ePBL, stoch_cs, dsv_dt, dsv_ds, ctke, skinbuoyflux, bbl_buoyflux, waves=waves)

1580

1581 call energetic_pbl_get_mld(cs%ePBL, bld(:,:), g, us)

1582 ! If visc%MLD or visc%h_ML exist, copy ePBL's BLD into them with appropriate conversions.

1583 if (associated(visc%h_ML)) call convert_mld_to_ml_thickness(bld, h, visc%h_ML, tv, g, gv)

1584 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

1585 if (cs%MLD_param_ePBL) visc%h_ML_param = visc%h_ML

1586 if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = skinbuoyflux(:,:)

1587

1588 ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.

1589 do k=2,nz ; do j=js,je ; do i=is,ie

1590 if (cs%ePBL_is_additive) then

1591 kd_add_here = kd_epbl(i,j,k)

1592 visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%ePBL_Prandtl*kd_epbl(i,j,k)

1593 else

1594 kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)

1595 visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), cs%ePBL_Prandtl*kd_epbl(i,j,k))

1596 endif

1597

1598 kd_heat(i,j,k) = kd_heat(i,j,k) + kd_add_here

1599 kd_salt(i,j,k) = kd_salt(i,j,k) + kd_add_here

1600 enddo ; enddo ; enddo

1601

1602 if (cs%debug) then

1603 call hchksum(ent_t, "after ePBL ent_t", g%HI, haloshift=0, unscale=gv%H_to_MKS)

1604 call hchksum(ent_s, "after ePBL ent_s", g%HI, haloshift=0, unscale=gv%H_to_MKS)

1605 call hchksum(kd_epbl, "after ePBL Kd_ePBL", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

1606 endif

1607

1608 else

1609 call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &

1610 optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, &

1611 cs%evap_CFL_limit, cs%minimum_forcing_depth, mld_h=visc%h_ML_param)

1612

1613 endif ! endif for CS%use_energetic_PBL

1614

1615 ! diagnose the tendencies due to boundary forcing

1616 ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme

1617 ! so all tendency diagnostics need to be posted on h_orig, and grids rebuilt afterwards

1618 if (cs%boundary_forcing_tendency_diag) then

1619 call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_orig, dt, g, gv, us, cs)

1620 if (cs%id_boundary_forcing_h > 0) call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h=h_orig)

1621 endif

1622 ! Boundary fluxes may have changed T, S, and h

1623 call diag_update_remap_grids(cs%diag)

1624 call cpu_clock_end(id_clock_remap)

1625 if (cs%debug) then

1626 call mom_forcing_chksum("after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)

1627 call mom_thermovar_chksum("after applyBoundaryFluxes ", tv, g, us)

1628 call mom_state_chksum("after applyBoundaryFluxes ", u, v, h, g, gv, us, haloshift=0)

1629 endif

1630 if (showcalltree) call calltree_waypoint("done with applyBoundaryFluxes (diabatic)")

1631 if (cs%debugConservation) call mom_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, g, gv, us)

1632

1633 if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")

1634 if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)

1635

1636 ! calculate change in temperature & salinity due to dia-coordinate surface diffusion

1637 if (associated(tv%T)) then

1638

1639 if (cs%debug) then

1640 call hchksum(ent_t, "before triDiagTS ent_t ", g%HI, haloshift=0, unscale=gv%H_to_MKS)

1641 call hchksum(ent_s, "before triDiagTS ent_s ", g%HI, haloshift=0, unscale=gv%H_to_MKS)

1642 endif

1643

1644 call cpu_clock_begin(id_clock_tridiag)

1645 ! Keep salinity from falling below a small but positive threshold.

1646 ! This constraint is needed for SIS1 ice model, which can extract

1647 ! more salt than is present in the ocean. SIS2 does not suffer

1648 ! from this limitation, in which case we can let salinity=0 and still

1649 ! have salt conserved with SIS2 ice. So for SIS2, we can run with

1650 ! BOUND_SALINITY=False in MOM.F90.

1651 if (associated(tv%S) .and. associated(tv%salt_deficit)) &

1652 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

1653

1654 if (cs%diabatic_diff_tendency_diag) then

1655 do k=1,nz ; do j=js,je ; do i=is,ie

1656 temp_diag(i,j,k) = tv%T(i,j,k)

1657 saln_diag(i,j,k) = tv%S(i,j,k)

1658 enddo ; enddo ; enddo

1659 endif

1660

1661 ! Find the vertical distances across layers, which may have been modified by the net surface flux

1662 call thickness_to_dz(h, tv, dz, g, gv, us)

1663

1664 ! set ent_t=dt*Kd_heat/h_int and est_s=dt*Kd_salt/h_int on interfaces for use in the tridiagonal solver.

1665 do j=js,je ; do i=is,ie

1666 ent_t(i,j,1) = 0. ; ent_t(i,j,nz+1) = 0.

1667 ent_s(i,j,1) = 0. ; ent_s(i,j,nz+1) = 0.

1668 enddo ; enddo

1669

1670 !$OMP parallel do default(shared) private(I_dzval)

1671 do k=2,nz ; do j=js,je ; do i=is,ie

1672 i_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))

1673 ent_t(i,j,k) = dt * i_dzval * kd_heat(i,j,k)

1674 ent_s(i,j,k) = dt * i_dzval * kd_salt(i,j,k)

1675 enddo ; enddo ; enddo

1676 if (showcalltree) call calltree_waypoint("done setting ent_t and ent_t from Kd_heat and " //&

1677 "Kd_salt (diabatic_ALE)")

1678

1679 ! Changes T and S via the tridiagonal solver; no change to h

1680 call tracer_vertdiff_eulerian(h, ent_t, dt, tv%T, g, gv)

1681 call tracer_vertdiff_eulerian(h, ent_s, dt, tv%S, g, gv)

1682

1683 ! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below

1684 if (cs%diabatic_diff_tendency_diag) then

1685 call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, g, gv, us, cs)

1686 endif

1687 call cpu_clock_end(id_clock_tridiag)

1688

1689 if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")

1690

1691 endif ! endif corresponding to if (associated(tv%T))

1692

1693 if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)

1694

1695 if (cs%debug) then

1696 call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)

1697 call mom_thermovar_chksum("after mixed layer ", tv, g, us)

1698 endif

1699

1700 ! Whenever thickness changes let the diag manager know, as the

1701 ! target grids for vertical remapping may need to be regenerated.

1702 call diag_update_remap_grids(cs%diag)

1703

1704 ! Set diffusivities for VBF diagnostics if enabled

1705 if (cs%use_energetic_PBL .and. associated(cs%VBF%Kd_ePBL)) cs%VBF%Kd_ePBL(:,:,:) = kd_epbl(:,:,:)

1706 if (associated(cs%VBF%Kd_salt)) cs%VBF%Kd_temp(:,:,:) = kd_heat(:,:,:)

1707 if (associated(cs%VBF%Kd_temp)) cs%VBF%Kd_salt(:,:,:) = kd_salt(:,:,:)

1708

1709 ! Diagnose the diapycnal diffusivities and other related quantities.

1710 if (cs%id_Kd_heat > 0) call post_data(cs%id_Kd_heat, kd_heat, cs%diag)

1711 if (cs%id_Kd_salt > 0) call post_data(cs%id_Kd_salt, kd_salt, cs%diag)

1712 if (cs%id_Kd_ePBL > 0) call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)

1713 if (cs%id_Kd_int > 0) then

1714 if (cs%double_diffuse .or. cs%useKPP) then

1715 ! Using this as a work array might cause confusion.

1716 do k=1,nz ; do j=js,je ; do i=is,ie

1717 kd_heat(i,j,k) = min(kd_heat(i,j,k), kd_salt(i,j,k))

1718 enddo ; enddo ; enddo

1719 endif

1720 call post_data(cs%id_Kd_int, kd_heat, cs%diag)

1721 endif

1722

1723 if (cs%id_ea_t > 0) call post_data(cs%id_ea_t, ent_t(:,:,1:nz), cs%diag)

1724 if (cs%id_eb_t > 0) call post_data(cs%id_eb_t, ent_t(:,:,2:nz+1), cs%diag)

1725 if (cs%id_ea_s > 0) call post_data(cs%id_ea_s, ent_s(:,:,1:nz), cs%diag)

1726 if (cs%id_eb_s > 0) call post_data(cs%id_eb_s, ent_s(:,:,2:nz+1), cs%diag)

1727

1728 idt = 1.0 / dt

1729 if (cs%id_Tdif > 0) then

1730 do j=js,je ; do i=is,ie

1731 tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0

1732 enddo ; enddo

1733 !$OMP parallel do default(shared)

1734 do k=2,nz ; do j=js,je ; do i=is,ie

1735 tdif_flx(i,j,k) = (idt * ent_t(i,j,k)) * (tv%T(i,j,k-1) - tv%T(i,j,k))

1736 enddo ; enddo ; enddo

1737 if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)

1738 endif

1739 if (cs%id_Sdif > 0) then

1740 do j=js,je ; do i=is,ie

1741 sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0

1742 enddo ; enddo

1743 !$OMP parallel do default(shared)

1744 do k=2,nz ; do j=js,je ; do i=is,ie

1745 sdif_flx(i,j,k) = (idt * ent_s(i,j,k)) * (tv%S(i,j,k-1) - tv%S(i,j,k))

1746 enddo ; enddo ; enddo

1747 if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)

1748 endif

1749

1750 if (cs%Use_KdWork_diag .or. cs%Use_N2_diag) then

1751 n2_salt(:,:,:) = 0.0

1752 n2_temp(:,:,:) = 0.0

1753 !Compute N2 and don't mask negatives here

1754 eosdom(:) = eos_domain(g%HI)

1755 if (nonbous) then

1756 !$OMP parallel do default(shared)

1757 do j=js,je

1758 if (associated(tv%p_surf)) then

1759 do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo

1760 else

1761 do i=is,ie ; p_i(i) = 0.0 ; enddo

1762 endif

1763 do k=2,nz

1764 do i=is,ie

1765 p_i(i) = p_i(i) + h_to_pres * h(i,j,k-1)

1766 enddo

1767 t_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))

1768 s_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))

1769 call calculate_specific_vol_derivs(t_i, s_i, p_i, dspv_dt, dspv_ds, tv%eqn_of_state, eosdom)

1770 do i=is,ie

1771 i_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))

1772 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))

1773 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))

1774 enddo

1775 enddo

1776 enddo

1777 else

1778 !$OMP parallel do default(shared)

1779 do j=js,je

1780 if (associated(tv%p_surf)) then

1781 do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo

1782 else

1783 do i=is,ie ; p_i(i) = 0.0 ; enddo

1784 endif

1785 do k=2,nz

1786 do i=is,ie

1787 p_i(i) = p_i(i) + h_to_pres* h(i,j,k-1)

1788 enddo

1789 t_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))

1790 s_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))

1791 call calculate_density_derivs(t_i, s_i, p_i, drhodt, drhods, tv%eqn_of_state, eosdom)

1792 do i=is,ie

1793 i_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))

1794 n2_salt(i,j,k) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (drhods(i) * (g_rho0 * i_h))

1795 n2_temp(i,j,k) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (drhodt(i) * (g_rho0 * i_h))

1796 enddo

1797 enddo

1798 enddo

1799 endif

1800 if (cs%id_N2_dd>0) call post_data(cs%id_N2_dd, n2_salt(:,:,:)+n2_temp(:,:,:), cs%diag)

1801 if (cs%id_N2_salt_dd>0) call post_data(cs%id_N2_salt_dd, n2_salt, cs%diag)

1802 if (cs%id_N2_temp_dd>0) call post_data(cs%id_N2_temp_dd, n2_temp, cs%diag)

1803

1804 if (cs%Use_KdWork_diag) then

1805 call kdwork_diagnostics(g,gv,us,cs%diag,cs%VBF,n2_salt,n2_temp,dz)

1806 endif

1807

1808 call deallocate_vbf_cs(cs%VBF)

1809

1810 endif

1811

1812 ! mixing of passive tracers from massless boundary layers to interior

1813 call cpu_clock_begin(id_clock_tracers)

1814

1815 if (cs%mix_boundary_tracer_ALE) then

1816 tr_ea_bbl = sqrt(dt * cs%Kd_BBL_tr)

1817 !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)

1818 do j=js,je

1819 do i=is,ie

1820 htot(i) = 0.0

1821 in_boundary(i) = (g%mask2dT(i,j) > 0.0)

1822 enddo

1823 do k=nz,2,-1 ; do i=is,ie

1824 if (in_boundary(i)) then

1825 htot(i) = htot(i) + h(i,j,k)

1826 ! If diapycnal mixing has been suppressed because this is a massless layer near the

1827 ! bottom, add some mixing of tracers between these layers. This flux is based on the

1828 ! harmonic mean of the two thicknesses, following what is done in layered mode. Kd_min_tr

1829 ! should be much less than the values in Kd_salt, perhaps a molecular diffusivity.

1830 add_ent = (dt * cs%Kd_min_tr) * &

1831 ((dz(i,j,k-1)+dz(i,j,k) + dz_neglect) / (dz(i,j,k-1)*dz(i,j,k) + dz_neglect2)) - &

1832 ent_s(i,j,k)

1833 if (htot(i) < tr_ea_bbl) then

1834 add_ent = max(0.0, add_ent, (tr_ea_bbl - htot(i)) - ent_s(i,j,k))

1835 elseif (add_ent < 0.0) then

1836 add_ent = 0.0 ; in_boundary(i) = .false.

1837 endif

1838

1839 ent_s(i,j,k) = ent_s(i,j,k) + add_ent

1840 endif

1841 enddo ; enddo

1842 enddo

1843 endif ! (CS%mix_boundary_tracer_ALE)

1844

1845 ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied

1846 call call_tracer_column_fns(h_orig, h, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), fluxes, bld, dt, &

1847 g, gv, us, tv, cs%optics, cs%tracer_flow_CSp, cs%debug, &

1848 kpp_csp=cs%KPP_CSp, &

1849 nonlocaltrans=kpp_nltscalar, &

1850 evap_cfl_limit=cs%evap_CFL_limit, &

1851 minimum_forcing_depth=cs%minimum_forcing_depth, h_bl=visc%h_ML)

1852

1853 call cpu_clock_end(id_clock_tracers)

1854

1855 ! Apply ALE sponge

1856 if (cs%use_sponge .and. associated(cs%ALE_sponge_CSp)) then

1857 call cpu_clock_begin(id_clock_sponge)

1858 call apply_ale_sponge(h, tv, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)

1859 call cpu_clock_end(id_clock_sponge)

1860 if (cs%debug) then

1861 call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)

1862 call mom_thermovar_chksum("apply_sponge ", tv, g, us)

1863 endif

1864 endif ! CS%use_sponge

1865

1866 ! Apply data assimilation incremental update -oda_incupd-

1867 if (cs%use_oda_incupd .and. associated(cs%oda_incupd_CSp)) then

1868 call mom_mesg("Starting ODA_INCUPD ", 5)

1869 call cpu_clock_begin(id_clock_oda_incupd)

1870 call apply_oda_incupd(h, tv, u, v, dt, g, gv, us, cs%oda_incupd_CSp)

1871 call cpu_clock_end(id_clock_oda_incupd)

1872 if (cs%debug) then

1873 call mom_state_chksum("apply_oda_incupd ", u, v, h, g, gv, us, haloshift=0)

1874 call mom_thermovar_chksum("apply_oda_incupd ", tv, g, us)

1875 endif

1876 endif ! CS%use_oda_incupd

1877

1878

1879 call cpu_clock_begin(id_clock_pass)

1880 ! visc%Kv_slow is not in the group pass because it has larger vertical extent.

1881 if (associated(visc%Kv_slow)) &

1882 call pass_var(visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)

1883 call cpu_clock_end(id_clock_pass)

1884

1885 call disable_averaging(cs%diag)

1886

1887 if (showcalltree) call calltree_leave("diabatic_ALE()")

1888

1889end subroutine diabatic_ale

1890

1891!> Imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers

1892!! using the original MOM6 algorithms.

1893subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &

1894 G, GV, US, CS, Waves)

1895 type(ocean_grid_type), intent(inout) :: G !< ocean grid structure

1896 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

1897 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

1898 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]

1899 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]

1900 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2]

1901 type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields

1902 !! unused have NULL ptrs

1903 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Active mixed layer depth [Z ~> m]

1904 type(forcing), intent(inout) :: fluxes !< points to forcing fields

1905 !! unused fields have NULL ptrs

1906 type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities,

1907 !! BBL properties and related fields

1908 type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum

1909 !! equations, to enable the later derived

1910 !! diagnostics, like energy budgets

1911 type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations

1912 real, intent(in) :: dt !< time increment [T ~> s]

1913 type(time_type), intent(in) :: Time_end !< Time at the end of the interval

1914 type(diabatic_cs), pointer :: CS !< module control structure

1915 type(wave_parameters_cs), pointer :: Waves !< Surface gravity waves

1916

1917 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &

1918 ea, & ! amount of fluid entrained from the layer above within

1919 ! one time step [H ~> m or kg m-2]

1920 eb, & ! amount of fluid entrained from the layer below within

1921 ! one time step [H ~> m or kg m-2]

1922 kd_lay, & ! diapycnal diffusivity of layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1923 h_orig, & ! initial layer thicknesses [H ~> m or kg m-2]

1924 dz, & ! The vertical distance between interfaces around a layer [Z ~> m]

1925 hold, & ! layer thickness before diapycnal entrainment, and later the initial

1926 ! layer thicknesses (if a mixed layer is used) [H ~> m or kg m-2]

1927 dz_old, & ! The initial vertical distance between interfaces around a layer

1928 ! or the distance before entrainment [Z ~> m]

1929 u_h, & ! Zonal velocities at thickness points after entrainment [L T-1 ~> m s-1]

1930 v_h, & ! Meridional velocities at thickness points after entrainment [L T-1 ~> m s-1]

1931 temp_diag, & ! Diagnostic array of previous temperatures [C ~> degC]

1932 saln_diag ! Diagnostic array of previous salinity [S ~> ppt]

1933 real, dimension(SZI_(G),SZJ_(G)) :: &

1934 h_MLD, & ! Active mixed layer thickness [H ~> m or kg m-2].

1935 U_star, & ! The friction velocity [Z T-1 ~> m s-1].

1936 KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

1937 KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

1938 Rcv_ml ! Coordinate density of mixed layer [R ~> kg m-3], used for applying sponges

1939

1940 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &

1941 ! These are targets so that the space can be shared with eaml & ebml.

1942 eatr, & ! The equivalent of ea for tracers, which differs from ea in that it tends to

1943 ! homogenize tracers in massless layers near the boundaries [H ~> m or kg m-2]

1944 ebtr ! The equivalent of eb for tracers, which differs from eb in that it tends to

1945 ! homogenize tracers in massless layers near the boundaries [H ~> m or kg m-2]

1946

1947 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &

1948 Kd_int, & ! diapycnal diffusivity of interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1949 Kd_heat, & ! diapycnal diffusivity of heat [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1950 Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1951 Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to

1952 ! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1953 kd_extra_s , & ! The extra diffusivity of salinity due to double diffusion relative to

1954 ! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1955 kpp_nltheat, & ! KPP non-local transport for heat [nondim]

1956 kpp_nltscalar, & ! KPP non-local transport for scalars [nondim]

1957 kpp_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]

1958 tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

1959 tadv_flx, & ! advective diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

1960 sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

1961 sadv_flx ! advective diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

1962

1963 ! The following 3 variables are only used with a bulk mixed layer.

1964 real, pointer, dimension(:,:,:) :: &

1965 eaml, & ! The equivalent of ea due to mixed layer processes [H ~> m or kg m-2].

1966 ebml ! The equivalent of eb due to mixed layer processes [H ~> m or kg m-2].

1967 ! eaml and ebml are pointers to eatr and ebtr so as to reuse the memory as

1968 ! the arrays are not needed at the same time.

1969

1970 integer :: kb(SZI_(G),SZJ_(G)) ! index of the lightest layer denser

1971 ! than the buffer layer [nondim]

1972

1973 real :: p_ref_cv(SZI_(G)) ! Reference pressure for the potential density that defines the

1974 ! coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].

1975

1976 logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below,

1977 ! where massive is defined as sufficiently thick that

1978 ! the no-flux boundary conditions have not restricted

1979 ! the entrainment - usually sqrt(Kd*dt).

1980

1981 real :: h_neglect ! A thickness that is so small it is usually lost

1982 ! in roundoff and can be neglected [H ~> m or kg m-2]

1983 real :: dz_neglect ! A vertical distance that is so small it is usually lost

1984 ! in roundoff and can be neglected [Z ~> m]

1985 real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2]

1986 real :: net_ent ! The net of ea-eb at an interface [H ~> m or kg m-2]

1987 real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2]

1988 real :: eaval ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]

1989 real :: hval ! hval is 2*h at velocity grid points [H ~> m or kg m-2]

1990 real :: h_tr ! h_tr is h at tracer points with a tiny thickness

1991 ! added to ensure positive definiteness [H ~> m or kg m-2]

1992 real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is

1993 ! coupled to the bottom within a timestep [H ~> m or kg m-2]

1994

1995 real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].

1996 real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]

1997 real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2]

1998 real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim]

1999 real :: c1(SZIB_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim]

2000

2001 real :: dt_mix ! The amount of time over which to apply mixing [T ~> s]

2002 real :: Idt ! The inverse time step [T-1 ~> s-1]

2003

2004 integer :: dir_flag ! An integer encoding the directions in which to do halo updates.

2005 logical :: showCallTree ! If true, show the call tree

2006 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state

2007 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, halo

2008

2009 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

2010 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

2011 nkmb = gv%nk_rho_varies

2012 h_neglect = gv%H_subroundoff

2013 dz_neglect = gv%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect

2014 kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0

2015

2016 showcalltree = calltree_showquery()

2017 if (showcalltree) call calltree_enter("layered_diabatic(), MOM_diabatic_driver.F90")

2018

2019 ! set equivalence between the same bits of memory for these arrays

2020 eaml => eatr ; ebml => ebtr

2021

2022 ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep

2023 call enable_averages(dt, time_end, cs%diag)

2024

2025 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then

2026 halo = cs%halo_TS_diff

2027 !$OMP parallel do default(shared)

2028 do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo

2029 h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0

2030 enddo ; enddo ; enddo

2031 endif

2032

2033 if (cs%use_geothermal) then

2034 call cpu_clock_begin(id_clock_geothermal)

2035 call geothermal_entraining(h, tv, dt, eaml, ebml, g, gv, us, cs%geothermal, halo=cs%halo_TS_diff)

2036 call cpu_clock_end(id_clock_geothermal)

2037 if (showcalltree) call calltree_waypoint("geothermal (diabatic)")

2038 if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)

2039 endif

2040

2041 ! Whenever thickness changes let the diag manager know, target grids

2042 ! for vertical remapping may need to be regenerated.

2043 call diag_update_remap_grids(cs%diag)

2044

2045 ! Set_pen_shortwave estimates the optical properties of the water column.

2046 ! It will need to be modified later to include information about the

2047 ! biological properties and layer thicknesses.

2048 if (associated(cs%optics)) &

2049 call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity, cs%tracer_flow_CSp)

2050

2051 if (cs%use_bulkmixedlayer) then

2052 if (cs%debug) call mom_forcing_chksum("Before mixedlayer", fluxes, g, us, haloshift=0)

2053

2054 if (cs%ML_mix_first > 0.0) then

2055! This subroutine

2056! (1) Cools the mixed layer.

2057! (2) Performs convective adjustment by mixed layer entrainment.

2058! (3) Heats the mixed layer and causes it to detrain to

2059! Monin-Obukhov depth or minimum mixed layer depth.

2060! (4) Uses any remaining TKE to drive mixed layer entrainment.

2061! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate)

2062 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)

2063

2064 call cpu_clock_begin(id_clock_mixedlayer)

2065 if (cs%ML_mix_first < 1.0) then

2066 ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???)

2067 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*cs%ML_mix_first, &

2068 eaml, ebml, g, gv, us, cs%bulkmixedlayer, cs%optics, &

2069 bld, h_mld, cs%aggregate_FW_forcing, dt, last_call=.false.)

2070 else

2071 ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???)

2072 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, &

2073 g, gv, us, cs%bulkmixedlayer, cs%optics, &

2074 bld, h_mld, cs%aggregate_FW_forcing, dt, last_call=.true.)

2075 if (associated(visc%h_ML)) visc%h_ML(:,:) = h_mld(:,:)

2076 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

2077 endif

2078

2079 ! Keep salinity from falling below a small but positive threshold.

2080 ! This constraint is needed for SIS1 ice model, which can extract

2081 ! more salt than is present in the ocean. SIS2 does not suffer

2082 ! from this limitation, in which case we can let salinity=0 and still

2083 ! have salt conserved with SIS2 ice. So for SIS2, we can run with

2084 ! BOUND_SALINITY=False in MOM.F90.

2085 if (associated(tv%S) .and. associated(tv%salt_deficit)) &

2086 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

2087 call cpu_clock_end(id_clock_mixedlayer)

2088 if (cs%debug) then

2089 call mom_state_chksum("After mixedlayer ", u, v, h, g, gv, us, haloshift=0)

2090 call mom_forcing_chksum("After mixedlayer", fluxes, g, us, haloshift=0)

2091 endif

2092 if (showcalltree) call calltree_waypoint("done with 1st bulkmixedlayer (diabatic)")

2093 if (cs%debugConservation) call mom_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, g, gv, us)

2094 endif

2095 endif

2096

2097 if (cs%debug) &

2098 call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)

2099 if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then

2100 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then

2101 call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, us, eaml, ebml)

2102 if (cs%debug) then

2103 call hchksum(eaml, "after find_uv_at_h eaml", g%HI, unscale=gv%H_to_MKS)

2104 call hchksum(ebml, "after find_uv_at_h ebml", g%HI, unscale=gv%H_to_MKS)

2105 endif

2106 else

2107 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)

2108 endif

2109 if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")

2110 endif

2111

2112 call cpu_clock_begin(id_clock_set_diffusivity)

2113 ! Sets: Kd_lay, Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb

2114 ! Also changes: visc%Kd_shear and visc%Kv_shear

2115 if ((cs%halo_TS_diff > 0) .and. (cs%ML_mix_first > 0.0)) then

2116 if (associated(tv%T)) call pass_var(tv%T, g%Domain, halo=cs%halo_TS_diff, complete=.false.)

2117 if (associated(tv%S)) call pass_var(tv%S, g%Domain, halo=cs%halo_TS_diff, complete=.false.)

2118 call pass_var(h, g%domain, halo=cs%halo_TS_diff, complete=.true.)

2119 endif

2120

2121 ! Update derived thermodynamic quantities.

2122 if ((cs%ML_mix_first > 0.0) .and. allocated(tv%SpV_avg)) then

2123 call calc_derived_thermo(tv, h, g, gv, us, halo=cs%halo_TS_diff)

2124 endif

2125

2126 if (cs%debug) &

2127 call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)

2128 if (cs%double_diffuse) then

2129 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, kd_int, g, gv, us, &

2130 cs%set_diff_CSp, cs%VBF, kd_lay=kd_lay, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)

2131 else

2132 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, kd_int, g, gv, us, &

2133 cs%set_diff_CSp, cs%VBF, kd_lay=kd_lay)

2134 endif

2135 call cpu_clock_end(id_clock_set_diffusivity)

2136 if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")

2137

2138 if (cs%debug) then

2139 call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)

2140 call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)

2141 call mom_thermovar_chksum("after set_diffusivity ", tv, g, us)

2142 call hchksum(kd_lay, "after set_diffusivity Kd_lay", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

2143 call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

2144 endif

2145

2146

2147 if (cs%useKPP) then

2148 call cpu_clock_begin(id_clock_kpp)

2149

2150 ! NOTE: The following do not require initialization, but their checksums do

2151 ! require initialization, and past versions were initialized to zero.

2152 kpp_nltheat(:,:,:) = 0.

2153 kpp_nltscalar(:,:,:) = 0.

2154 kpp_buoy_flux(:,:,:) = 0.

2155 kpp_temp_flux(:,:) = 0.

2156 kpp_salt_flux(:,:) = 0.

2157

2158 ! KPP needs the surface buoyancy flux but does not update state variables.

2159 ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.

2160 ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux

2161 ! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second)

2162 ! unlike other instances where the fluxes are integrated in time over a time-step.

2163 call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &

2164 kpp_buoy_flux, kpp_temp_flux, kpp_salt_flux)

2165 ! The KPP scheme calculates boundary layer diffusivities and non-local transport.

2166

2167 ! Set diffusivities for heat and salt separately

2168

2169 if (cs%double_diffuse) then

2170 ! Add contribution from double diffusion

2171 !$OMP parallel do default(shared)

2172 do k=1,nz+1 ; do j=js,je ; do i=is,ie

2173 kd_salt(i,j,k) = kd_int(i,j,k) + kd_extra_s(i,j,k)

2174 kd_heat(i,j,k) = kd_int(i,j,k) + kd_extra_t(i,j,k)

2175 enddo ; enddo ; enddo

2176 else

2177 !$OMP parallel do default(shared)

2178 do k=1,nz+1 ; do j=js,je ; do i=is,ie

2179 kd_salt(i,j,k) = kd_int(i,j,k)

2180 kd_heat(i,j,k) = kd_int(i,j,k)

2181 enddo ; enddo ; enddo

2182 endif

2183

2184 ! Determine the friction velocity, perhaps using the evovling surface density.

2185 call find_ustar(fluxes, tv, u_star, g, gv, us)

2186

2187 if ( associated(fluxes%lamult) ) then

2188 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &

2189 u_star, kpp_buoy_flux, waves=waves, lamult=fluxes%lamult)

2190

2191 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, tv, u_star, kpp_buoy_flux, kd_heat, &

2192 kd_salt, visc%Kv_shear, kpp_nltheat, kpp_nltscalar, waves=waves, lamult=fluxes%lamult)

2193 else

2194 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &

2195 u_star, kpp_buoy_flux, waves=waves)

2196

2197 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, tv, u_star, kpp_buoy_flux, kd_heat, &

2198 kd_salt, visc%Kv_shear, kpp_nltheat, kpp_nltscalar, waves=waves)

2199 endif

2200

2201 call kpp_get_bld(cs%KPP_CSp, bld(:,:), g, us)

2202 ! If visc%MLD or visc%h_ML exist, copy KPP's BLD into them with appropriate conversions.

2203 if (associated(visc%h_ML)) call convert_mld_to_ml_thickness(bld, h, visc%h_ML, tv, g, gv)

2204 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

2205 if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = kpp_buoy_flux(:,:,1)

2206

2207 if (.not. cs%KPPisPassive) then

2208 !$OMP parallel do default(shared)

2209 do k=1,nz+1 ; do j=js,je ; do i=is,ie

2210 kd_int(i,j,k) = min( kd_salt(i,j,k), kd_heat(i,j,k) )

2211 enddo ; enddo ; enddo

2212 if (cs%double_diffuse) then

2213 !$OMP parallel do default(shared)

2214 do k=1,nz+1 ; do j=js,je ; do i=is,ie

2215 kd_extra_s(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))

2216 kd_extra_t(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))

2217 enddo ; enddo ; enddo

2218 endif

2219 endif ! not passive

2220

2221 call cpu_clock_end(id_clock_kpp)

2222 if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")

2223 if (cs%debug) then

2224 call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)

2225 call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)

2226 call mom_thermovar_chksum("after KPP", tv, g, us)

2227 call hchksum(kd_lay, "after KPP Kd_lay", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

2228 call hchksum(kd_int, "after KPP Kd_Int", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

2229 endif

2230 endif ! endif for KPP

2231

2232 ! Add vertical diff./visc. due to convection (computed via CVMix)

2233 if (cs%use_CVMix_conv) then

2234 call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv, bld, kd_int, visc%Kv_shear)

2235 endif

2236

2237 if (cs%useKPP) then

2238 call cpu_clock_begin(id_clock_kpp)

2239 if (cs%debug) then

2240 call hchksum(kpp_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, &

2241 unscale=us%C_to_degC*gv%H_to_m*us%s_to_T)

2242 call hchksum(kpp_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, &

2243 unscale=us%S_to_ppt*gv%H_to_m*us%s_to_T)

2244 call hchksum(kpp_nltheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)

2245 call hchksum(kpp_nltscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)

2246 endif

2247 ! Apply non-local transport of heat and salt

2248 ! Changes: tv%T, tv%S

2249 call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, kpp_nltheat, kpp_temp_flux, &

2250 dt, tv%tr_T, tv%T, tv%C_p)

2251 call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, kpp_nltscalar, kpp_salt_flux, &

2252 dt, tv%tr_S, tv%S)

2253 call cpu_clock_end(id_clock_kpp)

2254 if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")

2255 if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)

2256

2257 if (cs%debug) then

2258 call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)

2259 call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)

2260 call mom_thermovar_chksum("after KPP_applyNLT ", tv, g, us)

2261 endif

2262 endif ! endif for KPP

2263

2264 ! Differential diffusion done here.

2265 ! Changes: tv%T, tv%S

2266 if (cs%double_diffuse .and. associated(tv%T)) then

2267

2268 call cpu_clock_begin(id_clock_differential_diff)

2269 call differential_diffuse_t_s(h, tv%T, tv%S, kd_extra_t, kd_extra_s, tv, dt, g, gv)

2270 call cpu_clock_end(id_clock_differential_diff)

2271 if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")

2272 if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g, gv, us)

2273

2274 ! increment heat and salt diffusivity.

2275 ! CS%useKPP==.true. already has extra_T and extra_S included

2276 if (.not. cs%useKPP) then

2277 !$OMP parallel do default(shared)

2278 do k=2,nz ; do j=js,je ; do i=is,ie

2279 kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)

2280 kd_salt(i,j,k) = kd_salt(i,j,k) + kd_extra_s(i,j,k)

2281 enddo ; enddo ; enddo

2282 endif

2283

2284 endif

2285

2286 ! Calculate layer entrainments and detrainments from diffusivities and differences between

2287 ! layer and target densities (i.e. do remapping as well as diffusion).

2288 call cpu_clock_begin(id_clock_entrain)

2289 ! Calculate appropriately limited diapycnal mass fluxes to account

2290 ! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb

2291 call entrainment_diffusive(h, tv, fluxes, dt, g, gv, us, cs%entrain_diffusive, &

2292 ea, eb, kb, kd_lay=kd_lay, kd_int=kd_int)

2293 call cpu_clock_end(id_clock_entrain)

2294 if (showcalltree) call calltree_waypoint("done with Entrainment_diffusive (diabatic)")

2295

2296 if (cs%debug) then

2297 call mom_forcing_chksum("after calc_entrain ", fluxes, g, us, haloshift=0)

2298 call mom_thermovar_chksum("after calc_entrain ", tv, g, us)

2299 call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, us, haloshift=0)

2300 call hchksum(ea, "after calc_entrain ea", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2301 call hchksum(eb, "after calc_entrain eb", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2302 endif

2303

2304 ! Save fields before boundary forcing is applied for tendency diagnostics

2305 if (cs%boundary_forcing_tendency_diag) then

2306 do k=1,nz ; do j=js,je ; do i=is,ie

2307 temp_diag(i,j,k) = tv%T(i,j,k)

2308 saln_diag(i,j,k) = tv%S(i,j,k)

2309 enddo ; enddo ; enddo

2310 endif

2311

2312 ! Update h according to divergence of the difference between

2313 ! ea and eb. We keep a record of the original h in hold.

2314 ! In the following, the checks for negative values are to guard

2315 ! against instances where entrainment drives a layer to

2316 ! negative thickness. This situation will never happen if

2317 ! enough iterations are permitted in Calculate_Entrainment.

2318 ! Even if too few iterations are allowed, it is still guarded

2319 ! against. In other words the checks are probably unnecessary.

2320 !$OMP parallel do default(shared)

2321 do j=js,je

2322 do i=is,ie

2323 hold(i,j,1) = h(i,j,1)

2324 h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2))

2325 hold(i,j,nz) = h(i,j,nz)

2326 h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1))

2327 if (h(i,j,1) <= 0.0) then

2328 h(i,j,1) = gv%Angstrom_H

2329 endif

2330 if (h(i,j,nz) <= 0.0) then

2331 h(i,j,nz) = gv%Angstrom_H

2332 endif

2333 enddo

2334 do k=2,nz-1 ; do i=is,ie

2335 hold(i,j,k) = h(i,j,k)

2336 h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &

2337 (eb(i,j,k) - ea(i,j,k+1)))

2338 if (h(i,j,k) <= 0.0) then

2339 h(i,j,k) = gv%Angstrom_H

2340 endif

2341 enddo ; enddo

2342 enddo

2343 ! Checks for negative thickness may have changed layer thicknesses

2344 call diag_update_remap_grids(cs%diag)

2345

2346 if (cs%debug) then

2347 call mom_state_chksum("after negative check ", u, v, h, g, gv, us, haloshift=0)

2348 call mom_forcing_chksum("after negative check ", fluxes, g, us, haloshift=0)

2349 call mom_thermovar_chksum("after negative check ", tv, g, us)

2350 endif

2351 if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")

2352 if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)

2353

2354 ! Here, T and S are updated according to ea and eb.

2355 ! If using the bulk mixed layer, T and S are also updated

2356 ! by surface fluxes (in fluxes%*).

2357 ! This is a very long block.

2358 if (cs%use_bulkmixedlayer) then

2359

2360 if (associated(tv%T)) then

2361 call cpu_clock_begin(id_clock_tridiag)

2362 ! Temperature and salinity (as state variables) are treated

2363 ! differently from other tracers to insure massless layers that

2364 ! are lighter than the mixed layer have temperatures and salinities

2365 ! that correspond to their prescribed densities.

2366 if (cs%massless_match_targets) then

2367 !$OMP parallel do default (shared) private(h_tr,b1,d1,c1,b_denom_1)

2368 do j=js,je

2369 do i=is,ie

2370 h_tr = hold(i,j,1) + h_neglect

2371 b1(i) = 1.0 / (h_tr + eb(i,j,1))

2372 d1(i) = h_tr * b1(i)

2373 tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1))

2374 tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1))

2375 enddo

2376 do k=2,nkmb ; do i=is,ie

2377 c1(i,k) = eb(i,j,k-1) * b1(i)

2378 h_tr = hold(i,j,k) + h_neglect

2379 b_denom_1 = h_tr + d1(i)*ea(i,j,k)

2380 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))

2381 if (k<nkmb) d1(i) = b_denom_1 * b1(i)

2382 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))

2383 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))

2384 enddo ; enddo

2385

2386 do k=nkmb+1,nz ; do i=is,ie

2387 if (k == kb(i,j)) then

2388 c1(i,k) = eb(i,j,k-1) * b1(i)

2389 d1(i) = (((eb(i,j,nkmb)-eb(i,j,k-1)) + hold(i,j,nkmb) + h_neglect) + &

2390 d1(i)*ea(i,j,nkmb)) * b1(i)

2391 h_tr = hold(i,j,k) + h_neglect

2392 b_denom_1 = h_tr + d1(i)*ea(i,j,k)

2393 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))

2394 d1(i) = b_denom_1 * b1(i)

2395 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,nkmb))

2396 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,nkmb))

2397 elseif (k > kb(i,j)) then

2398 c1(i,k) = eb(i,j,k-1) * b1(i)

2399 h_tr = hold(i,j,k) + h_neglect

2400 b_denom_1 = h_tr + d1(i)*ea(i,j,k)

2401 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))

2402 d1(i) = b_denom_1 * b1(i)

2403 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))

2404 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))

2405 elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j))

2406 ! The bottommost buffer layer might entrain all the mass from some

2407 ! of the interior layers that are thin and lighter in the coordinate

2408 ! density than that buffer layer. The T and S of these newly

2409 ! massless interior layers are unchanged.

2410 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)

2411 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)

2412 endif

2413 enddo ; enddo

2414

2415 do k=nz-1,nkmb,-1 ; do i=is,ie

2416 if (k >= kb(i,j)) then

2417 tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)

2418 tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)

2419 endif

2420 enddo ; enddo

2421 do i=is,ie ; if (kb(i,j) <= nz) then

2422 tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j))

2423 tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j))

2424 endif ; enddo

2425 do k=nkmb-1,1,-1 ; do i=is,ie

2426 tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)

2427 tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)

2428 enddo ; enddo

2429 enddo ! end of j loop

2430 else ! .not. massless_match_targets

2431 ! This simpler form allows T & S to be too dense for the layers

2432 ! between the buffer layers and the interior.

2433 ! Changes: T, S

2434 if (cs%tracer_tridiag) then

2435 call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)

2436 call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)

2437 else

2438 call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)

2439 endif

2440 endif ! massless_match_targets

2441 call cpu_clock_end(id_clock_tridiag)

2442

2443 endif ! endif for associated(T)

2444 if (cs%debugConservation) call mom_state_stats('BML tridiag', u, v, h, tv%T, tv%S, g, gv, us)

2445

2446 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then

2447 ! The mixed layer code has already been called, but there is some needed

2448 ! bookkeeping.

2449 !$OMP parallel do default(shared)

2450 do k=1,nz ; do j=js,je ; do i=is,ie

2451 hold(i,j,k) = h_orig(i,j,k)

2452 ea(i,j,k) = ea(i,j,k) + eaml(i,j,k)

2453 eb(i,j,k) = eb(i,j,k) + ebml(i,j,k)

2454 enddo ; enddo ; enddo

2455 if (cs%debug) then

2456 call hchksum(ea, "after ea = ea + eaml", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2457 call hchksum(eb, "after eb = eb + ebml", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2458 endif

2459 endif

2460

2461 if (cs%ML_mix_first < 1.0) then

2462 ! Call the mixed layer code now, perhaps for a second time.

2463 ! This subroutine (1) Cools the mixed layer.

2464 ! (2) Performs convective adjustment by mixed layer entrainment.

2465 ! (3) Heats the mixed layer and causes it to detrain to

2466 ! Monin-Obukhov depth or minimum mixed layer depth.

2467 ! (4) Uses any remaining TKE to drive mixed layer entrainment.

2468 ! (5) Possibly splits the buffer layer into two isopycnal layers.

2469

2470 call find_uv_at_h(u, v, hold, u_h, v_h, g, gv, us, ea, eb)

2471 if (cs%debug) call mom_state_chksum("find_uv_at_h1 ", u, v, h, g, gv, us, haloshift=0)

2472

2473 dt_mix = min(dt, dt*(1.0 - cs%ML_mix_first))

2474 call cpu_clock_begin(id_clock_mixedlayer)

2475 ! Changes: h, tv%T, tv%S, ea and eb (G is also inout???)

2476 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, &

2477 g, gv, us, cs%bulkmixedlayer, cs%optics, &

2478 bld, h_mld, cs%aggregate_FW_forcing, dt, last_call=.true.)

2479 if (associated(visc%h_ML)) visc%h_ML(:,:) = h_mld(:,:)

2480 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

2481

2482 ! Keep salinity from falling below a small but positive threshold.

2483 ! This constraint is needed for SIS1 ice model, which can extract

2484 ! more salt than is present in the ocean. SIS2 does not suffer

2485 ! from this limitation, in which case we can let salinity=0 and still

2486 ! have salt conserved with SIS2 ice. So for SIS2, we can run with

2487 ! BOUND_SALINITY=False in MOM.F90.

2488 if (associated(tv%S) .and. associated(tv%salt_deficit)) &

2489 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

2490

2491 call cpu_clock_end(id_clock_mixedlayer)

2492 if (showcalltree) call calltree_waypoint("done with 2nd bulkmixedlayer (diabatic)")

2493 if (cs%debugConservation) call mom_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, g, gv, us)

2494 endif

2495

2496 else ! following block for when NOT using BULKMIXEDLAYER

2497

2498 ! calculate change in temperature & salinity due to dia-coordinate surface diffusion

2499 if (associated(tv%T)) then

2500

2501 if (cs%debug) then

2502 call hchksum(ea, "before triDiagTS ea ", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2503 call hchksum(eb, "before triDiagTS eb ", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2504 endif

2505 call cpu_clock_begin(id_clock_tridiag)

2506

2507 ! Keep salinity from falling below a small but positive threshold.

2508 ! This constraint is needed for SIS1 ice model, which can extract

2509 ! more salt than is present in the ocean. SIS2 does not suffer

2510 ! from this limitation, in which case we can let salinity=0 and still

2511 ! have salt conserved with SIS2 ice. So for SIS2, we can run with

2512 ! BOUND_SALINITY=False in MOM.F90.

2513 if (associated(tv%S) .and. associated(tv%salt_deficit)) &

2514 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

2515

2516 if (cs%diabatic_diff_tendency_diag) then

2517 do k=1,nz ; do j=js,je ; do i=is,ie

2518 temp_diag(i,j,k) = tv%T(i,j,k)

2519 saln_diag(i,j,k) = tv%S(i,j,k)

2520 enddo ; enddo ; enddo

2521 endif

2522

2523 ! Changes T and S via the tridiagonal solver; no change to h

2524 if (cs%tracer_tridiag) then

2525 call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)

2526 call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)

2527 else

2528 call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)

2529 endif

2530

2531 ! diagnose temperature, salinity, heat, and salt tendencies

2532 ! Note: hold here refers to the thicknesses from before the dual-entrainment when using

2533 ! the bulk mixed layer scheme, so tendencies should be posted on hold.

2534 if (cs%diabatic_diff_tendency_diag) then

2535 call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, us, cs)

2536 if (cs%id_diabatic_diff_h > 0) call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h=hold)

2537 endif

2538

2539 call cpu_clock_end(id_clock_tridiag)

2540 if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")

2541

2542 endif ! endif corresponding to if (associated(tv%T))

2543 if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)

2544

2545 endif ! endif for the BULKMIXEDLAYER block

2546

2547 if (cs%debug) then

2548 call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)

2549 call mom_thermovar_chksum("after mixed layer ", tv, g, us)

2550 call hchksum(ea, "after mixed layer ea", g%HI, unscale=gv%H_to_MKS)

2551 call hchksum(eb, "after mixed layer eb", g%HI, unscale=gv%H_to_MKS)

2552 endif

2553

2554 call cpu_clock_begin(id_clock_remap)

2555 call regularize_layers(h, tv, dt, ea, eb, g, gv, us, cs%regularize_layers)

2556 call cpu_clock_end(id_clock_remap)

2557 if (showcalltree) call calltree_waypoint("done with regularize_layers (diabatic)")

2558 if (cs%debugConservation) call mom_state_stats('regularize_layers', u, v, h, tv%T, tv%S, g, gv, us)

2559

2560 ! Whenever thickness changes let the diag manager know, as the

2561 ! target grids for vertical remapping may need to be regenerated.

2562 if (associated(adp%du_dt_dia) .or. associated(adp%dv_dt_dia)) &

2563 ! Remapped d[uv]dt_dia require east/north halo updates of h

2564 call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)

2565 call diag_update_remap_grids(cs%diag)

2566

2567 ! diagnostics

2568 idt = 1.0 / dt

2569 if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0)) then

2570 do j=js,je ; do i=is,ie

2571 tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0

2572 tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0

2573 enddo ; enddo

2574 !$OMP parallel do default(shared)

2575 do k=2,nz ; do j=js,je ; do i=is,ie

2576 tdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &

2577 (tv%T(i,j,k-1) - tv%T(i,j,k))

2578 tadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &

2579 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))

2580 enddo ; enddo ; enddo

2581 endif

2582 if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0)) then

2583 do j=js,je ; do i=is,ie

2584 sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0

2585 sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0

2586 enddo ; enddo

2587 !$OMP parallel do default(shared)

2588 do k=2,nz ; do j=js,je ; do i=is,ie

2589 sdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &

2590 (tv%S(i,j,k-1) - tv%S(i,j,k))

2591 sadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &

2592 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))

2593 enddo ; enddo ; enddo

2594 endif

2595

2596 ! mixing of passive tracers from massless boundary layers to interior

2597 call cpu_clock_begin(id_clock_tracers)

2598

2599 ! Find the vertical distances across layers.

2600 if (cs%mix_boundary_tracers .or. cs%double_diffuse) &

2601 call thickness_to_dz(h, tv, dz, g, gv, us)

2602 if (cs%double_diffuse) &

2603 call thickness_to_dz(hold, tv, dz_old, g, gv, us)

2604

2605 if (cs%mix_boundary_tracers) then

2606 tr_ea_bbl = sqrt(dt * cs%Kd_BBL_tr)

2607 !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)

2608 do j=js,je

2609 do i=is,ie

2610 ebtr(i,j,nz) = eb(i,j,nz)

2611 htot(i) = 0.0

2612 in_boundary(i) = (g%mask2dT(i,j) > 0.0)

2613 enddo

2614 do k=nz,2,-1 ; do i=is,ie

2615 if (in_boundary(i)) then

2616 htot(i) = htot(i) + h(i,j,k)

2617 ! If diapycnal mixing has been suppressed because this is a massless

2618 ! layer near the bottom, add some mixing of tracers between these

2619 ! layers. This flux is based on the harmonic mean of the two

2620 ! thicknesses, as this corresponds pretty closely (to within

2621 ! differences in the density jumps between layers) with what is done

2622 ! in the calculation of the fluxes in the first place. Kd_min_tr

2623 ! should be much less than the values that have been set in Kd_lay,

2624 ! perhaps a molecular diffusivity.

2625 add_ent = (dt * cs%Kd_min_tr) * &

2626 ((dz(i,j,k-1) + dz(i,j,k) + dz_neglect) / &

2627 (dz(i,j,k-1)*dz(i,j,k) + dz_neglect2)) - &

2628 0.5*(ea(i,j,k) + eb(i,j,k-1))

2629 if (htot(i) < tr_ea_bbl) then

2630 add_ent = max(0.0, add_ent, &

2631 (tr_ea_bbl - htot(i)) - min(ea(i,j,k), eb(i,j,k-1)))

2632 elseif (add_ent < 0.0) then

2633 add_ent = 0.0 ; in_boundary(i) = .false.

2634 endif

2635

2636 ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent

2637 eatr(i,j,k) = ea(i,j,k) + add_ent

2638 else

2639 ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k)

2640 endif

2641 if (cs%double_diffuse) then ; if (kd_extra_s(i,j,k) > 0.0) then

2642 add_ent = (dt * kd_extra_s(i,j,k)) / &

2643 (0.25 * ((dz(i,j,k-1) + dz(i,j,k)) + (dz_old(i,j,k-1) + dz_old(i,j,k))) + dz_neglect)

2644 ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent

2645 eatr(i,j,k) = eatr(i,j,k) + add_ent

2646 endif ; endif

2647 enddo ; enddo

2648 do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo

2649

2650 enddo

2651

2652 call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, bld, dt, g, gv, us, tv, &

2653 cs%optics, cs%tracer_flow_CSp, cs%debug, &

2654 kpp_csp=cs%KPP_CSp, nonlocaltrans=kpp_nltscalar, h_bl=visc%h_ML)

2655

2656 elseif (cs%double_diffuse) then ! extra diffusivity for passive tracers

2657

2658 do j=js,je ; do i=is,ie

2659 ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1)

2660 enddo ; enddo

2661 !$OMP parallel do default(shared) private(add_ent)

2662 do k=nz,2,-1 ; do j=js,je ; do i=is,ie

2663 if (kd_extra_s(i,j,k) > 0.0) then

2664 add_ent = (dt * kd_extra_s(i,j,k)) / &

2665 (0.25 * ((dz(i,j,k-1) + dz(i,j,k)) + (dz_old(i,j,k-1) + dz_old(i,j,k))) + dz_neglect)

2666 else

2667 add_ent = 0.0

2668 endif

2669 ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent

2670 eatr(i,j,k) = ea(i,j,k) + add_ent

2671 enddo ; enddo ; enddo

2672

2673 call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, bld, dt, g, gv, us, tv, &

2674 cs%optics, cs%tracer_flow_CSp, cs%debug, &

2675 kpp_csp=cs%KPP_CSp, nonlocaltrans=kpp_nltscalar, h_bl=visc%h_ML)

2676

2677 else

2678 call call_tracer_column_fns(hold, h, ea, eb, fluxes, bld, dt, g, gv, us, tv, &

2679 cs%optics, cs%tracer_flow_CSp, cs%debug, &

2680 kpp_csp=cs%KPP_CSp, nonlocaltrans=kpp_nltscalar, h_bl=visc%h_ML)

2681

2682 endif ! (CS%mix_boundary_tracers)

2683

2684 call cpu_clock_end(id_clock_tracers)

2685

2686 ! sponges

2687 if (cs%use_sponge) then

2688 call cpu_clock_begin(id_clock_sponge)

2689 ! Layer mode sponge

2690 if (cs%use_bulkmixedlayer .and. associated(tv%eqn_of_state)) then

2691 do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo

2692 eosdom(:) = eos_domain(g%HI)

2693 !$OMP parallel do default(shared)

2694 do j=js,je

2695 call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, rcv_ml(:,j), &

2696 tv%eqn_of_state, eosdom)

2697 enddo

2698 call apply_sponge(h, tv, dt, g, gv, us, ea, eb, cs%sponge_CSp, rcv_ml)

2699 else

2700 call apply_sponge(h, tv, dt, g, gv, us, ea, eb, cs%sponge_CSp)

2701 endif

2702 call cpu_clock_end(id_clock_sponge)

2703 if (cs%debug) then

2704 call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)

2705 call mom_thermovar_chksum("apply_sponge ", tv, g, us)

2706 endif

2707 endif ! CS%use_sponge

2708

2709 ! Apply data assimilation incremental update -oda_incupd-

2710 if (cs%use_oda_incupd .and. associated(cs%oda_incupd_CSp)) then

2711 call cpu_clock_begin(id_clock_oda_incupd)

2712 call apply_oda_incupd(h, tv, u, v, dt, g, gv, us, cs%oda_incupd_CSp)

2713 call cpu_clock_end(id_clock_oda_incupd)

2714 if (cs%debug) then

2715 call mom_state_chksum("apply_oda_incupd ", u, v, h, g, gv, us, haloshift=0)

2716 call mom_thermovar_chksum("apply_oda_incupd ", tv, g, us)

2717 endif

2718 endif ! CS%use_oda_incupd

2719

2720

2721

2722! Save the diapycnal mass fluxes as a diagnostic field.

2723 if (associated(cdp%diapyc_vel)) then

2724 !$OMP parallel do default(shared)

2725 do j=js,je

2726 do k=2,nz ; do i=is,ie

2727 cdp%diapyc_vel(i,j,k) = idt * (ea(i,j,k) - eb(i,j,k-1))

2728 enddo ; enddo

2729 do i=is,ie

2730 cdp%diapyc_vel(i,j,1) = 0.0

2731 cdp%diapyc_vel(i,j,nz+1) = 0.0

2732 enddo

2733 enddo

2734 endif

2735

2736! For momentum, it is only the net flux that homogenizes within

2737! the mixed layer. Vertical viscosity that is proportional to the

2738! mixed layer turbulence is applied elsewhere.

2739 if (cs%use_bulkmixedlayer) then

2740 if (cs%debug) then

2741 call hchksum(ea, "before net flux rearrangement ea", g%HI, unscale=gv%H_to_MKS)

2742 call hchksum(eb, "before net flux rearrangement eb", g%HI, unscale=gv%H_to_MKS)

2743 endif

2744 !$OMP parallel do default(shared) private(net_ent)

2745 do j=js,je

2746 do k=2,gv%nkml ; do i=is,ie

2747 net_ent = ea(i,j,k) - eb(i,j,k-1)

2748 ea(i,j,k) = max(net_ent, 0.0)

2749 eb(i,j,k-1) = max(-net_ent, 0.0)

2750 enddo ; enddo

2751 enddo

2752 if (cs%debug) then

2753 call hchksum(ea, "after net flux rearrangement ea", g%HI, unscale=gv%H_to_MKS)

2754 call hchksum(eb, "after net flux rearrangement eb", g%HI, unscale=gv%H_to_MKS)

2755 endif

2756 endif

2757

2758! Initialize halo regions of ea, eb, and hold to default values.

2759 !$OMP parallel do default(shared)

2760 do k=1,nz

2761 do i=is-1,ie+1

2762 hold(i,js-1,k) = gv%Angstrom_H ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0

2763 hold(i,je+1,k) = gv%Angstrom_H ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0

2764 enddo

2765 do j=js,je

2766 hold(is-1,j,k) = gv%Angstrom_H ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0

2767 hold(ie+1,j,k) = gv%Angstrom_H ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0

2768 enddo

2769 enddo

2770

2771 call cpu_clock_begin(id_clock_pass)

2772 if (g%symmetric) then ; dir_flag = to_all+omit_corners

2773 else ; dir_flag = to_west+to_south+omit_corners ; endif

2774 call create_group_pass(cs%pass_hold_eb_ea, hold, g%Domain, dir_flag, halo=1)

2775 call create_group_pass(cs%pass_hold_eb_ea, eb, g%Domain, dir_flag, halo=1)

2776 call create_group_pass(cs%pass_hold_eb_ea, ea, g%Domain, dir_flag, halo=1)

2777 call do_group_pass(cs%pass_hold_eb_ea, g%Domain)

2778 call cpu_clock_end(id_clock_pass)

2779

2780 ! Use a tridiagonal solver to determine effect of the diapycnal

2781 ! advection on velocity field. It is assumed that water leaves

2782 ! or enters the ocean with the surface velocity.

2783 if (cs%debug) then

2784 call mom_state_chksum("before u/v tridiag ", u, v, h, g, gv, us, haloshift=0)

2785 call hchksum(ea, "before u/v tridiag ea", g%HI, unscale=gv%H_to_MKS)

2786 call hchksum(eb, "before u/v tridiag eb", g%HI, unscale=gv%H_to_MKS)

2787 call hchksum(hold, "before u/v tridiag hold", g%HI, unscale=gv%H_to_MKS)

2788 endif

2789 call cpu_clock_begin(id_clock_tridiag)

2790

2791 !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)

2792 do j=js,je

2793 do i=isq,ieq

2794 if (associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,1) = u(i,j,1)

2795 hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect

2796 b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))

2797 d1(i) = hval * b1(i)

2798 u(i,j,1) = b1(i) * (hval * u(i,j,1))

2799 enddo

2800 do k=2,nz ; do i=isq,ieq

2801 if (associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,k) = u(i,j,k)

2802 c1(i,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(i)

2803 eaval = ea(i,j,k) + ea(i+1,j,k)

2804 hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect

2805 b1(i) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(i)*eaval))

2806 d1(i) = (hval + d1(i)*eaval) * b1(i)

2807 u(i,j,k) = (hval*u(i,j,k) + eaval*u(i,j,k-1))*b1(i)

2808 enddo ; enddo

2809 do k=nz-1,1,-1 ; do i=isq,ieq

2810 u(i,j,k) = u(i,j,k) + c1(i,k+1)*u(i,j,k+1)

2811 if (associated(adp%du_dt_dia)) &

2812 adp%du_dt_dia(i,j,k) = (u(i,j,k) - adp%du_dt_dia(i,j,k)) * idt

2813 enddo ; enddo

2814 if (associated(adp%du_dt_dia)) then

2815 do i=isq,ieq

2816 adp%du_dt_dia(i,j,nz) = (u(i,j,nz)-adp%du_dt_dia(i,j,nz)) * idt

2817 enddo

2818 endif

2819 enddo

2820 if (cs%debug) then

2821 call mom_state_chksum("aft 1st loop tridiag ", u, v, h, g, gv, us, haloshift=0)

2822 endif

2823 !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)

2824 do j=jsq,jeq

2825 do i=is,ie

2826 if (associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,1) = v(i,j,1)

2827 hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect

2828 b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))

2829 d1(i) = hval * b1(i)

2830 v(i,j,1) = b1(i) * (hval * v(i,j,1))

2831 enddo

2832 do k=2,nz ; do i=is,ie

2833 if (associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,k) = v(i,j,k)

2834 c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)

2835 eaval = ea(i,j,k) + ea(i,j+1,k)

2836 hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect

2837 b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))

2838 d1(i) = (hval + d1(i)*eaval) * b1(i)

2839 v(i,j,k) = (hval*v(i,j,k) + eaval*v(i,j,k-1))*b1(i)

2840 enddo ; enddo

2841 do k=nz-1,1,-1 ; do i=is,ie

2842 v(i,j,k) = v(i,j,k) + c1(i,k+1)*v(i,j,k+1)

2843 if (associated(adp%dv_dt_dia)) &

2844 adp%dv_dt_dia(i,j,k) = (v(i,j,k) - adp%dv_dt_dia(i,j,k)) * idt

2845 enddo ; enddo

2846 if (associated(adp%dv_dt_dia)) then

2847 do i=is,ie

2848 adp%dv_dt_dia(i,j,nz) = (v(i,j,nz)-adp%dv_dt_dia(i,j,nz)) * idt

2849 enddo

2850 endif

2851 enddo

2852 call cpu_clock_end(id_clock_tridiag)

2853 if (cs%debug) then

2854 call mom_state_chksum("after u/v tridiag ", u, v, h, g, gv, us, haloshift=0)

2855 endif

2856

2857 ! Diagnose the diapycnal diffusivities and other related quantities.

2858 if (cs%id_Kd_int > 0) call post_data(cs%id_Kd_int, kd_int, cs%diag)

2859 if (cs%id_Kd_heat > 0) call post_data(cs%id_Kd_heat, kd_heat, cs%diag)

2860 if (cs%id_Kd_salt > 0) call post_data(cs%id_Kd_salt, kd_salt, cs%diag)

2861

2862 if (cs%id_ea > 0) call post_data(cs%id_ea, ea, cs%diag)

2863 if (cs%id_eb > 0) call post_data(cs%id_eb, eb, cs%diag)

2864

2865 if (cs%id_dudt_dia > 0) call post_data(cs%id_dudt_dia, adp%du_dt_dia, cs%diag)

2866 if (cs%id_dvdt_dia > 0) call post_data(cs%id_dvdt_dia, adp%dv_dt_dia, cs%diag)

2867 if (cs%id_wd > 0) call post_data(cs%id_wd, cdp%diapyc_vel, cs%diag)

2868

2869 if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)

2870 if (cs%id_Tadv > 0) call post_data(cs%id_Tadv, tadv_flx, cs%diag)

2871 if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)

2872 if (cs%id_Sadv > 0) call post_data(cs%id_Sadv, sadv_flx, cs%diag)

2873

2874 ! Diagnostics for thickness-weighted vertically averaged diapycnal accelerations

2875 if (cs%id_hf_dudt_dia_2d > 0) &

2876 call post_product_sum_u(cs%id_hf_dudt_dia_2d, adp%du_dt_dia, adp%diag_hfrac_u, g, nz, cs%diag)

2877 if (cs%id_hf_dvdt_dia_2d > 0) &

2878 call post_product_sum_v(cs%id_hf_dvdt_dia_2d, adp%dv_dt_dia, adp%diag_hfrac_v, g, nz, cs%diag)

2879

2880 call disable_averaging(cs%diag)

2881

2882 if (showcalltree) call calltree_leave("layered_diabatic()")

2883

2884end subroutine layered_diabatic

2885

2886!> Returns pointers or values of members within the diabatic_CS type. For extensibility,

2887!! each returned argument is an optional argument

2888subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit, minimum_forcing_depth, &

2889 KPP_CSp, energetic_PBL_CSp, diabatic_aux_CSp, diabatic_halo, use_KPP)

2890 type(diabatic_cs), target, intent(in) :: cs !< module control structure

2891 ! All output arguments are optional

2892 type(opacity_cs), optional, pointer :: opacity_csp !< A pointer to be set to the opacity control structure

2893 type(optics_type), optional, pointer :: optics_csp !< A pointer to be set to the optics control structure

2894 type(kpp_cs), optional, pointer :: kpp_csp !< A pointer to be set to the KPP CS

2895 type(energetic_pbl_cs), optional, pointer :: energetic_pbl_csp !< A pointer to be set to the ePBL CS

2896 real, optional, intent( out) :: evap_cfl_limit !<The largest fraction of a layer that can be

2897 !! evaporated in one time-step [nondim].

2898 real, optional, intent( out) :: minimum_forcing_depth !< The smallest depth over which heat

2899 !! and freshwater fluxes are applied [H ~> m or kg m-2].

2900 type(diabatic_aux_cs), optional, pointer :: diabatic_aux_csp !< A pointer to be set to the diabatic_aux

2901 !! control structure

2902 integer, optional, intent( out) :: diabatic_halo !< The halo size where the diabatic algorithms

2903 !! assume thermodynamics properties are valid.

2904 logical, optional, intent( out) :: use_kpp !< If true, diabatic is using KPP vertical mixing

2905

2906 ! Pointers to control structures

2907 if (present(opacity_csp)) opacity_csp => cs%opacity

2908 if (present(optics_csp)) optics_csp => cs%optics

2909 if (present(kpp_csp)) kpp_csp => cs%KPP_CSp

2910 if (present(energetic_pbl_csp) .and. cs%use_energetic_PBL) energetic_pbl_csp => cs%ePBL

2911 if (present(diabatic_aux_csp)) diabatic_aux_csp => cs%diabatic_aux_CSp

2912

2913 ! Constants within diabatic_CS

2914 if (present(evap_cfl_limit)) evap_cfl_limit = cs%evap_CFL_limit

2915 if (present(minimum_forcing_depth)) minimum_forcing_depth = cs%minimum_forcing_depth

2916 if (present(diabatic_halo)) diabatic_halo = cs%halo_diabatic

2917 if (present(use_kpp)) use_kpp = cs%use_KPP

2918end subroutine extract_diabatic_member

2919

2920!> Routine called for adiabatic physics

2921subroutine adiabatic(h, tv, fluxes, dt, G, GV, US, CS)

2922 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure

2923 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure

2924 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

2925 intent(inout) :: h !< thickness [H ~> m or kg m-2]

2926 type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields

2927 type(forcing), intent(inout) :: fluxes !< boundary fluxes

2928 real, intent(in) :: dt !< time step [T ~> s]

2929 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

2930 type(diabatic_cs), pointer :: cs !< module control structure

2931

2932 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: zeros ! An array of zeros with units of [H ~> m or kg m-2]

2933

2934 zeros(:,:,:) = 0.0

2935

2936 call call_tracer_column_fns(h, h, zeros, zeros, fluxes, zeros(:,:,1), dt, g, gv, us, tv, &

2937 cs%optics, cs%tracer_flow_CSp, cs%debug)

2938

2939end subroutine adiabatic

2940

2941

2942!> This routine diagnoses tendencies from application of diabatic diffusion

2943!! using ALE algorithm. Note that layer thickness is not altered by

2944!! diabatic diffusion.

2945subroutine diagnose_diabatic_diff_tendency(tv, h, temp_old, saln_old, dt, G, GV, US, CS)

2946 type(ocean_grid_type), intent(in) :: G !< ocean grid structure

2947 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

2948 type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields

2949 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< thickness [H ~> m or kg m-2]

2950 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: temp_old !< temperature prior to diabatic

2951 !! physics [C ~> degC]

2952 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: saln_old !< salinity prior to diabatic physics [S ~> ppt]

2953 real, intent(in) :: dt !< time step [T ~> s]

2954 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

2955 type(diabatic_cs), pointer :: CS !< module control structure

2956

2957 ! Local variables

2958 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_3d ! A 3-d work array for diagnostics [various]

2959 real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array for diagnostics [various]

2960 real :: Idt ! The inverse of the timestep [T-1 ~> s-1]

2961 integer :: i, j, k, is, ie, js, je, nz

2962 logical :: do_saln_tend ! Calculate salinity-based tendency diagnostics

2963

2964 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

2965 idt = 0.0 ; if (dt > 0.0) idt = 1. / dt

2966 work_3d(:,:,:) = 0.0

2967 work_2d(:,:) = 0.0

2968

2969

2970 ! temperature tendency

2971 do k=1,nz ; do j=js,je ; do i=is,ie

2972 work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt

2973 enddo ; enddo ; enddo

2974 if (cs%id_diabatic_diff_temp_tend > 0) then

2975 call post_data(cs%id_diabatic_diff_temp_tend, work_3d, cs%diag, alt_h=h)

2976 endif

2977

2978 ! heat tendency

2979 if (cs%id_diabatic_diff_heat_tend > 0 .or. cs%id_diabatic_diff_heat_tend_2d > 0) then

2980 do k=1,nz ; do j=js,je ; do i=is,ie

2981 work_3d(i,j,k) = h(i,j,k)*gv%H_to_RZ * tv%C_p * work_3d(i,j,k)

2982 enddo ; enddo ; enddo

2983 if (cs%id_diabatic_diff_heat_tend > 0) then

2984 call post_data(cs%id_diabatic_diff_heat_tend, work_3d, cs%diag, alt_h=h)

2985 endif

2986 if (cs%id_diabatic_diff_heat_tend_2d > 0) then

2987 work_2d(:,:) = 0.0

2988 do k=1,nz ; do j=js,je ; do i=is,ie

2989 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

2990 enddo ; enddo ; enddo

2991 call post_data(cs%id_diabatic_diff_heat_tend_2d, work_2d, cs%diag)

2992 endif

2993 endif

2994

2995 ! salinity tendency

2996 do_saln_tend = cs%id_diabatic_diff_saln_tend > 0 &

2997 .or. cs%id_diabatic_diff_salt_tend > 0 &

2998 .or. cs%id_diabatic_diff_salt_tend_2d > 0

2999

3000 if (do_saln_tend) then

3001 do k=1,nz ; do j=js,je ; do i=is,ie

3002 work_3d(i,j,k) = (tv%S(i,j,k) - saln_old(i,j,k)) * idt

3003 enddo ; enddo ; enddo

3004

3005 if (cs%id_diabatic_diff_saln_tend > 0) &

3006 call post_data(cs%id_diabatic_diff_saln_tend, work_3d, cs%diag, alt_h=h)

3007

3008 ! salt tendency

3009 if (cs%id_diabatic_diff_salt_tend > 0 .or. cs%id_diabatic_diff_salt_tend_2d > 0) then

3010 do k=1,nz ; do j=js,je ; do i=is,ie

3011 work_3d(i,j,k) = h(i,j,k) * work_3d(i,j,k)

3012 enddo ; enddo ; enddo

3013 if (cs%id_diabatic_diff_salt_tend > 0) then

3014 call post_data(cs%id_diabatic_diff_salt_tend, work_3d, cs%diag, alt_h=h)

3015 endif

3016 if (cs%id_diabatic_diff_salt_tend_2d > 0) then

3017 work_2d(:,:) = 0.0

3018 do k=1,nz ; do j=js,je ; do i=is,ie

3019 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

3020 enddo ; enddo ; enddo

3021 call post_data(cs%id_diabatic_diff_salt_tend_2d, work_2d, cs%diag)

3022 endif

3023 endif

3024 endif

3025

3026end subroutine diagnose_diabatic_diff_tendency

3027

3028

3029!> This routine diagnoses tendencies from application of boundary fluxes.

3030!! These impacts are generally 3d, in particular for penetrative shortwave.

3031!! Other fluxes contribute 3d in cases when the layers vanish or are very thin,

3032!! in which case we distribute the flux into k > 1 layers.

3033subroutine diagnose_boundary_forcing_tendency(tv, h, temp_old, saln_old, h_old, &

3034 dt, G, GV, US, CS)

3035 type(ocean_grid_type), intent(in) :: G !< ocean grid structure

3036 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

3037 type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields

3038 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

3039 intent(in) :: h !< thickness after boundary flux application [H ~> m or kg m-2]

3040 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

3041 intent(in) :: temp_old !< temperature prior to boundary flux application [C ~> degC]

3042 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

3043 intent(in) :: saln_old !< salinity prior to boundary flux application [S ~> ppt]

3044 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

3045 intent(in) :: h_old !< thickness prior to boundary flux application [H ~> m or kg m-2]

3046 real, intent(in) :: dt !< time step [T ~> s]

3047 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

3048 type(diabatic_cs), pointer :: CS !< module control structure

3049

3050 ! Local variables

3051 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_3d ! A 3-d work array for diagnostics [various]

3052 real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array for diagnostics [various]

3053 real :: Idt ! The inverse of the timestep [T-1 ~> s-1]

3054 integer :: i, j, k, is, ie, js, je, nz

3055

3056 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

3057 idt = 0.0 ; if (dt > 0.0) idt = 1. / dt

3058 work_3d(:,:,:) = 0.0

3059 work_2d(:,:) = 0.0

3060

3061 ! Thickness tendency

3062 if (cs%id_boundary_forcing_h_tendency > 0) then

3063 do k=1,nz ; do j=js,je ; do i=is,ie

3064 work_3d(i,j,k) = (h(i,j,k) - h_old(i,j,k))*idt

3065 enddo ; enddo ; enddo

3066 call post_data(cs%id_boundary_forcing_h_tendency, work_3d, cs%diag, alt_h=h_old)

3067 endif

3068

3069 ! temperature tendency

3070 if (cs%id_boundary_forcing_temp_tend > 0) then

3071 do k=1,nz ; do j=js,je ; do i=is,ie

3072 work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt

3073 enddo ; enddo ; enddo

3074 call post_data(cs%id_boundary_forcing_temp_tend, work_3d, cs%diag, alt_h=h_old)

3075 endif

3076

3077 ! heat tendency

3078 if (cs%id_boundary_forcing_heat_tend > 0 .or. cs%id_boundary_forcing_heat_tend_2d > 0) then

3079 do k=1,nz ; do j=js,je ; do i=is,ie

3080 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))

3081 enddo ; enddo ; enddo

3082 if (cs%id_boundary_forcing_heat_tend > 0) then

3083 call post_data(cs%id_boundary_forcing_heat_tend, work_3d, cs%diag, alt_h=h_old)

3084 endif

3085 if (cs%id_boundary_forcing_heat_tend_2d > 0) then

3086 work_2d(:,:) = 0.0

3087 do k=1,nz ; do j=js,je ; do i=is,ie

3088 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

3089 enddo ; enddo ; enddo

3090 call post_data(cs%id_boundary_forcing_heat_tend_2d, work_2d, cs%diag)

3091 endif

3092 endif

3093

3094 ! salinity tendency

3095 if (cs%id_boundary_forcing_saln_tend > 0) then

3096 do k=1,nz ; do j=js,je ; do i=is,ie

3097 work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*idt

3098 enddo ; enddo ; enddo

3099 call post_data(cs%id_boundary_forcing_saln_tend, work_3d, cs%diag, alt_h=h_old)

3100 endif

3101

3102 ! salt tendency

3103 if (cs%id_boundary_forcing_salt_tend > 0 .or. cs%id_boundary_forcing_salt_tend_2d > 0) then

3104 do k=1,nz ; do j=js,je ; do i=is,ie

3105 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))

3106 enddo ; enddo ; enddo

3107 if (cs%id_boundary_forcing_salt_tend > 0) then

3108 call post_data(cs%id_boundary_forcing_salt_tend, work_3d, cs%diag, alt_h=h_old)

3109 endif

3110 if (cs%id_boundary_forcing_salt_tend_2d > 0) then

3111 work_2d(:,:) = 0.0

3112 do k=1,nz ; do j=js,je ; do i=is,ie

3113 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

3114 enddo ; enddo ; enddo

3115 call post_data(cs%id_boundary_forcing_salt_tend_2d, work_2d, cs%diag)

3116 endif

3117 endif

3118

3119end subroutine diagnose_boundary_forcing_tendency

3120

3121

3122!> This routine diagnoses tendencies for temperature and heat from frazil formation.

3123!! This routine is called twice from within subroutine diabatic; at start and at

3124!! end of the diabatic processes. The impacts from frazil are generally a function

3125!! of depth. Hence, when checking heat budget, be sure to remove HFSIFRAZIL from HFDS in k=1.

3126subroutine diagnose_frazil_tendency(tv, h, temp_old, dt, G, GV, US, CS)

3127 type(ocean_grid_type), intent(in) :: G !< ocean grid structure

3128 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

3129 type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields

3130 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< thickness [H ~> m or kg m-2]

3131 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: temp_old !< temperature prior to frazil formation [C ~> degC]

3132 real, intent(in) :: dt !< time step [T ~> s]

3133 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

3134 type(diabatic_cs), pointer :: CS !< module control structure

3135

3136 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_3d ! A 3-d work array for diagnostics [various]

3137 real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array for diagnostics [various]

3138 real :: Idt ! The inverse of the timestep [T-1 ~> s-1]

3139 integer :: i, j, k, is, ie, js, je, nz

3140

3141 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

3142 idt = 0.0 ; if (dt > 0.0) idt = 1. / dt

3143

3144 ! temperature tendency

3145 if (cs%id_frazil_temp_tend > 0) then

3146 do k=1,nz ; do j=js,je ; do i=is,ie

3147 work_3d(i,j,k) = idt * (tv%T(i,j,k)-temp_old(i,j,k))

3148 enddo ; enddo ; enddo

3149 call post_data(cs%id_frazil_temp_tend, work_3d, cs%diag)

3150 endif

3151

3152 ! heat tendency

3153 if (cs%id_frazil_heat_tend > 0 .or. cs%id_frazil_heat_tend_2d > 0) then

3154 do k=1,nz ; do j=js,je ; do i=is,ie

3155 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))

3156 enddo ; enddo ; enddo

3157 if (cs%id_frazil_heat_tend > 0) call post_data(cs%id_frazil_heat_tend, work_3d, cs%diag)

3158

3159 ! As a consistency check, we must have

3160 ! FRAZIL_HEAT_TENDENCY_2d = HFSIFRAZIL

3161 if (cs%id_frazil_heat_tend_2d > 0) then

3162 work_2d(:,:) = 0.0

3163 do k=1,nz ; do j=js,je ; do i=is,ie

3164 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

3165 enddo ; enddo ; enddo

3166 call post_data(cs%id_frazil_heat_tend_2d, work_2d, cs%diag)

3167 endif

3168 endif

3169

3170end subroutine diagnose_frazil_tendency

3171

3172

3173!> A simplified version of diabatic_driver_init that will allow

3174!! tracer column functions to be called without allowing any

3175!! of the diabatic processes to be used.

3176subroutine adiabatic_driver_init(Time, G, param_file, diag, CS, &

3177 tracer_flow_CSp)

3178 type(time_type), intent(in) :: time !< current model time

3179 type(ocean_grid_type), intent(in) :: g !< model grid structure

3180 type(param_file_type), intent(in) :: param_file !< the file to parse for parameter values

3181 type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output

3182 type(diabatic_cs), pointer :: cs !< module control structure

3183 type(tracer_flow_control_cs), pointer :: tracer_flow_csp !< pointer to control structure of the

3184 !! tracer flow control module

3185

3186 ! This "include" declares and sets the variable "version".

3187# include "version_variable.h"

3188 character(len=40) :: mdl = "MOM_diabatic_driver" ! This module's name.

3189

3190 if (associated(cs)) then

3191 call mom_error(warning, "adiabatic_driver_init called with an "// &

3192 "associated control structure.")

3193 return

3194 else ; allocate(cs) ; endif

3195

3196 cs%diag => diag

3197 if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp

3198

3199 ! Set default, read and log parameters

3200 call log_version(param_file, mdl, version, &

3201 "The following parameters are used for diabatic processes.")

3202

3203 ! Check for any subsidiary parameters that are inconsistent with the adiabatic mode.

3204 call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &

3205 "If true, sponges may be applied anywhere in the domain. "//&

3206 "The exact location and properties of those sponges are "//&

3207 "specified via calls to initialize_sponge and possibly "//&

3208 "set_up_sponge_field.", default=.false., do_not_log=.true.)

3209 call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", cs%use_energetic_PBL, &

3210 "If true, use an implied energetics planetary boundary "//&

3211 "layer scheme to determine the diffusivity and viscosity "//&

3212 "in the surface boundary layer.", default=.false., do_not_log=.true.)

3213 call get_param(param_file, mdl, "USE_KPP", cs%use_KPP, &

3214 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&

3215 "to calculate diffusivities and non-local transport in the OBL.", &

3216 default=.false., do_not_log=.true.)

3217

3218 if (cs%use_sponge) call mom_error(warning, &

3219 "When ADIABATIC = True, it is inconsistent to set SPONGE = True.")

3220 if (cs%use_energetic_PBL) call mom_error(warning, &

3221 "When ADIABATIC = True, it is inconsistent to set ENERGETICS_SFC_PBL = True.")

3222 if (cs%use_KPP) call mom_error(warning, &

3223 "When ADIABATIC = True, it is inconsistent to set USE_KPP = True.")

3224

3225 if (cs%use_sponge .or. cs%use_energetic_PBL .or. cs%use_KPP) &

3226 call mom_error(fatal, "adiabatic_driver_init is aborting due to inconsistent parameter settings.")

3227

3228end subroutine adiabatic_driver_init

3229

3230

3231!> This routine initializes the diabatic driver module.

3232subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, diag, &

3233 ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, &

3234 ALE_sponge_CSp, oda_incupd_CSp, int_tide_CSp)

3235 type(time_type), target :: time !< model time

3236 type(ocean_grid_type), intent(inout) :: g !< model grid structure

3237 type(verticalgrid_type), intent(in) :: gv !< model vertical grid structure

3238 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

3239 type(param_file_type), intent(in) :: param_file !< file to parse for parameter values

3240 logical, intent(in) :: usealealgorithm !< logical for whether to use ALE remapping

3241 type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output

3242 type(accel_diag_ptrs), intent(inout) :: adp !< pointers to accelerations in momentum equations,

3243 !! to enable diagnostics, like energy budgets

3244 type(cont_diag_ptrs), intent(inout) :: cdp !< pointers to terms in continuity equations

3245 type(diabatic_cs), pointer :: cs !< module control structure

3246 type(tracer_flow_control_cs), pointer :: tracer_flow_csp !< pointer to control structure of the

3247 !! tracer flow control module

3248 type(sponge_cs), pointer :: sponge_csp !< pointer to the sponge module control structure

3249 type(ale_sponge_cs), pointer :: ale_sponge_csp !< pointer to the ALE sponge module control structure

3250 type(oda_incupd_cs), pointer :: oda_incupd_csp !< pointer to the ocean data assimilation incremental

3251 !! update module control structure

3252 type(int_tide_cs), pointer :: int_tide_csp !< pointer to the internal tide structure

3253

3254 ! Local variables

3255 real :: kd ! A diffusivity used in the default for other tracer diffusivities [Z2 T-1 ~> m2 s-1]

3256 logical :: use_temperature

3257 character(len=20) :: en1, en2, en3

3258

3259 ! This "include" declares and sets the variable "version".

3260# include "version_variable.h"

3261 character(len=40) :: mdl = "MOM_diabatic_driver" ! This module's name.

3262 character(len=48) :: thickness_units

3263 character(len=20) :: brine_plume_mld_def

3264 logical :: physical_obl_scheme, do_brine_plume

3265 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nbands

3266 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke

3267 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

3268

3269 cs%initialized = .true.

3270

3271 cs%diag => diag

3272 cs%Time => time

3273

3274 if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp

3275 if (associated(sponge_csp)) cs%sponge_CSp => sponge_csp

3276 if (associated(ale_sponge_csp)) cs%ALE_sponge_CSp => ale_sponge_csp

3277 if (associated(oda_incupd_csp)) cs%oda_incupd_CSp => oda_incupd_csp

3278 if (associated(int_tide_csp)) cs%int_tide_CSp => int_tide_csp

3279

3280 cs%useALEalgorithm = usealealgorithm

3281 cs%use_bulkmixedlayer = (gv%nkml > 0)

3282

3283 ! Set default, read and log parameters

3284 call log_version(param_file, mdl, version, &

3285 "The following parameters are used for diabatic processes.", &

3286 log_to_all=.true., debugging=.true.)

3287 call get_param(param_file, mdl, "USE_LEGACY_DIABATIC_DRIVER", cs%use_legacy_diabatic, &

3288 "If true, use a legacy version of the diabatic subroutine. "//&

3289 "This is temporary and is needed to avoid change in answers.", &

3290 default=.true.)

3291 call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &

3292 "If true, sponges may be applied anywhere in the domain. "//&

3293 "The exact location and properties of those sponges are "//&

3294 "specified via calls to initialize_sponge and possibly "//&

3295 "set_up_sponge_field.", default=.false.)

3296 call get_param(param_file, mdl, "ODA_INCUPD", cs%use_oda_incupd, &

3297 "If true, oda incremental updates will be applied "//&

3298 "everywhere in the domain.", default=.false.)

3299 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &

3300 "If true, temperature and salinity are used as state "//&

3301 "variables.", default=.true.)

3302 call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", cs%use_energetic_PBL, &

3303 "If true, use an implied energetics planetary boundary "//&

3304 "layer scheme to determine the diffusivity and viscosity "//&

3305 "in the surface boundary layer.", default=.false.)

3306 call get_param(param_file, mdl, "EPBL_IS_ADDITIVE", cs%ePBL_is_additive, &

3307 "If true, the diffusivity from ePBL is added to all "//&

3308 "other diffusivities. Otherwise, the larger of kappa-shear "//&

3309 "and ePBL diffusivities are used.", default=.true.)

3310 call get_param(param_file, mdl, "PRANDTL_EPBL", cs%ePBL_Prandtl, &

3311 "The Prandtl number used by ePBL to convert vertical diffusivities into "//&

3312 "viscosities.", default=1.0, units="nondim", do_not_log=.not.cs%use_energetic_PBL)

3313 call get_param(param_file, mdl, "USE_KPP", cs%use_KPP, &

3314 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&

3315 "to calculate diffusivities and non-local transport in the OBL.", &

3316 default=.false., do_not_log=.true.)

3317 cs%use_CVMix_ddiff = cvmix_ddiff_is_used(param_file)

3318

3319 cs%use_kappa_shear = kappa_shear_is_used(param_file)

3320 cs%use_CVMix_shear = cvmix_shear_is_used(param_file)

3321

3322 if (cs%use_bulkmixedlayer) then

3323 call get_param(param_file, mdl, "ML_MIX_FIRST", cs%ML_mix_first, &

3324 "The fraction of the mixed layer mixing that is applied "//&

3325 "before interior diapycnal mixing. 0 by default.", &

3326 units="nondim", default=0.0)

3327 call get_param(param_file, mdl, "NKBL", cs%nkbl, default=2, do_not_log=.true.)

3328 else

3329 cs%ML_mix_first = 0.0

3330 endif

3331 if (use_temperature) then

3332 call get_param(param_file, mdl, "DO_GEOTHERMAL", cs%use_geothermal, &

3333 "If true, apply geothermal heating.", default=.false.)

3334 else

3335 cs%use_geothermal = .false.

3336 endif

3337 call get_param(param_file, mdl, "INTERNAL_TIDES", cs%use_int_tides, &

3338 "If true, use the code that advances a separate set of "//&

3339 "equations for the internal tide energy density.", default=.false.)

3340

3341 call get_param(param_file, mdl, "MASSLESS_MATCH_TARGETS", cs%massless_match_targets, &

3342 "If true, the temperature and salinity of massless layers "//&

3343 "are kept consistent with their target densities. "//&

3344 "Otherwise the properties of massless layers evolve "//&

3345 "diffusively to match massive neighboring layers.", &

3346 default=.true.)

3347

3348 call get_param(param_file, mdl, "AGGREGATE_FW_FORCING", cs%aggregate_FW_forcing, &

3349 "If true, the net incoming and outgoing fresh water fluxes are combined "//&

3350 "and applied as either incoming or outgoing depending on the sign of the net. "//&

3351 "If false, the net incoming fresh water flux is added to the model and "//&

3352 "thereafter the net outgoing is removed from the topmost non-vanished "//&

3353 "layers of the updated state.", default=.true.)

3354

3355 call get_param(param_file, mdl, "DEBUG", cs%debug, &

3356 "If true, write out verbose debugging data.", &

3357 default=.false., debuggingparam=.true.)

3358 call get_param(param_file, mdl, "DEBUG_CONSERVATION", cs%debugConservation, &

3359 "If true, monitor conservation and extrema.", &

3360 default=.false., debuggingparam=.true.)

3361

3362 call get_param(param_file, mdl, "DEBUG_ENERGY_REQ", cs%debug_energy_req, &

3363 "If true, debug the energy requirements.", default=.false., do_not_log=.true.)

3364 call get_param(param_file, mdl, "MIX_BOUNDARY_TRACERS", cs%mix_boundary_tracers, &

3365 "If true, mix the passive tracers in massless layers at "//&

3366 "the bottom into the interior as though a diffusivity of "//&

3367 "KD_MIN_TR were operating.", default=.true.)

3368 call get_param(param_file, mdl, "MIX_BOUNDARY_TRACER_ALE", cs%mix_boundary_tracer_ALE, &

3369 "If true and in ALE mode, mix the passive tracers in massless layers at "//&

3370 "the bottom into the interior as though a diffusivity of "//&

3371 "KD_MIN_TR were operating.", default=.false., do_not_log=.not.cs%useALEalgorithm)

3372

3373 if (cs%mix_boundary_tracers .or. cs%mix_boundary_tracer_ALE) then

3374 call get_param(param_file, mdl, "KD", kd, units="m2 s-1", default=0.0, scale=us%m2_s_to_Z2_T)

3375 call get_param(param_file, mdl, "KD_MIN_TR", cs%Kd_min_tr, &

3376 "A minimal diffusivity that should always be applied to "//&

3377 "tracers, especially in massless layers near the bottom. "//&

3378 "The default is 0.1*KD.", &

3379 units="m2 s-1", default=0.1*kd*us%Z2_T_to_m2_s, scale=gv%m2_s_to_HZ_T)

3380 call get_param(param_file, mdl, "KD_BBL_TR", cs%Kd_BBL_tr, &

3381 "A bottom boundary layer tracer diffusivity that will "//&

3382 "allow for explicitly specified bottom fluxes. The "//&

3383 "entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) "//&

3384 "over the same distance.", &

3385 units="m2 s-1", default=0., scale=gv%m2_s_to_HZ_T*(us%Z_to_m*gv%m_to_H))

3386 ! The scaling factor here is usually equivalent to GV%m2_s_to_HZ_T*GV%Z_to_H.

3387 endif

3388

3389 call get_param(param_file, mdl, "TRACER_TRIDIAG", cs%tracer_tridiag, &

3390 "If true, use the passive tracer tridiagonal solver for T and S", &

3391 default=.false.)

3392

3393 call get_param(param_file, mdl, "MINIMUM_FORCING_DEPTH", cs%minimum_forcing_depth, &

3394 "The smallest depth over which forcing can be applied. This "//&

3395 "only takes effect when near-surface layers become thin "//&

3396 "relative to this scale, in which case the forcing tendencies "//&

3397 "scaled down by distributing the forcing over this depth scale.", &

3398 units="m", default=0.001, scale=gv%m_to_H)

3399 call get_param(param_file, mdl, "EVAP_CFL_LIMIT", cs%evap_CFL_limit, &

3400 "The largest fraction of a layer than can be lost to forcing "//&

3401 "(e.g. evaporation, sea-ice formation) in one time-step. The unused "//&

3402 "mass loss is passed down through the column.", &

3403 units="nondim", default=0.8)

3404

3405 call get_param(param_file, mdl, "DO_BRINE_PLUME", do_brine_plume, &

3406 "If true, enables a brine plume parameterizations (not logged here)", &

3407 do_not_log=.true.,default=.false.)

3408 if (do_brine_plume) then

3409 call get_param(param_file, mdl, "BRINE_PLUME_MLD_DEF", brine_plume_mld_def, &

3410 "A string that determines which mixed/mixing depth is used in setting "//&

3411 "the brine plume depth, \n Valid options are MLD_003, MLD_EN1, and H_ePBL",&

3412 default='MLD_EN1')

3413 select case (trim(brine_plume_mld_def))

3414 case ('MLD_003')

3415 cs%MLD_param_003 = .true.

3416 case ('MLD_EN1')

3417 cs%MLD_param_EN1 = .true.

3418 case ('H_ePBL')

3419 cs%MLD_param_ePBL = .true.

3420 case default

3421 call mom_error(fatal,"Invalid choice for BRINE_PLUME_MLD_DEF. Valid options are"//&

3422 "MLD_003, MLD_EN1, or H_ePBL.")

3423 end select

3424

3425 endif

3426

3427 if (cs%use_energetic_PBL .and. .not.cs%useALEalgorithm) &

3428 call mom_error(fatal, "diabatic_driver_init: "//&

3429 "ENERGETICS_SFC_PBL = True is only coded to work when USE_REGRIDDING = True.")

3430

3431 ! Register all available diagnostics for this module.

3432 thickness_units = get_thickness_units(gv)

3433

3434 cs%id_ea_t = register_diag_field('ocean_model', 'ea_t', diag%axesTL, time, &

3435 'Layer (heat) entrainment from above per timestep', 'm', conversion=gv%H_to_m)

3436 cs%id_eb_t = register_diag_field('ocean_model', 'eb_t', diag%axesTL, time, &

3437 'Layer (heat) entrainment from below per timestep', 'm', conversion=gv%H_to_m)

3438 cs%id_ea_s = register_diag_field('ocean_model', 'ea_s', diag%axesTL, time, &

3439 'Layer (salt) entrainment from above per timestep', 'm', conversion=gv%H_to_m)

3440 cs%id_eb_s = register_diag_field('ocean_model', 'eb_s', diag%axesTL, time, &

3441 'Layer (salt) entrainment from below per timestep', 'm', conversion=gv%H_to_m)

3442 ! used by layer diabatic

3443 cs%id_ea = register_diag_field('ocean_model', 'ea', diag%axesTL, time, &

3444 'Layer entrainment from above per timestep', 'm', conversion=gv%H_to_m)

3445 cs%id_eb = register_diag_field('ocean_model', 'eb', diag%axesTL, time, &

3446 'Layer entrainment from below per timestep', 'm', conversion=gv%H_to_m)

3447 if (.not.cs%useALEalgorithm) then

3448 cs%id_wd = register_diag_field('ocean_model', 'wd', diag%axesTi, time, &

3449 'Diapycnal velocity', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3450 if (cs%id_wd > 0) call safe_alloc_ptr(cdp%diapyc_vel,isd,ied,jsd,jed,nz+1)

3451

3452 cs%id_dudt_dia = register_diag_field('ocean_model', 'dudt_dia', diag%axesCuL, time, &

3453 'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)

3454 cs%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, time, &

3455 'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)

3456

3457 cs%id_hf_dudt_dia_2d = register_diag_field('ocean_model', 'hf_dudt_dia_2d', diag%axesCu1, time, &

3458 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Diapycnal Mixing', &

3459 'm s-2', conversion=us%L_T2_to_m_s2)

3460 if (cs%id_hf_dudt_dia_2d > 0) call safe_alloc_ptr(adp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)

3461

3462 cs%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, time, &

3463 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', &

3464 'm s-2', conversion=us%L_T2_to_m_s2)

3465 if (cs%id_hf_dvdt_dia_2d > 0) call safe_alloc_ptr(adp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)

3466

3467 if ((cs%id_dudt_dia > 0) .or. (cs%id_hf_dudt_dia_2d > 0)) &

3468 call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)

3469 if ((cs%id_dvdt_dia > 0) .or. (cs%id_hf_dvdt_dia_2d > 0)) &

3470 call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)

3471 endif

3472

3473 if (use_temperature) then

3474 cs%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff", diag%axesTi, &

3475 time, "Diffusive diapycnal temperature flux across interfaces", &

3476 units="degC m s-1", conversion=us%C_to_degC*gv%H_to_m*us%s_to_T)

3477 if (.not.cs%useALEalgorithm) then

3478 cs%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv", diag%axesTi, &

3479 time, "Advective diapycnal temperature flux across interfaces", &

3480 units="degC m s-1", conversion=us%C_to_degC*gv%H_to_m*us%s_to_T)

3481 endif

3482 cs%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff", diag%axesTi, &

3483 time, "Diffusive diapycnal salnity flux across interfaces", &

3484 units="psu m s-1", conversion=us%S_to_ppt*gv%H_to_m*us%s_to_T)

3485 if (.not.cs%useALEalgorithm) then

3486 cs%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv", diag%axesTi, &

3487 time, "Advective diapycnal salnity flux across interfaces", &

3488 units="psu m s-1", conversion=us%S_to_ppt*gv%H_to_m*us%s_to_T)

3489 endif

3490

3491 cs%id_N2_dd = register_diag_field('ocean_model',"N2_diabatic", diag%axesTi, &

3492 time, "Squared buoyancy frequency diagnosed after diffusion applied in thermodynamic timestep.", &

3493 "s-2", conversion=us%s_to_T**2)

3494 cs%id_N2_temp_dd = register_diag_field('ocean_model',"N2_temp_diabatic", diag%axesTi, &

3495 time, "Squared buoyancy frequency due to temperature stratification diagnosed after diffusion applied "//&

3496 "in thermodynamic timestep.", "s-2", conversion=us%s_to_T**2)

3497 cs%id_N2_salt_dd = register_diag_field('ocean_model',"N2_salt_diabatic", diag%axesTi, &

3498 time, "Squared buoyancy frequency due to salinity stratification diagnosed after diffusion applied "//&

3499 "in thermodynamic timestep.", "s-2", conversion=us%s_to_T**2)

3500 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.

3501

3502 cs%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, time, &

3503 'Mixed layer depth (delta rho = 0.03)', units='m', conversion=us%Z_to_m, &

3504 cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', &

3505 cmor_standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t')

3506 cs%id_mlotstsq = register_diag_field('ocean_model', 'mlotstsq', diag%axesT1, time, &

3507 long_name='Square of Ocean Mixed Layer Thickness Defined by Sigma T', &

3508 standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t', &

3509 units='m2', conversion=us%Z_to_m**2)

3510 cs%id_MLD_0125 = register_diag_field('ocean_model', 'MLD_0125', diag%axesT1, time, &

3511 'Mixed layer depth (delta rho = 0.125)', 'm', conversion=us%Z_to_m)

3512 call get_param(param_file, mdl, "MLD_EN_VALS", cs%MLD_En_vals, &

3513 "The energy values used to compute MLDs. If not set (or all set to 0.), the "//&

3514 "default will overwrite to 25., 2500., 250000.", units='J/m2', &

3515 defaults=(/25., 2500., 250000./), scale=us%W_m2_to_RZ3_T3*us%s_to_T)

3516 write(en1,'(F10.2)') cs%MLD_En_vals(1)*us%RZ3_T3_to_W_m2*us%T_to_s

3517 write(en2,'(F10.2)') cs%MLD_En_vals(2)*us%RZ3_T3_to_W_m2*us%T_to_s

3518 write(en3,'(F10.2)') cs%MLD_En_vals(3)*us%RZ3_T3_to_W_m2*us%T_to_s

3519 cs%id_MLD_EN1 = register_diag_field('ocean_model', 'MLD_EN1', diag%axesT1, time, &

3520 'Mixed layer depth for energy value set to '//trim(en1)//' J/m2 (Energy set by 1st MLD_EN_VALS)', &

3521 units='m', conversion=us%Z_to_m)

3522 cs%id_MLD_EN2 = register_diag_field('ocean_model', 'MLD_EN2', diag%axesT1, time, &

3523 'Mixed layer depth for energy value set to '//trim(en2)//' J/m2 (Energy set by 2nd MLD_EN_VALS)', &

3524 units='m', conversion=us%Z_to_m)

3525 cs%id_MLD_EN3 = register_diag_field('ocean_model', 'MLD_EN3', diag%axesT1, time, &

3526 'Mixed layer depth for energy value set to '//trim(en3)//' J/m2 (Energy set by 3rd MLD_EN_VALS)', &

3527 units='m', conversion=us%Z_to_m)

3528 call get_param(param_file, mdl, "BMLD_EN_VALS", cs%BMLD_En_vals, &

3529 "The energy values used to compute Bottom MLDs. If not set (or all set to 0.), the "//&

3530 "default will overwrite to 25., 2500., 250000.", units='J/m2', &

3531 defaults=(/25., 2500., 250000./), scale=us%W_m2_to_RZ3_T3*us%s_to_T)

3532 write(en1,'(F10.2)') cs%BMLD_En_vals(1)*us%RZ3_T3_to_W_m2*us%T_to_s

3533 write(en2,'(F10.2)') cs%BMLD_En_vals(2)*us%RZ3_T3_to_W_m2*us%T_to_s

3534 write(en3,'(F10.2)') cs%BMLD_En_vals(3)*us%RZ3_T3_to_W_m2*us%T_to_s

3535 cs%id_BMLD_EN1 = register_diag_field('ocean_model', 'BMLD_EN1', diag%axesT1, time, &

3536 'Bottom mixed layer depth for energy value set to '//trim(en1)//' J/m2 (Energy set by 1st MLD_EN_VALS)', &

3537 units='m', conversion=us%Z_to_m)

3538 cs%id_BMLD_EN2 = register_diag_field('ocean_model', 'BMLD_EN2', diag%axesT1, time, &

3539 'Bottom mixed layer depth for energy value set to '//trim(en2)//' J/m2 (Energy set by 2nd MLD_EN_VALS)', &

3540 units='m', conversion=us%Z_to_m)

3541 cs%id_BMLD_EN3 = register_diag_field('ocean_model', 'BMLD_EN3', diag%axesT1, time, &

3542 'Bottom mixed layer depth for energy value set to '//trim(en3)//' J/m2 (Energy set by 3rd MLD_EN_VALS)', &

3543 units='m', conversion=us%Z_to_m)

3544 if ((cs%id_MLD_EN1>0).or. (cs%id_MLD_EN2>0).or. (cs%id_MLD_EN3>0).or. &

3545 (cs%id_BMLD_EN1>0).or.(cs%id_BMLD_EN2>0).or.(cs%id_BMLD_EN3>0)) then

3546 call get_param(param_file, mdl, "USE_OM4_MLD_EN_ITER", cs%use_OM4_MLD_En_iter, &

3547 "If true, uses an older set of iteration coefficients in computing the PE based "//&

3548 "surface MLD to reproduce OM4 era models. False uses an updated (general) method.",&

3549 default=.true.)

3550 endif

3551 cs%id_subMLN2 = register_diag_field('ocean_model', 'subML_N2', diag%axesT1, time, &

3552 'Squared buoyancy frequency below mixed layer', units='s-2', conversion=us%s_to_T**2)

3553 cs%id_MLD_user = register_diag_field('ocean_model', 'MLD_user', diag%axesT1, time, &

3554 'Mixed layer depth (used defined)', units='m', conversion=us%Z_to_m)

3555 if (cs%id_MLD_003 > 0) then

3556 call get_param(param_file, mdl, "HREF_FOR_MLD", cs%ref_h_mld, &

3557 "Reference depth used to calculate the potential density used to find the mixed layer depth "//&

3558 "based on a delta rho = 0.03 kg/m3.", units='m', default=0.0, scale=us%m_to_Z)

3559 cs%id_MLD_003_zr = register_diag_field('ocean_model', 'MLD_003_refZ', diag%axesT1, time, &

3560 'Depth of reference density for MLD (delta rho = 0.03)', units='m', conversion=us%Z_to_m)

3561 cs%id_MLD_003_rr = register_diag_field('ocean_model', 'MLD_003_refRho', diag%axesT1, time, &

3562 'Reference density for MLD (delta rho = 0.03)', units='kg/m3', conversion=us%R_to_kg_m3)

3563 endif

3564 endif

3565

3566 call kdwork_init(time, g,gv,us,diag,cs%VBF,cs%Use_KdWork_diag)

3567 if (cs%Use_KdWork_diag.and.(.not.usealealgorithm)) &

3568 call mom_error(warning,"The KdWork diagnostics are not fully implemented for use in layer mode.")

3569 if (cs%Use_KdWork_diag.and.(cs%use_legacy_diabatic)) &

3570 call mom_error(warning,"The KdWork diagnostics are only approximate with the legacy diabatic driver.")

3571

3572 call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", cs%MLDdensityDifference, &

3573 "The density difference used to determine a diagnostic mixed "//&

3574 "layer depth, MLD_user, following the definition of Levitus 1982. "//&

3575 "The MLD is the depth at which the density is larger than the "//&

3576 "surface density by the specified amount.", &

3577 units='kg/m3', default=0.1, scale=us%kg_m3_to_R)

3578 call get_param(param_file, mdl, "DIAG_DEPTH_SUBML_N2", cs%dz_subML_N2, &

3579 "The distance over which to calculate a diagnostic of the "//&

3580 "stratification at the base of the mixed layer.", &

3581 units='m', default=50.0, scale=us%m_to_Z)

3582

3583 ! diagnostics for values prior to diabatic and prior to ALE

3584 cs%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, time, &

3585 'Zonal velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)

3586 cs%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, time, &

3587 'Meridional velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)

3588 cs%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, time, &

3589 'Layer Thickness before diabatic forcing', &

3590 trim(thickness_units), conversion=gv%H_to_MKS, v_extensive=.true.)

3591 cs%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, time, &

3592 'Interface Heights before diabatic forcing', 'm', conversion=us%Z_to_m)

3593 if (use_temperature) then

3594 cs%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, time, &

3595 'Potential Temperature', 'degC', conversion=us%C_to_degC)

3596 cs%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, time, &

3597 'Salinity', 'PSU', conversion=us%S_to_ppt)

3598 endif

3599

3600 cs%id_Kd_int = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, time, &

3601 'Total diapycnal diffusivity at interfaces', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)

3602 if (cs%use_energetic_PBL) then

3603 cs%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, time, &

3604 'ePBL diapycnal diffusivity at interfaces', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)

3605 endif

3606

3607 cs%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, time, &

3608 'Total diapycnal diffusivity for heat at interfaces', 'm2 s-1', conversion=gv%HZ_T_to_m2_s, &

3609 cmor_field_name='difvho', &

3610 cmor_standard_name='ocean_vertical_heat_diffusivity', &

3611 cmor_long_name='Ocean vertical heat diffusivity')

3612 cs%id_Kd_salt = register_diag_field('ocean_model', 'Kd_salt', diag%axesTi, time, &

3613 'Total diapycnal diffusivity for salt at interfaces', 'm2 s-1', conversion=gv%HZ_T_to_m2_s, &

3614 cmor_field_name='difvso', &

3615 cmor_standard_name='ocean_vertical_salt_diffusivity', &

3616 cmor_long_name='Ocean vertical salt diffusivity')

3617

3618 ! CS%useKPP is set to True if KPP-scheme is to be used, False otherwise.

3619 ! KPP_init() allocated CS%KPP_Csp and also sets CS%KPPisPassive

3620 cs%useKPP = kpp_init(param_file, g, gv, us, diag, time, cs%KPP_CSp, passive=cs%KPPisPassive)

3621

3622 ! Diagnostics for tendencies of temperature and salinity due to diabatic processes,

3623 ! available only for ALE algorithm.

3624 ! Diagnostics for tendencies of temperature and heat due to frazil

3625 cs%id_diabatic_diff_h = register_diag_field('ocean_model', 'diabatic_diff_h', diag%axesTL, time, &

3626 'Cell thickness used during diabatic diffusion', &

3627 thickness_units, conversion=gv%H_to_MKS, v_extensive=.true.)

3628

3629 if (cs%useALEalgorithm) then

3630 cs%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', &

3631 'diabatic_diff_temp_tendency', diag%axesTL, time, &

3632 'Diabatic diffusion temperature tendency', 'degC s-1', conversion=us%C_to_degC*us%s_to_T)

3633 if (cs%id_diabatic_diff_temp_tend > 0) then

3634 cs%diabatic_diff_tendency_diag = .true.

3635 endif

3636

3637 cs%id_diabatic_diff_saln_tend = register_diag_field('ocean_model',&

3638 'diabatic_diff_saln_tendency', diag%axesTL, time, &

3639 'Diabatic diffusion salinity tendency', 'psu s-1', conversion=us%S_to_ppt*us%s_to_T)

3640 if (cs%id_diabatic_diff_saln_tend > 0) then

3641 cs%diabatic_diff_tendency_diag = .true.

3642 endif

3643

3644 cs%id_diabatic_diff_heat_tend = register_diag_field('ocean_model', &

3645 'diabatic_heat_tendency', diag%axesTL, time, &

3646 'Diabatic diffusion heat tendency', &

3647 'W m-2', conversion=us%QRZ_T_to_W_m2, cmor_field_name='opottempdiff', &

3648 cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'// &

3649 'due_to_parameterized_dianeutral_mixing', &

3650 cmor_long_name='Tendency of sea water potential temperature expressed as heat content '// &

3651 'due to parameterized dianeutral mixing', &

3652 v_extensive=.true.)

3653 if (cs%id_diabatic_diff_heat_tend > 0) then

3654 cs%diabatic_diff_tendency_diag = .true.

3655 endif

3656

3657 cs%id_diabatic_diff_salt_tend = register_diag_field('ocean_model', &

3658 'diabatic_salt_tendency', diag%axesTL, time, &

3659 'Diabatic diffusion of salt tendency', &

3660 'kg m-2 s-1', conversion=us%S_to_ppt*0.001*gv%H_to_RZ*us%RZ_T_to_kg_m2s, &

3661 cmor_field_name='osaltdiff', &

3662 cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &

3663 'due_to_parameterized_dianeutral_mixing', &

3664 cmor_long_name='Tendency of sea water salinity expressed as salt content '// &

3665 'due to parameterized dianeutral mixing', &

3666 v_extensive=.true.)

3667 if (cs%id_diabatic_diff_salt_tend > 0) then

3668 cs%diabatic_diff_tendency_diag = .true.

3669 endif

3670

3671 ! This diagnostic should equal to roundoff if all is working well.

3672 cs%id_diabatic_diff_heat_tend_2d = register_diag_field('ocean_model', &

3673 'diabatic_heat_tendency_2d', diag%axesT1, time, &

3674 'Depth integrated diabatic diffusion heat tendency', &

3675 'W m-2', conversion=us%QRZ_T_to_W_m2, cmor_field_name='opottempdiff_2d', &

3676 cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'//&

3677 'due_to_parameterized_dianeutral_mixing_depth_integrated', &

3678 cmor_long_name='Tendency of sea water potential temperature expressed as heat content '//&

3679 'due to parameterized dianeutral mixing depth integrated')

3680 if (cs%id_diabatic_diff_heat_tend_2d > 0) then

3681 cs%diabatic_diff_tendency_diag = .true.

3682 endif

3683

3684 ! This diagnostic should equal to roundoff if all is working well.

3685 cs%id_diabatic_diff_salt_tend_2d = register_diag_field('ocean_model', &

3686 'diabatic_salt_tendency_2d', diag%axesT1, time, &

3687 'Depth integrated diabatic diffusion salt tendency', &

3688 'kg m-2 s-1', conversion=us%S_to_ppt*0.001*gv%H_to_RZ*us%RZ_T_to_kg_m2s, &

3689 cmor_field_name='osaltdiff_2d', &

3690 cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &

3691 'due_to_parameterized_dianeutral_mixing_depth_integrated', &

3692 cmor_long_name='Tendency of sea water salinity expressed as salt content '// &

3693 'due to parameterized dianeutral mixing depth integrated')

3694 if (cs%id_diabatic_diff_salt_tend_2d > 0) then

3695 cs%diabatic_diff_tendency_diag = .true.

3696 endif

3697

3698 ! Diagnostics for tendencies of thickness temperature and salinity due to boundary forcing,

3699 ! available only for ALE algorithm.

3700 ! Diagnostics for tendencies of temperature and heat due to frazil

3701 cs%id_boundary_forcing_h = register_diag_field('ocean_model', 'boundary_forcing_h', diag%axesTL, time, &

3702 'Cell thickness after applying boundary forcing', &

3703 thickness_units, conversion=gv%H_to_MKS, v_extensive=.true.)

3704 cs%id_boundary_forcing_h_tendency = register_diag_field('ocean_model', &

3705 'boundary_forcing_h_tendency', diag%axesTL, time, &

3706 'Cell thickness tendency due to boundary forcing', &

3707 trim(thickness_units)//" s-1", conversion=gv%H_to_MKS*us%s_to_T, v_extensive=.true.)

3708 if (cs%id_boundary_forcing_h_tendency > 0) then

3709 cs%boundary_forcing_tendency_diag = .true.

3710 endif

3711

3712 cs%id_boundary_forcing_temp_tend = register_diag_field('ocean_model',&

3713 'boundary_forcing_temp_tendency', diag%axesTL, time, &

3714 'Boundary forcing temperature tendency', 'degC s-1', conversion=us%C_to_degC*us%s_to_T)

3715 if (cs%id_boundary_forcing_temp_tend > 0) then

3716 cs%boundary_forcing_tendency_diag = .true.

3717 endif

3718

3719 cs%id_boundary_forcing_saln_tend = register_diag_field('ocean_model',&

3720 'boundary_forcing_saln_tendency', diag%axesTL, time, &

3721 'Boundary forcing saln tendency', 'psu s-1', conversion=us%S_to_ppt*us%s_to_T)

3722 if (cs%id_boundary_forcing_saln_tend > 0) then

3723 cs%boundary_forcing_tendency_diag = .true.

3724 endif

3725

3726 cs%id_boundary_forcing_heat_tend = register_diag_field('ocean_model',&

3727 'boundary_forcing_heat_tendency', diag%axesTL, time, &

3728 'Boundary forcing heat tendency', &

3729 'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive=.true.)

3730 if (cs%id_boundary_forcing_heat_tend > 0) then

3731 cs%boundary_forcing_tendency_diag = .true.

3732 endif

3733

3734 cs%id_boundary_forcing_salt_tend = register_diag_field('ocean_model',&

3735 'boundary_forcing_salt_tendency', diag%axesTL, time, &

3736 'Boundary forcing salt tendency', &

3737 'kg m-2 s-1', conversion=us%S_to_ppt*0.001*gv%H_to_RZ*us%RZ_T_to_kg_m2s, &

3738 v_extensive = .true.)

3739 if (cs%id_boundary_forcing_salt_tend > 0) then

3740 cs%boundary_forcing_tendency_diag = .true.

3741 endif

3742

3743 ! This diagnostic should equal to surface heat flux if all is working well.

3744 cs%id_boundary_forcing_heat_tend_2d = register_diag_field('ocean_model',&

3745 'boundary_forcing_heat_tendency_2d', diag%axesT1, time, &

3746 'Depth integrated boundary forcing of ocean heat', &

3747 'W m-2', conversion=us%QRZ_T_to_W_m2)

3748 if (cs%id_boundary_forcing_heat_tend_2d > 0) then

3749 cs%boundary_forcing_tendency_diag = .true.

3750 endif

3751

3752 ! This diagnostic should equal to surface salt flux if all is working well.

3753 cs%id_boundary_forcing_salt_tend_2d = register_diag_field('ocean_model',&

3754 'boundary_forcing_salt_tendency_2d', diag%axesT1, time, &

3755 'Depth integrated boundary forcing of ocean salt', &

3756 'kg m-2 s-1', conversion=us%S_to_ppt*0.001*gv%H_to_RZ*us%RZ_T_to_kg_m2s)

3757 if (cs%id_boundary_forcing_salt_tend_2d > 0) then

3758 cs%boundary_forcing_tendency_diag = .true.

3759 endif

3760 endif

3761

3762 ! diagnostics for tendencies of temp and heat due to frazil

3763 cs%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, time, &

3764 long_name='Cell Thickness', standard_name='cell_thickness', &

3765 units=thickness_units, conversion=gv%H_to_MKS, v_extensive=.true.)

3766

3767 ! diagnostic for tendency of temp due to frazil

3768 cs%id_frazil_temp_tend = register_diag_field('ocean_model',&

3769 'frazil_temp_tendency', diag%axesTL, time, &

3770 'Temperature tendency due to frazil formation', 'degC s-1', conversion=us%C_to_degC*us%s_to_T)

3771 if (cs%id_frazil_temp_tend > 0) then

3772 cs%frazil_tendency_diag = .true.

3773 endif

3774

3775 ! diagnostic for tendency of heat due to frazil

3776 cs%id_frazil_heat_tend = register_diag_field('ocean_model',&

3777 'frazil_heat_tendency', diag%axesTL, time, &

3778 'Heat tendency due to frazil formation', &

3779 'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive=.true.)

3780 if (cs%id_frazil_heat_tend > 0) then

3781 cs%frazil_tendency_diag = .true.

3782 endif

3783

3784 ! If all is working properly, this diagnostic should equal to hfsifrazil.

3785 cs%id_frazil_heat_tend_2d = register_diag_field('ocean_model',&

3786 'frazil_heat_tendency_2d', diag%axesT1, time, &

3787 'Depth integrated heat tendency due to frazil formation', &

3788 'W m-2', conversion=us%QRZ_T_to_W_m2)

3789 if (cs%id_frazil_heat_tend_2d > 0) then

3790 cs%frazil_tendency_diag = .true.

3791 endif

3792

3793 ! CS%use_CVMix_conv is set to True if CVMix convection will be used, otherwise it is False.

3794 cs%use_CVMix_conv = cvmix_conv_init(time, g, gv, us, param_file, diag, cs%CVMix_conv)

3795

3796 call entrain_diffusive_init(time, g, gv, us, param_file, diag, cs%entrain_diffusive, &

3797 just_read_params=cs%useALEalgorithm)

3798

3799 ! initialize the geothermal heating module

3800 if (cs%use_geothermal) &

3801 call geothermal_init(time, g, gv, us, param_file, diag, cs%geothermal, usealealgorithm)

3802

3803 ! initialize module for internal tide induced mixing

3804 if (cs%use_int_tides) then

3805 call int_tide_input_init(time, g, gv, us, param_file, diag, cs%int_tide_input_CSp, &

3806 cs%int_tide_input)

3807 call internal_tides_init(time, g, gv, us, param_file, diag, int_tide_csp)

3808 endif

3809

3810 !if (associated(int_tide_CSp)) CS%int_tide_CSp => int_tide_CSp

3811

3812 physical_obl_scheme = (cs%use_bulkmixedlayer .or. cs%use_KPP .or. cs%use_energetic_PBL)

3813 ! initialize module for setting diffusivities

3814 call set_diffusivity_init(time, g, gv, us, param_file, diag, cs%set_diff_CSp, cs%int_tide_CSp, &

3815 halo_ts=cs%halo_TS_diff, double_diffuse=cs%double_diffuse, &

3816 physical_obl_scheme=physical_obl_scheme)

3817

3818 cs%halo_diabatic = cs%halo_TS_diff

3819 if (cs%use_int_tides) cs%halo_diabatic = max(cs%halo_TS_diff, 2)

3820

3821 if (cs%useKPP .and. (cs%double_diffuse .and. .not.cs%use_CVMix_ddiff)) &

3822 call mom_error(fatal, 'diabatic_driver_init: DOUBLE_DIFFUSION (old method) does not work '//&

3823 'with KPP. Please set DOUBLE_DIFFUSION=False and USE_CVMIX_DDIFF=True.')

3824

3825 ! set up the clocks for this module

3826 id_clock_entrain = cpu_clock_id('(Ocean diabatic entrain)', grain=clock_module)

3827 if (cs%use_bulkmixedlayer) &

3828 id_clock_mixedlayer = cpu_clock_id('(Ocean mixed layer)', grain=clock_module)

3829 id_clock_remap = cpu_clock_id('(Ocean vert remap)', grain=clock_module)

3830 if (cs%use_geothermal) &

3831 id_clock_geothermal = cpu_clock_id('(Ocean geothermal)', grain=clock_routine)

3832 id_clock_set_diffusivity = cpu_clock_id('(Ocean set_diffusivity)', grain=clock_module)

3833 id_clock_kpp = cpu_clock_id('(Ocean KPP)', grain=clock_module)

3834 id_clock_tracers = cpu_clock_id('(Ocean tracer_columns)', grain=clock_module_driver+5)

3835 if (cs%use_sponge) &

3836 id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=clock_module)

3837 if (cs%use_oda_incupd) &

3838 id_clock_oda_incupd = cpu_clock_id('(Ocean inc. update data assimilation)', grain=clock_module)

3839 id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=clock_routine)

3840 id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=clock_routine)

3841 id_clock_differential_diff = -1 ; if (cs%double_diffuse .and. .not.cs%use_CVMix_ddiff) &

3842 id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=clock_routine)

3843

3844 ! initialize the auxiliary diabatic driver module

3845 call diabatic_aux_init(time, g, gv, us, param_file, diag, cs%diabatic_aux_CSp, &

3846 cs%useALEalgorithm, cs%use_energetic_PBL)

3847

3848 ! initialize the boundary layer modules

3849 if (cs%use_bulkmixedlayer) &

3850 call bulkmixedlayer_init(time, g, gv, us, param_file, diag, cs%bulkmixedlayer)

3851 if (cs%use_energetic_PBL) &

3852 call energetic_pbl_init(time, g, gv, us, param_file, diag, cs%ePBL)

3853

3854 call regularize_layers_init(time, g, gv, param_file, diag, cs%regularize_layers)

3855

3856 if (cs%debug_energy_req) &

3857 call diapyc_energy_req_init(time, g, gv, us, param_file, diag, cs%diapyc_en_rec_CSp)

3858

3859 ! obtain information about the number of bands for penetrative shortwave

3860 if (use_temperature) then

3861 call get_param(param_file, mdl, "PEN_SW_NBANDS", nbands, default=1)

3862 if (nbands > 0) then

3863 allocate(cs%optics)

3864 call opacity_init(time, g, gv, us, param_file, diag, cs%opacity, cs%optics)

3865 endif

3866 endif

3867

3868 ! Initialize the diagnostic grid storage

3869 call diag_grid_storage_init(cs%diag_grids_prev, g, gv, diag)

3870

3871end subroutine diabatic_driver_init

3872

3873!> Routine to register restarts, pass-through to children modules

3874subroutine register_diabatic_restarts(G, GV, US, param_file, int_tide_CSp, restart_CSp, CS)

3875 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure

3876 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure

3877 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

3878 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters

3879 type(int_tide_cs), pointer :: int_tide_csp !< Internal tide control structure

3880 type(mom_restart_cs), pointer :: restart_csp !< MOM restart control structure

3881 type(diabatic_cs), pointer :: cs !< module control structure

3882

3883 logical :: use_int_tides

3884

3885 if (associated(cs)) then

3886 call mom_error(warning, "diabatic_driver_init called with an "// &

3887 "associated control structure.")

3888 return

3889 else

3890 allocate(cs)

3891 endif

3892

3893 use_int_tides=.false.

3894

3895 call read_param(param_file, "INTERNAL_TIDES", use_int_tides)

3896

3897 if (use_int_tides) then

3898 call register_int_tide_restarts(g, gv, us, param_file, int_tide_csp, restart_csp)

3899 endif

3900

3901 call register_kpp_restarts(g, param_file, restart_csp, cs%KPP_CSp)

3902

3903end subroutine register_diabatic_restarts

3904

3905!> Routine to close the diabatic driver module

3906subroutine diabatic_driver_end(CS)

3907 type(diabatic_cs), intent(inout) :: cs !< module control structure

3908

3909 if (associated(cs%VBF)) then

3910 call kdwork_end(cs%VBF)

3911 endif

3912

3913 if (associated(cs%optics)) then

3914 call opacity_end(cs%opacity, cs%optics)

3915 deallocate(cs%optics)

3916 endif

3917

3918 if (cs%debug_energy_req) &

3919 call diapyc_energy_req_end(cs%diapyc_en_rec_CSp)

3920

3921 if (cs%use_energetic_PBL) &

3922 call energetic_pbl_end(cs%ePBL)

3923

3924 call diabatic_aux_end(cs%diabatic_aux_CSp)

3925

3926 call set_diffusivity_end(cs%set_diff_CSp)

3927

3928 deallocate(cs%set_diff_CSp)

3929

3930 if (cs%use_geothermal) &

3931 call geothermal_end(cs%geothermal)

3932

3933 if (cs%useKPP) &

3934 call kpp_end(cs%KPP_CSp)

3935

3936 ! GMM, the following is commented out because arrays in

3937 ! CS%diag_grids_prev are neither pointers or allocatables

3938 ! and, therefore, cannot be deallocated.

3939

3940 !call diag_grid_storage_end(CS%diag_grids_prev)

3941end subroutine diabatic_driver_end

3942

3943

3944!> \namespace mom_diabatic_driver

3945!!

3946!! By Robert Hallberg, Alistair Adcroft, and Stephen Griffies

3947!!

3948!! This program contains the subroutine that, along with the

3949!! subroutines that it calls, implements diapycnal mass and momentum

3950!! fluxes and a bulk mixed layer. The diapycnal diffusion can be

3951!! used without the bulk mixed layer.

3952!!

3953!! \section section_diabatic Outline of MOM diabatic

3954!!

3955!! * diabatic first determines the (diffusive) diapycnal mass fluxes

3956!! based on the convergence of the buoyancy fluxes within each layer.

3957!!

3958!! * The dual-stream entrainment scheme of MacDougall and Dewar (JPO,

3959!! 1997) is used for combined diapycnal advection and diffusion,

3960!! calculated implicitly and potentially with the Richardson number

3961!! dependent mixing, as described by Hallberg (MWR, 2000).

3962!!

3963!! * Diapycnal advection is the residual of diapycnal diffusion,

3964!! so the fully implicit upwind differencing scheme that is used is

3965!! entirely appropriate.

3966!!

3967!! * The downward buoyancy flux in each layer is determined from

3968!! an implicit calculation based on the previously

3969!! calculated flux of the layer above and an estimated flux in the

3970!! layer below. This flux is subject to the following conditions:

3971!! (1) the flux in the top and bottom layers are set by the boundary

3972!! conditions, and (2) no layer may be driven below a minimal thickness.

3973!! If there is a bulk mixed layer, the buffer layer is treated

3974!! as a fixed density layer with vanishingly small diffusivity.

3975!!

3976!! diabatic takes 5 arguments: the two velocities (u and v), the

3977!! thicknesses (h), a structure containing the forcing fields, and

3978!! the length of time over which to act (dt). The velocities and

3979!! thickness are taken as inputs and modified within the subroutine.

3980!! There is no limit on the time step.

3981

3982end module mom_diabatic_driver