MOM6
Data Types | Functions/Subroutines
mom_diabatic_driver Module Reference

Detailed Description

This routine drives the diabatic/dianeutral physics for MOM.

By Robert Hallberg, Alistair Adcroft, and Stephen Griffies.

This program contains the subroutine that, along with the subroutines that it calls, implements diapycnal mass and momentum fluxes and a bulk mixed layer. The diapycnal diffusion can be used without the bulk mixed layer.

Outline of MOM diabatic

  • diabatic first determines the (diffusive) diapycnal mass fluxes based on the convergence of the buoyancy fluxes within each layer.
  • The dual-stream entrainment scheme of MacDougall and Dewar (JPO, 1997) is used for combined diapycnal advection and diffusion, calculated implicitly and potentially with the Richardson number dependent mixing, as described by Hallberg (MWR, 2000).
  • Diapycnal advection is the residual of diapycnal diffusion, so the fully implicit upwind differencing scheme that is used is entirely appropriate.

  • The downward buoyancy flux in each layer is determined from an implicit calculation based on the previously calculated flux of the layer above and an estimated flux in the layer below. This flux is subject to the following conditions: (1) the flux in the top and bottom layers are set by the boundary conditions, and (2) no layer may be driven below an Angstrom thick- ness. If there is a bulk mixed layer, the buffer layer is treated as a fixed density layer with vanishingly small diffusivity.

    diabatic takes 5 arguments: the two velocities (u and v), the thicknesses (h), a structure containing the forcing fields, and the length of time over which to act (dt). The velocities and thickness are taken as inputs and modified within the subroutine. There is no limit on the time step.

Data Types

type  diabatic_cs
 Control structure for this module. More...
 

Functions/Subroutines

subroutine, public diabatic (u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, GV, US, CS, OBC, WAVES)
 This subroutine imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers. More...
 
subroutine diabatic_ale_legacy (u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, GV, US, CS, Waves)
 Applies diabatic forcing and diapycnal mixing of temperature, salinity and other tracers for use with an ALE algorithm. This version uses an older set of algorithms compared with diabatic_ALE. More...
 
subroutine diabatic_ale (u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, GV, US, CS, Waves)
 This subroutine imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers. More...
 
subroutine layered_diabatic (u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, GV, US, CS, WAVES)
 Imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers using the original MOM6 algorithms. More...
 
subroutine, public extract_diabatic_member (CS, opacity_CSp, optics_CSp, evap_CFL_limit, minimum_forcing_depth, KPP_CSp, energetic_PBL_CSp, diabatic_aux_CSp, diabatic_halo)
 Returns pointers or values of members within the diabatic_CS type. For extensibility, each returned argument is an optional argument. More...
 
subroutine, public adiabatic (h, tv, fluxes, dt, G, GV, US, CS)
 Routine called for adiabatic physics. More...
 
subroutine diagnose_diabatic_diff_tendency (tv, h, temp_old, saln_old, dt, G, GV, US, CS)
 This routine diagnoses tendencies from application of diabatic diffusion using ALE algorithm. Note that layer thickness is not altered by diabatic diffusion. More...
 
subroutine diagnose_boundary_forcing_tendency (tv, h, temp_old, saln_old, h_old, dt, G, GV, US, CS)
 This routine diagnoses tendencies from application of boundary fluxes. These impacts are generally 3d, in particular for penetrative shortwave. Other fluxes contribute 3d in cases when the layers vanish or are very thin, in which case we distribute the flux into k > 1 layers. More...
 
subroutine diagnose_frazil_tendency (tv, h, temp_old, dt, G, GV, US, CS)
 This routine diagnoses tendencies for temperature and heat from frazil formation. This routine is called twice from within subroutine diabatic; at start and at end of the diabatic processes. The impacts from frazil are generally a function of depth. Hence, when checking heat budget, be sure to remove HFSIFRAZIL from HFDS in k=1. More...
 
subroutine, public adiabatic_driver_init (Time, G, param_file, diag, CS, tracer_flow_CSp)
 A simplified version of diabatic_driver_init that will allow tracer column functions to be called without allowing any of the diabatic processes to be used. More...
 
subroutine, public diabatic_driver_init (Time, G, GV, US, param_file, useALEalgorithm, diag, ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, ALE_sponge_CSp)
 This routine initializes the diabatic driver module. More...
 
subroutine, public diabatic_driver_end (CS)
 Routine to close the diabatic driver module. More...
 

Variables

integer id_clock_entrain
 clock ids
 
integer id_clock_mixedlayer
 clock ids
 
integer id_clock_set_diffusivity
 clock ids
 
integer id_clock_tracers
 clock ids
 
integer id_clock_tridiag
 clock ids
 
integer id_clock_pass
 clock ids
 
integer id_clock_sponge
 clock ids
 
integer id_clock_geothermal
 clock ids
 
integer id_clock_differential_diff
 clock ids
 
integer id_clock_remap
 clock ids
 
integer id_clock_kpp
 clock ids
 

Function/Subroutine Documentation

◆ adiabatic()

subroutine, public mom_diabatic_driver::adiabatic ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  h,
type(thermo_var_ptrs), intent(inout)  tv,
type(forcing), intent(inout)  fluxes,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diabatic_cs), pointer  CS 
)

Routine called for adiabatic physics.

Parameters
[in,out]gocean grid structure
[in,out]hthickness [H ~> m or kg m-2]
[in,out]tvpoints to thermodynamic fields
[in,out]fluxesboundary fluxes
[in]dttime step [T ~> s]
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
csmodule control structure

Definition at line 2849 of file MOM_diabatic_driver.F90.

2849  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
2850  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2851  intent(inout) :: h !< thickness [H ~> m or kg m-2]
2852  type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
2853  type(forcing), intent(inout) :: fluxes !< boundary fluxes
2854  real, intent(in) :: dt !< time step [T ~> s]
2855  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2856  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2857  type(diabatic_CS), pointer :: CS !< module control structure
2858 
2859  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: zeros ! An array of zeros.
2860 
2861  zeros(:,:,:) = 0.0
2862 
2863  call call_tracer_column_fns(h, h, zeros, zeros, fluxes, zeros(:,:,1), dt, g, gv, us, tv, &
2864  cs%optics, cs%tracer_flow_CSp, cs%debug)
2865 

◆ adiabatic_driver_init()

subroutine, public mom_diabatic_driver::adiabatic_driver_init ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(diabatic_cs), pointer  CS,
type(tracer_flow_control_cs), pointer  tracer_flow_CSp 
)

A simplified version of diabatic_driver_init that will allow tracer column functions to be called without allowing any of the diabatic processes to be used.

Parameters
[in]timecurrent model time
[in]gmodel grid structure
[in]param_filethe file to parse for parameter values
[in,out]diagregulates diagnostic output
csmodule control structure
tracer_flow_csppointer to control structure of the tracer flow control module

Definition at line 3115 of file MOM_diabatic_driver.F90.

3115  type(time_type), intent(in) :: Time !< current model time
3116  type(ocean_grid_type), intent(in) :: G !< model grid structure
3117  type(param_file_type), intent(in) :: param_file !< the file to parse for parameter values
3118  type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
3119  type(diabatic_CS), pointer :: CS !< module control structure
3120  type(tracer_flow_control_CS), pointer :: tracer_flow_CSp !< pointer to control structure of the
3121  !! tracer flow control module
3122 
3123 ! This "include" declares and sets the variable "version".
3124 #include "version_variable.h"
3125  character(len=40) :: mdl = "MOM_diabatic_driver" ! This module's name.
3126 
3127  if (associated(cs)) then
3128  call mom_error(warning, "adiabatic_driver_init called with an "// &
3129  "associated control structure.")
3130  return
3131  else ; allocate(cs) ; endif
3132 
3133  cs%diag => diag
3134  if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
3135 
3136 ! Set default, read and log parameters
3137  call log_version(param_file, mdl, version, &
3138  "The following parameters are used for diabatic processes.")
3139 

◆ diabatic()

subroutine, public mom_diabatic_driver::diabatic ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(inout)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h,
type(thermo_var_ptrs), intent(inout)  tv,
real, dimension(:,:), pointer  Hml,
type(forcing), intent(inout)  fluxes,
type(vertvisc_type), intent(inout)  visc,
type(accel_diag_ptrs), intent(inout)  ADp,
type(cont_diag_ptrs), intent(inout)  CDp,
real, intent(in)  dt,
type(time_type), intent(in)  Time_end,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diabatic_cs), pointer  CS,
type(ocean_obc_type), optional, pointer  OBC,
type(wave_parameters_cs), optional, pointer  WAVES 
)

This subroutine imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers.

Parameters
[in,out]gocean grid structure
[in]gvocean vertical grid structure
[in,out]uzonal velocity [L T-1 ~> m s-1]
[in,out]vmeridional velocity [L T-1 ~> m s-1]
[in,out]hthickness [H ~> m or kg m-2]
[in,out]tvpoints to thermodynamic fields unused have NULL ptrs
hmlActive mixed layer depth [Z ~> m]
[in,out]fluxespoints to forcing fields unused fields have NULL ptrs
[in,out]viscvertical viscosities, BBL properies, and
[in,out]adprelated points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]cdppoints to terms in continuity equations
[in]dttime increment [T ~> s]
[in]time_endTime at the end of the interval
[in]usA dimensional unit scaling type
csmodule control structure
obcOpen boundaries control structure.
wavesSurface gravity waves

Definition at line 263 of file MOM_diabatic_driver.F90.

263  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
264  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
265  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
266  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
267  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2]
268  type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
269  !! unused have NULL ptrs
270  real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m]
271  type(forcing), intent(inout) :: fluxes !< points to forcing fields
272  !! unused fields have NULL ptrs
273  type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and
274  type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum
275  !! equations, to enable the later derived
276  !! diagnostics, like energy budgets
277  type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations
278  real, intent(in) :: dt !< time increment [T ~> s]
279  type(time_type), intent(in) :: Time_end !< Time at the end of the interval
280  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
281  type(diabatic_CS), pointer :: CS !< module control structure
282  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
283  type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves
284 
285  ! local variables
286  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
287  eta ! Interface heights before diapycnal mixing [m].
288  real, dimension(SZI_(G),SZJ_(G),CS%nMode) :: &
289  cn_IGW ! baroclinic internal gravity wave speeds
290  real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag ! diagnostic array for temp
291  integer :: i, j, k, m, is, ie, js, je, nz
292  logical :: showCallTree ! If true, show the call tree
293 
294  if (g%ke == 1) return
295 
296  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
297 
298  if (.not. associated(cs)) call mom_error(fatal, "MOM_diabatic_driver: "// &
299  "Module must be initialized before it is used.")
300  if (dt == 0.0) call mom_error(fatal, "MOM_diabatic_driver: "// &
301  "diabatic was called with a zero length timestep.")
302  if (dt < 0.0) call mom_error(fatal, "MOM_diabatic_driver: "// &
303  "diabatic was called with a negative timestep.")
304 
305  showcalltree = calltree_showquery()
306 
307  ! Offer diagnostics of various state varables at the start of diabatic
308  ! these are mostly for debugging purposes.
309  if (cs%id_u_predia > 0) call post_data(cs%id_u_predia, u, cs%diag)
310  if (cs%id_v_predia > 0) call post_data(cs%id_v_predia, v, cs%diag)
311  if (cs%id_h_predia > 0) call post_data(cs%id_h_predia, h, cs%diag)
312  if (cs%id_T_predia > 0) call post_data(cs%id_T_predia, tv%T, cs%diag)
313  if (cs%id_S_predia > 0) call post_data(cs%id_S_predia, tv%S, cs%diag)
314  if (cs%id_e_predia > 0) then
315  call find_eta(h, tv, g, gv, us, eta, eta_to_m=1.0)
316  call post_data(cs%id_e_predia, eta, cs%diag)
317  endif
318 
319  if (cs%debug) then
320  call mom_state_chksum("Start of diabatic ", u, v, h, g, gv, us, haloshift=0)
321  call mom_forcing_chksum("Start of diabatic", fluxes, g, us, haloshift=0)
322  endif
323  if (cs%debugConservation) call mom_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, g, gv, us)
324 
325  if (cs%debug_energy_req) &
326  call diapyc_energy_req_test(h, dt, tv, g, gv, us, cs%diapyc_en_rec_CSp)
327 
328  call cpu_clock_begin(id_clock_set_diffusivity)
329  call set_bbl_tke(u, v, h, fluxes, visc, g, gv, us, cs%set_diff_CSp, obc=obc)
330  call cpu_clock_end(id_clock_set_diffusivity)
331 
332  ! Frazil formation keeps the temperature above the freezing point.
333  ! make_frazil is deliberately called at both the beginning and at
334  ! the end of the diabatic processes.
335  if (associated(tv%T) .AND. associated(tv%frazil)) then
336  ! For frazil diagnostic, the first call covers the first half of the time step
337  call enable_averages(0.5*dt, time_end - real_to_time(0.5*us%T_to_s*dt), cs%diag)
338  if (cs%frazil_tendency_diag) then
339  do k=1,nz ; do j=js,je ; do i=is,ie
340  temp_diag(i,j,k) = tv%T(i,j,k)
341  enddo ; enddo ; enddo
342  endif
343 
344  if (associated(fluxes%p_surf_full)) then
345  call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp, fluxes%p_surf_full, halo=cs%halo_TS_diff)
346  else
347  call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp, halo=cs%halo_TS_diff)
348  endif
349  if (showcalltree) call calltree_waypoint("done with 1st make_frazil (diabatic)")
350 
351  if (cs%frazil_tendency_diag) then
352  call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, g, gv, us, cs)
353  if (cs%id_frazil_h > 0) call post_data(cs%id_frazil_h, h, cs%diag)
354  endif
355  call disable_averaging(cs%diag)
356  endif ! associated(tv%T) .AND. associated(tv%frazil)
357  if (cs%debugConservation) call mom_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, g, gv, us)
358 
359  if (cs%use_int_tides) then
360  ! This block provides an interface for the unresolved low-mode internal tide module.
361  call set_int_tide_input(u, v, h, tv, fluxes, cs%int_tide_input, dt, g, gv, us, &
362  cs%int_tide_input_CSp)
363  cn_igw(:,:,:) = 0.0
364  if (cs%uniform_test_cg > 0.0) then
365  do m=1,cs%nMode ; cn_igw(:,:,m) = cs%uniform_test_cg ; enddo
366  else
367  call wave_speeds(h, tv, g, gv, us, cs%nMode, cn_igw, full_halos=.true.)
368  endif
369 
370  call propagate_int_tide(h, tv, cn_igw, cs%int_tide_input%TKE_itidal_input, cs%int_tide_input%tideamp, &
371  cs%int_tide_input%Nb, dt, g, gv, us, cs%int_tide_CSp)
372  if (showcalltree) call calltree_waypoint("done with propagate_int_tide (diabatic)")
373  endif ! end CS%use_int_tides
374 
375  if (cs%useALEalgorithm .and. cs%use_legacy_diabatic) then
376  call diabatic_ale_legacy(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
377  g, gv, us, cs, waves)
378  elseif (cs%useALEalgorithm) then
379  call diabatic_ale(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
380  g, gv, us, cs, waves)
381  else
382  call layered_diabatic(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
383  g, gv, us, cs, waves)
384  endif
385 
386 
387  call cpu_clock_begin(id_clock_pass)
388  if (associated(visc%Kv_shear)) &
389  call pass_var(visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)
390  call cpu_clock_end(id_clock_pass)
391 
392  call disable_averaging(cs%diag)
393  ! Frazil formation keeps temperature above the freezing point.
394  ! make_frazil is deliberately called at both the beginning and at
395  ! the end of the diabatic processes.
396  if (associated(tv%T) .AND. associated(tv%frazil)) then
397  call enable_averages(0.5*dt, time_end, cs%diag)
398  if (cs%frazil_tendency_diag) then
399  do k=1,nz ; do j=js,je ; do i=is,ie
400  temp_diag(i,j,k) = tv%T(i,j,k)
401  enddo ; enddo ; enddo
402  endif
403 
404  if (associated(fluxes%p_surf_full)) then
405  call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp, fluxes%p_surf_full)
406  else
407  call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp)
408  endif
409 
410  if (cs%frazil_tendency_diag) then
411  call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, g, gv, us, cs)
412  if (cs%id_frazil_h > 0 ) call post_data(cs%id_frazil_h, h, cs%diag)
413  endif
414 
415  if (showcalltree) call calltree_waypoint("done with 2nd make_frazil (diabatic)")
416  if (cs%debugConservation) call mom_state_stats('2nd make_frazil', u, v, h, tv%T, tv%S, g, gv, us)
417  call disable_averaging(cs%diag)
418 
419  endif ! endif for frazil
420 
421 
422  ! Diagnose mixed layer depths.
423  call enable_averages(dt, time_end, cs%diag)
424  if (cs%id_MLD_003 > 0 .or. cs%id_subMLN2 > 0 .or. cs%id_mlotstsq > 0) then
425  call diagnosemldbydensitydifference(cs%id_MLD_003, h, tv, 0.03*us%kg_m3_to_R, g, gv, us, cs%diag, &
426  id_n2subml=cs%id_subMLN2, id_mldsq=cs%id_mlotstsq, dz_subml=cs%dz_subML_N2)
427  endif
428  if (cs%id_MLD_0125 > 0) then
429  call diagnosemldbydensitydifference(cs%id_MLD_0125, h, tv, 0.125*us%kg_m3_to_R, g, gv, us, cs%diag)
430  endif
431  if (cs%id_MLD_user > 0) then
432  call diagnosemldbydensitydifference(cs%id_MLD_user, h, tv, cs%MLDdensityDifference, g, gv, us, cs%diag)
433  endif
434  if ((cs%id_MLD_EN1 > 0) .or. (cs%id_MLD_EN2 > 0) .or. (cs%id_MLD_EN3 > 0)) then
435  call diagnosemldbyenergy((/cs%id_MLD_EN1, cs%id_MLD_EN2, cs%id_MLD_EN3/),&
436  h, tv, g, gv, us, cs%MLD_EN_VALS, cs%diag)
437  endif
438  if (cs%use_int_tides) then
439  if (cs%id_cg1 > 0) call post_data(cs%id_cg1, cn_igw(:,:,1),cs%diag)
440  do m=1,cs%nMode ; if (cs%id_cn(m) > 0) call post_data(cs%id_cn(m), cn_igw(:,:,m), cs%diag) ; enddo
441  endif
442  call disable_averaging(cs%diag)
443 
444  if (cs%debugConservation) call mom_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, g, gv, us)
445 

◆ diabatic_ale()

subroutine mom_diabatic_driver::diabatic_ale ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(inout)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h,
type(thermo_var_ptrs), intent(inout)  tv,
real, dimension(:,:), pointer  Hml,
type(forcing), intent(inout)  fluxes,
type(vertvisc_type), intent(inout)  visc,
type(accel_diag_ptrs), intent(inout)  ADp,
type(cont_diag_ptrs), intent(inout)  CDp,
real, intent(in)  dt,
type(time_type), intent(in)  Time_end,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diabatic_cs), pointer  CS,
type(wave_parameters_cs), optional, pointer  Waves 
)
private

This subroutine imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers.

Parameters
[in,out]gocean grid structure
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]uzonal velocity [L T-1 ~> m s-1]
[in,out]vmeridional velocity [L T-1 ~> m s-1]
[in,out]hthickness [H ~> m or kg m-2]
[in,out]tvpoints to thermodynamic fields unused have NULL ptrs
hmlActive mixed layer depth [Z ~> m]
[in,out]fluxespoints to forcing fields unused fields have NULL ptrs
[in,out]viscvertical viscosities, BBL properies, and
[in,out]adprelated points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]cdppoints to terms in continuity equations
[in]dttime increment [T ~> s]
[in]time_endTime at the end of the interval
csmodule control structure
wavesSurface gravity waves

Definition at line 1173 of file MOM_diabatic_driver.F90.

