MOM6
Data Types | Functions/Subroutines | Variables
mom_kappa_shear Module Reference

Data Types

type  kappa_shear_cs
 

Functions/Subroutines

subroutine, public calculate_kappa_shear (u_in, v_in, h, tv, p_surf, kappa_io, tke_io, kv_io, dt, G, GV, CS, initialize_all)
 Subroutine for calculating diffusivity and TKE. More...
 
subroutine calculate_projected_state (kappa, u0, v0, T0, S0, dt, nz, dz, I_dz_int, dbuoy_dT, dbuoy_dS, u, v, T, Sal, N2, S2, ks_int, ke_int)
 
subroutine find_kappa_tke (N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, nz, CS, K_Q, tke, kappa, kappa_src, local_src)
 This subroutine calculates new, consistent estimates of TKE and kappa. More...
 
logical function, public kappa_shear_init (Time, G, GV, param_file, diag, CS)
 
logical function, public kappa_shear_is_used (param_file)
 

Variables

character(len=40) mdl = "MOM_kappa_shear"
 

Function/Subroutine Documentation

◆ calculate_kappa_shear()

subroutine, public mom_kappa_shear::calculate_kappa_shear ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  u_in,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  v_in,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension(:,:), pointer  p_surf,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout)  kappa_io,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout)  tke_io,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout)  kv_io,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(kappa_shear_cs), pointer  CS,
logical, intent(in), optional  initialize_all 
)

Subroutine for calculating diffusivity and TKE.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]u_inInitial zonal velocity, in m s-1. (Intent in)
[in]v_inInitial meridional velocity, in m s-1.
[in]hLayer thicknesses, in H (usually m or kg m-2).
[in]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
p_surfThe pressure at the ocean surface in Pa (or NULL).
[in,out]kappa_ioThe diapycnal diffusivity at each interface
[in,out]tke_ioThe turbulent kinetic energy per unit mass at
[in,out]kv_ioThe vertical viscosity at each interface
[in]dtTime increment, in s.
csThe control structure returned by a previous call to kappa_shear_init.
[in]initialize_allIf present and false, the previous value of kappa is used to start the iterations

Definition at line 122 of file MOM_kappa_shear.F90.

References mom_eos::calculate_density_derivs(), calculate_projected_state(), and find_kappa_tke().

