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

191 !>@{ Diagnostic IDs

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

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

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

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

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

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

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

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

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

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

202

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

204 integer :: id_wd = -1

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

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

207

208 ! diagnostic for fields prior to applying diapycnal physics

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

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

211

212 integer :: id_diabatic_diff_temp_tend = -1

213 integer :: id_diabatic_diff_saln_tend = -1

214 integer :: id_diabatic_diff_heat_tend = -1

215 integer :: id_diabatic_diff_salt_tend = -1

216 integer :: id_diabatic_diff_heat_tend_2d = -1

217 integer :: id_diabatic_diff_salt_tend_2d = -1

218 integer :: id_diabatic_diff_h = -1

219

220 integer :: id_boundary_forcing_h = -1

221 integer :: id_boundary_forcing_h_tendency = -1

222 integer :: id_boundary_forcing_temp_tend = -1

223 integer :: id_boundary_forcing_saln_tend = -1

224 integer :: id_boundary_forcing_heat_tend = -1

225 integer :: id_boundary_forcing_salt_tend = -1

226 integer :: id_boundary_forcing_heat_tend_2d = -1

227 integer :: id_boundary_forcing_salt_tend_2d = -1

228

229 integer :: id_frazil_h = -1

230 integer :: id_frazil_temp_tend = -1

231 integer :: id_frazil_heat_tend = -1

232 integer :: id_frazil_heat_tend_2d = -1

233 !>@}

234

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

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

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

238

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

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

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

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

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

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

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

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

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

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

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

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

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

252

253

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

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

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

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

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

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

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

261

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

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

264

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

266end type diabatic_cs

267

268!>@{ clock ids

269integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity

270integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge

271integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap

272integer :: id_clock_kpp, id_clock_oda_incupd

273!>@}

274

275contains

276

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

278!! accompanying diapycnal advection of momentum and tracers.

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

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

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

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

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

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

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

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

287 !! unused have NULL ptrs

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

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

290 !! unused fields have NULL ptrs

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

292 !! BBL properties and related fields

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

294 !! equations, to enable the later derived

295 !! diagnostics, like energy budgets

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

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

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

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

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

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

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

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

304

305 ! local variables

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

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

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

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

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

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

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

313 real :: h_to_rl2_t2 ! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]

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

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

316

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

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

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

320

321 if (gv%ke == 1) return

322

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

324

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

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

327

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

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

330

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

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

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

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

335

336 showcalltree = calltree_showquery()

337

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

339 ! these are mostly for debugging purposes.

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

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

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

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

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

345 if (cs%id_e_predia > 0) then

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

347 call post_data(cs%id_e_predia, eta, cs%diag)

348 endif

349

350 ! Save a copy of the initial state if stochastic perturbations are active.

351 if (stoch_cs%do_sppt) then

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

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

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

355 endif

356

357 if (cs%debug) then

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

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

360 endif

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

362

363 if (cs%debug_energy_req) &

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

365

366 call cpu_clock_begin(id_clock_set_diffusivity)

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

368 call cpu_clock_end(id_clock_set_diffusivity)

369

370 ! Frazil formation keeps the temperature above the freezing point.

371 ! make_frazil is deliberately called at both the beginning and at

372 ! the end of the diabatic processes.

373 if (associated(tv%T) .AND. associated(tv%frazil)) then

374 ! For frazil diagnostic, the first call covers the first half of the time step

375 call enable_averages(0.5*dt, time_end - real_to_time(0.5*dt, unscale=us%T_to_s), cs%diag)

376 if (cs%frazil_tendency_diag) then

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

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

379 enddo ; enddo ; enddo

380 endif

381

382 if (associated(fluxes%p_surf_full)) then

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

384 else

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

386 endif

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

388

389 if (cs%frazil_tendency_diag) then

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

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

392 endif

393 call disable_averaging(cs%diag)

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

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

396

397 if (cs%use_int_tides) then

398 ! This block provides an interface for the unresolved low-mode internal tide module.

399 call set_int_tide_input(u, v, h, tv, fluxes, cs%int_tide_input, dt, g, gv, us, &

400 cs%int_tide_input_CSp)

401

402 call propagate_int_tide(h, tv, cs%int_tide_input%Nb, cs%int_tide_input%Rho_bot, dt, &

403 g, gv, us, cs%int_tide_input_CSp, cs%int_tide_CSp)

404

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

406 endif ! end CS%use_int_tides

407

408 if (cs%useALEalgorithm .and. cs%use_legacy_diabatic) then

409 call diabatic_ale_legacy(u, v, h, tv, bld, fluxes, visc, adp, cdp, dt, time_end, &

410 g, gv, us, cs, stoch_cs, waves)

411 elseif (cs%useALEalgorithm) then

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

413 g, gv, us, cs, stoch_cs, waves)

414 else

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

416 g, gv, us, cs, waves)

417 endif

418

419 call cpu_clock_begin(id_clock_pass)

420 if (associated(visc%sfc_buoy_flx)) &

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

422 if (associated(visc%h_ML)) &

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

424 if (associated(visc%MLD)) &

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

426 if (associated(visc%Kv_shear)) &

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

428 call cpu_clock_end(id_clock_pass)

429

430 call disable_averaging(cs%diag)

431 ! Frazil formation keeps temperature above the freezing point.

432 ! make_frazil is deliberately called at both the beginning and at

433 ! the end of the diabatic processes.

434 if (associated(tv%T) .AND. associated(tv%frazil)) then

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

436 if (cs%frazil_tendency_diag) then

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

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

439 enddo ; enddo ; enddo

440 endif

441

442 if (associated(fluxes%p_surf_full)) then

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

444 else

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

446 endif

447

448 if (cs%frazil_tendency_diag) then

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

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

451 endif

452

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

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

455 call disable_averaging(cs%diag)

456

457 endif ! endif for frazil

458

459 if (stoch_cs%do_sppt) then

460 ! perturb diabatic tendencies.

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

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

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

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

465 enddo ; enddo ; enddo

466 ! now that we have updated thickness and salinity, calculate freeing point

467 h_to_rl2_t2 = gv%H_to_RZ * gv%g_Earth

468 do j=js,je

469 ps(:) = 0.0

470 if (associated(fluxes%p_surf)) then

471 do i=is,ie

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

473 enddo

474 endif

475

476 do i=is,ie

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

478 enddo

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

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

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

482 enddo ; enddo

483 do k=1,nz

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

485 tv%eqn_of_state)

486 do i=is,ie

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

488 enddo

489 enddo

490 enddo

491

492 deallocate(h_in, t_in, s_in)

493 endif

494

495 ! Diagnose mixed layer depths.

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

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

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

499 cs%ref_h_mld, cs%id_MLD_003_zr, cs%id_MLD_003_rr, &

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

501 endif

502 if (cs%id_MLD_0125 > 0) then

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

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

505 endif

506 if (cs%id_MLD_user > 0) then

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

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

509 endif

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

511 ! Surface Mixed Layer diagnostic

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

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

514 endif

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

516 ! Bottom Mixed Layer diagnostic

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

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

519 endif

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

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

522

523 call disable_averaging(cs%diag)

524

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

526

527end subroutine diabatic

528

529

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

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

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

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

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

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

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

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

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

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

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

541 !! unused have NULL ptrs

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

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

544 !! unused fields have NULL ptrs

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

546 !! BBL properties and related fields

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

548 !! equations, to enable the later derived

549 !! diagnostics, like energy budgets

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

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

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

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

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

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

556

557 ! local variables

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

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

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

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

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

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

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

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

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

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

568

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

589

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

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

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

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

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

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

596

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

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

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

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

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

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

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

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

605

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

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

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

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

610

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

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

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

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

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

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

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

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

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

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

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

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

623

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

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

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

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

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

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

630

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

632

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

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

635

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

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

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

639 h_neglect = gv%H_subroundoff

640

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

642 g_rho0 = gv%g_Earth_Z_T2 / gv%H_to_RZ

643 h_to_pres = gv%H_to_RZ * gv%g_Earth

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

645

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

647

648 showcalltree = calltree_showquery()

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

650

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

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

653

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

655

656 if (cs%use_geothermal) then

657 call cpu_clock_begin(id_clock_geothermal)

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

659 call cpu_clock_end(id_clock_geothermal)

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

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

662 endif

663

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

665 ! for vertical remapping may need to be regenerated.

666 call diag_update_remap_grids(cs%diag)

667

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

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

670 ! biological properties and layer thicknesses.

671 if (associated(cs%optics)) &

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

673

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

675

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

677 if (cs%use_geothermal) then

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

679 else

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

681 endif

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

683 endif

684

685 call cpu_clock_begin(id_clock_set_diffusivity)

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

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

688 if (cs%debug) &

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

690 if (cs%double_diffuse) then

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

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

693 else

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

695 cs%set_diff_CSp, cs%VBF)

696 endif

697 call cpu_clock_end(id_clock_set_diffusivity)

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

699

700

701 if (cs%debug) then

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

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

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

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

706 endif

707

708 ! Set diffusivities for heat and salt separately

709 if (cs%useKPP) then

710 ! Add contribution from double diffusion

711 if (cs%double_diffuse) then

712 !$OMP parallel do default(shared)

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

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

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

716 enddo ; enddo ; enddo

717 else

718 !$OMP parallel do default(shared)

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

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

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

722 enddo ; enddo ; enddo

723 endif

724

725 if (cs%debug) then

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

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

728 endif

729

730 call cpu_clock_begin(id_clock_kpp)

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

732

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

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

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

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

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

738 kpp_temp_flux(:,:) = 0.

739 kpp_salt_flux(:,:) = 0.

740

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

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

743 ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux

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

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

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

747 kpp_buoy_flux, kpp_temp_flux, kpp_salt_flux)

748

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

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

751

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

753 if ( associated(fluxes%lamult) ) then

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

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

756

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

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

759 else

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

761 u_star, kpp_buoy_flux, waves=waves)

762

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

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

765 endif

766

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

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

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

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

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

772

773 if (.not.cs%KPPisPassive) then

774 !$OMP parallel do default(shared)

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

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

777 enddo ; enddo ; enddo

778 if (cs%double_diffuse) then

779 !$OMP parallel do default(shared)

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

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

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

783 enddo ; enddo ; enddo

784 endif

785 endif ! not passive

786

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

788 if (cs%debug) then

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

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

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

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

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

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

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

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

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

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

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

800 endif

801 ! Apply non-local transport of heat and salt

802 ! Changes: tv%T, tv%S

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

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

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

806 dt, tv%tr_S, tv%S)

807 call cpu_clock_end(id_clock_kpp)

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

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

810

811 if (cs%debug) then

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

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

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

815 endif

816 endif ! endif for KPP

817

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

819 ! Changes: tv%T, tv%S

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

821

822 call cpu_clock_begin(id_clock_differential_diff)

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

824 call cpu_clock_end(id_clock_differential_diff)

825

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

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

828

829 ! increment heat and salt diffusivity.

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

831 if (.not. cs%useKPP) then

832 !$OMP parallel do default(shared)

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

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

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

836 enddo ; enddo ; enddo

837 endif

838

839 endif

840

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

842 if (cs%use_CVMix_conv) then

843 ! Increment vertical diffusion and viscosity due to convection

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

845 endif

846

847 ! Find the vertical distances across layers.

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

849

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

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

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

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

854 enddo ; enddo

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

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

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

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

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

860 enddo ; enddo ; enddo

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

862

863 if (cs%debug) then

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

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

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

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

868 endif

869

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

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

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

873 enddo ; enddo ; enddo

874 if (cs%boundary_forcing_tendency_diag) then

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

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

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

878 enddo ; enddo ; enddo

879 endif

880

881 ! Apply forcing

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

883 call cpu_clock_begin(id_clock_remap)

884

885 if (cs%use_energetic_PBL) then

886

887 skinbuoyflux(:,:) = 0.0

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

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

890 cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux, mld_h=visc%h_ML)

891

892 if (cs%debug) then

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

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

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

896 unscale=us%RZ3_T3_to_W_m2*us%T_to_s)

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

898 unscale=us%kg_m3_to_R*us%degC_to_C)

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

900 unscale=us%kg_m3_to_R*us%ppt_to_S)

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

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

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

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

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

906 endif

907

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

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

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

911

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

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

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

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

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

917

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

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

920

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

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

923 if (cs%ePBL_is_additive) then

924 kd_add_here = kd_epbl(i,j,k)

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

926 else

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

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

929 endif

930

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

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

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

934

935 ! for diagnostics

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

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

938

939 enddo ; enddo ; enddo

940

941 if (cs%debug) then

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

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

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

945 endif

946

947 else

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

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

950 cs%evap_CFL_limit, cs%minimum_forcing_depth, mld_h=visc%h_ML)

951

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

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

954

955 endif ! endif for CS%use_energetic_PBL

956

957 ! diagnose the tendencies due to boundary forcing

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

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

960 if (cs%boundary_forcing_tendency_diag) then

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

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

963 endif

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

965 call diag_update_remap_grids(cs%diag)

966 call cpu_clock_end(id_clock_remap)

967 if (cs%debug) then

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

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

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

971 endif

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

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

974

975 if (cs%debug) then

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

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

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

979 endif

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

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

982

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

984 if (associated(tv%T)) then

985

986 if (cs%debug) then

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

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

989 endif

990

991 call cpu_clock_begin(id_clock_tridiag)

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

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

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

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

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

997 ! BOUND_SALINITY=False in MOM.F90.

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

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

1000

1001 if (cs%diabatic_diff_tendency_diag) then

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

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

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

1005 enddo ; enddo ; enddo

1006 endif

1007

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

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

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

1011 enddo ; enddo ; enddo

1012 if (cs%tracer_tridiag) then

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

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

1015 else

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

1017 endif

1018

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

1020 if (cs%diabatic_diff_tendency_diag) then

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

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

1023 endif

1024

1025 call cpu_clock_end(id_clock_tridiag)

1026

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

1028

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

1030

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

1032

1033 if (cs%debug) then

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

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

1036 endif

1037

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

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

1040 call diag_update_remap_grids(cs%diag)

1041

1042 ! Set diffusivities for VBF diagnostics if enabled

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

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

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

1046

1047

1048 ! Diagnose the diapycnal diffusivities and other related quantities.

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

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

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

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

1053

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

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

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

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

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

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

1060 idt = 1.0 / dt

1061 if (cs%id_Tdif > 0) then

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

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

1064 enddo ; enddo

1065 !$OMP parallel do default(shared)

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

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

1068 enddo ; enddo ; enddo

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

1070 endif

1071 if (cs%id_Sdif > 0) then

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

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

1074 enddo ; enddo

1075 !$OMP parallel do default(shared)

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

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

1078 enddo ; enddo ; enddo

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

1080 endif

1081

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

1083 n2_salt(:,:,:) = 0.0

1084 n2_temp(:,:,:) = 0.0

1085 !Compute N2 and don't mask negatives here

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

1087 if (nonbous) then

1088 !$OMP parallel do default(shared)

1089 do j=js,je

1090 if (associated(tv%p_surf)) then

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

1092 else

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

1094 endif

1095 do k=2,nz

1096 do i=is,ie

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

1098 enddo

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

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

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

1102 do i=is,ie

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

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

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

1106 enddo

1107 enddo

1108 enddo

1109 else

1110 !$OMP parallel do default(shared)

1111 do j=js,je

1112 if (associated(tv%p_surf)) then

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

1114 else

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

1116 endif

1117 do k=2,nz

1118 do i=is,ie

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

1120 enddo

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

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

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

1124 do i=is,ie

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

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

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

1128 enddo

1129 enddo

1130 enddo

1131 endif

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

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

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

1135

1136 if (cs%Use_KdWork_diag) then

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

1138 endif

1139

1140 call deallocate_vbf_cs(cs%VBF)

1141

1142 endif

1143

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

1145 call cpu_clock_begin(id_clock_tracers)

1146

1147 if (cs%mix_boundary_tracer_ALE) then

1148 tr_ea_bbl = sqrt(dt * cs%Kd_BBL_tr)

1149

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

1151 do j=js,je

1152 do i=is,ie

1153 htot(i) = 0.0

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

1155 enddo

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

1157 if (in_boundary(i)) then

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

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

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

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

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

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

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

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

1166 ! perhaps a molecular diffusivity.

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

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

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

1170 if (htot(i) < tr_ea_bbl) then

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

1172 elseif (add_ent < 0.0) then

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

1174 endif

1175

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

1177 endif

1178

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

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

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

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

1183 endif ; endif

1184 enddo ; enddo

1185

1186 enddo

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

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

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

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

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

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

1193 else

1194 add_ent = 0.0

1195 endif

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

1197 enddo ; enddo ; enddo

1198 endif ! (CS%mix_boundary_tracers)

1199

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

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

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

1203 kpp_csp=cs%KPP_CSp, &

1204 nonlocaltrans=kpp_nltscalar, &

1205 evap_cfl_limit=cs%evap_CFL_limit, &

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

1207

1208 call cpu_clock_end(id_clock_tracers)

1209

1210 ! Apply ALE sponge

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

1212 call cpu_clock_begin(id_clock_sponge)

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

1214 call cpu_clock_end(id_clock_sponge)

1215 if (cs%debug) then

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

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

1218 endif

1219 endif ! CS%use_sponge

1220

1221 ! Apply data assimilation incremental update -oda_incupd-

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

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

1224 call cpu_clock_begin(id_clock_oda_incupd)

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

1226 call cpu_clock_end(id_clock_oda_incupd)

1227 if (cs%debug) then

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

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

1230 endif

1231 endif ! CS%use_oda_incupd

1232

1233

1234

1235 call disable_averaging(cs%diag)

1236

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

1238

1239end subroutine diabatic_ale_legacy

1240

1241

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

1243!! accompanying diapycnal advection of momentum and tracers.

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

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

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

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

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

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

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

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

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

1253 !! unused have NULL ptrs

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

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

1256 !! unused fields have NULL ptrs

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

1258 !! BBL properties and related fields

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

1260 !! equations, to enable the later derived

1261 !! diagnostics, like energy budgets

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

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

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

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

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

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