1173  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1174  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1175  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1176  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1177  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
1178  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2]
1179  type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
1180  !! unused have NULL ptrs
1181  real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m]
1182  type(forcing), intent(inout) :: fluxes !< points to forcing fields
1183  !! unused fields have NULL ptrs
1184  type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and
1185  type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum
1186  !! equations, to enable the later derived
1187  !! diagnostics, like energy budgets
1188  type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations
1189  real, intent(in) :: dt !< time increment [T ~> s]
1190  type(time_type), intent(in) :: Time_end !< Time at the end of the interval
1191  type(diabatic_CS), pointer :: CS !< module control structure
1192  type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves
1193 
1194  ! local variables
1195  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1196  ea_s, & ! amount of fluid entrained from the layer above within
1197  ! one time step [H ~> m or kg m-2]
1198  eb_s, & ! amount of fluid entrained from the layer below within
1199  ! one time step [H ~> m or kg m-2]
1200  ea_t, & ! amount of fluid entrained from the layer above within
1201  ! one time step [H ~> m or kg m-2]
1202  eb_t, & ! amount of fluid entrained from the layer below within
1203  ! one time step [H ~> m or kg m-2]
1204  h_prebound, & ! initial layer thicknesses [H ~> m or kg m-2]
1205  dsv_dt, & ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
1206  dsv_ds, & ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
1207  ctke, & ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
1208  u_h, & ! zonal and meridional velocities at thickness points after
1209  v_h ! entrainment [L T-1 ~> m s-1]
1210  real, dimension(SZI_(G),SZJ_(G)) :: &
1211  SkinBuoyFlux! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
1212  real, dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag ! diagnostic array for thickness
1213  real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag ! diagnostic array for temp
1214  real, dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag ! diagnostic array for salinity
1215  real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn
1216 
1217  real :: net_ent ! The net of ea-eb at an interface.
1218 
1219  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1220  eatr, & ! The equivalent of ea and eb for tracers, which differ from ea and
1221  ebtr ! eb in that they tend to homogenize tracers in massless layers
1222  ! near the boundaries [H ~> m or kg m-2] (for Bous or non-Bouss)
1223 
1224  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
1225  Kd_int, & ! diapycnal diffusivity of interfaces [Z2 T-1 ~> m2 s-1]
1226  Kd_heat, & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1]
1227  Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1]
1228  Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to
1229  ! Kd_int [Z2 T-1 ~> m2 s-1].
1230  kd_extra_s , & ! The extra diffusivity of salinity due to double diffusion relative to
1231  ! Kd_int [Z2 T-1 ~> m2 s-1].
1232  kd_epbl, & ! boundary layer or convective diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]
1233  tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1234  tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1235  sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1236  sadv_flx ! advective diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1237 
1238  real, allocatable, dimension(:,:) :: &
1239  hf_dudt_dia_2d, hf_dvdt_dia_2d ! Depth sum of diapycnal mixing accelaration * fract. thickness [L T-2 ~> m s-2].
1240 
1241  logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below,
1242  ! where massive is defined as sufficiently thick that
1243  ! the no-flux boundary conditions have not restricted
1244  ! the entrainment - usually sqrt(Kd*dt).
1245 
1246  real :: b_denom_1 ! The first term in the denominator of b1
1247  ! [H ~> m or kg m-2]
1248  real :: h_neglect ! A thickness that is so small it is usually lost
1249  ! in roundoff and can be neglected
1250  ! [H ~> m or kg m-2]
1251  real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4]
1252  real :: add_ent ! Entrainment that needs to be added when mixing tracers
1253  ! [H ~> m or kg m-2]
1254  real :: eaval ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
1255  real :: hval ! hval is 2*h at velocity grid points [H ~> m or kg m-2]
1256  real :: h_tr ! h_tr is h at tracer points with a tiny thickness
1257  ! added to ensure positive definiteness [H ~> m or kg m-2]
1258  real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is
1259  ! coupled to the bottom within a timestep [H ~> m or kg m-2]
1260 
1261  real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].
1262  real :: b1(SZIB_(G)), d1(SZIB_(G)) ! b1, c1, and d1 are variables used by the
1263  real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver.
1264 
1265  real :: Kd_add_here ! An added diffusivity [Z2 T-1 ~> m2 s-1].
1266  real :: Idt ! The inverse time step [T-1 ~> s-1]
1267 
1268  integer :: dir_flag ! An integer encoding the directions in which to do halo updates.
1269  logical :: showCallTree ! If true, show the call tree
1270  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m
1271 
1272  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1273  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1274  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
1275  kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0
1276  ea_s(:,:,:) = 0.0; eb_s(:,:,:) = 0.0; ea_t(:,:,:) = 0.0; eb_t(:,:,:) = 0.0
1277 
1278  showcalltree = calltree_showquery()
1279  if (showcalltree) call calltree_enter("diabatic_ALE(), MOM_diabatic_driver.F90")
1280 
1281  if (.not. (cs%useALEalgorithm)) call mom_error(fatal, "MOM_diabatic_driver: "// &
1282  "The ALE algorithm must be enabled when using MOM_diabatic_driver.")
1283 
1284  ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
1285  call enable_averages(dt, time_end, cs%diag)
1286 
1287  if (cs%use_geothermal) then
1288  call cpu_clock_begin(id_clock_geothermal)
1289  call geothermal_in_place(h, tv, dt, g, gv, us, cs%geothermal_CSp, halo=cs%halo_TS_diff)
1290  call cpu_clock_end(id_clock_geothermal)
1291  if (showcalltree) call calltree_waypoint("geothermal (diabatic)")
1292  if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)
1293  endif
1294 
1295  ! Whenever thickness changes let the diag manager know, target grids
1296  ! for vertical remapping may need to be regenerated.
1297  call diag_update_remap_grids(cs%diag)
1298 
1299  ! Set_pen_shortwave estimates the optical properties of the water column.
1300  ! It will need to be modified later to include information about the
1301  ! biological properties and layer thicknesses.
1302  if (associated(cs%optics)) &
1303  call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
1304 
1305  if (cs%debug) call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)
1306 
1307  if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then
1308  if (cs%use_geothermal) then
1309  ! The presence of eatr and ebtr causes find_uv_at_h to use a tridiagonal solver,
1310  ! which changes answers at the level of roundoff because ((A*B / A) /= B).
1311  eatr(:,:,:) = 0.0 ; ebtr(:,:,:) = 0.0
1312  call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us, eatr, ebtr)
1313  else
1314  call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
1315  endif
1316  if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
1317  endif
1318 
1319  call cpu_clock_begin(id_clock_set_diffusivity)
1320  ! Sets: Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb
1321  ! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow
1322  if (cs%debug) &
1323  call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)
1324  if (cs%double_diffuse) then
1325  call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, cs%set_diff_CSp, &
1326  kd_int=kd_int, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)
1327  else
1328  call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, &
1329  cs%set_diff_CSp, kd_int=kd_int)
1330  endif
1331  call cpu_clock_end(id_clock_set_diffusivity)
1332  if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")
1333 
1334  if (cs%debug) then
1335  call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)
1336  call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)
1337  call mom_thermovar_chksum("after set_diffusivity ", tv, g)
1338  call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1339  endif
1340 
1341  ! Set diffusivities for heat and salt separately
1342 
1343  ! Add contribution from double diffusion
1344  if (cs%double_diffuse) then
1345  !$OMP parallel do default(shared)
1346  do k=1,nz+1 ; do j=js,je ; do i=is,ie
1347  kd_salt(i,j,k) = kd_int(i,j,k) + kd_extra_s(i,j,k)
1348  kd_heat(i,j,k) = kd_int(i,j,k) + kd_extra_t(i,j,k)
1349  enddo ; enddo ; enddo
1350  else
1351  !$OMP parallel do default(shared)
1352  do k=1,nz+1 ; do j=js,je ; do i=is,ie
1353  kd_salt(i,j,k) = kd_int(i,j,k)
1354  kd_heat(i,j,k) = kd_int(i,j,k)
1355  enddo ; enddo ; enddo
1356  endif
1357 
1358  if (cs%debug) then
1359  call hchksum(kd_heat, "after set_diffusivity Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1360  call hchksum(kd_salt, "after set_diffusivity Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1361  endif
1362 
1363  if (cs%useKPP) then
1364  call cpu_clock_begin(id_clock_kpp)
1365  ! total vertical viscosity in the interior is represented via visc%Kv_shear
1366  do k=1,nz+1 ; do j=js,je ; do i=is,ie
1367  visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k)
1368  enddo ; enddo ; enddo
1369 
1370  ! KPP needs the surface buoyancy flux but does not update state variables.
1371  ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
1372  ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux
1373  ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second)
1374  ! unlike other instances where the fluxes are integrated in time over a time-step.
1375  call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
1376  cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
1377 
1378  ! The KPP scheme calculates boundary layer diffusivities and non-local transport.
1379  call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &
1380  fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
1381 
1382  call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
1383  kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
1384 
1385  if (associated(hml)) then
1386  call kpp_get_bld(cs%KPP_CSp, hml(:,:), g, us)
1387  call pass_var(hml, g%domain, halo=1)
1388  ! If visc%MLD exists, copy KPP's BLD into it
1389  if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
1390  endif
1391 
1392  call cpu_clock_end(id_clock_kpp)
1393  if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")
1394  if (cs%debug) then
1395  call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)
1396  call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)
1397  call mom_thermovar_chksum("after KPP", tv, g)
1398  call hchksum(kd_heat, "after KPP Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1399  call hchksum(kd_salt, "after KPP Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1400  endif
1401 
1402  endif ! endif for KPP
1403 
1404 
1405  if (cs%useKPP) then
1406  call cpu_clock_begin(id_clock_kpp)
1407  if (cs%debug) then
1408  call hchksum(cs%KPP_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
1409  call hchksum(cs%KPP_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
1410  call hchksum(cs%KPP_NLTheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)
1411  call hchksum(cs%KPP_NLTscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
1412  endif
1413  ! Apply non-local transport of heat and salt
1414  ! Changes: tv%T, tv%S
1415  call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat, cs%KPP_temp_flux, &
1416  us%T_to_s*dt, tv%T, us%Q_to_J_kg*tv%C_p)
1417  call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, &
1418  us%T_to_s*dt, tv%S)
1419  call cpu_clock_end(id_clock_kpp)
1420  if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
1421  if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)
1422 
1423  if (cs%debug) then
1424  call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)
1425  call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)
1426  call mom_thermovar_chksum("after KPP_applyNLT ", tv, g)
1427  endif
1428  endif ! endif for KPP
1429 
1430  ! This is the "old" method for applying differential diffusion.
1431  ! Changes: tv%T, tv%S
1432  if (cs%double_diffuse .and. associated(tv%T) .and. (.not.cs%use_CVMix_ddiff)) then
1433 
1434  call cpu_clock_begin(id_clock_differential_diff)
1435  call differential_diffuse_t_s(h, tv%T, tv%S, kd_extra_t, kd_extra_s, dt, g, gv)
1436  call cpu_clock_end(id_clock_differential_diff)
1437 
1438  if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")
1439  if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g, gv, us)
1440 
1441  ! increment heat and salt diffusivity.
1442  ! CS%useKPP==.true. already has extra_T and extra_S included
1443  if (.not. cs%useKPP) then
1444  !$OMP parallel do default(shared)
1445  do k=2,nz ; do j=js,je ; do i=is,ie
1446  kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)
1447  kd_salt(i,j,k) = kd_salt(i,j,k) + kd_extra_s(i,j,k)
1448  enddo ; enddo ; enddo
1449  endif
1450 
1451  endif
1452 
1453  ! Calculate vertical mixing due to convection (computed via CVMix)
1454  if (cs%use_CVMix_conv) then
1455  ! Increment vertical diffusion and viscosity due to convection
1456  if (cs%useKPP) then
1457  call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml, kd=kd_heat, kv=visc%Kv_shear, kd_aux=kd_salt)
1458  else
1459  call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml, kd=kd_heat, kv=visc%Kv_slow, kd_aux=kd_salt)
1460  endif
1461  endif
1462 
1463  ! Save fields before boundary forcing is applied for tendency diagnostics
1464  if (cs%boundary_forcing_tendency_diag) then
1465  do k=1,nz ; do j=js,je ; do i=is,ie
1466  h_diag(i,j,k) = h(i,j,k)
1467  temp_diag(i,j,k) = tv%T(i,j,k)
1468  saln_diag(i,j,k) = tv%S(i,j,k)
1469  enddo ; enddo ; enddo
1470  endif
1471 
1472  ! Apply forcing
1473  call cpu_clock_begin(id_clock_remap)
1474 
1475  ! Changes made to following fields: h, tv%T and tv%S.
1476  do k=1,nz ; do j=js,je ; do i=is,ie
1477  h_prebound(i,j,k) = h(i,j,k)
1478  enddo ; enddo ; enddo
1479  if (cs%use_energetic_PBL) then
1480 
1481  skinbuoyflux(:,:) = 0.0
1482  call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
1483  optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, cs%evap_CFL_limit, &
1484  cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux)
1485 
1486  if (cs%debug) then
1487  call hchksum(ea_t, "after applyBoundaryFluxes ea_t", g%HI, haloshift=0, scale=gv%H_to_m)
1488  call hchksum(eb_t, "after applyBoundaryFluxes eb_t", g%HI, haloshift=0, scale=gv%H_to_m)
1489  call hchksum(ea_s, "after applyBoundaryFluxes ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
1490  call hchksum(eb_s, "after applyBoundaryFluxes eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
1491  call hchksum(ctke, "after applyBoundaryFluxes cTKE", g%HI, haloshift=0, &
1492  scale=us%RZ3_T3_to_W_m2*us%T_to_s)
1493  call hchksum(dsv_dt, "after applyBoundaryFluxes dSV_dT", g%HI, haloshift=0, scale=us%kg_m3_to_R)
1494  call hchksum(dsv_ds, "after applyBoundaryFluxes dSV_dS", g%HI, haloshift=0, scale=us%kg_m3_to_R)
1495  endif
1496 
1497  call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
1498  call energetic_pbl(h, u_h, v_h, tv, fluxes, dt, kd_epbl, g, gv, us, &
1499  cs%energetic_PBL_CSp, dsv_dt, dsv_ds, ctke, skinbuoyflux, waves=waves)
1500 
1501  if (associated(hml)) then
1502  call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hml(:,:), g, us)
1503  call pass_var(hml, g%domain, halo=1)
1504  ! If visc%MLD exists, copy ePBL's MLD into it
1505  if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
1506  elseif (associated(visc%MLD)) then
1507  call energetic_pbl_get_mld(cs%energetic_PBL_CSp, visc%MLD, g, us)
1508  call pass_var(visc%MLD, g%domain, halo=1)
1509  endif
1510 
1511  ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
1512  do k=2,nz ; do j=js,je ; do i=is,ie
1513  if (cs%ePBL_is_additive) then
1514  kd_add_here = kd_epbl(i,j,k)
1515  visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%ePBL_Prandtl*kd_epbl(i,j,k)
1516  else
1517  kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)
1518  visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), cs%ePBL_Prandtl*kd_epbl(i,j,k))
1519  endif
1520 
1521  kd_heat(i,j,k) = kd_heat(i,j,k) + kd_add_here
1522  kd_salt(i,j,k) = kd_salt(i,j,k) + kd_add_here
1523 
1524  enddo ; enddo ; enddo
1525 
1526  if (cs%debug) then
1527  call hchksum(ea_t, "after ePBL ea_t", g%HI, haloshift=0, scale=gv%H_to_m)
1528  call hchksum(eb_t, "after ePBL eb_t", g%HI, haloshift=0, scale=gv%H_to_m)
1529  call hchksum(ea_s, "after ePBL ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
1530  call hchksum(eb_s, "after ePBL eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
1531  call hchksum(kd_epbl, "after ePBL Kd_ePBL", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1532  endif
1533 
1534  else
1535  call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
1536  optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, &
1537  cs%evap_CFL_limit, cs%minimum_forcing_depth)
1538 
1539  endif ! endif for CS%use_energetic_PBL
1540 
1541  ! diagnose the tendencies due to boundary forcing
1542  ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme
1543  ! so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards
1544  if (cs%boundary_forcing_tendency_diag) then
1545  call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, g, gv, us, cs)
1546  if (cs%id_boundary_forcing_h > 0) call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h = h_diag)
1547  endif
1548  ! Boundary fluxes may have changed T, S, and h
1549  call diag_update_remap_grids(cs%diag)
1550  call cpu_clock_end(id_clock_remap)
1551  if (cs%debug) then
1552  call mom_forcing_chksum("after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)
1553  call mom_thermovar_chksum("after applyBoundaryFluxes ", tv, g)
1554  call mom_state_chksum("after applyBoundaryFluxes ", u, v, h, g, gv, us, haloshift=0)
1555  endif
1556  if (showcalltree) call calltree_waypoint("done with applyBoundaryFluxes (diabatic)")
1557  if (cs%debugConservation) call mom_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, g, gv, us)
1558 
1559  if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")
1560  if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)
1561 
1562  ! calculate change in temperature & salinity due to dia-coordinate surface diffusion
1563  if (associated(tv%T)) then
1564 
1565  if (cs%debug) then
1566  call hchksum(ea_t, "before triDiagTS ea_t ", g%HI, haloshift=0, scale=gv%H_to_m)
1567  call hchksum(eb_t, "before triDiagTS eb_t ", g%HI, haloshift=0, scale=gv%H_to_m)
1568  call hchksum(ea_s, "before triDiagTS ea_s ", g%HI, haloshift=0, scale=gv%H_to_m)
1569  call hchksum(eb_s, "before triDiagTS eb_s ", g%HI, haloshift=0, scale=gv%H_to_m)
1570  endif
1571 
1572  call cpu_clock_begin(id_clock_tridiag)
1573  ! Keep salinity from falling below a small but positive threshold.
1574  ! This constraint is needed for SIS1 ice model, which can extract
1575  ! more salt than is present in the ocean. SIS2 does not suffer
1576  ! from this limitation, in which case we can let salinity=0 and still
1577  ! have salt conserved with SIS2 ice. So for SIS2, we can run with
1578  ! BOUND_SALINITY=False in MOM.F90.
1579  if (associated(tv%S) .and. associated(tv%salt_deficit)) &
1580  call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
1581 
1582  if (cs%diabatic_diff_tendency_diag) then
1583  do k=1,nz ; do j=js,je ; do i=is,ie
1584  temp_diag(i,j,k) = tv%T(i,j,k)
1585  saln_diag(i,j,k) = tv%S(i,j,k)
1586  enddo ; enddo ; enddo
1587  endif
1588 
1589  ! set ea_t=eb_t=Kd_heat and ea_s=eb_s=Kd_salt on interfaces for use in the tridiagonal solver.
1590  do j=js,je ; do i=is,ie
1591  ea_t(i,j,1) = 0.; ea_s(i,j,1) = 0.
1592  enddo ; enddo
1593 
1594  !$OMP parallel do default(shared) private(hval)
1595  do k=2,nz ; do j=js,je ; do i=is,ie
1596  hval = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
1597  ea_t(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_heat(i,j,k)
1598  eb_t(i,j,k-1) = ea_t(i,j,k)
1599  ea_s(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_salt(i,j,k)
1600  eb_s(i,j,k-1) = ea_s(i,j,k)
1601  enddo ; enddo ; enddo
1602  do j=js,je ; do i=is,ie
1603  eb_t(i,j,nz) = 0. ; eb_s(i,j,nz) = 0.
1604  enddo ; enddo
1605  if (showcalltree) call calltree_waypoint("done setting ea_t,ea_s,eb_t,eb_s from Kd_heat" //&
1606  "and Kd_salt (diabatic)")
1607 
1608  ! Initialize halo regions of ea_t, eb_t, ea_s and eb_s to default values.
1609  !$OMP parallel do default(shared)
1610  do k=1,nz
1611  do i=is-1,ie+1
1612  ea_t(i,js-1,k) = 0.0 ; eb_t(i,js-1,k) = 0.0
1613  ea_s(i,js-1,k) = 0.0 ; eb_s(i,js-1,k) = 0.0
1614  ea_t(i,je+1,k) = 0.0 ; eb_t(i,je+1,k) = 0.0
1615  ea_s(i,je+1,k) = 0.0 ; eb_s(i,je+1,k) = 0.0
1616  enddo
1617  do j=js,je
1618  ea_t(is-1,j,k) = 0.0 ; eb_t(is-1,j,k) = 0.0
1619  ea_s(is-1,j,k) = 0.0 ; eb_s(is-1,j,k) = 0.0
1620  ea_t(ie+1,j,k) = 0.0 ; eb_t(ie+1,j,k) = 0.0
1621  ea_s(ie+1,j,k) = 0.0 ; eb_s(ie+1,j,k) = 0.0
1622  enddo
1623  enddo
1624 
1625  ! Changes T and S via the tridiagonal solver; no change to h
1626  call tracer_vertdiff(h, ea_t, eb_t, dt, tv%T, g, gv)
1627  call tracer_vertdiff(h, ea_s, eb_s, dt, tv%S, g, gv)
1628 
1629 
1630  ! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below
1631  if (cs%diabatic_diff_tendency_diag) then
1632  call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, g, gv, us, cs)
1633  endif
1634  call cpu_clock_end(id_clock_tridiag)
1635 
1636  if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")
1637 
1638  endif ! endif corresponding to if (associated(tv%T))
1639 
1640  if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)
1641 
1642  if (cs%debug) then
1643  call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)
1644  call mom_thermovar_chksum("after mixed layer ", tv, g)
1645  endif
1646 
1647  ! Whenever thickness changes let the diag manager know, as the
1648  ! target grids for vertical remapping may need to be regenerated.
1649  call diag_update_remap_grids(cs%diag)
1650 
1651  ! diagnostics
1652  idt = 1.0 / dt
1653  if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0)) then
1654  do j=js,je ; do i=is,ie
1655  tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
1656  tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
1657  enddo ; enddo
1658  !$OMP parallel do default(shared)
1659  do k=2,nz ; do j=js,je ; do i=is,ie
1660  tdif_flx(i,j,k) = (idt * 0.5*(ea_t(i,j,k) + eb_t(i,j,k-1))) * &
1661  (tv%T(i,j,k-1) - tv%T(i,j,k))
1662  tadv_flx(i,j,k) = (idt * (ea_t(i,j,k) - eb_t(i,j,k-1))) * &
1663  0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
1664  enddo ; enddo ; enddo
1665  endif
1666  if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0)) then
1667  do j=js,je ; do i=is,ie
1668  sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
1669  sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
1670  enddo ; enddo
1671  !$OMP parallel do default(shared)
1672  do k=2,nz ; do j=js,je ; do i=is,ie
1673  sdif_flx(i,j,k) = (idt * 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))) * &
1674  (tv%S(i,j,k-1) - tv%S(i,j,k))
1675  sadv_flx(i,j,k) = (idt * (ea_s(i,j,k) - eb_s(i,j,k-1))) * &
1676  0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
1677  enddo ; enddo ; enddo
1678  endif
1679 
1680  ! mixing of passive tracers from massless boundary layers to interior
1681  call cpu_clock_begin(id_clock_tracers)
1682 
1683  if (cs%mix_boundary_tracers) then
1684  tr_ea_bbl = gv%Z_to_H * sqrt(dt*cs%Kd_BBL_tr)
1685  !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
1686  do j=js,je
1687  do i=is,ie
1688  ebtr(i,j,nz) = eb_s(i,j,nz)
1689  htot(i) = 0.0
1690  in_boundary(i) = (g%mask2dT(i,j) > 0.0)
1691  enddo
1692  do k=nz,2,-1 ; do i=is,ie
1693  if (in_boundary(i)) then
1694  htot(i) = htot(i) + h(i,j,k)
1695  ! If diapycnal mixing has been suppressed because this is a massless
1696  ! layer near the bottom, add some mixing of tracers between these
1697  ! layers. This flux is based on the harmonic mean of the two
1698  ! thicknesses, as this corresponds pretty closely (to within
1699  ! differences in the density jumps between layers) with what is done
1700  ! in the calculation of the fluxes in the first place. Kd_min_tr
1701  ! should be much less than the values that have been set in Kd_int,
1702  ! perhaps a molecular diffusivity.
1703  add_ent = ((dt * cs%Kd_min_tr) * gv%Z_to_H**2) * &
1704  ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
1705  (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
1706  0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))
1707  if (htot(i) < tr_ea_bbl) then
1708  add_ent = max(0.0, add_ent, &
1709  (tr_ea_bbl - htot(i)) - min(ea_s(i,j,k),eb_s(i,j,k-1)))
1710  elseif (add_ent < 0.0) then
1711  add_ent = 0.0 ; in_boundary(i) = .false.
1712  endif
1713 
1714  ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
1715  eatr(i,j,k) = ea_s(i,j,k) + add_ent
1716  else
1717  ebtr(i,j,k-1) = eb_s(i,j,k-1) ; eatr(i,j,k) = ea_s(i,j,k)
1718  endif
1719 
1720  if (cs%double_diffuse) then ; if (kd_extra_s(i,j,k) > 0.0) then
1721  add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
1722  (0.5 * (h(i,j,k-1) + h(i,j,k)) + &
1723  h_neglect)
1724  ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
1725  eatr(i,j,k) = eatr(i,j,k) + add_ent
1726  endif ; endif
1727  enddo ; enddo
1728  do i=is,ie ; eatr(i,j,1) = ea_s(i,j,1) ; enddo
1729 
1730  enddo
1731 
1732  ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
1733  ! so use h_prebound as the old value.
1734  !### I think that in the following, ea_s and eb_s should be eatr and ebtr.
1735  call call_tracer_column_fns(h_prebound, h, ea_s, eb_s, fluxes, hml, dt, g, gv, us, tv, &
1736  cs%optics, cs%tracer_flow_CSp, cs%debug, &
1737  evap_cfl_limit = cs%evap_CFL_limit, &
1738  minimum_forcing_depth=cs%minimum_forcing_depth)
1739 
1740  elseif (cs%double_diffuse) then ! extra diffusivity for passive tracers
1741 
1742  do j=js,je ; do i=is,ie
1743  ebtr(i,j,nz) = eb_s(i,j,nz) ; eatr(i,j,1) = ea_s(i,j,1)
1744  enddo ; enddo
1745  !$OMP parallel do default(shared) private(add_ent)
1746  do k=nz,2,-1 ; do j=js,je ; do i=is,ie
1747  if (kd_extra_s(i,j,k) > 0.0) then
1748  add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
1749  (0.5 * (h(i,j,k-1) + h(i,j,k)) + &
1750  h_neglect)
1751  else
1752  add_ent = 0.0
1753  endif
1754  ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
1755  eatr(i,j,k) = ea_s(i,j,k) + add_ent
1756  enddo ; enddo ; enddo
1757 
1758  ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
1759  call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
1760  cs%optics, cs%tracer_flow_CSp, cs%debug,&
1761  evap_cfl_limit = cs%evap_CFL_limit, &
1762  minimum_forcing_depth=cs%minimum_forcing_depth)
1763  else
1764  ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
1765  !### eatr and ebtr may not be initialized or may be 0, depending on CS%geothermal.
1766  call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
1767  cs%optics, cs%tracer_flow_CSp, cs%debug, &
1768  evap_cfl_limit = cs%evap_CFL_limit, &
1769  minimum_forcing_depth=cs%minimum_forcing_depth)
1770  endif ! (CS%mix_boundary_tracers)
1771 
1772  call cpu_clock_end(id_clock_tracers)
1773 
1774  ! Apply ALE sponge
1775  if (cs%use_sponge) then
1776  call cpu_clock_begin(id_clock_sponge)
1777  if (associated(cs%ALE_sponge_CSp)) then
1778  call apply_ale_sponge(h, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)
1779  endif
1780 
1781  call cpu_clock_end(id_clock_sponge)
1782  if (cs%debug) then
1783  call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)
1784  call mom_thermovar_chksum("apply_sponge ", tv, g)
1785  endif
1786  endif ! CS%use_sponge
1787 
1788  call cpu_clock_begin(id_clock_pass)
1789  if (g%symmetric) then ; dir_flag = to_all+omit_corners
1790  else ; dir_flag = to_west+to_south+omit_corners ; endif
1791  call create_group_pass(cs%pass_hold_eb_ea, eb_t, g%Domain, dir_flag, halo=1)
1792  call create_group_pass(cs%pass_hold_eb_ea, eb_s, g%Domain, dir_flag, halo=1)
1793  call create_group_pass(cs%pass_hold_eb_ea, ea_t, g%Domain, dir_flag, halo=1)
1794  call create_group_pass(cs%pass_hold_eb_ea, ea_s, g%Domain, dir_flag, halo=1)
1795  call do_group_pass(cs%pass_hold_eb_ea, g%Domain)
1796  ! visc%Kv_slow is not in the group pass because it has larger vertical extent.
1797  if (associated(visc%Kv_slow)) &
1798  call pass_var(visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)
1799  call cpu_clock_end(id_clock_pass)
1800 
1801  call disable_averaging(cs%diag)
1802 
1803  ! Diagnose the diapycnal diffusivities and other related quantities.
1804  call enable_averages(dt, time_end, cs%diag)
1805 
1806  if (cs%id_Kd_interface > 0) call post_data(cs%id_Kd_interface, kd_int, cs%diag)
1807  if (cs%id_Kd_heat > 0) call post_data(cs%id_Kd_heat, kd_heat, cs%diag)
1808  if (cs%id_Kd_salt > 0) call post_data(cs%id_Kd_salt, kd_salt, cs%diag)
1809  if (cs%id_Kd_ePBL > 0) call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)
1810 
1811  if (cs%id_ea_t > 0) call post_data(cs%id_ea_t, ea_t, cs%diag)
1812  if (cs%id_eb_t > 0) call post_data(cs%id_eb_t, eb_t, cs%diag)
1813  if (cs%id_ea_s > 0) call post_data(cs%id_ea_s, ea_s, cs%diag)
1814  if (cs%id_eb_s > 0) call post_data(cs%id_eb_s, eb_s, cs%diag)
1815 
1816  if (cs%id_dudt_dia > 0) call post_data(cs%id_dudt_dia, adp%du_dt_dia, cs%diag)
1817  if (cs%id_dvdt_dia > 0) call post_data(cs%id_dvdt_dia, adp%dv_dt_dia, cs%diag)
1818 
1819  if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)
1820  if (cs%id_Tadv > 0) call post_data(cs%id_Tadv, tadv_flx, cs%diag)
1821  if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)
1822  if (cs%id_Sadv > 0) call post_data(cs%id_Sadv, sadv_flx, cs%diag)
1823 
1824  !! Diagnostics for terms multiplied by fractional thicknesses
1825  if (cs%id_hf_dudt_dia_2d > 0) then
1826  allocate(hf_dudt_dia_2d(g%IsdB:g%IedB,g%jsd:g%jed))
1827  hf_dudt_dia_2d(:,:) = 0.0
1828  do k=1,nz ; do j=js,je ; do i=isq,ieq
1829  hf_dudt_dia_2d(i,j) = hf_dudt_dia_2d(i,j) + adp%du_dt_dia(i,j,k) * adp%diag_hfrac_u(i,j,k)
1830  enddo ; enddo ; enddo
1831  call post_data(cs%id_hf_dudt_dia_2d, hf_dudt_dia_2d, cs%diag)
1832  deallocate(hf_dudt_dia_2d)
1833  endif
1834 
1835  if (cs%id_hf_dvdt_dia_2d > 0) then
1836  allocate(hf_dvdt_dia_2d(g%isd:g%ied,g%JsdB:g%JedB))
1837  hf_dvdt_dia_2d(:,:) = 0.0
1838  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1839  hf_dvdt_dia_2d(i,j) = hf_dvdt_dia_2d(i,j) + adp%dv_dt_dia(i,j,k) * adp%diag_hfrac_v(i,j,k)
1840  enddo ; enddo ; enddo
1841  call post_data(cs%id_hf_dvdt_dia_2d, hf_dvdt_dia_2d, cs%diag)
1842  deallocate(hf_dvdt_dia_2d)
1843  endif
1844 
1845  call disable_averaging(cs%diag)
1846 
1847  if (showcalltree) call calltree_leave("diabatic_ALE()")
1848 

