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 2856 of file MOM_diabatic_driver.F90.

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

◆ 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 3122 of file MOM_diabatic_driver.F90.

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

◆ 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 262 of file MOM_diabatic_driver.F90.

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

◆ 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 1175 of file MOM_diabatic_driver.F90.

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

◆ 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 452 of file MOM_diabatic_driver.F90.

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

◆ 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 3750 of file MOM_diabatic_driver.F90.

3750  type(diabatic_cs), pointer :: cs !< module control structure
3751 
3752  if (.not.associated(cs)) return
3753 
3754  call diabatic_aux_end(cs%diabatic_aux_CSp)
3755 
3756  call entrain_diffusive_end(cs%entrain_diffusive_CSp)
3757  call set_diffusivity_end(cs%set_diff_CSp)
3758 
3759  if (cs%useKPP) then
3760  deallocate( cs%KPP_buoy_flux )
3761  deallocate( cs%KPP_temp_flux )
3762  deallocate( cs%KPP_salt_flux )
3763  deallocate( cs%KPP_NLTheat )
3764  deallocate( cs%KPP_NLTscalar )
3765  call kpp_end(cs%KPP_CSp)
3766  endif
3767 
3768  if (cs%use_CVMix_conv) call cvmix_conv_end(cs%CVMix_conv_csp)
3769 
3770  if (cs%use_energetic_PBL) &
3771  call energetic_pbl_end(cs%energetic_PBL_CSp)
3772  if (cs%debug_energy_req) &
3773  call diapyc_energy_req_end(cs%diapyc_en_rec_CSp)
3774 
3775  if (associated(cs%optics)) then
3776  call opacity_end(cs%opacity_CSp, cs%optics)
3777  deallocate(cs%optics)
3778  endif
3779 
3780  ! GMM, the following is commented out because arrays in
3781  ! CS%diag_grids_prev are neither pointers or allocatables
3782  ! and, therefore, cannot be deallocated.
3783 
3784  !call diag_grid_storage_end(CS%diag_grids_prev)
3785 
3786  deallocate(cs)
3787 

◆ 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 3154 of file MOM_diabatic_driver.F90.

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

◆ 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 2973 of file MOM_diabatic_driver.F90.

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

◆ 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 2880 of file MOM_diabatic_driver.F90.

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

◆ 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 3070 of file MOM_diabatic_driver.F90.

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

◆ 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 2826 of file MOM_diabatic_driver.F90.

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

◆ 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 1857 of file MOM_diabatic_driver.F90.

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