MOM_diapyc_energy_req.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!> Calculates the energy requirements of mixing.
7
8!! \author By Robert Hallberg, May 2015
9
10use mom_diag_mediator, only : diag_ctrl, time_type, post_data, register_diag_field
11use mom_eos, only : calculate_specific_vol_derivs, calculate_density
12use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
13use mom_file_parser, only : get_param, log_version, param_file_type
14use mom_grid, only : ocean_grid_type
15use mom_interface_heights, only : thickness_to_dz
19
20implicit none ; private
21
23
24! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
25! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
26! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
27! vary with the Boussinesq approximation, the Boussinesq variant is given first.
28
29!> This control structure holds parameters for the MOM_diapyc_energy_req module
30type, public :: diapyc_energy_req_cs ; private
31 logical :: initialized = .false. !< A variable that is here because empty
32 !! structures are not permitted by some compilers.
33 real :: test_kh_scaling !< A scaling factor for the diapycnal diffusivity [nondim]
34 real :: colht_scaling !< A scaling factor for the column height change correction term [nondim]
35 real :: vonkar !< The von Karman coefficient as used in this module [nondim]
36 logical :: use_test_kh_profile !< If true, use the internal test diffusivity profile in place of
37 !! any that might be passed in as an argument.
38 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
39 !! regulate the timing of diagnostic output.
40
41 !>@{ Diagnostic IDs
42 integer :: id_ert=-1, id_erb=-1, id_erc=-1, id_erh=-1, id_kddt=-1, id_kd=-1
43 integer :: id_chct=-1, id_chcb=-1, id_chcc=-1, id_chch=-1
44 integer :: id_t0=-1, id_tf=-1, id_s0=-1, id_sf=-1, id_n2_0=-1, id_n2_f=-1
45 integer :: id_h=-1, id_zint=-1
46 !>@}
48
49contains
50
51!> This subroutine helps test the accuracy of the diapycnal mixing energy requirement code
52!! by writing diagnostics, possibly using an intensely mixing test profile of diffusivity
53subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int)
54 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
55 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
56 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
57 real, dimension(G%isd:G%ied,G%jsd:G%jed,GV%ke), &
58 intent(in) :: h_3d !< Layer thickness before entrainment [H ~> m or kg m-2].
59 type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
60 !! available thermodynamic fields.
61 !! Absent fields have NULL ptrs.
62 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s].
63 type(diapyc_energy_req_cs), pointer :: cs !< This module's control structure.
64 real, dimension(G%isd:G%ied,G%jsd:G%jed,GV%ke+1), &
65 optional, intent(in) :: kd_int !< Interface diffusivities [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
66
67 ! Local variables
68 real, dimension(GV%ke) :: &
69 t0, s0, & ! T0 & S0 are columns of initial temperatures and salinities [C ~> degC] and [S ~> ppt].
70 h_col, & ! h_col is a column of thicknesses h at tracer points [H ~> m or kg m-2].
71 dz_col ! dz_col is a column of vertical distances across layers at tracer points [Z ~> m]
72 real, dimension( G%isd:G%ied,GV%ke) :: &
73 dz_2d ! A 2-d slice of the vertical distance across layers [Z ~> m]
74 real, dimension(GV%ke+1) :: &
75 kd, & ! A column of diapycnal diffusivities at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
76 h_top, h_bot ! Distances from the top or bottom [H ~> m or kg m-2].
77 real :: dz_h_int ! The ratio of the vertical distances across the layers surrounding an interface
78 ! over the layer thicknesses [H Z-1 ~> nondim or kg m-3]
79 real :: ustar ! The local friction velocity [Z T-1 ~> m s-1]
80 real :: absf ! The absolute value of the Coriolis parameter [T-1 ~> s-1]
81 real :: htot ! The sum of the thicknesses [H ~> m or kg m-2].
82 real :: energy_kd ! The energy used by diapycnal mixing [R Z L2 T-3 ~> W m-2].
83 real :: tmp1 ! A temporary array [H2 ~> m2 or kg2 m-4]
84 integer :: i, j, k, is, ie, js, je, nz
85 logical :: may_print
86 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
87
88 if (.not. associated(cs)) call mom_error(fatal, "diapyc_energy_req_test: "// &
89 "Module must be initialized before it is used.")
90
91 if (.not. cs%initialized) call mom_error(fatal, "diapyc_energy_req_test: "// &
92 "Module must be initialized before it is used.")
93
94!$OMP do
95 do j=js,je
96 call thickness_to_dz(h_3d, tv, dz_2d, j, g, gv)
97
98 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
99
100 do k=1,nz
101 t0(k) = tv%T(i,j,k) ; s0(k) = tv%S(i,j,k)
102 h_col(k) = h_3d(i,j,k)
103 dz_col(k) = dz_2d(i,k)
104 enddo
105
106 if (present(kd_int) .and. .not.cs%use_test_Kh_profile) then
107 do k=1,nz+1 ; kd(k) = cs%test_Kh_scaling*kd_int(i,j,k) ; enddo
108 else
109 htot = 0.0 ; h_top(1) = 0.0
110 do k=1,nz
111 h_top(k+1) = h_top(k) + h_col(k)
112 enddo
113 htot = h_top(nz+1)
114
115 h_bot(nz+1) = 0.0
116 do k=nz,1,-1
117 h_bot(k) = h_bot(k+1) + h_col(k)
118 enddo
119
120 ustar = 0.01*us%m_to_Z*us%T_to_s ! Change this to being an input parameter?
121 absf = 0.25*((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
122 (abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))))
123 kd(1) = 0.0 ; kd(nz+1) = 0.0
124 if (gv%Boussinesq) then
125 do k=2,nz
126 tmp1 = h_top(k) * h_bot(k)
127 kd(k) = cs%test_Kh_scaling * &
128 ustar * cs%VonKar * (tmp1*ustar) / (absf*gv%H_to_Z*tmp1 + htot*ustar)
129 enddo
130 else
131 do k=2,nz
132 tmp1 = h_top(k) * h_bot(k)
133 dz_h_int = (dz_2d(j,k-1) + dz_2d(j,k) + gv%dz_subroundoff) / &
134 (h_3d(i,j,k-1) + h_3d(i,j,k) + gv%H_subroundoff)
135 kd(k) = cs%test_Kh_scaling * &
136 ustar * cs%VonKar * (tmp1*ustar) / (dz_h_int*absf*tmp1 + htot*ustar)
137 enddo
138 endif
139 endif
140 may_print = is_root_pe() .and. (i==ie) .and. (j==je)
141 call diapyc_energy_req_calc(h_col, dz_col, t0, s0, kd, energy_kd, dt, tv, g, gv, us, &
142 may_print=may_print, cs=cs)
143 endif ; enddo
144 enddo
145
146end subroutine diapyc_energy_req_test
147
148!> This subroutine uses a substantially refactored tridiagonal equation for
149!! diapycnal mixing of temperature and salinity to estimate the potential energy
150!! change due to diapycnal mixing in a column of water. It does this estimate
151!! 4 different ways, all of which should be equivalent, but reports only one.
152!! The various estimates are taken because they will later be used as templates
153!! for other bits of code
154subroutine diapyc_energy_req_calc(h_in, dz_in, T_in, S_in, Kd, energy_Kd, dt, tv, &
155 G, GV, US, may_print, CS)
156 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
157 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
158 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
159 real, dimension(GV%ke), intent(in) :: h_in !< Layer thickness before entrainment,
160 !! [H ~> m or kg m-2]
161 real, dimension(GV%ke), intent(in) :: dz_in !< Vertical distance across layers before
162 !! entrainment [Z ~> m]
163 real, dimension(GV%ke), intent(in) :: t_in !< The layer temperatures [C ~> degC].
164 real, dimension(GV%ke), intent(in) :: s_in !< The layer salinities [S ~> ppt].
165 real, dimension(GV%ke+1), intent(in) :: kd !< The interfaces diapycnal diffusivities
166 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
167 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s].
168 real, intent(out) :: energy_kd !< The column-integrated rate of energy
169 !! consumption by diapycnal diffusion [R Z L2 T-3 ~> W m-2].
170 type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
171 !! available thermodynamic fields.
172 !! Absent fields have NULL ptrs.
173 logical, optional, intent(in) :: may_print !< If present and true, write out diagnostics
174 !! of energy use.
175 type(diapyc_energy_req_cs), &
176 optional, pointer :: cs !< This module's control structure.
177
178! This subroutine uses a substantially refactored tridiagonal equation for
179! diapycnal mixing of temperature and salinity to estimate the potential energy
180! change due to diapycnal mixing in a column of water. It does this estimate
181! 4 different ways, all of which should be equivalent, but reports only one.
182! The various estimates are taken because they will later be used as templates
183! for other bits of code.
184
185 real, dimension(GV%ke) :: &
186 p_lay, & ! Average pressure of a layer [R L2 T-2 ~> Pa].
187 dsv_dt, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1].
188 dsv_ds, & ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
189 t0, s0, & ! Initial temperatures and salinities [C ~> degC] and [S ~> ppt].
190 tf, sf, & ! New final values of the temperatures and salinities [C ~> degC] and [S ~> ppt].
191 th_a, & ! An effective temperature times a thickness in the layer above, including implicit
192 ! mixing effects with other yet higher layers [C H ~> degC m or degC kg m-2].
193 sh_a, & ! An effective salinity times a thickness in the layer above, including implicit
194 ! mixing effects with other yet higher layers [S H ~> ppt m or ppt kg m-2].
195 th_b, & ! An effective temperature times a thickness in the layer below, including implicit
196 ! mixing effects with other yet lower layers [C H ~> degC m or degC kg m-2].
197 sh_b, & ! An effective salinity times a thickness in the layer below, including implicit
198 ! mixing effects with other yet lower layers [S H ~> ppt m or ppt kg m-2].
199 dt_to_dpe, & ! Partial derivative of column potential energy with the temperature changes within
200 ! a layer [R Z L2 T-2 C-1 ~> J m-2 degC-1]
201 ds_to_dpe, & ! Partial derivative of column potential energy with the salinity changes within
202 ! a layer [R Z L2 T-2 S-1 ~> J m-2 ppt-1]
203 dt_to_dcolht, & ! Partial derivative of the total column height with the temperature
204 ! changes within a layer [Z C-1 ~> m degC-1]
205 ds_to_dcolht, & ! Partial derivative of the total column height with the
206 ! salinity changes within a layer [Z S-1 ~> m ppt-1].
207 dt_to_dcolht_a, & ! Partial derivative of the total column height with the temperature changes
208 ! within a layer, including the implicit effects of mixing with layers higher
209 ! in the water column [Z C-1 ~> m degC-1].
210 ds_to_dcolht_a, & ! Partial derivative of the total column height with the salinity changes
211 ! within a layer, including the implicit effects of mixing with layers higher
212 ! of mixing with layers higher in the water column [Z S-1 ~> m ppt-1].
213 dt_to_dcolht_b, & ! Partial derivative of the total column height with the temperature changes
214 ! within a layer, including the implicit effects of mixing with layers lower
215 ! in the water column [Z C-1 ~> m degC-1].
216 ds_to_dcolht_b, & ! Partial derivative of the total column height with the salinity changes
217 ! within a layer, including the implicit effects of mixing with layers lower
218 ! in the water column [Z S-1 ~> m ppt-1].
219 dt_to_dpe_a, & ! Partial derivative of column potential energy with the temperature changes
220 ! within a layer, including the implicit effects of mixing with layers higher
221 ! in the water column, in units of [R Z L2 T-2 C-1 ~> J m-2 degC-1].
222 ds_to_dpe_a, & ! Partial derivative of column potential energy with the salinity changes
223 ! within a layer, including the implicit effects of mixing with layers higher
224 ! in the water column, in units of [R Z L2 T-2 S-1 ~> J m-2 ppt-1].
225 dt_to_dpe_b, & ! Partial derivative of column potential energy with the temperature changes
226 ! within a layer, including the implicit effects of mixing with layers lower
227 ! in the water column, in units of [R Z L2 T-2 C-1 ~> J m-2 degC-1].
228 ds_to_dpe_b, & ! Partial derivative of column potential energy with the salinity changes
229 ! within a layer, including the implicit effects of mixing with layers lower
230 ! in the water column, in units of [R Z L2 T-2 S-1 ~> J m-2 ppt-1].
231 hp_a, & ! An effective pivot thickness of the layer including the effects
232 ! of coupling with layers above [H ~> m or kg m-2]. This is the first term
233 ! in the denominator of b1 in a downward-oriented tridiagonal solver.
234 hp_b, & ! An effective pivot thickness of the layer including the effects
235 ! of coupling with layers below [H ~> m or kg m-2]. This is the first term
236 ! in the denominator of b1 in an upward-oriented tridiagonal solver.
237 c1_a, & ! c1_a is used by a downward-oriented tridiagonal solver [nondim].
238 c1_b, & ! c1_b is used by an upward-oriented tridiagonal solver [nondim].
239 h_tr, & ! h_tr is h at tracer points with a h_neglect added to
240 ! ensure positive definiteness [H ~> m or kg m-2].
241 dz_tr ! dz_tr is dz at tracer points with dz_neglect added to
242 ! ensure positive definiteness [Z ~> m]
243 ! Note that the following arrays have extra (ficticious) layers above or below the
244 ! water column for code convenience
245 real, dimension(0:GV%ke+1) :: &
246 te, se ! Running incomplete estimates of the new temperatures and salinities [C ~> degC] and [S ~> ppt]
247 real, dimension(0:GV%ke) :: &
248 te_a, se_a ! Running incomplete estimates of the new temperatures and salinities [C ~> degC] and [S ~> ppt]
249 real, dimension(GV%ke+1) :: &
250 te_b, se_b ! Running incomplete estimates of the new temperatures and salinities [C ~> degC] and [S ~> ppt]
251 real, dimension(GV%ke+1) :: &
252 pres, & ! Interface pressures [R L2 T-2 ~> Pa].
253 pres_z, & ! The hydrostatic interface pressure, which is used to relate
254 ! the changes in column thickness to the energy that is radiated
255 ! as gravity waves and unavailable to drive mixing [R L2 T-2 ~> J m-3].
256 z_int, & ! Interface heights relative to the surface [H ~> m or kg m-2].
257 n2, & ! An estimate of the buoyancy frequency [T-2 ~> s-2].
258 kddt_h, & ! The diapycnal diffusivity times a timestep divided by the
259 ! average thicknesses around a layer [H ~> m or kg m-2].
260 kddt_h_a, & ! The value of Kddt_h for layers above the central point in the
261 ! tridiagonal solver [H ~> m or kg m-2].
262 kddt_h_b, & ! The value of Kddt_h for layers below the central point in the
263 ! tridiagonal solver [H ~> m or kg m-2].
264 kd_so_far ! The value of Kddt_h that has been applied already in
265 ! calculating the energy changes [H ~> m or kg m-2].
266 real, dimension(GV%ke+1,4) :: &
267 pe_chg_k, & ! The integrated potential energy change within a timestep due
268 ! to the diffusivity at interface K for 4 different orders of
269 ! accumulating the diffusivities [R Z L2 T-2 ~> J m-2].
270 colht_cor_k ! The correction to the potential energy change due to
271 ! changes in the net column height [R Z L2 T-2 ~> J m-2].
272 real :: b1 ! b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
273 real :: kd0 ! The value of Kddt_h that has already been applied [H ~> m or kg m-2].
274 real :: dkd ! The change in the value of Kddt_h [H ~> m or kg m-2].
275 real :: h_neglect ! A thickness that is so small it is usually lost
276 ! in roundoff and can be neglected [H ~> m or kg m-2].
277 real :: kddt_h_guess ! A guess of the final value of Kddt_h [H ~> m or kg m-2].
278 real :: dmass ! The mass per unit area within a layer [R Z ~> kg m-2].
279 real :: dpres ! The hydrostatic pressure change across a layer [R L2 T-2 ~> Pa].
280 real :: rho_here ! The in-situ density [R ~> kg m-3].
281 real :: pe_change ! The change in column potential energy from applying Kddt_h at the
282 ! present interface [R L2 Z T-2 ~> J m-2].
283 real :: colht_cor ! The correction to PE_chg that is made due to a net
284 ! change in the column height [R L2 Z T-2 ~> J m-2].
285 real :: htot ! A running sum of thicknesses [H ~> m or kg m-2].
286 real :: dztot ! A running sum of vertical distances across layers [Z ~> m]
287 logical :: do_print
288
289 ! The following are a bunch of diagnostic arrays for debugging purposes.
290 real, dimension(GV%ke) :: &
291 ta, tb, & ! Copies of temperature profiles for debugging [C ~> degC]
292 sa, sb ! Copies of salinity profiles for debugging [S ~> ppt]
293 real, dimension(GV%ke+1) :: &
294 dpea_dkd, dpea_dkd_est, & ! Estimates of the partial derivative of the column potential energy
295 ! change with Kddt_h [R Z L2 T-2 H-1 ~> J m-3 or J kg-1].
296 dpeb_dkd, dpeb_dkd_est, & ! Estimates of the partial derivative of the column potential energy
297 ! change with Kddt_h [R Z L2 T-2 H-1 ~> J m-3 or J kg-1].
298 dpea_dkd_err, dpeb_dkd_err, & ! Differences in estimates of the partial derivative of the column
299 ! potential energy change with Kddt_h [R Z L2 T-2 H-1 ~> J m-3 or J kg-1].
300 dpea_dkd_err_norm, dpeb_dkd_err_norm, & ! Normalized changes in sensitivities [nondim]
301 dpea_dkd_trunc, dpeb_dkd_trunc ! Estimates of the truncation error in estimates of the partial
302 ! derivative of the column potential energy change with
303 ! Kddt_h [R Z L2 T-2 H-1 ~> J m-3 or J kg-1].
304 real :: pe_chg_tot1a, pe_chg_tot2a ! Changes in column potential energy [R Z L2 T-2 ~> J m-2]
305 real :: pe_chg_tot1b, pe_chg_tot2b ! Changes in column potential energy [R Z L2 T-2 ~> J m-2]
306 real :: pe_chg_tot1c, pe_chg_tot2c ! Changes in column potential energy [R Z L2 T-2 ~> J m-2]
307 real :: pe_chg_tot1d, pe_chg_tot2d ! Changes in column potential energy [R Z L2 T-2 ~> J m-2]
308 real :: t_chg_tota, t_chg_totb ! Vertically integrated temperature changes [C H ~> degC m or degC kg m-2]
309 real :: t_chg_totc, t_chg_totd ! Vertically integrated temperature changes [C H ~> degC m or degC kg m-2]
310 real :: pe_chg(6) ! The potential energy change within the first few iterations [R Z L2 T-2 ~> J m-2]
311
312 integer :: k, nz, itt, k_cent
313 logical :: surface_bl, bottom_bl, central, halves, debug
314 nz = gv%ke
315 h_neglect = gv%H_subroundoff
316
317 debug = .true.
318
319 surface_bl = .true. ; bottom_bl = .true. ; halves = .true.
320 central = .true. ; k_cent = nz/2
321
322 do_print = .false. ; if (present(may_print) .and. present(cs)) do_print = may_print
323
324 dpea_dkd(:) = 0.0 ; dpea_dkd_est(:) = 0.0 ; dpea_dkd_err(:) = 0.0
325 dpea_dkd_err_norm(:) = 0.0 ; dpea_dkd_trunc(:) = 0.0
326 dpeb_dkd(:) = 0.0 ; dpeb_dkd_est(:) = 0.0 ; dpeb_dkd_err(:) = 0.0
327 dpeb_dkd_err_norm(:) = 0.0 ; dpeb_dkd_trunc(:) = 0.0
328
329 htot = 0.0 ; dztot = 0.0 ; pres(1) = 0.0 ; pres_z(1) = 0.0 ; z_int(1) = 0.0
330 do k=1,nz
331 t0(k) = t_in(k) ; s0(k) = s_in(k)
332 h_tr(k) = h_in(k)
333 dz_tr(k) = dz_in(k)
334 htot = htot + h_tr(k)
335 dztot = dztot + dz_tr(k)
336 pres(k+1) = pres(k) + (gv%g_Earth * gv%H_to_RZ) * h_tr(k)
337 pres_z(k+1) = pres(k+1)
338 p_lay(k) = 0.5*(pres(k) + pres(k+1))
339 z_int(k+1) = z_int(k) - h_tr(k)
340 enddo
341 do k=1,nz
342 h_tr(k) = max(h_tr(k), 1e-15*htot)
343 dz_tr(k) = max(dz_tr(k), 1e-15*dztot)
344 enddo
345
346 ! Introduce a diffusive flux variable, Kddt_h(K) = ea(k) = eb(k-1)
347
348 kddt_h(1) = 0.0 ; kddt_h(nz+1) = 0.0
349 do k=2,nz
350 kddt_h(k) = min(dt * kd(k) / (0.5*(dz_tr(k-1) + dz_tr(k))), 1e3*dztot)
351 enddo
352
353 ! Zero out the temperature and salinity estimates in the extra (ficticious) layers.
354 ! The actual values set here are irrelevant (so long as they are not NaNs) because they
355 ! are always multiplied by a zero value of Kddt_h reflecting the no-flux boundary condition.
356 te(0) = 0.0 ; se(0) = 0.0 ; te(nz+1) = 0.0 ; se(nz+1) = 0.0
357 te_a(0) = 0.0 ; se_a(0) = 0.0
358 te_b(nz+1) = 0.0 ; se_b(nz+1) = 0.0
359
360 ! Solve the tridiagonal equations for new temperatures.
361
362 call calculate_specific_vol_derivs(t0, s0, p_lay, dsv_dt, dsv_ds, tv%eqn_of_state)
363
364 do k=1,nz
365 dmass = gv%H_to_RZ * h_tr(k)
366 dpres = (gv%g_Earth * gv%H_to_RZ) * h_tr(k)
367 dt_to_dpe(k) = (dmass * (pres(k) + 0.5*dpres)) * dsv_dt(k)
368 ds_to_dpe(k) = (dmass * (pres(k) + 0.5*dpres)) * dsv_ds(k)
369 dt_to_dcolht(k) = dmass * dsv_dt(k) * cs%ColHt_scaling
370 ds_to_dcolht(k) = dmass * dsv_ds(k) * cs%ColHt_scaling
371 enddo
372
373! PE_chg_k(1) = 0.0 ; PE_chg_k(nz+1) = 0.0
374 ! PEchg(:) = 0.0
375 pe_chg_k(:,:) = 0.0 ; colht_cor_k(:,:) = 0.0
376
377 if (surface_bl) then ! This version is appropriate for a surface boundary layer.
378
379 ! Set up values appropriate for no diffusivity.
380 do k=1,nz
381 hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
382 dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
383 dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
384 dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
385 dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
386 enddo
387
388 do k=2,nz
389 ! At this point, Kddt_h(K) will be unknown because its value may depend
390 ! on how much energy is available.
391
392 ! Precalculate some temporary expressions that are independent of Kddt_h_guess.
393 th_a(k-1) = h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2)
394 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2)
395 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
396
397
398 ! Find the energy change due to a guess at the strength of diffusion at interface K.
399
400 kddt_h_guess = kddt_h(k)
401 call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
402 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
403 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
404 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
405 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
406 pe_chg=pe_chg_k(k,1), dpec_dkd=dpea_dkd(k), &
407 pe_colht_cor=colht_cor_k(k,1))
408
409 if (debug) then
410 do itt=1,5
411 kddt_h_guess = (1.0+0.01*(itt-3))*kddt_h(k)
412
413 call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
414 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
415 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
416 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
417 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
418 pe_chg=pe_chg(itt))
419 enddo
420 ! Compare with a 4th-order finite difference estimate.
421 dpea_dkd_est(k) = (4.0*(pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
422 (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))) / 3.0
423 dpea_dkd_trunc(k) = (pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
424 (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))
425 dpea_dkd_err(k) = (dpea_dkd_est(k) - dpea_dkd(k))
426 dpea_dkd_err_norm(k) = (dpea_dkd_est(k) - dpea_dkd(k)) / &
427 (abs(dpea_dkd_est(k)) + abs(dpea_dkd(k)) + 1e-100*us%RZ_to_kg_m2*us%L_T_to_m_s**2)
428 endif
429
430 ! At this point, the final value of Kddt_h(K) is known, so the estimated
431 ! properties for layer k-1 can be calculated.
432
433 b1 = 1.0 / (hp_a(k-1) + kddt_h(k))
434 c1_a(k) = kddt_h(k) * b1
435 te(k-1) = b1 * (h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2))
436 se(k-1) = b1 * (h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2))
437
438 hp_a(k) = h_tr(k) + (hp_a(k-1) * b1) * kddt_h(k)
439 dt_to_dpe_a(k) = dt_to_dpe(k) + c1_a(k)*dt_to_dpe_a(k-1)
440 ds_to_dpe_a(k) = ds_to_dpe(k) + c1_a(k)*ds_to_dpe_a(k-1)
441 dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1_a(k)*dt_to_dcolht_a(k-1)
442 ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1_a(k)*ds_to_dcolht_a(k-1)
443
444 enddo
445
446 b1 = 1.0 / (hp_a(nz))
447 tf(nz) = b1 * (h_tr(nz) * t0(nz) + kddt_h(nz) * te(nz-1))
448 sf(nz) = b1 * (h_tr(nz) * s0(nz) + kddt_h(nz) * se(nz-1))
449
450 do k=nz-1,1,-1
451 tf(k) = te(k) + c1_a(k+1)*tf(k+1)
452 sf(k) = se(k) + c1_a(k+1)*sf(k+1)
453 enddo
454
455 if (debug) then
456 do k=1,nz ; ta(k) = tf(k) ; sa(k) = sf(k) ; enddo
457 pe_chg_tot1a = 0.0 ; pe_chg_tot2a = 0.0 ; t_chg_tota = 0.0
458 do k=1,nz
459 pe_chg_tot1a = pe_chg_tot1a + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
460 ds_to_dpe(k) * (sf(k) - s0(k)))
461 t_chg_tota = t_chg_tota + h_tr(k) * (tf(k) - t0(k))
462 enddo
463 do k=2,nz
464 pe_chg_tot2a = pe_chg_tot2a + (pe_chg_k(k,1) - colht_cor_k(k,1))
465 enddo
466 endif
467
468 endif
469
470 if (bottom_bl) then ! This version is appropriate for a bottom boundary layer.
471
472 ! Set up values appropriate for no diffusivity.
473 do k=1,nz
474 hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
475 dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
476 dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
477 dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
478 dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
479 enddo
480
481 do k=nz,2,-1 ! Loop over interior interfaces.
482 ! At this point, Kddt_h(K) will be unknown because its value may depend
483 ! on how much energy is available.
484
485 ! Precalculate some temporary expressions that are independent of Kddt_h_guess.
486 th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
487 th_b(k) = h_tr(k) * t0(k) + kddt_h(k+1) * te(k+1)
488 sh_b(k) = h_tr(k) * s0(k) + kddt_h(k+1) * se(k+1)
489
490 ! Find the energy change due to a guess at the strength of diffusion at interface K.
491 kddt_h_guess = kddt_h(k)
492
493 call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
494 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
495 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
496 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
497 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
498 pe_chg=pe_chg_k(k,2), dpec_dkd=dpeb_dkd(k), &
499 pe_colht_cor=colht_cor_k(k,2))
500
501 if (debug) then
502 ! Compare with a 4th-order finite difference estimate.
503 do itt=1,5
504 kddt_h_guess = (1.0+0.01*(itt-3))*kddt_h(k)
505
506 call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
507 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
508 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
509 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
510 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
511 pe_chg=pe_chg(itt))
512 enddo
513
514 dpeb_dkd_est(k) = (4.0*(pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
515 (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))) / 3.0
516 dpeb_dkd_trunc(k) = (pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
517 (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))
518 dpeb_dkd_err(k) = (dpeb_dkd_est(k) - dpeb_dkd(k))
519 dpeb_dkd_err_norm(k) = (dpeb_dkd_est(k) - dpeb_dkd(k)) / &
520 (abs(dpeb_dkd_est(k)) + abs(dpeb_dkd(k)) + 1e-100*us%RZ_to_kg_m2*us%L_T_to_m_s**2)
521 endif
522
523 ! At this point, the final value of Kddt_h(K) is known, so the estimated
524 ! properties for layer k can be calculated.
525
526 b1 = 1.0 / (hp_b(k) + kddt_h(k))
527 c1_b(k) = kddt_h(k) * b1
528
529 te(k) = b1 * (h_tr(k) * t0(k) + kddt_h(k+1) * te(k+1))
530 se(k) = b1 * (h_tr(k) * s0(k) + kddt_h(k+1) * se(k+1))
531
532 hp_b(k-1) = h_tr(k-1) + (hp_b(k) * b1) * kddt_h(k)
533 dt_to_dpe_b(k-1) = dt_to_dpe(k-1) + c1_b(k)*dt_to_dpe_b(k)
534 ds_to_dpe_b(k-1) = ds_to_dpe(k-1) + c1_b(k)*ds_to_dpe_b(k)
535 dt_to_dcolht_b(k-1) = dt_to_dcolht(k-1) + c1_b(k)*dt_to_dcolht_b(k)
536 ds_to_dcolht_b(k-1) = ds_to_dcolht(k-1) + c1_b(k)*ds_to_dcolht_b(k)
537
538 enddo
539
540 b1 = 1.0 / (hp_b(1))
541 tf(1) = b1 * (h_tr(1) * t0(1) + kddt_h(2) * te(2))
542 sf(1) = b1 * (h_tr(1) * s0(1) + kddt_h(2) * se(2))
543
544 do k=2,nz
545 tf(k) = te(k) + c1_b(k)*tf(k-1)
546 sf(k) = se(k) + c1_b(k)*sf(k-1)
547 enddo
548
549 if (debug) then
550 do k=1,nz ; tb(k) = tf(k) ; sb(k) = sf(k) ; enddo
551 pe_chg_tot1b = 0.0 ; pe_chg_tot2b = 0.0 ; t_chg_totb = 0.0
552 do k=1,nz
553 pe_chg_tot1b = pe_chg_tot1b + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
554 ds_to_dpe(k) * (sf(k) - s0(k)))
555 t_chg_totb = t_chg_totb + h_tr(k) * (tf(k) - t0(k))
556 enddo
557 do k=2,nz
558 pe_chg_tot2b = pe_chg_tot2b + (pe_chg_k(k,2) - colht_cor_k(k,2))
559 enddo
560 endif
561
562 endif
563
564 if (central) then
565
566 ! Set up values appropriate for no diffusivity.
567 do k=1,nz
568 hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
569 dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
570 dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
571 dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
572 dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
573 enddo
574
575 ! Calculate the dependencies on layers above.
576 kddt_h_a(1) = 0.0
577 do k=2,nz ! Loop over interior interfaces.
578 ! First calculate some terms that are independent of the change in Kddt_h(K).
579 kd0 = 0.0 ! This might need to be changed - it is the already applied value of Kddt_h(K).
580
581 th_a(k-1) = h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te_a(k-2)
582 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se_a(k-2)
583 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
584
585 kddt_h_a(k) = 0.0 ; if (k < k_cent) kddt_h_a(k) = kddt_h(k)
586 dkd = kddt_h_a(k)
587
588 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
589 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
590 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
591 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
592 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
593 pe_chg=pe_change, pe_colht_cor=colht_cor)
594 pe_chg_k(k,3) = pe_change
595 colht_cor_k(k,3) = colht_cor
596
597 b1 = 1.0 / (hp_a(k-1) + kddt_h_a(k))
598 c1_a(k) = kddt_h_a(k) * b1
599
600 te_a(k-1) = b1 * (h_tr(k-1) * t0(k-1) + kddt_h_a(k-1) * te_a(k-2))
601 se_a(k-1) = b1 * (h_tr(k-1) * s0(k-1) + kddt_h_a(k-1) * se_a(k-2))
602
603 hp_a(k) = h_tr(k) + (hp_a(k-1) * b1) * kddt_h_a(k)
604 dt_to_dpe_a(k) = dt_to_dpe(k) + c1_a(k)*dt_to_dpe_a(k-1)
605 ds_to_dpe_a(k) = ds_to_dpe(k) + c1_a(k)*ds_to_dpe_a(k-1)
606 dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1_a(k)*dt_to_dcolht_a(k-1)
607 ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1_a(k)*ds_to_dcolht_a(k-1)
608 enddo
609
610 ! Calculate the dependencies on layers below.
611 kddt_h_b(nz+1) = 0.0
612 do k=nz,2,-1 ! Loop over interior interfaces.
613 ! First calculate some terms that are independent of the change in Kddt_h(K).
614 kd0 = 0.0 ! This might need to be changed - it is the already applied value of Kddt_h(K).
615
616 th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
617! Th_a(k-1) = h_tr(k-1) * T0(k-1) + Kddt_h(K-1) * Te_a(k-2)
618! Sh_a(k-1) = h_tr(k-1) * S0(k-1) + Kddt_h(K-1) * Se_a(k-2)
619
620 th_b(k) = h_tr(k) * t0(k) + kddt_h(k+1) * te_b(k+1)
621 sh_b(k) = h_tr(k) * s0(k) + kddt_h(k+1) * se_b(k+1)
622
623 kddt_h_b(k) = 0.0 ; if (k > k_cent) kddt_h_b(k) = kddt_h(k)
624 dkd = kddt_h_b(k)
625
626 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
627 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
628 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
629 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
630 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
631 pe_chg=pe_change, pe_colht_cor=colht_cor)
632 pe_chg_k(k,3) = pe_chg_k(k,3) + pe_change
633 colht_cor_k(k,3) = colht_cor_k(k,3) + colht_cor
634
635 b1 = 1.0 / (hp_b(k) + kddt_h_b(k))
636 c1_b(k) = kddt_h_b(k) * b1
637
638 te_b(k) = b1 * (h_tr(k) * t0(k) + kddt_h_b(k+1) * te_b(k+1))
639 se_b(k) = b1 * (h_tr(k) * s0(k) + kddt_h_b(k+1) * se_b(k+1))
640
641 hp_b(k-1) = h_tr(k-1) + (hp_b(k) * b1) * kddt_h_b(k)
642 dt_to_dpe_b(k-1) = dt_to_dpe(k-1) + c1_b(k)*dt_to_dpe_b(k)
643 ds_to_dpe_b(k-1) = ds_to_dpe(k-1) + c1_b(k)*ds_to_dpe_b(k)
644 dt_to_dcolht_b(k-1) = dt_to_dcolht(k-1) + c1_b(k)*dt_to_dcolht_b(k)
645 ds_to_dcolht_b(k-1) = ds_to_dcolht(k-1) + c1_b(k)*ds_to_dcolht_b(k)
646
647 enddo
648
649 ! Calculate the final solution for the layers surrounding interface K_cent
650 k = k_cent
651
652 ! First calculate some terms that are independent of the change in Kddt_h(K).
653 kd0 = 0.0 ! This might need to be changed - it is the already applied value of Kddt_h(K).
654
655 th_a(k-1) = h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te_a(k-2)
656 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se_a(k-2)
657 th_b(k) = h_tr(k) * t0(k) + kddt_h(k+1) * te_b(k+1)
658 sh_b(k) = h_tr(k) * s0(k) + kddt_h(k+1) * se_b(k+1)
659
660 dkd = kddt_h(k)
661
662 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
663 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
664 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
665 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
666 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
667 pe_chg=pe_change, pe_colht_cor=colht_cor)
668 pe_chg_k(k,3) = pe_chg_k(k,3) + pe_change
669 colht_cor_k(k,3) = colht_cor_k(k,3) + colht_cor
670
671
672 ! To derive the following, first to a partial update for the estimated
673 ! temperatures and salinities in the layers around this interface:
674 ! b1_a = 1.0 / (hp_a(k-1) + Kddt_h(K))
675 ! b1_b = 1.0 / (hp_b(k) + Kddt_h(K))
676 ! Te_up(k) = Th_b(k) * b1_b ; Se_up(k) = Sh_b(k) * b1_b
677 ! Te_up(k-1) = Th_a(k-1) * b1_a ; Se_up(k-1) = Sh_a(k-1) * b1_a
678 ! Find the final values of T & S for both layers around K_cent, using that
679 ! c1_a(K) = Kddt_h(K) * b1_a ; c1_b(K) = Kddt_h(K) * b1_b
680 ! Tf(K_cent-1) = Te_up(K_cent-1) + c1_a(K_cent)*Tf(K_cent)
681 ! Tf(K_cent) = Te_up(K_cent) + c1_b(K_cent)*Tf(K_cent-1)
682 ! but further exploiting the expressions for c1_a and c1_b to avoid
683 ! subtraction in the denominator, and use only a single division.
684 b1 = 1.0 / (hp_a(k-1)*hp_b(k) + kddt_h(k)*(hp_a(k-1) + hp_b(k)))
685 tf(k-1) = ((hp_b(k) + kddt_h(k)) * th_a(k-1) + kddt_h(k) * th_b(k) ) * b1
686 sf(k-1) = ((hp_b(k) + kddt_h(k)) * sh_a(k-1) + kddt_h(k) * sh_b(k) ) * b1
687 tf(k) = (kddt_h(k) * th_a(k-1) + (hp_a(k-1) + kddt_h(k)) * th_b(k) ) * b1
688 sf(k) = (kddt_h(k) * sh_a(k-1) + (hp_a(k-1) + kddt_h(k)) * sh_b(k) ) * b1
689
690 c1_a(k) = kddt_h(k) / (hp_a(k-1) + kddt_h(k))
691 c1_b(k) = kddt_h(k) / (hp_b(k) + kddt_h(k))
692
693 ! Now update the other layer working outward from k_cent to determine the final
694 ! temperatures and salinities.
695 do k=k_cent-2,1,-1
696 tf(k) = te_a(k) + c1_a(k+1)*tf(k+1)
697 sf(k) = se_a(k) + c1_a(k+1)*sf(k+1)
698 enddo
699 do k=k_cent+1,nz
700 tf(k) = te_b(k) + c1_b(k)*tf(k-1)
701 sf(k) = se_b(k) + c1_b(k)*sf(k-1)
702 enddo
703
704 if (debug) then
705 pe_chg_tot1c = 0.0 ; pe_chg_tot2c = 0.0 ; t_chg_totc = 0.0
706 do k=1,nz
707 pe_chg_tot1c = pe_chg_tot1c + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
708 ds_to_dpe(k) * (sf(k) - s0(k)))
709 t_chg_totc = t_chg_totc + h_tr(k) * (tf(k) - t0(k))
710 enddo
711 do k=2,nz
712 pe_chg_tot2c = pe_chg_tot2c + (pe_chg_k(k,3) - colht_cor_k(k,3))
713 enddo
714 endif
715
716 endif
717
718 if (halves) then
719
720 ! Set up values appropriate for no diffusivity.
721 do k=1,nz
722 hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
723 dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
724 dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
725 dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
726 dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
727 enddo
728 do k=1,nz+1
729 kd_so_far(k) = 0.0
730 enddo
731
732 ! Calculate the dependencies on layers above.
733 do k=2,nz ! Loop over interior interfaces.
734 ! First calculate some terms that are independent of the change in Kddt_h(K).
735 kd0 = kd_so_far(k)
736
737 th_a(k-1) = h_tr(k-1) * t0(k-1) + kd_so_far(k-1) * te(k-2)
738 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kd_so_far(k-1) * se(k-2)
739 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
740
741 dkd = 0.5 * kddt_h(k) - kd_so_far(k)
742
743 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
744 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
745 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
746 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
747 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
748 pe_chg=pe_change, pe_colht_cor=colht_cor)
749
750 pe_chg_k(k,4) = pe_change
751 colht_cor_k(k,4) = colht_cor
752
753 kd_so_far(k) = kd_so_far(k) + dkd
754
755 b1 = 1.0 / (hp_a(k-1) + kd_so_far(k))
756 c1_a(k) = kd_so_far(k) * b1
757
758 te(k-1) = b1 * (h_tr(k-1) * t0(k-1) + kd_so_far(k-1) * te(k-2))
759 se(k-1) = b1 * (h_tr(k-1) * s0(k-1) + kd_so_far(k-1) * se(k-2))
760
761 hp_a(k) = h_tr(k) + (hp_a(k-1) * b1) * kd_so_far(k)
762 dt_to_dpe_a(k) = dt_to_dpe(k) + c1_a(k)*dt_to_dpe_a(k-1)
763 ds_to_dpe_a(k) = ds_to_dpe(k) + c1_a(k)*ds_to_dpe_a(k-1)
764 dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1_a(k)*dt_to_dcolht_a(k-1)
765 ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1_a(k)*ds_to_dcolht_a(k-1)
766 enddo
767
768 ! Calculate the dependencies on layers below.
769 do k=nz,2,-1 ! Loop over interior interfaces.
770 ! First calculate some terms that are independent of the change in Kddt_h(K).
771 kd0 = kd_so_far(k)
772
773 th_a(k-1) = h_tr(k-1) * t0(k-1) + kd_so_far(k-1) * te(k-2)
774 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kd_so_far(k-1) * se(k-2)
775 th_b(k) = h_tr(k) * t0(k) + kd_so_far(k+1) * te(k+1)
776 sh_b(k) = h_tr(k) * s0(k) + kd_so_far(k+1) * se(k+1)
777
778 dkd = kddt_h(k) - kd_so_far(k)
779
780 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
781 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
782 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
783 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
784 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
785 pe_chg=pe_change, pe_colht_cor=colht_cor)
786
787 pe_chg_k(k,4) = pe_chg_k(k,4) + pe_change
788 colht_cor_k(k,4) = colht_cor_k(k,4) + colht_cor
789
790
791 kd_so_far(k) = kd_so_far(k) + dkd
792
793 b1 = 1.0 / (hp_b(k) + kd_so_far(k))
794 c1_b(k) = kd_so_far(k) * b1
795
796 te(k) = b1 * (h_tr(k) * t0(k) + kd_so_far(k+1) * te(k+1))
797 se(k) = b1 * (h_tr(k) * s0(k) + kd_so_far(k+1) * se(k+1))
798
799 hp_b(k-1) = h_tr(k-1) + (hp_b(k) * b1) * kd_so_far(k)
800 dt_to_dpe_b(k-1) = dt_to_dpe(k-1) + c1_b(k)*dt_to_dpe_b(k)
801 ds_to_dpe_b(k-1) = ds_to_dpe(k-1) + c1_b(k)*ds_to_dpe_b(k)
802 dt_to_dcolht_b(k-1) = dt_to_dcolht(k-1) + c1_b(k)*dt_to_dcolht_b(k)
803 ds_to_dcolht_b(k-1) = ds_to_dcolht(k-1) + c1_b(k)*ds_to_dcolht_b(k)
804
805 enddo
806
807 ! Now update the other layer working down from the top to determine the
808 ! final temperatures and salinities.
809 b1 = 1.0 / (hp_b(1))
810 tf(1) = b1 * (h_tr(1) * t0(1) + kddt_h(2) * te(2))
811 sf(1) = b1 * (h_tr(1) * s0(1) + kddt_h(2) * se(2))
812 do k=2,nz
813 tf(k) = te(k) + c1_b(k)*tf(k-1)
814 sf(k) = se(k) + c1_b(k)*sf(k-1)
815 enddo
816
817 if (debug) then
818 pe_chg_tot1d = 0.0 ; pe_chg_tot2d = 0.0 ; t_chg_totd = 0.0
819 do k=1,nz
820 pe_chg_tot1d = pe_chg_tot1d + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
821 ds_to_dpe(k) * (sf(k) - s0(k)))
822 t_chg_totd = t_chg_totd + h_tr(k) * (tf(k) - t0(k))
823 enddo
824 do k=2,nz
825 pe_chg_tot2d = pe_chg_tot2d + (pe_chg_k(k,4) - colht_cor_k(k,4))
826 enddo
827 endif
828
829 endif
830
831 energy_kd = 0.0 ; do k=2,nz ; energy_kd = energy_kd + pe_chg_k(k,1) ; enddo
832 energy_kd = energy_kd / dt
833
834 if (do_print) then
835 if (cs%id_ERt>0) call post_data(cs%id_ERt, pe_chg_k(:,1), cs%diag)
836 if (cs%id_ERb>0) call post_data(cs%id_ERb, pe_chg_k(:,2), cs%diag)
837 if (cs%id_ERc>0) call post_data(cs%id_ERc, pe_chg_k(:,3), cs%diag)
838 if (cs%id_ERh>0) call post_data(cs%id_ERh, pe_chg_k(:,4), cs%diag)
839 if (cs%id_Kddt>0) call post_data(cs%id_Kddt, kddt_h, cs%diag)
840 if (cs%id_Kd>0) call post_data(cs%id_Kd, kd, cs%diag)
841 if (cs%id_h>0) call post_data(cs%id_h, h_tr, cs%diag)
842 if (cs%id_zInt>0) call post_data(cs%id_zInt, z_int, cs%diag)
843 if (cs%id_CHCt>0) call post_data(cs%id_CHCt, colht_cor_k(:,1), cs%diag)
844 if (cs%id_CHCb>0) call post_data(cs%id_CHCb, colht_cor_k(:,2), cs%diag)
845 if (cs%id_CHCc>0) call post_data(cs%id_CHCc, colht_cor_k(:,3), cs%diag)
846 if (cs%id_CHCh>0) call post_data(cs%id_CHCh, colht_cor_k(:,4), cs%diag)
847 if (cs%id_T0>0) call post_data(cs%id_T0, t0, cs%diag)
848 if (cs%id_Tf>0) call post_data(cs%id_Tf, tf, cs%diag)
849 if (cs%id_S0>0) call post_data(cs%id_S0, s0, cs%diag)
850 if (cs%id_Sf>0) call post_data(cs%id_Sf, sf, cs%diag)
851 if (cs%id_N2_0>0) then
852 n2(1) = 0.0 ; n2(nz+1) = 0.0
853 do k=2,nz
854 call calculate_density(0.5*(t0(k-1) + t0(k)), 0.5*(s0(k-1) + s0(k)), &
855 pres(k), rho_here, tv%eqn_of_state)
856 n2(k) = (gv%g_Earth_Z_T2 * rho_here / (0.5*(dz_tr(k-1) + dz_tr(k)))) * &
857 ( 0.5*(dsv_dt(k-1) + dsv_dt(k)) * (t0(k-1) - t0(k)) + &
858 0.5*(dsv_ds(k-1) + dsv_ds(k)) * (s0(k-1) - s0(k)) )
859 enddo
860 call post_data(cs%id_N2_0, n2, cs%diag)
861 endif
862 if (cs%id_N2_f>0) then
863 n2(1) = 0.0 ; n2(nz+1) = 0.0
864 do k=2,nz
865 call calculate_density(0.5*(tf(k-1) + tf(k)), 0.5*(sf(k-1) + sf(k)), &
866 pres(k), rho_here, tv%eqn_of_state)
867 n2(k) = (gv%g_Earth_Z_T2 * rho_here / (0.5*(dz_tr(k-1) + dz_tr(k)))) * &
868 ( 0.5*(dsv_dt(k-1) + dsv_dt(k)) * (tf(k-1) - tf(k)) + &
869 0.5*(dsv_ds(k-1) + dsv_ds(k)) * (sf(k-1) - sf(k)) )
870 enddo
871 call post_data(cs%id_N2_f, n2, cs%diag)
872 endif
873 endif
874
875end subroutine diapyc_energy_req_calc
876
877!> This subroutine calculates the change in potential energy and or derivatives
878!! for several changes in an interface's diapycnal diffusivity times a timestep.
879subroutine find_pe_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, &
880 dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, &
881 pres_Z, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, &
882 PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0, PE_ColHt_cor)
883 real, intent(in) :: Kddt_h0 !< The previously used diffusivity at an interface times
884 !! the time step and divided by the average of the
885 !! thicknesses around the interface [H ~> m or kg m-2].
886 real, intent(in) :: dKddt_h !< The trial change in the diffusivity at an interface times
887 !! the time step and divided by the average of the
888 !! thicknesses around the interface [H ~> m or kg m-2].
889 real, intent(in) :: hp_a !< The effective pivot thickness of the layer above the
890 !! interface, given by h_k plus a term that
891 !! is a fraction (determined from the tridiagonal solver) of
892 !! Kddt_h for the interface above [H ~> m or kg m-2].
893 real, intent(in) :: hp_b !< The effective pivot thickness of the layer below the
894 !! interface, given by h_k plus a term that
895 !! is a fraction (determined from the tridiagonal solver) of
896 !! Kddt_h for the interface above [H ~> m or kg m-2].
897 real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer
898 !! above, including implicit mixing effects with other
899 !! yet higher layers [C H ~> degC m or degC kg m-2].
900 real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer
901 !! above, including implicit mixing effects with other
902 !! yet higher layers [S H ~> ppt m or ppt kg m-2].
903 real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer
904 !! below, including implicit mixing effects with other
905 !! yet lower layers [C H ~> degC m or degC kg m-2].
906 real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer
907 !! below, including implicit mixing effects with other
908 !! yet lower layers [S H ~> ppt m or ppt kg m-2].
909 real, intent(in) :: dT_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
910 !! a layer's temperature change to the change in column potential
911 !! energy, including all implicit diffusive changes in the
912 !! temperatures of all the layers above [R Z L2 T-2 C-1 ~> J m-2 degC-1].
913 real, intent(in) :: dS_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
914 !! a layer's salinity change to the change in column potential
915 !! energy, including all implicit diffusive changes in the
916 !! salinities of all the layers above [R Z L2 T-2 S-1 ~> J m-2 ppt-1].
917 real, intent(in) :: dT_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
918 !! a layer's temperature change to the change in column potential
919 !! energy, including all implicit diffusive changes in the
920 !! temperatures of all the layers below [R Z L2 T-2 C-1 ~> J m-2 degC-1].
921 real, intent(in) :: dS_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
922 !! a layer's salinity change to the change in column potential
923 !! energy, including all implicit diffusive changes in the
924 !! salinities of all the layers below [R Z L2 T-2 S-1 ~> J m-2 ppt-1].
925 real, intent(in) :: pres_Z !< The hydrostatic interface pressure, which relates
926 !! the changes in column thickness to the energy that is radiated
927 !! as gravity waves and unavailable to drive mixing [R L2 T-2 ~> J m-3].
928 real, intent(in) :: dT_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dT) relating
929 !! a layer's temperature change to the change in column
930 !! height, including all implicit diffusive changes
931 !! in the temperatures of all the layers above [Z C-1 ~> m degC-1].
932 real, intent(in) :: dS_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dS) relating
933 !! a layer's salinity change to the change in column
934 !! height, including all implicit diffusive changes
935 !! in the salinities of all the layers above [Z S-1 ~> m ppt-1].
936 real, intent(in) :: dT_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dT) relating
937 !! a layer's temperature change to the change in column
938 !! height, including all implicit diffusive changes
939 !! in the temperatures of all the layers below [Z C-1 ~> m degC-1].
940 real, intent(in) :: dS_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dS) relating
941 !! a layer's salinity change to the change in column
942 !! height, including all implicit diffusive changes
943 !! in the salinities of all the layers below [Z S-1 ~> m ppt-1].
944
945 real, intent(out) :: PE_chg !< The change in column potential energy from applying
946 !! Kddt_h at the present interface [R Z L2 T-2 ~> J m-2].
947 real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h,
948 !! [R Z L2 T-2 H-1 ~> J m-3 or J kg-1].
949 real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could
950 !! be realized by applying a huge value of Kddt_h at the
951 !! present interface [R Z L2 T-2 ~> J m-2].
952 real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the
953 !! limit where Kddt_h = 0 [R Z L2 T-2 H-1 ~> J m-3 or J kg-1].
954 real, optional, intent(out) :: PE_ColHt_cor !< The correction to PE_chg that is made due to a net
955 !! change in the column height [R Z L2 T-2 ~> J m-2].
956
957 ! Local variables
958 real :: hps ! The sum of the two effective pivot thicknesses [H ~> m or kg m-2].
959 real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term [H2 ~> m2 or kg2 m-4].
960 real :: dT_c ! The core term in the expressions for the temperature changes [C H2 ~> degC m2 or degC kg2 m-4].
961 real :: dS_c ! The core term in the expressions for the salinity changes [S H2 ~> ppt m2 or ppt kg2 m-4].
962 real :: PEc_core ! The diffusivity-independent core term in the expressions
963 ! for the potential energy changes [H3 R Z L2 T-2 ~> J m or J kg3 m-8].
964 real :: ColHt_core ! The diffusivity-independent core term in the expressions
965 ! for the column height changes [H3 Z ~> m4 or kg3 m-5].
966 real :: ColHt_chg ! The change in the column height [Z ~> m].
967 real :: y1_3 ! A local temporary term in [H-3 ~> m-3 or m6 kg-3].
968 real :: y1_4 ! A local temporary term in [H-4 ~> m-4 or m8 kg-4].
969
970 ! The expression for the change in potential energy used here is derived
971 ! from the expression for the final estimates of the changes in temperature
972 ! and salinities, and then extensively manipulated to get it into its most
973 ! succinct form. The derivation is not necessarily obvious, but it demonstrably
974 ! works by comparison with separate calculations of the energy changes after
975 ! the tridiagonal solver for the final changes in temperature and salinity are
976 ! applied.
977
978 hps = hp_a + hp_b
979 bdt1 = hp_a * hp_b + kddt_h0 * hps
980 dt_c = hp_a * th_b - hp_b * th_a
981 ds_c = hp_a * sh_b - hp_b * sh_a
982 pec_core = hp_b * (dt_to_dpe_a * dt_c + ds_to_dpe_a * ds_c) - &
983 hp_a * (dt_to_dpe_b * dt_c + ds_to_dpe_b * ds_c)
984 colht_core = hp_b * (dt_to_dcolht_a * dt_c + ds_to_dcolht_a * ds_c) - &
985 hp_a * (dt_to_dcolht_b * dt_c + ds_to_dcolht_b * ds_c)
986
987 ! Find the change in column potential energy due to the change in the
988 ! diffusivity at this interface by dKddt_h.
989 y1_3 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
990 pe_chg = pec_core * y1_3
991 colht_chg = colht_core * y1_3
992 if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
993
994 if (present(pe_colht_cor)) pe_colht_cor = -pres_z * min(colht_chg, 0.0)
995
996 if (present(dpec_dkd)) then
997 ! Find the derivative of the potential energy change with dKddt_h.
998 y1_4 = 1.0 / (bdt1 + dkddt_h * hps)**2
999 dpec_dkd = pec_core * y1_4
1000 colht_chg = colht_core * y1_4
1001 if (colht_chg < 0.0) dpec_dkd = dpec_dkd - pres_z * colht_chg
1002 endif
1003
1004 if (present(dpe_max)) then
1005 ! This expression is the limit of PE_chg for infinite dKddt_h.
1006 y1_3 = 1.0 / (bdt1 * hps)
1007 dpe_max = pec_core * y1_3
1008 colht_chg = colht_core * y1_3
1009 if (colht_chg < 0.0) dpe_max = dpe_max - pres_z * colht_chg
1010 endif
1011
1012 if (present(dpec_dkd_0)) then
1013 ! This expression is the limit of dPEc_dKd for dKddt_h = 0.
1014 y1_4 = 1.0 / bdt1**2
1015 dpec_dkd_0 = pec_core * y1_4
1016 colht_chg = colht_core * y1_4
1017 if (colht_chg < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z * colht_chg
1018 endif
1019
1020end subroutine find_pe_chg
1021
1022
1023!> Initialize parameters and allocate memory associated with the diapycnal energy requirement module.
1024subroutine diapyc_energy_req_init(Time, G, GV, US, param_file, diag, CS)
1025 type(time_type), intent(in) :: time !< model time
1026 type(ocean_grid_type), intent(in) :: g !< model grid structure
1027 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1028 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1029 type(param_file_type), intent(in) :: param_file !< file to parse for parameter values
1030 type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output
1031 type(diapyc_energy_req_cs), pointer :: cs !< module control structure
1032
1033! This include declares and sets the variable "version".
1034#include "version_variable.h"
1035 character(len=40) :: mdl = "MOM_diapyc_energy_req" ! This module's name.
1036
1037 if (.not.associated(cs)) then ; allocate(cs)
1038 else ; return ; endif
1039
1040 cs%initialized = .true.
1041 cs%diag => diag
1042
1043 ! Read all relevant parameters and write them to the model log.
1044 call log_version(param_file, mdl, version, "")
1045 call get_param(param_file, mdl, "ENERGY_REQ_KH_SCALING", cs%test_Kh_scaling, &
1046 "A scaling factor for the diapycnal diffusivity used in "//&
1047 "testing the energy requirements.", default=1.0, units="nondim")
1048 call get_param(param_file, mdl, "ENERGY_REQ_COL_HT_SCALING", cs%ColHt_scaling, &
1049 "A scaling factor for the column height change correction "//&
1050 "used in testing the energy requirements.", default=1.0, units="nondim")
1051 call get_param(param_file, mdl, "ENERGY_REQ_USE_TEST_PROFILE", cs%use_test_Kh_profile, &
1052 "If true, use the internal test diffusivity profile in "//&
1053 "place of any that might be passed in as an argument.", default=.false.)
1054 call get_param(param_file, mdl, 'VON_KARMAN_CONST', cs%vonKar, &
1055 'The value the von Karman constant as used for mixed layer viscosity.', &
1056 units='nondim', default=0.41)
1057
1058 cs%id_ERt = register_diag_field('ocean_model', 'EnReqTest_ERt', diag%axesZi, time, &
1059 "Diffusivity Energy Requirements, top-down", &
1060 "J m-2", conversion=us%RZ_to_kg_m2*us%L_T_to_m_s**2)
1061 cs%id_ERb = register_diag_field('ocean_model', 'EnReqTest_ERb', diag%axesZi, time, &
1062 "Diffusivity Energy Requirements, bottom-up", &
1063 "J m-2", conversion=us%RZ_to_kg_m2*us%L_T_to_m_s**2)
1064 cs%id_ERc = register_diag_field('ocean_model', 'EnReqTest_ERc', diag%axesZi, time, &
1065 "Diffusivity Energy Requirements, center-last", &
1066 "J m-2", conversion=us%RZ_to_kg_m2*us%L_T_to_m_s**2)
1067 cs%id_ERh = register_diag_field('ocean_model', 'EnReqTest_ERh', diag%axesZi, time, &
1068 "Diffusivity Energy Requirements, halves", &
1069 "J m-2", conversion=us%RZ_to_kg_m2*us%L_T_to_m_s**2)
1070 cs%id_Kddt = register_diag_field('ocean_model', 'EnReqTest_Kddt', diag%axesZi, time, &
1071 "Implicit diffusive coupling coefficient", "m", conversion=gv%H_to_m)
1072 cs%id_Kd = register_diag_field('ocean_model', 'EnReqTest_Kd', diag%axesZi, time, &
1073 "Diffusivity in test", "m2 s-1", conversion=us%Z2_T_to_m2_s)
1074 cs%id_h = register_diag_field('ocean_model', 'EnReqTest_h_lay', diag%axesZL, time, &
1075 "Test column layer thicknesses", "m", conversion=gv%H_to_m)
1076 cs%id_zInt = register_diag_field('ocean_model', 'EnReqTest_z_int', diag%axesZi, time, &
1077 "Test column layer interface heights", "m", conversion=gv%H_to_m)
1078 cs%id_CHCt = register_diag_field('ocean_model', 'EnReqTest_CHCt', diag%axesZi, time, &
1079 "Column Height Correction to Energy Requirements, top-down", &
1080 "J m-2", conversion=us%RZ_to_kg_m2*us%L_T_to_m_s**2)
1081 cs%id_CHCb = register_diag_field('ocean_model', 'EnReqTest_CHCb', diag%axesZi, time, &
1082 "Column Height Correction to Energy Requirements, bottom-up", &
1083 "J m-2", conversion=us%RZ_to_kg_m2*us%L_T_to_m_s**2)
1084 cs%id_CHCc = register_diag_field('ocean_model', 'EnReqTest_CHCc', diag%axesZi, time, &
1085 "Column Height Correction to Energy Requirements, center-last", &
1086 "J m-2", conversion=us%RZ_to_kg_m2*us%L_T_to_m_s**2)
1087 cs%id_CHCh = register_diag_field('ocean_model', 'EnReqTest_CHCh', diag%axesZi, time, &
1088 "Column Height Correction to Energy Requirements, halves", &
1089 "J m-2", conversion=us%RZ_to_kg_m2*us%L_T_to_m_s**2)
1090 cs%id_T0 = register_diag_field('ocean_model', 'EnReqTest_T0', diag%axesZL, time, &
1091 "Temperature before mixing", "deg C", conversion=us%C_to_degC)
1092 cs%id_Tf = register_diag_field('ocean_model', 'EnReqTest_Tf', diag%axesZL, time, &
1093 "Temperature after mixing", "deg C", conversion=us%C_to_degC)
1094 cs%id_S0 = register_diag_field('ocean_model', 'EnReqTest_S0', diag%axesZL, time, &
1095 "Salinity before mixing", "g kg-1", conversion=us%S_to_ppt)
1096 cs%id_Sf = register_diag_field('ocean_model', 'EnReqTest_Sf', diag%axesZL, time, &
1097 "Salinity after mixing", "g kg-1", conversion=us%S_to_ppt)
1098 cs%id_N2_0 = register_diag_field('ocean_model', 'EnReqTest_N2_0', diag%axesZi, time, &
1099 "Squared buoyancy frequency before mixing", "second-2", conversion=us%s_to_T**2)
1100 cs%id_N2_f = register_diag_field('ocean_model', 'EnReqTest_N2_f', diag%axesZi, time, &
1101 "Squared buoyancy frequency after mixing", "second-2", conversion=us%s_to_T**2)
1102
1103end subroutine diapyc_energy_req_init
1104
1105!> Clean up and deallocate memory associated with the diapycnal energy requirement module.
1106subroutine diapyc_energy_req_end(CS)
1107 type(diapyc_energy_req_cs), pointer :: cs !< Diapycnal energy requirement control structure that
1108 !! will be deallocated in this subroutine.
1109 if (associated(cs)) deallocate(cs)
1110end subroutine diapyc_energy_req_end
1111
1112end module mom_diapyc_energy_req