1268

1269 ! local variables

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

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

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

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

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

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

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

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

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

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

1280

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1301

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

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

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

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

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

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

1308

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

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

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

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

1313

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

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

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

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

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

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

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

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

1322

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1340

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

1342

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

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

1345

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

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

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

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

1350 h_neglect = gv%H_subroundoff

1351

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

1353 g_rho0 = gv%g_Earth_Z_T2 / gv%H_to_RZ

1354 h_to_pres = gv%H_to_RZ * gv%g_Earth

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

1356

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

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

1359

1360 showcalltree = calltree_showquery()

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

1362

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

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

1365

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

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

1368

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

1370

1371 if (cs%use_geothermal) then

1372 call cpu_clock_begin(id_clock_geothermal)

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

1374 call cpu_clock_end(id_clock_geothermal)

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

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

1377 endif

1378

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

1380 ! for vertical remapping may need to be regenerated.

1381 call diag_update_remap_grids(cs%diag)

1382

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

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

1385 ! biological properties and layer thicknesses.

1386 if (associated(cs%optics)) &

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

1388

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

1390

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

1392 if (cs%use_geothermal) then

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

1394 else

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

1396 endif

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

1398 endif

1399

1400 call cpu_clock_begin(id_clock_set_diffusivity)

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

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

1403 if (cs%debug) &

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

1405 if (cs%double_diffuse) then

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

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

1408 else

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

1410 cs%set_diff_CSp, cs%VBF)

1411 endif

1412 call cpu_clock_end(id_clock_set_diffusivity)

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

1414

1415 if (cs%debug) then

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

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

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

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

1420 endif

1421

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

1423 if (cs%double_diffuse) then

1424 ! Add contributions from double diffusion

1425 !$OMP parallel do default(shared)

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

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

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

1429 enddo ; enddo ; enddo

1430 else

1431 !$OMP parallel do default(shared)

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

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

1434 enddo ; enddo ; enddo

1435 endif

1436

1437 if (cs%debug) then

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

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

1440 endif

1441

1442 if (cs%useKPP) then

1443 call cpu_clock_begin(id_clock_kpp)

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

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

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

1447 enddo ; enddo ; enddo

1448

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

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

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

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

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

1454 kpp_temp_flux(:,:) = 0.

1455 kpp_salt_flux(:,:) = 0.

1456

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

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

1459 ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux

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

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

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

1463 kpp_buoy_flux, kpp_temp_flux, kpp_salt_flux)

1464

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

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

1467

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

1469 if ( associated(fluxes%lamult) ) then

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

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

1472

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

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

1475 else

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

1477 u_star, kpp_buoy_flux, waves=waves)

1478

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

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

1481 endif

1482

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

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

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

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

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

1488

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

1490 if (cs%debug) then

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

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

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

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

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

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

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

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

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

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

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

1502 endif

1503 ! Apply non-local transport of heat and salt

1504 ! Changes: tv%T, tv%S

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

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

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

1508 dt, tv%tr_S, tv%S)

1509 call cpu_clock_end(id_clock_kpp)

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

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

1512

1513 if (cs%debug) then

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

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

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

1517 endif

1518 endif ! endif for KPP

1519

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

1521 if (cs%use_CVMix_conv) then

1522 ! Increment vertical diffusion and viscosity due to convection

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

1524 endif

1525

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

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

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

1529 enddo ; enddo ; enddo

1530 if (cs%boundary_forcing_tendency_diag) then

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

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

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

1534 enddo ; enddo ; enddo

1535 endif

1536

1537 ! Apply forcing

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

1539 call cpu_clock_begin(id_clock_remap)

1540

1541 if (cs%use_energetic_PBL) then

1542

1543 skinbuoyflux(:,:) = 0.0

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

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

1546 cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux, mld_h=visc%h_ML)

1547

1548 if (cs%debug) then

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

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

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

1552 unscale=us%RZ3_T3_to_W_m2*us%T_to_s)

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

1554 unscale=us%kg_m3_to_R*us%degC_to_C)

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

1556 unscale=us%kg_m3_to_R*us%ppt_to_S)

1557 endif

1558

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

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

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

1562

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

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

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

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

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

1568

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

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

1571 if (cs%ePBL_is_additive) then

1572 kd_add_here = kd_epbl(i,j,k)

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

1574 else

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

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

1577 endif

1578

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

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

1581 enddo ; enddo ; enddo

1582

1583 if (cs%debug) then

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

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

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

1587 endif

1588

1589 else

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

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

1592 cs%evap_CFL_limit, cs%minimum_forcing_depth, mld_h=visc%h_ML)

1593

1594 endif ! endif for CS%use_energetic_PBL

1595

1596 ! diagnose the tendencies due to boundary forcing

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

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

1599 if (cs%boundary_forcing_tendency_diag) then

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

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

1602 endif

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

1604 call diag_update_remap_grids(cs%diag)

1605 call cpu_clock_end(id_clock_remap)

1606 if (cs%debug) then

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

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

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

1610 endif

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

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

1613

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

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

1616

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

1618 if (associated(tv%T)) then

1619

1620 if (cs%debug) then

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

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

1623 endif

1624

1625 call cpu_clock_begin(id_clock_tridiag)

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

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

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

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

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

1631 ! BOUND_SALINITY=False in MOM.F90.

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

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

1634

1635 if (cs%diabatic_diff_tendency_diag) then

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

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

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

1639 enddo ; enddo ; enddo

1640 endif

1641

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

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

1644

1645 ! set ent_t=dt*Kd_heat/h_int and est_s=dt*Kd_salt/h_int on interfaces for use in the tridiagonal solver.

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

1647 ent_t(i,j,1) = 0. ; ent_t(i,j,nz+1) = 0.

1648 ent_s(i,j,1) = 0. ; ent_s(i,j,nz+1) = 0.

1649 enddo ; enddo

1650

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

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

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

1654 ent_t(i,j,k) = dt * i_dzval * kd_heat(i,j,k)

1655 ent_s(i,j,k) = dt * i_dzval * kd_salt(i,j,k)

1656 enddo ; enddo ; enddo

1657 if (showcalltree) call calltree_waypoint("done setting ent_t and ent_t from Kd_heat and " //&

1658 "Kd_salt (diabatic_ALE)")

1659

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

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

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

1663

1664 ! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below

1665 if (cs%diabatic_diff_tendency_diag) then

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

1667 endif

1668 call cpu_clock_end(id_clock_tridiag)

1669

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

1671

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

1673

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

1675

1676 if (cs%debug) then

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

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

1679 endif

1680

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

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

1683 call diag_update_remap_grids(cs%diag)

1684

1685 ! Set diffusivities for VBF diagnostics if enabled

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

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

1688 if (associated(cs%VBF%Kd_temp)) cs%VBF%Kd_salt(:,:,:) = kd_salt(:,:,:)

1689

1690 ! Diagnose the diapycnal diffusivities and other related quantities.

1691 if (cs%id_Kd_heat > 0) call post_data(cs%id_Kd_heat, kd_heat, cs%diag)

1692 if (cs%id_Kd_salt > 0) call post_data(cs%id_Kd_salt, kd_salt, cs%diag)

1693 if (cs%id_Kd_ePBL > 0) call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)

1694 if (cs%id_Kd_int > 0) then

1695 if (cs%double_diffuse .or. cs%useKPP) then

1696 ! Using this as a work array might cause confusion.

1697 do k=1,nz ; do j=js,je ; do i=is,ie

1698 kd_heat(i,j,k) = min(kd_heat(i,j,k), kd_salt(i,j,k))

1699 enddo ; enddo ; enddo

1700 endif

1701 call post_data(cs%id_Kd_int, kd_heat, cs%diag)

1702 endif

1703

1704 if (cs%id_ea_t > 0) call post_data(cs%id_ea_t, ent_t(:,:,1:nz), cs%diag)

1705 if (cs%id_eb_t > 0) call post_data(cs%id_eb_t, ent_t(:,:,2:nz+1), cs%diag)

1706 if (cs%id_ea_s > 0) call post_data(cs%id_ea_s, ent_s(:,:,1:nz), cs%diag)

1707 if (cs%id_eb_s > 0) call post_data(cs%id_eb_s, ent_s(:,:,2:nz+1), cs%diag)

1708

1709 idt = 1.0 / dt

1710 if (cs%id_Tdif > 0) then

1711 do j=js,je ; do i=is,ie

1712 tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0

1713 enddo ; enddo

1714 !$OMP parallel do default(shared)

1715 do k=2,nz ; do j=js,je ; do i=is,ie

1716 tdif_flx(i,j,k) = (idt * ent_t(i,j,k)) * (tv%T(i,j,k-1) - tv%T(i,j,k))

1717 enddo ; enddo ; enddo

1718 if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)

1719 endif

1720 if (cs%id_Sdif > 0) then

1721 do j=js,je ; do i=is,ie

1722 sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0

1723 enddo ; enddo

1724 !$OMP parallel do default(shared)

1725 do k=2,nz ; do j=js,je ; do i=is,ie

1726 sdif_flx(i,j,k) = (idt * ent_s(i,j,k)) * (tv%S(i,j,k-1) - tv%S(i,j,k))

1727 enddo ; enddo ; enddo

1728 if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)

1729 endif

1730

1731 if (cs%Use_KdWork_diag .or. cs%Use_N2_diag) then

1732 n2_salt(:,:,:) = 0.0

1733 n2_temp(:,:,:) = 0.0

1734 !Compute N2 and don't mask negatives here

1735 eosdom(:) = eos_domain(g%HI)

1736 if (nonbous) then

1737 !$OMP parallel do default(shared)

1738 do j=js,je

1739 if (associated(tv%p_surf)) then

1740 do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo

1741 else

1742 do i=is,ie ; p_i(i) = 0.0 ; enddo

1743 endif

1744 do k=2,nz

1745 do i=is,ie

1746 p_i(i) = p_i(i) + h_to_pres * h(i,j,k-1)

1747 enddo

1748 t_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))

1749 s_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))

1750 call calculate_specific_vol_derivs(t_i, s_i, p_i, dspv_dt, dspv_ds, tv%eqn_of_state, eosdom)

1751 do i=is,ie

1752 i_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))

1753 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))

1754 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))

1755 enddo

1756 enddo

1757 enddo

1758 else

1759 !$OMP parallel do default(shared)

1760 do j=js,je

1761 if (associated(tv%p_surf)) then

1762 do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo

1763 else

1764 do i=is,ie ; p_i(i) = 0.0 ; enddo

1765 endif

1766 do k=2,nz

1767 do i=is,ie

1768 p_i(i) = p_i(i) + h_to_pres* h(i,j,k-1)

1769 enddo

1770 t_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))

1771 s_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))

1772 call calculate_density_derivs(t_i, s_i, p_i, drhodt, drhods, tv%eqn_of_state, eosdom)

1773 do i=is,ie

1774 i_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))

1775 n2_salt(i,j,k) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (drhods(i) * (g_rho0 * i_h))

1776 n2_temp(i,j,k) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (drhodt(i) * (g_rho0 * i_h))

1777 enddo

1778 enddo

1779 enddo

1780 endif

1781 if (cs%id_N2_dd>0) call post_data(cs%id_N2_dd, n2_salt(:,:,:)+n2_temp(:,:,:), cs%diag)

1782 if (cs%id_N2_salt_dd>0) call post_data(cs%id_N2_salt_dd, n2_salt, cs%diag)

1783 if (cs%id_N2_temp_dd>0) call post_data(cs%id_N2_temp_dd, n2_temp, cs%diag)

1784

1785 if (cs%Use_KdWork_diag) then

1786 call kdwork_diagnostics(g,gv,us,cs%diag,cs%VBF,n2_salt,n2_temp,dz)

1787 endif

1788

1789 call deallocate_vbf_cs(cs%VBF)

1790

1791 endif

1792

1793 ! mixing of passive tracers from massless boundary layers to interior

1794 call cpu_clock_begin(id_clock_tracers)

1795

1796 if (cs%mix_boundary_tracer_ALE) then

1797 tr_ea_bbl = sqrt(dt * cs%Kd_BBL_tr)

1798 !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)

1799 do j=js,je

1800 do i=is,ie

1801 htot(i) = 0.0

1802 in_boundary(i) = (g%mask2dT(i,j) > 0.0)

1803 enddo

1804 do k=nz,2,-1 ; do i=is,ie

1805 if (in_boundary(i)) then

1806 htot(i) = htot(i) + h(i,j,k)

1807 ! If diapycnal mixing has been suppressed because this is a massless layer near the

1808 ! bottom, add some mixing of tracers between these layers. This flux is based on the

1809 ! harmonic mean of the two thicknesses, following what is done in layered mode. Kd_min_tr

1810 ! should be much less than the values in Kd_salt, perhaps a molecular diffusivity.

1811 add_ent = (dt * cs%Kd_min_tr) * &

1812 ((dz(i,j,k-1)+dz(i,j,k) + dz_neglect) / (dz(i,j,k-1)*dz(i,j,k) + dz_neglect2)) - &

1813 ent_s(i,j,k)

1814 if (htot(i) < tr_ea_bbl) then

1815 add_ent = max(0.0, add_ent, (tr_ea_bbl - htot(i)) - ent_s(i,j,k))

1816 elseif (add_ent < 0.0) then

1817 add_ent = 0.0 ; in_boundary(i) = .false.

1818 endif

1819

1820 ent_s(i,j,k) = ent_s(i,j,k) + add_ent

1821 endif

1822 enddo ; enddo

1823 enddo

1824 endif ! (CS%mix_boundary_tracer_ALE)

1825

1826 ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied

1827 call call_tracer_column_fns(h_orig, h, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), fluxes, bld, dt, &

1828 g, gv, us, tv, cs%optics, cs%tracer_flow_CSp, cs%debug, &

1829 kpp_csp=cs%KPP_CSp, &

1830 nonlocaltrans=kpp_nltscalar, &

1831 evap_cfl_limit=cs%evap_CFL_limit, &

1832 minimum_forcing_depth=cs%minimum_forcing_depth, h_bl=visc%h_ML)

1833

1834 call cpu_clock_end(id_clock_tracers)

1835

1836 ! Apply ALE sponge

1837 if (cs%use_sponge .and. associated(cs%ALE_sponge_CSp)) then

1838 call cpu_clock_begin(id_clock_sponge)

1839 call apply_ale_sponge(h, tv, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)

1840 call cpu_clock_end(id_clock_sponge)

1841 if (cs%debug) then

1842 call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)

1843 call mom_thermovar_chksum("apply_sponge ", tv, g, us)

1844 endif

1845 endif ! CS%use_sponge

1846

1847 ! Apply data assimilation incremental update -oda_incupd-

1848 if (cs%use_oda_incupd .and. associated(cs%oda_incupd_CSp)) then

1849 call mom_mesg("Starting ODA_INCUPD ", 5)

1850 call cpu_clock_begin(id_clock_oda_incupd)

1851 call apply_oda_incupd(h, tv, u, v, dt, g, gv, us, cs%oda_incupd_CSp)

1852 call cpu_clock_end(id_clock_oda_incupd)

1853 if (cs%debug) then

1854 call mom_state_chksum("apply_oda_incupd ", u, v, h, g, gv, us, haloshift=0)

1855 call mom_thermovar_chksum("apply_oda_incupd ", tv, g, us)

1856 endif

1857 endif ! CS%use_oda_incupd

1858

1859

1860 call cpu_clock_begin(id_clock_pass)

1861 ! visc%Kv_slow is not in the group pass because it has larger vertical extent.

1862 if (associated(visc%Kv_slow)) &

1863 call pass_var(visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)

1864 call cpu_clock_end(id_clock_pass)

1865

1866 call disable_averaging(cs%diag)

1867

1868 if (showcalltree) call calltree_leave("diabatic_ALE()")

1869

1870end subroutine diabatic_ale

1871

1872!> Imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers

1873!! using the original MOM6 algorithms.

1874subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &

1875 G, GV, US, CS, Waves)

1876 type(ocean_grid_type), intent(inout) :: G !< ocean grid structure

1877 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

1878 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

1879 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]

1880 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]

1881 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2]

1882 type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields

1883 !! unused have NULL ptrs

1884 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Active mixed layer depth [Z ~> m]

1885 type(forcing), intent(inout) :: fluxes !< points to forcing fields

1886 !! unused fields have NULL ptrs

1887 type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities,

1888 !! BBL properties and related fields

1889 type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum

1890 !! equations, to enable the later derived

1891 !! diagnostics, like energy budgets

1892 type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations

1893 real, intent(in) :: dt !< time increment [T ~> s]

1894 type(time_type), intent(in) :: Time_end !< Time at the end of the interval

1895 type(diabatic_cs), pointer :: CS !< module control structure

1896 type(wave_parameters_cs), pointer :: Waves !< Surface gravity waves

1897

1898 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &

1899 ea, & ! amount of fluid entrained from the layer above within

1900 ! one time step [H ~> m or kg m-2]

1901 eb, & ! amount of fluid entrained from the layer below within

1902 ! one time step [H ~> m or kg m-2]

1903 kd_lay, & ! diapycnal diffusivity of layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1904 h_orig, & ! initial layer thicknesses [H ~> m or kg m-2]

1905 dz, & ! The vertical distance between interfaces around a layer [Z ~> m]

1906 hold, & ! layer thickness before diapycnal entrainment, and later the initial

1907 ! layer thicknesses (if a mixed layer is used) [H ~> m or kg m-2]

1908 dz_old, & ! The initial vertical distance between interfaces around a layer

1909 ! or the distance before entrainment [Z ~> m]

1910 u_h, & ! Zonal velocities at thickness points after entrainment [L T-1 ~> m s-1]

1911 v_h, & ! Meridional velocities at thickness points after entrainment [L T-1 ~> m s-1]

1912 temp_diag, & ! Diagnostic array of previous temperatures [C ~> degC]

1913 saln_diag ! Diagnostic array of previous salinity [S ~> ppt]

1914 real, dimension(SZI_(G),SZJ_(G)) :: &