Referenced by mom_set_diffusivity::set_diffusivity().

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  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
125  intent(in) :: u_in !< Initial zonal velocity, in m s-1. (Intent in)
126  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
127  intent(in) :: v_in !< Initial meridional velocity, in m s-1.
128  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
129  intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2).
130  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
131  !! available thermodynamic fields. Absent fields
132  !! have NULL ptrs.
133  real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface in Pa
134  !! (or NULL).
135  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
136  intent(inout) :: kappa_io !< The diapycnal diffusivity at each interface
137  !! (not layer!) in m2 s-1. Initially this is the
138  !! value from the previous timestep, which may
139  !! accelerate the iteration toward convergence.
140  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
141  intent(inout) :: tke_io !< The turbulent kinetic energy per unit mass at
142  !! each interface (not layer!) in m2 s-2.
143  !! Initially this is the value from the previous
144  !! timestep, which may accelerate the iteration
145  !! toward convergence.
146  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
147  intent(inout) :: kv_io !< The vertical viscosity at each interface
148  !! (not layer!) in m2 s-1. This discards any
149  !! previous value i.e. intent(out) and simply
150  !! sets Kv = Prandtl * Kd_turb
151  real, intent(in) :: dt !< Time increment, in s.
152  type(kappa_shear_cs), pointer :: cs !< The control structure returned by a previous
153  !! call to kappa_shear_init.
154  logical, optional, intent(in) :: initialize_all !< If present and false, the previous
155  !! value of kappa is used to start the iterations
156 !
157 ! ----------------------------------------------
158 ! Subroutine for calculating diffusivity and TKE
159 ! ----------------------------------------------
160 ! Arguments: u_in - Initial zonal velocity, in m s-1. (Intent in)
161 ! (in) v_in - Initial meridional velocity, in m s-1.
162 ! (in) h - Layer thickness, in m or kg m-2.
163 ! (in) tv - A structure containing pointers to any available
164 ! thermodynamic fields. Absent fields have NULL ptrs.
165 ! (in) p_surf - The pressure at the ocean surface in Pa (or NULL).
166 ! (in/out) kappa_io - The diapycnal diffusivity at each interface
167 ! (not layer!) in m2 s-1. Initially this is the value
168 ! from the previous timestep, which may accelerate the
169 ! iteration toward convergence.
170 ! (in/out) tke_io - The turbulent kinetic energy per unit mass at each
171 ! interface (not layer!) in m2 s-2. Initially this is the
172 ! value from the previous timestep, which may accelerate
173 ! the iteration toward convergence.
174 ! (in/out) kv_io - The vertical viscosity at each interface
175 ! (not layer!) in m2 s-1. This discards any previous value
176 ! i.e. intent(out) and simply sets Kv = Prandtl * Kd_turb
177 ! (in) dt - Time increment, in s.
178 ! (in) G - The ocean's grid structure.
179 ! (in) GV - The ocean's vertical grid structure.
180 ! (in) CS - The control structure returned by a previous call to
181 ! kappa_shear_init.
182 ! (in,opt) initialize_all - If present and false, the previous value of
183 ! kappa is used to start the iterations.
184  real, dimension(SZI_(G),SZK_(G)) :: &
185  h_2d, & ! A 2-D version of h, but converted to m.
186  u_2d, v_2d, t_2d, s_2d, rho_2d ! 2-D versions of u_in, v_in, T, S, and rho.
187  real, dimension(SZI_(G),SZK_(G)+1) :: &
188  kappa_2d, tke_2d ! 2-D versions of various kappa_io and tke_io.
189  real, dimension(SZK_(G)) :: &
190  u, & ! The zonal velocity after a timestep of mixing, in m s-1.
191  v, & ! The meridional velocity after a timestep of mixing, in m s-1.
192  idz, & ! The inverse of the distance between TKE points, in m.
193  t, & ! The potential temperature after a timestep of mixing, in C.
194  sal, & ! The salinity after a timestep of mixing, in psu.
195  dz, & ! The layer thickness, in m.
196  u0xdz, & ! The initial zonal velocity times dz, in m2 s-1.
197  v0xdz, & ! The initial meridional velocity times dz, in m2 s-1.
198  t0xdz, & ! The initial temperature times dz, in C m.
199  s0xdz, & ! The initial salinity times dz, in PSU m.
200  u_test, v_test, t_test, s_test
201  real, dimension(SZK_(G)+1) :: &
202  n2, & ! The squared buoyancy frequency at an interface, in s-2.
203  dz_int, & ! The extent of a finite-volume space surrounding an interface,
204  ! as used in calculating kappa and TKE, in m.
205  i_dz_int, & ! The inverse of the distance between velocity & density points
206  ! above and below an interface, in m-1. This is used to
207  ! calculate N2, shear, and fluxes, and it might differ from
208  ! 1/dz_Int, as they have different uses.
209  s2, & ! The squared shear at an interface, in s-2.
210  a1, & ! a1 is the coupling between adjacent interfaces in the TKE,
211  ! velocity, and density equations, in m s-1 or m.
212  c1, & ! c1 is used in the tridiagonal (and similar) solvers.
213  k_src, & ! The shear-dependent source term in the kappa equation, in s-1.
214  kappa_src, & ! The shear-dependent source term in the kappa equation in s-1.
215  kappa, & ! The shear-driven diapycnal diffusivity at an interface, in
216  ! units of m2 s-1.
217  tke, & ! The Turbulent Kinetic Energy per unit mass at an interface,
218  ! in units of m2 s-2.
219  kappa_avg, & ! The time-weighted average of kappa, in m2 s-1.
220  tke_avg, & ! The time-weighted average of TKE, in m2 s-2.
221  kappa_out, & ! The kappa that results from the kappa equation, in m2 s-1.
222  kappa_mid, & ! The average of the initial and predictor estimates of kappa,
223  ! in units of m2 s-1.
224  tke_pred, & ! The value of TKE from a predictor step, in m2 s-2.
225  kappa_pred, & ! The value of kappa from a predictor step, in m2 s-1.
226  pressure, & ! The pressure at an interface, in Pa.
227  t_int, & ! The temperature interpolated to an interface, in C.
228  sal_int, & ! The salinity interpolated to an interface, in psu.
229  dbuoy_dt, & ! The partial derivatives of buoyancy with changes in
230  dbuoy_ds, & ! temperature and salinity, in m s-2 K-1 and m s-2 psu-1.
231  i_l2_bdry, & ! The inverse of the square of twice the harmonic mean
232  ! distance to the top and bottom boundaries, in m-2.
233  k_q, & ! Diffusivity divided by TKE, in s.
234  k_q_tmp, & ! Diffusivity divided by TKE, in s.
235  local_src_avg, & ! The time-integral of the local source, nondim.
236  tol_min, & ! Minimum tolerated ksrc for the corrector step, in s-1.
237  tol_max, & ! Maximum tolerated ksrc for the corrector step, in s-1.
238  tol_chg, & ! The tolerated change integrated in time, nondim.
239  dist_from_top, & ! The distance from the top surface, in m.
240  local_src ! The sum of all sources of kappa, including kappa_src and
241  ! sources from the elliptic term, in s-1.
242  real :: f2 ! The squared Coriolis parameter of each column, in s-2.
243  real :: dist_from_bot ! The distance from the bottom surface, in m.
244 
245  real :: b1 ! The inverse of the pivot in the tridiagonal equations.
246  real :: bd1 ! A term in the denominator of b1.
247  real :: d1 ! 1 - c1 in the tridiagonal equations.
248  real :: gr0 ! Rho_0 times g in kg m-2 s-2.
249  real :: g_r0 ! g_R0 is g/Rho in m4 kg-1 s-2.
250  real :: norm ! A factor that normalizes two weights to 1, in m-2.
251  real :: tol_dksrc, tol2 ! ### Tolerances that need to be set better later.
252  real :: tol_dksrc_low ! The tolerance for the fractional decrease in ksrc
253  ! within an iteration. 0 < tol_dksrc_low < 1.
254  real :: ri_crit ! The critical shear Richardson number for shear-
255  ! driven mixing. The theoretical value is 0.25.
256  real :: dt_rem ! The remaining time to advance the solution, in s.
257  real :: dt_now ! The time step used in the current iteration, in s.
258  real :: dt_wt ! The fractional weight of the current iteration, ND.
259  real :: dt_test ! A time-step that is being tested for whether it
260  ! gives acceptably small changes in k_src, in s.
261  real :: idtt ! Idtt = 1 / dt_test, in s-1.
262  real :: dt_inc ! An increment to dt_test that is being tested, in s.
263 
264  real :: dz_in_lay ! The running sum of the thickness in a layer, in m.
265  real :: k0dt ! The background diffusivity times the timestep, in m2.
266  real :: dz_massless ! A layer thickness that is considered massless, in m.
267  logical :: valid_dt ! If true, all levels so far exhibit acceptably small
268  ! changes in k_src.
269  logical :: use_temperature ! If true, temperature and salinity have been
270  ! allocated and are being used as state variables.
271  logical :: new_kappa = .true. ! If true, ignore the value of kappa from the
272  ! last call to this subroutine.
273 
274  integer, dimension(SZK_(G)+1) :: kc ! The index map between the original
275  ! interfaces and the interfaces with massless layers
276  ! merged into nearby massive layers.
277  real, dimension(SZK_(G)+1) :: kf ! The fractional weight of interface kc+1 for
278  ! interpolating back to the original index space, ND.
279  integer :: ks_kappa, ke_kappa ! The k-range with nonzero kappas.
280  integer :: dt_halvings ! The number of times that the time-step is halved
281  ! in seeking an acceptable timestep. If none is
282  ! found, dt_rem*0.5^dt_halvings is used.
283  integer :: dt_refinements ! The number of 2-fold refinements that will be used
284  ! to estimate the maximum permitted time step. I.e.,
285  ! the resolution is 1/2^dt_refinements.
286  integer :: is, ie, js, je, i, j, k, nz, nzc, itt, itt_dt
287 
288  ! Arrays that are used in find_kappa_tke if it is included in this subroutine.
289  real, dimension(SZK_(G)) :: &
290  aq, & ! aQ is the coupling between adjacent interfaces in the TKE
291  ! equations, in m s-1.
292  dqdz ! Half the partial derivative of TKE with depth, m s-2.
293  real, dimension(SZK_(G)+1) :: &
294  dk, & ! The change in kappa, in m2 s-1.
295  dq, & ! The change in TKE, in m2 s-1.
296  cq, ck, & ! cQ and cK are the upward influences in the tridiagonal and
297  ! hexadiagonal solvers for the TKE and kappa equations, ND.
298  i_ld2, & ! 1/Ld^2, where Ld is the effective decay length scale
299  ! for kappa, in units of m-2.
300  tke_decay, & ! The local TKE decay rate in s-1.
301  dqmdk, & ! With Newton's method the change in dQ(K-1) due to dK(K), s.
302  dkdq, & ! With Newton's method the change in dK(K) due to dQ(K), s-1.
303  e1 ! The fractional change in a layer TKE due to a change in the
304  ! TKE of the layer above when all the kappas below are 0.
305  ! e1 is nondimensional, and 0 < e1 < 1.
306 
307 
308  ! Diagnostics that should be deleted?
309 #ifdef ADD_DIAGNOSTICS
310  real, dimension(SZK_(G)+1) :: & ! Additional diagnostics.
311  i_ld2_1d
312  real, dimension(SZI_(G),SZK_(G)+1) :: & ! 2-D versions of diagnostics.
313  i_ld2_2d, dz_int_2d
314  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: & ! 3-D versions of diagnostics.
315  i_ld2_3d, dz_int_3d
316 #endif
317 #ifdef DEBUG
318  integer :: max_debug_itt ; parameter(max_debug_itt=20)
319  real :: wt(szk_(g)+1), wt_tot, i_wt_tot, wt_itt
320  real, dimension(SZK_(G)+1) :: &
321  ri_k, tke_prev, dtke, dkappa, dtke_norm, &
322  ksrc_av ! The average through the iterations of k_src, in s-1.
323  real, dimension(SZK_(G)+1,0:max_debug_itt) :: &
324  tke_it1, n2_it1, sh2_it1, ksrc_it1, kappa_it1, kprev_it1
325  real, dimension(SZK_(G)+1,1:max_debug_itt) :: &
326  dkappa_it1, wt_it1, k_q_it1, d_dkappa_it1, dkappa_norm
327  real, dimension(SZK_(G),0:max_debug_itt) :: &
328  u_it1, v_it1, rho_it1, t_it1, s_it1
329  real, dimension(0:max_debug_itt) :: &
330  dk_wt_it1, dkpos_wt_it1, dkneg_wt_it1, k_mag
331  real, dimension(max_debug_itt) :: dt_it1
332 #endif
333  is = g%isc ; ie = g%iec; js = g%jsc ; je = g%jec ; nz = g%ke
334 
335  ! These are hard-coded for now. Perhaps these could be made dynamic later?
336  ! tol_dksrc = 0.5*tol_ksrc_chg ; tol_dksrc_low = 1.0 - 1.0/tol_ksrc_chg ?
337  tol_dksrc = 10.0 ; tol_dksrc_low = 0.95 ; tol2 = 2.0*cs%kappa_tol_err
338  dt_refinements = 5 ! Selected so that 1/2^dt_refinements < 1-tol_dksrc_low
339 
340  use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
341  new_kappa = .true. ; if (present(initialize_all)) new_kappa = initialize_all
342 
343  ri_crit = cs%Rino_crit
344  gr0 = gv%Rho0*gv%g_Earth ; g_r0 = gv%g_Earth/gv%Rho0
345 
346  k0dt = dt*cs%kappa_0
347  dz_massless = 0.1*sqrt(k0dt)
348 
349 !$OMP parallel do default(none) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,new_kappa, &
350 !$OMP tv,G,GV,CS,kappa_io,dz_massless,k0dt,p_surf,gR0,g_R0, &
351 !$OMP dt,tol_dksrc,tol_dksrc_low,tol2,dt_refinements, &
352 #ifdef ADD_DIAGNOSTICS
353 !$OMP I_Ld2_3d,dz_Int_3d, &
354 #endif
355 !$OMP Ri_crit,tke_io,kv_io) &
356 !$OMP private(h_2d,u_2d,v_2d,T_2d,S_2d,rho_2d,kappa_2d,nzc,dz, &
357 !$OMP u0xdz,v0xdz,T0xdz,S0xdz,kc,Idz,kf,I_dz_int,dz_in_lay, &
358 !$OMP dist_from_top,a1,b1,u,v,T,Sal,c1,d1,bd1,dz_Int,Norm, &
359 !$OMP dist_from_bot,I_L2_bdry,f2,pressure,T_int,Sal_int, &
360 !$OMP dbuoy_dT,dbuoy_dS,kappa,K_Q,N2,S2,dt_rem,kappa_avg, &
361 !$OMP tke_avg,local_src_avg,tke,kappa_out,kappa_src, &
362 !$OMP local_src,ks_kappa,ke_kappa,dt_now,dt_test,tol_max, &
363 !$OMP tol_min,tol_chg,u_test, v_test, T_test, S_test, &
364 !$OMP valid_dt,Idtt,k_src,dt_inc,dt_wt,kappa_mid,K_Q_tmp, &
365 #ifdef ADD_DIAGNOSTICS
366 !$OMP I_Ld2_1d,I_Ld2_2d, dz_Int_2d, &
367 #endif
368 !$OMP tke_pred,kappa_pred,tke_2d)
369 
370  do j=js,je
371  do k=1,nz ; do i=is,ie
372  h_2d(i,k) = h(i,j,k)*gv%H_to_m
373  u_2d(i,k) = u_in(i,j,k) ; v_2d(i,k) = v_in(i,j,k)
374  enddo ; enddo
375  if (use_temperature) then ; do k=1,nz ; do i=is,ie
376  t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
377  enddo ; enddo ; else ; do k=1,nz ; do i=is,ie
378  rho_2d(i,k) = gv%Rlay(k) ! Could be tv%Rho(i,j,k) ?
379  enddo ; enddo ; endif
380  if (.not.new_kappa) then ; do k=1,nz+1 ; do i=is,ie
381  kappa_2d(i,k) = kappa_io(i,j,k)
382  enddo ; enddo ; endif
383 
384 !---------------------------------------
385 ! Work on each column.
386 !---------------------------------------
387  do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
388  ! call cpu_clock_begin(id_clock_setup)
389  ! Store a transposed version of the initial arrays.
390  ! Any elimination of massless layers would occur here.
391  if (cs%eliminate_massless) then
392  nzc = 1
393  do k=1,nz
394  ! Zero out the thicknesses of all layers, even if they are unused.
395  dz(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
396  t0xdz(k) = 0.0 ; s0xdz(k) = 0.0
397 
398  ! Add a new layer if this one has mass.
399 ! if ((dz(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless)) nzc = nzc+1
400  if ((k>cs%nkml) .and. (dz(nzc) > 0.0) .and. &
401  (h_2d(i,k) > dz_massless)) nzc = nzc+1
402 
403  ! Only merge clusters of massless layers.
404 ! if ((dz(nzc) > dz_massless) .or. &
405 ! ((dz(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless))) nzc = nzc+1
406 
407  kc(k) = nzc
408  dz(nzc) = dz(nzc) + h_2d(i,k)
409  u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
410  v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
411  if (use_temperature) then
412  t0xdz(nzc) = t0xdz(nzc) + t_2d(i,k)*h_2d(i,k)
413  s0xdz(nzc) = s0xdz(nzc) + s_2d(i,k)*h_2d(i,k)
414  else
415  t0xdz(nzc) = t0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
416  s0xdz(nzc) = s0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
417  endif
418  enddo
419  kc(nz+1) = nzc+1
420  ! Set up Idz as the inverse of layer thicknesses.
421  do k=1,nzc ; idz(k) = 1.0 / dz(k) ; enddo
422 
423  ! Now determine kf, the fractional weight of interface kc when
424  ! interpolating between interfaces kc and kc+1.
425  kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
426  do k=2,nz
427  if (kc(k) > kc(k-1)) then
428  kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
429  else
430  kf(k) = dz_in_lay*idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
431  endif
432  enddo
433  kf(nz+1) = 0.0
434  else
435  do k=1,nz
436  dz(k) = h_2d(i,k)
437  u0xdz(k) = u_2d(i,k)*dz(k) ; v0xdz(k) = v_2d(i,k)*dz(k)
438  enddo
439  if (use_temperature) then
440  do k=1,nz
441  t0xdz(k) = t_2d(i,k)*dz(k) ; s0xdz(k) = s_2d(i,k)*dz(k)
442  enddo
443  else
444  do k=1,nz
445  t0xdz(k) = rho_2d(i,k)*dz(k) ; s0xdz(k) = rho_2d(i,k)*dz(k)
446  enddo
447  endif
448  nzc = nz
449  do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
450  ! Set up Idz as the inverse of layer thicknesses.
451  do k=1,nzc ; idz(k) = 1.0 / dz(k) ; enddo
452  endif
453 
454  ! Set up I_dz_int as the inverse of the distance between
455  ! adjacent layer centers.
456  i_dz_int(1) = 2.0 / dz(1)
457  dist_from_top(1) = 0.0
458  do k=2,nzc
459  i_dz_int(k) = 2.0 / (dz(k-1) + dz(k))
460  dist_from_top(k) = dist_from_top(k-1) + dz(k-1)
461  enddo
462  i_dz_int(nzc+1) = 2.0 / dz(nzc)
463 
464  ! Determine the velocities and thicknesses after eliminating massless
465  ! layers and applying a time-step of background diffusion.
466  if (nzc > 1) then
467  a1(2) = k0dt*i_dz_int(2)
468  b1 = 1.0 / (dz(1)+a1(2))
469  u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
470  t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
471  c1(2) = a1(2) * b1 ; d1 = dz(1) * b1 ! = 1 - c1
472  do k=2,nzc-1
473  bd1 = dz(k) + d1*a1(k)
474  a1(k+1) = k0dt*i_dz_int(k+1)
475  b1 = 1.0 / (bd1 + a1(k+1))
476  u(k) = b1 * (u0xdz(k) + a1(k)*u(k-1))
477  v(k) = b1 * (v0xdz(k) + a1(k)*v(k-1))
478  t(k) = b1 * (t0xdz(k) + a1(k)*t(k-1))
479  sal(k) = b1 * (s0xdz(k) + a1(k)*sal(k-1))
480  c1(k+1) = a1(k+1) * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
481  enddo
482  ! rho or T and S have insulating boundary conditions, u & v use no-slip
483  ! bottom boundary conditions (if kappa0 > 0).
484  ! For no-slip bottom boundary conditions
485  b1 = 1.0 / ((dz(nzc) + d1*a1(nzc)) + k0dt*i_dz_int(nzc+1))
486  u(nzc) = b1 * (u0xdz(nzc) + a1(nzc)*u(nzc-1))
487  v(nzc) = b1 * (v0xdz(nzc) + a1(nzc)*v(nzc-1))
488  ! For insulating boundary conditions
489  b1 = 1.0 / (dz(nzc) + d1*a1(nzc))
490  t(nzc) = b1 * (t0xdz(nzc) + a1(nzc)*t(nzc-1))
491  sal(nzc) = b1 * (s0xdz(nzc) + a1(nzc)*sal(nzc-1))
492  do k=nzc-1,1,-1
493  u(k) = u(k) + c1(k+1)*u(k+1) ; v(k) = v(k) + c1(k+1)*v(k+1)
494  t(k) = t(k) + c1(k+1)*t(k+1) ; sal(k) = sal(k) + c1(k+1)*sal(k+1)
495  enddo
496  else
497  ! This is correct, but probably unnecessary.
498  b1 = 1.0 / (dz(1) + k0dt*i_dz_int(2))
499  u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
500  b1 = 1.0 / dz(1)
501  t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
502  endif
503 
504  ! This uses half the harmonic mean of thicknesses to provide two estimates
505  ! of the boundary between cells, and the inverse of the harmonic mean to
506  ! weight the two estimates. The net effect is that interfaces around thin
507  ! layers have thin cells, and the total thickness adds up properly.
508  ! The top- and bottom- interfaces have zero thickness, consistent with
509  ! adding additional zero thickness layers.
510  dz_int(1) = 0.0 ; dz_int(2) = dz(1)
511  do k=2,nzc-1
512  norm = 1.0 / (dz(k)*(dz(k-1)+dz(k+1)) + 2.0*dz(k-1)*dz(k+1))
513  dz_int(k) = dz_int(k) + dz(k) * ( ((dz(k)+dz(k+1)) * dz(k-1)) * norm)
514  dz_int(k+1) = dz(k) * ( ((dz(k-1)+dz(k)) * dz(k+1)) * norm)
515  enddo
516  dz_int(nzc) = dz_int(nzc) + dz(nzc) ; dz_int(nzc+1) = 0.0
517 
518 #ifdef ADD_DIAGNOSTICS
519  do k=1,nzc+1 ; i_ld2_1d(k) = 0.0 ; enddo
520 #endif
521 
522  dist_from_bot = 0.0
523  do k=nzc,2,-1
524  dist_from_bot = dist_from_bot + dz(k)
525  i_l2_bdry(k) = (dist_from_top(k) + dist_from_bot)**2 / &
526  (dist_from_top(k) * dist_from_bot)**2
527  enddo
528  f2 = 0.25*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
529  (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
530 
531  ! Calculate thermodynamic coefficients and an initial estimate of N2.
532  if (use_temperature) then
533  pressure(1) = 0.0
534  if (associated(p_surf)) pressure(1) = p_surf(i,j)
535  do k=2,nzc
536  pressure(k) = pressure(k-1) + gr0*dz(k-1)
537  t_int(k) = 0.5*(t(k-1) + t(k))
538  sal_int(k) = 0.5*(sal(k-1) + sal(k))
539  enddo
540  call calculate_density_derivs(t_int, sal_int, pressure, dbuoy_dt, &
541  dbuoy_ds, 2, nzc-1, tv%eqn_of_state)
542  do k=2,nzc
543  dbuoy_dt(k) = -g_r0*dbuoy_dt(k)
544  dbuoy_ds(k) = -g_r0*dbuoy_ds(k)
545  enddo
546  else
547  do k=1,nzc+1 ; dbuoy_dt(k) = -g_r0 ; dbuoy_ds(k) = 0.0 ; enddo
548  endif
549 
550  ! ----------------------------------------------------
551  ! Calculate kappa, here defined at interfaces.
552  ! ----------------------------------------------------
553  if (new_kappa) then
554  do k=1,nzc+1 ; kappa(k) = 1.0 ; k_q(k) = 0.0 ; enddo
555  else
556  do k=1,nzc+1 ; kappa(k) = kappa_2d(i,k) ; k_q(k) = 0.0 ; enddo
557  endif
558 
559 #ifdef DEBUG
560  n2(1) = 0.0 ; n2(nzc+1) = 0.0
561  do k=2,nzc
562  n2(k) = max((dbuoy_dt(k) * (t0xdz(k-1)*idz(k-1) - t0xdz(k)*idz(k)) + &
563  dbuoy_ds(k) * (s0xdz(k-1)*idz(k-1) - s0xdz(k)*idz(k))) * &
564  i_dz_int(k), 0.0)
565  enddo
566  do k=1,nzc
567  u_it1(k,0) = u0xdz(k)*idz(k) ; v_it1(k,0) = v0xdz(k)*idz(k)
568  t_it1(k,0) = t0xdz(k)*idz(k) ; s_it1(k,0) = s0xdz(k)*idz(k)
569  enddo
570  do k=1,nzc+1
571  kprev_it1(k,0) = kappa(k) ; kappa_it1(k,0) = kappa(k)
572  tke_it1(k,0) = tke(k)
573  n2_it1(k,0) = n2(k) ; sh2_it1(k,0) = s2(k) ; ksrc_it1(k,0) = k_src(k)
574  enddo
575  do k=nzc+1,nz
576  u_it1(k,0) = 0.0 ; v_it1(k,0) = 0.0
577  t_it1(k,0) = 0.0 ; s_it1(k,0) = 0.0
578  kprev_it1(k+1,0) = 0.0 ; kappa_it1(k+1,0) = 0.0 ; tke_it1(k+1,0) = 0.0
579  n2_it1(k+1,0) = 0.0 ; sh2_it1(k+1,0) = 0.0 ; ksrc_it1(k+1,0) = 0.0
580  enddo
581  do itt=1,max_debug_itt
582  dt_it1(itt) = 0.0
583  do k=1,nz
584  u_it1(k,itt) = 0.0 ; v_it1(k,itt) = 0.0
585  t_it1(k,itt) = 0.0 ; s_it1(k,itt) = 0.0
586  rho_it1(k,itt) = 0.0
587  enddo
588  do k=1,nz+1
589  kprev_it1(k,itt) = 0.0 ; kappa_it1(k,itt) = 0.0 ; tke_it1(k,itt) = 0.0
590  n2_it1(k,itt) = 0.0 ; sh2_it1(k,itt) = 0.0
591  ksrc_it1(k,itt) = 0.0
592  dkappa_it1(k,itt) = 0.0 ; wt_it1(k,itt) = 0.0
593  k_q_it1(k,itt) = 0.0 ; d_dkappa_it1(k,itt) = 0.0
594  enddo
595  enddo
596  do k=1,nz+1 ; ksrc_av(k) = 0.0 ; enddo
597 #endif
598 
599  ! This call just calculates N2 and S2.
600  call calculate_projected_state(kappa, u, v, t, sal, 0.0, nzc, &
601  dz, i_dz_int, dbuoy_dt, dbuoy_ds, &
602  u, v, t, sal, n2=n2, s2=s2)
603  ! ----------------------------------------------------
604  ! Iterate
605  ! ----------------------------------------------------
606  dt_rem = dt
607  do k=1,nzc+1
608  kappa_avg(k) = 0.0 ; tke_avg(k) = 0.0
609  local_src_avg(k) = 0.0
610  ! Use the grid spacings to scale errors in the source.
611  if ( dz_int(k) > 0.0 ) &
612  local_src_avg(k) = 0.1 * k0dt * i_dz_int(k) / dz_int(k)
613  enddo
614 
615  ! call cpu_clock_end(id_clock_setup)
616 
617 ! do itt=1,CS%max_RiNo_it
618  do itt=1,cs%max_KS_it
619 
620  ! ----------------------------------------------------
621  ! Calculate new values of u, v, rho, N^2 and S.
622  ! ----------------------------------------------------
623 #ifdef DEBUG
624  do k=1,nzc+1
625  ri_k(k) = 1e3 ; if (s2(k) > 1e-3*n2(k)) ri_k(k) = n2(k) / s2(k)
626  tke_prev(k) = tke(k)
627  enddo
628 #endif
629 
630  ! call cpu_clock_begin(id_clock_KQ)
631  call find_kappa_tke(n2, s2, kappa, idz, dz_int, i_l2_bdry, f2, &
632  nzc, cs, k_q, tke, kappa_out, kappa_src, local_src)
633  ! call cpu_clock_end(id_clock_KQ)
634 
635  ! call cpu_clock_begin(id_clock_avg)
636  ! Determine the range of non-zero values of kappa_out.
637  ks_kappa = nz+1 ; ke_kappa = 0
638  do k=2,nzc ; if (kappa_out(k) > 0.0) then
639  ks_kappa = k ; exit
640  endif ; enddo
641  do k=nzc,ks_kappa,-1 ; if (kappa_out(k) > 0.0) then
642  ke_kappa = k ; exit
643  endif ; enddo
644  if (ke_kappa == nzc) kappa_out(nzc+1) = 0.0
645  ! call cpu_clock_end(id_clock_avg)
646 
647  ! Determine how long to use this value of kappa (dt_now).
648 
649  ! call cpu_clock_begin(id_clock_project)
650  if ((ke_kappa < ks_kappa) .or. (itt==cs%max_RiNo_it)) then
651  dt_now = dt_rem
652  else
653  ! Limit dt_now so that |k_src(k)-kappa_src(k)| < tol * local_src(k)
654  dt_test = dt_rem
655  do k=2,nzc
656  tol_max(k) = kappa_src(k) + tol_dksrc * local_src(k)
657  tol_min(k) = kappa_src(k) - tol_dksrc_low * local_src(k)
658  tol_chg(k) = tol2 * local_src_avg(k)
659  enddo
660 
661  do itt_dt=1,(cs%max_KS_it+1-itt)/2
662  ! The maximum number of times that the time-step is halved in
663  ! seeking an acceptable timestep is reduced with each iteration,
664  ! so that as the maximum number of iterations is approached, the
665  ! whole remaining timestep is used. Typically, an acceptable
666  ! timestep is found long before the minimum is reached, so the
667  ! value of max_KS_it may be unimportant, especially if it is large
668  ! enough.
669  call calculate_projected_state(kappa_out, u, v, t, sal, 0.5*dt_test, nzc, &
670  dz, i_dz_int, dbuoy_dt, dbuoy_ds, &
671  u_test, v_test, t_test, s_test, n2, s2, &
672  ks_int = ks_kappa, ke_int = ke_kappa)
673  valid_dt = .true.
674  idtt = 1.0 / dt_test
675  do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
676  if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
677  k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
678  ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
679  if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
680  (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))) then
681  valid_dt = .false. ; exit
682  endif
683  else
684  if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))) then
685  valid_dt = .false. ; k_src(k) = 0.0 ; exit
686  endif
687  endif
688  enddo
689 
690  if (valid_dt) exit
691  dt_test = 0.5*dt_test
692  enddo
693  if ((dt_test < dt_rem) .and. valid_dt) then
694  dt_inc = 0.5*dt_test
695  do itt_dt=1,dt_refinements
696  call calculate_projected_state(kappa_out, u, v, t, sal, &
697  0.5*(dt_test+dt_inc), nzc, dz, i_dz_int, dbuoy_dt, &
698  dbuoy_ds, u_test, v_test, t_test, s_test, n2, s2, &
699  ks_int = ks_kappa, ke_int = ke_kappa)
700  valid_dt = .true.
701  idtt = 1.0 / (dt_test+dt_inc)
702  do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
703  if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
704  k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
705  ((ri_crit*s2(k) - n2(k)) / &
706  (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
707  if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
708  (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))) then
709  valid_dt = .false. ; exit
710  endif
711  else
712  if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))) then
713  valid_dt = .false. ; k_src(k) = 0.0 ; exit
714  endif
715  endif
716  enddo
717 
718  if (valid_dt) dt_test = dt_test + dt_inc
719  dt_inc = 0.5*dt_inc
720  enddo
721  else
722  dt_inc = 0.0
723  endif
724 
725  dt_now = min(dt_test*(1.0+cs%kappa_tol_err)+dt_inc,dt_rem)
726  do k=2,nzc
727  local_src_avg(k) = local_src_avg(k) + dt_now * local_src(k)
728  enddo
729  endif ! Are all the values of kappa_out 0?
730  ! call cpu_clock_end(id_clock_project)
731 
732  ! The state has already been projected forward. Now find new values of kappa.
733 
734  if (ke_kappa < ks_kappa) then
735  ! There is no mixing now, and will not be again.
736  ! call cpu_clock_begin(id_clock_avg)
737  dt_wt = dt_rem / dt ; dt_rem = 0.0
738  do k=1,nzc+1
739  kappa_mid(k) = 0.0
740  ! This would be here but does nothing.
741  ! kappa_avg(K) = kappa_avg(K) + kappa_mid(K)*dt_wt
742  tke_avg(k) = tke_avg(k) + dt_wt*tke(k)
743 #ifdef DEBUG
744  tke_pred(k) = tke(k) ; kappa_pred(k) = 0.0 ; kappa(k) = 0.0
745 #endif
746  enddo
747  ! call cpu_clock_end(id_clock_avg)
748  else
749  ! call cpu_clock_begin(id_clock_project)
750  call calculate_projected_state(kappa_out, u, v, t, sal, dt_now, nzc, &
751  dz, i_dz_int, dbuoy_dt, dbuoy_ds, &
752  u_test, v_test, t_test, s_test, n2=n2, s2=s2, &
753  ks_int = ks_kappa, ke_int = ke_kappa)
754  ! call cpu_clock_end(id_clock_project)
755 
756  ! call cpu_clock_begin(id_clock_KQ)
757  do k=1,nzc+1 ; k_q_tmp(k) = k_q(k) ; enddo
758  call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
759  nzc, cs, k_q_tmp, tke_pred, kappa_pred)
760  ! call cpu_clock_end(id_clock_KQ)
761 
762  ks_kappa = nz+1 ; ke_kappa = 0
763  do k=1,nzc+1
764  kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
765  if ((kappa_mid(k) > 0.0) .and. (k<ks_kappa)) ks_kappa = k
766  if (kappa_mid(k) > 0.0) ke_kappa = k
767  enddo
768 
769  ! call cpu_clock_begin(id_clock_project)
770  call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, &
771  dz, i_dz_int, dbuoy_dt, dbuoy_ds, &
772  u_test, v_test, t_test, s_test, n2=n2, s2=s2, &
773  ks_int = ks_kappa, ke_int = ke_kappa)
774  ! call cpu_clock_end(id_clock_project)
775 
776  ! call cpu_clock_begin(id_clock_KQ)
777  call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
778  nzc, cs, k_q, tke_pred, kappa_pred)
779  ! call cpu_clock_end(id_clock_KQ)
780 
781  ! call cpu_clock_begin(id_clock_avg)
782  dt_wt = dt_now / dt ; dt_rem = dt_rem - dt_now
783  do k=1,nzc+1
784  kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
785  kappa_avg(k) = kappa_avg(k) + kappa_mid(k)*dt_wt
786  tke_avg(k) = tke_avg(k) + dt_wt*0.5*(tke_pred(k) + tke(k))
787  kappa(k) = kappa_pred(k) ! First guess for the next iteration.
788  enddo
789  ! call cpu_clock_end(id_clock_avg)
790  endif
791 
792  if (dt_rem > 0.0) then
793  ! Update the values of u, v, T, Sal, N2, and S2 for the next iteration.
794  ! call cpu_clock_begin(id_clock_project)
795  call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, &
796  dz, i_dz_int, dbuoy_dt, dbuoy_ds, &
797  u, v, t, sal, n2, s2)
798  ! call cpu_clock_end(id_clock_project)
799  endif
800 
801 #ifdef DEBUG
802  if (itt <= max_debug_itt) then
803  dt_it1(itt) = dt_now
804  dk_wt_it1(itt) = 0.0 ; dkpos_wt_it1(itt) = 0.0 ; dkneg_wt_it1(itt) = 0.0
805  k_mag(itt) = 0.0
806  wt_itt = 1.0/real(itt) ; wt_tot = 0.0
807  do k=1,nzc+1
808  ksrc_av(k) = (1.0-wt_itt)*ksrc_av(k) + wt_itt*k_src(k)
809  wt_tot = wt_tot + dz_int(k) * ksrc_av(k)
810  enddo
811  ! Use the 1/0=0 convention.
812  i_wt_tot = 0.0 ; if (wt_tot > 0.0) i_wt_tot = 1.0/wt_tot
813 
814  do k=1,nzc+1
815  wt(k) = (dz_int(k)*ksrc_av(k)) * i_wt_tot
816  k_mag(itt) = k_mag(itt) + wt(k)*kappa_mid(k)
817  dkappa_it1(k,itt) = kappa_pred(k) - kappa_out(k)
818  dk_wt_it1(itt) = dk_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
819  if (dk > 0.0) then
820  dkpos_wt_it1(itt) = dkpos_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
821  else
822  dkneg_wt_it1(itt) = dkneg_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
823  endif
824  wt_it1(k,itt) = wt(k)
825  enddo
826  endif
827  do k=1,nzc+1
828  ri_k(k) = 1e3 ; if (n2(k) < 1e3 * s2(k)) ri_k(k) = n2(k) / s2(k)
829  dtke(k) = tke_pred(k) - tke(k)
830  dtke_norm(k) = dtke(k) / (0.5*(tke(k) + tke_pred(k)))
831  dkappa(k) = kappa_pred(k) - kappa_out(k)
832  enddo
833  if (itt <= max_debug_itt) then
834  do k=1,nzc
835  u_it1(k,itt) = u(k) ; v_it1(k,itt) = v(k)
836  t_it1(k,itt) = t(k) ; s_it1(k,itt) = sal(k)
837  enddo
838  do k=1,nzc+1
839  kprev_it1(k,itt)=kappa_out(k)
840  kappa_it1(k,itt)=kappa_mid(k) ; tke_it1(k,itt) = 0.5*(tke(k)+tke_pred(k))
841  n2_it1(k,itt)=n2(k) ; sh2_it1(k,itt)=s2(k)
842  ksrc_it1(k,itt) = kappa_src(k)
843  k_q_it1(k,itt) = kappa_out(k) / (tke(k))
844  if (itt > 1) then
845  if (abs(dkappa_it1(k,itt-1)) > 1e-20) &
846  d_dkappa_it1(k,itt) = dkappa_it1(k,itt) / dkappa_it1(k,itt-1)
847  endif
848  dkappa_norm(k,itt) = dkappa(k) / max(0.5*(kappa_pred(k) + kappa_out(k)), 1e-100)
849  enddo
850  endif
851 #endif
852 
853  if (dt_rem <= 0.0) exit
854 
855  enddo ! end itt loop
856 
857  ! call cpu_clock_begin(id_clock_setup)
858  ! Extrapolate from the vertically reduced grid back to the original layers.
859  if (nz == nzc) then
860  do k=1,nz+1
861  kappa_2d(i,k) = kappa_avg(k)
862  tke_2d(i,k) = tke(k)
863  enddo
864  else
865  do k=1,nz+1
866  if (kf(k) == 0.0) then
867  kappa_2d(i,k) = kappa_avg(kc(k))
868  tke_2d(i,k) = tke_avg(kc(k))
869  else
870  kappa_2d(i,k) = (1.0-kf(k)) * kappa_avg(kc(k)) + &
871  kf(k) * kappa_avg(kc(k)+1)
872  tke_2d(i,k) = (1.0-kf(k)) * tke_avg(kc(k)) + &
873  kf(k) * tke_avg(kc(k)+1)
874  endif
875  enddo
876  endif
877 #ifdef ADD_DIAGNOSTICS
878  i_ld2_2d(i,1) = 0.0 ; dz_int_2d(i,1) = dz_int(1)
879  do k=2,nzc
880  i_ld2_2d(i,k) = (n2(k) / cs%lambda**2 + f2) / &
881  max(tke(k),1e-30) + i_l2_bdry(k)
882  dz_int_2d(i,k) = dz_int(k)
883  enddo
884  i_ld2_2d(i,nzc+1) = 0.0 ; dz_int_2d(i,nzc+1) = dz_int(nzc+1)
885  do k=nzc+2,nz+1
886  i_ld2_2d(i,k) = 0.0 ; dz_int_2d(i,k) = 0.0
887  enddo
888 #endif
889  ! call cpu_clock_end(id_clock_setup)
890  else ! Land points, still inside the i-loop.
891  do k=1,nz+1
892  kappa_2d(i,k) = 0.0 ; tke_2d(i,k) = 0.0
893 #ifdef ADD_DIAGNOSTICS
894  i_ld2_2d(i,k) = 0.0
895  dz_int_2d(i,k) = dz_int(k)
896 #endif
897  enddo
898  endif ; enddo ! i-loop
899 
900  do k=1,nz+1 ; do i=is,ie
901  kappa_io(i,j,k) = g%mask2dT(i,j) * kappa_2d(i,k)
902  tke_io(i,j,k) = g%mask2dT(i,j) * tke_2d(i,k)
903  kv_io(i,j,k) = ( g%mask2dT(i,j) * kappa_2d(i,k) ) * cs%Prandtl_turb
904 #ifdef ADD_DIAGNOSTICS
905  i_ld2_3d(i,j,k) = i_ld2_2d(i,k)
906  dz_int_3d(i,j,k) = dz_int_2d(i,k)
907 #endif
908  enddo ; enddo
909 
910  enddo ! end of j-loop
911 
912  if (cs%debug) then
913  call hchksum(kappa_io,"kappa",g%HI)
914  call hchksum(tke_io,"tke",g%HI)
915  endif
916 
917  if (cs%id_Kd_shear > 0) call post_data(cs%id_Kd_shear, kappa_io, cs%diag)
918  if (cs%id_TKE > 0) call post_data(cs%id_TKE, tke_io, cs%diag)
919 #ifdef ADD_DIAGNOSTICS
920  if (cs%id_ILd2 > 0) call post_data(cs%id_ILd2, i_ld2_3d, cs%diag)
921  if (cs%id_dz_Int > 0) call post_data(cs%id_dz_Int, dz_int_3d, cs%diag)
922 #endif
923 
924  return
925 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ calculate_projected_state()

