MOM6
MOM_diapyc_energy_req.F90
Go to the documentation of this file.
1 !> Calculates the energy requirements of mixing.
2 module mom_diapyc_energy_req
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 !! \author By Robert Hallberg, May 2015
7 
8 use mom_diag_mediator, only : diag_ctrl, time_type, post_data, register_diag_field
9 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
10 use mom_file_parser, only : get_param, log_version, param_file_type
11 use mom_grid, only : ocean_grid_type
12 use mom_unit_scaling, only : unit_scale_type
13 use mom_variables, only : thermo_var_ptrs
14 use mom_verticalgrid, only : verticalgrid_type
15 use mom_eos, only : calculate_specific_vol_derivs, calculate_density
16 
17 implicit none ; private
18 
19 public diapyc_energy_req_init, diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_end
20 
21 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
22 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
23 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
24 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
25 
26 !> This control structure holds parameters for the MOM_diapyc_energy_req module
27 type, public :: diapyc_energy_req_cs ; private
28  logical :: initialized = .false. !< A variable that is here because empty
29  !! structures are not permitted by some compilers.
30  real :: test_kh_scaling !< A scaling factor for the diapycnal diffusivity.
31  real :: colht_scaling !< A scaling factor for the column height change correction term.
32  logical :: use_test_kh_profile !< If true, use the internal test diffusivity profile in place of
33  !! any that might be passed in as an argument.
34  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
35  !! regulate the timing of diagnostic output.
36 
37  !>@{ Diagnostic IDs
38  integer :: id_ert=-1, id_erb=-1, id_erc=-1, id_erh=-1, id_kddt=-1, id_kd=-1
39  integer :: id_chct=-1, id_chcb=-1, id_chcc=-1, id_chch=-1
40  integer :: id_t0=-1, id_tf=-1, id_s0=-1, id_sf=-1, id_n2_0=-1, id_n2_f=-1
41  integer :: id_h=-1, id_zint=-1
42  !!@}
43 end type diapyc_energy_req_cs
44 
45 contains
46 
47 !> This subroutine helps test the accuracy of the diapycnal mixing energy requirement code
48 !! by writing diagnostics, possibly using an intensely mixing test profile of diffusivity
49 subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int)
50  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
51  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
52  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
53  real, dimension(G%isd:G%ied,G%jsd:G%jed,GV%ke), &
54  intent(in) :: h_3d !< Layer thickness before entrainment [H ~> m or kg m-2].
55  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
56  !! available thermodynamic fields.
57  !! Absent fields have NULL ptrs.
58  real, intent(in) :: dt !< The amount of time covered by this call [T ~> s].
59  type(diapyc_energy_req_cs), pointer :: cs !< This module's control structure.
60  real, dimension(G%isd:G%ied,G%jsd:G%jed,GV%ke+1), &
61  optional, intent(in) :: kd_int !< Interface diffusivities [Z2 T-1 ~> m2 s-1].
62 
63  ! Local variables
64  real, dimension(GV%ke) :: &
65  t0, s0, & ! T0 & S0 are columns of initial temperatures and salinities [degC] and g/kg.
66  h_col ! h_col is a column of thicknesses h at tracer points [H ~> m or kg m-2].
67  real, dimension(GV%ke+1) :: &
68  kd, & ! A column of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1].
69  h_top, h_bot ! Distances from the top or bottom [H ~> m or kg m-2].
70  real :: ustar, absf, htot
71  real :: energy_kd ! The energy used by diapycnal mixing [W m-2].
72  real :: tmp1 ! A temporary array.
73  integer :: i, j, k, is, ie, js, je, nz, itt
74  logical :: may_print
75  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
76 
77  if (.not. associated(cs)) call mom_error(fatal, "diapyc_energy_req_test: "// &
78  "Module must be initialized before it is used.")
79 
80 !$OMP do
81  do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
82  if (present(kd_int) .and. .not.cs%use_test_Kh_profile) then
83  do k=1,nz+1 ; kd(k) = cs%test_Kh_scaling*kd_int(i,j,k) ; enddo
84  else
85  htot = 0.0 ; h_top(1) = 0.0
86  do k=1,nz
87  t0(k) = tv%T(i,j,k) ; s0(k) = tv%S(i,j,k)
88  h_col(k) = h_3d(i,j,k)
89  h_top(k+1) = h_top(k) + h_col(k)
90  enddo
91  htot = h_top(nz+1)
92  h_bot(nz+1) = 0.0
93  do k=nz,1,-1
94  h_bot(k) = h_bot(k+1) + h_col(k)
95  enddo
96 
97  ustar = 0.01*us%m_to_Z*us%T_to_s ! Change this to being an input parameter?
98  absf = 0.25*((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
99  (abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))))
100  kd(1) = 0.0 ; kd(nz+1) = 0.0
101  do k=2,nz
102  tmp1 = h_top(k) * h_bot(k) * gv%H_to_Z
103  kd(k) = cs%test_Kh_scaling * &
104  ustar * 0.41 * (tmp1*ustar) / (absf*tmp1 + htot*ustar)
105  enddo
106  endif
107  may_print = is_root_pe() .and. (i==ie) .and. (j==je)
108  call diapyc_energy_req_calc(h_col, t0, s0, kd, energy_kd, dt, tv, g, gv, us, &
109  may_print=may_print, cs=cs)
110  endif ; enddo ; enddo
111 
112 end subroutine diapyc_energy_req_test
113 
114 !> This subroutine uses a substantially refactored tridiagonal equation for
115 !! diapycnal mixing of temperature and salinity to estimate the potential energy
116 !! change due to diapycnal mixing in a column of water. It does this estimate
117 !! 4 different ways, all of which should be equivalent, but reports only one.
118 !! The various estimates are taken because they will later be used as templates
119 !! for other bits of code
120 subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, &
121  G, GV, US, may_print, CS)
122  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
123  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
124  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
125  real, dimension(GV%ke), intent(in) :: h_in !< Layer thickness before entrainment,
126  !! [H ~> m or kg m-2].
127  real, dimension(GV%ke), intent(in) :: t_in !< The layer temperatures [degC].
128  real, dimension(GV%ke), intent(in) :: s_in !< The layer salinities [ppt].
129  real, dimension(GV%ke+1), intent(in) :: kd !< The interfaces diapycnal diffusivities
130  !! [Z2 T-1 ~> m2 s-1].
131  real, intent(in) :: dt !< The amount of time covered by this call [T ~> s].
132  real, intent(out) :: energy_kd !< The column-integrated rate of energy
133  !! consumption by diapycnal diffusion [W m-2].
134  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
135  !! available thermodynamic fields.
136  !! Absent fields have NULL ptrs.
137  logical, optional, intent(in) :: may_print !< If present and true, write out diagnostics
138  !! of energy use.
139  type(diapyc_energy_req_cs), &
140  optional, pointer :: cs !< This module's control structure.
141 
142 ! This subroutine uses a substantially refactored tridiagonal equation for
143 ! diapycnal mixing of temperature and salinity to estimate the potential energy
144 ! change due to diapycnal mixing in a column of water. It does this estimate
145 ! 4 different ways, all of which should be equivalent, but reports only one.
146 ! The various estimates are taken because they will later be used as templates
147 ! for other bits of code.
148 
149  real, dimension(GV%ke) :: &
150  p_lay, & ! Average pressure of a layer [Pa].
151  dsv_dt, & ! Partial derivative of specific volume with temperature [m3 kg-1 degC-1].
152  dsv_ds, & ! Partial derivative of specific volume with salinity [m3 kg-1 ppt-1].
153  t0, s0, & ! Initial temperatures and salinities [degC] and [ppt].
154  te, se, & ! Running incomplete estimates of the new temperatures and salinities [degC] and [ppt].
155  te_a, se_a, & ! Running incomplete estimates of the new temperatures and salinities [degC] and [ppt].
156  te_b, se_b, & ! Running incomplete estimates of the new temperatures and salinities [degC] and [ppt].
157  tf, sf, & ! New final values of the temperatures and salinities [degC] and [ppt].
158  dte, dse, & ! Running (1-way) estimates of temperature and salinity change [degC] and [ppt].
159  dte_a, dse_a, & ! Running (1-way) estimates of temperature and salinity change [degC] and [ppt].
160  dte_b, dse_b, & ! Running (1-way) estimates of temperature and salinity change [degC] and [ppt].
161  th_a, & ! An effective temperature times a thickness in the layer above, including implicit
162  ! mixing effects with other yet higher layers [degC H ~> degC m or degC kg m-2].
163  sh_a, & ! An effective salinity times a thickness in the layer above, including implicit
164  ! mixing effects with other yet higher layers [ppt H ~> ppt m or ppt kg m-2].
165  th_b, & ! An effective temperature times a thickness in the layer below, including implicit
166  ! mixing effects with other yet lower layers [degC H ~> degC m or degC kg m-2].
167  sh_b, & ! An effective salinity times a thickness in the layer below, including implicit
168  ! mixing effects with other yet lower layers [ppt H ~> ppt m or ppt kg m-2].
169  dt_to_dpe, & ! Partial derivative of column potential energy with the temperature
170  ds_to_dpe, & ! and salinity changes within a layer [J m-2 degC-1] and [J m-2 ppt-1].
171  dt_to_dcolht, & ! Partial derivatives of the total column height with the temperature
172  ds_to_dcolht, & ! and salinity changes within a layer [Z degC-1 ~> m degC-1] and [Z ppt-1 ~> m ppt-1].
173  dt_to_dcolht_a, & ! Partial derivatives of the total column height with the temperature
174  ds_to_dcolht_a, & ! and salinity changes within a layer, including the implicit effects
175  ! of mixing with layers higher in the water colun [Z degC-1 ~> m degC-1] and [Z ppt-1 ~> m ppt-1].
176  dt_to_dcolht_b, & ! Partial derivatives of the total column height with the temperature
177  ds_to_dcolht_b, & ! and salinity changes within a layer, including the implicit effects
178  ! of mixing with layers lower in the water colun [Z degC-1 ~> m degC-1] and [Z ppt-1 ~> m ppt-1].
179  dt_to_dpe_a, & ! Partial derivatives of column potential energy with the temperature
180  ds_to_dpe_a, & ! and salinity changes within a layer, including the implicit effects
181  ! of mixing with layers higher in the water column, in
182  ! units of [J m-2 degC-1] and [J m-2 ppt-1].
183  dt_to_dpe_b, & ! Partial derivatives of column potential energy with the temperature
184  ds_to_dpe_b, & ! and salinity changes within a layer, including the implicit effects
185  ! of mixing with layers lower in the water column, in
186  ! units of [J m-2 degC-1] and [J m-2 ppt-1].
187  hp_a, & ! An effective pivot thickness of the layer including the effects
188  ! of coupling with layers above [H ~> m or kg m-2]. This is the first term
189  ! in the denominator of b1 in a downward-oriented tridiagonal solver.
190  hp_b, & ! An effective pivot thickness of the layer including the effects
191  ! of coupling with layers below [H ~> m or kg m-2]. This is the first term
192  ! in the denominator of b1 in an upward-oriented tridiagonal solver.
193  c1_a, & ! c1_a is used by a downward-oriented tridiagonal solver [nondim].
194  c1_b, & ! c1_b is used by an upward-oriented tridiagonal solver [nondim].
195  h_tr ! h_tr is h at tracer points with a h_neglect added to
196  ! ensure positive definiteness [H ~> m or kg m-2].
197  real, dimension(GV%ke+1) :: &
198  pres, & ! Interface pressures [Pa].
199  pres_z, & ! Interface pressures with a rescaling factor to convert interface height
200  ! movements into changes in column potential energy [J m-2 Z-1 ~> J m-3].
201  z_int, & ! Interface heights relative to the surface [H ~> m or kg m-2].
202  n2, & ! An estimate of the buoyancy frequency [T-2 ~> s-2].
203  kddt_h, & ! The diapycnal diffusivity times a timestep divided by the
204  ! average thicknesses around a layer [H ~> m or kg m-2].
205  kddt_h_a, & ! The value of Kddt_h for layers above the central point in the
206  ! tridiagonal solver [H ~> m or kg m-2].
207  kddt_h_b, & ! The value of Kddt_h for layers below the central point in the
208  ! tridiagonal solver [H ~> m or kg m-2].
209  kd_so_far ! The value of Kddt_h that has been applied already in
210  ! calculating the energy changes [H ~> m or kg m-2].
211  real, dimension(GV%ke+1,4) :: &
212  pe_chg_k, & ! The integrated potential energy change within a timestep due
213  ! to the diffusivity at interface K for 4 different orders of
214  ! accumulating the diffusivities [J m-2].
215  colht_cor_k ! The correction to the potential energy change due to
216  ! changes in the net column height [J m-2].
217  real :: &
218  b1 ! b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
219  real :: &
220  i_b1 ! The inverse of b1 [H ~> m or kg m-2].
221  real :: kd0 ! The value of Kddt_h that has already been applied [H ~> m or kg m-2].
222  real :: dkd ! The change in the value of Kddt_h [H ~> m or kg m-2].
223  real :: h_neglect ! A thickness that is so small it is usually lost
224  ! in roundoff and can be neglected [H ~> m or kg m-2].
225  real :: dte_term ! A diffusivity-independent term related to the temperature
226  ! change in the layer below the interface [degC H ~> degC m or degC kg m-2].
227  real :: dse_term ! A diffusivity-independent term related to the salinity
228  ! change in the layer below the interface [ppt H ~> ppt m or ppt kg m-2].
229  real :: kddt_h_guess ! A guess of the final value of Kddt_h [H ~> m or kg m-2].
230  real :: dmass ! The mass per unit area within a layer [kg m-2].
231  real :: dpres ! The hydrostatic pressure change across a layer [Pa].
232  real :: dmke_max ! The maximum amount of mean kinetic energy that could be
233  ! converted to turbulent kinetic energy if the velocity in
234  ! the layer below an interface were homogenized with all of
235  ! the water above the interface [J m-2 = kg s-2].
236  real :: rho_here ! The in-situ density [kg m-3].
237  real :: pe_change ! The change in column potential energy from applying Kddt_h at the
238  ! present interface [J m-2].
239  real :: colht_cor ! The correction to PE_chg that is made due to a net
240  ! change in the column height [J m-2].
241  real :: htot ! A running sum of thicknesses [H ~> m or kg m-2].
242  real :: dte_t2, dt_km1_t2, dt_k_t2 ! Temporary arrays describing temperature changes [degC].
243  real :: dse_t2, ds_km1_t2, ds_k_t2 ! Temporary arrays describing salinity changes [ppt].
244  logical :: do_print
245 
246  ! The following are a bunch of diagnostic arrays for debugging purposes.
247  real, dimension(GV%ke) :: &
248  ta, sa, tb, sb
249  real, dimension(GV%ke+1) :: &
250  dpea_dkd, dpea_dkd_est, dpea_dkd_err, dpea_dkd_trunc, dpea_dkd_err_norm, &
251  dpeb_dkd, dpeb_dkd_est, dpeb_dkd_err, dpeb_dkd_trunc, dpeb_dkd_err_norm
252  real :: pe_chg_tot1a, pe_chg_tot2a, t_chg_tota
253  real :: pe_chg_tot1b, pe_chg_tot2b, t_chg_totb
254  real :: pe_chg_tot1c, pe_chg_tot2c, t_chg_totc
255  real :: pe_chg_tot1d, pe_chg_tot2d, t_chg_totd
256  real, dimension(GV%ke+1) :: dpechg_dkd
257  real :: pe_chg(6)
258  real, dimension(6) :: dt_k_itt, ds_k_itt, dt_km1_itt, ds_km1_itt
259 
260  integer :: k, nz, itt, max_itt, k_cent
261  logical :: surface_bl, bottom_bl, central, halves, debug
262  logical :: old_pe_calc
263  nz = g%ke
264  h_neglect = gv%H_subroundoff
265 
266  debug = .true.
267 
268  surface_bl = .true. ; bottom_bl = .true. ; halves = .true.
269  central = .true. ; k_cent = nz/2
270 
271  do_print = .false. ; if (present(may_print) .and. present(cs)) do_print = may_print
272 
273  dpea_dkd(:) = 0.0 ; dpea_dkd_est(:) = 0.0 ; dpea_dkd_err(:) = 0.0
274  dpea_dkd_err_norm(:) = 0.0 ; dpea_dkd_trunc(:) = 0.0
275  dpeb_dkd(:) = 0.0 ; dpeb_dkd_est(:) = 0.0 ; dpeb_dkd_err(:) = 0.0
276  dpeb_dkd_err_norm(:) = 0.0 ; dpeb_dkd_trunc(:) = 0.0
277 
278  htot = 0.0 ; pres(1) = 0.0 ; pres_z(1) = 0.0 ; z_int(1) = 0.0
279  do k=1,nz
280  t0(k) = t_in(k) ; s0(k) = s_in(k)
281  h_tr(k) = h_in(k)
282  htot = htot + h_tr(k)
283  pres(k+1) = pres(k) + gv%H_to_Pa * h_tr(k)
284  pres_z(k+1) = us%Z_to_m * pres(k+1)
285  p_lay(k) = 0.5*(pres(k) + pres(k+1))
286  z_int(k+1) = z_int(k) - h_tr(k)
287  enddo
288  do k=1,nz
289  h_tr(k) = max(h_tr(k), 1e-15*htot)
290  enddo
291 
292  ! Introduce a diffusive flux variable, Kddt_h(K) = ea(k) = eb(k-1)
293 
294  kddt_h(1) = 0.0 ; kddt_h(nz+1) = 0.0
295  do k=2,nz
296  kddt_h(k) = min((gv%Z_to_H**2*dt)*kd(k) / (0.5*(h_tr(k-1) + h_tr(k))), 1e3*htot)
297  enddo
298 
299  ! Solve the tridiagonal equations for new temperatures.
300 
301  call calculate_specific_vol_derivs(t0, s0, p_lay, dsv_dt, dsv_ds, 1, nz, tv%eqn_of_state)
302 
303  do k=1,nz
304  dmass = gv%H_to_kg_m2 * h_tr(k)
305  dpres = gv%H_to_Pa * h_tr(k)
306  dt_to_dpe(k) = (dmass * (pres(k) + 0.5*dpres)) * dsv_dt(k)
307  ds_to_dpe(k) = (dmass * (pres(k) + 0.5*dpres)) * dsv_ds(k)
308  dt_to_dcolht(k) = dmass * us%m_to_Z * dsv_dt(k) * cs%ColHt_scaling
309  ds_to_dcolht(k) = dmass * us%m_to_Z * dsv_ds(k) * cs%ColHt_scaling
310  enddo
311 
312 ! PE_chg_k(1) = 0.0 ; PE_chg_k(nz+1) = 0.0
313  ! PEchg(:) = 0.0
314  pe_chg_k(:,:) = 0.0 ; colht_cor_k(:,:) = 0.0
315  dpechg_dkd(:) = 0.0
316 
317  if (surface_bl) then ! This version is appropriate for a surface boundary layer.
318  old_pe_calc = .false.
319 
320  ! Set up values appropriate for no diffusivity.
321  do k=1,nz
322  hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
323  dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
324  dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
325  dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
326  dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
327  enddo
328 
329  do k=2,nz
330  ! At this point, Kddt_h(K) will be unknown because its value may depend
331  ! on how much energy is available.
332 
333  ! Precalculate some temporary expressions that are independent of Kddt_h_guess.
334  if (old_pe_calc) then
335  if (k==2) then
336  dt_km1_t2 = (t0(k)-t0(k-1))
337  ds_km1_t2 = (s0(k)-s0(k-1))
338  dte_t2 = 0.0 ; dse_t2 = 0.0
339  else
340  dte_t2 = kddt_h(k-1) * ((t0(k-2) - t0(k-1)) + dte(k-2))
341  dse_t2 = kddt_h(k-1) * ((s0(k-2) - s0(k-1)) + dse(k-2))
342  dt_km1_t2 = (t0(k)-t0(k-1)) - &
343  (kddt_h(k-1) / hp_a(k-1)) * ((t0(k-2) - t0(k-1)) + dte(k-2))
344  ds_km1_t2 = (s0(k)-s0(k-1)) - &
345  (kddt_h(k-1) / hp_a(k-1)) * ((s0(k-2) - s0(k-1)) + dse(k-2))
346  endif
347  dte_term = dte_t2 + hp_a(k-1) * (t0(k-1)-t0(k))
348  dse_term = dse_t2 + hp_a(k-1) * (s0(k-1)-s0(k))
349  else
350  if (k<=2) then
351  th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
352  else
353  th_a(k-1) = h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2)
354  sh_a(k-1) = h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2)
355  endif
356  th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
357  endif
358 
359 
360  ! Find the energy change due to a guess at the strength of diffusion at interface K.
361 
362  kddt_h_guess = kddt_h(k)
363  if (old_pe_calc) then
364  call find_pe_chg_orig(kddt_h_guess, h_tr(k), hp_a(k-1), &
365  dte_term, dse_term, dt_km1_t2, ds_km1_t2, &
366  dt_to_dpe(k), ds_to_dpe(k), dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
367  pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
368  dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
369  pe_chg_k(k,1), dpea_dkd(k))
370  else
371  call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
372  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
373  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
374  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
375  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
376  pe_chg=pe_chg_k(k,1), dpec_dkd=dpea_dkd(k), &
377  colht_cor=colht_cor_k(k,1))
378  endif
379 
380  if (debug) then
381  do itt=1,5
382  kddt_h_guess = (1.0+0.01*(itt-3))*kddt_h(k)
383 
384  if (old_pe_calc) then
385  call find_pe_chg_orig(kddt_h_guess, h_tr(k), hp_a(k-1), &
386  dte_term, dse_term, dt_km1_t2, ds_km1_t2, &
387  dt_to_dpe(k), ds_to_dpe(k), dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
388  pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
389  dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
390  pe_chg=pe_chg(itt))
391  else
392  call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
393  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
394  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
395  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
396  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
397  pe_chg=pe_chg(itt))
398  endif
399  enddo
400  ! Compare with a 4th-order finite difference estimate.
401  dpea_dkd_est(k) = (4.0*(pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
402  (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))) / 3.0
403  dpea_dkd_trunc(k) = (pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
404  (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))
405  dpea_dkd_err(k) = (dpea_dkd_est(k) - dpea_dkd(k))
406  dpea_dkd_err_norm(k) = (dpea_dkd_est(k) - dpea_dkd(k)) / &
407  (abs(dpea_dkd_est(k)) + abs(dpea_dkd(k)) + 1e-100)
408  endif
409 
410  ! At this point, the final value of Kddt_h(K) is known, so the estimated
411  ! properties for layer k-1 can be calculated.
412 
413  b1 = 1.0 / (hp_a(k-1) + kddt_h(k))
414  c1_a(k) = kddt_h(k) * b1
415  if (k==2) then
416  te(1) = b1*(h_tr(1)*t0(1))
417  se(1) = b1*(h_tr(1)*s0(1))
418  else
419  te(k-1) = b1 * (h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2))
420  se(k-1) = b1 * (h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2))
421  endif
422  if (old_pe_calc) then
423  dte(k-1) = b1 * ( kddt_h(k)*(t0(k)-t0(k-1)) + dte_t2 )
424  dse(k-1) = b1 * ( kddt_h(k)*(s0(k)-s0(k-1)) + dse_t2 )
425  endif
426 
427  hp_a(k) = h_tr(k) + (hp_a(k-1) * b1) * kddt_h(k)
428  dt_to_dpe_a(k) = dt_to_dpe(k) + c1_a(k)*dt_to_dpe_a(k-1)
429  ds_to_dpe_a(k) = ds_to_dpe(k) + c1_a(k)*ds_to_dpe_a(k-1)
430  dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1_a(k)*dt_to_dcolht_a(k-1)
431  ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1_a(k)*ds_to_dcolht_a(k-1)
432 
433  enddo
434 
435  b1 = 1.0 / (hp_a(nz))
436  tf(nz) = b1 * (h_tr(nz) * t0(nz) + kddt_h(nz) * te(nz-1))
437  sf(nz) = b1 * (h_tr(nz) * s0(nz) + kddt_h(nz) * se(nz-1))
438  if (old_pe_calc) then
439  dte(nz) = b1 * kddt_h(nz) * ((t0(nz-1)-t0(nz)) + dte(nz-1))
440  dse(nz) = b1 * kddt_h(nz) * ((s0(nz-1)-s0(nz)) + dse(nz-1))
441  endif
442 
443  do k=nz-1,1,-1
444  tf(k) = te(k) + c1_a(k+1)*tf(k+1)
445  sf(k) = se(k) + c1_a(k+1)*sf(k+1)
446  enddo
447 
448  if (debug) then
449  do k=1,nz ; ta(k) = tf(k) ; sa(k) = sf(k) ; enddo
450  pe_chg_tot1a = 0.0 ; pe_chg_tot2a = 0.0 ; t_chg_tota = 0.0
451  do k=1,nz
452  pe_chg_tot1a = pe_chg_tot1a + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
453  ds_to_dpe(k) * (sf(k) - s0(k)))
454  t_chg_tota = t_chg_tota + h_tr(k) * (tf(k) - t0(k))
455  enddo
456  do k=2,nz
457  pe_chg_tot2a = pe_chg_tot2a + (pe_chg_k(k,1) - colht_cor_k(k,1))
458  enddo
459  endif
460 
461  endif
462 
463  if (bottom_bl) then ! This version is appropriate for a bottom boundary layer.
464  old_pe_calc = .false.
465 
466  ! Set up values appropriate for no diffusivity.
467  do k=1,nz
468  hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
469  dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
470  dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
471  dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
472  dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
473  enddo
474 
475  do k=nz,2,-1 ! Loop over interior interfaces.
476  ! At this point, Kddt_h(K) will be unknown because its value may depend
477  ! on how much energy is available.
478 
479  ! Precalculate some temporary expressions that are independent of Kddt_h_guess.
480  if (old_pe_calc) then
481  if (k==nz) then
482  dt_k_t2 = (t0(k-1)-t0(k))
483  ds_k_t2 = (s0(k-1)-s0(k))
484  dte_t2 = 0.0 ; dse_t2 = 0.0
485  else
486  dte_t2 = kddt_h(k+1) * ((t0(k+1) - t0(k)) + dte(k+1))
487  dse_t2 = kddt_h(k+1) * ((s0(k+1) - s0(k)) + dse(k+1))
488  dt_k_t2 = (t0(k-1)-t0(k)) - &
489  (kddt_h(k+1)/ hp_b(k)) * ((t0(k+1) - t0(k)) + dte(k+1))
490  ds_k_t2 = (s0(k-1)-s0(k)) - &
491  (kddt_h(k+1)/ hp_b(k)) * ((s0(k+1) - s0(k)) + dse(k+1))
492  endif
493  dte_term = dte_t2 + hp_b(k) * (t0(k)-t0(k-1))
494  dse_term = dse_t2 + hp_b(k) * (s0(k)-s0(k-1))
495  else
496  th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
497  if (k>=nz) then
498  th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
499  else
500  th_b(k) = h_tr(k) * t0(k) + kddt_h(k+1) * te(k+1)
501  sh_b(k) = h_tr(k) * s0(k) + kddt_h(k+1) * se(k+1)
502  endif
503  endif
504 
505  ! Find the energy change due to a guess at the strength of diffusion at interface K.
506  kddt_h_guess = kddt_h(k)
507 
508  if (old_pe_calc) then
509  call find_pe_chg_orig(kddt_h_guess, h_tr(k-1), hp_b(k), &
510  dte_term, dse_term, dt_k_t2, ds_k_t2, &
511  dt_to_dpe(k-1), ds_to_dpe(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
512  pres_z(k), dt_to_dcolht(k-1), ds_to_dcolht(k-1), &
513  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
514  pe_chg=pe_chg_k(k,2), dpec_dkd=dpeb_dkd(k))
515  else
516  call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
517  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
518  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
519  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
520  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
521  pe_chg=pe_chg_k(k,2), dpec_dkd=dpeb_dkd(k), &
522  colht_cor=colht_cor_k(k,2))
523  endif
524 
525  if (debug) then
526  ! Compare with a 4th-order finite difference estimate.
527  do itt=1,5
528  kddt_h_guess = (1.0+0.01*(itt-3))*kddt_h(k)
529 
530  if (old_pe_calc) then
531  call find_pe_chg_orig(kddt_h_guess, h_tr(k-1), hp_b(k), &
532  dte_term, dse_term, dt_k_t2, ds_k_t2, &
533  dt_to_dpe(k-1), ds_to_dpe(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
534  pres_z(k), dt_to_dcolht(k-1), ds_to_dcolht(k-1), &
535  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
536  pe_chg=pe_chg(itt))
537  else
538  call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
539  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
540  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
541  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
542  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
543  pe_chg=pe_chg(itt))
544  endif
545  enddo
546 
547  dpeb_dkd_est(k) = (4.0*(pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
548  (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))) / 3.0
549  dpeb_dkd_trunc(k) = (pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
550  (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))
551  dpeb_dkd_err(k) = (dpeb_dkd_est(k) - dpeb_dkd(k))
552  dpeb_dkd_err_norm(k) = (dpeb_dkd_est(k) - dpeb_dkd(k)) / &
553  (abs(dpeb_dkd_est(k)) + abs(dpeb_dkd(k)) + 1e-100)
554  endif
555 
556  ! At this point, the final value of Kddt_h(K) is known, so the estimated
557  ! properties for layer k can be calculated.
558 
559  b1 = 1.0 / (hp_b(k) + kddt_h(k))
560  c1_b(k) = kddt_h(k) * b1
561  if (k==nz) then
562  te(nz) = b1* (h_tr(nz)*t0(nz))
563  se(nz) = b1* (h_tr(nz)*s0(nz))
564  else
565  te(k) = b1 * (h_tr(k) * t0(k) + kddt_h(k+1) * te(k+1))
566  se(k) = b1 * (h_tr(k) * s0(k) + kddt_h(k+1) * se(k+1))
567  endif
568  if (old_pe_calc) then
569  dte(k) = b1 * ( kddt_h(k)*(t0(k-1)-t0(k)) + dte_t2 )
570  dse(k) = b1 * ( kddt_h(k)*(s0(k-1)-s0(k)) + dse_t2 )
571  endif
572 
573  hp_b(k-1) = h_tr(k-1) + (hp_b(k) * b1) * kddt_h(k)
574  dt_to_dpe_b(k-1) = dt_to_dpe(k-1) + c1_b(k)*dt_to_dpe_b(k)
575  ds_to_dpe_b(k-1) = ds_to_dpe(k-1) + c1_b(k)*ds_to_dpe_b(k)
576  dt_to_dcolht_b(k-1) = dt_to_dcolht(k-1) + c1_b(k)*dt_to_dcolht_b(k)
577  ds_to_dcolht_b(k-1) = ds_to_dcolht(k-1) + c1_b(k)*ds_to_dcolht_b(k)
578 
579  enddo
580 
581  b1 = 1.0 / (hp_b(1))
582  tf(1) = b1 * (h_tr(1) * t0(1) + kddt_h(2) * te(2))
583  sf(1) = b1 * (h_tr(1) * s0(1) + kddt_h(2) * se(2))
584  if (old_pe_calc) then
585  dte(1) = b1 * kddt_h(2) * ((t0(2)-t0(1)) + dte(2))
586  dse(1) = b1 * kddt_h(2) * ((s0(2)-s0(1)) + dse(2))
587  endif
588 
589  do k=2,nz
590  tf(k) = te(k) + c1_b(k)*tf(k-1)
591  sf(k) = se(k) + c1_b(k)*sf(k-1)
592  enddo
593 
594  if (debug) then
595  do k=1,nz ; tb(k) = tf(k) ; sb(k) = sf(k) ; enddo
596  pe_chg_tot1b = 0.0 ; pe_chg_tot2b = 0.0 ; t_chg_totb = 0.0
597  do k=1,nz
598  pe_chg_tot1b = pe_chg_tot1b + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
599  ds_to_dpe(k) * (sf(k) - s0(k)))
600  t_chg_totb = t_chg_totb + h_tr(k) * (tf(k) - t0(k))
601  enddo
602  do k=2,nz
603  pe_chg_tot2b = pe_chg_tot2b + (pe_chg_k(k,2) - colht_cor_k(k,2))
604  enddo
605  endif
606 
607  endif
608 
609  if (central) then
610 
611  ! Set up values appropriate for no diffusivity.
612  do k=1,nz
613  hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
614  dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
615  dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
616  dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
617  dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
618  enddo
619 
620  ! Calculate the dependencies on layers above.
621  kddt_h_a(1) = 0.0
622  do k=2,nz ! Loop over interior interfaces.
623  ! First calculate some terms that are independent of the change in Kddt_h(K).
624  kd0 = 0.0 ! This might need to be changed - it is the already applied value of Kddt_h(K).
625  if (k<=2) then
626  th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
627  else
628  th_a(k-1) = h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te_a(k-2)
629  sh_a(k-1) = h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se_a(k-2)
630  endif
631  th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
632 
633  kddt_h_a(k) = 0.0 ; if (k < k_cent) kddt_h_a(k) = kddt_h(k)
634  dkd = kddt_h_a(k)
635 
636  call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
637  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
638  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
639  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
640  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
641  pe_chg=pe_change, colht_cor=colht_cor)
642  pe_chg_k(k,3) = pe_change
643  colht_cor_k(k,3) = colht_cor
644 
645  b1 = 1.0 / (hp_a(k-1) + kddt_h_a(k))
646  c1_a(k) = kddt_h_a(k) * b1
647  if (k==2) then
648  te_a(1) = b1*(h_tr(1)*t0(1))
649  se_a(1) = b1*(h_tr(1)*s0(1))
650  else
651  te_a(k-1) = b1 * (h_tr(k-1) * t0(k-1) + kddt_h_a(k-1) * te_a(k-2))
652  se_a(k-1) = b1 * (h_tr(k-1) * s0(k-1) + kddt_h_a(k-1) * se_a(k-2))
653  endif
654 
655  hp_a(k) = h_tr(k) + (hp_a(k-1) * b1) * kddt_h_a(k)
656  dt_to_dpe_a(k) = dt_to_dpe(k) + c1_a(k)*dt_to_dpe_a(k-1)
657  ds_to_dpe_a(k) = ds_to_dpe(k) + c1_a(k)*ds_to_dpe_a(k-1)
658  dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1_a(k)*dt_to_dcolht_a(k-1)
659  ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1_a(k)*ds_to_dcolht_a(k-1)
660  enddo
661 
662  ! Calculate the dependencies on layers below.
663  kddt_h_b(nz+1) = 0.0
664  do k=nz,2,-1 ! Loop over interior interfaces.
665  ! First calculate some terms that are independent of the change in Kddt_h(K).
666  kd0 = 0.0 ! This might need to be changed - it is the already applied value of Kddt_h(K).
667 ! if (K<=2) then
668  th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
669 ! else
670 ! Th_a(k-1) = h_tr(k-1) * T0(k-1) + Kddt_h(K-1) * Te_a(k-2)
671 ! Sh_a(k-1) = h_tr(k-1) * S0(k-1) + Kddt_h(K-1) * Se_a(k-2)
672 ! endif
673  if (k>=nz) then
674  th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
675  else
676  th_b(k) = h_tr(k) * t0(k) + kddt_h(k+1) * te_b(k+1)
677  sh_b(k) = h_tr(k) * s0(k) + kddt_h(k+1) * se_b(k+1)
678  endif
679 
680  kddt_h_b(k) = 0.0 ; if (k > k_cent) kddt_h_b(k) = kddt_h(k)
681  dkd = kddt_h_b(k)
682 
683  call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
684  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
685  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
686  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
687  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
688  pe_chg=pe_change, colht_cor=colht_cor)
689  pe_chg_k(k,3) = pe_chg_k(k,3) + pe_change
690  colht_cor_k(k,3) = colht_cor_k(k,3) + colht_cor
691 
692  b1 = 1.0 / (hp_b(k) + kddt_h_b(k))
693  c1_b(k) = kddt_h_b(k) * b1
694  if (k==nz) then
695  te_b(k) = b1 * (h_tr(k)*t0(k))
696  se_b(k) = b1 * (h_tr(k)*s0(k))
697  else
698  te_b(k) = b1 * (h_tr(k) * t0(k) + kddt_h_b(k+1) * te_b(k+1))
699  se_b(k) = b1 * (h_tr(k) * s0(k) + kddt_h_b(k+1) * se_b(k+1))
700  endif
701 
702  hp_b(k-1) = h_tr(k-1) + (hp_b(k) * b1) * kddt_h_b(k)
703  dt_to_dpe_b(k-1) = dt_to_dpe(k-1) + c1_b(k)*dt_to_dpe_b(k)
704  ds_to_dpe_b(k-1) = ds_to_dpe(k-1) + c1_b(k)*ds_to_dpe_b(k)
705  dt_to_dcolht_b(k-1) = dt_to_dcolht(k-1) + c1_b(k)*dt_to_dcolht_b(k)
706  ds_to_dcolht_b(k-1) = ds_to_dcolht(k-1) + c1_b(k)*ds_to_dcolht_b(k)
707 
708  enddo
709 
710  ! Calculate the final solution for the layers surrounding interface K_cent
711  k = k_cent
712 
713  ! First calculate some terms that are independent of the change in Kddt_h(K).
714  kd0 = 0.0 ! This might need to be changed - it is the already applied value of Kddt_h(K).
715  if (k<=2) then
716  th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
717  else
718  th_a(k-1) = h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te_a(k-2)
719  sh_a(k-1) = h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se_a(k-2)
720  endif
721  if (k>=nz) then
722  th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
723  else
724  th_b(k) = h_tr(k) * t0(k) + kddt_h(k+1) * te_b(k+1)
725  sh_b(k) = h_tr(k) * s0(k) + kddt_h(k+1) * se_b(k+1)
726  endif
727 
728  dkd = kddt_h(k)
729 
730  call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
731  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
732  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
733  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
734  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
735  pe_chg=pe_change, colht_cor=colht_cor)
736  pe_chg_k(k,3) = pe_chg_k(k,3) + pe_change
737  colht_cor_k(k,3) = colht_cor_k(k,3) + colht_cor
738 
739 
740  ! To derive the following, first to a partial update for the estimated
741  ! temperatures and salinities in the layers around this interface:
742  ! b1_a = 1.0 / (hp_a(k-1) + Kddt_h(K))
743  ! b1_b = 1.0 / (hp_b(k) + Kddt_h(K))
744  ! Te_up(k) = Th_b(k) * b1_b ; Se_up(k) = Sh_b(k) * b1_b
745  ! Te_up(k-1) = Th_a(k-1) * b1_a ; Se_up(k-1) = Sh_a(k-1) * b1_a
746  ! Find the final values of T & S for both layers around K_cent, using that
747  ! c1_a(K) = Kddt_h(K) * b1_a ; c1_b(K) = Kddt_h(K) * b1_b
748  ! Tf(K_cent-1) = Te_up(K_cent-1) + c1_a(K_cent)*Tf(K_cent)
749  ! Tf(K_cent) = Te_up(K_cent) + c1_b(K_cent)*Tf(K_cent-1)
750  ! but further exploiting the expressions for c1_a and c1_b to avoid
751  ! subtraction in the denominator, and use only a single division.
752  b1 = 1.0 / (hp_a(k-1)*hp_b(k) + kddt_h(k)*(hp_a(k-1) + hp_b(k)))
753  tf(k-1) = ((hp_b(k) + kddt_h(k)) * th_a(k-1) + kddt_h(k) * th_b(k) ) * b1
754  sf(k-1) = ((hp_b(k) + kddt_h(k)) * sh_a(k-1) + kddt_h(k) * sh_b(k) ) * b1
755  tf(k) = (kddt_h(k) * th_a(k-1) + (hp_a(k-1) + kddt_h(k)) * th_b(k) ) * b1
756  sf(k) = (kddt_h(k) * sh_a(k-1) + (hp_a(k-1) + kddt_h(k)) * sh_b(k) ) * b1
757 
758  c1_a(k) = kddt_h(k) / (hp_a(k-1) + kddt_h(k))
759  c1_b(k) = kddt_h(k) / (hp_b(k) + kddt_h(k))
760 
761  ! Now update the other layer working outward from k_cent to determine the final
762  ! temperatures and salinities.
763  do k=k_cent-2,1,-1
764  tf(k) = te_a(k) + c1_a(k+1)*tf(k+1)
765  sf(k) = se_a(k) + c1_a(k+1)*sf(k+1)
766  enddo
767  do k=k_cent+1,nz
768  tf(k) = te_b(k) + c1_b(k)*tf(k-1)
769  sf(k) = se_b(k) + c1_b(k)*sf(k-1)
770  enddo
771 
772  if (debug) then
773  pe_chg_tot1c = 0.0 ; pe_chg_tot2c = 0.0 ; t_chg_totc = 0.0
774  do k=1,nz
775  pe_chg_tot1c = pe_chg_tot1c + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
776  ds_to_dpe(k) * (sf(k) - s0(k)))
777  t_chg_totc = t_chg_totc + h_tr(k) * (tf(k) - t0(k))
778  enddo
779  do k=2,nz
780  pe_chg_tot2c = pe_chg_tot2c + (pe_chg_k(k,3) - colht_cor_k(k,3))
781  enddo
782  endif
783 
784  endif
785 
786  if (halves) then
787 
788  ! Set up values appropriate for no diffusivity.
789  do k=1,nz
790  hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
791  dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
792  dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
793  dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
794  dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
795  enddo
796  do k=1,nz+1
797  kd_so_far(k) = 0.0
798  enddo
799 
800  ! Calculate the dependencies on layers above.
801  kddt_h_a(1) = 0.0
802  do k=2,nz ! Loop over interior interfaces.
803  ! First calculate some terms that are independent of the change in Kddt_h(K).
804  kd0 = kd_so_far(k)
805  if (k<=2) then
806  th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
807  else
808  th_a(k-1) = h_tr(k-1) * t0(k-1) + kd_so_far(k-1) * te(k-2)
809  sh_a(k-1) = h_tr(k-1) * s0(k-1) + kd_so_far(k-1) * se(k-2)
810  endif
811  th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
812 
813  dkd = 0.5 * kddt_h(k) - kd_so_far(k)
814 
815  call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
816  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
817  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
818  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
819  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
820  pe_chg=pe_change, colht_cor=colht_cor)
821 
822  pe_chg_k(k,4) = pe_change
823  colht_cor_k(k,4) = colht_cor
824 
825  kd_so_far(k) = kd_so_far(k) + dkd
826 
827  b1 = 1.0 / (hp_a(k-1) + kd_so_far(k))
828  c1_a(k) = kd_so_far(k) * b1
829  if (k==2) then
830  te(1) = b1*(h_tr(1)*t0(1))
831  se(1) = b1*(h_tr(1)*s0(1))
832  else
833  te(k-1) = b1 * (h_tr(k-1) * t0(k-1) + kd_so_far(k-1) * te(k-2))
834  se(k-1) = b1 * (h_tr(k-1) * s0(k-1) + kd_so_far(k-1) * se(k-2))
835  endif
836 
837  hp_a(k) = h_tr(k) + (hp_a(k-1) * b1) * kd_so_far(k)
838  dt_to_dpe_a(k) = dt_to_dpe(k) + c1_a(k)*dt_to_dpe_a(k-1)
839  ds_to_dpe_a(k) = ds_to_dpe(k) + c1_a(k)*ds_to_dpe_a(k-1)
840  dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1_a(k)*dt_to_dcolht_a(k-1)
841  ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1_a(k)*ds_to_dcolht_a(k-1)
842  enddo
843 
844  ! Calculate the dependencies on layers below.
845  do k=nz,2,-1 ! Loop over interior interfaces.
846  ! First calculate some terms that are independent of the change in Kddt_h(K).
847  kd0 = kd_so_far(k)
848  if (k<=2) then
849  th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
850  else
851  th_a(k-1) = h_tr(k-1) * t0(k-1) + kd_so_far(k-1) * te(k-2)
852  sh_a(k-1) = h_tr(k-1) * s0(k-1) + kd_so_far(k-1) * se(k-2)
853  endif
854  if (k>=nz) then
855  th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
856  else
857  th_b(k) = h_tr(k) * t0(k) + kd_so_far(k+1) * te(k+1)
858  sh_b(k) = h_tr(k) * s0(k) + kd_so_far(k+1) * se(k+1)
859  endif
860 
861  dkd = kddt_h(k) - kd_so_far(k)
862 
863  call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
864  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
865  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
866  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
867  dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
868  pe_chg=pe_change, colht_cor=colht_cor)
869 
870  pe_chg_k(k,4) = pe_chg_k(k,4) + pe_change
871  colht_cor_k(k,4) = colht_cor_k(k,4) + colht_cor
872 
873 
874  kd_so_far(k) = kd_so_far(k) + dkd
875 
876  b1 = 1.0 / (hp_b(k) + kd_so_far(k))
877  c1_b(k) = kd_so_far(k) * b1
878  if (k==nz) then
879  te(k) = b1 * (h_tr(k)*t0(k))
880  se(k) = b1 * (h_tr(k)*s0(k))
881  else
882  te(k) = b1 * (h_tr(k) * t0(k) + kd_so_far(k+1) * te(k+1))
883  se(k) = b1 * (h_tr(k) * s0(k) + kd_so_far(k+1) * se(k+1))
884  endif
885 
886  hp_b(k-1) = h_tr(k-1) + (hp_b(k) * b1) * kd_so_far(k)
887  dt_to_dpe_b(k-1) = dt_to_dpe(k-1) + c1_b(k)*dt_to_dpe_b(k)
888  ds_to_dpe_b(k-1) = ds_to_dpe(k-1) + c1_b(k)*ds_to_dpe_b(k)
889  dt_to_dcolht_b(k-1) = dt_to_dcolht(k-1) + c1_b(k)*dt_to_dcolht_b(k)
890  ds_to_dcolht_b(k-1) = ds_to_dcolht(k-1) + c1_b(k)*ds_to_dcolht_b(k)
891 
892  enddo
893 
894  ! Now update the other layer working down from the top to determine the
895  ! final temperatures and salinities.
896  b1 = 1.0 / (hp_b(1))
897  tf(1) = b1 * (h_tr(1) * t0(1) + kddt_h(2) * te(2))
898  sf(1) = b1 * (h_tr(1) * s0(1) + kddt_h(2) * se(2))
899  do k=2,nz
900  tf(k) = te(k) + c1_b(k)*tf(k-1)
901  sf(k) = se(k) + c1_b(k)*sf(k-1)
902  enddo
903 
904  if (debug) then
905  pe_chg_tot1d = 0.0 ; pe_chg_tot2d = 0.0 ; t_chg_totd = 0.0
906  do k=1,nz
907  pe_chg_tot1d = pe_chg_tot1d + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
908  ds_to_dpe(k) * (sf(k) - s0(k)))
909  t_chg_totd = t_chg_totd + h_tr(k) * (tf(k) - t0(k))
910  enddo
911  do k=2,nz
912  pe_chg_tot2d = pe_chg_tot2d + (pe_chg_k(k,4) - colht_cor_k(k,4))
913  enddo
914  endif
915 
916  endif
917 
918  energy_kd = 0.0 ; do k=2,nz ; energy_kd = energy_kd + pe_chg_k(k,1) ; enddo
919  energy_kd = energy_kd / dt
920  k=nz
921 
922  if (do_print) then
923  if (cs%id_ERt>0) call post_data(cs%id_ERt, pe_chg_k(:,1), cs%diag)
924  if (cs%id_ERb>0) call post_data(cs%id_ERb, pe_chg_k(:,2), cs%diag)
925  if (cs%id_ERc>0) call post_data(cs%id_ERc, pe_chg_k(:,3), cs%diag)
926  if (cs%id_ERh>0) call post_data(cs%id_ERh, pe_chg_k(:,4), cs%diag)
927  if (cs%id_Kddt>0) call post_data(cs%id_Kddt, kddt_h, cs%diag)
928  if (cs%id_Kd>0) call post_data(cs%id_Kd, kd, cs%diag)
929  if (cs%id_h>0) call post_data(cs%id_h, h_tr, cs%diag)
930  if (cs%id_zInt>0) call post_data(cs%id_zInt, z_int, cs%diag)
931  if (cs%id_CHCt>0) call post_data(cs%id_CHCt, colht_cor_k(:,1), cs%diag)
932  if (cs%id_CHCb>0) call post_data(cs%id_CHCb, colht_cor_k(:,2), cs%diag)
933  if (cs%id_CHCc>0) call post_data(cs%id_CHCc, colht_cor_k(:,3), cs%diag)
934  if (cs%id_CHCh>0) call post_data(cs%id_CHCh, colht_cor_k(:,4), cs%diag)
935  if (cs%id_T0>0) call post_data(cs%id_T0, t0, cs%diag)
936  if (cs%id_Tf>0) call post_data(cs%id_Tf, tf, cs%diag)
937  if (cs%id_S0>0) call post_data(cs%id_S0, s0, cs%diag)
938  if (cs%id_Sf>0) call post_data(cs%id_Sf, sf, cs%diag)
939  if (cs%id_N2_0>0) then
940  n2(1) = 0.0 ; n2(nz+1) = 0.0
941  do k=2,nz
942  call calculate_density(0.5*(t0(k-1) + t0(k)), 0.5*(s0(k-1) + s0(k)), &
943  pres(k), rho_here, tv%eqn_of_state)
944  n2(k) = ((us%L_to_Z**2*gv%g_Earth) * rho_here / (0.5*gv%H_to_Z*(h_tr(k-1) + h_tr(k)))) * &
945  ( 0.5*(dsv_dt(k-1) + dsv_dt(k)) * (t0(k-1) - t0(k)) + &
946  0.5*(dsv_ds(k-1) + dsv_ds(k)) * (s0(k-1) - s0(k)) )
947  enddo
948  call post_data(cs%id_N2_0, n2, cs%diag)
949  endif
950  if (cs%id_N2_f>0) then
951  n2(1) = 0.0 ; n2(nz+1) = 0.0
952  do k=2,nz
953  call calculate_density(0.5*(tf(k-1) + tf(k)), 0.5*(sf(k-1) + sf(k)), &
954  pres(k), rho_here, tv%eqn_of_state)
955  n2(k) = ((us%L_to_Z**2*gv%g_Earth) * rho_here / (0.5*gv%H_to_Z*(h_tr(k-1) + h_tr(k)))) * &
956  ( 0.5*(dsv_dt(k-1) + dsv_dt(k)) * (tf(k-1) - tf(k)) + &
957  0.5*(dsv_ds(k-1) + dsv_ds(k)) * (sf(k-1) - sf(k)) )
958  enddo
959  call post_data(cs%id_N2_f, n2, cs%diag)
960  endif
961  endif
962 
963 end subroutine diapyc_energy_req_calc
964 
965 !> This subroutine calculates the change in potential energy and or derivatives
966 !! for several changes in an interfaces's diapycnal diffusivity times a timestep.
967 subroutine find_pe_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, &
968  dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, &
969  pres_Z, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, &
970  PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0, ColHt_cor)
971  real, intent(in) :: Kddt_h0 !< The previously used diffusivity at an interface times
972  !! the time step and divided by the average of the
973  !! thicknesses around the interface [H ~> m or kg m-2].
974  real, intent(in) :: dKddt_h !< The trial change in the diffusivity at an interface times
975  !! the time step and divided by the average of the
976  !! thicknesses around the interface [H ~> m or kg m-2].
977  real, intent(in) :: hp_a !< The effective pivot thickness of the layer above the
978  !! interface, given by h_k plus a term that
979  !! is a fraction (determined from the tridiagonal solver) of
980  !! Kddt_h for the interface above [H ~> m or kg m-2].
981  real, intent(in) :: hp_b !< The effective pivot thickness of the layer below the
982  !! interface, given by h_k plus a term that
983  !! is a fraction (determined from the tridiagonal solver) of
984  !! Kddt_h for the interface above [H ~> m or kg m-2].
985  real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer
986  !! above, including implicit mixing effects with other
987  !! yet higher layers [degC H ~> degC m or degC kg m-2].
988  real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer
989  !! above, including implicit mixing effects with other
990  !! yet higher layers [ppt H ~> ppt m or ppt kg m-2].
991  real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer
992  !! below, including implicit mixing effects with other
993  !! yet lower layers [degC H ~> degC m or degC kg m-2].
994  real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer
995  !! below, including implicit mixing effects with other
996  !! yet lower layers [ppt H ~> ppt m or ppt kg m-2].
997  real, intent(in) :: dT_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
998  !! a layer's temperature change to the change in column
999  !! potential energy, including all implicit diffusive changes
1000  !! in the temperatures of all the layers above [J m-2 degC-1].
1001  real, intent(in) :: dS_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
1002  !! a layer's salinity change to the change in column
1003  !! potential energy, including all implicit diffusive changes
1004  !! in the salinities of all the layers above [J m-2 ppt-1].
1005  real, intent(in) :: dT_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
1006  !! a layer's temperature change to the change in column
1007  !! potential energy, including all implicit diffusive changes
1008  !! in the temperatures of all the layers below [J m-2 degC-1].
1009  real, intent(in) :: dS_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
1010  !! a layer's salinity change to the change in column
1011  !! potential energy, including all implicit diffusive changes
1012  !! in the salinities of all the layers below [J m-2 ppt-1].
1013  real, intent(in) :: pres_Z !< The hydrostatic interface pressure, which is used to relate
1014  !! the changes in column thickness to the energy that is radiated
1015  !! as gravity waves and unavailable to drive mixing [J m-2 Z-1 ~> J m-3].
1016  real, intent(in) :: dT_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dT) relating
1017  !! a layer's temperature change to the change in column
1018  !! height, including all implicit diffusive changes
1019  !! in the temperatures of all the layers above [Z degC-1 ~> m degC-1].
1020  real, intent(in) :: dS_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dS) relating
1021  !! a layer's salinity change to the change in column
1022  !! height, including all implicit diffusive changes
1023  !! in the salinities of all the layers above [Z ppt-1 ~> m ppt-1].
1024  real, intent(in) :: dT_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dT) relating
1025  !! a layer's temperature change to the change in column
1026  !! height, including all implicit diffusive changes
1027  !! in the temperatures of all the layers below [Z degC-1 ~> m degC-1].
1028  real, intent(in) :: dS_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dS) relating
1029  !! a layer's salinity change to the change in column
1030  !! height, including all implicit diffusive changes
1031  !! in the salinities of all the layers below [Z ppt-1 ~> m ppt-1].
1032 
1033  real, optional, intent(out) :: PE_chg !< The change in column potential energy from applying
1034  !! Kddt_h at the present interface [J m-2].
1035  real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h,
1036  !! [J m-2 H-1 ~> J m-3 or J kg-1].
1037  real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could
1038  !! be realizedd by applying a huge value of Kddt_h at the
1039  !! present interface [J m-2].
1040  real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the
1041  !! limit where Kddt_h = 0 [J m-2 H-1 ~> J m-3 or J kg-1].
1042  real, optional, intent(out) :: ColHt_cor !< The correction to PE_chg that is made due to a net
1043  !! change in the column height [J m-2].
1044 
1045  real :: hps ! The sum of the two effective pivot thicknesses [H ~> m or kg m-2].
1046  real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term [H2 ~> m2 or kg2 m-4].
1047  real :: dT_c ! The core term in the expressions for the temperature changes [degC H2 ~> degC m2 or degC kg2 m-4].
1048  real :: dS_c ! The core term in the expressions for the salinity changes [psu H2 ~> psu m2 or psu kg2 m-4].
1049  real :: PEc_core ! The diffusivity-independent core term in the expressions
1050  ! for the potential energy changes [J m-3].
1051  real :: ColHt_core ! The diffusivity-independent core term in the expressions
1052  ! for the column height changes [J m-3].
1053  real :: ColHt_chg ! The change in the column height [Z ~> m].
1054  real :: y1 ! A local temporary term, in [H-3] or [H-4] in various contexts.
1055 
1056  ! The expression for the change in potential energy used here is derived
1057  ! from the expression for the final estimates of the changes in temperature
1058  ! and salinities, and then extensively manipulated to get it into its most
1059  ! succint form. The derivation is not necessarily obvious, but it demonstrably
1060  ! works by comparison with separate calculations of the energy changes after
1061  ! the tridiagonal solver for the final changes in temperature and salinity are
1062  ! applied.
1063 
1064  hps = hp_a + hp_b
1065  bdt1 = hp_a * hp_b + kddt_h0 * hps
1066  dt_c = hp_a * th_b - hp_b * th_a
1067  ds_c = hp_a * sh_b - hp_b * sh_a
1068  pec_core = hp_b * (dt_to_dpe_a * dt_c + ds_to_dpe_a * ds_c) - &
1069  hp_a * (dt_to_dpe_b * dt_c + ds_to_dpe_b * ds_c)
1070  colht_core = hp_b * (dt_to_dcolht_a * dt_c + ds_to_dcolht_a * ds_c) - &
1071  hp_a * (dt_to_dcolht_b * dt_c + ds_to_dcolht_b * ds_c)
1072 
1073  if (present(pe_chg)) then
1074  ! Find the change in column potential energy due to the change in the
1075  ! diffusivity at this interface by dKddt_h.
1076  y1 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
1077  pe_chg = pec_core * y1
1078  colht_chg = colht_core * y1
1079  if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
1080  if (present(colht_cor)) colht_cor = -pres_z * min(colht_chg, 0.0)
1081  elseif (present(colht_cor)) then
1082  y1 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
1083  colht_cor = -pres_z * min(colht_core * y1, 0.0)
1084  endif
1085 
1086  if (present(dpec_dkd)) then
1087  ! Find the derivative of the potential energy change with dKddt_h.
1088  y1 = 1.0 / (bdt1 + dkddt_h * hps)**2
1089  dpec_dkd = pec_core * y1
1090  colht_chg = colht_core * y1
1091  if (colht_chg < 0.0) dpec_dkd = dpec_dkd - pres_z * colht_chg
1092  endif
1093 
1094  if (present(dpe_max)) then
1095  ! This expression is the limit of PE_chg for infinite dKddt_h.
1096  y1 = 1.0 / (bdt1 * hps)
1097  dpe_max = pec_core * y1
1098  colht_chg = colht_core * y1
1099  if (colht_chg < 0.0) dpe_max = dpe_max - pres_z * colht_chg
1100  endif
1101 
1102  if (present(dpec_dkd_0)) then
1103  ! This expression is the limit of dPEc_dKd for dKddt_h = 0.
1104  y1 = 1.0 / bdt1**2
1105  dpec_dkd_0 = pec_core * y1
1106  colht_chg = colht_core * y1
1107  if (colht_chg < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z * colht_chg
1108  endif
1109 
1110 end subroutine find_pe_chg
1111 
1112 
1113 !> This subroutine calculates the change in potential energy and or derivatives
1114 !! for several changes in an interfaces's diapycnal diffusivity times a timestep
1115 !! using the original form used in the first version of ePBL.
1116 subroutine find_pe_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, &
1117  dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, &
1118  dT_to_dPEa, dS_to_dPEa, pres_Z, dT_to_dColHt_k, &
1119  dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, &
1120  PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0)
1121  real, intent(in) :: Kddt_h !< The diffusivity at an interface times the time step and
1122  !! divided by the average of the thicknesses around the
1123  !! interface [H ~> m or kg m-2].
1124  real, intent(in) :: h_k !< The thickness of the layer below the interface [H ~> m or kg m-2].
1125  real, intent(in) :: b_den_1 !< The first term in the denominator of the pivot
1126  !! for the tridiagonal solver, given by h_k plus a term that
1127  !! is a fraction (determined from the tridiagonal solver) of
1128  !! Kddt_h for the interface above [H ~> m or kg m-2].
1129  real, intent(in) :: dTe_term !< A diffusivity-independent term related to the temperature change
1130  !! in the layer below the interface [degC H ~> degC m or degC kg m-2].
1131  real, intent(in) :: dSe_term !< A diffusivity-independent term related to the salinity change
1132  !! in the layer below the interface [ppt H ~> ppt m or ppt kg m-2].
1133  real, intent(in) :: dT_km1_t2 !< A diffusivity-independent term related to the
1134  !! temperature change in the layer above the interface [degC].
1135  real, intent(in) :: dS_km1_t2 !< A diffusivity-independent term related to the
1136  !! salinity change in the layer above the interface [ppt].
1137  real, intent(in) :: pres_Z !< The hydrostatic interface pressure, which is used to relate
1138  !! the changes in column thickness to the energy that is radiated
1139  !! as gravity waves and unavailable to drive mixing [J m-2 Z-1 ~> J m-3].
1140  real, intent(in) :: dT_to_dPE_k !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
1141  !! a layer's temperature change to the change in column
1142  !! potential energy, including all implicit diffusive changes
1143  !! in the temperatures of all the layers below [J m-2 degC-1].
1144  real, intent(in) :: dS_to_dPE_k !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
1145  !! a layer's salinity change to the change in column
1146  !! potential energy, including all implicit diffusive changes
1147  !! in the salinities of all the layers below [J m-2 ppt-1].
1148  real, intent(in) :: dT_to_dPEa !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
1149  !! a layer's temperature change to the change in column
1150  !! potential energy, including all implicit diffusive changes
1151  !! in the temperatures of all the layers above [J m-2 degC-1].
1152  real, intent(in) :: dS_to_dPEa !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
1153  !! a layer's salinity change to the change in column
1154  !! potential energy, including all implicit diffusive changes
1155  !! in the salinities of all the layers above [J m-2 ppt-1].
1156  real, intent(in) :: dT_to_dColHt_k !< A factor (mass_lay*dSColHtc_vol/dT) relating
1157  !! a layer's temperature change to the change in column
1158  !! height, including all implicit diffusive changes
1159  !! in the temperatures of all the layers below [Z degC-1 ~> m degC-1].
1160  real, intent(in) :: dS_to_dColHt_k !< A factor (mass_lay*dSColHtc_vol/dS) relating
1161  !! a layer's salinity change to the change in column
1162  !! height, including all implicit diffusive changes
1163  !! in the salinities of all the layers below [Z ppt-1 ~> m ppt-1].
1164  real, intent(in) :: dT_to_dColHta !< A factor (mass_lay*dSColHtc_vol/dT) relating
1165  !! a layer's temperature change to the change in column
1166  !! height, including all implicit diffusive changes
1167  !! in the temperatures of all the layers above [Z degC-1 ~> m degC-1].
1168  real, intent(in) :: dS_to_dColHta !< A factor (mass_lay*dSColHtc_vol/dS) relating
1169  !! a layer's salinity change to the change in column
1170  !! height, including all implicit diffusive changes
1171  !! in the salinities of all the layers above [Z ppt-1 ~> m ppt-1].
1172 
1173  real, optional, intent(out) :: PE_chg !< The change in column potential energy from applying
1174  !! Kddt_h at the present interface [J m-2].
1175  real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h,
1176  !! [J m-2 H-1 ~> J m-3 or J kg-1].
1177  real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could
1178  !! be realized by applying a huge value of Kddt_h at the
1179  !! present interface [J m-2].
1180  real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the
1181  !! limit where Kddt_h = 0 [J m-2 H-1 ~> J m-3 or J kg-1].
1182 
1183 ! This subroutine determines the total potential energy change due to mixing
1184 ! at an interface, including all of the implicit effects of the prescribed
1185 ! mixing at interfaces above. Everything here is derived by careful manipulation
1186 ! of the robust tridiagonal solvers used for tracers by MOM6. The results are
1187 ! positive for mixing in a stably stratified environment.
1188 ! The comments describing these arguments are for a downward mixing pass, but
1189 ! this routine can also be used for an upward pass with the sense of direction
1190 ! reversed.
1191 
1192  real :: b1 ! b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
1193  real :: b1Kd ! Temporary array [nondim]
1194  real :: ColHt_chg ! The change in column thickness [Z ~> m].
1195  real :: dColHt_max ! The change in column thickness for infinite diffusivity [Z ~> m].
1196  real :: dColHt_dKd ! The partial derivative of column thickness with diffusivity [s Z-1 ~> s m-1].
1197  real :: dT_k, dT_km1 ! Temporary arrays [degC].
1198  real :: dS_k, dS_km1 ! Temporary arrays [ppt].
1199  real :: I_Kr_denom ! Temporary arrays [H-2 ~> m-2 or m4 kg-2].
1200  real :: dKr_dKd ! Nondimensional temporary array [nondim].
1201  real :: ddT_k_dKd, ddT_km1_dKd ! Temporary arrays [degC H-1 ~> m-1 or m2 kg-1].
1202  real :: ddS_k_dKd, ddS_km1_dKd ! Temporary arrays [ppt H-1 ~> ppt m-1 or ppt m2 kg-1].
1203 
1204  b1 = 1.0 / (b_den_1 + kddt_h)
1205  b1kd = kddt_h*b1
1206 
1207  ! Start with the temperature change in layer k-1 due to the diffusivity at
1208  ! interface K without considering the effects of changes in layer k.
1209 
1210  ! Calculate the change in PE due to the diffusion at interface K
1211  ! if Kddt_h(K+1) = 0.
1212  i_kr_denom = 1.0 / (h_k*b_den_1 + (b_den_1 + h_k)*kddt_h)
1213 
1214  dt_k = (kddt_h*i_kr_denom) * dte_term
1215  ds_k = (kddt_h*i_kr_denom) * dse_term
1216 
1217  if (present(pe_chg)) then
1218  ! Find the change in energy due to diffusion with strength Kddt_h at this interface.
1219  ! Increment the temperature changes in layer k-1 due the changes in layer k.
1220  dt_km1 = b1kd * ( dt_k + dt_km1_t2 )
1221  ds_km1 = b1kd * ( ds_k + ds_km1_t2 )
1222 
1223  pe_chg = (dt_to_dpe_k * dt_k + dt_to_dpea * dt_km1) + &
1224  (ds_to_dpe_k * ds_k + ds_to_dpea * ds_km1)
1225  colht_chg = (dt_to_dcolht_k * dt_k + dt_to_dcolhta * dt_km1) + &
1226  (ds_to_dcolht_k * ds_k + ds_to_dcolhta * ds_km1)
1227  if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
1228  endif
1229 
1230  if (present(dpec_dkd)) then
1231  ! Find the derivatives of the temperature and salinity changes with Kddt_h.
1232  dkr_dkd = (h_k*b_den_1) * i_kr_denom**2
1233 
1234  ddt_k_dkd = dkr_dkd * dte_term
1235  dds_k_dkd = dkr_dkd * dse_term
1236  ddt_km1_dkd = (b1**2 * b_den_1) * ( dt_k + dt_km1_t2 ) + b1kd * ddt_k_dkd
1237  dds_km1_dkd = (b1**2 * b_den_1) * ( ds_k + ds_km1_t2 ) + b1kd * dds_k_dkd
1238 
1239  ! Calculate the partial derivative of Pe_chg with Kddt_h.
1240  dpec_dkd = (dt_to_dpe_k * ddt_k_dkd + dt_to_dpea * ddt_km1_dkd) + &
1241  (ds_to_dpe_k * dds_k_dkd + ds_to_dpea * dds_km1_dkd)
1242  dcolht_dkd = (dt_to_dcolht_k * ddt_k_dkd + dt_to_dcolhta * ddt_km1_dkd) + &
1243  (ds_to_dcolht_k * dds_k_dkd + ds_to_dcolhta * dds_km1_dkd)
1244  if (dcolht_dkd < 0.0) dpec_dkd = dpec_dkd - pres_z * dcolht_dkd
1245  endif
1246 
1247  if (present(dpe_max)) then
1248  ! This expression is the limit of PE_chg for infinite Kddt_h.
1249  dpe_max = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) + &
1250  ((dt_to_dpe_k + dt_to_dpea) * dte_term + &
1251  (ds_to_dpe_k + ds_to_dpea) * dse_term) / (b_den_1 + h_k)
1252  dcolht_max = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) + &
1253  ((dt_to_dcolht_k + dt_to_dcolhta) * dte_term + &
1254  (ds_to_dcolht_k + ds_to_dcolhta) * dse_term) / (b_den_1 + h_k)
1255  if (dcolht_max < 0.0) dpe_max = dpe_max - pres_z*dcolht_max
1256  endif
1257 
1258  if (present(dpec_dkd_0)) then
1259  ! This expression is the limit of dPEc_dKd for Kddt_h = 0.
1260  dpec_dkd_0 = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) / (b_den_1) + &
1261  (dt_to_dpe_k * dte_term + ds_to_dpe_k * dse_term) / (h_k*b_den_1)
1262  dcolht_dkd = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) / (b_den_1) + &
1263  (dt_to_dcolht_k * dte_term + ds_to_dcolht_k * dse_term) / (h_k*b_den_1)
1264  if (dcolht_dkd < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z*dcolht_dkd
1265  endif
1266 
1267 end subroutine find_pe_chg_orig
1268 
1269 !> Initialize parameters and allocate memory associated with the diapycnal energy requirement module.
1270 subroutine diapyc_energy_req_init(Time, G, GV, US, param_file, diag, CS)
1271  type(time_type), intent(in) :: time !< model time
1272  type(ocean_grid_type), intent(in) :: g !< model grid structure
1273  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1274  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1275  type(param_file_type), intent(in) :: param_file !< file to parse for parameter values
1276  type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output
1277  type(diapyc_energy_req_cs), pointer :: cs !< module control structure
1278 
1279  integer, save :: init_calls = 0
1280 ! This include declares and sets the variable "version".
1281 #include "version_variable.h"
1282  character(len=40) :: mdl = "MOM_diapyc_energy_req" ! This module's name.
1283  character(len=256) :: mesg ! Message for error messages.
1284 
1285  if (.not.associated(cs)) then ; allocate(cs)
1286  else ; return ; endif
1287 
1288  cs%initialized = .true.
1289  cs%diag => diag
1290 
1291  ! Read all relevant parameters and write them to the model log.
1292  call log_version(param_file, mdl, version, "")
1293  call get_param(param_file, mdl, "ENERGY_REQ_KH_SCALING", cs%test_Kh_scaling, &
1294  "A scaling factor for the diapycnal diffusivity used in "//&
1295  "testing the energy requirements.", default=1.0, units="nondim")
1296  call get_param(param_file, mdl, "ENERGY_REQ_COL_HT_SCALING", cs%ColHt_scaling, &
1297  "A scaling factor for the column height change correction "//&
1298  "used in testing the energy requirements.", default=1.0, units="nondim")
1299  call get_param(param_file, mdl, "ENERGY_REQ_USE_TEST_PROFILE", &
1300  cs%use_test_Kh_profile, &
1301  "If true, use the internal test diffusivity profile in "//&
1302  "place of any that might be passed in as an argument.", default=.false.)
1303 
1304  cs%id_ERt = register_diag_field('ocean_model', 'EnReqTest_ERt', diag%axesZi, time, &
1305  "Diffusivity Energy Requirements, top-down", "J m-2")
1306  cs%id_ERb = register_diag_field('ocean_model', 'EnReqTest_ERb', diag%axesZi, time, &
1307  "Diffusivity Energy Requirements, bottom-up", "J m-2")
1308  cs%id_ERc = register_diag_field('ocean_model', 'EnReqTest_ERc', diag%axesZi, time, &
1309  "Diffusivity Energy Requirements, center-last", "J m-2")
1310  cs%id_ERh = register_diag_field('ocean_model', 'EnReqTest_ERh', diag%axesZi, time, &
1311  "Diffusivity Energy Requirements, halves", "J m-2")
1312  cs%id_Kddt = register_diag_field('ocean_model', 'EnReqTest_Kddt', diag%axesZi, time, &
1313  "Implicit diffusive coupling coefficient", "m", conversion=gv%H_to_m)
1314  cs%id_Kd = register_diag_field('ocean_model', 'EnReqTest_Kd', diag%axesZi, time, &
1315  "Diffusivity in test", "m2 s-1", conversion=us%Z_to_m**2)
1316  cs%id_h = register_diag_field('ocean_model', 'EnReqTest_h_lay', diag%axesZL, time, &
1317  "Test column layer thicknesses", "m", conversion=gv%H_to_m)
1318  cs%id_zInt = register_diag_field('ocean_model', 'EnReqTest_z_int', diag%axesZi, time, &
1319  "Test column layer interface heights", "m", conversion=gv%H_to_m)
1320  cs%id_CHCt = register_diag_field('ocean_model', 'EnReqTest_CHCt', diag%axesZi, time, &
1321  "Column Height Correction to Energy Requirements, top-down", "J m-2")
1322  cs%id_CHCb = register_diag_field('ocean_model', 'EnReqTest_CHCb', diag%axesZi, time, &
1323  "Column Height Correction to Energy Requirements, bottom-up", "J m-2")
1324  cs%id_CHCc = register_diag_field('ocean_model', 'EnReqTest_CHCc', diag%axesZi, time, &
1325  "Column Height Correction to Energy Requirements, center-last", "J m-2")
1326  cs%id_CHCh = register_diag_field('ocean_model', 'EnReqTest_CHCh', diag%axesZi, time, &
1327  "Column Height Correction to Energy Requirements, halves", "J m-2")
1328  cs%id_T0 = register_diag_field('ocean_model', 'EnReqTest_T0', diag%axesZL, time, &
1329  "Temperature before mixing", "deg C")
1330  cs%id_Tf = register_diag_field('ocean_model', 'EnReqTest_Tf', diag%axesZL, time, &
1331  "Temperature after mixing", "deg C")
1332  cs%id_S0 = register_diag_field('ocean_model', 'EnReqTest_S0', diag%axesZL, time, &
1333  "Salinity before mixing", "g kg-1")
1334  cs%id_Sf = register_diag_field('ocean_model', 'EnReqTest_Sf', diag%axesZL, time, &
1335  "Salinity after mixing", "g kg-1")
1336  cs%id_N2_0 = register_diag_field('ocean_model', 'EnReqTest_N2_0', diag%axesZi, time, &
1337  "Squared buoyancy frequency before mixing", "second-2", conversion=us%s_to_T**2)
1338  cs%id_N2_f = register_diag_field('ocean_model', 'EnReqTest_N2_f', diag%axesZi, time, &
1339  "Squared buoyancy frequency after mixing", "second-2", conversion=us%s_to_T**2)
1340 
1341 end subroutine diapyc_energy_req_init
1342 
1343 !> Clean up and deallocate memory associated with the diapycnal energy requirement module.
1344 subroutine diapyc_energy_req_end(CS)
1345  type(diapyc_energy_req_cs), pointer :: cs !< Diapycnal energy requirement control structure that
1346  !! will be deallocated in this subroutine.
1347  if (associated(cs)) deallocate(cs)
1348 end subroutine diapyc_energy_req_end
1349 
1350 end module mom_diapyc_energy_req
mom_diapyc_energy_req::diapyc_energy_req_calc
subroutine, public diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV, US, may_print, CS)
This subroutine uses a substantially refactored tridiagonal equation for diapycnal mixing of temperat...
Definition: MOM_diapyc_energy_req.F90:122
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_diapyc_energy_req::find_pe_chg_orig
subroutine find_pe_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, dT_to_dPEa, dS_to_dPEa, pres_Z, dT_to_dColHt_k, dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0)
This subroutine calculates the change in potential energy and or derivatives for several changes in a...
Definition: MOM_diapyc_energy_req.F90:1121
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_diapyc_energy_req::find_pe_chg
subroutine find_pe_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, pres_Z, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0, ColHt_cor)
This subroutine calculates the change in potential energy and or derivatives for several changes in a...
Definition: MOM_diapyc_energy_req.F90:971
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_diapyc_energy_req::diapyc_energy_req_init
subroutine, public diapyc_energy_req_init(Time, G, GV, US, param_file, diag, CS)
Initialize parameters and allocate memory associated with the diapycnal energy requirement module.
Definition: MOM_diapyc_energy_req.F90:1271
mom_diapyc_energy_req::diapyc_energy_req_test
subroutine, public diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int)
This subroutine helps test the accuracy of the diapycnal mixing energy requirement code by writing di...
Definition: MOM_diapyc_energy_req.F90:50
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_diapyc_energy_req::diapyc_energy_req_end
subroutine, public diapyc_energy_req_end(CS)
Clean up and deallocate memory associated with the diapycnal energy requirement module.
Definition: MOM_diapyc_energy_req.F90:1345
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_diapyc_energy_req::diapyc_energy_req_cs
This control structure holds parameters for the MOM_diapyc_energy_req module.
Definition: MOM_diapyc_energy_req.F90:27
mom_eos::calculate_specific_vol_derivs
subroutine, public calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS, scale)
Calls the appropriate subroutine to calculate specific volume derivatives for an array.
Definition: MOM_EOS.F90:546
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_diapyc_energy_req
Calculates the energy requirements of mixing.
Definition: MOM_diapyc_energy_req.F90:2
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60