1915 h_MLD, & ! Active mixed layer thickness [H ~> m or kg m-2].

1916 U_star, & ! The friction velocity [Z T-1 ~> m s-1].

1917 KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

1918 KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

1919 Rcv_ml ! Coordinate density of mixed layer [R ~> kg m-3], used for applying sponges

1920

1921 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &

1922 ! These are targets so that the space can be shared with eaml & ebml.

1923 eatr, & ! The equivalent of ea for tracers, which differs from ea in that it tends to

1924 ! homogenize tracers in massless layers near the boundaries [H ~> m or kg m-2]

1925 ebtr ! The equivalent of eb for tracers, which differs from eb in that it tends to

1926 ! homogenize tracers in massless layers near the boundaries [H ~> m or kg m-2]

1927

1928 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &

1929 Kd_int, & ! diapycnal diffusivity of interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1930 Kd_heat, & ! diapycnal diffusivity of heat [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1931 Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1932 Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to

1933 ! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1934 kd_extra_s , & ! The extra diffusivity of salinity due to double diffusion relative to

1935 ! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

1936 kpp_nltheat, & ! KPP non-local transport for heat [nondim]

1937 kpp_nltscalar, & ! KPP non-local transport for scalars [nondim]

1938 kpp_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]

1939 tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

1940 tadv_flx, & ! advective diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]

1941 sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

1942 sadv_flx ! advective diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

1943

1944 ! The following 3 variables are only used with a bulk mixed layer.

1945 real, pointer, dimension(:,:,:) :: &

1946 eaml, & ! The equivalent of ea due to mixed layer processes [H ~> m or kg m-2].

1947 ebml ! The equivalent of eb due to mixed layer processes [H ~> m or kg m-2].

1948 ! eaml and ebml are pointers to eatr and ebtr so as to reuse the memory as

1949 ! the arrays are not needed at the same time.

1950

1951 integer :: kb(SZI_(G),SZJ_(G)) ! index of the lightest layer denser

1952 ! than the buffer layer [nondim]

1953

1954 real :: p_ref_cv(SZI_(G)) ! Reference pressure for the potential density that defines the

1955 ! coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].

1956

1957 logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below,

1958 ! where massive is defined as sufficiently thick that

1959 ! the no-flux boundary conditions have not restricted

1960 ! the entrainment - usually sqrt(Kd*dt).

1961

1962 real :: h_neglect ! A thickness that is so small it is usually lost

1963 ! in roundoff and can be neglected [H ~> m or kg m-2]

1964 real :: dz_neglect ! A vertical distance that is so small it is usually lost

1965 ! in roundoff and can be neglected [Z ~> m]

1966 real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2]

1967 real :: net_ent ! The net of ea-eb at an interface [H ~> m or kg m-2]

1968 real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2]

1969 real :: eaval ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]

1970 real :: hval ! hval is 2*h at velocity grid points [H ~> m or kg m-2]

1971 real :: h_tr ! h_tr is h at tracer points with a tiny thickness

1972 ! added to ensure positive definiteness [H ~> m or kg m-2]

1973 real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is

1974 ! coupled to the bottom within a timestep [H ~> m or kg m-2]

1975

1976 real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].

1977 real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]

1978 real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2]

1979 real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim]

1980 real :: c1(SZIB_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim]

1981

1982 real :: dt_mix ! The amount of time over which to apply mixing [T ~> s]

1983 real :: Idt ! The inverse time step [T-1 ~> s-1]

1984

1985 integer :: dir_flag ! An integer encoding the directions in which to do halo updates.

1986 logical :: showCallTree ! If true, show the call tree

1987 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state

1988 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, halo

1989

1990 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

1991 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

1992 nkmb = gv%nk_rho_varies

1993 h_neglect = gv%H_subroundoff

1994 dz_neglect = gv%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect

1995 kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0

1996

1997 showcalltree = calltree_showquery()

1998 if (showcalltree) call calltree_enter("layered_diabatic(), MOM_diabatic_driver.F90")

1999

2000 ! set equivalence between the same bits of memory for these arrays

2001 eaml => eatr ; ebml => ebtr

2002

2003 ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep

2004 call enable_averages(dt, time_end, cs%diag)

2005

2006 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then

2007 halo = cs%halo_TS_diff

2008 !$OMP parallel do default(shared)

2009 do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo

2010 h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0

2011 enddo ; enddo ; enddo

2012 endif

2013

2014 if (cs%use_geothermal) then

2015 call cpu_clock_begin(id_clock_geothermal)

2016 call geothermal_entraining(h, tv, dt, eaml, ebml, g, gv, us, cs%geothermal, halo=cs%halo_TS_diff)

2017 call cpu_clock_end(id_clock_geothermal)

2018 if (showcalltree) call calltree_waypoint("geothermal (diabatic)")

2019 if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)

2020 endif

2021

2022 ! Whenever thickness changes let the diag manager know, target grids

2023 ! for vertical remapping may need to be regenerated.

2024 call diag_update_remap_grids(cs%diag)

2025

2026 ! Set_pen_shortwave estimates the optical properties of the water column.

2027 ! It will need to be modified later to include information about the

2028 ! biological properties and layer thicknesses.

2029 if (associated(cs%optics)) &

2030 call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity, cs%tracer_flow_CSp)

2031

2032 if (cs%use_bulkmixedlayer) then

2033 if (cs%debug) call mom_forcing_chksum("Before mixedlayer", fluxes, g, us, haloshift=0)

2034

2035 if (cs%ML_mix_first > 0.0) then

2036! This subroutine

2037! (1) Cools the mixed layer.

2038! (2) Performs convective adjustment by mixed layer entrainment.

2039! (3) Heats the mixed layer and causes it to detrain to

2040! Monin-Obukhov depth or minimum mixed layer depth.

2041! (4) Uses any remaining TKE to drive mixed layer entrainment.

2042! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate)

2043 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)

2044

2045 call cpu_clock_begin(id_clock_mixedlayer)

2046 if (cs%ML_mix_first < 1.0) then

2047 ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???)

2048 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*cs%ML_mix_first, &

2049 eaml, ebml, g, gv, us, cs%bulkmixedlayer, cs%optics, &

2050 bld, h_mld, cs%aggregate_FW_forcing, dt, last_call=.false.)

2051 else

2052 ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???)

2053 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, &

2054 g, gv, us, cs%bulkmixedlayer, cs%optics, &

2055 bld, h_mld, cs%aggregate_FW_forcing, dt, last_call=.true.)

2056 if (associated(visc%h_ML)) visc%h_ML(:,:) = h_mld(:,:)

2057 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

2058 endif

2059

2060 ! Keep salinity from falling below a small but positive threshold.

2061 ! This constraint is needed for SIS1 ice model, which can extract

2062 ! more salt than is present in the ocean. SIS2 does not suffer

2063 ! from this limitation, in which case we can let salinity=0 and still

2064 ! have salt conserved with SIS2 ice. So for SIS2, we can run with

2065 ! BOUND_SALINITY=False in MOM.F90.

2066 if (associated(tv%S) .and. associated(tv%salt_deficit)) &

2067 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

2068 call cpu_clock_end(id_clock_mixedlayer)

2069 if (cs%debug) then

2070 call mom_state_chksum("After mixedlayer ", u, v, h, g, gv, us, haloshift=0)

2071 call mom_forcing_chksum("After mixedlayer", fluxes, g, us, haloshift=0)

2072 endif

2073 if (showcalltree) call calltree_waypoint("done with 1st bulkmixedlayer (diabatic)")

2074 if (cs%debugConservation) call mom_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, g, gv, us)

2075 endif

2076 endif

2077

2078 if (cs%debug) &

2079 call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)

2080 if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then

2081 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then

2082 call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, us, eaml, ebml)

2083 if (cs%debug) then

2084 call hchksum(eaml, "after find_uv_at_h eaml", g%HI, unscale=gv%H_to_MKS)

2085 call hchksum(ebml, "after find_uv_at_h ebml", g%HI, unscale=gv%H_to_MKS)

2086 endif

2087 else

2088 call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)

2089 endif

2090 if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")

2091 endif

2092

2093 call cpu_clock_begin(id_clock_set_diffusivity)

2094 ! Sets: Kd_lay, Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb

2095 ! Also changes: visc%Kd_shear and visc%Kv_shear

2096 if ((cs%halo_TS_diff > 0) .and. (cs%ML_mix_first > 0.0)) then

2097 if (associated(tv%T)) call pass_var(tv%T, g%Domain, halo=cs%halo_TS_diff, complete=.false.)

2098 if (associated(tv%S)) call pass_var(tv%S, g%Domain, halo=cs%halo_TS_diff, complete=.false.)

2099 call pass_var(h, g%domain, halo=cs%halo_TS_diff, complete=.true.)

2100 endif

2101

2102 ! Update derived thermodynamic quantities.

2103 if ((cs%ML_mix_first > 0.0) .and. allocated(tv%SpV_avg)) then

2104 call calc_derived_thermo(tv, h, g, gv, us, halo=cs%halo_TS_diff)

2105 endif

2106

2107 if (cs%debug) &

2108 call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)

2109 if (cs%double_diffuse) then

2110 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, kd_int, g, gv, us, &

2111 cs%set_diff_CSp, cs%VBF, kd_lay=kd_lay, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)

2112 else

2113 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, kd_int, g, gv, us, &

2114 cs%set_diff_CSp, cs%VBF, kd_lay=kd_lay)

2115 endif

2116 call cpu_clock_end(id_clock_set_diffusivity)

2117 if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")

2118

2119 if (cs%debug) then

2120 call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)

2121 call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)

2122 call mom_thermovar_chksum("after set_diffusivity ", tv, g, us)

2123 call hchksum(kd_lay, "after set_diffusivity Kd_lay", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

2124 call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

2125 endif

2126

2127

2128 if (cs%useKPP) then

2129 call cpu_clock_begin(id_clock_kpp)

2130

2131 ! NOTE: The following do not require initialization, but their checksums do

2132 ! require initialization, and past versions were initialized to zero.

2133 kpp_nltheat(:,:,:) = 0.

2134 kpp_nltscalar(:,:,:) = 0.

2135 kpp_buoy_flux(:,:,:) = 0.

2136 kpp_temp_flux(:,:) = 0.

2137 kpp_salt_flux(:,:) = 0.

2138

2139 ! KPP needs the surface buoyancy flux but does not update state variables.

2140 ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.

2141 ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux

2142 ! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second)

2143 ! unlike other instances where the fluxes are integrated in time over a time-step.

2144 call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &

2145 kpp_buoy_flux, kpp_temp_flux, kpp_salt_flux)

2146 ! The KPP scheme calculates boundary layer diffusivities and non-local transport.

2147

2148 ! Set diffusivities for heat and salt separately

2149

2150 if (cs%double_diffuse) then

2151 ! Add contribution from double diffusion

2152 !$OMP parallel do default(shared)

2153 do k=1,nz+1 ; do j=js,je ; do i=is,ie

2154 kd_salt(i,j,k) = kd_int(i,j,k) + kd_extra_s(i,j,k)

2155 kd_heat(i,j,k) = kd_int(i,j,k) + kd_extra_t(i,j,k)

2156 enddo ; enddo ; enddo

2157 else

2158 !$OMP parallel do default(shared)

2159 do k=1,nz+1 ; do j=js,je ; do i=is,ie

2160 kd_salt(i,j,k) = kd_int(i,j,k)

2161 kd_heat(i,j,k) = kd_int(i,j,k)

2162 enddo ; enddo ; enddo

2163 endif

2164

2165 ! Determine the friction velocity, perhaps using the evovling surface density.

2166 call find_ustar(fluxes, tv, u_star, g, gv, us)

2167

2168 if ( associated(fluxes%lamult) ) then

2169 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &

2170 u_star, kpp_buoy_flux, waves=waves, lamult=fluxes%lamult)

2171

2172 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, tv, u_star, kpp_buoy_flux, kd_heat, &

2173 kd_salt, visc%Kv_shear, kpp_nltheat, kpp_nltscalar, waves=waves, lamult=fluxes%lamult)

2174 else

2175 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &

2176 u_star, kpp_buoy_flux, waves=waves)

2177

2178 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, tv, u_star, kpp_buoy_flux, kd_heat, &

2179 kd_salt, visc%Kv_shear, kpp_nltheat, kpp_nltscalar, waves=waves)

2180 endif

2181

2182 call kpp_get_bld(cs%KPP_CSp, bld(:,:), g, us)

2183 ! If visc%MLD or visc%h_ML exist, copy KPP's BLD into them with appropriate conversions.

2184 if (associated(visc%h_ML)) call convert_mld_to_ml_thickness(bld, h, visc%h_ML, tv, g, gv)

2185 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

2186 if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = kpp_buoy_flux(:,:,1)

2187

2188 if (.not. cs%KPPisPassive) then

2189 !$OMP parallel do default(shared)

2190 do k=1,nz+1 ; do j=js,je ; do i=is,ie

2191 kd_int(i,j,k) = min( kd_salt(i,j,k), kd_heat(i,j,k) )

2192 enddo ; enddo ; enddo

2193 if (cs%double_diffuse) then

2194 !$OMP parallel do default(shared)

2195 do k=1,nz+1 ; do j=js,je ; do i=is,ie

2196 kd_extra_s(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))

2197 kd_extra_t(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))

2198 enddo ; enddo ; enddo

2199 endif

2200 endif ! not passive

2201

2202 call cpu_clock_end(id_clock_kpp)

2203 if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")

2204 if (cs%debug) then

2205 call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)

2206 call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)

2207 call mom_thermovar_chksum("after KPP", tv, g, us)

2208 call hchksum(kd_lay, "after KPP Kd_lay", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

2209 call hchksum(kd_int, "after KPP Kd_Int", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)

2210 endif

2211 endif ! endif for KPP

2212

2213 ! Add vertical diff./visc. due to convection (computed via CVMix)

2214 if (cs%use_CVMix_conv) then

2215 call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv, bld, kd_int, visc%Kv_shear)

2216 endif

2217

2218 if (cs%useKPP) then

2219 call cpu_clock_begin(id_clock_kpp)

2220 if (cs%debug) then

2221 call hchksum(kpp_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, &

2222 unscale=us%C_to_degC*gv%H_to_m*us%s_to_T)

2223 call hchksum(kpp_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, &

2224 unscale=us%S_to_ppt*gv%H_to_m*us%s_to_T)

2225 call hchksum(kpp_nltheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)

2226 call hchksum(kpp_nltscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)

2227 endif

2228 ! Apply non-local transport of heat and salt

2229 ! Changes: tv%T, tv%S

2230 call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, kpp_nltheat, kpp_temp_flux, &

2231 dt, tv%tr_T, tv%T, tv%C_p)

2232 call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, kpp_nltscalar, kpp_salt_flux, &

2233 dt, tv%tr_S, tv%S)

2234 call cpu_clock_end(id_clock_kpp)

2235 if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")

2236 if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)

2237

2238 if (cs%debug) then

2239 call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)

2240 call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)

2241 call mom_thermovar_chksum("after KPP_applyNLT ", tv, g, us)

2242 endif

2243 endif ! endif for KPP

2244

2245 ! Differential diffusion done here.

2246 ! Changes: tv%T, tv%S

2247 if (cs%double_diffuse .and. associated(tv%T)) then

2248

2249 call cpu_clock_begin(id_clock_differential_diff)

2250 call differential_diffuse_t_s(h, tv%T, tv%S, kd_extra_t, kd_extra_s, tv, dt, g, gv)

2251 call cpu_clock_end(id_clock_differential_diff)

2252 if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")

2253 if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g, gv, us)

2254

2255 ! increment heat and salt diffusivity.

2256 ! CS%useKPP==.true. already has extra_T and extra_S included

2257 if (.not. cs%useKPP) then

2258 !$OMP parallel do default(shared)

2259 do k=2,nz ; do j=js,je ; do i=is,ie

2260 kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)

2261 kd_salt(i,j,k) = kd_salt(i,j,k) + kd_extra_s(i,j,k)

2262 enddo ; enddo ; enddo

2263 endif

2264

2265 endif

2266

2267 ! Calculate layer entrainments and detrainments from diffusivities and differences between

2268 ! layer and target densities (i.e. do remapping as well as diffusion).

2269 call cpu_clock_begin(id_clock_entrain)

2270 ! Calculate appropriately limited diapycnal mass fluxes to account

2271 ! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb

2272 call entrainment_diffusive(h, tv, fluxes, dt, g, gv, us, cs%entrain_diffusive, &

2273 ea, eb, kb, kd_lay=kd_lay, kd_int=kd_int)

2274 call cpu_clock_end(id_clock_entrain)

2275 if (showcalltree) call calltree_waypoint("done with Entrainment_diffusive (diabatic)")

2276

2277 if (cs%debug) then

2278 call mom_forcing_chksum("after calc_entrain ", fluxes, g, us, haloshift=0)

2279 call mom_thermovar_chksum("after calc_entrain ", tv, g, us)

2280 call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, us, haloshift=0)

2281 call hchksum(ea, "after calc_entrain ea", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2282 call hchksum(eb, "after calc_entrain eb", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2283 endif

2284

2285 ! Save fields before boundary forcing is applied for tendency diagnostics

2286 if (cs%boundary_forcing_tendency_diag) then

2287 do k=1,nz ; do j=js,je ; do i=is,ie

2288 temp_diag(i,j,k) = tv%T(i,j,k)

2289 saln_diag(i,j,k) = tv%S(i,j,k)

2290 enddo ; enddo ; enddo

2291 endif

2292

2293 ! Update h according to divergence of the difference between

2294 ! ea and eb. We keep a record of the original h in hold.

2295 ! In the following, the checks for negative values are to guard

2296 ! against instances where entrainment drives a layer to

2297 ! negative thickness. This situation will never happen if

2298 ! enough iterations are permitted in Calculate_Entrainment.

2299 ! Even if too few iterations are allowed, it is still guarded

2300 ! against. In other words the checks are probably unnecessary.

2301 !$OMP parallel do default(shared)

2302 do j=js,je

2303 do i=is,ie

2304 hold(i,j,1) = h(i,j,1)

2305 h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2))