◆ diabatic_ale_legacy()

subroutine mom_diabatic_driver::diabatic_ale_legacy ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(inout)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h,
type(thermo_var_ptrs), intent(inout)  tv,
real, dimension(:,:), pointer  Hml,
type(forcing), intent(inout)  fluxes,
type(vertvisc_type), intent(inout)  visc,
type(accel_diag_ptrs), intent(inout)  ADp,
type(cont_diag_ptrs), intent(inout)  CDp,
real, intent(in)  dt,
type(time_type), intent(in)  Time_end,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diabatic_cs), pointer  CS,
type(wave_parameters_cs), optional, pointer  Waves 
)
private

Applies diabatic forcing and diapycnal mixing of temperature, salinity and other tracers for use with an ALE algorithm. This version uses an older set of algorithms compared with diabatic_ALE.

Parameters
[in,out]gocean grid structure
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]uzonal velocity [L T-1 ~> m s-1]
[in,out]vmeridional velocity [L T-1 ~> m s-1]
[in,out]hthickness [H ~> m or kg m-2]
[in,out]tvpoints to thermodynamic fields unused have NULL ptrs
hmlActive mixed layer depth [Z ~> m]
[in,out]fluxespoints to forcing fields unused fields have NULL ptrs
[in,out]viscvertical viscosities, BBL properies, and
[in,out]adprelated points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]cdppoints to terms in continuity equations
[in]dttime increment [T ~> s]
[in]time_endTime at the end of the interval
csmodule control structure
wavesSurface gravity waves

Definition at line 453 of file MOM_diabatic_driver.F90.