subroutine mom_kappa_shear::calculate_projected_state ( real, dimension(nz+1), intent(in)  kappa,
real, dimension(nz), intent(in)  u0,
real, dimension(nz), intent(in)  v0,
real, dimension(nz), intent(in)  T0,
real, dimension(nz), intent(in)  S0,
real, intent(in)  dt,
integer, intent(in)  nz,
real, dimension(nz), intent(in)  dz,
real, dimension(nz+1), intent(in)  I_dz_int,
real, dimension(nz+1), intent(in)  dbuoy_dT,
real, dimension(nz+1), intent(in)  dbuoy_dS,
real, dimension(nz), intent(inout)  u,
real, dimension(nz), intent(inout)  v,
real, dimension(nz), intent(inout)  T,
real, dimension(nz), intent(inout)  Sal,
real, dimension(nz+1), intent(inout), optional  N2,
real, dimension(nz+1), intent(inout), optional  S2,
integer, intent(in), optional  ks_int,
integer, intent(in), optional  ke_int 
)
private
Parameters
[in]nzThe number of layers (after eliminating massless layers?).
[in]kappaThe diapycnal diffusivity at interfaces, in m2 s-1.
[in]u0The initial zonal velocity, in m s-1.
[in]v0The initial meridional velocity, in m s-1.
[in]t0The initial temperature, in C.
[in]s0The initial salinity, in PSU.
[in]dzThe grid spacing of layers, in m.
[in]i_dz_intThe inverse of the layer's thicknesses, in m-1.
[in]dbuoy_dtThe partial derivative of buoyancy with temperature, in m s-2 C-1.
[in]dbuoy_dsThe partial derivative of buoyancy with salinity, in m s-2 PSU-1.
[in]dtThe time step in s.
[in,out]uThe zonal velocity after dt, in m s-1.
[in,out]vThe meridional velocity after dt, in m s-1.
[in,out]tThe temperature after dt, in C.
[in,out]salThe salinity after dt, in PSU.
[in,out]n2The buoyancy frequency squared at interfaces,
[in,out]s2The squared shear at interfaces, in s-2.
[in]ks_intThe topmost k-index with a non-zero diffusivity.
[in]ke_intThe bottommost k-index with a non-zero diffusivity.