2306 hold(i,j,nz) = h(i,j,nz)

2307 h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1))

2308 if (h(i,j,1) <= 0.0) then

2309 h(i,j,1) = gv%Angstrom_H

2310 endif

2311 if (h(i,j,nz) <= 0.0) then

2312 h(i,j,nz) = gv%Angstrom_H

2313 endif

2314 enddo

2315 do k=2,nz-1 ; do i=is,ie

2316 hold(i,j,k) = h(i,j,k)

2317 h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &

2318 (eb(i,j,k) - ea(i,j,k+1)))

2319 if (h(i,j,k) <= 0.0) then

2320 h(i,j,k) = gv%Angstrom_H

2321 endif

2322 enddo ; enddo

2323 enddo

2324 ! Checks for negative thickness may have changed layer thicknesses

2325 call diag_update_remap_grids(cs%diag)

2326

2327 if (cs%debug) then

2328 call mom_state_chksum("after negative check ", u, v, h, g, gv, us, haloshift=0)

2329 call mom_forcing_chksum("after negative check ", fluxes, g, us, haloshift=0)

2330 call mom_thermovar_chksum("after negative check ", tv, g, us)

2331 endif

2332 if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")

2333 if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)

2334

2335 ! Here, T and S are updated according to ea and eb.

2336 ! If using the bulk mixed layer, T and S are also updated

2337 ! by surface fluxes (in fluxes%*).

2338 ! This is a very long block.

2339 if (cs%use_bulkmixedlayer) then

2340

2341 if (associated(tv%T)) then

2342 call cpu_clock_begin(id_clock_tridiag)

2343 ! Temperature and salinity (as state variables) are treated

2344 ! differently from other tracers to insure massless layers that

2345 ! are lighter than the mixed layer have temperatures and salinities

2346 ! that correspond to their prescribed densities.

2347 if (cs%massless_match_targets) then

2348 !$OMP parallel do default (shared) private(h_tr,b1,d1,c1,b_denom_1)

2349 do j=js,je

2350 do i=is,ie

2351 h_tr = hold(i,j,1) + h_neglect

2352 b1(i) = 1.0 / (h_tr + eb(i,j,1))

2353 d1(i) = h_tr * b1(i)

2354 tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1))

2355 tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1))

2356 enddo

2357 do k=2,nkmb ; do i=is,ie

2358 c1(i,k) = eb(i,j,k-1) * b1(i)

2359 h_tr = hold(i,j,k) + h_neglect

2360 b_denom_1 = h_tr + d1(i)*ea(i,j,k)

2361 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))

2362 if (k<nkmb) d1(i) = b_denom_1 * b1(i)

2363 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))

2364 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))

2365 enddo ; enddo

2366

2367 do k=nkmb+1,nz ; do i=is,ie

2368 if (k == kb(i,j)) then

2369 c1(i,k) = eb(i,j,k-1) * b1(i)

2370 d1(i) = (((eb(i,j,nkmb)-eb(i,j,k-1)) + hold(i,j,nkmb) + h_neglect) + &

2371 d1(i)*ea(i,j,nkmb)) * b1(i)

2372 h_tr = hold(i,j,k) + h_neglect

2373 b_denom_1 = h_tr + d1(i)*ea(i,j,k)

2374 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))

2375 d1(i) = b_denom_1 * b1(i)

2376 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,nkmb))

2377 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,nkmb))

2378 elseif (k > kb(i,j)) then

2379 c1(i,k) = eb(i,j,k-1) * b1(i)

2380 h_tr = hold(i,j,k) + h_neglect

2381 b_denom_1 = h_tr + d1(i)*ea(i,j,k)

2382 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))

2383 d1(i) = b_denom_1 * b1(i)

2384 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))

2385 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))

2386 elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j))

2387 ! The bottommost buffer layer might entrain all the mass from some

2388 ! of the interior layers that are thin and lighter in the coordinate

2389 ! density than that buffer layer. The T and S of these newly

2390 ! massless interior layers are unchanged.

2391 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)

2392 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)

2393 endif

2394 enddo ; enddo

2395

2396 do k=nz-1,nkmb,-1 ; do i=is,ie

2397 if (k >= kb(i,j)) then

2398 tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)

2399 tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)

2400 endif

2401 enddo ; enddo

2402 do i=is,ie ; if (kb(i,j) <= nz) then

2403 tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j))

2404 tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j))

2405 endif ; enddo

2406 do k=nkmb-1,1,-1 ; do i=is,ie

2407 tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)

2408 tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)

2409 enddo ; enddo

2410 enddo ! end of j loop

2411 else ! .not. massless_match_targets

2412 ! This simpler form allows T & S to be too dense for the layers

2413 ! between the buffer layers and the interior.

2414 ! Changes: T, S

2415 if (cs%tracer_tridiag) then

2416 call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)

2417 call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)

2418 else

2419 call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)

2420 endif

2421 endif ! massless_match_targets

2422 call cpu_clock_end(id_clock_tridiag)

2423

2424 endif ! endif for associated(T)

2425 if (cs%debugConservation) call mom_state_stats('BML tridiag', u, v, h, tv%T, tv%S, g, gv, us)

2426

2427 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then

2428 ! The mixed layer code has already been called, but there is some needed

2429 ! bookkeeping.

2430 !$OMP parallel do default(shared)

2431 do k=1,nz ; do j=js,je ; do i=is,ie

2432 hold(i,j,k) = h_orig(i,j,k)

2433 ea(i,j,k) = ea(i,j,k) + eaml(i,j,k)

2434 eb(i,j,k) = eb(i,j,k) + ebml(i,j,k)

2435 enddo ; enddo ; enddo

2436 if (cs%debug) then

2437 call hchksum(ea, "after ea = ea + eaml", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2438 call hchksum(eb, "after eb = eb + ebml", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2439 endif

2440 endif

2441

2442 if (cs%ML_mix_first < 1.0) then

2443 ! Call the mixed layer code now, perhaps for a second time.

2444 ! This subroutine (1) Cools the mixed layer.

2445 ! (2) Performs convective adjustment by mixed layer entrainment.

2446 ! (3) Heats the mixed layer and causes it to detrain to

2447 ! Monin-Obukhov depth or minimum mixed layer depth.

2448 ! (4) Uses any remaining TKE to drive mixed layer entrainment.

2449 ! (5) Possibly splits the buffer layer into two isopycnal layers.

2450

2451 call find_uv_at_h(u, v, hold, u_h, v_h, g, gv, us, ea, eb)

2452 if (cs%debug) call mom_state_chksum("find_uv_at_h1 ", u, v, h, g, gv, us, haloshift=0)

2453

2454 dt_mix = min(dt, dt*(1.0 - cs%ML_mix_first))

2455 call cpu_clock_begin(id_clock_mixedlayer)

2456 ! Changes: h, tv%T, tv%S, ea and eb (G is also inout???)

2457 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, &

2458 g, gv, us, cs%bulkmixedlayer, cs%optics, &

2459 bld, h_mld, cs%aggregate_FW_forcing, dt, last_call=.true.)

2460 if (associated(visc%h_ML)) visc%h_ML(:,:) = h_mld(:,:)

2461 if (associated(visc%MLD)) visc%MLD(:,:) = bld(:,:)

2462

2463 ! Keep salinity from falling below a small but positive threshold.

2464 ! This constraint is needed for SIS1 ice model, which can extract

2465 ! more salt than is present in the ocean. SIS2 does not suffer

2466 ! from this limitation, in which case we can let salinity=0 and still

2467 ! have salt conserved with SIS2 ice. So for SIS2, we can run with

2468 ! BOUND_SALINITY=False in MOM.F90.

2469 if (associated(tv%S) .and. associated(tv%salt_deficit)) &

2470 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

2471

2472 call cpu_clock_end(id_clock_mixedlayer)

2473 if (showcalltree) call calltree_waypoint("done with 2nd bulkmixedlayer (diabatic)")

2474 if (cs%debugConservation) call mom_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, g, gv, us)

2475 endif

2476

2477 else ! following block for when NOT using BULKMIXEDLAYER

2478

2479 ! calculate change in temperature & salinity due to dia-coordinate surface diffusion

2480 if (associated(tv%T)) then

2481

2482 if (cs%debug) then

2483 call hchksum(ea, "before triDiagTS ea ", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2484 call hchksum(eb, "before triDiagTS eb ", g%HI, haloshift=0, unscale=gv%H_to_MKS)

2485 endif

2486 call cpu_clock_begin(id_clock_tridiag)

2487

2488 ! Keep salinity from falling below a small but positive threshold.

2489 ! This constraint is needed for SIS1 ice model, which can extract

2490 ! more salt than is present in the ocean. SIS2 does not suffer

2491 ! from this limitation, in which case we can let salinity=0 and still

2492 ! have salt conserved with SIS2 ice. So for SIS2, we can run with

2493 ! BOUND_SALINITY=False in MOM.F90.

2494 if (associated(tv%S) .and. associated(tv%salt_deficit)) &

2495 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

2496

2497 if (cs%diabatic_diff_tendency_diag) then

2498 do k=1,nz ; do j=js,je ; do i=is,ie

2499 temp_diag(i,j,k) = tv%T(i,j,k)

2500 saln_diag(i,j,k) = tv%S(i,j,k)

2501 enddo ; enddo ; enddo

2502 endif

2503

2504 ! Changes T and S via the tridiagonal solver; no change to h

2505 if (cs%tracer_tridiag) then

2506 call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)

2507 call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)

2508 else

2509 call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)

2510 endif

2511

2512 ! diagnose temperature, salinity, heat, and salt tendencies

2513 ! Note: hold here refers to the thicknesses from before the dual-entrainment when using

2514 ! the bulk mixed layer scheme, so tendencies should be posted on hold.

2515 if (cs%diabatic_diff_tendency_diag) then

2516 call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, us, cs)

2517 if (cs%id_diabatic_diff_h > 0) call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h=hold)

2518 endif

2519

2520 call cpu_clock_end(id_clock_tridiag)

2521 if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")

2522

2523 endif ! endif corresponding to if (associated(tv%T))

2524 if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)

2525

2526 endif ! endif for the BULKMIXEDLAYER block

2527

2528 if (cs%debug) then

2529 call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)

2530 call mom_thermovar_chksum("after mixed layer ", tv, g, us)

2531 call hchksum(ea, "after mixed layer ea", g%HI, unscale=gv%H_to_MKS)

2532 call hchksum(eb, "after mixed layer eb", g%HI, unscale=gv%H_to_MKS)

2533 endif

2534

2535 call cpu_clock_begin(id_clock_remap)

2536 call regularize_layers(h, tv, dt, ea, eb, g, gv, us, cs%regularize_layers)

2537 call cpu_clock_end(id_clock_remap)

2538 if (showcalltree) call calltree_waypoint("done with regularize_layers (diabatic)")

2539 if (cs%debugConservation) call mom_state_stats('regularize_layers', u, v, h, tv%T, tv%S, g, gv, us)

2540

2541 ! Whenever thickness changes let the diag manager know, as the

2542 ! target grids for vertical remapping may need to be regenerated.

2543 if (associated(adp%du_dt_dia) .or. associated(adp%dv_dt_dia)) &

2544 ! Remapped d[uv]dt_dia require east/north halo updates of h

2545 call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)

2546 call diag_update_remap_grids(cs%diag)

2547

2548 ! diagnostics

2549 idt = 1.0 / dt

2550 if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0)) then

2551 do j=js,je ; do i=is,ie

2552 tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0

2553 tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0

2554 enddo ; enddo

2555 !$OMP parallel do default(shared)

2556 do k=2,nz ; do j=js,je ; do i=is,ie

2557 tdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &

2558 (tv%T(i,j,k-1) - tv%T(i,j,k))

2559 tadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &

2560 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))

2561 enddo ; enddo ; enddo

2562 endif

2563 if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0)) then

2564 do j=js,je ; do i=is,ie

2565 sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0

2566 sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0

2567 enddo ; enddo

2568 !$OMP parallel do default(shared)

2569 do k=2,nz ; do j=js,je ; do i=is,ie

2570 sdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &

2571 (tv%S(i,j,k-1) - tv%S(i,j,k))

2572 sadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &

2573 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))

2574 enddo ; enddo ; enddo

2575 endif

2576

2577 ! mixing of passive tracers from massless boundary layers to interior

2578 call cpu_clock_begin(id_clock_tracers)

2579

2580 ! Find the vertical distances across layers.

2581 if (cs%mix_boundary_tracers .or. cs%double_diffuse) &

2582 call thickness_to_dz(h, tv, dz, g, gv, us)

2583 if (cs%double_diffuse) &

2584 call thickness_to_dz(hold, tv, dz_old, g, gv, us)

2585

2586 if (cs%mix_boundary_tracers) then

2587 tr_ea_bbl = sqrt(dt * cs%Kd_BBL_tr)

2588 !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)

2589 do j=js,je

2590 do i=is,ie

2591 ebtr(i,j,nz) = eb(i,j,nz)

2592 htot(i) = 0.0

2593 in_boundary(i) = (g%mask2dT(i,j) > 0.0)

2594 enddo

2595 do k=nz,2,-1 ; do i=is,ie

2596 if (in_boundary(i)) then

2597 htot(i) = htot(i) + h(i,j,k)

2598 ! If diapycnal mixing has been suppressed because this is a massless

2599 ! layer near the bottom, add some mixing of tracers between these

2600 ! layers. This flux is based on the harmonic mean of the two

2601 ! thicknesses, as this corresponds pretty closely (to within

2602 ! differences in the density jumps between layers) with what is done

2603 ! in the calculation of the fluxes in the first place. Kd_min_tr

2604 ! should be much less than the values that have been set in Kd_lay,

2605 ! perhaps a molecular diffusivity.

2606 add_ent = (dt * cs%Kd_min_tr) * &

2607 ((dz(i,j,k-1) + dz(i,j,k) + dz_neglect) / &

2608 (dz(i,j,k-1)*dz(i,j,k) + dz_neglect2)) - &

2609 0.5*(ea(i,j,k) + eb(i,j,k-1))

2610 if (htot(i) < tr_ea_bbl) then

2611 add_ent = max(0.0, add_ent, &

2612 (tr_ea_bbl - htot(i)) - min(ea(i,j,k), eb(i,j,k-1)))

2613 elseif (add_ent < 0.0) then

2614 add_ent = 0.0 ; in_boundary(i) = .false.

2615 endif

2616

2617 ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent

2618 eatr(i,j,k) = ea(i,j,k) + add_ent

2619 else

2620 ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k)

2621 endif

2622 if (cs%double_diffuse) then ; if (kd_extra_s(i,j,k) > 0.0) then

2623 add_ent = (dt * kd_extra_s(i,j,k)) / &

2624 (0.25 * ((dz(i,j,k-1) + dz(i,j,k)) + (dz_old(i,j,k-1) + dz_old(i,j,k))) + dz_neglect)

2625 ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent

2626 eatr(i,j,k) = eatr(i,j,k) + add_ent

2627 endif ; endif

2628 enddo ; enddo

2629 do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo

2630

2631 enddo

2632

2633 call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, bld, dt, g, gv, us, tv, &

2634 cs%optics, cs%tracer_flow_CSp, cs%debug, &

2635 kpp_csp=cs%KPP_CSp, nonlocaltrans=kpp_nltscalar, h_bl=visc%h_ML)

2636

2637 elseif (cs%double_diffuse) then ! extra diffusivity for passive tracers

2638

2639 do j=js,je ; do i=is,ie

2640 ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1)

2641 enddo ; enddo

2642 !$OMP parallel do default(shared) private(add_ent)

2643 do k=nz,2,-1 ; do j=js,je ; do i=is,ie

2644 if (kd_extra_s(i,j,k) > 0.0) then

2645 add_ent = (dt * kd_extra_s(i,j,k)) / &

2646 (0.25 * ((dz(i,j,k-1) + dz(i,j,k)) + (dz_old(i,j,k-1) + dz_old(i,j,k))) + dz_neglect)

2647 else

2648 add_ent = 0.0

2649 endif

2650 ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent

2651 eatr(i,j,k) = ea(i,j,k) + add_ent

2652 enddo ; enddo ; enddo

2653

2654 call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, bld, dt, g, gv, us, tv, &

2655 cs%optics, cs%tracer_flow_CSp, cs%debug, &

2656 kpp_csp=cs%KPP_CSp, nonlocaltrans=kpp_nltscalar, h_bl=visc%h_ML)

2657

2658 else

2659 call call_tracer_column_fns(hold, h, ea, eb, fluxes, bld, dt, g, gv, us, tv, &

2660 cs%optics, cs%tracer_flow_CSp, cs%debug, &

2661 kpp_csp=cs%KPP_CSp, nonlocaltrans=kpp_nltscalar, h_bl=visc%h_ML)

2662

2663 endif ! (CS%mix_boundary_tracers)

2664

2665 call cpu_clock_end(id_clock_tracers)

2666

2667 ! sponges

2668 if (cs%use_sponge) then

2669 call cpu_clock_begin(id_clock_sponge)

2670 ! Layer mode sponge

2671 if (cs%use_bulkmixedlayer .and. associated(tv%eqn_of_state)) then

2672 do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo

2673 eosdom(:) = eos_domain(g%HI)

2674 !$OMP parallel do default(shared)

2675 do j=js,je

2676 call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, rcv_ml(:,j), &

2677 tv%eqn_of_state, eosdom)

2678 enddo

2679 call apply_sponge(h, tv, dt, g, gv, us, ea, eb, cs%sponge_CSp, rcv_ml)

2680 else

2681 call apply_sponge(h, tv, dt, g, gv, us, ea, eb, cs%sponge_CSp)

2682 endif

2683 call cpu_clock_end(id_clock_sponge)

2684 if (cs%debug) then

2685 call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)

2686 call mom_thermovar_chksum("apply_sponge ", tv, g, us)

2687 endif

2688 endif ! CS%use_sponge

2689

2690 ! Apply data assimilation incremental update -oda_incupd-

2691 if (cs%use_oda_incupd .and. associated(cs%oda_incupd_CSp)) then