453  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
454  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
455  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
456  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
457  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
458  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2]
459  type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
460  !! unused have NULL ptrs
461  real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m]
462  type(forcing), intent(inout) :: fluxes !< points to forcing fields
463  !! unused fields have NULL ptrs
464  type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and
465  type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum
466  !! equations, to enable the later derived
467  !! diagnostics, like energy budgets
468  type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations
469  real, intent(in) :: dt !< time increment [T ~> s]
470  type(time_type), intent(in) :: Time_end !< Time at the end of the interval
471  type(diabatic_CS), pointer :: CS !< module control structure
472  type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves
473 
474  ! local variables
475  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
476  ea_s, & ! amount of fluid entrained from the layer above within
477  ! one time step [H ~> m or kg m-2]
478  eb_s, & ! amount of fluid entrained from the layer below within
479  ! one time step [H ~> m or kg m-2]
480  ea_t, & ! amount of fluid entrained from the layer above within
481  ! one time step [H ~> m or kg m-2]
482  eb_t, & ! amount of fluid entrained from the layer below within
483  ! one time step [H ~> m or kg m-2]
484  h_orig, & ! initial layer thicknesses [H ~> m or kg m-2]
485  h_prebound, & ! initial layer thicknesses [H ~> m or kg m-2]
486  hold, & ! layer thickness before diapycnal entrainment, and later
487  ! the initial layer thicknesses (if a mixed layer is used),
488  ! [H ~> m or kg m-2]
489  dsv_dt, & ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
490  dsv_ds, & ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
491  ctke, & ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
492  u_h, & ! zonal and meridional velocities at thickness points after
493  v_h ! entrainment [L T-1 ~> m s-1]
494  real, dimension(SZI_(G),SZJ_(G)) :: &
495  SkinBuoyFlux! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
496  real, dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag ! diagnostic array for thickness
497  real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag ! diagnostic array for temp
498  real, dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag ! diagnostic array for salinity
499  real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn
500 
501  real :: net_ent ! The net of ea-eb at an interface.
502 
503  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
504  eatr, & ! The equivalent of ea and eb for tracers, which differ from ea and
505  ebtr ! eb in that they tend to homogenize tracers in massless layers
506  ! near the boundaries [H ~> m or kg m-2] (for Bous or non-Bouss)
507 
508  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
509  Kd_int, & ! diapycnal diffusivity of interfaces [Z2 T-1 ~> m2 s-1]
510  Kd_heat, & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1]
511  Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1]
512  Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to
513  ! Kd_int [Z2 T-1 ~> m2 s-1].
514  kd_extra_s , & ! The extra diffusivity of salinity due to double diffusion relative to
515  ! Kd_int [Z2 T-1 ~> m2 s-1].
516  kd_epbl, & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]
517  tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
518  tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
519  sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
520  sadv_flx ! advective diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
521 
522  real, allocatable, dimension(:,:) :: &
523  hf_dudt_dia_2d, hf_dvdt_dia_2d ! Depth sum of diapycnal mixing accelaration * fract. thickness [L T-2 ~> m s-2].
524 
525  logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below,
526  ! where massive is defined as sufficiently thick that
527  ! the no-flux boundary conditions have not restricted
528  ! the entrainment - usually sqrt(Kd*dt).
529 
530  real :: b_denom_1 ! The first term in the denominator of b1
531  ! [H ~> m or kg m-2]
532  real :: h_neglect ! A thickness that is so small it is usually lost
533  ! in roundoff and can be neglected
534  ! [H ~> m or kg m-2]
535  real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4]
536  real :: add_ent ! Entrainment that needs to be added when mixing tracers
537  ! [H ~> m or kg m-2]
538  real :: eaval ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
539  real :: hval ! hval is 2*h at velocity grid points [H ~> m or kg m-2]
540  real :: h_tr ! h_tr is h at tracer points with a tiny thickness
541  ! added to ensure positive definiteness [H ~> m or kg m-2]
542  real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is
543  ! coupled to the bottom within a timestep [H ~> m or kg m-2]
544 
545  real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].
546  real :: b1(SZIB_(G)), d1(SZIB_(G)) ! b1, c1, and d1 are variables used by the
547  real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver.
548 
549  real :: Ent_int ! The diffusive entrainment rate at an interface [H ~> m or kg m-2]
550  real :: Idt ! The inverse time step [T-1 ~> s-1]
551 
552  integer :: dir_flag ! An integer encoding the directions in which to do halo updates.
553  logical :: showCallTree ! If true, show the call tree
554  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m
555 
556  integer :: ig, jg ! global indices for testing testing itide point source (BDM)
557  real :: Kd_add_here ! An added diffusivity [Z2 T-1 ~> m2 s-1].
558 
559  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
560  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
561  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
562  kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0
563 
564  showcalltree = calltree_showquery()
565  if (showcalltree) call calltree_enter("diabatic_ALE_legacy(), MOM_diabatic_driver.F90")
566 
567  ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
568  call enable_averages(dt, time_end, cs%diag)
569 
570  if (cs%use_geothermal) then
571  call cpu_clock_begin(id_clock_geothermal)
572  call geothermal_in_place(h, tv, dt, g, gv, us, cs%geothermal_CSp, halo=cs%halo_TS_diff)
573  call cpu_clock_end(id_clock_geothermal)
574  if (showcalltree) call calltree_waypoint("geothermal (diabatic)")
575  if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)
576  endif
577 
578  ! Whenever thickness changes let the diag manager know, target grids
579  ! for vertical remapping may need to be regenerated.
580  call diag_update_remap_grids(cs%diag)
581 
582  ! Set_pen_shortwave estimates the optical properties of the water column.
583  ! It will need to be modified later to include information about the
584  ! biological properties and layer thicknesses.
585  if (associated(cs%optics)) &
586  call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
587 
588  if (cs%debug) call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)
589 
590  if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then
591  if (cs%use_geothermal) then
592  ! The presence of eatr and ebtr causes find_uv_at_h to use a tridiagonal solver,
593  ! which changes answers at the level of roundoff because ((A*B / A) /= B).
594  eatr(:,:,:) = 0.0 ; ebtr(:,:,:) = 0.0
595  call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us, eatr, ebtr)
596  else
597  call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
598  endif
599  if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
600  endif
601 
602  call cpu_clock_begin(id_clock_set_diffusivity)
603  ! Sets: Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb
604  ! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow
605  if (cs%debug) &
606  call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)
607  if (cs%double_diffuse) then
608  call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, cs%set_diff_CSp, &
609  kd_int=kd_int, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)
610  else
611  call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, &
612  cs%set_diff_CSp, kd_int=kd_int)
613  endif
614  call cpu_clock_end(id_clock_set_diffusivity)
615  if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")
616 
617  ! Set diffusivities for heat and salt separately
618 
619  if (cs%useKPP) then
620  ! Add contribution from double diffusion
621  if (cs%double_diffuse) then
622  !$OMP parallel do default(shared)
623  do k=1,nz+1 ; do j=js,je ; do i=is,ie
624  kd_salt(i,j,k) = kd_int(i,j,k) + kd_extra_s(i,j,k)
625  kd_heat(i,j,k) = kd_int(i,j,k) + kd_extra_t(i,j,k)
626  enddo ; enddo ; enddo
627  else
628  !$OMP parallel do default(shared)
629  do k=1,nz+1 ; do j=js,je ; do i=is,ie
630  kd_salt(i,j,k) = kd_int(i,j,k)
631  kd_heat(i,j,k) = kd_int(i,j,k)
632  enddo ; enddo ; enddo
633  endif
634  endif
635 
636  if (cs%debug) then
637  call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)
638  call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)
639  call mom_thermovar_chksum("after set_diffusivity ", tv, g)
640  call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
641  call hchksum(kd_heat, "after set_diffusivity Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
642  call hchksum(kd_salt, "after set_diffusivity Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
643  endif
644 
645  if (cs%useKPP) then
646  call cpu_clock_begin(id_clock_kpp)
647  ! total vertical viscosity in the interior is represented via visc%Kv_shear
648 
649  ! KPP needs the surface buoyancy flux but does not update state variables.
650  ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
651  ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux
652  ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second)
653  ! unlike other instances where the fluxes are integrated in time over a time-step.
654  call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
655  cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
656  ! The KPP scheme calculates boundary layer diffusivities and non-local transport.
657 
658  call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &
659  fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
660 
661  call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
662  kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
663 
664  if (associated(hml)) then
665  call kpp_get_bld(cs%KPP_CSp, hml(:,:), g, us)
666  call pass_var(hml, g%domain, halo=1)
667  ! If visc%MLD exists, copy KPP's BLD into it
668  if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
669  endif
670 
671  if (.not.cs%KPPisPassive) then
672  !$OMP parallel do default(shared)
673  do k=1,nz+1 ; do j=js,je ; do i=is,ie
674  kd_int(i,j,k) = min( kd_salt(i,j,k), kd_heat(i,j,k) )
675  enddo ; enddo ; enddo
676  if (cs%double_diffuse) then
677  !$OMP parallel do default(shared)
678  do k=1,nz+1 ; do j=js,je ; do i=is,ie
679  kd_extra_s(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))
680  kd_extra_t(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))
681  enddo ; enddo ; enddo
682  endif
683  endif ! not passive
684 
685  call cpu_clock_end(id_clock_kpp)
686  if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")
687  if (cs%debug) then
688  call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)
689  call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)
690  call mom_thermovar_chksum("after KPP", tv, g)
691  call hchksum(kd_heat, "after KPP Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
692  call hchksum(kd_salt, "after KPP Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
693  endif
694 
695  endif ! endif for KPP
696 
697 
698  if (cs%useKPP) then
699  call cpu_clock_begin(id_clock_kpp)
700  if (cs%debug) then
701  call hchksum(cs%KPP_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
702  call hchksum(cs%KPP_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
703  call hchksum(cs%KPP_NLTheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)
704  call hchksum(cs%KPP_NLTscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
705  endif
706  ! Apply non-local transport of heat and salt
707  ! Changes: tv%T, tv%S
708  call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat, cs%KPP_temp_flux, &
709  us%T_to_s*dt, tv%T, us%Q_to_J_kg*tv%C_p)
710  call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, &
711  us%T_to_s*dt, tv%S)
712  call cpu_clock_end(id_clock_kpp)
713  if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
714  if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)
715 
716  if (cs%debug) then
717  call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)
718  call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)
719  call mom_thermovar_chksum("after KPP_applyNLT ", tv, g)
720  endif
721  endif ! endif for KPP
722 
723  ! This is the "old" method for applying differential diffusion.
724  ! Changes: tv%T, tv%S
725  if (cs%double_diffuse .and. associated(tv%T)) then
726 
727  call cpu_clock_begin(id_clock_differential_diff)
728  call differential_diffuse_t_s(h, tv%T, tv%S, kd_extra_t, kd_extra_s, dt, g, gv)
729  call cpu_clock_end(id_clock_differential_diff)
730 
731  if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")
732  if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g, gv, us)
733 
734  ! increment heat and salt diffusivity.
735  ! CS%useKPP==.true. already has extra_T and extra_S included
736  if (.not. cs%useKPP) then
737  !$OMP parallel do default(shared)
738  do k=2,nz ; do j=js,je ; do i=is,ie
739  kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)
740  kd_salt(i,j,k) = kd_salt(i,j,k) + kd_extra_s(i,j,k)
741  enddo ; enddo ; enddo
742  endif
743 
744  endif
745 
746  ! Calculate vertical mixing due to convection (computed via CVMix)
747  if (cs%use_CVMix_conv) then
748  ! Increment vertical diffusion and viscosity due to convection
749  call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml, kd=kd_int, kv=visc%Kv_slow)
750  endif
751 
752  ! This block sets ea, eb from h and Kd_int.
753  do j=js,je ; do i=is,ie
754  ea_s(i,j,1) = 0.0
755  enddo ; enddo
756  !$OMP parallel do default(shared) private(hval)
757  do k=2,nz ; do j=js,je ; do i=is,ie
758  hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
759  ea_s(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_int(i,j,k)
760  eb_s(i,j,k-1) = ea_s(i,j,k)
761  ea_t(i,j,k-1) = ea_s(i,j,k-1) ; eb_t(i,j,k-1) = eb_s(i,j,k-1)
762  enddo ; enddo ; enddo
763  do j=js,je ; do i=is,ie
764  eb_s(i,j,nz) = 0.0
765  ea_t(i,j,nz) = ea_s(i,j,nz) ; eb_t(i,j,nz) = eb_s(i,j,nz)
766  enddo ; enddo
767  if (showcalltree) call calltree_waypoint("done setting ea,eb from Kd_int (diabatic)")
768 
769  if (cs%debug) then
770  call mom_forcing_chksum("after calc_entrain ", fluxes, g, us, haloshift=0)
771  call mom_thermovar_chksum("after calc_entrain ", tv, g)
772  call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, us, haloshift=0)
773  call hchksum(ea_s, "after calc_entrain ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
774  call hchksum(eb_s, "after calc_entrain eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
775  endif
776 
777  ! Save fields before boundary forcing is applied for tendency diagnostics
778  if (cs%boundary_forcing_tendency_diag) then
779  do k=1,nz ; do j=js,je ; do i=is,ie
780  h_diag(i,j,k) = h(i,j,k)
781  temp_diag(i,j,k) = tv%T(i,j,k)
782  saln_diag(i,j,k) = tv%S(i,j,k)
783  enddo ; enddo ; enddo
784  endif
785 
786  ! Apply forcing
787  call cpu_clock_begin(id_clock_remap)
788 
789  ! Changes made to following fields: h, tv%T and tv%S.
790  do k=1,nz ; do j=js,je ; do i=is,ie
791  h_prebound(i,j,k) = h(i,j,k)
792  enddo ; enddo ; enddo
793  if (cs%use_energetic_PBL) then
794 
795  skinbuoyflux(:,:) = 0.0
796  call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
797  optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, cs%evap_CFL_limit, &
798  cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux)
799 
800  if (cs%debug) then
801  call hchksum(ea_t, "after applyBoundaryFluxes ea_t", g%HI, haloshift=0, scale=gv%H_to_m)
802  call hchksum(eb_t, "after applyBoundaryFluxes eb_t", g%HI, haloshift=0, scale=gv%H_to_m)
803  call hchksum(ea_s, "after applyBoundaryFluxes ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
804  call hchksum(eb_s, "after applyBoundaryFluxes eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
805  call hchksum(ctke, "after applyBoundaryFluxes cTKE", g%HI, haloshift=0, &
806  scale=us%RZ3_T3_to_W_m2*us%T_to_s)
807  call hchksum(dsv_dt, "after applyBoundaryFluxes dSV_dT", g%HI, haloshift=0, scale=us%kg_m3_to_R)
808  call hchksum(dsv_ds, "after applyBoundaryFluxes dSV_dS", g%HI, haloshift=0, scale=us%kg_m3_to_R)
809  endif
810 
811  call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
812  call energetic_pbl(h, u_h, v_h, tv, fluxes, dt, kd_epbl, g, gv, us, &
813  cs%energetic_PBL_CSp, dsv_dt, dsv_ds, ctke, skinbuoyflux, waves=waves)
814 
815  if (associated(hml)) then
816  call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hml(:,:), g, us)
817  call pass_var(hml, g%domain, halo=1)
818  ! If visc%MLD exists, copy ePBL's MLD into it
819  if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
820  elseif (associated(visc%MLD)) then
821  call energetic_pbl_get_mld(cs%energetic_PBL_CSp, visc%MLD, g, us)
822  call pass_var(visc%MLD, g%domain, halo=1)
823  endif
824 
825  ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
826  do k=2,nz ; do j=js,je ; do i=is,ie
827  if (cs%ePBL_is_additive) then
828  kd_add_here = kd_epbl(i,j,k)
829  visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%ePBL_Prandtl*kd_epbl(i,j,k)
830  else
831  kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)
832  visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), cs%ePBL_Prandtl*kd_epbl(i,j,k))
833  endif
834 
835  ent_int = kd_add_here * (gv%Z_to_H**2 * dt) / &
836  (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect)
837  eb_s(i,j,k-1) = eb_s(i,j,k-1) + ent_int
838  ea_s(i,j,k) = ea_s(i,j,k) + ent_int
839  kd_int(i,j,k) = kd_int(i,j,k) + kd_add_here
840 
841  ! for diagnostics
842  kd_heat(i,j,k) = kd_heat(i,j,k) + kd_int(i,j,k)
843  kd_salt(i,j,k) = kd_salt(i,j,k) + kd_int(i,j,k)
844 
845  enddo ; enddo ; enddo
846 
847  if (cs%debug) then
848  call hchksum(ea_t, "after ePBL ea_t", g%HI, haloshift=0, scale=gv%H_to_m)
849  call hchksum(eb_t, "after ePBL eb_t", g%HI, haloshift=0, scale=gv%H_to_m)
850  call hchksum(ea_s, "after ePBL ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
851  call hchksum(eb_s, "after ePBL eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
852  call hchksum(kd_epbl, "after ePBL Kd_ePBL", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
853  endif
854 
855  else
856  call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
857  optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, &
858  cs%evap_CFL_limit, cs%minimum_forcing_depth)
859 
860  endif ! endif for CS%use_energetic_PBL
861 
862  ! diagnose the tendencies due to boundary forcing
863  ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme
864  ! so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards
865  if (cs%boundary_forcing_tendency_diag) then
866  call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, g, gv, us, cs)
867  if (cs%id_boundary_forcing_h > 0) call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h = h_diag)
868  endif
869  ! Boundary fluxes may have changed T, S, and h
870  call diag_update_remap_grids(cs%diag)
871  call cpu_clock_end(id_clock_remap)
872  if (cs%debug) then
873  call mom_forcing_chksum("after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)
874  call mom_thermovar_chksum("after applyBoundaryFluxes ", tv, g)
875  call mom_state_chksum("after applyBoundaryFluxes ", u, v, h, g, gv, us, haloshift=0)
876  endif
877  if (showcalltree) call calltree_waypoint("done with applyBoundaryFluxes (diabatic)")
878  if (cs%debugConservation) call mom_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, g, gv, us)
879 
880  ! Update h according to divergence of the difference between
881  ! ea and eb. We keep a record of the original h in hold.
882  ! In the following, the checks for negative values are to guard against
883  ! instances where entrainment drives a layer to negative thickness.
884  !### This code may be unnecessary, but the negative-thickness checks do appear to change
885  ! answers slightly in some cases.
886  !$OMP parallel do default(shared)
887  do j=js,je
888  do i=is,ie
889  hold(i,j,1) = h(i,j,1)
890  ! Does nothing with ALE: h(i,j,1) = h(i,j,1) + (eb_s(i,j,1) - ea_s(i,j,2))
891  hold(i,j,nz) = h(i,j,nz)
892  ! Does nothing with ALE: h(i,j,nz) = h(i,j,nz) + (ea_s(i,j,nz) - eb_s(i,j,nz-1))
893  if (h(i,j,1) <= 0.0) h(i,j,1) = gv%Angstrom_H
894  if (h(i,j,nz) <= 0.0) h(i,j,nz) = gv%Angstrom_H
895  enddo
896  do k=2,nz-1 ; do i=is,ie
897  hold(i,j,k) = h(i,j,k)
898  ! Does nothing with ALE: h(i,j,k) = h(i,j,k) + ((ea_s(i,j,k) - eb_s(i,j,k-1)) + &
899  ! (eb_s(i,j,k) - ea_s(i,j,k+1)))
900  if (h(i,j,k) <= 0.0) h(i,j,k) = gv%Angstrom_H
901  enddo ; enddo
902  enddo
903  ! Checks for negative thickness may have changed layer thicknesses
904  call diag_update_remap_grids(cs%diag)
905 
906  if (cs%debug) then
907  call mom_state_chksum("after negative check ", u, v, h, g, gv, us, haloshift=0)
908  call mom_forcing_chksum("after negative check ", fluxes, g, us, haloshift=0)
909  call mom_thermovar_chksum("after negative check ", tv, g)
910  endif
911  if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")
912  if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)
913 
914  ! calculate change in temperature & salinity due to dia-coordinate surface diffusion
915  if (associated(tv%T)) then
916 
917  if (cs%debug) then
918  call hchksum(ea_t, "before triDiagTS ea_t ", g%HI, haloshift=0, scale=gv%H_to_m)
919  call hchksum(eb_t, "before triDiagTS eb_t ", g%HI, haloshift=0, scale=gv%H_to_m)
920  call hchksum(ea_s, "before triDiagTS ea_s ", g%HI, haloshift=0, scale=gv%H_to_m)
921  call hchksum(eb_s, "before triDiagTS eb_s ", g%HI, haloshift=0, scale=gv%H_to_m)
922  endif
923 
924  call cpu_clock_begin(id_clock_tridiag)
925  ! Keep salinity from falling below a small but positive threshold.
926  ! This constraint is needed for SIS1 ice model, which can extract
927  ! more salt than is present in the ocean. SIS2 does not suffer
928  ! from this limitation, in which case we can let salinity=0 and still
929  ! have salt conserved with SIS2 ice. So for SIS2, we can run with
930  ! BOUND_SALINITY=False in MOM.F90.
931  if (associated(tv%S) .and. associated(tv%salt_deficit)) &
932  call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
933 
934  if (cs%diabatic_diff_tendency_diag) then
935  do k=1,nz ; do j=js,je ; do i=is,ie
936  temp_diag(i,j,k) = tv%T(i,j,k)
937  saln_diag(i,j,k) = tv%S(i,j,k)
938  enddo ; enddo ; enddo
939  endif
940 
941  ! Changes T and S via the tridiagonal solver; no change to h
942  do k=1,nz ; do j=js,je ; do i=is,ie
943  ea_t(i,j,k) = ea_s(i,j,k) ; eb_t(i,j,k) = eb_s(i,j,k)
944  enddo ; enddo ; enddo
945  if (cs%tracer_tridiag) then
946  call tracer_vertdiff(hold, ea_t, eb_t, dt, tv%T, g, gv)
947  call tracer_vertdiff(hold, ea_s, eb_s, dt, tv%S, g, gv)
948  else
949  call tridiagts(g, gv, is, ie, js, je, hold, ea_s, eb_s, tv%T, tv%S)
950  endif
951 
952  ! diagnose temperature, salinity, heat, and salt tendencies
953  ! Note: hold here refers to the thicknesses from before the dual-entraintment when using
954  ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will (not?) have changed
955  ! In either case, tendencies should be posted on hold
956  if (cs%diabatic_diff_tendency_diag) then
957  call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, us, cs)
958  if (cs%id_diabatic_diff_h > 0) call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h = hold)
959  endif
960 
961  call cpu_clock_end(id_clock_tridiag)
962 
963  if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")
964 
965  endif ! endif corresponding to if (associated(tv%T))
966 
967  if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)
968 
969  if (cs%debug) then
970  call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)
971  call mom_thermovar_chksum("after mixed layer ", tv, g)
972  endif
973 
974  ! Whenever thickness changes let the diag manager know, as the
975  ! target grids for vertical remapping may need to be regenerated.
976  if (associated(adp%du_dt_dia) .or. associated(adp%dv_dt_dia)) &
977  ! Remapped d[uv]dt_dia require east/north halo updates of h
978  call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
979  call diag_update_remap_grids(cs%diag)
980 
981  ! diagnostics
982  idt = 1.0 / dt
983  if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0)) then
984  do j=js,je ; do i=is,ie
985  tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
986  tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
987  enddo ; enddo
988  !$OMP parallel do default(shared)
989  do k=2,nz ; do j=js,je ; do i=is,ie
990  tdif_flx(i,j,k) = (idt * 0.5*(ea_t(i,j,k) + eb_t(i,j,k-1))) * &
991  (tv%T(i,j,k-1) - tv%T(i,j,k))
992  tadv_flx(i,j,k) = (idt * (ea_t(i,j,k) - eb_t(i,j,k-1))) * &
993  0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
994  enddo ; enddo ; enddo
995  endif
996  if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0)) then
997  do j=js,je ; do i=is,ie
998  sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
999  sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
1000  enddo ; enddo
1001  !$OMP parallel do default(shared)
1002  do k=2,nz ; do j=js,je ; do i=is,ie
1003  sdif_flx(i,j,k) = (idt * 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))) * &
1004  (tv%S(i,j,k-1) - tv%S(i,j,k))
1005  sadv_flx(i,j,k) = (idt * (ea_s(i,j,k) - eb_s(i,j,k-1))) * &
1006  0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
1007  enddo ; enddo ; enddo
1008  endif
1009 
1010  ! mixing of passive tracers from massless boundary layers to interior
1011  call cpu_clock_begin(id_clock_tracers)
1012 
1013  if (cs%mix_boundary_tracers) then
1014  tr_ea_bbl = gv%Z_to_H * sqrt(dt*cs%Kd_BBL_tr)
1015  !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
1016  do j=js,je
1017  do i=is,ie
1018  ebtr(i,j,nz) = eb_s(i,j,nz)
1019  htot(i) = 0.0
1020  in_boundary(i) = (g%mask2dT(i,j) > 0.0)
1021  enddo
1022  do k=nz,2,-1 ; do i=is,ie
1023  if (in_boundary(i)) then
1024  htot(i) = htot(i) + h(i,j,k)
1025  ! If diapycnal mixing has been suppressed because this is a massless
1026  ! layer near the bottom, add some mixing of tracers between these
1027  ! layers. This flux is based on the harmonic mean of the two
1028  ! thicknesses, as this corresponds pretty closely (to within
1029  ! differences in the density jumps between layers) with what is done
1030  ! in the calculation of the fluxes in the first place. Kd_min_tr
1031  ! should be much less than the values that have been set in Kd_int,
1032  ! perhaps a molecular diffusivity.
1033  add_ent = ((dt * cs%Kd_min_tr) * gv%Z_to_H**2) * &
1034  ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
1035  (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
1036  0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))
1037  if (htot(i) < tr_ea_bbl) then
1038  add_ent = max(0.0, add_ent, &
1039  (tr_ea_bbl - htot(i)) - min(ea_s(i,j,k),eb_s(i,j,k-1)))
1040  elseif (add_ent < 0.0) then
1041  add_ent = 0.0 ; in_boundary(i) = .false.
1042  endif
1043 
1044  ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
1045  eatr(i,j,k) = ea_s(i,j,k) + add_ent
1046  else
1047  ebtr(i,j,k-1) = eb_s(i,j,k-1) ; eatr(i,j,k) = ea_s(i,j,k)
1048  endif
1049 
1050  if (cs%double_diffuse) then ; if (kd_extra_s(i,j,k) > 0.0) then
1051  add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
1052  (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + h_neglect)
1053  ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
1054  eatr(i,j,k) = eatr(i,j,k) + add_ent
1055  endif ; endif
1056  enddo ; enddo
1057  do i=is,ie ; eatr(i,j,1) = ea_s(i,j,1) ; enddo
1058 
1059  enddo
1060 
1061  ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
1062  ! so h_prebound is used for the old thickness.
1063  !### I think that in the following, ea_s and eb_s should be eatr and ebtr.
1064  call call_tracer_column_fns(h_prebound, h, ea_s, eb_s, fluxes, hml, dt, g, gv, us, tv, &
1065  cs%optics, cs%tracer_flow_CSp, cs%debug, &
1066  evap_cfl_limit = cs%evap_CFL_limit, &
1067  minimum_forcing_depth=cs%minimum_forcing_depth)
1068 
1069  elseif (cs%double_diffuse) then ! extra diffusivity for passive tracers
1070 
1071  do j=js,je ; do i=is,ie
1072  ebtr(i,j,nz) = eb_s(i,j,nz) ; eatr(i,j,1) = ea_s(i,j,1)
1073  enddo ; enddo
1074  !$OMP parallel do default(shared) private(add_ent)
1075  do k=nz,2,-1 ; do j=js,je ; do i=is,ie
1076  if (kd_extra_s(i,j,k) > 0.0) then
1077  add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
1078  (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + h_neglect)
1079  else
1080  add_ent = 0.0
1081  endif
1082  ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
1083  eatr(i,j,k) = ea_s(i,j,k) + add_ent
1084  enddo ; enddo ; enddo
1085 
1086  ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
1087  call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
1088  cs%optics, cs%tracer_flow_CSp, cs%debug,&
1089  evap_cfl_limit = cs%evap_CFL_limit, &
1090  minimum_forcing_depth=cs%minimum_forcing_depth)
1091  else
1092  ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
1093  !### eatr and ebtr may not be initialized or may be 0, depending on CS%geothermal.
1094  call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
1095  cs%optics, cs%tracer_flow_CSp, cs%debug, &
1096  evap_cfl_limit = cs%evap_CFL_limit, &
1097  minimum_forcing_depth=cs%minimum_forcing_depth)
1098  endif ! (CS%mix_boundary_tracers)
1099 
1100  call cpu_clock_end(id_clock_tracers)
1101 
1102  ! Apply ALE sponge
1103  if (cs%use_sponge) then
1104  call cpu_clock_begin(id_clock_sponge)
1105  if (associated(cs%ALE_sponge_CSp)) then
1106  call apply_ale_sponge(h, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)
1107  endif
1108 
1109  call cpu_clock_end(id_clock_sponge)
1110  if (cs%debug) then
1111  call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)
1112  call mom_thermovar_chksum("apply_sponge ", tv, g)
1113  endif
1114  endif ! CS%use_sponge
1115 
1116  call disable_averaging(cs%diag)
1117  ! Diagnose the diapycnal diffusivities and other related quantities.
1118  call enable_averages(dt, time_end, cs%diag)
1119 
1120  if (cs%id_Kd_interface > 0) call post_data(cs%id_Kd_interface, kd_int, cs%diag)
1121  if (cs%id_Kd_heat > 0) call post_data(cs%id_Kd_heat, kd_heat, cs%diag)
1122  if (cs%id_Kd_salt > 0) call post_data(cs%id_Kd_salt, kd_salt, cs%diag)
1123  if (cs%id_Kd_ePBL > 0) call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)
1124 
1125  if (cs%id_ea > 0) call post_data(cs%id_ea, ea_s, cs%diag)
1126  if (cs%id_eb > 0) call post_data(cs%id_eb, eb_s, cs%diag)
1127  if (cs%id_ea_t > 0) call post_data(cs%id_ea_t, ea_t, cs%diag)
1128  if (cs%id_eb_t > 0) call post_data(cs%id_eb_t, eb_t, cs%diag)
1129  if (cs%id_ea_s > 0) call post_data(cs%id_ea_s, ea_s, cs%diag)
1130  if (cs%id_eb_s > 0) call post_data(cs%id_eb_s, eb_s, cs%diag)
1131 
1132  if (cs%id_dudt_dia > 0) call post_data(cs%id_dudt_dia, adp%du_dt_dia, cs%diag)
1133  if (cs%id_dvdt_dia > 0) call post_data(cs%id_dvdt_dia, adp%dv_dt_dia, cs%diag)
1134  if (cs%id_wd > 0) call post_data(cs%id_wd, cdp%diapyc_vel, cs%diag)
1135 
1136  if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)
1137  if (cs%id_Tadv > 0) call post_data(cs%id_Tadv, tadv_flx, cs%diag)
1138  if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)
1139  if (cs%id_Sadv > 0) call post_data(cs%id_Sadv, sadv_flx, cs%diag)
1140 
1141  !! Diagnostics for terms multiplied by fractional thicknesses
1142  if (cs%id_hf_dudt_dia_2d > 0) then
1143  allocate(hf_dudt_dia_2d(g%IsdB:g%IedB,g%jsd:g%jed))
1144  hf_dudt_dia_2d(:,:) = 0.0
1145  do k=1,nz ; do j=js,je ; do i=isq,ieq
1146  hf_dudt_dia_2d(i,j) = hf_dudt_dia_2d(i,j) + adp%du_dt_dia(i,j,k) * adp%diag_hfrac_u(i,j,k)
1147  enddo ; enddo ; enddo
1148  call post_data(cs%id_hf_dudt_dia_2d, hf_dudt_dia_2d, cs%diag)
1149  deallocate(hf_dudt_dia_2d)
1150  endif
1151 
1152  if (cs%id_hf_dvdt_dia_2d > 0) then
1153  allocate(hf_dvdt_dia_2d(g%isd:g%ied,g%JsdB:g%JedB))
1154  hf_dvdt_dia_2d(:,:) = 0.0
1155  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1156  hf_dvdt_dia_2d(i,j) = hf_dvdt_dia_2d(i,j) + adp%dv_dt_dia(i,j,k) * adp%diag_hfrac_v(i,j,k)
1157  enddo ; enddo ; enddo
1158  call post_data(cs%id_hf_dvdt_dia_2d, hf_dvdt_dia_2d, cs%diag)
1159  deallocate(hf_dvdt_dia_2d)
1160  endif
1161 
1162  call disable_averaging(cs%diag)
1163 
1164  if (showcalltree) call calltree_leave("diabatic_ALE_legacy()")
1165 

◆ diabatic_driver_end()

subroutine, public mom_diabatic_driver::diabatic_driver_end ( type(diabatic_cs), pointer  CS)

Routine to close the diabatic driver module.

Parameters
csmodule control structure

Definition at line 3733 of file MOM_diabatic_driver.F90.

3733  type(diabatic_CS), pointer :: CS !< module control structure
3734 
3735  if (.not.associated(cs)) return
3736 
3737  call diabatic_aux_end(cs%diabatic_aux_CSp)
3738 
3739  call entrain_diffusive_end(cs%entrain_diffusive_CSp)
3740  call set_diffusivity_end(cs%set_diff_CSp)
3741 
3742  if (cs%useKPP) then
3743  deallocate( cs%KPP_buoy_flux )
3744  deallocate( cs%KPP_temp_flux )
3745  deallocate( cs%KPP_salt_flux )
3746  deallocate( cs%KPP_NLTheat )
3747  deallocate( cs%KPP_NLTscalar )
3748  call kpp_end(cs%KPP_CSp)
3749  endif
3750 
3751  if (cs%use_CVMix_conv) call cvmix_conv_end(cs%CVMix_conv_csp)
3752 
3753  if (cs%use_energetic_PBL) &
3754  call energetic_pbl_end(cs%energetic_PBL_CSp)
3755  if (cs%debug_energy_req) &
3756  call diapyc_energy_req_end(cs%diapyc_en_rec_CSp)
3757 
3758  if (associated(cs%optics)) then
3759  call opacity_end(cs%opacity_CSp, cs%optics)
3760  deallocate(cs%optics)
3761  endif
3762 
3763  ! GMM, the following is commented out because arrays in
3764  ! CS%diag_grids_prev are neither pointers or allocatables
3765  ! and, therefore, cannot be deallocated.
3766 
3767  !call diag_grid_storage_end(CS%diag_grids_prev)
3768 
3769  deallocate(cs)
3770 

◆ diabatic_driver_init()

subroutine, public mom_diabatic_driver::diabatic_driver_init ( type(time_type), target  Time,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
logical, intent(in)  useALEalgorithm,
type(diag_ctrl), intent(inout), target  diag,
type(accel_diag_ptrs), intent(inout)  ADp,
type(cont_diag_ptrs), intent(inout)  CDp,
type(diabatic_cs), pointer  CS,
type(tracer_flow_control_cs), pointer  tracer_flow_CSp,
type(sponge_cs), pointer  sponge_CSp,
type(ale_sponge_cs), pointer  ALE_sponge_CSp 
)

This routine initializes the diabatic driver module.

Parameters
timemodel time
[in,out]gmodel grid structure
[in]gvmodel vertical grid structure
[in]usA dimensional unit scaling type
[in]param_filefile to parse for parameter values
[in]usealealgorithmlogical for whether to use ALE remapping
[in,out]diagstructure to regulate diagnostic output
[in,out]adppointers to accelerations in momentum equations, to enable diagnostics, like energy budgets
[in,out]cdppointers to terms in continuity equations
csmodule control structure
tracer_flow_csppointer to control structure of the tracer flow control module
sponge_csppointer to the sponge module control structure
ale_sponge_csppointer to the ALE sponge module control structure

Definition at line 3147 of file MOM_diabatic_driver.F90.