Definition at line 931 of file MOM_kappa_shear.F90.

Referenced by calculate_kappa_shear().

931 !< This subroutine calculates the velocities, temperature and salinity that
932 !! the water column will have after mixing for dt with diffusivities kappa. It
933 !! may also calculate the projected buoyancy frequency and shear.
934  integer, intent(in) :: nz !< The number of layers (after eliminating massless
935  !! layers?).
936  real, dimension(nz+1), intent(in) :: kappa !< The diapycnal diffusivity at interfaces,
937  !! in m2 s-1.
938  real, dimension(nz), intent(in) :: u0 !< The initial zonal velocity, in m s-1.
939  real, dimension(nz), intent(in) :: v0 !< The initial meridional velocity, in m s-1.
940  real, dimension(nz), intent(in) :: t0 !< The initial temperature, in C.
941  real, dimension(nz), intent(in) :: s0 !< The initial salinity, in PSU.
942  real, dimension(nz), intent(in) :: dz !< The grid spacing of layers, in m.
943  real, dimension(nz+1), intent(in) :: i_dz_int !< The inverse of the layer's thicknesses,
944  !! in m-1.
945  real, dimension(nz+1), intent(in) :: dbuoy_dt !< The partial derivative of buoyancy with
946  !! temperature, in m s-2 C-1.
947  real, dimension(nz+1), intent(in) :: dbuoy_ds !< The partial derivative of buoyancy with
948  !! salinity, in m s-2 PSU-1.
949  real, intent(in) :: dt !< The time step in s.
950  real, dimension(nz), intent(inout) :: u !< The zonal velocity after dt, in m s-1.
951  real, dimension(nz), intent(inout) :: v !< The meridional velocity after dt, in m s-1.
952  real, dimension(nz), intent(inout) :: t !< The temperature after dt, in C.
953  real, dimension(nz), intent(inout) :: sal !< The salinity after dt, in PSU.
954  real, dimension(nz+1), optional, &
955  intent(inout) :: n2 !< The buoyancy frequency squared at interfaces,
956  !! in s-2.
957  real, dimension(nz+1), optional, &
958  intent(inout) :: s2 !< The squared shear at interfaces, in s-2.
959  integer, optional, intent(in) :: ks_int !< The topmost k-index with a non-zero diffusivity.
960  integer, optional, intent(in) :: ke_int !< The bottommost k-index with a non-zero
961  !! diffusivity.
962 
963  ! Arguments: kappa - The diapycnal diffusivity at interfaces, in m2 s-1.
964  ! (in) Sh - The shear at interfaces, in s-1.
965  ! (in) u0 - The initial zonal velocity, in m s-1.
966  ! (in) v0 - The initial meridional velocity, in m s-1.
967  ! (in) T0 - The initial temperature, in C.
968  ! (in) S0 - The initial salinity, in PSU.
969  ! (in) nz - The number of layers (after eliminating massless layers?).
970  ! (in) dz - The grid spacing of layers, in m.
971  ! (in) I_dz_int - The inverse of the layer's thicknesses, in m-1.
972  ! (in) dbuoy_dT - The partial derivative of buoyancy with temperature,
973  ! in m s-2 C-1.
974  ! (in) dbuoy_dS - The partial derivative of buoyancy with salinity,
975  ! in m s-2 PSU-1.
976  ! (in) dt - The time step in s.
977  ! (in) nz - The number of layers to work on.
978  ! (out) u - The zonal velocity after dt, in m s-1.
979  ! (out) v - The meridional velocity after dt, in m s-1.
980  ! (in) T - The temperature after dt, in C.
981  ! (in) Sal - The salinity after dt, in PSU.
982  ! (out) N2 - The buoyancy frequency squared at interfaces, in s-2.
983  ! (out) S2 - The squared shear at interfaces, in s-2.
984  ! (in,opt) ks_int - The topmost k-index with a non-zero diffusivity.
985  ! (in,opt) ke_int - The bottommost k-index with a non-zero diffusivity.
986 
987  ! UNCOMMENT THE FOLLOWING IF NOT CONTAINED IN THE OUTER SUBROUTINE.
988  real, dimension(nz+1) :: c1
989  real :: a_a, a_b, b1, d1, bd1, b1nz_0
990  integer :: k, ks, ke
991 
992  ks = 1 ; ke = nz
993  if (present(ks_int)) ks = max(ks_int-1,1)
994  if (present(ke_int)) ke = min(ke_int,nz)
995 
996  if (ks > ke) return
997 
998  if (dt > 0.0) then
999  a_b = dt*(kappa(ks+1)*i_dz_int(ks+1))
1000  b1 = 1.0 / (dz(ks) + a_b)
1001  c1(ks+1) = a_b * b1 ; d1 = dz(ks) * b1 ! = 1 - c1
1002 
1003  u(ks) = (b1 * dz(ks))*u0(ks) ; v(ks) = (b1 * dz(ks))*v0(ks)
1004  t(ks) = (b1 * dz(ks))*t0(ks) ; sal(ks) = (b1 * dz(ks))*s0(ks)
1005  do k=ks+1,ke-1
1006  a_a = a_b
1007  a_b = dt*(kappa(k+1)*i_dz_int(k+1))
1008  bd1 = dz(k) + d1*a_a
1009  b1 = 1.0 / (bd1 + a_b)
1010  c1(k+1) = a_b * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
1011 
1012  u(k) = b1 * (dz(k)*u0(k) + a_a*u(k-1))
1013  v(k) = b1 * (dz(k)*v0(k) + a_a*v(k-1))
1014  t(k) = b1 * (dz(k)*t0(k) + a_a*t(k-1))
1015  sal(k) = b1 * (dz(k)*s0(k) + a_a*sal(k-1))
1016  enddo
1017  ! T and S have insulating boundary conditions, u & v use no-slip
1018  ! bottom boundary conditions at the solid bottom.
1019 
1020  ! For insulating boundary conditions or mixing simply stopping, use...
1021  a_a = a_b
1022  b1 = 1.0 / (dz(ke) + d1*a_a)
1023  t(ke) = b1 * (dz(ke)*t0(ke) + a_a*t(ke-1))
1024  sal(ke) = b1 * (dz(ke)*s0(ke) + a_a*sal(ke-1))
1025 
1026  ! There is no distinction between the effective boundary conditions for
1027  ! tracers and velocities if the mixing is separated from the bottom, but if
1028  ! the mixing goes all the way to the bottom, use no-slip BCs for velocities.
1029  if (ke == nz) then
1030  a_b = dt*(kappa(nz+1)*i_dz_int(nz+1))
1031  b1nz_0 = 1.0 / ((dz(nz) + d1*a_a) + a_b)
1032  else
1033  b1nz_0 = b1
1034  endif
1035  u(ke) = b1nz_0 * (dz(ke)*u0(ke) + a_a*u(ke-1))
1036  v(ke) = b1nz_0 * (dz(ke)*v0(ke) + a_a*v(ke-1))
1037 
1038  do k=ke-1,ks,-1
1039  u(k) = u(k) + c1(k+1)*u(k+1)
1040  v(k) = v(k) + c1(k+1)*v(k+1)
1041  t(k) = t(k) + c1(k+1)*t(k+1)
1042  sal(k) = sal(k) + c1(k+1)*sal(k+1)
1043  enddo
1044  else ! dt <= 0.0
1045  do k=1,nz
1046  u(k) = u0(k) ; v(k) = v0(k) ; t(k) = t0(k) ; sal(k) = s0(k)
1047  enddo
1048  endif
1049 
1050  if (present(s2)) then
1051  s2(1) = 0.0 ; s2(nz+1) = 0.0
1052  if (ks > 1) &
1053  s2(ks) = ((u(ks)-u0(ks-1))**2 + (v(ks)-v0(ks-1))**2) * i_dz_int(ks)**2
1054  do k=ks+1,ke
1055  s2(k) = ((u(k)-u(k-1))**2 + (v(k)-v(k-1))**2) * i_dz_int(k)**2
1056  enddo
1057  if (ke<nz) &
1058  s2(ke+1) = ((u0(ke+1)-u(ke))**2 + (v0(ke+1)-v(ke))**2) * i_dz_int(ke+1)**2
1059  endif
1060 
1061  if (present(n2)) then
1062  n2(1) = 0.0 ; n2(nz+1) = 0.0
1063  if (ks > 1) &
1064  n2(ks) = max(0.0, i_dz_int(ks) * &
1065  (dbuoy_dt(ks) * (t0(ks-1)-t(ks)) + dbuoy_ds(ks) * (s0(ks-1)-sal(ks))))
1066  do k=ks+1,ke
1067  n2(k) = max(0.0, i_dz_int(k) * &
1068  (dbuoy_dt(k) * (t(k-1)-t(k)) + dbuoy_ds(k) * (sal(k-1)-sal(k))))
1069  enddo
1070  if (ke<nz) &
1071  n2(ke+1) = max(0.0, i_dz_int(ke+1) * &
1072  (dbuoy_dt(ke+1) * (t(ke)-t0(ke+1)) + dbuoy_ds(ke+1) * (sal(ke)-s0(ke+1))))
1073  endif
1074 
Here is the caller graph for this function:

◆ find_kappa_tke()

subroutine mom_kappa_shear::find_kappa_tke ( real, dimension(nz+1), intent(in)  N2,
real, dimension(nz+1), intent(in)  S2,
real, dimension(nz+1), intent(in)  kappa_in,
real, dimension(nz), intent(in)  Idz,
real, dimension(nz+1), intent(in)  dz_Int,
real, dimension(nz+1), intent(in)  I_L2_bdry,
real, intent(in)  f2,
integer, intent(in)  nz,
type(kappa_shear_cs), pointer  CS,
real, dimension(nz+1), intent(inout)  K_Q,
real, dimension(nz+1), intent(out)  tke,
real, dimension(nz+1), intent(out)  kappa,
real, dimension(nz+1), intent(out), optional  kappa_src,
real, dimension(nz+1), intent(out), optional  local_src 
)
private

This subroutine calculates new, consistent estimates of TKE and kappa.

Parameters
[in]nzThe number of layers to work on.
[in]n2The buoyancy frequency squared at interfaces, in s-2.
[in]s2The squared shear at interfaces, in s-2.
[in]kappa_inThe initial guess at the diffusivity, in m2 s-1.
[in]dz_intThe thicknesses associated with interfaces, in m.
[in]i_l2_bdryThe inverse of the squared distance to boundaries, m2.
[in]idzThe inverse grid spacing of layers, in m-1.
[in]f2The squared Coriolis parameter, in s-2.
csA pointer to this module's control structure.
[in,out]k_qThe shear-driven diapycnal diffusivity divided by the turbulent kinetic energy per unit mass at interfaces, in s.
[out]tkeThe turbulent kinetic energy per unit mass at interfaces, in units of m2 s-2.
[out]kappaThe diapycnal diffusivity at interfaces, in m2 s-1.
[out]kappa_srcThe source term for kappa, in s-1.
[out]local_srcThe sum of all local sources for kappa,

Definition at line 1080 of file MOM_kappa_shear.F90.

Referenced by calculate_kappa_shear().