2692 call cpu_clock_begin(id_clock_oda_incupd)

2693 call apply_oda_incupd(h, tv, u, v, dt, g, gv, us, cs%oda_incupd_CSp)

2694 call cpu_clock_end(id_clock_oda_incupd)

2695 if (cs%debug) then

2696 call mom_state_chksum("apply_oda_incupd ", u, v, h, g, gv, us, haloshift=0)

2697 call mom_thermovar_chksum("apply_oda_incupd ", tv, g, us)

2698 endif

2699 endif ! CS%use_oda_incupd

2700

2701

2702

2703! Save the diapycnal mass fluxes as a diagnostic field.

2704 if (associated(cdp%diapyc_vel)) then

2705 !$OMP parallel do default(shared)

2706 do j=js,je

2707 do k=2,nz ; do i=is,ie

2708 cdp%diapyc_vel(i,j,k) = idt * (ea(i,j,k) - eb(i,j,k-1))

2709 enddo ; enddo

2710 do i=is,ie

2711 cdp%diapyc_vel(i,j,1) = 0.0

2712 cdp%diapyc_vel(i,j,nz+1) = 0.0

2713 enddo

2714 enddo

2715 endif

2716

2717! For momentum, it is only the net flux that homogenizes within

2718! the mixed layer. Vertical viscosity that is proportional to the

2719! mixed layer turbulence is applied elsewhere.

2720 if (cs%use_bulkmixedlayer) then

2721 if (cs%debug) then

2722 call hchksum(ea, "before net flux rearrangement ea", g%HI, unscale=gv%H_to_MKS)

2723 call hchksum(eb, "before net flux rearrangement eb", g%HI, unscale=gv%H_to_MKS)

2724 endif

2725 !$OMP parallel do default(shared) private(net_ent)

2726 do j=js,je

2727 do k=2,gv%nkml ; do i=is,ie

2728 net_ent = ea(i,j,k) - eb(i,j,k-1)

2729 ea(i,j,k) = max(net_ent, 0.0)

2730 eb(i,j,k-1) = max(-net_ent, 0.0)

2731 enddo ; enddo

2732 enddo

2733 if (cs%debug) then

2734 call hchksum(ea, "after net flux rearrangement ea", g%HI, unscale=gv%H_to_MKS)

2735 call hchksum(eb, "after net flux rearrangement eb", g%HI, unscale=gv%H_to_MKS)

2736 endif

2737 endif

2738

2739! Initialize halo regions of ea, eb, and hold to default values.

2740 !$OMP parallel do default(shared)

2741 do k=1,nz

2742 do i=is-1,ie+1

2743 hold(i,js-1,k) = gv%Angstrom_H ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0

2744 hold(i,je+1,k) = gv%Angstrom_H ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0

2745 enddo

2746 do j=js,je

2747 hold(is-1,j,k) = gv%Angstrom_H ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0

2748 hold(ie+1,j,k) = gv%Angstrom_H ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0

2749 enddo

2750 enddo

2751

2752 call cpu_clock_begin(id_clock_pass)

2753 if (g%symmetric) then ; dir_flag = to_all+omit_corners

2754 else ; dir_flag = to_west+to_south+omit_corners ; endif

2755 call create_group_pass(cs%pass_hold_eb_ea, hold, g%Domain, dir_flag, halo=1)

2756 call create_group_pass(cs%pass_hold_eb_ea, eb, g%Domain, dir_flag, halo=1)

2757 call create_group_pass(cs%pass_hold_eb_ea, ea, g%Domain, dir_flag, halo=1)

2758 call do_group_pass(cs%pass_hold_eb_ea, g%Domain)

2759 call cpu_clock_end(id_clock_pass)

2760

2761 ! Use a tridiagonal solver to determine effect of the diapycnal

2762 ! advection on velocity field. It is assumed that water leaves

2763 ! or enters the ocean with the surface velocity.

2764 if (cs%debug) then

2765 call mom_state_chksum("before u/v tridiag ", u, v, h, g, gv, us, haloshift=0)

2766 call hchksum(ea, "before u/v tridiag ea", g%HI, unscale=gv%H_to_MKS)

2767 call hchksum(eb, "before u/v tridiag eb", g%HI, unscale=gv%H_to_MKS)

2768 call hchksum(hold, "before u/v tridiag hold", g%HI, unscale=gv%H_to_MKS)

2769 endif

2770 call cpu_clock_begin(id_clock_tridiag)

2771

2772 !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)

2773 do j=js,je

2774 do i=isq,ieq

2775 if (associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,1) = u(i,j,1)

2776 hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect

2777 b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))

2778 d1(i) = hval * b1(i)

2779 u(i,j,1) = b1(i) * (hval * u(i,j,1))

2780 enddo

2781 do k=2,nz ; do i=isq,ieq

2782 if (associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,k) = u(i,j,k)

2783 c1(i,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(i)

2784 eaval = ea(i,j,k) + ea(i+1,j,k)

2785 hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect

2786 b1(i) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(i)*eaval))

2787 d1(i) = (hval + d1(i)*eaval) * b1(i)

2788 u(i,j,k) = (hval*u(i,j,k) + eaval*u(i,j,k-1))*b1(i)

2789 enddo ; enddo

2790 do k=nz-1,1,-1 ; do i=isq,ieq

2791 u(i,j,k) = u(i,j,k) + c1(i,k+1)*u(i,j,k+1)

2792 if (associated(adp%du_dt_dia)) &

2793 adp%du_dt_dia(i,j,k) = (u(i,j,k) - adp%du_dt_dia(i,j,k)) * idt

2794 enddo ; enddo

2795 if (associated(adp%du_dt_dia)) then

2796 do i=isq,ieq

2797 adp%du_dt_dia(i,j,nz) = (u(i,j,nz)-adp%du_dt_dia(i,j,nz)) * idt

2798 enddo

2799 endif

2800 enddo

2801 if (cs%debug) then

2802 call mom_state_chksum("aft 1st loop tridiag ", u, v, h, g, gv, us, haloshift=0)

2803 endif

2804 !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)

2805 do j=jsq,jeq

2806 do i=is,ie

2807 if (associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,1) = v(i,j,1)

2808 hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect

2809 b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))

2810 d1(i) = hval * b1(i)

2811 v(i,j,1) = b1(i) * (hval * v(i,j,1))

2812 enddo

2813 do k=2,nz ; do i=is,ie

2814 if (associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,k) = v(i,j,k)

2815 c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)

2816 eaval = ea(i,j,k) + ea(i,j+1,k)

2817 hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect

2818 b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))

2819 d1(i) = (hval + d1(i)*eaval) * b1(i)

2820 v(i,j,k) = (hval*v(i,j,k) + eaval*v(i,j,k-1))*b1(i)

2821 enddo ; enddo

2822 do k=nz-1,1,-1 ; do i=is,ie

2823 v(i,j,k) = v(i,j,k) + c1(i,k+1)*v(i,j,k+1)

2824 if (associated(adp%dv_dt_dia)) &

2825 adp%dv_dt_dia(i,j,k) = (v(i,j,k) - adp%dv_dt_dia(i,j,k)) * idt

2826 enddo ; enddo

2827 if (associated(adp%dv_dt_dia)) then

2828 do i=is,ie

2829 adp%dv_dt_dia(i,j,nz) = (v(i,j,nz)-adp%dv_dt_dia(i,j,nz)) * idt

2830 enddo

2831 endif

2832 enddo

2833 call cpu_clock_end(id_clock_tridiag)

2834 if (cs%debug) then

2835 call mom_state_chksum("after u/v tridiag ", u, v, h, g, gv, us, haloshift=0)

2836 endif

2837

2838 ! Diagnose the diapycnal diffusivities and other related quantities.

2839 if (cs%id_Kd_int > 0) call post_data(cs%id_Kd_int, kd_int, cs%diag)

2840 if (cs%id_Kd_heat > 0) call post_data(cs%id_Kd_heat, kd_heat, cs%diag)

2841 if (cs%id_Kd_salt > 0) call post_data(cs%id_Kd_salt, kd_salt, cs%diag)

2842

2843 if (cs%id_ea > 0) call post_data(cs%id_ea, ea, cs%diag)

2844 if (cs%id_eb > 0) call post_data(cs%id_eb, eb, cs%diag)

2845

2846 if (cs%id_dudt_dia > 0) call post_data(cs%id_dudt_dia, adp%du_dt_dia, cs%diag)

2847 if (cs%id_dvdt_dia > 0) call post_data(cs%id_dvdt_dia, adp%dv_dt_dia, cs%diag)

2848 if (cs%id_wd > 0) call post_data(cs%id_wd, cdp%diapyc_vel, cs%diag)

2849

2850 if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)

2851 if (cs%id_Tadv > 0) call post_data(cs%id_Tadv, tadv_flx, cs%diag)

2852 if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)

2853 if (cs%id_Sadv > 0) call post_data(cs%id_Sadv, sadv_flx, cs%diag)

2854

2855 ! Diagnostics for thickness-weighted vertically averaged diapycnal accelerations

2856 if (cs%id_hf_dudt_dia_2d > 0) &

2857 call post_product_sum_u(cs%id_hf_dudt_dia_2d, adp%du_dt_dia, adp%diag_hfrac_u, g, nz, cs%diag)

2858 if (cs%id_hf_dvdt_dia_2d > 0) &

2859 call post_product_sum_v(cs%id_hf_dvdt_dia_2d, adp%dv_dt_dia, adp%diag_hfrac_v, g, nz, cs%diag)

2860

2861 call disable_averaging(cs%diag)

2862

2863 if (showcalltree) call calltree_leave("layered_diabatic()")

2864

2865end subroutine layered_diabatic

2866

2867!> Returns pointers or values of members within the diabatic_CS type. For extensibility,

2868!! each returned argument is an optional argument

2869subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit, minimum_forcing_depth, &

2870 KPP_CSp, energetic_PBL_CSp, diabatic_aux_CSp, diabatic_halo, use_KPP)

2871 type(diabatic_cs), target, intent(in) :: cs !< module control structure

2872 ! All output arguments are optional

2873 type(opacity_cs), optional, pointer :: opacity_csp !< A pointer to be set to the opacity control structure

2874 type(optics_type), optional, pointer :: optics_csp !< A pointer to be set to the optics control structure

2875 type(kpp_cs), optional, pointer :: kpp_csp !< A pointer to be set to the KPP CS

2876 type(energetic_pbl_cs), optional, pointer :: energetic_pbl_csp !< A pointer to be set to the ePBL CS

2877 real, optional, intent( out) :: evap_cfl_limit !<The largest fraction of a layer that can be

2878 !! evaporated in one time-step [nondim].

2879 real, optional, intent( out) :: minimum_forcing_depth !< The smallest depth over which heat

2880 !! and freshwater fluxes are applied [H ~> m or kg m-2].

2881 type(diabatic_aux_cs), optional, pointer :: diabatic_aux_csp !< A pointer to be set to the diabatic_aux

2882 !! control structure

2883 integer, optional, intent( out) :: diabatic_halo !< The halo size where the diabatic algorithms

2884 !! assume thermodynamics properties are valid.

2885 logical, optional, intent( out) :: use_kpp !< If true, diabatic is using KPP vertical mixing

2886

2887 ! Pointers to control structures

2888 if (present(opacity_csp)) opacity_csp => cs%opacity

2889 if (present(optics_csp)) optics_csp => cs%optics

2890 if (present(kpp_csp)) kpp_csp => cs%KPP_CSp

2891 if (present(energetic_pbl_csp) .and. cs%use_energetic_PBL) energetic_pbl_csp => cs%ePBL

2892 if (present(diabatic_aux_csp)) diabatic_aux_csp => cs%diabatic_aux_CSp

2893

2894 ! Constants within diabatic_CS

2895 if (present(evap_cfl_limit)) evap_cfl_limit = cs%evap_CFL_limit

2896 if (present(minimum_forcing_depth)) minimum_forcing_depth = cs%minimum_forcing_depth

2897 if (present(diabatic_halo)) diabatic_halo = cs%halo_diabatic

2898 if (present(use_kpp)) use_kpp = cs%use_KPP

2899end subroutine extract_diabatic_member

2900

2901!> Routine called for adiabatic physics

2902subroutine adiabatic(h, tv, fluxes, dt, G, GV, US, CS)

2903 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure

2904 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure

2905 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

2906 intent(inout) :: h !< thickness [H ~> m or kg m-2]

2907 type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields

2908 type(forcing), intent(inout) :: fluxes !< boundary fluxes

2909 real, intent(in) :: dt !< time step [T ~> s]

2910 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

2911 type(diabatic_cs), pointer :: cs !< module control structure

2912

2913 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: zeros ! An array of zeros with units of [H ~> m or kg m-2]

2914

2915 zeros(:,:,:) = 0.0

2916

2917 call call_tracer_column_fns(h, h, zeros, zeros, fluxes, zeros(:,:,1), dt, g, gv, us, tv, &

2918 cs%optics, cs%tracer_flow_CSp, cs%debug)

2919

2920end subroutine adiabatic

2921

2922

2923!> This routine diagnoses tendencies from application of diabatic diffusion

2924!! using ALE algorithm. Note that layer thickness is not altered by

2925!! diabatic diffusion.

2926subroutine diagnose_diabatic_diff_tendency(tv, h, temp_old, saln_old, dt, G, GV, US, CS)

2927 type(ocean_grid_type), intent(in) :: G !< ocean grid structure

2928 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

2929 type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields

2930 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< thickness [H ~> m or kg m-2]

2931 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: temp_old !< temperature prior to diabatic

2932 !! physics [C ~> degC]

2933 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: saln_old !< salinity prior to diabatic physics [S ~> ppt]

2934 real, intent(in) :: dt !< time step [T ~> s]

2935 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

2936 type(diabatic_cs), pointer :: CS !< module control structure

2937

2938 ! Local variables

2939 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_3d ! A 3-d work array for diagnostics [various]

2940 real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array for diagnostics [various]

2941 real :: Idt ! The inverse of the timestep [T-1 ~> s-1]

2942 integer :: i, j, k, is, ie, js, je, nz

2943 logical :: do_saln_tend ! Calculate salinity-based tendency diagnostics

2944

2945 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

2946 idt = 0.0 ; if (dt > 0.0) idt = 1. / dt

2947 work_3d(:,:,:) = 0.0

2948 work_2d(:,:) = 0.0

2949

2950

2951 ! temperature tendency

2952 do k=1,nz ; do j=js,je ; do i=is,ie

2953 work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt

2954 enddo ; enddo ; enddo

2955 if (cs%id_diabatic_diff_temp_tend > 0) then

2956 call post_data(cs%id_diabatic_diff_temp_tend, work_3d, cs%diag, alt_h=h)

2957 endif

2958

2959 ! heat tendency

2960 if (cs%id_diabatic_diff_heat_tend > 0 .or. cs%id_diabatic_diff_heat_tend_2d > 0) then

2961 do k=1,nz ; do j=js,je ; do i=is,ie

2962 work_3d(i,j,k) = h(i,j,k)*gv%H_to_RZ * tv%C_p * work_3d(i,j,k)

2963 enddo ; enddo ; enddo

2964 if (cs%id_diabatic_diff_heat_tend > 0) then

2965 call post_data(cs%id_diabatic_diff_heat_tend, work_3d, cs%diag, alt_h=h)

2966 endif

2967 if (cs%id_diabatic_diff_heat_tend_2d > 0) then

2968 work_2d(:,:) = 0.0

2969 do k=1,nz ; do j=js,je ; do i=is,ie

2970 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

2971 enddo ; enddo ; enddo

2972 call post_data(cs%id_diabatic_diff_heat_tend_2d, work_2d, cs%diag)

2973 endif

2974 endif

2975

2976 ! salinity tendency

2977 do_saln_tend = cs%id_diabatic_diff_saln_tend > 0 &

2978 .or. cs%id_diabatic_diff_salt_tend > 0 &

2979 .or. cs%id_diabatic_diff_salt_tend_2d > 0

2980

2981 if (do_saln_tend) then

2982 do k=1,nz ; do j=js,je ; do i=is,ie

2983 work_3d(i,j,k) = (tv%S(i,j,k) - saln_old(i,j,k)) * idt

2984 enddo ; enddo ; enddo

2985

2986 if (cs%id_diabatic_diff_saln_tend > 0) &

2987 call post_data(cs%id_diabatic_diff_saln_tend, work_3d, cs%diag, alt_h=h)

2988

2989 ! salt tendency

2990 if (cs%id_diabatic_diff_salt_tend > 0 .or. cs%id_diabatic_diff_salt_tend_2d > 0) then

2991 do k=1,nz ; do j=js,je ; do i=is,ie

2992 work_3d(i,j,k) = h(i,j,k) * work_3d(i,j,k)

2993 enddo ; enddo ; enddo

2994 if (cs%id_diabatic_diff_salt_tend > 0) then

2995 call post_data(cs%id_diabatic_diff_salt_tend, work_3d, cs%diag, alt_h=h)

2996 endif

2997 if (cs%id_diabatic_diff_salt_tend_2d > 0) then

2998 work_2d(:,:) = 0.0

2999 do k=1,nz ; do j=js,je ; do i=is,ie

3000 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

3001 enddo ; enddo ; enddo

3002 call post_data(cs%id_diabatic_diff_salt_tend_2d, work_2d, cs%diag)

3003 endif

3004 endif

3005 endif

3006

3007end subroutine diagnose_diabatic_diff_tendency

3008

3009

3010!> This routine diagnoses tendencies from application of boundary fluxes.

3011!! These impacts are generally 3d, in particular for penetrative shortwave.

3012!! Other fluxes contribute 3d in cases when the layers vanish or are very thin,

3013!! in which case we distribute the flux into k > 1 layers.

3014subroutine diagnose_boundary_forcing_tendency(tv, h, temp_old, saln_old, h_old, &

3015 dt, G, GV, US, CS)

3016 type(ocean_grid_type), intent(in) :: G !< ocean grid structure

3017 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

3018 type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields

3019 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

3020 intent(in) :: h !< thickness after boundary flux application [H ~> m or kg m-2]

3021 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

3022 intent(in) :: temp_old !< temperature prior to boundary flux application [C ~> degC]

3023 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

3024 intent(in) :: saln_old !< salinity prior to boundary flux application [S ~> ppt]

3025 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

3026 intent(in) :: h_old !< thickness prior to boundary flux application [H ~> m or kg m-2]

3027 real, intent(in) :: dt !< time step [T ~> s]

3028 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

3029 type(diabatic_cs), pointer :: CS !< module control structure

3030

3031 ! Local variables

3032 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_3d ! A 3-d work array for diagnostics [various]

3033 real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array for diagnostics [various]

3034 real :: Idt ! The inverse of the timestep [T-1 ~> s-1]

3035 integer :: i, j, k, is, ie, js, je, nz

3036

3037 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

3038 idt = 0.0 ; if (dt > 0.0) idt = 1. / dt

3039 work_3d(:,:,:) = 0.0

3040 work_2d(:,:) = 0.0

3041

3042 ! Thickness tendency

3043 if (cs%id_boundary_forcing_h_tendency > 0) then

3044 do k=1,nz ; do j=js,je ; do i=is,ie

3045 work_3d(i,j,k) = (h(i,j,k) - h_old(i,j,k))*idt

3046 enddo ; enddo ; enddo

3047 call post_data(cs%id_boundary_forcing_h_tendency, work_3d, cs%diag, alt_h=h_old)

3048 endif

3049

3050 ! temperature tendency

3051 if (cs%id_boundary_forcing_temp_tend > 0) then

3052 do k=1,nz ; do j=js,je ; do i=is,ie

3053 work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt

3054 enddo ; enddo ; enddo

3055 call post_data(cs%id_boundary_forcing_temp_tend, work_3d, cs%diag, alt_h=h_old)

3056 endif

3057

3058 ! heat tendency

3059 if (cs%id_boundary_forcing_heat_tend > 0 .or. cs%id_boundary_forcing_heat_tend_2d > 0) then

3060 do k=1,nz ; do j=js,je ; do i=is,ie

3061 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))

3062 enddo ; enddo ; enddo

3063 if (cs%id_boundary_forcing_heat_tend > 0) then

3064 call post_data(cs%id_boundary_forcing_heat_tend, work_3d, cs%diag, alt_h=h_old)

3065 endif

3066 if (cs%id_boundary_forcing_heat_tend_2d > 0) then

3067 work_2d(:,:) = 0.0

3068 do k=1,nz ; do j=js,je ; do i=is,ie

3069 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

3070 enddo ; enddo ; enddo

3071 call post_data(cs%id_boundary_forcing_heat_tend_2d, work_2d, cs%diag)

3072 endif

3073 endif

3074

3075 ! salinity tendency

3076 if (cs%id_boundary_forcing_saln_tend > 0) then

3077 do k=1,nz ; do j=js,je ; do i=is,ie

3078 work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*idt

3079 enddo ; enddo ; enddo

3080 call post_data(cs%id_boundary_forcing_saln_tend, work_3d, cs%diag, alt_h=h_old)

3081 endif

3082

3083 ! salt tendency

3084 if (cs%id_boundary_forcing_salt_tend > 0 .or. cs%id_boundary_forcing_salt_tend_2d > 0) then

3085 do k=1,nz ; do j=js,je ; do i=is,ie

3086 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))

3087 enddo ; enddo ; enddo

3088 if (cs%id_boundary_forcing_salt_tend > 0) then

3089 call post_data(cs%id_boundary_forcing_salt_tend, work_3d, cs%diag, alt_h=h_old)

3090 endif

3091 if (cs%id_boundary_forcing_salt_tend_2d > 0) then

3092 work_2d(:,:) = 0.0

3093 do k=1,nz ; do j=js,je ; do i=is,ie

3094 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

3095 enddo ; enddo ; enddo

3096 call post_data(cs%id_boundary_forcing_salt_tend_2d, work_2d, cs%diag)

3097 endif

3098 endif

3099

3100end subroutine diagnose_boundary_forcing_tendency

3101

3102

3103!> This routine diagnoses tendencies for temperature and heat from frazil formation.

3104!! This routine is called twice from within subroutine diabatic; at start and at

3105!! end of the diabatic processes. The impacts from frazil are generally a function

3106!! of depth. Hence, when checking heat budget, be sure to remove HFSIFRAZIL from HFDS in k=1.

3107subroutine diagnose_frazil_tendency(tv, h, temp_old, dt, G, GV, US, CS)

3108 type(ocean_grid_type), intent(in) :: G !< ocean grid structure

3109 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure

3110 type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields

3111 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< thickness [H ~> m or kg m-2]

3112 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: temp_old !< temperature prior to frazil formation [C ~> degC]

3113 real, intent(in) :: dt !< time step [T ~> s]

3114 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

3115 type(diabatic_cs), pointer :: CS !< module control structure

3116

3117 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_3d ! A 3-d work array for diagnostics [various]

3118 real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array for diagnostics [various]

3119 real :: Idt ! The inverse of the timestep [T-1 ~> s-1]

3120 integer :: i, j, k, is, ie, js, je, nz

3121

3122 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

3123 idt = 0.0 ; if (dt > 0.0) idt = 1. / dt

3124

3125 ! temperature tendency

3126 if (cs%id_frazil_temp_tend > 0) then

3127 do k=1,nz ; do j=js,je ; do i=is,ie

3128 work_3d(i,j,k) = idt * (tv%T(i,j,k)-temp_old(i,j,k))

3129 enddo ; enddo ; enddo

3130 call post_data(cs%id_frazil_temp_tend, work_3d, cs%diag)

3131 endif

3132

3133 ! heat tendency

3134 if (cs%id_frazil_heat_tend > 0 .or. cs%id_frazil_heat_tend_2d > 0) then

3135 do k=1,nz ; do j=js,je ; do i=is,ie

3136 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))

3137 enddo ; enddo ; enddo

3138 if (cs%id_frazil_heat_tend > 0) call post_data(cs%id_frazil_heat_tend, work_3d, cs%diag)

3139

3140 ! As a consistency check, we must have

3141 ! FRAZIL_HEAT_TENDENCY_2d = HFSIFRAZIL

3142 if (cs%id_frazil_heat_tend_2d > 0) then

3143 work_2d(:,:) = 0.0

3144 do k=1,nz ; do j=js,je ; do i=is,ie

3145 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)

3146 enddo ; enddo ; enddo

3147 call post_data(cs%id_frazil_heat_tend_2d, work_2d, cs%diag)

3148 endif

3149 endif

3150

3151end subroutine diagnose_frazil_tendency

3152

3153

3154!> A simplified version of diabatic_driver_init that will allow

3155!! tracer column functions to be called without allowing any

3156!! of the diabatic processes to be used.

3157subroutine adiabatic_driver_init(Time, G, param_file, diag, CS, &

3158 tracer_flow_CSp)

3159 type(time_type), intent(in) :: time !< current model time

3160 type(ocean_grid_type), intent(in) :: g !< model grid structure

3161 type(param_file_type), intent(in) :: param_file !< the file to parse for parameter values

3162 type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output

3163 type(diabatic_cs), pointer :: cs !< module control structure

3164 type(tracer_flow_control_cs), pointer :: tracer_flow_csp !< pointer to control structure of the

3165 !! tracer flow control module

3166

3167 ! This "include" declares and sets the variable "version".

3168# include "version_variable.h"

3169 character(len=40) :: mdl = "MOM_diabatic_driver" ! This module's name.

3170

3171 if (associated(cs)) then

3172 call mom_error(warning, "adiabatic_driver_init called with an "// &

3173 "associated control structure.")

3174 return

3175 else ; allocate(cs) ; endif

3176

3177 cs%diag => diag

3178 if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp

3179

3180 ! Set default, read and log parameters

3181 call log_version(param_file, mdl, version, &

3182 "The following parameters are used for diabatic processes.")

3183

3184 ! Check for any subsidiary parameters that are inconsistent with the adiabatic mode.

3185 call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &

3186 "If true, sponges may be applied anywhere in the domain. "//&

3187 "The exact location and properties of those sponges are "//&

3188 "specified via calls to initialize_sponge and possibly "//&

3189 "set_up_sponge_field.", default=.false., do_not_log=.true.)

3190 call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", cs%use_energetic_PBL, &

3191 "If true, use an implied energetics planetary boundary "//&

3192 "layer scheme to determine the diffusivity and viscosity "//&

3193 "in the surface boundary layer.", default=.false., do_not_log=.true.)

3194 call get_param(param_file, mdl, "USE_KPP", cs%use_KPP, &

3195 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&

3196 "to calculate diffusivities and non-local transport in the OBL.", &

3197 default=.false., do_not_log=.true.)

3198

3199 if (cs%use_sponge) call mom_error(warning, &

3200 "When ADIABATIC = True, it is inconsistent to set SPONGE = True.")

3201 if (cs%use_energetic_PBL) call mom_error(warning, &

3202 "When ADIABATIC = True, it is inconsistent to set ENERGETICS_SFC_PBL = True.")

3203 if (cs%use_KPP) call mom_error(warning, &

3204 "When ADIABATIC = True, it is inconsistent to set USE_KPP = True.")

3205

3206 if (cs%use_sponge .or. cs%use_energetic_PBL .or. cs%use_KPP) &

3207 call mom_error(fatal, "adiabatic_driver_init is aborting due to inconsistent parameter settings.")

3208

3209end subroutine adiabatic_driver_init

3210

3211

3212!> This routine initializes the diabatic driver module.

3213subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, diag, &

3214 ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, &

3215 ALE_sponge_CSp, oda_incupd_CSp, int_tide_CSp)

3216 type(time_type), target :: time !< model time

3217 type(ocean_grid_type), intent(inout) :: g !< model grid structure

3218 type(verticalgrid_type), intent(in) :: gv !< model vertical grid structure

3219 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

3220 type(param_file_type), intent(in) :: param_file !< file to parse for parameter values

3221 logical, intent(in) :: usealealgorithm !< logical for whether to use ALE remapping

3222 type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output

3223 type(accel_diag_ptrs), intent(inout) :: adp !< pointers to accelerations in momentum equations,

3224 !! to enable diagnostics, like energy budgets

3225 type(cont_diag_ptrs), intent(inout) :: cdp !< pointers to terms in continuity equations

3226 type(diabatic_cs), pointer :: cs !< module control structure

3227 type(tracer_flow_control_cs), pointer :: tracer_flow_csp !< pointer to control structure of the

3228 !! tracer flow control module

3229 type(sponge_cs), pointer :: sponge_csp !< pointer to the sponge module control structure

3230 type(ale_sponge_cs), pointer :: ale_sponge_csp !< pointer to the ALE sponge module control structure

3231 type(oda_incupd_cs), pointer :: oda_incupd_csp !< pointer to the ocean data assimilation incremental

3232 !! update module control structure

3233 type(int_tide_cs), pointer :: int_tide_csp !< pointer to the internal tide structure

3234

3235 ! Local variables

3236 real :: kd ! A diffusivity used in the default for other tracer diffusivities [Z2 T-1 ~> m2 s-1]

3237 logical :: use_temperature

3238 character(len=20) :: en1, en2, en3

3239

3240 ! This "include" declares and sets the variable "version".

3241# include "version_variable.h"

3242 character(len=40) :: mdl = "MOM_diabatic_driver" ! This module's name.

3243 character(len=48) :: thickness_units

3244 logical :: physical_obl_scheme

3245 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nbands

3246 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke

3247 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

3248

3249 cs%initialized = .true.

3250

3251 cs%diag => diag

3252 cs%Time => time

3253

3254 if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp

3255 if (associated(sponge_csp)) cs%sponge_CSp => sponge_csp

3256 if (associated(ale_sponge_csp)) cs%ALE_sponge_CSp => ale_sponge_csp

3257 if (associated(oda_incupd_csp)) cs%oda_incupd_CSp => oda_incupd_csp

3258 if (associated(int_tide_csp)) cs%int_tide_CSp => int_tide_csp

3259

3260 cs%useALEalgorithm = usealealgorithm

3261 cs%use_bulkmixedlayer = (gv%nkml > 0)

3262

3263 ! Set default, read and log parameters

3264 call log_version(param_file, mdl, version, &

3265 "The following parameters are used for diabatic processes.", &

3266 log_to_all=.true., debugging=.true.)

3267 call get_param(param_file, mdl, "USE_LEGACY_DIABATIC_DRIVER", cs%use_legacy_diabatic, &

3268 "If true, use a legacy version of the diabatic subroutine. "//&

3269 "This is temporary and is needed to avoid change in answers.", &

3270 default=.true.)

3271 call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &

3272 "If true, sponges may be applied anywhere in the domain. "//&

3273 "The exact location and properties of those sponges are "//&

3274 "specified via calls to initialize_sponge and possibly "//&

3275 "set_up_sponge_field.", default=.false.)

3276 call get_param(param_file, mdl, "ODA_INCUPD", cs%use_oda_incupd, &

3277 "If true, oda incremental updates will be applied "//&

3278 "everywhere in the domain.", default=.false.)

3279 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &

3280 "If true, temperature and salinity are used as state "//&

3281 "variables.", default=.true.)

3282 call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", cs%use_energetic_PBL, &

3283 "If true, use an implied energetics planetary boundary "//&

3284 "layer scheme to determine the diffusivity and viscosity "//&

3285 "in the surface boundary layer.", default=.false.)

3286 call get_param(param_file, mdl, "EPBL_IS_ADDITIVE", cs%ePBL_is_additive, &

3287 "If true, the diffusivity from ePBL is added to all "//&

3288 "other diffusivities. Otherwise, the larger of kappa-shear "//&

3289 "and ePBL diffusivities are used.", default=.true.)

3290 call get_param(param_file, mdl, "PRANDTL_EPBL", cs%ePBL_Prandtl, &

3291 "The Prandtl number used by ePBL to convert vertical diffusivities into "//&

3292 "viscosities.", default=1.0, units="nondim", do_not_log=.not.cs%use_energetic_PBL)

3293 call get_param(param_file, mdl, "USE_KPP", cs%use_KPP, &

3294 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&

3295 "to calculate diffusivities and non-local transport in the OBL.", &

3296 default=.false., do_not_log=.true.)

3297 cs%use_CVMix_ddiff = cvmix_ddiff_is_used(param_file)

3298

3299 cs%use_kappa_shear = kappa_shear_is_used(param_file)

3300 cs%use_CVMix_shear = cvmix_shear_is_used(param_file)

3301

3302 if (cs%use_bulkmixedlayer) then

3303 call get_param(param_file, mdl, "ML_MIX_FIRST", cs%ML_mix_first, &

3304 "The fraction of the mixed layer mixing that is applied "//&

3305 "before interior diapycnal mixing. 0 by default.", &

3306 units="nondim", default=0.0)

3307 call get_param(param_file, mdl, "NKBL", cs%nkbl, default=2, do_not_log=.true.)

3308 else

3309 cs%ML_mix_first = 0.0

3310 endif

3311 if (use_temperature) then

3312 call get_param(param_file, mdl, "DO_GEOTHERMAL", cs%use_geothermal, &

3313 "If true, apply geothermal heating.", default=.false.)

3314 else

3315 cs%use_geothermal = .false.

3316 endif

3317 call get_param(param_file, mdl, "INTERNAL_TIDES", cs%use_int_tides, &

3318 "If true, use the code that advances a separate set of "//&

3319 "equations for the internal tide energy density.", default=.false.)

3320

3321 call get_param(param_file, mdl, "MASSLESS_MATCH_TARGETS", cs%massless_match_targets, &

3322 "If true, the temperature and salinity of massless layers "//&

3323 "are kept consistent with their target densities. "//&

3324 "Otherwise the properties of massless layers evolve "//&

3325 "diffusively to match massive neighboring layers.", &

3326 default=.true.)

3327

3328 call get_param(param_file, mdl, "AGGREGATE_FW_FORCING", cs%aggregate_FW_forcing, &

3329 "If true, the net incoming and outgoing fresh water fluxes are combined "//&

3330 "and applied as either incoming or outgoing depending on the sign of the net. "//&

3331 "If false, the net incoming fresh water flux is added to the model and "//&

3332 "thereafter the net outgoing is removed from the topmost non-vanished "//&

3333 "layers of the updated state.", default=.true.)

3334

3335 call get_param(param_file, mdl, "DEBUG", cs%debug, &

3336 "If true, write out verbose debugging data.", &

3337 default=.false., debuggingparam=.true.)

3338 call get_param(param_file, mdl, "DEBUG_CONSERVATION", cs%debugConservation, &

3339 "If true, monitor conservation and extrema.", &

3340 default=.false., debuggingparam=.true.)

3341

3342 call get_param(param_file, mdl, "DEBUG_ENERGY_REQ", cs%debug_energy_req, &

3343 "If true, debug the energy requirements.", default=.false., do_not_log=.true.)

3344 call get_param(param_file, mdl, "MIX_BOUNDARY_TRACERS", cs%mix_boundary_tracers, &

3345 "If true, mix the passive tracers in massless layers at "//&

3346 "the bottom into the interior as though a diffusivity of "//&

3347 "KD_MIN_TR were operating.", default=.true.)

3348 call get_param(param_file, mdl, "MIX_BOUNDARY_TRACER_ALE", cs%mix_boundary_tracer_ALE, &

3349 "If true and in ALE mode, mix the passive tracers in massless layers at "//&

3350 "the bottom into the interior as though a diffusivity of "//&

3351 "KD_MIN_TR were operating.", default=.false., do_not_log=.not.cs%useALEalgorithm)

3352

3353 if (cs%mix_boundary_tracers .or. cs%mix_boundary_tracer_ALE) then

3354 call get_param(param_file, mdl, "KD", kd, units="m2 s-1", default=0.0, scale=us%m2_s_to_Z2_T)