3147  type(time_type), target :: Time !< model time
3148  type(ocean_grid_type), intent(inout) :: G !< model grid structure
3149  type(verticalGrid_type), intent(in) :: GV !< model vertical grid structure
3150  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3151  type(param_file_type), intent(in) :: param_file !< file to parse for parameter values
3152  logical, intent(in) :: useALEalgorithm !< logical for whether to use ALE remapping
3153  type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output
3154  type(accel_diag_ptrs), intent(inout) :: ADp !< pointers to accelerations in momentum equations,
3155  !! to enable diagnostics, like energy budgets
3156  type(cont_diag_ptrs), intent(inout) :: CDp !< pointers to terms in continuity equations
3157  type(diabatic_CS), pointer :: CS !< module control structure
3158  type(tracer_flow_control_CS), pointer :: tracer_flow_CSp !< pointer to control structure of the
3159  !! tracer flow control module
3160  type(sponge_CS), pointer :: sponge_CSp !< pointer to the sponge module control structure
3161  type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< pointer to the ALE sponge module control structure
3162 
3163  real :: Kd ! A diffusivity used in the default for other tracer diffusivities, in MKS units [m2 s-1]
3164  integer :: num_mode
3165  logical :: use_temperature
3166  character(len=20) :: EN1, EN2, EN3
3167 
3168 ! This "include" declares and sets the variable "version".
3169 #include "version_variable.h"
3170  character(len=40) :: mdl = "MOM_diabatic_driver" ! This module's name.
3171  character(len=48) :: thickness_units
3172  character(len=40) :: var_name
3173  character(len=160) :: var_descript
3174  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nbands, m
3175  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
3176  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3177 
3178  if (associated(cs)) then
3179  call mom_error(warning, "diabatic_driver_init called with an "// &
3180  "associated control structure.")
3181  return
3182  else
3183  allocate(cs)
3184  endif
3185 
3186  cs%diag => diag
3187  cs%Time => time
3188 
3189  if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
3190  if (associated(sponge_csp)) cs%sponge_CSp => sponge_csp
3191  if (associated(ale_sponge_csp)) cs%ALE_sponge_CSp => ale_sponge_csp
3192 
3193  cs%useALEalgorithm = usealealgorithm
3194  cs%bulkmixedlayer = (gv%nkml > 0)
3195 
3196  ! Set default, read and log parameters
3197  call log_version(param_file, mdl, version, &
3198  "The following parameters are used for diabatic processes.", &
3199  log_to_all=.true., debugging=.true.)
3200  call get_param(param_file, mdl, "USE_LEGACY_DIABATIC_DRIVER", cs%use_legacy_diabatic, &
3201  "If true, use a legacy version of the diabatic subroutine. "//&
3202  "This is temporary and is needed to avoid change in answers.", &
3203  default=.true.)
3204  call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &
3205  "If true, sponges may be applied anywhere in the domain. "//&
3206  "The exact location and properties of those sponges are "//&
3207  "specified via calls to initialize_sponge and possibly "//&
3208  "set_up_sponge_field.", default=.false.)
3209  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
3210  "If true, temperature and salinity are used as state "//&
3211  "variables.", default=.true.)
3212  call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", cs%use_energetic_PBL, &
3213  "If true, use an implied energetics planetary boundary "//&
3214  "layer scheme to determine the diffusivity and viscosity "//&
3215  "in the surface boundary layer.", default=.false.)
3216  call get_param(param_file, mdl, "EPBL_IS_ADDITIVE", cs%ePBL_is_additive, &
3217  "If true, the diffusivity from ePBL is added to all "//&
3218  "other diffusivities. Otherwise, the larger of kappa-shear "//&
3219  "and ePBL diffusivities are used.", default=.true.)
3220  call get_param(param_file, mdl, "PRANDTL_EPBL", cs%ePBL_Prandtl, &
3221  "The Prandtl number used by ePBL to convert vertical diffusivities into "//&
3222  "viscosities.", default=1.0, units="nondim", do_not_log=.not.cs%use_energetic_PBL)
3223  call get_param(param_file, mdl, "USE_KPP", cs%use_KPP, &
3224  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
3225  "to calculate diffusivities and non-local transport in the OBL.", &
3226  default=.false., do_not_log=.true.)
3227  cs%use_CVMix_ddiff = cvmix_ddiff_is_used(param_file)
3228 
3229  cs%use_kappa_shear = kappa_shear_is_used(param_file)
3230  cs%use_CVMix_shear = cvmix_shear_is_used(param_file)
3231 
3232  if (cs%bulkmixedlayer) then
3233  call get_param(param_file, mdl, "ML_MIX_FIRST", cs%ML_mix_first, &
3234  "The fraction of the mixed layer mixing that is applied "//&
3235  "before interior diapycnal mixing. 0 by default.", &
3236  units="nondim", default=0.0)
3237  call get_param(param_file, mdl, "NKBL", cs%nkbl, default=2, do_not_log=.true.)
3238  else
3239  cs%ML_mix_first = 0.0
3240  endif
3241  if (use_temperature) then
3242  call get_param(param_file, mdl, "DO_GEOTHERMAL", cs%use_geothermal, &
3243  "If true, apply geothermal heating.", default=.false.)
3244  else
3245  cs%use_geothermal = .false.
3246  endif
3247  call get_param(param_file, mdl, "INTERNAL_TIDES", cs%use_int_tides, &
3248  "If true, use the code that advances a separate set of "//&
3249  "equations for the internal tide energy density.", default=.false.)
3250  cs%nMode = 1
3251  if (cs%use_int_tides) then
3252  call get_param(param_file, mdl, "INTERNAL_TIDE_MODES", cs%nMode, &
3253  "The number of distinct internal tide modes "//&
3254  "that will be calculated.", default=1, do_not_log=.true.)
3255  call get_param(param_file, mdl, "UNIFORM_TEST_CG", cs%uniform_test_cg, &
3256  "If positive, a uniform group velocity of internal tide for test case", &
3257  default=-1., units="m s-1", scale=us%m_s_to_L_T)
3258  endif
3259 
3260  call get_param(param_file, mdl, "MASSLESS_MATCH_TARGETS", &
3261  cs%massless_match_targets, &
3262  "If true, the temperature and salinity of massless layers "//&
3263  "are kept consistent with their target densities. "//&
3264  "Otherwise the properties of massless layers evolve "//&
3265  "diffusively to match massive neighboring layers.", &
3266  default=.true.)
3267 
3268  call get_param(param_file, mdl, "AGGREGATE_FW_FORCING", cs%aggregate_FW_forcing, &
3269  "If true, the net incoming and outgoing fresh water fluxes are combined "//&
3270  "and applied as either incoming or outgoing depending on the sign of the net. "//&
3271  "If false, the net incoming fresh water flux is added to the model and "//&
3272  "thereafter the net outgoing is removed from the topmost non-vanished "//&
3273  "layers of the updated state.", default=.true.)
3274 
3275  call get_param(param_file, mdl, "DEBUG", cs%debug, &
3276  "If true, write out verbose debugging data.", &
3277  default=.false., debuggingparam=.true.)
3278  call get_param(param_file, mdl, "DEBUG_CONSERVATION", cs%debugConservation, &
3279  "If true, monitor conservation and extrema.", &
3280  default=.false., debuggingparam=.true.)
3281 
3282  call get_param(param_file, mdl, "DEBUG_ENERGY_REQ", cs%debug_energy_req, &
3283  "If true, debug the energy requirements.", default=.false., do_not_log=.true.)
3284  call get_param(param_file, mdl, "MIX_BOUNDARY_TRACERS", cs%mix_boundary_tracers, &
3285  "If true, mix the passive tracers in massless layers at "//&
3286  "the bottom into the interior as though a diffusivity of "//&
3287  "KD_MIN_TR were operating.", default=.true.)
3288 
3289  if (cs%mix_boundary_tracers) then
3290  call get_param(param_file, mdl, "KD", kd, default=0.0)
3291  call get_param(param_file, mdl, "KD_MIN_TR", cs%Kd_min_tr, &
3292  "A minimal diffusivity that should always be applied to "//&
3293  "tracers, especially in massless layers near the bottom. "//&
3294  "The default is 0.1*KD.", units="m2 s-1", default=0.1*kd, scale=us%m2_s_to_Z2_T)
3295  call get_param(param_file, mdl, "KD_BBL_TR", cs%Kd_BBL_tr, &
3296  "A bottom boundary layer tracer diffusivity that will "//&
3297  "allow for explicitly specified bottom fluxes. The "//&
3298  "entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) "//&
3299  "over the same distance.", units="m2 s-1", default=0., scale=us%m2_s_to_Z2_T)
3300  endif
3301 
3302  call get_param(param_file, mdl, "TRACER_TRIDIAG", cs%tracer_tridiag, &
3303  "If true, use the passive tracer tridiagonal solver for T and S", &
3304  default=.false.)
3305 
3306  call get_param(param_file, mdl, "MINIMUM_FORCING_DEPTH", cs%minimum_forcing_depth, &
3307  "The smallest depth over which forcing can be applied. This "//&
3308  "only takes effect when near-surface layers become thin "//&
3309  "relative to this scale, in which case the forcing tendencies "//&
3310  "scaled down by distributing the forcing over this depth scale.", &
3311  units="m", default=0.001, scale=gv%m_to_H)
3312  call get_param(param_file, mdl, "EVAP_CFL_LIMIT", cs%evap_CFL_limit, &
3313  "The largest fraction of a layer than can be lost to forcing "//&
3314  "(e.g. evaporation, sea-ice formation) in one time-step. The unused "//&
3315  "mass loss is passed down through the column.", &
3316  units="nondim", default=0.8)
3317 
3318 
3319  ! Register all available diagnostics for this module.
3320  thickness_units = get_thickness_units(gv)
3321 
3322  cs%id_ea_t = register_diag_field('ocean_model', 'ea_t', diag%axesTL, time, &
3323  'Layer (heat) entrainment from above per timestep', 'm', conversion=gv%H_to_m)
3324  cs%id_eb_t = register_diag_field('ocean_model', 'eb_t', diag%axesTL, time, &
3325  'Layer (heat) entrainment from below per timestep', 'm', conversion=gv%H_to_m)
3326  cs%id_ea_s = register_diag_field('ocean_model', 'ea_s', diag%axesTL, time, &
3327  'Layer (salt) entrainment from above per timestep', 'm', conversion=gv%H_to_m)
3328  cs%id_eb_s = register_diag_field('ocean_model', 'eb_s', diag%axesTL, time, &
3329  'Layer (salt) entrainment from below per timestep', 'm', conversion=gv%H_to_m)
3330  ! used by layer diabatic
3331  cs%id_ea = register_diag_field('ocean_model', 'ea', diag%axesTL, time, &
3332  'Layer entrainment from above per timestep', 'm', conversion=gv%H_to_m)
3333  cs%id_eb = register_diag_field('ocean_model', 'eb', diag%axesTL, time, &
3334  'Layer entrainment from below per timestep', 'm', conversion=gv%H_to_m)
3335  cs%id_wd = register_diag_field('ocean_model', 'wd', diag%axesTi, time, &
3336  'Diapycnal velocity', 'm s-1', conversion=gv%H_to_m)
3337  if (cs%id_wd > 0) call safe_alloc_ptr(cdp%diapyc_vel,isd,ied,jsd,jed,nz+1)
3338 
3339  cs%id_dudt_dia = register_diag_field('ocean_model', 'dudt_dia', diag%axesCuL, time, &
3340  'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)
3341  cs%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, time, &
3342  'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)
3343 
3344  cs%id_hf_dudt_dia_2d = register_diag_field('ocean_model', 'hf_dudt_dia_2d', diag%axesCu1, time, &
3345  'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Diapycnal Mixing', 'm s-2', &
3346  conversion=us%L_T2_to_m_s2)
3347  if (cs%id_hf_dudt_dia_2d > 0) then
3348  call safe_alloc_ptr(adp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
3349  call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
3350  endif
3351 
3352  cs%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, time, &
3353  'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', 'm s-2', &
3354  conversion=us%L_T2_to_m_s2)
3355  if (cs%id_hf_dvdt_dia_2d > 0) then
3356  call safe_alloc_ptr(adp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
3357  call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
3358  endif
3359 
3360  if (cs%use_int_tides) then
3361  cs%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, &
3362  time, 'First baroclinic mode (eigen) speed', 'm s-1', conversion=us%L_T_to_m_s)
3363  allocate(cs%id_cn(cs%nMode)) ; cs%id_cn(:) = -1
3364  do m=1,cs%nMode
3365  write(var_name, '("cn_mode",i1)') m
3366  write(var_descript, '("Baroclinic (eigen) speed of mode ",i1)') m
3367  cs%id_cn(m) = register_diag_field('ocean_model',var_name, diag%axesT1, &
3368  time, var_descript, 'm s-1', conversion=us%L_T_to_m_s)
3369  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
3370  enddo
3371  endif
3372 
3373  if (use_temperature) then
3374  cs%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff", diag%axesTi, &
3375  time, "Diffusive diapycnal temperature flux across interfaces", &
3376  "degC m s-1", conversion=gv%H_to_m*us%s_to_T)
3377  cs%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv", diag%axesTi, &
3378  time, "Advective diapycnal temperature flux across interfaces", &
3379  "degC m s-1", conversion=gv%H_to_m*us%s_to_T)
3380  cs%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff", diag%axesTi, &
3381  time, "Diffusive diapycnal salnity flux across interfaces", &
3382  "psu m s-1", conversion=gv%H_to_m*us%s_to_T)
3383  cs%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv", diag%axesTi, &
3384  time, "Advective diapycnal salnity flux across interfaces", &
3385  "psu m s-1", conversion=gv%H_to_m*us%s_to_T)
3386  cs%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, time, &
3387  'Mixed layer depth (delta rho = 0.03)', 'm', conversion=us%Z_to_m, &
3388  cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', &
3389  cmor_standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t')
3390  cs%id_mlotstsq = register_diag_field('ocean_model', 'mlotstsq', diag%axesT1, time, &
3391  long_name='Square of Ocean Mixed Layer Thickness Defined by Sigma T', &
3392  standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t', &
3393  units='m2', conversion=us%Z_to_m**2)
3394  cs%id_MLD_0125 = register_diag_field('ocean_model', 'MLD_0125', diag%axesT1, time, &
3395  'Mixed layer depth (delta rho = 0.125)', 'm', conversion=us%Z_to_m)
3396  call get_param(param_file, mdl, "MLD_EN_VALS", cs%MLD_EN_VALS, &
3397  "The energy values used to compute MLDs. If not set (or all set to 0.), the "//&
3398  "default will overwrite to 25., 2500., 250000.",units='J/m2', default=0., &
3399  scale=us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**2)
3400  if ((cs%MLD_EN_VALS(1)==0.).and.(cs%MLD_EN_VALS(2)==0.).and.(cs%MLD_EN_VALS(3)==0.)) then
3401  cs%MLD_EN_VALS = (/25.*us%kg_m3_to_R*us%m_to_Z*us%m_to_L**2*us%T_to_s**2,&
3402  2500.*us%kg_m3_to_R*us%m_to_Z*us%m_to_L**2*us%T_to_s**2,&
3403  250000.*us%kg_m3_to_R*us%m_to_Z*us%m_to_L**2*us%T_to_s**2/)
3404  endif
3405  write(en1,'(F10.2)') cs%MLD_EN_VALS(1)*us%R_to_kg_m3*us%Z_to_m*us%L_to_m**2*us%s_to_T**2
3406  write(en2,'(F10.2)') cs%MLD_EN_VALS(2)*us%R_to_kg_m3*us%Z_to_m*us%L_to_m**2*us%s_to_T**2
3407  write(en3,'(F10.2)') cs%MLD_EN_VALS(3)*us%R_to_kg_m3*us%Z_to_m*us%L_to_m**2*us%s_to_T**2
3408  cs%id_MLD_EN1 = register_diag_field('ocean_model', 'MLD_EN1', diag%axesT1, time, &
3409  'Mixed layer depth for energy value set to '//trim(en1)//' J/m2 (Energy set by 1st MLD_EN_VALS)', &
3410  'm', conversion=us%Z_to_m)
3411  cs%id_MLD_EN2 = register_diag_field('ocean_model', 'MLD_EN2', diag%axesT1, time, &
3412  'Mixed layer depth for energy value set to '//trim(en2)//' J/m2 (Energy set by 2nd MLD_EN_VALS)', &
3413  'm', conversion=us%Z_to_m)
3414  cs%id_MLD_EN3 = register_diag_field('ocean_model', 'MLD_EN3', diag%axesT1, time, &
3415  'Mixed layer depth for energy value set to '//trim(en3)//' J/m2 (Energy set by 3rd MLD_EN_VALS)', &
3416  'm', conversion=us%Z_to_m)
3417  cs%id_subMLN2 = register_diag_field('ocean_model', 'subML_N2', diag%axesT1, time, &
3418  'Squared buoyancy frequency below mixed layer', 's-2', conversion=us%s_to_T**2)
3419  cs%id_MLD_user = register_diag_field('ocean_model', 'MLD_user', diag%axesT1, time, &
3420  'Mixed layer depth (used defined)', 'm', conversion=us%Z_to_m)
3421  endif
3422  call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", cs%MLDdensityDifference, &
3423  "The density difference used to determine a diagnostic mixed "//&
3424  "layer depth, MLD_user, following the definition of Levitus 1982. "//&
3425  "The MLD is the depth at which the density is larger than the "//&
3426  "surface density by the specified amount.", &
3427  units='kg/m3', default=0.1, scale=us%kg_m3_to_R)
3428  call get_param(param_file, mdl, "DIAG_DEPTH_SUBML_N2", cs%dz_subML_N2, &
3429  "The distance over which to calculate a diagnostic of the "//&
3430  "stratification at the base of the mixed layer.", &
3431  units='m', default=50.0, scale=us%m_to_Z)
3432 
3433  if (cs%id_dudt_dia > 0) call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
3434  if (cs%id_dvdt_dia > 0) call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
3435 
3436  ! diagnostics for values prior to diabatic and prior to ALE
3437  cs%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, time, &
3438  'Zonal velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)
3439  cs%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, time, &
3440  'Meridional velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)
3441  cs%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, time, &
3442  'Layer Thickness before diabatic forcing', &
3443  trim(thickness_units), conversion=gv%H_to_MKS, v_extensive=.true.)
3444  cs%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, time, &
3445  'Interface Heights before diabatic forcing', 'm')
3446  if (use_temperature) then
3447  cs%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, time, &
3448  'Potential Temperature', 'degC')
3449  cs%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, time, &
3450  'Salinity', 'PSU')
3451  endif
3452 
3453 
3454  !call set_diffusivity_init(Time, G, param_file, diag, CS%set_diff_CSp, CS%int_tide_CSp)
3455  cs%id_Kd_interface = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, time, &
3456  'Total diapycnal diffusivity at interfaces', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
3457  if (cs%use_energetic_PBL) then
3458  cs%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, time, &
3459  'ePBL diapycnal diffusivity at interfaces', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
3460  endif
3461 
3462  cs%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, time, &
3463  'Total diapycnal diffusivity for heat at interfaces', 'm2 s-1', conversion=us%Z2_T_to_m2_s, &
3464  cmor_field_name='difvho', &
3465  cmor_standard_name='ocean_vertical_heat_diffusivity', &
3466  cmor_long_name='Ocean vertical heat diffusivity')
3467  cs%id_Kd_salt = register_diag_field('ocean_model', 'Kd_salt', diag%axesTi, time, &
3468  'Total diapycnal diffusivity for salt at interfaces', 'm2 s-1', conversion=us%Z2_T_to_m2_s, &
3469  cmor_field_name='difvso', &
3470  cmor_standard_name='ocean_vertical_salt_diffusivity', &
3471  cmor_long_name='Ocean vertical salt diffusivity')
3472 
3473  ! CS%useKPP is set to True if KPP-scheme is to be used, False otherwise.
3474  ! KPP_init() allocated CS%KPP_Csp and also sets CS%KPPisPassive
3475  cs%useKPP = kpp_init(param_file, g, gv, us, diag, time, cs%KPP_CSp, passive=cs%KPPisPassive)
3476  if (cs%useKPP) then
3477  allocate( cs%KPP_NLTheat(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_NLTheat(:,:,:) = 0.
3478  allocate( cs%KPP_NLTscalar(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_NLTscalar(:,:,:) = 0.
3479  endif
3480  if (cs%useKPP) then
3481  allocate( cs%KPP_buoy_flux(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_buoy_flux(:,:,:) = 0.
3482  allocate( cs%KPP_temp_flux(isd:ied,jsd:jed) ) ; cs%KPP_temp_flux(:,:) = 0.
3483  allocate( cs%KPP_salt_flux(isd:ied,jsd:jed) ) ; cs%KPP_salt_flux(:,:) = 0.
3484  endif
3485 
3486 
3487  ! diagnostics for tendencies of temp and saln due to diabatic processes
3488  ! available only for ALE algorithm.
3489  ! diagnostics for tendencies of temp and heat due to frazil
3490  cs%id_diabatic_diff_h = register_diag_field('ocean_model', 'diabatic_diff_h', diag%axesTL, time, &
3491  long_name='Cell thickness used during diabatic diffusion', &
3492  units='m', conversion=gv%H_to_m, v_extensive=.true.)
3493  if (cs%useALEalgorithm) then
3494  cs%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', &
3495  'diabatic_diff_temp_tendency', diag%axesTL, time, &
3496  'Diabatic diffusion temperature tendency', 'degC s-1', conversion=us%s_to_T)
3497  if (cs%id_diabatic_diff_temp_tend > 0) then
3498  cs%diabatic_diff_tendency_diag = .true.
3499  endif
3500 
3501  cs%id_diabatic_diff_saln_tend = register_diag_field('ocean_model',&
3502  'diabatic_diff_saln_tendency', diag%axesTL, time, &
3503  'Diabatic diffusion salinity tendency', 'psu s-1', conversion=us%s_to_T)
3504  if (cs%id_diabatic_diff_saln_tend > 0) then
3505  cs%diabatic_diff_tendency_diag = .true.
3506  endif
3507 
3508  cs%id_diabatic_diff_heat_tend = register_diag_field('ocean_model', &
3509  'diabatic_heat_tendency', diag%axesTL, time, &
3510  'Diabatic diffusion heat tendency', &
3511  'W m-2', conversion=us%QRZ_T_to_W_m2, cmor_field_name='opottempdiff', &
3512  cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'// &
3513  'due_to_parameterized_dianeutral_mixing', &
3514  cmor_long_name='Tendency of sea water potential temperature expressed as heat content '// &
3515  'due to parameterized dianeutral mixing', &
3516  v_extensive=.true.)
3517  if (cs%id_diabatic_diff_heat_tend > 0) then
3518  cs%diabatic_diff_tendency_diag = .true.
3519  endif
3520 
3521  cs%id_diabatic_diff_salt_tend = register_diag_field('ocean_model', &
3522  'diabatic_salt_tendency', diag%axesTL, time, &
3523  'Diabatic diffusion of salt tendency', &
3524  'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s, cmor_field_name='osaltdiff', &
3525  cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &
3526  'due_to_parameterized_dianeutral_mixing', &
3527  cmor_long_name='Tendency of sea water salinity expressed as salt content '// &
3528  'due to parameterized dianeutral mixing', &
3529  v_extensive=.true.)
3530  if (cs%id_diabatic_diff_salt_tend > 0) then
3531  cs%diabatic_diff_tendency_diag = .true.
3532  endif
3533 
3534  ! This diagnostic should equal to roundoff if all is working well.
3535  cs%id_diabatic_diff_heat_tend_2d = register_diag_field('ocean_model', &
3536  'diabatic_heat_tendency_2d', diag%axesT1, time, &
3537  'Depth integrated diabatic diffusion heat tendency', &
3538  'W m-2', conversion=us%QRZ_T_to_W_m2, cmor_field_name='opottempdiff_2d', &
3539  cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'//&
3540  'due_to_parameterized_dianeutral_mixing_depth_integrated', &
3541  cmor_long_name='Tendency of sea water potential temperature expressed as heat content '//&
3542  'due to parameterized dianeutral mixing depth integrated')
3543  if (cs%id_diabatic_diff_heat_tend_2d > 0) then
3544  cs%diabatic_diff_tendency_diag = .true.
3545  endif
3546 
3547  ! This diagnostic should equal to roundoff if all is working well.
3548  cs%id_diabatic_diff_salt_tend_2d = register_diag_field('ocean_model', &
3549  'diabatic_salt_tendency_2d', diag%axesT1, time, &
3550  'Depth integrated diabatic diffusion salt tendency', &
3551  'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s, cmor_field_name='osaltdiff_2d', &
3552  cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &
3553  'due_to_parameterized_dianeutral_mixing_depth_integrated', &
3554  cmor_long_name='Tendency of sea water salinity expressed as salt content '// &
3555  'due to parameterized dianeutral mixing depth integrated')
3556  if (cs%id_diabatic_diff_salt_tend_2d > 0) then
3557  cs%diabatic_diff_tendency_diag = .true.
3558  endif
3559 
3560  ! diagnostics for tendencies of thickness temp and saln due to boundary forcing
3561  ! available only for ALE algorithm.
3562  ! diagnostics for tendencies of temp and heat due to frazil
3563  cs%id_boundary_forcing_h = register_diag_field('ocean_model', 'boundary_forcing_h', diag%axesTL, time, &
3564  long_name='Cell thickness after applying boundary forcing', &
3565  units='m', conversion=gv%H_to_m, v_extensive=.true.)
3566  cs%id_boundary_forcing_h_tendency = register_diag_field('ocean_model', &
3567  'boundary_forcing_h_tendency', diag%axesTL, time, &
3568  'Cell thickness tendency due to boundary forcing', 'm s-1', &
3569  conversion=gv%H_to_m*us%s_to_T, v_extensive=.true.)
3570  if (cs%id_boundary_forcing_h_tendency > 0) then
3571  cs%boundary_forcing_tendency_diag = .true.
3572  endif
3573 
3574  cs%id_boundary_forcing_temp_tend = register_diag_field('ocean_model',&
3575  'boundary_forcing_temp_tendency', diag%axesTL, time, &
3576  'Boundary forcing temperature tendency', 'degC s-1', conversion=us%s_to_T)
3577  if (cs%id_boundary_forcing_temp_tend > 0) then
3578  cs%boundary_forcing_tendency_diag = .true.
3579  endif
3580 
3581  cs%id_boundary_forcing_saln_tend = register_diag_field('ocean_model',&
3582  'boundary_forcing_saln_tendency', diag%axesTL, time, &
3583  'Boundary forcing saln tendency', 'psu s-1', conversion=us%s_to_T)
3584  if (cs%id_boundary_forcing_saln_tend > 0) then
3585  cs%boundary_forcing_tendency_diag = .true.
3586  endif
3587 
3588  cs%id_boundary_forcing_heat_tend = register_diag_field('ocean_model',&
3589  'boundary_forcing_heat_tendency', diag%axesTL, time, &
3590  'Boundary forcing heat tendency', &
3591  'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive = .true.)
3592  if (cs%id_boundary_forcing_heat_tend > 0) then
3593  cs%boundary_forcing_tendency_diag = .true.
3594  endif
3595 
3596  cs%id_boundary_forcing_salt_tend = register_diag_field('ocean_model',&
3597  'boundary_forcing_salt_tendency', diag%axesTL, time, &
3598  'Boundary forcing salt tendency', 'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s, &
3599  v_extensive = .true.)
3600  if (cs%id_boundary_forcing_salt_tend > 0) then
3601  cs%boundary_forcing_tendency_diag = .true.
3602  endif
3603 
3604  ! This diagnostic should equal to surface heat flux if all is working well.
3605  cs%id_boundary_forcing_heat_tend_2d = register_diag_field('ocean_model',&
3606  'boundary_forcing_heat_tendency_2d', diag%axesT1, time, &
3607  'Depth integrated boundary forcing of ocean heat', &
3608  'W m-2', conversion=us%QRZ_T_to_W_m2)
3609  if (cs%id_boundary_forcing_heat_tend_2d > 0) then
3610  cs%boundary_forcing_tendency_diag = .true.
3611  endif
3612 
3613  ! This diagnostic should equal to surface salt flux if all is working well.
3614  cs%id_boundary_forcing_salt_tend_2d = register_diag_field('ocean_model',&
3615  'boundary_forcing_salt_tendency_2d', diag%axesT1, time, &
3616  'Depth integrated boundary forcing of ocean salt', &
3617  'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s)
3618  if (cs%id_boundary_forcing_salt_tend_2d > 0) then
3619  cs%boundary_forcing_tendency_diag = .true.
3620  endif
3621  endif
3622 
3623  ! diagnostics for tendencies of temp and heat due to frazil
3624  cs%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, time, &
3625  long_name='Cell Thickness', standard_name='cell_thickness', &
3626  units='m', conversion=gv%H_to_m, v_extensive=.true.)
3627 
3628  ! diagnostic for tendency of temp due to frazil
3629  cs%id_frazil_temp_tend = register_diag_field('ocean_model',&
3630  'frazil_temp_tendency', diag%axesTL, time, &
3631  'Temperature tendency due to frazil formation', 'degC s-1', conversion=us%s_to_T)
3632  if (cs%id_frazil_temp_tend > 0) then
3633  cs%frazil_tendency_diag = .true.
3634  endif
3635 
3636  ! diagnostic for tendency of heat due to frazil
3637  cs%id_frazil_heat_tend = register_diag_field('ocean_model',&
3638  'frazil_heat_tendency', diag%axesTL, time, &
3639  'Heat tendency due to frazil formation', &
3640  'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive=.true.)
3641  if (cs%id_frazil_heat_tend > 0) then
3642  cs%frazil_tendency_diag = .true.
3643  endif
3644 
3645  ! if all is working propertly, this diagnostic should equal to hfsifrazil
3646  cs%id_frazil_heat_tend_2d = register_diag_field('ocean_model',&
3647  'frazil_heat_tendency_2d', diag%axesT1, time, &
3648  'Depth integrated heat tendency due to frazil formation', &
3649  'W m-2', conversion=us%QRZ_T_to_W_m2)
3650  if (cs%id_frazil_heat_tend_2d > 0) then
3651  cs%frazil_tendency_diag = .true.
3652  endif
3653 
3654  if (cs%frazil_tendency_diag) then
3655  allocate(cs%frazil_temp_diag(isd:ied,jsd:jed,nz) ) ; cs%frazil_temp_diag(:,:,:) = 0.
3656  allocate(cs%frazil_heat_diag(isd:ied,jsd:jed,nz) ) ; cs%frazil_heat_diag(:,:,:) = 0.
3657  endif
3658 
3659  ! CS%use_CVMix_conv is set to True if CVMix convection will be used, otherwise it is False.
3660  cs%use_CVMix_conv = cvmix_conv_init(time, g, gv, us, param_file, diag, cs%CVMix_conv_csp)
3661 
3662  call entrain_diffusive_init(time, g, gv, us, param_file, diag, cs%entrain_diffusive_CSp, &
3663  just_read_params=cs%useALEalgorithm)
3664 
3665  ! initialize the geothermal heating module
3666  if (cs%use_geothermal) &
3667  call geothermal_init(time, g, gv, us, param_file, diag, cs%geothermal_CSp, usealealgorithm)
3668 
3669  ! initialize module for internal tide induced mixing
3670  if (cs%use_int_tides) then
3671  call int_tide_input_init(time, g, gv, us, param_file, diag, cs%int_tide_input_CSp, &
3672  cs%int_tide_input)
3673  call internal_tides_init(time, g, gv, us, param_file, diag, cs%int_tide_CSp)
3674  endif
3675 
3676  ! initialize module for setting diffusivities
3677  call set_diffusivity_init(time, g, gv, us, param_file, diag, cs%set_diff_CSp, cs%int_tide_CSp, &
3678  halo_ts=cs%halo_TS_diff, double_diffuse=cs%double_diffuse)
3679 
3680  if (cs%useKPP .and. (cs%double_diffuse .and. .not.cs%use_CVMix_ddiff)) &
3681  call mom_error(fatal, 'diabatic_driver_init: DOUBLE_DIFFUSION (old method) does not work '//&
3682  'with KPP. Please set DOUBLE_DIFFUSION=False and USE_CVMIX_DDIFF=True.')
3683 
3684  ! set up the clocks for this module
3685  id_clock_entrain = cpu_clock_id('(Ocean diabatic entrain)', grain=clock_module)
3686  if (cs%bulkmixedlayer) &
3687  id_clock_mixedlayer = cpu_clock_id('(Ocean mixed layer)', grain=clock_module)
3688  id_clock_remap = cpu_clock_id('(Ocean vert remap)', grain=clock_module)
3689  if (cs%use_geothermal) &
3690  id_clock_geothermal = cpu_clock_id('(Ocean geothermal)', grain=clock_routine)
3691  id_clock_set_diffusivity = cpu_clock_id('(Ocean set_diffusivity)', grain=clock_module)
3692  id_clock_kpp = cpu_clock_id('(Ocean KPP)', grain=clock_module)
3693  id_clock_tracers = cpu_clock_id('(Ocean tracer_columns)', grain=clock_module_driver+5)
3694  if (cs%use_sponge) &
3695  id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=clock_module)
3696  id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=clock_routine)
3697  id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=clock_routine)
3698  id_clock_differential_diff = -1 ; if (cs%double_diffuse .and. .not.cs%use_CVMix_ddiff) &
3699  id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=clock_routine)
3700 
3701  ! initialize the auxiliary diabatic driver module
3702  call diabatic_aux_init(time, g, gv, us, param_file, diag, cs%diabatic_aux_CSp, &
3703  cs%useALEalgorithm, cs%use_energetic_PBL)
3704 
3705  ! initialize the boundary layer modules
3706  if (cs%bulkmixedlayer) &
3707  call bulkmixedlayer_init(time, g, gv, us, param_file, diag, cs%bulkmixedlayer_CSp)
3708  if (cs%use_energetic_PBL) &
3709  call energetic_pbl_init(time, g, gv, us, param_file, diag, cs%energetic_PBL_CSp)
3710 
3711  call regularize_layers_init(time, g, gv, param_file, diag, cs%regularize_layers_CSp)
3712 
3713  if (cs%debug_energy_req) &
3714  call diapyc_energy_req_init(time, g, gv, us, param_file, diag, cs%diapyc_en_rec_CSp)
3715 
3716  ! obtain information about the number of bands for penetrative shortwave
3717  if (use_temperature) then
3718  call get_param(param_file, mdl, "PEN_SW_NBANDS", nbands, default=1)
3719  if (nbands > 0) then
3720  allocate(cs%optics)
3721  call opacity_init(time, g, gv, us, param_file, diag, cs%opacity_CSp, cs%optics)
3722  endif
3723  endif
3724 
3725  ! Initialize the diagnostic grid storage
3726  call diag_grid_storage_init(cs%diag_grids_prev, g, diag)
3727 

◆ diagnose_boundary_forcing_tendency()

subroutine mom_diabatic_driver::diagnose_boundary_forcing_tendency ( type(thermo_var_ptrs), intent(in)  tv,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  temp_old,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  saln_old,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_old,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diabatic_cs), pointer  CS 
)
private

This routine diagnoses tendencies from application of boundary fluxes. These impacts are generally 3d, in particular for penetrative shortwave. Other fluxes contribute 3d in cases when the layers vanish or are very thin, in which case we distribute the flux into k > 1 layers.

Parameters
[in]gocean grid structure
[in]gvocean vertical grid structure
[in]tvpoints to updated thermodynamic fields
[in]hthickness after boundary flux application [H ~> m or kg m-2]
[in]temp_oldtemperature prior to boundary flux application [degC]
[in]saln_oldsalinity prior to boundary flux application [ppt]
[in]h_oldthickness prior to boundary flux application [H ~> m or kg m-2]
[in]dttime step [T ~> s]
[in]usA dimensional unit scaling type
csmodule control structure

Definition at line 2966 of file MOM_diabatic_driver.F90.

2966  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2967  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2968  type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields
2969  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2970  intent(in) :: h !< thickness after boundary flux application [H ~> m or kg m-2]
2971  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2972  intent(in) :: temp_old !< temperature prior to boundary flux application [degC]
2973  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2974  intent(in) :: saln_old !< salinity prior to boundary flux application [ppt]
2975  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2976  intent(in) :: h_old !< thickness prior to boundary flux application [H ~> m or kg m-2]
2977  real, intent(in) :: dt !< time step [T ~> s]
2978  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2979  type(diabatic_CS), pointer :: CS !< module control structure
2980 
2981  ! Local variables
2982  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
2983  real, dimension(SZI_(G),SZJ_(G)) :: work_2d
2984  real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
2985  real :: ppt2mks = 0.001 ! Conversion factor from g/kg to kg/kg.
2986  integer :: i, j, k, is, ie, js, je, nz
2987 
2988  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2989  idt = 0.0 ; if (dt > 0.0) idt = 1. / dt
2990  work_3d(:,:,:) = 0.0
2991  work_2d(:,:) = 0.0
2992 
2993  ! Thickness tendency
2994  if (cs%id_boundary_forcing_h_tendency > 0) then
2995  do k=1,nz ; do j=js,je ; do i=is,ie
2996  work_3d(i,j,k) = (h(i,j,k) - h_old(i,j,k))*idt
2997  enddo ; enddo ; enddo
2998  call post_data(cs%id_boundary_forcing_h_tendency, work_3d, cs%diag, alt_h=h_old)
2999  endif
3000 
3001  ! temperature tendency
3002  if (cs%id_boundary_forcing_temp_tend > 0) then
3003  do k=1,nz ; do j=js,je ; do i=is,ie
3004  work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt
3005  enddo ; enddo ; enddo
3006  call post_data(cs%id_boundary_forcing_temp_tend, work_3d, cs%diag, alt_h=h_old)
3007  endif
3008 
3009  ! heat tendency
3010  if (cs%id_boundary_forcing_heat_tend > 0 .or. cs%id_boundary_forcing_heat_tend_2d > 0) then
3011  do k=1,nz ; do j=js,je ; do i=is,ie
3012  work_3d(i,j,k) = gv%H_to_RZ * tv%C_p * idt * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * temp_old(i,j,k))
3013  enddo ; enddo ; enddo
3014  if (cs%id_boundary_forcing_heat_tend > 0) then
3015  call post_data(cs%id_boundary_forcing_heat_tend, work_3d, cs%diag, alt_h = h_old)
3016  endif
3017  if (cs%id_boundary_forcing_heat_tend_2d > 0) then
3018  do j=js,je ; do i=is,ie
3019  work_2d(i,j) = 0.0
3020  do k=1,nz
3021  work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
3022  enddo
3023  enddo ; enddo
3024  call post_data(cs%id_boundary_forcing_heat_tend_2d, work_2d, cs%diag)
3025  endif
3026  endif
3027 
3028  ! salinity tendency
3029  if (cs%id_boundary_forcing_saln_tend > 0) then
3030  do k=1,nz ; do j=js,je ; do i=is,ie
3031  work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*idt
3032  enddo ; enddo ; enddo
3033  call post_data(cs%id_boundary_forcing_saln_tend, work_3d, cs%diag, alt_h = h_old)
3034  endif
3035 
3036  ! salt tendency
3037  if (cs%id_boundary_forcing_salt_tend > 0 .or. cs%id_boundary_forcing_salt_tend_2d > 0) then
3038  do k=1,nz ; do j=js,je ; do i=is,ie
3039  work_3d(i,j,k) = gv%H_to_RZ * ppt2mks * idt * (h(i,j,k) * tv%S(i,j,k) - h_old(i,j,k) * saln_old(i,j,k))
3040  enddo ; enddo ; enddo
3041  if (cs%id_boundary_forcing_salt_tend > 0) then
3042  call post_data(cs%id_boundary_forcing_salt_tend, work_3d, cs%diag, alt_h = h_old)
3043  endif
3044  if (cs%id_boundary_forcing_salt_tend_2d > 0) then
3045  do j=js,je ; do i=is,ie
3046  work_2d(i,j) = 0.0
3047  do k=1,nz
3048  work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
3049  enddo
3050  enddo ; enddo
3051  call post_data(cs%id_boundary_forcing_salt_tend_2d, work_2d, cs%diag)
3052  endif
3053  endif
3054 

◆ diagnose_diabatic_diff_tendency()

subroutine mom_diabatic_driver::diagnose_diabatic_diff_tendency ( type(thermo_var_ptrs), intent(in)  tv,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  temp_old,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  saln_old,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diabatic_cs), pointer  CS 
)
private