1080  integer, intent(in) :: nz !< The number of layers to work on.
1081  real, dimension(nz+1), intent(in) :: n2 !< The buoyancy frequency squared at interfaces,
1082  !! in s-2.
1083  real, dimension(nz+1), intent(in) :: s2 !< The squared shear at interfaces, in s-2.
1084  real, dimension(nz+1), intent(in) :: kappa_in !< The initial guess at the diffusivity,
1085  !! in m2 s-1.
1086  real, dimension(nz+1), intent(in) :: dz_int !< The thicknesses associated with interfaces,
1087  !! in m.
1088  real, dimension(nz+1), intent(in) :: i_l2_bdry !< The inverse of the squared distance to
1089  !! boundaries, m2.
1090  real, dimension(nz), intent(in) :: idz !< The inverse grid spacing of layers, in m-1.
1091  real, intent(in) :: f2 !< The squared Coriolis parameter, in s-2.
1092  type(kappa_shear_cs), pointer :: cs !< A pointer to this module's control structure.
1093  real, dimension(nz+1), intent(inout) :: k_q !< The shear-driven diapycnal diffusivity divided by
1094  !! the turbulent kinetic energy per unit mass at
1095  !! interfaces, in s.
1096  real, dimension(nz+1), intent(out) :: tke !< The turbulent kinetic energy per unit mass at
1097  !! interfaces, in units of m2 s-2.
1098  real, dimension(nz+1), intent(out) :: kappa !< The diapycnal diffusivity at interfaces,
1099  !! in m2 s-1.
1100  real, dimension(nz+1), optional, &
1101  intent(out) :: kappa_src !< The source term for kappa, in s-1.
1102  real, dimension(nz+1), optional, &
1103  intent(out) :: local_src !< The sum of all local sources for kappa,
1104  !! in s-1.
1105 ! This subroutine calculates new, consistent estimates of TKE and kappa.
1106 
1107 ! Arguments: N2 - The buoyancy frequency squared at interfaces, in s-2.
1108 ! (in) S2 - The squared shear at interfaces, in s-2.
1109 ! (in) kappa_in - The initial guess at the diffusivity, in m2 s-1.
1110 ! (in) Idz - The inverse grid spacing of layers, in m-1.
1111 ! (in) dz_Int - The thicknesses associated with interfaces, in m.
1112 ! (in) I_L2_bdry - The inverse of the squared distance to boundaries, m2.
1113 ! (in) f2 - The squared Coriolis parameter, in s-2.
1114 ! (in) nz - The number of layers to work on.
1115 ! (in) CS - A pointer to this module's control structure.
1116 ! (inout) K_Q - The shear-driven diapycnal diffusivity divided by the
1117 ! turbulent kinetic energy per unit mass at interfaces, in s.
1118 ! (out) tke - The turbulent kinetic energy per unit mass at interfaces,
1119 ! in units of m2 s-2.
1120 ! (out) kappa - The diapycnal diffusivity at interfaces, in m2 s-1.
1121 ! (out,opt) kappa_src - The source term for kappa, in s-1.
1122 ! (out,opt) local_src - The sum of all local sources for kappa, in s-1.
1123 
1124 ! UNCOMMENT THE FOLLOWING IF NOT CONTAINED IN Calculate_kappa_shear
1125  real, dimension(nz) :: &
1126  aq, & ! aQ is the coupling between adjacent interfaces in the TKE
1127  ! equations, in m s-1.
1128  dqdz ! Half the partial derivative of TKE with depth, m s-2.
1129  real, dimension(nz+1) :: &
1130  dk, & ! The change in kappa, in m2 s-1.
1131  dq, & ! The change in TKE, in m2 s-1.
1132  cq, ck, & ! cQ and cK are the upward influences in the tridiagonal and
1133  ! hexadiagonal solvers for the TKE and kappa equations, ND.
1134  i_ld2, & ! 1/Ld^2, where Ld is the effective decay length scale
1135  ! for kappa, in units of m-2.
1136  tke_decay, & ! The local TKE decay rate in s-1.
1137  k_src, & ! The source term in the kappa equation, in s-1.
1138  dqmdk, & ! With Newton's method the change in dQ(k-1) due to dK(k), s.
1139  dkdq, & ! With Newton's method the change in dK(k) due to dQ(k), s-1.
1140  e1 ! The fractional change in a layer TKE due to a change in the
1141  ! TKE of the layer above when all the kappas below are 0.
1142  ! e1 is nondimensional, and 0 < e1 < 1.
1143  real :: tke_src ! The net source of TKE due to mixing against the shear
1144  ! and stratification, in m2 s-3. (For convenience,
1145  ! a term involving the non-dissipation of q0 is also
1146  ! included here.)
1147  real :: bq, bk ! The inverse of the pivot in the tridiagonal equations.
1148  real :: bd1 ! A term in the denominator of bQ or bK.
1149  real :: cqcomp, ckcomp ! 1 - cQ or 1 - cK in the tridiagonal equations.
1150  real :: c_s2 ! The coefficient for the decay of TKE due to
1151  ! shear (i.e. proportional to |S|*tke), nondimensional.
1152  real :: c_n2 ! The coefficient for the decay of TKE due to
1153  ! stratification (i.e. proportional to N*tke), nondim.
1154  real :: ri_crit ! The critical shear Richardson number for shear-
1155  ! driven mixing. The theoretical value is 0.25.
1156  real :: q0 ! The background level of TKE, in m2 s-2.
1157  real :: ilambda2 ! 1.0 / CS%lambda**2.
1158  real :: tke_min ! The minimum value of shear-driven TKE that can be
1159  ! solved for, in m2 s-2.
1160  real :: kappa0 ! The background diapycnal diffusivity, in m2 s-1.
1161  real :: max_err ! The maximum value of norm_err in a column, nondim.
1162  real :: kappa_trunc ! Diffusivities smaller than this are rounded to 0, m2 s-1.
1163 
1164  real :: eden1, eden2, i_eden, ome ! Variables used in calculating e1.
1165  real :: diffusive_src ! The diffusive source in the kappa equation, in m s-1.
1166  real :: chg_by_k0 ! The value of k_src that leads to an increase of
1167  ! kappa_0 if only the diffusive term is a sink, in s-1.
1168 
1169  real :: kappa_mean ! A mean value of kappa, in m2 s-1.
1170  real :: newton_test ! The value of relative error that will cause the next
1171  ! iteration to use Newton's method.
1172  ! Temporary variables used in the Newton's method iterations.
1173  real :: decay_term, i_q, kap_src, v1, v2
1174 
1175  real :: tol_err ! The tolerance for max_err that determines when to
1176  ! stop iterating.
1177  real :: newton_err ! The tolerance for max_err that determines when to
1178  ! start using Newton's method. Empirically, an initial
1179  ! value of about 0.2 seems to be most efficient.
1180  real, parameter :: roundoff = 1.0e-16 ! A negligible fractional change in TKE.
1181  ! This could be larger but performance gains are small.
1182 
1183  logical :: tke_noflux_bottom_bc = .false. ! Specify the boundary conditions
1184  logical :: tke_noflux_top_bc = .false. ! that are applied to the TKE eqns.
1185  logical :: do_newton ! If .true., use Newton's method for the next iteration.
1186  logical :: abort_newton ! If .true., an Newton's method has encountered a 0
1187  ! pivot, and should not have been used.
1188  logical :: was_newton ! The value of do_Newton before checking convergence.
1189  logical :: within_tolerance ! If .true., all points are within tolerance to
1190  ! enable this subroutine to return.
1191  integer :: ks_src, ke_src ! The range indices that have nonzero k_src.
1192  integer :: ks_kappa, ke_kappa, ke_tke ! The ranges of k-indices that are or
1193  integer :: ks_kappa_prev, ke_kappa_prev ! were being worked on.
1194  integer :: itt, k, k2
1195 #ifdef DEBUG
1196  integer :: max_debug_itt ; parameter(max_debug_itt=20)
1197  real :: k_err_lin, q_err_lin
1198  real, dimension(nz+1) :: &
1199  kappa_prev, & ! The value of kappa at the start of the current iteration, in m2 s-1.
1200  tke_prev ! The value of TKE at the start of the current iteration, in m2 s-2.
1201  real, dimension(nz+1,1:max_debug_itt) :: &
1202  tke_it1, kappa_it1, kprev_it1, & ! Various values from each iteration.
1203  dkappa_it1, k_q_it1, d_dkappa_it1, dkappa_norm_it1
1204  real :: norm_err ! The absolute change in kappa between iterations,
1205  ! normalized by the value of kappa, nondim.
1206  real :: max_tke_err, min_tke_err, tke_err(nz) ! Various normalized TKE changes.
1207  integer :: it2
1208 #endif
1209 
1210  c_n2 = cs%C_N**2 ; c_s2 = cs%C_S**2
1211  q0 = cs%TKE_bg ; kappa0 = cs%kappa_0 ; tke_min = max(cs%TKE_bg,1.0e-20)
1212  ri_crit = cs%Rino_crit
1213  ilambda2 = 1.0 / cs%lambda**2
1214  kappa_trunc = 0.01*kappa0 ! ### CHANGE THIS HARD-WIRING LATER?
1215  do_newton = .false. ; abort_newton = .false.
1216  tol_err = cs%kappa_tol_err
1217  newton_err = 0.2 ! This initial value may be automatically reduced later.
1218 
1219  ks_kappa = 2 ; ke_kappa = nz ; ks_kappa_prev = 2 ; ke_kappa_prev = nz
1220 
1221  ke_src = 0 ; ks_src = nz+1
1222  do k=2,nz
1223  if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
1224 ! Ri = N2(K) / S2(K)
1225 ! k_src(K) = (2.0 * CS%Shearmix_rate * sqrt(S2(K))) * &
1226 ! ((Ri_crit - Ri) / (Ri_crit + CS%FRi_curvature*Ri))
1227  k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1228  ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1229  ke_src = k
1230  if (ks_src > k) ks_src = k
1231  else
1232  k_src(k) = 0.0
1233  endif
1234  enddo
1235 
1236  ! If there is no source anywhere, return kappa(K) = 0.
1237  if (ks_src > ke_src) then
1238  do k=1,nz+1
1239  kappa(k) = 0.0 ; k_q(k) = 0.0 ; tke(k) = tke_min
1240  enddo
1241  if (present(kappa_src)) then ; do k=1,nz+1 ; kappa_src(k) = 0.0 ; enddo ; endif
1242  if (present(local_src)) then ; do k=1,nz+1 ; local_src(k) = 0.0 ; enddo ; endif
1243  return
1244  endif
1245 
1246  do k=1,nz+1
1247  kappa(k) = kappa_in(k)
1248 ! TKE_decay(K) = c_n*sqrt(N2(K)) + c_s*sqrt(S2(K)) ! The expression in JHL.
1249  tke_decay(k) = sqrt(c_n2*n2(k) + c_s2*s2(k))
1250  if ((kappa(k) > 0.0) .and. (k_q(k) > 0.0)) then
1251  tke(k) = kappa(k) / k_q(k)
1252  else
1253  tke(k) = tke_min
1254  endif
1255  enddo
1256  ! Apply boundary conditions to kappa.
1257  kappa(1) = 0.0 ; kappa(nz+1) = 0.0
1258 
1259  ! Calculate the term (e1) that allows changes in TKE to be calculated quickly
1260  ! below the deepest nonzero value of kappa. If kappa = 0, below interface
1261  ! k-1, the final changes in TKE are related by dQ(K+1) = e1(K+1)*dQ(K).
1262  eden2 = kappa0 * idz(nz)
1263  if (tke_noflux_bottom_bc) then
1264  eden1 = dz_int(nz+1)*tke_decay(nz+1)
1265  i_eden = 1.0 / (eden2 + eden1)
1266  e1(nz+1) = eden2 * i_eden ; ome = eden1 * i_eden
1267  else
1268  e1(nz+1) = 0.0 ; ome = 1.0
1269  endif
1270  do k=nz,2,-1
1271  eden1 = dz_int(k)*tke_decay(k) + ome * eden2
1272  eden2 = kappa0 * idz(k-1)
1273  i_eden = 1.0 / (eden2 + eden1)
1274  e1(k) = eden2 * i_eden ; ome = eden1 * i_eden ! = 1-e1
1275  enddo
1276  e1(1) = 0.0
1277 
1278 
1279  ! Iterate here to convergence to within some tolerance of order tol_err.
1280  do itt=1,cs%max_RiNo_it
1281 
1282  ! ----------------------------------------------------
1283  ! Calculate TKE
1284  ! ----------------------------------------------------
1285 
1286 #ifdef DEBUG
1287  do k=1,nz+1 ; kappa_prev(k) = kappa(k) ; tke_prev(k) = tke(k) ; enddo
1288 #endif
1289 
1290  if (.not.do_newton) then
1291  ! Use separate steps of the TKE and kappa equations, that are
1292  ! explicit in the nonlinear source terms, implicit in a linearized
1293  ! version of the nonlinear sink terms, and implicit in the linear
1294  ! terms.
1295 
1296  ke_tke = max(ke_kappa,ke_kappa_prev)+1
1297  ! aQ is the coupling between adjacent interfaces in m s-1.
1298  do k=1,min(ke_tke,nz)
1299  aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
1300  enddo
1301  dq(1) = -tke(1)
1302  if (tke_noflux_top_bc) then
1303  tke_src = kappa0*s2(1) + q0 * tke_decay(1) ! Uses that kappa(1) = 0
1304  bd1 = dz_int(1) * tke_decay(1)
1305  bq = 1.0 / (bd1 + aq(1))
1306  tke(1) = bq * (dz_int(1)*tke_src)
1307  cq(2) = aq(1) * bq ; cqcomp = bd1 * bq ! = 1 - cQ
1308  else
1309  tke(1) = q0 ; cq(2) = 0.0 ; cqcomp = 1.0
1310  endif
1311  do k=2,ke_tke-1
1312  dq(k) = -tke(k)
1313  tke_src = (kappa(k) + kappa0)*s2(k) + q0*tke_decay(k)
1314  bd1 = dz_int(k)*(tke_decay(k) + n2(k)*k_q(k)) + cqcomp*aq(k-1)
1315  bq = 1.0 / (bd1 + aq(k))
1316  tke(k) = bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1))
1317  cq(k+1) = aq(k) * bq ; cqcomp = bd1 * bq ! = 1 - cQ
1318  enddo
1319  if ((ke_tke == nz+1) .and. .not.(tke_noflux_bottom_bc)) then
1320  tke(nz+1) = tke_min
1321  dq(nz+1) = 0.0
1322  else
1323  k = ke_tke
1324  tke_src = kappa0*s2(k) + q0*tke_decay(k) ! Uses that kappa(ke_tke) = 0
1325  if (k == nz+1) then
1326  dq(k) = -tke(k)
1327  bq = 1.0 / (dz_int(k)*tke_decay(k) + cqcomp*aq(k-1))
1328  tke(k) = max(tke_min, bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)))
1329  dq(k) = tke(k) + dq(k)
1330  else
1331  bq = 1.0 / ((dz_int(k)*tke_decay(k) + cqcomp*aq(k-1)) + aq(k))
1332  cq(k+1) = aq(k) * bq
1333  ! Account for all changes deeper in the water column.
1334  dq(k) = -tke(k)
1335  tke(k) = max((bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)) + &
1336  cq(k+1)*(tke(k+1) - e1(k+1)*tke(k))) / &
1337  (1.0 - cq(k+1)*e1(k+1)), tke_min)
1338  dq(k) = tke(k) + dq(k)
1339 
1340  ! Adjust TKE deeper in the water column in case ke_tke increases.
1341  ! This might not be strictly necessary?
1342  do k=ke_tke+1,nz+1
1343  dq(k) = e1(k)*dq(k-1)
1344  tke(k) = max(tke(k) + dq(k), tke_min)
1345  if (abs(dq(k)) < roundoff*tke(k)) exit
1346  enddo
1347  do k2=k+1,nz
1348  if (dq(k2) == 0.0) exit
1349  dq(k2) = 0.0
1350  enddo
1351  endif
1352  endif
1353  do k=ke_tke-1,1,-1
1354  tke(k) = max(tke(k) + cq(k+1)*tke(k+1), tke_min)
1355  dq(k) = tke(k) + dq(k)
1356  enddo
1357 
1358  ! ----------------------------------------------------
1359  ! Calculate kappa, here defined at interfaces.
1360  ! ----------------------------------------------------
1361 
1362  ke_kappa_prev = ke_kappa ; ks_kappa_prev = ks_kappa
1363 
1364  dk(1) = 0.0 ! kappa takes boundary values of 0.
1365  ck(2) = 0.0 ; ckcomp = 1.0
1366  if (itt == 1) then ; dO k=2,nz
1367  i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
1368  enddo ; endif
1369  do k=2,nz
1370  dk(k) = -kappa(k)
1371  if (itt>1) &
1372  i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
1373  bd1 = dz_int(k)*i_ld2(k) + ckcomp*idz(k-1)
1374  bk = 1.0 / (bd1 + idz(k))
1375 
1376  kappa(k) = bk * (idz(k-1)*kappa(k-1) + dz_int(k) * k_src(k))
1377  ck(k+1) = idz(k) * bk ; ckcomp = bd1 * bk ! = 1 - cK(K+1)
1378 
1379  ! Neglect values that are smaller than kappa_trunc.
1380  if (kappa(k) < ckcomp*kappa_trunc) then
1381  kappa(k) = 0.0
1382  if (k > ke_src) then ; ke_kappa = k-1 ; k_q(k) = 0.0 ; exit ; endif
1383  elseif (kappa(k) < 2.0*ckcomp*kappa_trunc) then
1384  kappa(k) = 2.0 * (kappa(k) - ckcomp*kappa_trunc)
1385  endif
1386  enddo
1387  k_q(ke_kappa) = kappa(ke_kappa) / tke(ke_kappa)
1388  dk(ke_kappa) = dk(ke_kappa) + kappa(ke_kappa)
1389  do k=ke_kappa+2,ke_kappa_prev
1390  dk(k) = -kappa(k) ; kappa(k) = 0.0 ; k_q(k) = 0.0
1391  enddo
1392  do k=ke_kappa-1,2,-1
1393  kappa(k) = kappa(k) + ck(k+1)*kappa(k+1)
1394  ! Neglect values that are smaller than kappa_trunc.
1395  if (kappa(k) <= kappa_trunc) then
1396  kappa(k) = 0.0
1397  if (k < ks_src) then ; ks_kappa = k+1 ; k_q(k) = 0.0 ; exit ; endif
1398  elseif (kappa(k) < 2.0*kappa_trunc) then
1399  kappa(k) = 2.0 * (kappa(k) - kappa_trunc)
1400  endif
1401 
1402  dk(k) = dk(k) + kappa(k)
1403  k_q(k) = kappa(k) / tke(k)
1404  enddo
1405  do k=ks_kappa_prev,ks_kappa-2 ; kappa(k) = 0.0 ; k_q(k) = 0.0 ; enddo
1406 
1407  else ! do_Newton is .true.
1408 ! Once the solutions are close enough, use a Newton's method solver of the
1409 ! whole system to accelerate convergence.
1410  ks_kappa_prev = ks_kappa ; ke_kappa_prev = ke_kappa ; ke_kappa = nz
1411  ks_kappa = 2
1412  dk(1) = 0.0 ; ck(2) = 0.0 ; ckcomp = 1.0 ; dkdq(1) = 0.0
1413  aq(1) = (0.5*(kappa(1)+kappa(2))+kappa0) * idz(1)
1414  dqdz(1) = 0.5*(tke(1) - tke(2))*idz(1)
1415  if (tke_noflux_top_bc) then
1416  tke_src = dz_int(1) * (kappa0*s2(1) - (tke(1) - q0)*tke_decay(1)) - &
1417  aq(1) * (tke(1) - tke(2))
1418 
1419  bq = 1.0 / (aq(1) + dz_int(1)*tke_decay(1))
1420  cq(2) = aq(1) * bq
1421  cqcomp = (dz_int(1)*tke_decay(1)) * bq ! = 1 - cQ(2)
1422  dqmdk(2) = -dqdz(1) * bq
1423  dq(1) = bq * tke_src
1424  else
1425  dq(1) = 0.0 ; cq(2) = 0.0 ; cqcomp = 1.0 ; dqmdk(2) = 0.0
1426  endif
1427  do k=2,nz
1428  i_q = 1.0 / tke(k)
1429  i_ld2(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
1430 
1431  kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa(k)) + &
1432  idz(k-1)*(kappa(k-1)-kappa(k)) - idz(k)*(kappa(k)-kappa(k+1))
1433 
1434  ! Ensure that the pivot is always positive, and that 0 <= cK <= 1.
1435  ! Otherwise do not use Newton's method.
1436  decay_term = -idz(k-1)*dqmdk(k)*dkdq(k-1) + dz_int(k)*i_ld2(k)
1437  if (decay_term < 0.0) then ; abort_newton = .true. ; exit ; endif
1438  bk = 1.0 / (idz(k) + idz(k-1)*ckcomp + decay_term)
1439 
1440  ck(k+1) = bk * idz(k)
1441  ckcomp = bk * (idz(k-1)*ckcomp + decay_term) ! = 1-cK(K+1)
1442  dkdq(k) = bk * (idz(k-1)*dkdq(k-1)*cq(k) + &
1443  (n2(k)*ilambda2 + f2)*i_q**2*kappa(k))
1444  dk(k) = bk * (kap_src + idz(k-1)*dk(k-1) + idz(k-1)*dkdq(k-1)*dq(k-1))
1445 
1446  ! Truncate away negligibly small values of kappa.
1447  if (dk(k) <= ckcomp*(kappa_trunc - kappa(k))) then
1448  dk(k) = -ckcomp*kappa(k)
1449 ! if (K > ke_src) then ; ke_kappa = k-1 ; K_Q(K) = 0.0 ; exit ; endif
1450  elseif (dk(k) < ckcomp*(2.0*kappa_trunc - kappa(k))) then
1451  dk(k) = 2.0 * dk(k) - ckcomp*(2.0*kappa_trunc - kappa(k))
1452  endif
1453 
1454  ! Solve for dQ(K)...
1455  aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
1456  dqdz(k) = 0.5*(tke(k) - tke(k+1))*idz(k)
1457  tke_src = dz_int(k) * ((kappa(k) + kappa0)*s2(k) - kappa(k)*n2(k) - &
1458  (tke(k) - q0)*tke_decay(k)) - &
1459  (aq(k) * (tke(k) - tke(k+1)) - aq(k-1) * (tke(k-1) - tke(k)))
1460  v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
1461  v2 = (v1*dqmdk(k) + dqdz(k-1)*ck(k)) + &
1462  ((dqdz(k-1) - dqdz(k)) + dz_int(k)*(s2(k) - n2(k)))
1463 
1464  ! Ensure that the pivot is always positive, and that 0 <= cQ <= 1.
1465  ! Otherwise do not use Newton's method.
1466  decay_term = dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k) - v2*dkdq(k)
1467  if (decay_term < 0.0) then ; abort_newton = .true. ; exit ; endif
1468  bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term))
1469 
1470  cq(k+1) = aq(k) * bq
1471  cqcomp = (cqcomp*aq(k-1) + decay_term) * bq
1472  dqmdk(k+1) = (v2 * ck(k+1) - dqdz(k)) * bq
1473 
1474  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1475  dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + &
1476  (v2 * dk(k) + tke_src)), cqcomp*(-0.5*tke(k)))
1477 
1478  ! Check whether the next layer will be affected by any nonzero kappas.
1479  if ((itt > 1) .and. (k > ke_src) .and. (dk(k) == 0.0) .and. &
1480  ((kappa(k) + kappa(k+1)) == 0.0)) then
1481  ! Could also do .and. (bQ*abs(tke_src) < roundoff*TKE(K)) then
1482  ke_kappa = k-1 ; exit
1483  endif
1484  enddo
1485  if ((ke_kappa == nz) .and. (.not. abort_newton)) then
1486  dk(nz+1) = 0.0 ; dkdq(nz+1) = 0.0
1487  if (tke_noflux_bottom_bc) then
1488  k = nz+1
1489  tke_src = dz_int(k) * (kappa0*s2(k) - (tke(k) - q0)*tke_decay(k)) + &
1490  aq(k-1) * (tke(k-1) - tke(k))
1491 
1492  v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
1493  decay_term = max(0.0, dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k))
1494  if (decay_term < 0.0) then
1495  abort_newton = .true.
1496  else
1497  bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term))
1498  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1499  dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + tke_src), &
1500  -0.5*tke(k))
1501  tke(k) = max(tke(k) + dq(k), tke_min)
1502  endif
1503  else
1504  dq(nz+1) = 0.0
1505  endif
1506  elseif (.not. abort_newton) then
1507  ! Alter the first-guess determination of dQ(K).
1508  dq(ke_kappa+1) = dq(ke_kappa+1) / (1.0 - cq(ke_kappa+2)*e1(ke_kappa+2))
1509  tke(ke_kappa+1) = max(tke(ke_kappa+1) + dq(ke_kappa+1), tke_min)
1510  do k=ke_kappa+2,nz+1
1511 #ifdef DEBUG
1512  if (k < nz+1) then
1513  ! Ignore this source?
1514  aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
1515  tke_src = (dz_int(k) * (kappa0*s2(k) - (tke(k)-q0)*tke_decay(k)) - &
1516  (aq(k) * (tke(k) - tke(k+1)) - &
1517  aq(k-1) * (tke(k-1) - tke(k))) ) / &
1518  (aq(k) + (aq(k-1) + dz_int(k)*tke_decay(k)))
1519  endif
1520 #endif
1521  dk(k) = 0.0
1522  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1523  dq(k) = max(e1(k)*dq(k-1),-0.5*tke(k))
1524  tke(k) = max(tke(k) + dq(k), tke_min)
1525  if (abs(dq(k)) < roundoff*tke(k)) exit
1526  enddo
1527 #ifdef DEBUG
1528  do k2=k+1,ke_kappa_prev+1 ; dq(k2) = 0.0 ; dk(k2) = 0.0 ; enddo
1529  do k=k2,nz+1 ; if (dq(k) == 0.0) exit ; dq(k) = 0.0 ; dk(k) = 0.0 ; enddo
1530 #endif
1531  endif
1532  if (.not. abort_newton) then
1533  do k=ke_kappa,2,-1
1534  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1535  dq(k) = max(dq(k) + (cq(k+1)*dq(k+1) + dqmdk(k+1) * dk(k+1)), &
1536  -0.5*tke(k))
1537  tke(k) = max(tke(k) + dq(k), tke_min)
1538  dk(k) = dk(k) + (ck(k+1)*dk(k+1) + dkdq(k) * dq(k))
1539  ! Truncate away negligibly small values of kappa.
1540  if (dk(k) <= kappa_trunc - kappa(k)) then
1541  dk(k) = -kappa(k)
1542  kappa(k) = 0.0
1543  if ((k < ks_src) .and. (k+1 > ks_kappa)) ks_kappa = k+1
1544  elseif (dk(k) < 2.0*kappa_trunc - kappa(k)) then
1545  dk(k) = 2.0*dk(k) - (2.0*kappa_trunc - kappa(k))
1546  kappa(k) = max(kappa(k) + dk(k), 0.0) ! The max is for paranoia.
1547  if (k<=ks_kappa) ks_kappa = 2
1548  else
1549  kappa(k) = kappa(k) + dk(k)
1550  if (k<=ks_kappa) ks_kappa = 2
1551  endif
1552  enddo
1553  dq(1) = max(dq(1) + cq(2)*dq(2) + dqmdk(2) * dk(2), tke_min - tke(1))
1554  tke(1) = max(tke(1) + dq(1), tke_min)
1555  dk(1) = 0.0
1556  endif
1557 
1558 #ifdef DEBUG
1559  ! Check these solutions for consistency.
1560  do k=2,nz
1561  ! In these equations, K_err_lin and Q_err_lin should be at round-off levels
1562  ! compared with the dominant terms, perhaps, dz_Int*I_Ld2*kappa and
1563  ! dz_Int*TKE_decay*TKE. The exception is where, either 1) the decay term has been
1564  ! been increased to ensure a positive pivot, or 2) negative TKEs have been
1565  ! truncated, or 3) small or negative kappas have been rounded toward 0.
1566  i_q = 1.0 / tke(k)
1567  i_ld2(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
1568 
1569  kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa_prev(k)) + &
1570  (idz(k-1)*(kappa_prev(k-1)-kappa_prev(k)) - &
1571  idz(k)*(kappa_prev(k)-kappa_prev(k+1)))
1572  k_err_lin = -idz(k-1)*(dk(k-1)-dk(k)) + idz(k)*(dk(k)-dk(k+1)) + &
1573  dz_int(k)*i_ld2(k)*dk(k) - kap_src - &
1574  (n2(k)*ilambda2 + f2)*i_q**2*kappa_prev(k) * dq(k)
1575 
1576  tke_src = dz_int(k) * ((kappa_prev(k) + kappa0)*s2(k) - &
1577  kappa_prev(k)*n2(k) - (tke_prev(k) - q0)*tke_decay(k)) - &
1578  (aq(k) * (tke_prev(k) - tke_prev(k+1)) - &
1579  aq(k-1) * (tke_prev(k-1) - tke_prev(k)))
1580  q_err_lin = (aq(k-1) * (dq(k-1)-dq(k)) - aq(k) * (dq(k)-dq(k+1))) - &
1581  0.5*(tke_prev(k)-tke_prev(k+1))*idz(k) * (dk(k) + dk(k+1)) - &
1582  0.5*(tke_prev(k)-tke_prev(k-1))*idz(k-1)* (dk(k-1) + dk(k)) + &
1583  dz_int(k) * (dk(k) * (s2(k) - n2(k)) - dq(k)*tke_decay(k)) + tke_src
1584  enddo
1585 #endif
1586  endif ! End of the Newton's method solver.
1587 
1588  ! Test kappa for convergence...
1589 #ifdef DEBUG
1590  max_err = 0.0 ; max_tke_err = 0.0 ; min_tke_err = 0.0
1591  do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
1592  norm_err = abs(kappa(k) - kappa_prev(k)) / &
1593  (kappa0 + 0.5*(kappa(k) + kappa_prev(k)))
1594  if (max_err < norm_err) max_err = norm_err
1595 
1596  tke_err(k) = dq(k) / (tke(k) - 0.5*dq(k))
1597  if (tke_err(k) > max_tke_err) max_tke_err = tke_err(k)
1598  if (tke_err(k) < min_tke_err) min_tke_err = tke_err(k)
1599  enddo
1600  if (do_newton) then
1601  if (max(max_err,max_tke_err,-min_tke_err) >= 2.0*newton_err) then
1602  do_newton = .false. ; abort_newton = .true.
1603  endif
1604  else
1605  if (max(max_err,max_tke_err,-min_tke_err) < newton_err) do_newton = .true.
1606  endif
1607  within_tolerance = (max_err < tol_err)
1608 #else
1609  ! max_err = 0.0
1610  if ((tol_err < newton_err) .and. (.not.abort_newton)) then
1611  ! A lower tolerance is used to switch to Newton's method than to
1612  ! switch back.
1613  newton_test = newton_err ; if (do_newton) newton_test = 2.0*newton_err
1614  was_newton = do_newton
1615  within_tolerance = .true. ; do_newton = .true.
1616  do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
1617  kappa_mean = kappa0 + (kappa(k) - 0.5*dk(k))
1618  if (abs(dk(k)) > newton_test * kappa_mean) then
1619  if (do_newton) abort_newton = .true.
1620  within_tolerance = .false. ; do_newton = .false. ; exit
1621  elseif (abs(dk(k)) > tol_err * kappa_mean) then
1622  within_tolerance = .false. ; if (.not.do_newton) exit
1623  endif
1624  if (abs(dq(k)) > newton_test*(tke(k) - 0.5*dq(k))) then
1625  if (do_newton) abort_newton = .true.
1626  do_newton = .false. ; if (.not.within_tolerance) exit
1627  endif
1628  enddo
1629 
1630  else ! Newton's method will not be used again, so no need to check.
1631  within_tolerance = .true.
1632  do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
1633  if (abs(dk(k)) > tol_err * (kappa0 + (kappa(k) - 0.5*dk(k)))) then
1634  within_tolerance = .false. ; exit
1635  endif
1636  enddo
1637  endif
1638 #endif
1639 
1640  if (abort_newton) then
1641  do_newton = .false. ; abort_newton = .false.
1642  ! We went to Newton too quickly last time, so restrict the tolerance.
1643  newton_err = 0.5*newton_err
1644  ke_kappa_prev = nz
1645  do k=2,nz ; k_q(k) = kappa(k) / max(tke(k), tke_min) ; enddo
1646  endif
1647 
1648 #ifdef DEBUG
1649  if (itt <= max_debug_itt) then
1650  do k=1,nz+1
1651  kprev_it1(k,itt)=kappa_prev(k)
1652  kappa_it1(k,itt)=kappa(k) ; tke_it1(k,itt) = tke(k)
1653  dkappa_it1(k,itt) = kappa(k) - kappa_prev(k)
1654  dkappa_norm_it1(k,itt) = (kappa(k) - kappa_prev(k)) / &
1655  (kappa0 + 0.5*(kappa(k) + kappa_prev(k)))
1656  k_q_it1(k,itt) = kappa(k) / max(tke(k),tke_min)
1657  d_dkappa_it1(k,itt) = 0.0
1658  if (itt > 1) then ; if (abs(dkappa_it1(k,itt-1)) > 1e-20) &
1659  d_dkappa_it1(k,itt) = dkappa_it1(k,itt) / dkappa_it1(k,itt-1)
1660  endif
1661  enddo
1662  endif
1663 #endif
1664 
1665  if (within_tolerance) exit
1666 
1667  enddo
1668 
1669 #ifdef DEBUG
1670  do it2=itt+1,max_debug_itt ; do k=1,nz+1
1671  kprev_it1(k,it2) = 0.0 ; kappa_it1(k,it2) = 0.0 ; tke_it1(k,it2) = 0.0
1672  dkappa_it1(k,it2) = 0.0 ; k_q_it1(k,it2) = 0.0 ; d_dkappa_it1(k,it2) = 0.0
1673  enddo ; enddo
1674 #endif
1675 
1676  if (do_newton) then ! K_Q needs to be calculated.
1677  do k=1,ks_kappa-1 ; k_q(k) = 0.0 ; enddo
1678  do k=ks_kappa,ke_kappa ; k_q(k) = kappa(k) / tke(k) ; enddo
1679  do k=ke_kappa+1,nz+1 ; k_q(k) = 0.0 ; enddo
1680  endif
1681 
1682  if (present(local_src)) then
1683  local_src(1) = 0.0 ; local_src(nz+1) = 0.0
1684  do k=2,nz
1685  diffusive_src = idz(k-1)*(kappa(k-1)-kappa(k)) + &
1686  idz(k)*(kappa(k+1)-kappa(k))
1687  chg_by_k0 = kappa0 * ((idz(k-1)+idz(k)) / dz_int(k) + i_ld2(k))
1688  if (diffusive_src <= 0.0) then
1689  local_src(k) = k_src(k) + chg_by_k0
1690  else
1691  local_src(k) = (k_src(k) + chg_by_k0) + diffusive_src / dz_int(k)
1692  endif
1693  enddo
1694  endif
1695  if (present(kappa_src)) then
1696  kappa_src(1) = 0.0 ; kappa_src(nz+1) = 0.0
1697  do k=2,nz
1698  kappa_src(k) = k_src(k)
1699  enddo
1700  endif
1701 
Here is the caller graph for this function:

◆ kappa_shear_init()

logical function, public mom_kappa_shear::kappa_shear_init ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(kappa_shear_cs), pointer  CS 
)
Parameters
[in]timeThe current model time.
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]param_fileA structure to parse for run-time parameters.
[in,out]diagA structure that is used to regulate diagnostic output.
csA pointer that is set to point to the control structure for this module

Definition at line 1706 of file MOM_kappa_shear.F90.

References mdl, and mom_error_handler::mom_error().

1706  type(time_type), intent(in) :: time !< The current model time.
1707  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1708  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1709  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1710  !! parameters.
1711  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
1712  !! output.
1713  type(kappa_shear_cs), pointer :: cs !< A pointer that is set to point to the control
1714  !! structure for this module
1715 ! Arguments: Time - The current model time.
1716 ! (in) G - The ocean's grid structure.
1717 ! (in) GV - The ocean's vertical grid structure.
1718 ! (in) param_file - A structure indicating the open file to parse for
1719 ! model parameter values.
1720 ! (in) diag - A structure that is used to regulate diagnostic output.
1721 ! (in/out) CS - A pointer that is set to point to the control structure
1722 ! for this module
1723 ! (returns) kappa_shear_init - True if module is to be used, False otherwise
1724  logical :: merge_mixedlayer
1725 ! This include declares and sets the variable "version".
1726 #include "version_variable.h"
1727  real :: kd_normal ! The KD of the main model, read here only as a parameter
1728  ! for setting the default of KD_SMOOTH
1729  if (associated(cs)) then
1730  call mom_error(warning, "kappa_shear_init called with an associated "// &
1731  "control structure.")
1732  return
1733  endif
1734  allocate(cs)
1735 
1736  ! The Jackson-Hallberg-Legg shear mixing parameterization uses the following
1737  ! 6 nondimensional coefficients. That paper gives 3 best fit parameter sets.
1738  ! Ri_Crit Rate FRi_Curv K_buoy TKE_N TKE_Shear
1739  ! p1: 0.25 0.089 -0.97 0.82 0.24 0.14
1740  ! p2: 0.30 0.085 -0.94 0.86 0.26 0.13
1741  ! p3: 0.35 0.088 -0.89 0.81 0.28 0.12
1742  ! Future research will reveal how these should be modified to take
1743  ! subgridscale inhomogeneity into account.
1744 
1745 ! Set default, read and log parameters
1746  call log_version(param_file, mdl, version, &
1747  "Parameterization of shear-driven turbulence following Jackson, Hallberg and Legg, JPO 2008")
1748  call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_init, &
1749  "If true, use the Jackson-Hallberg-Legg (JPO 2008) \n"//&
1750  "shear mixing parameterization.", default=.false.)
1751  call get_param(param_file, mdl, "RINO_CRIT", cs%RiNo_crit, &
1752  "The critical Richardson number for shear mixing.", &
1753  units="nondim", default=0.25)
1754  call get_param(param_file, mdl, "SHEARMIX_RATE", cs%Shearmix_rate, &
1755  "A nondimensional rate scale for shear-driven entrainment.\n"//&
1756  "Jackson et al find values in the range of 0.085-0.089.", &
1757  units="nondim", default=0.089)
1758  call get_param(param_file, mdl, "MAX_RINO_IT", cs%max_RiNo_it, &
1759  "The maximum number of iterations that may be used to \n"//&
1760  "estimate the Richardson number driven mixing.", &
1761  units="nondim", default=50)
1762  call get_param(param_file, mdl, "KD", kd_normal, default=1.0e-7, do_not_log=.true.)
1763  call get_param(param_file, mdl, "KD_KAPPA_SHEAR_0", cs%kappa_0, &
1764  "The background diffusivity that is used to smooth the \n"//&
1765  "density and shear profiles before solving for the \n"//&
1766  "diffusivities. Defaults to value of KD.", units="m2 s-1", default=kd_normal)
1767  call get_param(param_file, mdl, "FRI_CURVATURE", cs%FRi_curvature, &
1768  "The nondimensional curvature of the function of the \n"//&
1769  "Richardson number in the kappa source term in the \n"//&
1770  "Jackson et al. scheme.", units="nondim", default=-0.97)
1771  call get_param(param_file, mdl, "TKE_N_DECAY_CONST", cs%C_N, &
1772  "The coefficient for the decay of TKE due to \n"//&
1773  "stratification (i.e. proportional to N*tke). \n"//&
1774  "The values found by Jackson et al. are 0.24-0.28.", &
1775  units="nondim", default=0.24)
1776 ! call get_param(param_file, mdl, "LAYER_KAPPA_STAGGER", CS%layer_stagger, &
1777 ! default=.false.)
1778  call get_param(param_file, mdl, "TKE_SHEAR_DECAY_CONST", cs%C_S, &
1779  "The coefficient for the decay of TKE due to shear (i.e. \n"//&
1780  "proportional to |S|*tke). The values found by Jackson \n"//&
1781  "et al. are 0.14-0.12.", units="nondim", default=0.14)
1782  call get_param(param_file, mdl, "KAPPA_BUOY_SCALE_COEF", cs%lambda, &
1783  "The coefficient for the buoyancy length scale in the \n"//&
1784  "kappa equation. The values found by Jackson et al. are \n"//&
1785  "in the range of 0.81-0.86.", units="nondim", default=0.82)
1786  call get_param(param_file, mdl, "KAPPA_N_OVER_S_SCALE_COEF2", cs%lambda2_N_S, &
1787  "The square of the ratio of the coefficients of the \n"//&
1788  "buoyancy and shear scales in the diffusivity equation, \n"//&
1789  "Set this to 0 (the default) to eliminate the shear scale. \n"//&
1790  "This is only used if USE_JACKSON_PARAM is true.", &
1791  units="nondim", default=0.0)
1792  call get_param(param_file, mdl, "KAPPA_SHEAR_TOL_ERR", cs%kappa_tol_err, &
1793  "The fractional error in kappa that is tolerated. \n"//&
1794  "Iteration stops when changes between subsequent \n"//&
1795  "iterations are smaller than this everywhere in a \n"//&
1796  "column. The peak diffusivities usually converge most \n"//&
1797  "rapidly, and have much smaller errors than this.", &
1798  units="nondim", default=0.1)
1799  call get_param(param_file, mdl, "TKE_BACKGROUND", cs%TKE_bg, &
1800  "A background level of TKE used in the first iteration \n"//&
1801  "of the kappa equation. TKE_BACKGROUND could be 0.", &
1802  units="m2 s-2", default=0.0)
1803  call get_param(param_file, mdl, "KAPPA_SHEAR_ELIM_MASSLESS", cs%eliminate_massless, &
1804  "If true, massless layers are merged with neighboring \n"//&
1805  "massive layers in this calculation. The default is \n"//&
1806  "true and I can think of no good reason why it should \n"//&
1807  "be false. This is only used if USE_JACKSON_PARAM is true.", &
1808  default=.true.)
1809  call get_param(param_file, mdl, "MAX_KAPPA_SHEAR_IT", cs%max_KS_it, &
1810  "The maximum number of iterations that may be used to \n"//&
1811  "estimate the time-averaged diffusivity.", units="nondim", &
1812  default=13)
1813  call get_param(param_file, mdl, "PRANDTL_TURB", cs%Prandtl_turb, &
1814  "The turbulent Prandtl number applied to shear \n"//&
1815  "instability.", units="nondim", default=1.0, do_not_log=.true.)
1816  call get_param(param_file, mdl, "DEBUG_KAPPA_SHEAR", cs%debug, &
1817  "If true, write debugging data for the kappa-shear code. \n"//&
1818  "Caution: this option is _very_ verbose and should only \n"//&
1819  "be used in single-column mode!", default=.false.)
1820 
1821 ! id_clock_KQ = cpu_clock_id('Ocean KS kappa_shear',grain=CLOCK_ROUTINE)
1822 ! id_clock_avg = cpu_clock_id('Ocean KS avg',grain=CLOCK_ROUTINE)
1823 ! id_clock_project = cpu_clock_id('Ocean KS project',grain=CLOCK_ROUTINE)
1824 ! id_clock_setup = cpu_clock_id('Ocean KS setup',grain=CLOCK_ROUTINE)
1825 
1826  cs%nkml = 1
1827  if (gv%nkml>0) then
1828  call get_param(param_file, mdl, "KAPPA_SHEAR_MERGE_ML",merge_mixedlayer, &
1829  "If true, combine the mixed layers together before \n"//&
1830  "solving the kappa-shear equations.", default=.true.)
1831  if (merge_mixedlayer) cs%nkml = gv%nkml
1832  endif
1833 
1834 ! Forego remainder of initialization if not using this scheme
1835  if (.not. kappa_shear_init) return
1836 
1837  cs%diag => diag
1838 
1839  cs%id_Kd_shear = register_diag_field('ocean_model','Kd_shear',diag%axesTi,time, &
1840  'Shear-driven Diapycnal Diffusivity', 'meter2 second-1')
1841  cs%id_TKE = register_diag_field('ocean_model','TKE_shear',diag%axesTi,time, &
1842  'Shear-driven Turbulent Kinetic Energy', 'meter2 second-2')
1843 #ifdef ADD_DIAGNOSTICS
1844  cs%id_ILd2 = register_diag_field('ocean_model','ILd2_shear',diag%axesTi,time, &
1845  'Inverse kappa decay scale at interfaces', 'meter-2')
1846  cs%id_dz_Int = register_diag_field('ocean_model','dz_Int_shear',diag%axesTi,time, &
1847  'Finite volume thickness of interfaces', 'meter')
1848 #endif
1849 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:

◆ kappa_shear_is_used()

logical function, public mom_kappa_shear::kappa_shear_is_used ( type(param_file_type), intent(in)  param_file)
Parameters
[in]param_fileA structure to parse for run-time parameters

Definition at line 1853 of file MOM_kappa_shear.F90.

References mdl.

Referenced by mom_cvmix_shear::cvmix_shear_init(), mom_diabatic_driver::diabatic_driver_init(), mom_set_visc::set_visc_init(), and mom_set_visc::set_visc_register_restarts().

1853 ! Reads the parameter "USE_JACKSON_PARAM" and returns state.
1854 ! This function allows other modules to know whether this parameterization will
1855 ! be used without needing to duplicate the log entry.
1856  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1857  call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_is_used, &
1858  default=.false., do_not_log = .true.)
Here is the caller graph for this function:

Variable Documentation

◆ mdl

character(len=40) mom_kappa_shear::mdl = "MOM_kappa_shear"

Definition at line 112 of file MOM_kappa_shear.F90.

Referenced by kappa_shear_init(), and kappa_shear_is_used().

112  character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.