3355 call get_param(param_file, mdl, "KD_MIN_TR", cs%Kd_min_tr, &

3356 "A minimal diffusivity that should always be applied to "//&

3357 "tracers, especially in massless layers near the bottom. "//&

3358 "The default is 0.1*KD.", &

3359 units="m2 s-1", default=0.1*kd*us%Z2_T_to_m2_s, scale=gv%m2_s_to_HZ_T)

3360 call get_param(param_file, mdl, "KD_BBL_TR", cs%Kd_BBL_tr, &

3361 "A bottom boundary layer tracer diffusivity that will "//&

3362 "allow for explicitly specified bottom fluxes. The "//&

3363 "entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) "//&

3364 "over the same distance.", &

3365 units="m2 s-1", default=0., scale=gv%m2_s_to_HZ_T*(us%Z_to_m*gv%m_to_H))

3366 ! The scaling factor here is usually equivalent to GV%m2_s_to_HZ_T*GV%Z_to_H.

3367 endif

3368

3369 call get_param(param_file, mdl, "TRACER_TRIDIAG", cs%tracer_tridiag, &

3370 "If true, use the passive tracer tridiagonal solver for T and S", &

3371 default=.false.)

3372

3373 call get_param(param_file, mdl, "MINIMUM_FORCING_DEPTH", cs%minimum_forcing_depth, &

3374 "The smallest depth over which forcing can be applied. This "//&

3375 "only takes effect when near-surface layers become thin "//&

3376 "relative to this scale, in which case the forcing tendencies "//&

3377 "scaled down by distributing the forcing over this depth scale.", &

3378 units="m", default=0.001, scale=gv%m_to_H)

3379 call get_param(param_file, mdl, "EVAP_CFL_LIMIT", cs%evap_CFL_limit, &

3380 "The largest fraction of a layer than can be lost to forcing "//&

3381 "(e.g. evaporation, sea-ice formation) in one time-step. The unused "//&

3382 "mass loss is passed down through the column.", &

3383 units="nondim", default=0.8)

3384

3385 if (cs%use_energetic_PBL .and. .not.cs%useALEalgorithm) &

3386 call mom_error(fatal, "diabatic_driver_init: "//&

3387 "ENERGETICS_SFC_PBL = True is only coded to work when USE_REGRIDDING = True.")

3388

3389 ! Register all available diagnostics for this module.

3390 thickness_units = get_thickness_units(gv)

3391

3392 cs%id_ea_t = register_diag_field('ocean_model', 'ea_t', diag%axesTL, time, &

3393 'Layer (heat) entrainment from above per timestep', 'm', conversion=gv%H_to_m)

3394 cs%id_eb_t = register_diag_field('ocean_model', 'eb_t', diag%axesTL, time, &

3395 'Layer (heat) entrainment from below per timestep', 'm', conversion=gv%H_to_m)

3396 cs%id_ea_s = register_diag_field('ocean_model', 'ea_s', diag%axesTL, time, &

3397 'Layer (salt) entrainment from above per timestep', 'm', conversion=gv%H_to_m)

3398 cs%id_eb_s = register_diag_field('ocean_model', 'eb_s', diag%axesTL, time, &

3399 'Layer (salt) entrainment from below per timestep', 'm', conversion=gv%H_to_m)

3400 ! used by layer diabatic

3401 cs%id_ea = register_diag_field('ocean_model', 'ea', diag%axesTL, time, &

3402 'Layer entrainment from above per timestep', 'm', conversion=gv%H_to_m)

3403 cs%id_eb = register_diag_field('ocean_model', 'eb', diag%axesTL, time, &

3404 'Layer entrainment from below per timestep', 'm', conversion=gv%H_to_m)

3405 if (.not.cs%useALEalgorithm) then

3406 cs%id_wd = register_diag_field('ocean_model', 'wd', diag%axesTi, time, &

3407 'Diapycnal velocity', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3408 if (cs%id_wd > 0) call safe_alloc_ptr(cdp%diapyc_vel,isd,ied,jsd,jed,nz+1)

3409

3410 cs%id_dudt_dia = register_diag_field('ocean_model', 'dudt_dia', diag%axesCuL, time, &

3411 'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)

3412 cs%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, time, &

3413 'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)

3414

3415 cs%id_hf_dudt_dia_2d = register_diag_field('ocean_model', 'hf_dudt_dia_2d', diag%axesCu1, time, &

3416 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Diapycnal Mixing', &

3417 'm s-2', conversion=us%L_T2_to_m_s2)

3418 if (cs%id_hf_dudt_dia_2d > 0) call safe_alloc_ptr(adp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)

3419

3420 cs%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, time, &

3421 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', &

3422 'm s-2', conversion=us%L_T2_to_m_s2)

3423 if (cs%id_hf_dvdt_dia_2d > 0) call safe_alloc_ptr(adp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)

3424

3425 if ((cs%id_dudt_dia > 0) .or. (cs%id_hf_dudt_dia_2d > 0)) &

3426 call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)

3427 if ((cs%id_dvdt_dia > 0) .or. (cs%id_hf_dvdt_dia_2d > 0)) &

3428 call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)

3429 endif

3430

3431 if (use_temperature) then

3432 cs%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff", diag%axesTi, &

3433 time, "Diffusive diapycnal temperature flux across interfaces", &

3434 units="degC m s-1", conversion=us%C_to_degC*gv%H_to_m*us%s_to_T)

3435 if (.not.cs%useALEalgorithm) then

3436 cs%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv", diag%axesTi, &

3437 time, "Advective diapycnal temperature flux across interfaces", &

3438 units="degC m s-1", conversion=us%C_to_degC*gv%H_to_m*us%s_to_T)

3439 endif

3440 cs%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff", diag%axesTi, &

3441 time, "Diffusive diapycnal salnity flux across interfaces", &

3442 units="psu m s-1", conversion=us%S_to_ppt*gv%H_to_m*us%s_to_T)

3443 if (.not.cs%useALEalgorithm) then

3444 cs%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv", diag%axesTi, &

3445 time, "Advective diapycnal salnity flux across interfaces", &

3446 units="psu m s-1", conversion=us%S_to_ppt*gv%H_to_m*us%s_to_T)

3447 endif

3448

3449 cs%id_N2_dd = register_diag_field('ocean_model',"N2_diabatic", diag%axesTi, &

3450 time, "Squared buoyancy frequency diagnosed after diffusion applied in thermodynamic timestep.", &

3451 "s-2", conversion=us%s_to_T**2)

3452 cs%id_N2_temp_dd = register_diag_field('ocean_model',"N2_temp_diabatic", diag%axesTi, &

3453 time, "Squared buoyancy frequency due to temperature stratification diagnosed after diffusion applied "//&

3454 "in thermodynamic timestep.", "s-2", conversion=us%s_to_T**2)

3455 cs%id_N2_salt_dd = register_diag_field('ocean_model',"N2_salt_diabatic", diag%axesTi, &

3456 time, "Squared buoyancy frequency due to salinity stratification diagnosed after diffusion applied "//&

3457 "in thermodynamic timestep.", "s-2", conversion=us%s_to_T**2)

3458 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.

3459

3460 cs%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, time, &

3461 'Mixed layer depth (delta rho = 0.03)', units='m', conversion=us%Z_to_m, &

3462 cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', &

3463 cmor_standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t')

3464 cs%id_mlotstsq = register_diag_field('ocean_model', 'mlotstsq', diag%axesT1, time, &

3465 long_name='Square of Ocean Mixed Layer Thickness Defined by Sigma T', &

3466 standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t', &

3467 units='m2', conversion=us%Z_to_m**2)

3468 cs%id_MLD_0125 = register_diag_field('ocean_model', 'MLD_0125', diag%axesT1, time, &

3469 'Mixed layer depth (delta rho = 0.125)', 'm', conversion=us%Z_to_m)

3470 call get_param(param_file, mdl, "MLD_EN_VALS", cs%MLD_En_vals, &

3471 "The energy values used to compute MLDs. If not set (or all set to 0.), the "//&

3472 "default will overwrite to 25., 2500., 250000.", units='J/m2', &

3473 defaults=(/25., 2500., 250000./), scale=us%W_m2_to_RZ3_T3*us%s_to_T)

3474 write(en1,'(F10.2)') cs%MLD_En_vals(1)*us%RZ3_T3_to_W_m2*us%T_to_s

3475 write(en2,'(F10.2)') cs%MLD_En_vals(2)*us%RZ3_T3_to_W_m2*us%T_to_s

3476 write(en3,'(F10.2)') cs%MLD_En_vals(3)*us%RZ3_T3_to_W_m2*us%T_to_s

3477 cs%id_MLD_EN1 = register_diag_field('ocean_model', 'MLD_EN1', diag%axesT1, time, &

3478 'Mixed layer depth for energy value set to '//trim(en1)//' J/m2 (Energy set by 1st MLD_EN_VALS)', &

3479 units='m', conversion=us%Z_to_m)

3480 cs%id_MLD_EN2 = register_diag_field('ocean_model', 'MLD_EN2', diag%axesT1, time, &

3481 'Mixed layer depth for energy value set to '//trim(en2)//' J/m2 (Energy set by 2nd MLD_EN_VALS)', &

3482 units='m', conversion=us%Z_to_m)

3483 cs%id_MLD_EN3 = register_diag_field('ocean_model', 'MLD_EN3', diag%axesT1, time, &

3484 'Mixed layer depth for energy value set to '//trim(en3)//' J/m2 (Energy set by 3rd MLD_EN_VALS)', &

3485 units='m', conversion=us%Z_to_m)

3486 call get_param(param_file, mdl, "BMLD_EN_VALS", cs%BMLD_En_vals, &

3487 "The energy values used to compute Bottom MLDs. If not set (or all set to 0.), the "//&

3488 "default will overwrite to 25., 2500., 250000.", units='J/m2', &

3489 defaults=(/25., 2500., 250000./), scale=us%W_m2_to_RZ3_T3*us%s_to_T)

3490 write(en1,'(F10.2)') cs%BMLD_En_vals(1)*us%RZ3_T3_to_W_m2*us%T_to_s

3491 write(en2,'(F10.2)') cs%BMLD_En_vals(2)*us%RZ3_T3_to_W_m2*us%T_to_s

3492 write(en3,'(F10.2)') cs%BMLD_En_vals(3)*us%RZ3_T3_to_W_m2*us%T_to_s

3493 cs%id_BMLD_EN1 = register_diag_field('ocean_model', 'BMLD_EN1', diag%axesT1, time, &

3494 'Bottom mixed layer depth for energy value set to '//trim(en1)//' J/m2 (Energy set by 1st MLD_EN_VALS)', &

3495 units='m', conversion=us%Z_to_m)

3496 cs%id_BMLD_EN2 = register_diag_field('ocean_model', 'BMLD_EN2', diag%axesT1, time, &

3497 'Bottom mixed layer depth for energy value set to '//trim(en2)//' J/m2 (Energy set by 2nd MLD_EN_VALS)', &

3498 units='m', conversion=us%Z_to_m)

3499 cs%id_BMLD_EN3 = register_diag_field('ocean_model', 'BMLD_EN3', diag%axesT1, time, &

3500 'Bottom mixed layer depth for energy value set to '//trim(en3)//' J/m2 (Energy set by 3rd MLD_EN_VALS)', &

3501 units='m', conversion=us%Z_to_m)

3502 if ((cs%id_MLD_EN1>0).or. (cs%id_MLD_EN2>0).or. (cs%id_MLD_EN3>0).or. &

3503 (cs%id_BMLD_EN1>0).or.(cs%id_BMLD_EN2>0).or.(cs%id_BMLD_EN3>0)) then

3504 call get_param(param_file, mdl, "USE_OM4_MLD_EN_ITER", cs%use_OM4_MLD_En_iter, &

3505 "If true, uses an older set of iteration coefficients in computing the PE based "//&

3506 "surface MLD to reproduce OM4 era models. False uses an updated (general) method.",&

3507 default=.true.)

3508 endif

3509 cs%id_subMLN2 = register_diag_field('ocean_model', 'subML_N2', diag%axesT1, time, &

3510 'Squared buoyancy frequency below mixed layer', units='s-2', conversion=us%s_to_T**2)

3511 cs%id_MLD_user = register_diag_field('ocean_model', 'MLD_user', diag%axesT1, time, &

3512 'Mixed layer depth (used defined)', units='m', conversion=us%Z_to_m)

3513 if (cs%id_MLD_003 > 0) then

3514 call get_param(param_file, mdl, "HREF_FOR_MLD", cs%ref_h_mld, &

3515 "Reference depth used to calculate the potential density used to find the mixed layer depth "//&

3516 "based on a delta rho = 0.03 kg/m3.", units='m', default=0.0, scale=us%m_to_Z)

3517 cs%id_MLD_003_zr = register_diag_field('ocean_model', 'MLD_003_refZ', diag%axesT1, time, &

3518 'Depth of reference density for MLD (delta rho = 0.03)', units='m', conversion=us%Z_to_m)

3519 cs%id_MLD_003_rr = register_diag_field('ocean_model', 'MLD_003_refRho', diag%axesT1, time, &

3520 'Reference density for MLD (delta rho = 0.03)', units='kg/m3', conversion=us%R_to_kg_m3)

3521 endif

3522 endif

3523

3524 call kdwork_init(time, g,gv,us,diag,cs%VBF,cs%Use_KdWork_diag)

3525 if (cs%Use_KdWork_diag.and.(.not.usealealgorithm)) &

3526 call mom_error(warning,"The KdWork diagnostics are not fully implemented for use in layer mode.")

3527 if (cs%Use_KdWork_diag.and.(cs%use_legacy_diabatic)) &

3528 call mom_error(warning,"The KdWork diagnostics are only approximate with the legacy diabatic driver.")

3529

3530 call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", cs%MLDdensityDifference, &

3531 "The density difference used to determine a diagnostic mixed "//&

3532 "layer depth, MLD_user, following the definition of Levitus 1982. "//&

3533 "The MLD is the depth at which the density is larger than the "//&

3534 "surface density by the specified amount.", &

3535 units='kg/m3', default=0.1, scale=us%kg_m3_to_R)

3536 call get_param(param_file, mdl, "DIAG_DEPTH_SUBML_N2", cs%dz_subML_N2, &

3537 "The distance over which to calculate a diagnostic of the "//&

3538 "stratification at the base of the mixed layer.", &

3539 units='m', default=50.0, scale=us%m_to_Z)

3540

3541 ! diagnostics for values prior to diabatic and prior to ALE

3542 cs%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, time, &

3543 'Zonal velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)

3544 cs%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, time, &

3545 'Meridional velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)

3546 cs%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, time, &

3547 'Layer Thickness before diabatic forcing', &

3548 trim(thickness_units), conversion=gv%H_to_MKS, v_extensive=.true.)

3549 cs%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, time, &

3550 'Interface Heights before diabatic forcing', 'm', conversion=us%Z_to_m)

3551 if (use_temperature) then

3552 cs%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, time, &

3553 'Potential Temperature', 'degC', conversion=us%C_to_degC)

3554 cs%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, time, &

3555 'Salinity', 'PSU', conversion=us%S_to_ppt)

3556 endif

3557

3558 cs%id_Kd_int = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, time, &

3559 'Total diapycnal diffusivity at interfaces', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)

3560 if (cs%use_energetic_PBL) then

3561 cs%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, time, &

3562 'ePBL diapycnal diffusivity at interfaces', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)

3563 endif

3564

3565 cs%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, time, &

3566 'Total diapycnal diffusivity for heat at interfaces', 'm2 s-1', conversion=gv%HZ_T_to_m2_s, &

3567 cmor_field_name='difvho', &

3568 cmor_standard_name='ocean_vertical_heat_diffusivity', &

3569 cmor_long_name='Ocean vertical heat diffusivity')

3570 cs%id_Kd_salt = register_diag_field('ocean_model', 'Kd_salt', diag%axesTi, time, &

3571 'Total diapycnal diffusivity for salt at interfaces', 'm2 s-1', conversion=gv%HZ_T_to_m2_s, &

3572 cmor_field_name='difvso', &

3573 cmor_standard_name='ocean_vertical_salt_diffusivity', &

3574 cmor_long_name='Ocean vertical salt diffusivity')

3575

3576 ! CS%useKPP is set to True if KPP-scheme is to be used, False otherwise.

3577 ! KPP_init() allocated CS%KPP_Csp and also sets CS%KPPisPassive

3578 cs%useKPP = kpp_init(param_file, g, gv, us, diag, time, cs%KPP_CSp, passive=cs%KPPisPassive)

3579

3580 ! Diagnostics for tendencies of temperature and salinity due to diabatic processes,

3581 ! available only for ALE algorithm.

3582 ! Diagnostics for tendencies of temperature and heat due to frazil

3583 cs%id_diabatic_diff_h = register_diag_field('ocean_model', 'diabatic_diff_h', diag%axesTL, time, &

3584 'Cell thickness used during diabatic diffusion', &

3585 thickness_units, conversion=gv%H_to_MKS, v_extensive=.true.)

3586

3587 if (cs%useALEalgorithm) then

3588 cs%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', &

3589 'diabatic_diff_temp_tendency', diag%axesTL, time, &

3590 'Diabatic diffusion temperature tendency', 'degC s-1', conversion=us%C_to_degC*us%s_to_T)

3591 if (cs%id_diabatic_diff_temp_tend > 0) then

3592 cs%diabatic_diff_tendency_diag = .true.

3593 endif

3594

3595 cs%id_diabatic_diff_saln_tend = register_diag_field('ocean_model',&

3596 'diabatic_diff_saln_tendency', diag%axesTL, time, &

3597 'Diabatic diffusion salinity tendency', 'psu s-1', conversion=us%S_to_ppt*us%s_to_T)

3598 if (cs%id_diabatic_diff_saln_tend > 0) then

3599 cs%diabatic_diff_tendency_diag = .true.

3600 endif

3601