This routine diagnoses tendencies from application of diabatic diffusion using ALE algorithm. Note that layer thickness is not altered by diabatic diffusion.

Parameters
[in]gocean grid structure
[in]gvocean vertical grid structure
[in]tvpoints to updated thermodynamic fields
[in]hthickness [H ~> m or kg m-2]
[in]temp_oldtemperature prior to diabatic physics
[in]saln_oldsalinity prior to diabatic physics [ppt]
[in]dttime step [T ~> s]
[in]usA dimensional unit scaling type
csmodule control structure

Definition at line 2873 of file MOM_diabatic_driver.F90.

2873  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2874  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2875  type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields
2876  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< thickness [H ~> m or kg m-2]
2877  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: temp_old !< temperature prior to diabatic physics
2878  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: saln_old !< salinity prior to diabatic physics [ppt]
2879  real, intent(in) :: dt !< time step [T ~> s]
2880  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2881  type(diabatic_CS), pointer :: CS !< module control structure
2882 
2883  ! Local variables
2884  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
2885  real, dimension(SZI_(G),SZJ_(G)) :: work_2d
2886  real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
2887  real :: ppt2mks = 0.001 ! Conversion factor from g/kg to kg/kg.
2888  integer :: i, j, k, is, ie, js, je, nz
2889  logical :: do_saln_tend ! Calculate salinity-based tendency diagnosics
2890 
2891  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2892  idt = 0.0 ; if (dt > 0.0) idt = 1. / dt
2893  work_3d(:,:,:) = 0.0
2894  work_2d(:,:) = 0.0
2895 
2896 
2897  ! temperature tendency
2898  do k=1,nz ; do j=js,je ; do i=is,ie
2899  work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt
2900  enddo ; enddo ; enddo
2901  if (cs%id_diabatic_diff_temp_tend > 0) then
2902  call post_data(cs%id_diabatic_diff_temp_tend, work_3d, cs%diag, alt_h=h)
2903  endif
2904 
2905  ! heat tendency
2906  if (cs%id_diabatic_diff_heat_tend > 0 .or. cs%id_diabatic_diff_heat_tend_2d > 0) then
2907  do k=1,nz ; do j=js,je ; do i=is,ie
2908  work_3d(i,j,k) = h(i,j,k)*gv%H_to_RZ * tv%C_p * work_3d(i,j,k)
2909  enddo ; enddo ; enddo
2910  if (cs%id_diabatic_diff_heat_tend > 0) then
2911  call post_data(cs%id_diabatic_diff_heat_tend, work_3d, cs%diag, alt_h=h)
2912  endif
2913  if (cs%id_diabatic_diff_heat_tend_2d > 0) then
2914  do j=js,je ; do i=is,ie
2915  work_2d(i,j) = 0.0
2916  do k=1,nz
2917  work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
2918  enddo
2919  enddo ; enddo
2920  call post_data(cs%id_diabatic_diff_heat_tend_2d, work_2d, cs%diag)
2921  endif
2922  endif
2923 
2924  ! salinity tendency
2925  do_saln_tend = cs%id_diabatic_diff_saln_tend > 0 &
2926  .or. cs%id_diabatic_diff_salt_tend > 0 &
2927  .or. cs%id_diabatic_diff_salt_tend_2d > 0
2928 
2929  if (do_saln_tend) then
2930  do k=1,nz ; do j=js,je ; do i=is,ie
2931  work_3d(i,j,k) = (tv%S(i,j,k) - saln_old(i,j,k)) * idt
2932  enddo ; enddo ; enddo
2933 
2934  if (cs%id_diabatic_diff_saln_tend > 0) &
2935  call post_data(cs%id_diabatic_diff_saln_tend, work_3d, cs%diag, alt_h=h)
2936 
2937  ! salt tendency
2938  if (cs%id_diabatic_diff_salt_tend > 0 .or. cs%id_diabatic_diff_salt_tend_2d > 0) then
2939  do k=1,nz ; do j=js,je ; do i=is,ie
2940  work_3d(i,j,k) = h(i,j,k)*gv%H_to_RZ * ppt2mks * work_3d(i,j,k)
2941  enddo ; enddo ; enddo
2942  if (cs%id_diabatic_diff_salt_tend > 0) then
2943  call post_data(cs%id_diabatic_diff_salt_tend, work_3d, cs%diag, alt_h=h)
2944  endif
2945  if (cs%id_diabatic_diff_salt_tend_2d > 0) then
2946  do j=js,je ; do i=is,ie
2947  work_2d(i,j) = 0.0
2948  do k=1,nz
2949  work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
2950  enddo
2951  enddo ; enddo
2952  call post_data(cs%id_diabatic_diff_salt_tend_2d, work_2d, cs%diag)
2953  endif
2954  endif
2955  endif
2956 

◆ diagnose_frazil_tendency()

subroutine mom_diabatic_driver::diagnose_frazil_tendency ( type(thermo_var_ptrs), intent(in)  tv,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  temp_old,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diabatic_cs), pointer  CS 
)
private

This routine diagnoses tendencies for temperature and heat from frazil formation. This routine is called twice from within subroutine diabatic; at start and at end of the diabatic processes. The impacts from frazil are generally a function of depth. Hence, when checking heat budget, be sure to remove HFSIFRAZIL from HFDS in k=1.

Parameters
[in]gocean grid structure
[in]gvocean vertical grid structure
[in]tvpoints to updated thermodynamic fields
[in]hthickness [H ~> m or kg m-2]
[in]temp_oldtemperature prior to frazil formation [degC]
[in]dttime step [T ~> s]
[in]usA dimensional unit scaling type
csmodule control structure

Definition at line 3063 of file MOM_diabatic_driver.F90.

3063  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
3064  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
3065  type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields
3066  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< thickness [H ~> m or kg m-2]
3067  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: temp_old !< temperature prior to frazil formation [degC]
3068  real, intent(in) :: dt !< time step [T ~> s]
3069  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3070  type(diabatic_CS), pointer :: CS !< module control structure
3071 
3072  real, dimension(SZI_(G),SZJ_(G)) :: work_2d
3073  real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
3074  integer :: i, j, k, is, ie, js, je, nz
3075 
3076  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3077  idt = 0.0 ; if (dt > 0.0) idt = 1. / dt
3078 
3079  ! temperature tendency
3080  if (cs%id_frazil_temp_tend > 0) then
3081  do k=1,nz ; do j=js,je ; do i=is,ie
3082  cs%frazil_temp_diag(i,j,k) = idt * (tv%T(i,j,k)-temp_old(i,j,k))
3083  enddo ; enddo ; enddo
3084  call post_data(cs%id_frazil_temp_tend, cs%frazil_temp_diag(:,:,:), cs%diag)
3085  endif
3086 
3087  ! heat tendency
3088  if (cs%id_frazil_heat_tend > 0 .or. cs%id_frazil_heat_tend_2d > 0) then
3089  do k=1,nz ; do j=js,je ; do i=is,ie
3090  cs%frazil_heat_diag(i,j,k) = gv%H_to_RZ * tv%C_p * h(i,j,k) * idt * (tv%T(i,j,k)-temp_old(i,j,k))
3091  enddo ; enddo ; enddo
3092  if (cs%id_frazil_heat_tend > 0) call post_data(cs%id_frazil_heat_tend, cs%frazil_heat_diag(:,:,:), cs%diag)
3093 
3094  ! As a consistency check, we must have
3095  ! FRAZIL_HEAT_TENDENCY_2d = HFSIFRAZIL
3096  if (cs%id_frazil_heat_tend_2d > 0) then
3097  do j=js,je ; do i=is,ie
3098  work_2d(i,j) = 0.0
3099  do k=1,nz
3100  work_2d(i,j) = work_2d(i,j) + cs%frazil_heat_diag(i,j,k)
3101  enddo
3102  enddo ; enddo
3103  call post_data(cs%id_frazil_heat_tend_2d, work_2d, cs%diag)
3104  endif
3105  endif
3106 

◆ extract_diabatic_member()

subroutine, public mom_diabatic_driver::extract_diabatic_member ( type(diabatic_cs), intent(in)  CS,
type(opacity_cs), optional, pointer  opacity_CSp,
type(optics_type), optional, pointer  optics_CSp,
real, intent(out), optional  evap_CFL_limit,
real, intent(out), optional  minimum_forcing_depth,
type(kpp_cs), optional, pointer  KPP_CSp,
type(energetic_pbl_cs), optional, pointer  energetic_PBL_CSp,
type(diabatic_aux_cs), optional, pointer  diabatic_aux_CSp,
integer, intent(out), optional  diabatic_halo 
)

Returns pointers or values of members within the diabatic_CS type. For extensibility, each returned argument is an optional argument.

Parameters
[in]csmodule control structure
opacity_cspA pointer to be set to the opacity control structure
optics_cspA pointer to be set to the optics control structure
kpp_cspA pointer to be set to the KPP CS
energetic_pbl_cspA pointer to be set to the ePBL CS
[out]evap_cfl_limitThe largest fraction of a layer that can be evaporated in one time-step [nondim].
[out]minimum_forcing_depthThe smallest depth over which heat and freshwater fluxes are applied [H ~> m or kg m-2].
diabatic_aux_cspA pointer to be set to the diabatic_aux control structure
[out]diabatic_haloThe halo size where the diabatic algorithms assume thermodynamics properties are valid.

Definition at line 2819 of file MOM_diabatic_driver.F90.