3602 cs%id_diabatic_diff_heat_tend = register_diag_field('ocean_model', &

3603 'diabatic_heat_tendency', diag%axesTL, time, &

3604 'Diabatic diffusion heat tendency', &

3605 'W m-2', conversion=us%QRZ_T_to_W_m2, cmor_field_name='opottempdiff', &

3606 cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'// &

3607 'due_to_parameterized_dianeutral_mixing', &

3608 cmor_long_name='Tendency of sea water potential temperature expressed as heat content '// &

3609 'due to parameterized dianeutral mixing', &

3610 v_extensive=.true.)

3611 if (cs%id_diabatic_diff_heat_tend > 0) then

3612 cs%diabatic_diff_tendency_diag = .true.

3613 endif

3614

3615 cs%id_diabatic_diff_salt_tend = register_diag_field('ocean_model', &

3616 'diabatic_salt_tendency', diag%axesTL, time, &

3617 'Diabatic diffusion of salt tendency', &

3618 'kg m-2 s-1', conversion=us%S_to_ppt*0.001*gv%H_to_RZ*us%RZ_T_to_kg_m2s, &

3619 cmor_field_name='osaltdiff', &

3620 cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &

3621 'due_to_parameterized_dianeutral_mixing', &

3622 cmor_long_name='Tendency of sea water salinity expressed as salt content '// &

3623 'due to parameterized dianeutral mixing', &

3624 v_extensive=.true.)

3625 if (cs%id_diabatic_diff_salt_tend > 0) then

3626 cs%diabatic_diff_tendency_diag = .true.

3627 endif

3628

3629 ! This diagnostic should equal to roundoff if all is working well.

3630 cs%id_diabatic_diff_heat_tend_2d = register_diag_field('ocean_model', &

3631 'diabatic_heat_tendency_2d', diag%axesT1, time, &

3632 'Depth integrated diabatic diffusion heat tendency', &

3633 'W m-2', conversion=us%QRZ_T_to_W_m2, cmor_field_name='opottempdiff_2d', &

3634 cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'//&

3635 'due_to_parameterized_dianeutral_mixing_depth_integrated', &

3636 cmor_long_name='Tendency of sea water potential temperature expressed as heat content '//&

3637 'due to parameterized dianeutral mixing depth integrated')

3638 if (cs%id_diabatic_diff_heat_tend_2d > 0) then

3639 cs%diabatic_diff_tendency_diag = .true.

3640 endif

3641

3642 ! This diagnostic should equal to roundoff if all is working well.

3643 cs%id_diabatic_diff_salt_tend_2d = register_diag_field('ocean_model', &

3644 'diabatic_salt_tendency_2d', diag%axesT1, time, &

3645 'Depth integrated diabatic diffusion salt tendency', &

3646 'kg m-2 s-1', conversion=us%S_to_ppt*0.001*gv%H_to_RZ*us%RZ_T_to_kg_m2s, &

3647 cmor_field_name='osaltdiff_2d', &

3648 cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &

3649 'due_to_parameterized_dianeutral_mixing_depth_integrated', &

3650 cmor_long_name='Tendency of sea water salinity expressed as salt content '// &

3651 'due to parameterized dianeutral mixing depth integrated')

3652 if (cs%id_diabatic_diff_salt_tend_2d > 0) then

3653 cs%diabatic_diff_tendency_diag = .true.

3654 endif

3655

3656 ! Diagnostics for tendencies of thickness temperature and salinity due to boundary forcing,

3657 ! available only for ALE algorithm.

3658 ! Diagnostics for tendencies of temperature and heat due to frazil

3659 cs%id_boundary_forcing_h = register_diag_field('ocean_model', 'boundary_forcing_h', diag%axesTL, time, &

3660 'Cell thickness after applying boundary forcing', &

3661 thickness_units, conversion=gv%H_to_MKS, v_extensive=.true.)

3662 cs%id_boundary_forcing_h_tendency = register_diag_field('ocean_model', &

3663 'boundary_forcing_h_tendency', diag%axesTL, time, &

3664 'Cell thickness tendency due to boundary forcing', &

3665 trim(thickness_units)//" s-1", conversion=gv%H_to_MKS*us%s_to_T, v_extensive=.true.)

3666 if (cs%id_boundary_forcing_h_tendency > 0) then

3667 cs%boundary_forcing_tendency_diag = .true.

3668 endif

3669

3670 cs%id_boundary_forcing_temp_tend = register_diag_field('ocean_model',&

3671 'boundary_forcing_temp_tendency', diag%axesTL, time, &

3672 'Boundary forcing temperature tendency', 'degC s-1', conversion=us%C_to_degC*us%s_to_T)

3673 if (cs%id_boundary_forcing_temp_tend > 0) then

3674 cs%boundary_forcing_tendency_diag = .true.

3675 endif

3676

3677 cs%id_boundary_forcing_saln_tend = register_diag_field('ocean_model',&

3678 'boundary_forcing_saln_tendency', diag%axesTL, time, &

3679 'Boundary forcing saln tendency', 'psu s-1', conversion=us%S_to_ppt*us%s_to_T)

3680 if (cs%id_boundary_forcing_saln_tend > 0) then

3681 cs%boundary_forcing_tendency_diag = .true.

3682 endif

3683

3684 cs%id_boundary_forcing_heat_tend = register_diag_field('ocean_model',&

3685 'boundary_forcing_heat_tendency', diag%axesTL, time, &

3686 'Boundary forcing heat tendency', &

3687 'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive=.true.)

3688 if (cs%id_boundary_forcing_heat_tend > 0) then

3689 cs%boundary_forcing_tendency_diag = .true.

3690 endif

3691

3692 cs%id_boundary_forcing_salt_tend = register_diag_field('ocean_model',&

3693 'boundary_forcing_salt_tendency', diag%axesTL, time, &

3694 'Boundary forcing salt tendency', &

3695 'kg m-2 s-1', conversion=us%S_to_ppt*0.001*gv%H_to_RZ*us%RZ_T_to_kg_m2s, &

3696 v_extensive = .true.)

3697 if (cs%id_boundary_forcing_salt_tend > 0) then

3698 cs%boundary_forcing_tendency_diag = .true.

3699 endif

3700

3701 ! This diagnostic should equal to surface heat flux if all is working well.

3702 cs%id_boundary_forcing_heat_tend_2d = register_diag_field('ocean_model',&

3703 'boundary_forcing_heat_tendency_2d', diag%axesT1, time, &

3704 'Depth integrated boundary forcing of ocean heat', &

3705 'W m-2', conversion=us%QRZ_T_to_W_m2)

3706 if (cs%id_boundary_forcing_heat_tend_2d > 0) then

3707 cs%boundary_forcing_tendency_diag = .true.

3708 endif

3709

3710 ! This diagnostic should equal to surface salt flux if all is working well.

3711 cs%id_boundary_forcing_salt_tend_2d = register_diag_field('ocean_model',&

3712 'boundary_forcing_salt_tendency_2d', diag%axesT1, time, &

3713 'Depth integrated boundary forcing of ocean salt', &

3714 'kg m-2 s-1', conversion=us%S_to_ppt*0.001*gv%H_to_RZ*us%RZ_T_to_kg_m2s)

3715 if (cs%id_boundary_forcing_salt_tend_2d > 0) then

3716 cs%boundary_forcing_tendency_diag = .true.

3717 endif

3718 endif

3719

3720 ! diagnostics for tendencies of temp and heat due to frazil

3721 cs%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, time, &

3722 long_name='Cell Thickness', standard_name='cell_thickness', &

3723 units=thickness_units, conversion=gv%H_to_MKS, v_extensive=.true.)

3724

3725 ! diagnostic for tendency of temp due to frazil

3726 cs%id_frazil_temp_tend = register_diag_field('ocean_model',&

3727 'frazil_temp_tendency', diag%axesTL, time, &

3728 'Temperature tendency due to frazil formation', 'degC s-1', conversion=us%C_to_degC*us%s_to_T)

3729 if (cs%id_frazil_temp_tend > 0) then

3730 cs%frazil_tendency_diag = .true.

3731 endif

3732

3733 ! diagnostic for tendency of heat due to frazil

3734 cs%id_frazil_heat_tend = register_diag_field('ocean_model',&

3735 'frazil_heat_tendency', diag%axesTL, time, &

3736 'Heat tendency due to frazil formation', &

3737 'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive=.true.)

3738 if (cs%id_frazil_heat_tend > 0) then

3739 cs%frazil_tendency_diag = .true.

3740 endif

3741

3742 ! If all is working properly, this diagnostic should equal to hfsifrazil.

3743 cs%id_frazil_heat_tend_2d = register_diag_field('ocean_model',&

3744 'frazil_heat_tendency_2d', diag%axesT1, time, &

3745 'Depth integrated heat tendency due to frazil formation', &

3746 'W m-2', conversion=us%QRZ_T_to_W_m2)

3747 if (cs%id_frazil_heat_tend_2d > 0) then

3748 cs%frazil_tendency_diag = .true.

3749 endif

3750

3751 ! CS%use_CVMix_conv is set to True if CVMix convection will be used, otherwise it is False.

3752 cs%use_CVMix_conv = cvmix_conv_init(time, g, gv, us, param_file, diag, cs%CVMix_conv)

3753

3754 call entrain_diffusive_init(time, g, gv, us, param_file, diag, cs%entrain_diffusive, &

3755 just_read_params=cs%useALEalgorithm)

3756

3757 ! initialize the geothermal heating module

3758 if (cs%use_geothermal) &

3759 call geothermal_init(time, g, gv, us, param_file, diag, cs%geothermal, usealealgorithm)

3760

3761 ! initialize module for internal tide induced mixing

3762 if (cs%use_int_tides) then

3763 call int_tide_input_init(time, g, gv, us, param_file, diag, cs%int_tide_input_CSp, &

3764 cs%int_tide_input)

3765 call internal_tides_init(time, g, gv, us, param_file, diag, int_tide_csp)

3766 endif

3767

3768 !if (associated(int_tide_CSp)) CS%int_tide_CSp => int_tide_CSp

3769

3770 physical_obl_scheme = (cs%use_bulkmixedlayer .or. cs%use_KPP .or. cs%use_energetic_PBL)

3771 ! initialize module for setting diffusivities

3772 call set_diffusivity_init(time, g, gv, us, param_file, diag, cs%set_diff_CSp, cs%int_tide_CSp, &

3773 halo_ts=cs%halo_TS_diff, double_diffuse=cs%double_diffuse, &

3774 physical_obl_scheme=physical_obl_scheme)

3775

3776 cs%halo_diabatic = cs%halo_TS_diff

3777 if (cs%use_int_tides) cs%halo_diabatic = max(cs%halo_TS_diff, 2)

3778

3779 if (cs%useKPP .and. (cs%double_diffuse .and. .not.cs%use_CVMix_ddiff)) &

3780 call mom_error(fatal, 'diabatic_driver_init: DOUBLE_DIFFUSION (old method) does not work '//&

3781 'with KPP. Please set DOUBLE_DIFFUSION=False and USE_CVMIX_DDIFF=True.')

3782

3783 ! set up the clocks for this module

3784 id_clock_entrain = cpu_clock_id('(Ocean diabatic entrain)', grain=clock_module)

3785 if (cs%use_bulkmixedlayer) &

3786 id_clock_mixedlayer = cpu_clock_id('(Ocean mixed layer)', grain=clock_module)

3787 id_clock_remap = cpu_clock_id('(Ocean vert remap)', grain=clock_module)

3788 if (cs%use_geothermal) &

3789 id_clock_geothermal = cpu_clock_id('(Ocean geothermal)', grain=clock_routine)

3790 id_clock_set_diffusivity = cpu_clock_id('(Ocean set_diffusivity)', grain=clock_module)

3791 id_clock_kpp = cpu_clock_id('(Ocean KPP)', grain=clock_module)

3792 id_clock_tracers = cpu_clock_id('(Ocean tracer_columns)', grain=clock_module_driver+5)

3793 if (cs%use_sponge) &

3794 id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=clock_module)

3795 if (cs%use_oda_incupd) &

3796 id_clock_oda_incupd = cpu_clock_id('(Ocean inc. update data assimilation)', grain=clock_module)

3797 id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=clock_routine)

3798 id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=clock_routine)

3799 id_clock_differential_diff = -1 ; if (cs%double_diffuse .and. .not.cs%use_CVMix_ddiff) &

3800 id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=clock_routine)

3801

3802 ! initialize the auxiliary diabatic driver module

3803 call diabatic_aux_init(time, g, gv, us, param_file, diag, cs%diabatic_aux_CSp, &

3804 cs%useALEalgorithm, cs%use_energetic_PBL)

3805

3806 ! initialize the boundary layer modules

3807 if (cs%use_bulkmixedlayer) &

3808 call bulkmixedlayer_init(time, g, gv, us, param_file, diag, cs%bulkmixedlayer)

3809 if (cs%use_energetic_PBL) &

3810 call energetic_pbl_init(time, g, gv, us, param_file, diag, cs%ePBL)

3811

3812 call regularize_layers_init(time, g, gv, param_file, diag, cs%regularize_layers)

3813

3814 if (cs%debug_energy_req) &

3815 call diapyc_energy_req_init(time, g, gv, us, param_file, diag, cs%diapyc_en_rec_CSp)

3816

3817 ! obtain information about the number of bands for penetrative shortwave

3818 if (use_temperature) then

3819 call get_param(param_file, mdl, "PEN_SW_NBANDS", nbands, default=1)

3820 if (nbands > 0) then

3821 allocate(cs%optics)

3822 call opacity_init(time, g, gv, us, param_file, diag, cs%opacity, cs%optics)

3823 endif

3824 endif

3825

3826 ! Initialize the diagnostic grid storage

3827 call diag_grid_storage_init(cs%diag_grids_prev, g, gv, diag)

3828

3829end subroutine diabatic_driver_init

3830

3831!> Routine to register restarts, pass-through to children modules

3832subroutine register_diabatic_restarts(G, GV, US, param_file, int_tide_CSp, restart_CSp, CS)

3833 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure

3834 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure

3835 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

3836 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters

3837 type(int_tide_cs), pointer :: int_tide_csp !< Internal tide control structure

3838 type(mom_restart_cs), pointer :: restart_csp !< MOM restart control structure

3839 type(diabatic_cs), pointer :: cs !< module control structure

3840

3841 logical :: use_int_tides

3842

3843 if (associated(cs)) then

3844 call mom_error(warning, "diabatic_driver_init called with an "// &

3845 "associated control structure.")

3846 return

3847 else

3848 allocate(cs)

3849 endif

3850

3851 use_int_tides=.false.

3852

3853 call read_param(param_file, "INTERNAL_TIDES", use_int_tides)

3854

3855 if (use_int_tides) then

3856 call register_int_tide_restarts(g, gv, us, param_file, int_tide_csp, restart_csp)

3857 endif

3858

3859 call register_kpp_restarts(g, param_file, restart_csp, cs%KPP_CSp)

3860

3861end subroutine register_diabatic_restarts

3862

3863!> Routine to close the diabatic driver module

3864subroutine diabatic_driver_end(CS)

3865 type(diabatic_cs), intent(inout) :: cs !< module control structure

3866

3867 if (associated(cs%VBF)) then

3868 call kdwork_end(cs%VBF)

3869 endif

3870

3871 if (associated(cs%optics)) then

3872 call opacity_end(cs%opacity, cs%optics)

3873 deallocate(cs%optics)

3874 endif

3875

3876 if (cs%debug_energy_req) &

3877 call diapyc_energy_req_end(cs%diapyc_en_rec_CSp)

3878

3879 if (cs%use_energetic_PBL) &

3880 call energetic_pbl_end(cs%ePBL)

3881

3882 call diabatic_aux_end(cs%diabatic_aux_CSp)

3883

3884 call set_diffusivity_end(cs%set_diff_CSp)

3885

3886 deallocate(cs%set_diff_CSp)

3887

3888 if (cs%use_geothermal) &

3889 call geothermal_end(cs%geothermal)

3890

3891 if (cs%useKPP) &

3892 call kpp_end(cs%KPP_CSp)

3893

3894 ! GMM, the following is commented out because arrays in

3895 ! CS%diag_grids_prev are neither pointers or allocatables

3896 ! and, therefore, cannot be deallocated.

3897

3898 !call diag_grid_storage_end(CS%diag_grids_prev)

3899end subroutine diabatic_driver_end

3900

3901

3902!> \namespace mom_diabatic_driver

3903!!

3904!! By Robert Hallberg, Alistair Adcroft, and Stephen Griffies

3905!!

3906!! This program contains the subroutine that, along with the

3907!! subroutines that it calls, implements diapycnal mass and momentum

3908!! fluxes and a bulk mixed layer. The diapycnal diffusion can be

3909!! used without the bulk mixed layer.

3910!!

3911!! \section section_diabatic Outline of MOM diabatic

3912!!

3913!! * diabatic first determines the (diffusive) diapycnal mass fluxes

3914!! based on the convergence of the buoyancy fluxes within each layer.

3915!!

3916!! * The dual-stream entrainment scheme of MacDougall and Dewar (JPO,

3917!! 1997) is used for combined diapycnal advection and diffusion,

3918!! calculated implicitly and potentially with the Richardson number

3919!! dependent mixing, as described by Hallberg (MWR, 2000).

3920!!

3921!! * Diapycnal advection is the residual of diapycnal diffusion,

3922!! so the fully implicit upwind differencing scheme that is used is

3923!! entirely appropriate.

3924!!

3925!! * The downward buoyancy flux in each layer is determined from

3926!! an implicit calculation based on the previously

3927!! calculated flux of the layer above and an estimated flux in the

3928!! layer below. This flux is subject to the following conditions:

3929!! (1) the flux in the top and bottom layers are set by the boundary

3930!! conditions, and (2) no layer may be driven below a minimal thickness.

3931!! If there is a bulk mixed layer, the buffer layer is treated

3932!! as a fixed density layer with vanishingly small diffusivity.

3933!!

3934!! diabatic takes 5 arguments: the two velocities (u and v), the

3935!! thicknesses (h), a structure containing the forcing fields, and

3936!! the length of time over which to act (dt). The velocities and

3937!! thickness are taken as inputs and modified within the subroutine.

3938!! There is no limit on the time step.

3939

3940end module mom_diabatic_driver