2819  type(diabatic_CS), intent(in ) :: CS !< module control structure
2820  ! All output arguments are optional
2821  type(opacity_CS), optional, pointer :: opacity_CSp !< A pointer to be set to the opacity control structure
2822  type(optics_type), optional, pointer :: optics_CSp !< A pointer to be set to the optics control structure
2823  type(KPP_CS), optional, pointer :: KPP_CSp !< A pointer to be set to the KPP CS
2824  type(energetic_PBL_CS), optional, pointer :: energetic_PBL_CSp !< A pointer to be set to the ePBL CS
2825  real, optional, intent( out) :: evap_CFL_limit !<The largest fraction of a layer that can be
2826  !! evaporated in one time-step [nondim].
2827  real, optional, intent( out) :: minimum_forcing_depth !< The smallest depth over which heat
2828  !! and freshwater fluxes are applied [H ~> m or kg m-2].
2829  type(diabatic_aux_CS), optional, pointer :: diabatic_aux_CSp !< A pointer to be set to the diabatic_aux
2830  !! control structure
2831  integer, optional, intent( out) :: diabatic_halo !< The halo size where the diabatic algorithms
2832  !! assume thermodynamics properties are valid.
2833 
2834  ! Pointers to control structures
2835  if (present(opacity_csp)) opacity_csp => cs%opacity_CSp
2836  if (present(optics_csp)) optics_csp => cs%optics
2837  if (present(kpp_csp)) kpp_csp => cs%KPP_CSp
2838  if (present(energetic_pbl_csp)) energetic_pbl_csp => cs%energetic_PBL_CSp
2839 
2840  ! Constants within diabatic_CS
2841  if (present(evap_cfl_limit)) evap_cfl_limit = cs%evap_CFL_limit
2842  if (present(minimum_forcing_depth)) minimum_forcing_depth = cs%minimum_forcing_depth
2843  if (present(diabatic_halo)) diabatic_halo = cs%halo_TS_diff
2844 

◆ layered_diabatic()

subroutine mom_diabatic_driver::layered_diabatic ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(inout)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h,
type(thermo_var_ptrs), intent(inout)  tv,
real, dimension(:,:), pointer  Hml,
type(forcing), intent(inout)  fluxes,
type(vertvisc_type), intent(inout)  visc,
type(accel_diag_ptrs), intent(inout)  ADp,
type(cont_diag_ptrs), intent(inout)  CDp,
real, intent(in)  dt,
type(time_type), intent(in)  Time_end,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diabatic_cs), pointer  CS,
type(wave_parameters_cs), optional, pointer  WAVES 
)
private

Imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers using the original MOM6 algorithms.

Parameters
[in,out]gocean grid structure
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]uzonal velocity [L T-1 ~> m s-1]
[in,out]vmeridional velocity [L T-1 ~> m s-1]
[in,out]hthickness [H ~> m or kg m-2]
[in,out]tvpoints to thermodynamic fields unused have NULL ptrs
hmlActive mixed layer depth [Z ~> m]
[in,out]fluxespoints to forcing fields unused fields have NULL ptrs
[in,out]viscvertical viscosities, BBL properies, and
[in,out]adprelated points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]cdppoints to terms in continuity equations
[in]dttime increment [T ~> s]
[in]time_endTime at the end of the interval
csmodule control structure
wavesSurface gravity waves

Definition at line 1855 of file MOM_diabatic_driver.F90.

1855  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1856  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1857  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1858  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1859  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
1860  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2]
1861  type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields
1862  !! unused have NULL ptrs
1863  real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m]
1864  type(forcing), intent(inout) :: fluxes !< points to forcing fields
1865  !! unused fields have NULL ptrs
1866  type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and
1867  type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum
1868  !! equations, to enable the later derived
1869  !! diagnostics, like energy budgets
1870  type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations
1871  real, intent(in) :: dt !< time increment [T ~> s]
1872  type(time_type), intent(in) :: Time_end !< Time at the end of the interval
1873  type(diabatic_CS), pointer :: CS !< module control structure
1874  type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves
1875 
1876  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1877  ea, & ! amount of fluid entrained from the layer above within
1878  ! one time step [H ~> m or kg m-2]
1879  eb, & ! amount of fluid entrained from the layer below within
1880  ! one time step [H ~> m or kg m-2]
1881  kd_lay, & ! diapycnal diffusivity of layers [Z2 T-1 ~> m2 s-1]
1882  h_orig, & ! initial layer thicknesses [H ~> m or kg m-2]
1883  h_prebound, & ! initial layer thicknesses [H ~> m or kg m-2]
1884  hold, & ! layer thickness before diapycnal entrainment, and later
1885  ! the initial layer thicknesses (if a mixed layer is used),
1886  ! [H ~> m or kg m-2]
1887  dsv_dt, & ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
1888  dsv_ds, & ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
1889  u_h, & ! zonal and meridional velocities at thickness points after
1890  v_h ! entrainment [L T-1 ~> m s-1]
1891  real, dimension(SZI_(G),SZJ_(G)) :: &
1892  Rcv_ml, & ! coordinate density of mixed layer, used for applying sponges
1893  SkinBuoyFlux! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
1894  real, dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag ! diagnostic array for thickness
1895  real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag ! diagnostic array for temp
1896  real, dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag ! diagnostic array for salinity
1897  real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn
1898 
1899  real :: net_ent ! The net of ea-eb at an interface.
1900 
1901  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
1902  ! These are targets so that the space can be shared with eaml & ebml.
1903  eatr, & ! The equivalent of ea and eb for tracers, which differ from ea and
1904  ebtr ! eb in that they tend to homogenize tracers in massless layers
1905  ! near the boundaries [H ~> m or kg m-2] (for Bous or non-Bouss)
1906 
1907  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
1908  Kd_int, & ! diapycnal diffusivity of interfaces [Z2 T-1 ~> m2 s-1]
1909  Kd_heat, & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1]
1910  Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1]
1911  Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to
1912  ! Kd_int [Z2 T-1 ~> m2 s-1].
1913  kd_extra_s , & ! The extra diffusivity of salinity due to double diffusion relative to
1914  ! Kd_int [Z2 T-1 ~> m2 s-1].
1915  kd_epbl, & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]
1916  tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1917  tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1918  sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1919  sadv_flx ! advective diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1920 
1921  real, allocatable, dimension(:,:) :: &
1922  hf_dudt_dia_2d, hf_dvdt_dia_2d ! Depth sum of diapycnal mixing accelaration * fract. thickness [L T-2 ~> m s-2].
1923 
1924  ! The following 5 variables are only used with a bulk mixed layer.
1925  real, pointer, dimension(:,:,:) :: &
1926  eaml, & ! The equivalent of ea due to mixed layer processes [H ~> m or kg m-2].
1927  ebml ! The equivalent of eb due to mixed layer processes [H ~> m or kg m-2].
1928  ! eaml and ebml are pointers to eatr and ebtr so as to reuse the memory as
1929  ! the arrays are not needed at the same time.
1930 
1931  integer :: kb(SZI_(G),SZJ_(G)) ! index of the lightest layer denser
1932  ! than the buffer layer [nondim]
1933 
1934  real :: p_ref_cv(SZI_(G)) ! Reference pressure for the potential density that defines the
1935  ! coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].
1936 
1937  logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below,
1938  ! where massive is defined as sufficiently thick that
1939  ! the no-flux boundary conditions have not restricted
1940  ! the entrainment - usually sqrt(Kd*dt).
1941 
1942  real :: b_denom_1 ! The first term in the denominator of b1
1943  ! [H ~> m or kg m-2]
1944  real :: h_neglect ! A thickness that is so small it is usually lost
1945  ! in roundoff and can be neglected
1946  ! [H ~> m or kg m-2]
1947  real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4]
1948  real :: add_ent ! Entrainment that needs to be added when mixing tracers
1949  ! [H ~> m or kg m-2]
1950  real :: eaval ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
1951  real :: hval ! hval is 2*h at velocity grid points [H ~> m or kg m-2]
1952  real :: h_tr ! h_tr is h at tracer points with a tiny thickness
1953  ! added to ensure positive definiteness [H ~> m or kg m-2]
1954  real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is
1955  ! coupled to the bottom within a timestep [H ~> m or kg m-2]
1956 
1957  real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].
1958  real :: b1(SZIB_(G)), d1(SZIB_(G)) ! b1, c1, and d1 are variables used by the
1959  real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver.
1960 
1961  real :: dt_mix ! The amount of time over which to apply mixing [T ~> s]
1962  real :: Idt ! The inverse time step [T-1 ~> s-1]
1963  real :: Idt_accel ! The inverse time step times rescaling factors [T-1 ~> s-1]
1964 
1965  integer :: dir_flag ! An integer encoding the directions in which to do halo updates.
1966  logical :: showCallTree ! If true, show the call tree
1967  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1968  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo
1969 
1970  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1971  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1972  nkmb = gv%nk_rho_varies
1973  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
1974  kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0
1975 
1976 
1977  showcalltree = calltree_showquery()
1978  if (showcalltree) call calltree_enter("layered_diabatic(), MOM_diabatic_driver.F90")
1979 
1980  ! set equivalence between the same bits of memory for these arrays
1981  eaml => eatr ; ebml => ebtr
1982 
1983  ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
1984  call enable_averages(dt, time_end, cs%diag)
1985 
1986  if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
1987  halo = cs%halo_TS_diff
1988  !$OMP parallel do default(shared)
1989  do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo
1990  h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0
1991  enddo ; enddo ; enddo
1992  endif
1993 
1994  if (cs%use_geothermal) then
1995  call cpu_clock_begin(id_clock_geothermal)
1996  call geothermal_entraining(h, tv, dt, eaml, ebml, g, gv, us, cs%geothermal_CSp, halo=cs%halo_TS_diff)
1997  call cpu_clock_end(id_clock_geothermal)
1998  if (showcalltree) call calltree_waypoint("geothermal (diabatic)")
1999  if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)
2000  endif
2001 
2002  ! Whenever thickness changes let the diag manager know, target grids
2003  ! for vertical remapping may need to be regenerated.
2004  call diag_update_remap_grids(cs%diag)
2005 
2006  ! Set_pen_shortwave estimates the optical properties of the water column.
2007  ! It will need to be modified later to include information about the
2008  ! biological properties and layer thicknesses.
2009  if (associated(cs%optics)) &
2010  call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
2011 
2012  if (cs%bulkmixedlayer) then
2013  if (cs%debug) call mom_forcing_chksum("Before mixedlayer", fluxes, g, us, haloshift=0)
2014 
2015  if (cs%ML_mix_first > 0.0) then
2016 ! This subroutine
2017 ! (1) Cools the mixed layer.
2018 ! (2) Performs convective adjustment by mixed layer entrainment.
2019 ! (3) Heats the mixed layer and causes it to detrain to
2020 ! Monin-Obukhov depth or minimum mixed layer depth.
2021 ! (4) Uses any remaining TKE to drive mixed layer entrainment.
2022 ! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate)
2023  call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
2024 
2025  call cpu_clock_begin(id_clock_mixedlayer)
2026  if (cs%ML_mix_first < 1.0) then
2027  ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???)
2028  call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*cs%ML_mix_first, &
2029  eaml,ebml, g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
2030  hml, cs%aggregate_FW_forcing, dt, last_call=.false.)
2031  ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???)
2032  call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, &
2033  g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
2034  hml, cs%aggregate_FW_forcing, dt, last_call=.true.)
2035  endif
2036 
2037  ! Keep salinity from falling below a small but positive threshold.
2038  ! This constraint is needed for SIS1 ice model, which can extract
2039  ! more salt than is present in the ocean. SIS2 does not suffer
2040  ! from this limitation, in which case we can let salinity=0 and still
2041  ! have salt conserved with SIS2 ice. So for SIS2, we can run with
2042  ! BOUND_SALINITY=False in MOM.F90.
2043  if (associated(tv%S) .and. associated(tv%salt_deficit)) &
2044  call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
2045  call cpu_clock_end(id_clock_mixedlayer)
2046  if (cs%debug) then
2047  call mom_state_chksum("After mixedlayer ", u, v, h, g, gv, us, haloshift=0)
2048  call mom_forcing_chksum("After mixedlayer", fluxes, g, us, haloshift=0)
2049  endif
2050  if (showcalltree) call calltree_waypoint("done with 1st bulkmixedlayer (diabatic)")
2051  if (cs%debugConservation) call mom_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, g, gv, us)
2052  endif
2053  endif
2054 
2055  if (cs%debug) &
2056  call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)
2057  if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then
2058  if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
2059  call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, us, eaml, ebml)
2060  if (cs%debug) then
2061  call hchksum(eaml, "after find_uv_at_h eaml", g%HI, scale=gv%H_to_m)
2062  call hchksum(ebml, "after find_uv_at_h ebml", g%HI, scale=gv%H_to_m)
2063  endif
2064  else
2065  call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
2066  endif
2067  if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
2068  endif
2069 
2070  call cpu_clock_begin(id_clock_set_diffusivity)
2071  ! Sets: Kd_lay, Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb
2072  ! Also changes: visc%Kd_shear and visc%Kv_shear
2073  if ((cs%halo_TS_diff > 0) .and. (cs%ML_mix_first > 0.0)) then
2074  if (associated(tv%T)) call pass_var(tv%T, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
2075  if (associated(tv%T)) call pass_var(tv%S, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
2076  call pass_var(h, g%domain, halo=cs%halo_TS_diff, complete=.true.)
2077  endif
2078  if (cs%debug) &
2079  call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)
2080  if (cs%double_diffuse) then
2081  call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, cs%set_diff_CSp, &
2082  kd_lay=kd_lay, kd_int=kd_int, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)
2083  else
2084  call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, &
2085  cs%set_diff_CSp, kd_lay=kd_lay, kd_int=kd_int)
2086  endif
2087  call cpu_clock_end(id_clock_set_diffusivity)
2088  if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")
2089 
2090  if (cs%debug) then
2091  call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)
2092  call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)
2093  call mom_thermovar_chksum("after set_diffusivity ", tv, g)
2094  call hchksum(kd_lay, "after set_diffusivity Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2095  call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2096  endif
2097 
2098 
2099  if (cs%useKPP) then
2100  call cpu_clock_begin(id_clock_kpp)
2101  ! KPP needs the surface buoyancy flux but does not update state variables.
2102  ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
2103  ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux
2104  ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second)
2105  ! unlike other instances where the fluxes are integrated in time over a time-step.
2106  call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
2107  cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
2108  ! The KPP scheme calculates boundary layer diffusivities and non-local transport.
2109 
2110  ! Set diffusivities for heat and salt separately
2111 
2112  if (cs%double_diffuse) then
2113  ! Add contribution from double diffusion
2114  !$OMP parallel do default(shared)
2115  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2116  kd_salt(i,j,k) = kd_int(i,j,k) + kd_extra_s(i,j,k)
2117  kd_heat(i,j,k) = kd_int(i,j,k) + kd_extra_t(i,j,k)
2118  enddo ; enddo ; enddo
2119  else
2120  !$OMP parallel do default(shared)
2121  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2122  kd_salt(i,j,k) = kd_int(i,j,k)
2123  kd_heat(i,j,k) = kd_int(i,j,k)
2124  enddo ; enddo ; enddo
2125  endif
2126 
2127  call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &
2128  fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
2129 
2130  call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
2131  kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
2132 
2133  if (associated(hml)) then
2134  call kpp_get_bld(cs%KPP_CSp, hml(:,:), g, us)
2135  call pass_var(hml, g%domain, halo=1)
2136  ! If visc%MLD exists, copy KPP's BLD into it
2137  if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
2138  endif
2139 
2140  if (.not. cs%KPPisPassive) then
2141  !$OMP parallel do default(shared)
2142  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2143  kd_int(i,j,k) = min( kd_salt(i,j,k), kd_heat(i,j,k) )
2144  enddo ; enddo ; enddo
2145  if (cs%double_diffuse) then
2146  !$OMP parallel do default(shared)
2147  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2148  kd_extra_s(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))
2149  kd_extra_t(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))
2150  enddo ; enddo ; enddo
2151  endif
2152  endif ! not passive
2153 
2154  call cpu_clock_end(id_clock_kpp)
2155  if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")
2156  if (cs%debug) then
2157  call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)
2158  call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)
2159  call mom_thermovar_chksum("after KPP", tv, g)
2160  call hchksum(kd_lay, "after KPP Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2161  call hchksum(kd_int, "after KPP Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2162  endif
2163 
2164  endif ! endif for KPP
2165 
2166  ! Add vertical diff./visc. due to convection (computed via CVMix)
2167  if (cs%use_CVMix_conv) then
2168  call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml, kd=kd_int, kv=visc%Kv_slow)
2169  endif
2170 
2171  if (cs%useKPP) then
2172  call cpu_clock_begin(id_clock_kpp)
2173  if (cs%debug) then
2174  call hchksum(cs%KPP_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
2175  call hchksum(cs%KPP_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
2176  call hchksum(cs%KPP_NLTheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)
2177  call hchksum(cs%KPP_NLTscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
2178  endif
2179  ! Apply non-local transport of heat and salt
2180  ! Changes: tv%T, tv%S
2181  call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat, cs%KPP_temp_flux, &
2182  us%T_to_s*dt, tv%T, us%Q_to_J_kg*tv%C_p)
2183  call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, &
2184  us%T_to_s*dt, tv%S)
2185  call cpu_clock_end(id_clock_kpp)
2186  if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
2187  if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)
2188 
2189  if (cs%debug) then
2190  call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)
2191  call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)
2192  call mom_thermovar_chksum("after KPP_applyNLT ", tv, g)
2193  endif
2194  endif ! endif for KPP
2195 
2196  ! Differential diffusion done here.
2197  ! Changes: tv%T, tv%S
2198  if (cs%double_diffuse .and. associated(tv%T)) then
2199 
2200  call cpu_clock_begin(id_clock_differential_diff)
2201  call differential_diffuse_t_s(h, tv%T, tv%S, kd_extra_t, kd_extra_s, dt, g, gv)
2202  call cpu_clock_end(id_clock_differential_diff)
2203  if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")
2204  if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g, gv, us)
2205 
2206  ! increment heat and salt diffusivity.
2207  ! CS%useKPP==.true. already has extra_T and extra_S included
2208  if (.not. cs%useKPP) then
2209  !$OMP parallel do default(shared)
2210  do k=2,nz ; do j=js,je ; do i=is,ie
2211  kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)
2212  kd_salt(i,j,k) = kd_salt(i,j,k) + kd_extra_s(i,j,k)
2213  enddo ; enddo ; enddo
2214  endif
2215 
2216  endif
2217 
2218  ! Calculate layer entrainments and detrainments from diffusivities and differences between
2219  ! layer and target densities (i.e. do remapping as well as diffusion).
2220  call cpu_clock_begin(id_clock_entrain)
2221  ! Calculate appropriately limited diapycnal mass fluxes to account
2222  ! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb
2223  call entrainment_diffusive(h, tv, fluxes, dt, g, gv, us, cs%entrain_diffusive_CSp, &
2224  ea, eb, kb, kd_lay=kd_lay, kd_int=kd_int)
2225  call cpu_clock_end(id_clock_entrain)
2226  if (showcalltree) call calltree_waypoint("done with Entrainment_diffusive (diabatic)")
2227 
2228  if (cs%debug) then
2229  call mom_forcing_chksum("after calc_entrain ", fluxes, g, us, haloshift=0)
2230  call mom_thermovar_chksum("after calc_entrain ", tv, g)
2231  call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, us, haloshift=0)
2232  call hchksum(ea, "after calc_entrain ea", g%HI, haloshift=0, scale=gv%H_to_m)
2233  call hchksum(eb, "after calc_entrain eb", g%HI, haloshift=0, scale=gv%H_to_m)
2234  endif
2235 
2236  ! Save fields before boundary forcing is applied for tendency diagnostics
2237  if (cs%boundary_forcing_tendency_diag) then
2238  do k=1,nz ; do j=js,je ; do i=is,ie
2239  h_diag(i,j,k) = h(i,j,k)
2240  temp_diag(i,j,k) = tv%T(i,j,k)
2241  saln_diag(i,j,k) = tv%S(i,j,k)
2242  enddo ; enddo ; enddo
2243  endif
2244 
2245  ! Update h according to divergence of the difference between
2246  ! ea and eb. We keep a record of the original h in hold.
2247  ! In the following, the checks for negative values are to guard
2248  ! against instances where entrainment drives a layer to
2249  ! negative thickness. This situation will never happen if
2250  ! enough iterations are permitted in Calculate_Entrainment.
2251  ! Even if too few iterations are allowed, it is still guarded
2252  ! against. In other words the checks are probably unnecessary.
2253  !$OMP parallel do default(shared)
2254  do j=js,je
2255  do i=is,ie
2256  hold(i,j,1) = h(i,j,1)
2257  h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2))
2258  hold(i,j,nz) = h(i,j,nz)
2259  h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1))
2260  if (h(i,j,1) <= 0.0) then
2261  h(i,j,1) = gv%Angstrom_H
2262  endif
2263  if (h(i,j,nz) <= 0.0) then
2264  h(i,j,nz) = gv%Angstrom_H
2265  endif
2266  enddo
2267  do k=2,nz-1 ; do i=is,ie
2268  hold(i,j,k) = h(i,j,k)
2269  h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
2270  (eb(i,j,k) - ea(i,j,k+1)))
2271  if (h(i,j,k) <= 0.0) then
2272  h(i,j,k) = gv%Angstrom_H
2273  endif
2274  enddo ; enddo
2275  enddo
2276  ! Checks for negative thickness may have changed layer thicknesses
2277  call diag_update_remap_grids(cs%diag)
2278 
2279  if (cs%debug) then
2280  call mom_state_chksum("after negative check ", u, v, h, g, gv, us, haloshift=0)
2281  call mom_forcing_chksum("after negative check ", fluxes, g, us, haloshift=0)
2282  call mom_thermovar_chksum("after negative check ", tv, g)
2283  endif
2284  if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")
2285  if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)
2286 
2287  ! Here, T and S are updated according to ea and eb.
2288  ! If using the bulk mixed layer, T and S are also updated
2289  ! by surface fluxes (in fluxes%*).
2290  ! This is a very long block.
2291  if (cs%bulkmixedlayer) then
2292 
2293  if (associated(tv%T)) then
2294  call cpu_clock_begin(id_clock_tridiag)
2295  ! Temperature and salinity (as state variables) are treated
2296  ! differently from other tracers to insure massless layers that
2297  ! are lighter than the mixed layer have temperatures and salinities
2298  ! that correspond to their prescribed densities.
2299  if (cs%massless_match_targets) then
2300  !$OMP parallel do default (shared) private(h_tr,b1,d1,c1,b_denom_1)
2301  do j=js,je
2302  do i=is,ie
2303  h_tr = hold(i,j,1) + h_neglect
2304  b1(i) = 1.0 / (h_tr + eb(i,j,1))
2305  d1(i) = h_tr * b1(i)
2306  tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1))
2307  tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1))
2308  enddo
2309  do k=2,nkmb ; do i=is,ie
2310  c1(i,k) = eb(i,j,k-1) * b1(i)
2311  h_tr = hold(i,j,k) + h_neglect
2312  b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2313  b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2314  if (k<nkmb) d1(i) = b_denom_1 * b1(i)
2315  tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
2316  tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
2317  enddo ; enddo
2318 
2319  do k=nkmb+1,nz ; do i=is,ie
2320  if (k == kb(i,j)) then
2321  c1(i,k) = eb(i,j,k-1) * b1(i)
2322  d1(i) = (((eb(i,j,nkmb)-eb(i,j,k-1)) + hold(i,j,nkmb) + h_neglect) + &
2323  d1(i)*ea(i,j,nkmb)) * b1(i)
2324  h_tr = hold(i,j,k) + h_neglect
2325  b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2326  b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2327  d1(i) = b_denom_1 * b1(i)
2328  tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,nkmb))
2329  tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,nkmb))
2330  elseif (k > kb(i,j)) then
2331  c1(i,k) = eb(i,j,k-1) * b1(i)
2332  h_tr = hold(i,j,k) + h_neglect
2333  b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2334  b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2335  d1(i) = b_denom_1 * b1(i)
2336  tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
2337  tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
2338  elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j))
2339  ! The bottommost buffer layer might entrain all the mass from some
2340  ! of the interior layers that are thin and lighter in the coordinate
2341  ! density than that buffer layer. The T and S of these newly
2342  ! massless interior layers are unchanged.
2343  tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k)
2344  tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k)
2345  endif
2346  enddo ; enddo
2347 
2348  do k=nz-1,nkmb,-1 ; do i=is,ie
2349  if (k >= kb(i,j)) then
2350  tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
2351  tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
2352  endif
2353  enddo ; enddo
2354  do i=is,ie ; if (kb(i,j) <= nz) then
2355  tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j))
2356  tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j))
2357  endif ; enddo
2358  do k=nkmb-1,1,-1 ; do i=is,ie
2359  tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
2360  tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
2361  enddo ; enddo
2362  enddo ! end of j loop
2363  else ! .not. massless_match_targets
2364  ! This simpler form allows T & S to be too dense for the layers
2365  ! between the buffer layers and the interior.
2366  ! Changes: T, S
2367  if (cs%tracer_tridiag) then
2368  call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
2369  call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
2370  else
2371  call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
2372  endif
2373  endif ! massless_match_targets
2374  call cpu_clock_end(id_clock_tridiag)
2375 
2376  endif ! endif for associated(T)
2377  if (cs%debugConservation) call mom_state_stats('BML tridiag', u, v, h, tv%T, tv%S, g, gv, us)
2378 
2379  if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
2380  ! The mixed layer code has already been called, but there is some needed
2381  ! bookkeeping.
2382  !$OMP parallel do default(shared)
2383  do k=1,nz ; do j=js,je ; do i=is,ie
2384  hold(i,j,k) = h_orig(i,j,k)
2385  ea(i,j,k) = ea(i,j,k) + eaml(i,j,k)
2386  eb(i,j,k) = eb(i,j,k) + ebml(i,j,k)
2387  enddo ; enddo ; enddo
2388  if (cs%debug) then
2389  call hchksum(ea, "after ea = ea + eaml", g%HI, haloshift=0, scale=gv%H_to_m)
2390  call hchksum(eb, "after eb = eb + ebml", g%HI, haloshift=0, scale=gv%H_to_m)
2391  endif
2392  endif
2393 
2394  if (cs%ML_mix_first < 1.0) then
2395  ! Call the mixed layer code now, perhaps for a second time.
2396  ! This subroutine (1) Cools the mixed layer.
2397  ! (2) Performs convective adjustment by mixed layer entrainment.
2398  ! (3) Heats the mixed layer and causes it to detrain to
2399  ! Monin-Obukhov depth or minimum mixed layer depth.
2400  ! (4) Uses any remaining TKE to drive mixed layer entrainment.
2401  ! (5) Possibly splits the buffer layer into two isopycnal layers.
2402 
2403  call find_uv_at_h(u, v, hold, u_h, v_h, g, gv, us, ea, eb)
2404  if (cs%debug) call mom_state_chksum("find_uv_at_h1 ", u, v, h, g, gv, us, haloshift=0)
2405 
2406  dt_mix = min(dt, dt*(1.0 - cs%ML_mix_first))
2407  call cpu_clock_begin(id_clock_mixedlayer)
2408  ! Changes: h, tv%T, tv%S, ea and eb (G is also inout???)
2409  call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, &
2410  g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
2411  hml, cs%aggregate_FW_forcing, dt, last_call=.true.)
2412 
2413  ! Keep salinity from falling below a small but positive threshold.
2414  ! This constraint is needed for SIS1 ice model, which can extract
2415  ! more salt than is present in the ocean. SIS2 does not suffer
2416  ! from this limitation, in which case we can let salinity=0 and still
2417  ! have salt conserved with SIS2 ice. So for SIS2, we can run with
2418  ! BOUND_SALINITY=False in MOM.F90.
2419  if (associated(tv%S) .and. associated(tv%salt_deficit)) &
2420  call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
2421 
2422  call cpu_clock_end(id_clock_mixedlayer)
2423  if (showcalltree) call calltree_waypoint("done with 2nd bulkmixedlayer (diabatic)")
2424  if (cs%debugConservation) call mom_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, g, gv, us)
2425  endif
2426 
2427  else ! following block for when NOT using BULKMIXEDLAYER
2428 
2429  ! calculate change in temperature & salinity due to dia-coordinate surface diffusion
2430  if (associated(tv%T)) then
2431 
2432  if (cs%debug) then
2433  call hchksum(ea, "before triDiagTS ea ", g%HI, haloshift=0, scale=gv%H_to_m)
2434  call hchksum(eb, "before triDiagTS eb ", g%HI, haloshift=0, scale=gv%H_to_m)
2435  endif
2436  call cpu_clock_begin(id_clock_tridiag)
2437 
2438  ! Keep salinity from falling below a small but positive threshold.
2439  ! This constraint is needed for SIS1 ice model, which can extract
2440  ! more salt than is present in the ocean. SIS2 does not suffer
2441  ! from this limitation, in which case we can let salinity=0 and still
2442  ! have salt conserved with SIS2 ice. So for SIS2, we can run with
2443  ! BOUND_SALINITY=False in MOM.F90.
2444  if (associated(tv%S) .and. associated(tv%salt_deficit)) &
2445  call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
2446 
2447  if (cs%diabatic_diff_tendency_diag) then
2448  do k=1,nz ; do j=js,je ; do i=is,ie
2449  temp_diag(i,j,k) = tv%T(i,j,k)
2450  saln_diag(i,j,k) = tv%S(i,j,k)
2451  enddo ; enddo ; enddo
2452  endif
2453 
2454  ! Changes T and S via the tridiagonal solver; no change to h
2455  if (cs%tracer_tridiag) then
2456  call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
2457  call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
2458  else
2459  call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
2460  endif
2461 
2462  ! diagnose temperature, salinity, heat, and salt tendencies
2463  ! Note: hold here refers to the thicknesses from before the dual-entraintment when using
2464  ! the bulk mixed layer scheme, so tendencies should be posted on hold.
2465  if (cs%diabatic_diff_tendency_diag) then
2466  call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, us, cs)
2467  if (cs%id_diabatic_diff_h > 0) call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h=hold)
2468  endif
2469 
2470  call cpu_clock_end(id_clock_tridiag)
2471  if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")
2472 
2473  endif ! endif corresponding to if (associated(tv%T))
2474  if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)
2475 
2476  endif ! endif for the BULKMIXEDLAYER block
2477 
2478  if (cs%debug) then
2479  call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)
2480  call mom_thermovar_chksum("after mixed layer ", tv, g)
2481  call hchksum(ea, "after mixed layer ea", g%HI, scale=gv%H_to_m)
2482  call hchksum(eb, "after mixed layer eb", g%HI, scale=gv%H_to_m)
2483  endif
2484 
2485  call cpu_clock_begin(id_clock_remap)
2486  call regularize_layers(h, tv, dt, ea, eb, g, gv, us, cs%regularize_layers_CSp)
2487  call cpu_clock_end(id_clock_remap)
2488  if (showcalltree) call calltree_waypoint("done with regularize_layers (diabatic)")
2489  if (cs%debugConservation) call mom_state_stats('regularize_layers', u, v, h, tv%T, tv%S, g, gv, us)
2490 
2491  ! Whenever thickness changes let the diag manager know, as the
2492  ! target grids for vertical remapping may need to be regenerated.
2493  if (associated(adp%du_dt_dia) .or. associated(adp%dv_dt_dia)) &
2494  ! Remapped d[uv]dt_dia require east/north halo updates of h
2495  call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
2496  call diag_update_remap_grids(cs%diag)
2497 
2498  ! diagnostics
2499  idt = 1.0 / dt
2500  if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0)) then
2501  do j=js,je ; do i=is,ie
2502  tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
2503  tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
2504  enddo ; enddo
2505  !$OMP parallel do default(shared)
2506  do k=2,nz ; do j=js,je ; do i=is,ie
2507  tdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
2508  (tv%T(i,j,k-1) - tv%T(i,j,k))
2509  tadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
2510  0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
2511  enddo ; enddo ; enddo
2512  endif
2513  if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0)) then
2514  do j=js,je ; do i=is,ie
2515  sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
2516  sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
2517  enddo ; enddo
2518  !$OMP parallel do default(shared)
2519  do k=2,nz ; do j=js,je ; do i=is,ie
2520  sdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
2521  (tv%S(i,j,k-1) - tv%S(i,j,k))
2522  sadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
2523  0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
2524  enddo ; enddo ; enddo
2525  endif
2526 
2527  ! mixing of passive tracers from massless boundary layers to interior
2528  call cpu_clock_begin(id_clock_tracers)
2529  if (cs%mix_boundary_tracers) then
2530  tr_ea_bbl = gv%Z_to_H * sqrt(dt*cs%Kd_BBL_tr)
2531  !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
2532  do j=js,je
2533  do i=is,ie
2534  ebtr(i,j,nz) = eb(i,j,nz)
2535  htot(i) = 0.0
2536  in_boundary(i) = (g%mask2dT(i,j) > 0.0)
2537  enddo
2538  do k=nz,2,-1 ; do i=is,ie
2539  if (in_boundary(i)) then
2540  htot(i) = htot(i) + h(i,j,k)
2541  ! If diapycnal mixing has been suppressed because this is a massless
2542  ! layer near the bottom, add some mixing of tracers between these
2543  ! layers. This flux is based on the harmonic mean of the two
2544  ! thicknesses, as this corresponds pretty closely (to within
2545  ! differences in the density jumps between layers) with what is done
2546  ! in the calculation of the fluxes in the first place. Kd_min_tr
2547  ! should be much less than the values that have been set in Kd_lay,
2548  ! perhaps a molecular diffusivity.
2549  add_ent = ((dt * cs%Kd_min_tr) * gv%Z_to_H**2) * &
2550  ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
2551  (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
2552  0.5*(ea(i,j,k) + eb(i,j,k-1))
2553  if (htot(i) < tr_ea_bbl) then
2554  add_ent = max(0.0, add_ent, &
2555  (tr_ea_bbl - htot(i)) - min(ea(i,j,k), eb(i,j,k-1)))
2556  elseif (add_ent < 0.0) then
2557  add_ent = 0.0 ; in_boundary(i) = .false.
2558  endif
2559 
2560  ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
2561  eatr(i,j,k) = ea(i,j,k) + add_ent
2562  else
2563  ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k)
2564  endif
2565  if (cs%double_diffuse) then ; if (kd_extra_s(i,j,k) > 0.0) then
2566  add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
2567  (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
2568  h_neglect)
2569  ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
2570  eatr(i,j,k) = eatr(i,j,k) + add_ent
2571  endif ; endif
2572  enddo ; enddo
2573  do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo
2574 
2575  enddo
2576 
2577  call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
2578  cs%optics, cs%tracer_flow_CSp, cs%debug)
2579 
2580  elseif (cs%double_diffuse) then ! extra diffusivity for passive tracers
2581 
2582  do j=js,je ; do i=is,ie
2583  ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1)
2584  enddo ; enddo
2585  !$OMP parallel do default(shared) private(add_ent)
2586  do k=nz,2,-1 ; do j=js,je ; do i=is,ie
2587  if (kd_extra_s(i,j,k) > 0.0) then
2588  add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
2589  (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
2590  h_neglect)
2591  else
2592  add_ent = 0.0
2593  endif
2594  ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
2595  eatr(i,j,k) = ea(i,j,k) + add_ent
2596  enddo ; enddo ; enddo
2597 
2598  call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
2599  cs%optics, cs%tracer_flow_CSp, cs%debug)
2600 
2601  else
2602  call call_tracer_column_fns(hold, h, ea, eb, fluxes, hml, dt, g, gv, us, tv, &
2603  cs%optics, cs%tracer_flow_CSp, cs%debug)
2604 
2605  endif ! (CS%mix_boundary_tracers)
2606 
2607  call cpu_clock_end(id_clock_tracers)
2608 
2609  ! sponges
2610  if (cs%use_sponge) then
2611  call cpu_clock_begin(id_clock_sponge)
2612  ! Layer mode sponge
2613  if (cs%bulkmixedlayer .and. associated(tv%eqn_of_state)) then
2614  do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo
2615  eosdom(:) = eos_domain(g%HI)
2616  !$OMP parallel do default(shared)
2617  do j=js,je
2618  call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, rcv_ml(:,j), &
2619  tv%eqn_of_state, eosdom)
2620  enddo
2621  call apply_sponge(h, dt, g, gv, us, ea, eb, cs%sponge_CSp, rcv_ml)
2622  else
2623  call apply_sponge(h, dt, g, gv, us, ea, eb, cs%sponge_CSp)
2624  endif
2625  call cpu_clock_end(id_clock_sponge)
2626  if (cs%debug) then
2627  call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)
2628  call mom_thermovar_chksum("apply_sponge ", tv, g)
2629  endif
2630  endif ! CS%use_sponge
2631 
2632 ! Save the diapycnal mass fluxes as a diagnostic field.
2633  if (associated(cdp%diapyc_vel)) then
2634  !$OMP parallel do default(shared)
2635  do j=js,je
2636  do k=2,nz ; do i=is,ie
2637  cdp%diapyc_vel(i,j,k) = us%s_to_T*idt * (ea(i,j,k) - eb(i,j,k-1))
2638  enddo ; enddo
2639  do i=is,ie
2640  cdp%diapyc_vel(i,j,1) = 0.0
2641  cdp%diapyc_vel(i,j,nz+1) = 0.0
2642  enddo
2643  enddo
2644  endif
2645 
2646 ! For momentum, it is only the net flux that homogenizes within
2647 ! the mixed layer. Vertical viscosity that is proportional to the
2648 ! mixed layer turbulence is applied elsewhere.
2649  if (cs%bulkmixedlayer) then
2650  if (cs%debug) then
2651  call hchksum(ea, "before net flux rearrangement ea", g%HI, scale=gv%H_to_m)
2652  call hchksum(eb, "before net flux rearrangement eb", g%HI, scale=gv%H_to_m)
2653  endif
2654  !$OMP parallel do default(shared) private(net_ent)
2655  do j=js,je
2656  do k=2,gv%nkml ; do i=is,ie
2657  net_ent = ea(i,j,k) - eb(i,j,k-1)
2658  ea(i,j,k) = max(net_ent, 0.0)
2659  eb(i,j,k-1) = max(-net_ent, 0.0)
2660  enddo ; enddo
2661  enddo
2662  if (cs%debug) then
2663  call hchksum(ea, "after net flux rearrangement ea", g%HI, scale=gv%H_to_m)
2664  call hchksum(eb, "after net flux rearrangement eb", g%HI, scale=gv%H_to_m)
2665  endif
2666  endif
2667 
2668 ! Initialize halo regions of ea, eb, and hold to default values.
2669  !$OMP parallel do default(shared)
2670  do k=1,nz
2671  do i=is-1,ie+1
2672  hold(i,js-1,k) = gv%Angstrom_H ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0
2673  hold(i,je+1,k) = gv%Angstrom_H ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0
2674  enddo
2675  do j=js,je
2676  hold(is-1,j,k) = gv%Angstrom_H ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0
2677  hold(ie+1,j,k) = gv%Angstrom_H ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0
2678  enddo
2679  enddo
2680 
2681  call cpu_clock_begin(id_clock_pass)
2682  if (g%symmetric) then ; dir_flag = to_all+omit_corners
2683  else ; dir_flag = to_west+to_south+omit_corners ; endif
2684  call create_group_pass(cs%pass_hold_eb_ea, hold, g%Domain, dir_flag, halo=1)
2685  call create_group_pass(cs%pass_hold_eb_ea, eb, g%Domain, dir_flag, halo=1)
2686  call create_group_pass(cs%pass_hold_eb_ea, ea, g%Domain, dir_flag, halo=1)
2687  call do_group_pass(cs%pass_hold_eb_ea, g%Domain)
2688  call cpu_clock_end(id_clock_pass)
2689 
2690  ! Use a tridiagonal solver to determine effect of the diapycnal
2691  ! advection on velocity field. It is assumed that water leaves
2692  ! or enters the ocean with the surface velocity.
2693  if (cs%debug) then
2694  call mom_state_chksum("before u/v tridiag ", u, v, h, g, gv, us, haloshift=0)
2695  call hchksum(ea, "before u/v tridiag ea", g%HI, scale=gv%H_to_m)
2696  call hchksum(eb, "before u/v tridiag eb", g%HI, scale=gv%H_to_m)
2697  call hchksum(hold, "before u/v tridiag hold", g%HI, scale=gv%H_to_m)
2698  endif
2699  call cpu_clock_begin(id_clock_tridiag)
2700  idt_accel = 1.0 / dt
2701  !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)
2702  do j=js,je
2703  do i=isq,ieq
2704  if (associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,1) = u(i,j,1)
2705  hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect
2706  b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))
2707  d1(i) = hval * b1(i)
2708  u(i,j,1) = b1(i) * (hval * u(i,j,1))
2709  enddo
2710  do k=2,nz ; do i=isq,ieq
2711  if (associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,k) = u(i,j,k)
2712  c1(i,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(i)
2713  eaval = ea(i,j,k) + ea(i+1,j,k)
2714  hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect
2715  b1(i) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(i)*eaval))
2716  d1(i) = (hval + d1(i)*eaval) * b1(i)
2717  u(i,j,k) = (hval*u(i,j,k) + eaval*u(i,j,k-1))*b1(i)
2718  enddo ; enddo
2719  do k=nz-1,1,-1 ; do i=isq,ieq
2720  u(i,j,k) = u(i,j,k) + c1(i,k+1)*u(i,j,k+1)
2721  if (associated(adp%du_dt_dia)) &
2722  adp%du_dt_dia(i,j,k) = (u(i,j,k) - adp%du_dt_dia(i,j,k)) * idt_accel
2723  enddo ; enddo
2724  if (associated(adp%du_dt_dia)) then
2725  do i=isq,ieq
2726  adp%du_dt_dia(i,j,nz) = (u(i,j,nz)-adp%du_dt_dia(i,j,nz)) * idt_accel
2727  enddo
2728  endif
2729  enddo
2730  if (cs%debug) then
2731  call mom_state_chksum("aft 1st loop tridiag ", u, v, h, g, gv, us, haloshift=0)
2732  endif
2733  !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)
2734  do j=jsq,jeq
2735  do i=is,ie
2736  if (associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,1) = v(i,j,1)
2737  hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect
2738  b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))
2739  d1(i) = hval * b1(i)
2740  v(i,j,1) = b1(i) * (hval * v(i,j,1))
2741  enddo
2742  do k=2,nz ; do i=is,ie
2743  if (associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,k) = v(i,j,k)
2744  c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)
2745  eaval = ea(i,j,k) + ea(i,j+1,k)
2746  hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect
2747  b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))
2748  d1(i) = (hval + d1(i)*eaval) * b1(i)
2749  v(i,j,k) = (hval*v(i,j,k) + eaval*v(i,j,k-1))*b1(i)
2750  enddo ; enddo
2751  do k=nz-1,1,-1 ; do i=is,ie
2752  v(i,j,k) = v(i,j,k) + c1(i,k+1)*v(i,j,k+1)
2753  if (associated(adp%dv_dt_dia)) &
2754  adp%dv_dt_dia(i,j,k) = (v(i,j,k) - adp%dv_dt_dia(i,j,k)) * idt_accel
2755  enddo ; enddo
2756  if (associated(adp%dv_dt_dia)) then
2757  do i=is,ie
2758  adp%dv_dt_dia(i,j,nz) = (v(i,j,nz)-adp%dv_dt_dia(i,j,nz)) * idt_accel
2759  enddo
2760  endif
2761  enddo
2762  call cpu_clock_end(id_clock_tridiag)
2763  if (cs%debug) then
2764  call mom_state_chksum("after u/v tridiag ", u, v, h, g, gv, us, haloshift=0)
2765  endif
2766 
2767  call disable_averaging(cs%diag)
2768  ! Diagnose the diapycnal diffusivities and other related quantities.
2769  call enable_averages(dt, time_end, cs%diag)
2770 
2771  if (cs%id_Kd_interface > 0) call post_data(cs%id_Kd_interface, kd_int, cs%diag)
2772  if (cs%id_Kd_heat > 0) call post_data(cs%id_Kd_heat, kd_heat, cs%diag)
2773  if (cs%id_Kd_salt > 0) call post_data(cs%id_Kd_salt, kd_salt, cs%diag)
2774  if (cs%id_Kd_ePBL > 0) call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)
2775 
2776  if (cs%id_ea > 0) call post_data(cs%id_ea, ea, cs%diag)
2777  if (cs%id_eb > 0) call post_data(cs%id_eb, eb, cs%diag)
2778 
2779  if (cs%id_dudt_dia > 0) call post_data(cs%id_dudt_dia, adp%du_dt_dia, cs%diag)
2780  if (cs%id_dvdt_dia > 0) call post_data(cs%id_dvdt_dia, adp%dv_dt_dia, cs%diag)
2781  if (cs%id_wd > 0) call post_data(cs%id_wd, cdp%diapyc_vel, cs%diag)
2782 
2783  if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)
2784  if (cs%id_Tadv > 0) call post_data(cs%id_Tadv, tadv_flx, cs%diag)
2785  if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)
2786  if (cs%id_Sadv > 0) call post_data(cs%id_Sadv, sadv_flx, cs%diag)
2787 
2788  !! Diagnostics for terms multiplied by fractional thicknesses
2789  if (cs%id_hf_dudt_dia_2d > 0) then
2790  allocate(hf_dudt_dia_2d(g%IsdB:g%IedB,g%jsd:g%jed))
2791  hf_dudt_dia_2d(:,:) = 0.0
2792  do k=1,nz ; do j=js,je ; do i=isq,ieq
2793  hf_dudt_dia_2d(i,j) = hf_dudt_dia_2d(i,j) + adp%du_dt_dia(i,j,k) * adp%diag_hfrac_u(i,j,k)
2794  enddo ; enddo ; enddo
2795  call post_data(cs%id_hf_dudt_dia_2d, hf_dudt_dia_2d, cs%diag)
2796  deallocate(hf_dudt_dia_2d)
2797  endif
2798 
2799  if (cs%id_hf_dvdt_dia_2d > 0) then
2800  allocate(hf_dvdt_dia_2d(g%isd:g%ied,g%JsdB:g%JedB))
2801  hf_dvdt_dia_2d(:,:) = 0.0
2802  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
2803  hf_dvdt_dia_2d(i,j) = hf_dvdt_dia_2d(i,j) + adp%dv_dt_dia(i,j,k) * adp%diag_hfrac_v(i,j,k)
2804  enddo ; enddo ; enddo
2805  call post_data(cs%id_hf_dvdt_dia_2d, hf_dvdt_dia_2d, cs%diag)
2806  deallocate(hf_dvdt_dia_2d)
2807  endif
2808 
2809  call disable_averaging(cs%diag)
2810 
2811  if (showcalltree) call calltree_leave("layered_diabatic()")
2812