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, G, GV, CS)
	This subroutine imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers. More...
subroutine, public	extract_diabatic_member (CS, opacity_CSp, optics_CSp, evap_CFL_limit, minimum_forcing_depth)
	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, CS)
	Routine called for adiabatic physics. More...
subroutine	diagnose_diabatic_diff_tendency (tv, h, temp_old, saln_old, dt, G, GV, 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, 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, CS, ncall)
	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, diag_to_Z_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, param_file, useALEalgorithm, diag, ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, ALE_sponge_CSp, diag_to_Z_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
integer	id_clock_mixedlayer
integer	id_clock_set_diffusivity
integer	id_clock_tracers
integer	id_clock_tridiag
integer	id_clock_pass
integer	id_clock_sponge
integer	id_clock_geothermal
integer	id_clock_differential_diff
integer	id_clock_remap
integer	id_clock_kpp

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(diabatic_cs), pointer	CS
	)

Routine called for adiabatic physics.

Parameters

[in,out]	g	ocean grid structure
[in,out]	h	thickness (m for Bouss or kg/m2 for non-Bouss)
[in,out]	tv	points to thermodynamic fields
[in,out]	fluxes	boundary fluxes
[in]	dt	time step (seconds)
[in]	gv	ocean vertical grid structure
	cs	module control structure

Definition at line 1503 of file MOM_diabatic_driver.F90.

References mom_tracer_flow_control::call_tracer_column_fns().

Referenced by mom::step_mom_thermo().

  type(ocean_grid_type),                    intent(inout) :: g      !< ocean grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h      !< thickness (m for Bouss or kg/m2 for non-Bouss)
  type(thermo_var_ptrs),                    intent(inout) :: tv     !< points to thermodynamic fields
  type(forcing),                            intent(inout) :: fluxes !< boundary fluxes
  real,                                     intent(in)    :: dt     !< time step (seconds)
  type(verticalgrid_type),                  intent(in)    :: gv     !< ocean vertical grid structure
  type(diabatic_cs),                        pointer       :: cs     !< module control structure

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: zeros  ! An array of zeros.

  zeros(:,:,:) = 0.0

  call call_tracer_column_fns(h, h, zeros, zeros, fluxes, zeros(:,:,1), dt, g, gv, tv, &
                              cs%optics, cs%tracer_flow_CSp, cs%debug)

Here is the call graph for this function:

Here is the caller graph for this function:

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,
		type(diag_to_z_cs), pointer	diag_to_Z_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]	time	current model time
[in]	g	model grid structure
[in]	param_file	the file to parse for parameter values
[in,out]	diag	regulates diagnostic output
	cs	module control structure
	tracer_flow_csp	points to control structure of tracer flow control module
	diag_to_z_csp	pointer to Z-diagnostics control structure

Definition at line 1756 of file MOM_diabatic_driver.F90.

Referenced by mom::initialize_mom().

  type(time_type),         intent(in)    :: time              !< current model time
  type(ocean_grid_type),   intent(in)    :: g                 !< model grid structure
  type(param_file_type),   intent(in)    :: param_file        !< the file to parse for parameter values
  type(diag_ctrl), target, intent(inout) :: diag              !< regulates diagnostic output
  type(diabatic_cs),       pointer       :: cs                !< module control structure
  type(tracer_flow_control_cs), pointer  :: tracer_flow_csp   !< points to control structure of tracer flow control module
  type(diag_to_z_cs),      pointer       :: diag_to_z_csp     !< pointer to Z-diagnostics control structure

! This "include" declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mod  = "MOM_diabatic_driver" ! This module's name.

  if (associated(cs)) then
    call mom_error(warning, "adiabatic_driver_init called with an "// &
                            "associated control structure.")
    return
  else ; allocate(cs) ; endif

  cs%diag => diag
  if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
  if (associated(diag_to_z_csp)) cs%diag_to_Z_CSp => diag_to_z_csp

! Set default, read and log parameters
  call log_version(param_file, mod, version, &
                   "The following parameters are used for diabatic processes.")

Here is the caller graph for this function:

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(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(diabatic_cs), pointer	CS
	)

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

Parameters

[in,out]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in,out]	u	zonal velocity (m/s)
[in,out]	v	meridional velocity (m/s)
[in,out]	h	thickness (m for Bouss / kg/m2 for non-Bouss)
[in,out]	tv	points to thermodynamic fields; unused have NULL ptrs
	hml	active mixed layer depth
[in,out]	fluxes	points to forcing fields; unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and related
[in,out]	adp	points to accelerations in momentum equations, to enable the later derived diagn, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment (seconds)
	cs	module control structure

Definition at line 235 of file MOM_diabatic_driver.F90.

References mom_ale_sponge::apply_ale_sponge(), mom_sponge::apply_sponge(), mom_diabatic_aux::applyboundaryfluxesinout(), mom_bulk_mixed_layer::bulkmixedlayer(), mom_diag_to_z::calc_zint_diags(), mom_forcing_type::calculatebuoyancyflux2d(), mom_tracer_flow_control::call_tracer_column_fns(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), diagnose_boundary_forcing_tendency(), diagnose_diabatic_diff_tendency(), diagnose_frazil_tendency(), mom_diabatic_aux::diagnosemldbydensitydifference(), mom_diapyc_energy_req::diapyc_energy_req_test(), mom_diffconvection::diffconvection_calculate(), mom_domains::do_group_pass(), mom_energetic_pbl::energetic_pbl_get_mld(), mom_diabatic_aux::find_uv_at_h(), mom_geothermal::geothermal(), id_clock_differential_diff, id_clock_entrain, id_clock_geothermal, id_clock_kpp, id_clock_mixedlayer, id_clock_pass, id_clock_remap, id_clock_set_diffusivity, id_clock_sponge, id_clock_tracers, id_clock_tridiag, mom_kpp::kpp_nonlocaltransport_saln(), mom_kpp::kpp_nonlocaltransport_temp(), mom_checksum_packages::mom_state_stats(), mom_variables::mom_thermovar_chksum(), mom_diag_mediator::query_averaging_enabled(), mom_regularize_layers::regularize_layers(), mom_opacity::set_opacity(), mom_tracer_diabatic::tracer_vertdiff(), and mom_wave_speed::wave_speeds().

  type(ocean_grid_type),                     intent(inout) :: g      !< ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: gv     !< ocean vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u      !< zonal velocity (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v      !< meridional velocity (m/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h      !< thickness (m for Bouss / kg/m2 for non-Bouss)
  type(thermo_var_ptrs),                     intent(inout) :: tv     !< points to thermodynamic fields; unused have NULL ptrs
  real, dimension(:,:),                      pointer       :: hml    !< active mixed layer depth
  type(forcing),                             intent(inout) :: fluxes !< points to forcing fields; unused fields have NULL ptrs
  type(vertvisc_type),                       intent(inout) :: visc   !< vertical viscosities, BBL properies, and related
  type(accel_diag_ptrs),                     intent(inout) :: adp    !< points to accelerations in momentum equations,
                                                                     !! to enable the later derived diagn, like energy budgets
  type(cont_diag_ptrs),                      intent(inout) :: cdp    !< points to terms in continuity equations
  real,                                      intent(in)    :: dt     !< time increment (seconds)
  type(diabatic_cs),                         pointer       :: cs     !< module control structure

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    ea,     &    ! amount of fluid entrained from the layer above within
                 ! one time step  (m for Bouss, kg/m^2 for non-Bouss)
    eb,     &    ! amount of fluid entrained from the layer below within
                 ! one time step  (m for Bouss, kg/m^2 for non-Bouss)
    kd,     &    ! diapycnal diffusivity of layers (m^2/sec)
    h_orig, &    ! initial layer thicknesses (m for Bouss, kg/m^2 for non-Bouss)
    h_prebound, &    ! initial layer thicknesses (m for Bouss, kg/m^2 for non-Bouss)
    hold,   &    ! layer thickness before diapycnal entrainment, and later
                 ! the initial layer thicknesses (if a mixed layer is used),
                 ! (m for Bouss, kg/m^2 for non-Bouss)
    dsv_dt, &    ! The partial derivatives of specific volume with temperature
    dsv_ds, &    ! and salinity in m^3/(kg K) and m^3/(kg ppt).
    ctke,   &    ! convective TKE requirements for each layer in J/m^2.
    u_h,    &    ! zonal and meridional velocities at thickness points after
    v_h          ! entrainment (m/s)

  real, dimension(SZI_(G),SZJ_(G),CS%nMode) :: &
    cn       ! baroclinic gravity wave speeds (formerly cg1 - BDM)

  real, dimension(SZI_(G),SZJ_(G)) :: &
    rcv_ml, &   ! coordinate density of mixed layer, used for applying sponges
    skinbuoyflux! 2d surface buoyancy flux (m2/s3), used by ePBL
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag                ! diagnostic array for thickness
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag             ! diagnostic array for temp
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag             ! diagnostic array for salinity
  real, dimension(SZI_(G),SZJ_(G))      :: tendency_2d           ! depth integrated content tendency for diagn
  real, dimension(SZI_(G),SZJ_(G))      :: tke_itidal_input_test ! override of energy input for testing (BDM)

  real :: net_ent  ! The net of ea-eb at an interface.

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
             ! These are targets so that the space can be shared with eaml & ebml.
    eatr, &  ! The equivalent of ea and eb for tracers, which differ from ea and
    ebtr     ! eb in that they tend to homogenize tracers in massless layers
             ! near the boundaries (m for Bouss and kg/m^2 for non-Bouss)

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), target :: &
    kd_int,   & ! diapycnal diffusivity of interfaces (m^2/s)
    kd_heat,  & ! diapycnal diffusivity of heat (m^2/s)
    kd_salt,  & ! diapycnal diffusivity of salt and passive tracers (m^2/s)
    kd_epbl,  & ! test array of diapycnal diffusivities at interfaces (m^2/s)
    tdif_flx, & ! diffusive diapycnal heat flux across interfaces (K m/s)
    tadv_flx, & ! advective diapycnal heat flux across interfaces (K m/s)
    sdif_flx, & ! diffusive diapycnal salt flux across interfaces (ppt m/s)
    sadv_flx    ! advective diapycnal salt flux across interfaces (ppt m/s)

  ! The following 5 variables are only used with a bulk mixed layer.
  real, pointer, dimension(:,:,:) :: &
    eaml, &  ! The equivalent of ea and eb due to mixed layer processes,
    ebml     ! (m for Bouss and kg/m^2 for non-Bouss).  These will be
             ! pointers to eatr and ebtr so as to reuse the memory as
             ! the arrays are not needed at the same time.

  integer :: kb(szi_(g),szj_(g)) ! index of the lightest layer denser
                                 ! than the buffer laye (nondimensional)

  real :: p_ref_cv(szi_(g))      ! Reference pressure for the potential
                                 ! density which defines the coordinate
                                 ! variable, set to P_Ref, in Pa.

  logical :: in_boundary(szi_(g)) ! True if there are no massive layers below,
                                  ! where massive is defined as sufficiently thick that
                                  ! the no-flux boundary conditions have not restricted
                                  ! the entrainment - usually sqrt(Kd*dt).

  real :: b_denom_1    ! The first term in the denominator of b1
                       ! (m for Bouss, kg/m^2 for non-Bouss)
  real :: h_neglect    ! A thickness that is so small it is usually lost
                       ! in roundoff and can be neglected
                       ! (m for Bouss and kg/m^2 for non-Bouss)
  real :: h_neglect2   ! h_neglect^2  (m^2 for Bouss, kg^2/m^4 for non-Bouss)
  real :: add_ent      ! Entrainment that needs to be added when mixing tracers
                       ! (m for Bouss and kg/m^2 for non-Bouss)
  real :: eaval        ! eaval is 2*ea at velocity grid points (m for Bouss, kg/m^2 for non-Bouss)
  real :: hval         ! hval is 2*h at velocity grid points (m for Bouss, kg/m^2 for non-Bouss)
  real :: h_tr         ! h_tr is h at tracer points with a tiny thickness
                       ! added to ensure positive definiteness (m for Bouss, kg/m^2 for non-Bouss)
  real :: tr_ea_bbl    ! The diffusive tracer thickness in the BBL that is
                       ! coupled to the bottom within a timestep (m)

  real :: htot(szib_(g))             ! The summed thickness from the bottom, in m.
  real :: b1(szib_(g)), d1(szib_(g)) ! b1, c1, and d1 are variables used by the
  real :: c1(szib_(g),szk_(g))       ! tridiagonal solver.

  real :: ent_int ! The diffusive entrainment rate at an interface
                  ! (H units = m for Bouss, kg/m^2 for non-Bouss).
  real :: dt_mix  ! amount of time over which to apply mixing (seconds)
  real :: idt     ! inverse time step (1/s)

  type(p3d) :: z_ptrs(7)  ! pointers to diagnostics to be interpolated to depth
  integer :: num_z_diags  ! number of diagnostics to be interpolated to depth
  integer :: z_ids(7)     ! id numbers of diagnostics to be interpolated to depth
  logical :: showcalltree ! If true, show the call tree
  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb, m

  integer :: ig, jg      ! global indices for testing testing itide point source (BDM)
  logical :: avg_enabled ! for testing internal tides (BDM)
  real :: kd_add_here    ! An added diffusivity in m2/s

  is   = g%isc  ; ie  = g%iec  ; js  = g%jsc  ; je  = g%jec ; nz = g%ke
  isq  = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  nkmb = gv%nk_rho_varies
  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
  kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0


  if (nz == 1) return
  showcalltree = calltree_showquery()
  if (showcalltree) call calltree_enter("diabatic(), MOM_diabatic_driver.F90")

  ! set equivalence between the same bits of memory for these arrays
  eaml => eatr ; ebml => ebtr

  ! inverse time step
  idt = 1.0 / dt

  if (.not. associated(cs)) call mom_error(fatal, "MOM_diabatic_driver: "// &
         "Module must be initialized before it is used.")

  if (cs%debug) then
    call mom_state_chksum("Start of diabatic ", u, v, h, g, gv, haloshift=0)
    call mom_forcing_chksum("Start of diabatic", fluxes, g, haloshift=0)
  endif
  if (cs%debugConservation) call mom_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, g)

  if (cs%debug_energy_req) &
    call diapyc_energy_req_test(h, dt, tv, g, gv, cs%diapyc_en_rec_CSp)


  call cpu_clock_begin(id_clock_set_diffusivity)
  call set_bbl_tke(u, v, h, fluxes, visc, g, gv, cs%set_diff_CSp)
  call cpu_clock_end(id_clock_set_diffusivity)

  ! Frazil formation keeps the temperature above the freezing point.
  ! make_frazil is deliberately called at both the beginning and at
  ! the end of the diabatic processes.
  if (ASSOCIATED(tv%T) .AND. ASSOCIATED(tv%frazil)) then

    if(cs%frazil_tendency_diag) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        temp_diag(i,j,k) = tv%T(i,j,k)
      enddo ; enddo ; enddo
    endif

    if (ASSOCIATED(fluxes%p_surf_full)) then
      call make_frazil(h, tv, g, gv, cs%diabatic_aux_CSp, fluxes%p_surf_full)
    else
      call make_frazil(h, tv, g, gv, cs%diabatic_aux_CSp)
    endif
    if (showcalltree) call calltree_waypoint("done with 1st make_frazil (diabatic)")

    if (cs%frazil_tendency_diag) then
      call diagnose_frazil_tendency(tv, h, temp_diag, dt, g, gv, cs, 1)
    endif

  endif

  if (cs%debugConservation) call mom_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, g)

  if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_orig,h,eaml,ebml)
    do k=1,nz ; do j=js,je ; do i=is,ie
      h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0
    enddo ; enddo ; enddo
  endif

  if (cs%use_geothermal) then
    call cpu_clock_begin(id_clock_geothermal)
    call geothermal(h, tv, dt, eaml, ebml, g, gv, cs%geothermal_CSp)
    call cpu_clock_end(id_clock_geothermal)
    if (showcalltree) call calltree_waypoint("geothermal (diabatic)")
    if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g)
  endif


  ! Whenever thickness changes let the diag manager know, target grids
  ! for vertical remapping may need to be regenerated.
  call diag_update_remap_grids(cs%diag)

  ! Set_opacity estimates the optical properties of the water column.
  ! It will need to be modified later to include information about the
  ! biological properties and layer thicknesses.
  if (associated(cs%optics)) &
    call set_opacity(cs%optics, fluxes, g, gv, cs%opacity_CSp)

  if (cs%bulkmixedlayer) then
    if (cs%debug) then
      call mom_forcing_chksum("Before mixedlayer", fluxes, g, haloshift=0)
    endif

    if (cs%ML_mix_first > 0.0) then
!  This subroutine
!    (1) Cools the mixed layer.
!    (2) Performs convective adjustment by mixed layer entrainment.
!    (3) Heats the mixed layer and causes it to detrain to
!        Monin-Obukhov depth or minimum mixed layer depth.
!    (4) Uses any remaining TKE to drive mixed layer entrainment.
!    (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate)
      call find_uv_at_h(u, v, h, u_h, v_h, g, gv)

      call cpu_clock_begin(id_clock_mixedlayer)
      if (cs%ML_mix_first < 1.0) then
        ! Changes: h, tv%T, tv%S, eaml and ebml  (G is also inout???)
        call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*cs%ML_mix_first, &
                            eaml,ebml, g, gv, cs%bulkmixedlayer_CSp, cs%optics, &
                            hml, cs%aggregate_FW_forcing, dt, last_call=.false.)
        if (cs%salt_reject_below_ML) &
          call insert_brine(h, tv, g, gv, fluxes, nkmb, cs%diabatic_aux_CSp, &
                            dt*cs%ML_mix_first, cs%id_brine_lay)
      else
        ! Changes: h, tv%T, tv%S, eaml and ebml  (G is also inout???)
        call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, &
                        g, gv, cs%bulkmixedlayer_CSp, cs%optics, &
                        hml, cs%aggregate_FW_forcing, dt, last_call=.true.)
      endif

      !  Keep salinity from falling below a small but positive threshold.
      !  This constraint is needed for SIS1 ice model, which can extract
      !  more salt than is present in the ocean. SIS2 does not suffer
      !  from this limitation, in which case we can let salinity=0 and still
      !  have salt conserved with SIS2 ice. So for SIS2, we can run with
      !  BOUND_SALINITY=False in MOM.F90.
      if (ASSOCIATED(tv%S) .and. ASSOCIATED(tv%salt_deficit)) &
        call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
      call cpu_clock_end(id_clock_mixedlayer)
      if (cs%debug) then
        call mom_state_chksum("After mixedlayer ", u, v, h, g, gv, haloshift=0)
        call mom_forcing_chksum("After mixedlayer", fluxes, g, haloshift=0)
      endif
      if (showcalltree) call calltree_waypoint("done with 1st bulkmixedlayer (diabatic)")
      if (cs%debugConservation) call mom_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, g)
    endif
  endif

  if (cs%debug) then
    call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, haloshift=0)
  endif
  if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then
    if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
      call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, eaml, ebml)
      if (cs%debug) then
        call hchksum(eaml, "after find_uv_at_h eaml",g%HI, scale=gv%H_to_m)
        call hchksum(ebml, "after find_uv_at_h ebml",g%HI, scale=gv%H_to_m)
      endif
    else
      call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
    endif
    if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
  endif

  if (cs%use_int_tides) then
    !   This block provides an interface for the unresolved low-mode internal
    ! tide module (BDM).

    ! PROVIDE ENERGY DISTRIBUTION (calculate time-varying energy source)
    call set_int_tide_input(u, v, h, tv, fluxes, cs%int_tide_input, dt, g, gv, &
                            cs%int_tide_input_CSp)
    ! CALCULATE MODAL VELOCITIES
    cn(:,:,:) = 0.0
    if (cs%uniform_cg) then
       ! SET TO CONSTANT VALUE TO TEST PROPAGATE CODE
       do m=1,cs%nMode ; cn(:,:,m) = cs%cg_test ; enddo
    else
       call wave_speeds(h, tv, g, gv, cs%nMode, cn, full_halos=.true.)
       ! uncomment the lines below for a hard-coded cn that changes linearly with latitude
       !do j=G%jsd,G%jed ; do i=G%isd,G%ied
       !  cn(i,j,:) = ((7.-1.)/14000000.)*G%geoLatBu(i,j) + (1.-((7.-1.)/14000000.)*-7000000.)
       !enddo ; enddo
    endif

    if (cs%int_tide_source_test) then
      ! BUILD 2D ARRAY WITH POINT SOURCE FOR TESTING
      !  This block of code should be moved into set_int_tide_input. -RWH
      tke_itidal_input_test(:,:) = 0.0
      avg_enabled = query_averaging_enabled(cs%diag,time_end=cs%time_end)
      if (cs%time_end <= cs%time_max_source) then
        do j=g%jsc,g%jec ; do i=g%isc,g%iec
          !INPUT ARBITRARY ENERGY POINT SOURCE
          if ((g%idg_offset + i == cs%int_tide_source_x) .and. &
              (g%jdg_offset + j == cs%int_tide_source_y)) then
            tke_itidal_input_test(i,j) = 1.0
          endif
        enddo ; enddo
      endif
      ! CALL ROUTINE USING PRESCRIBED KE FOR TESTING
      call propagate_int_tide(h, tv, cn, tke_itidal_input_test, &
                            cs%int_tide_input%tideamp, cs%int_tide_input%Nb, dt, g, gv, cs%int_tide_CSp)
    else
      ! CALL ROUTINE USING CALCULATED KE INPUT
      call propagate_int_tide(h, tv, cn, cs%int_tide_input%TKE_itidal_input, &
                              cs%int_tide_input%tideamp, cs%int_tide_input%Nb, dt, g, gv, cs%int_tide_CSp)
    endif
    if (showcalltree) call calltree_waypoint("done with propagate_int_tide (diabatic)")
  endif

  call cpu_clock_begin(id_clock_set_diffusivity)
  ! Sets: Kd, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S
  ! Also changes: visc%Kd_turb, visc%TKE_turb (not clear that TKE_turb is used as input ????)
  ! And sets visc%Kv_turb
  call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, cs%set_diff_CSp, kd, kd_int)
  call cpu_clock_end(id_clock_set_diffusivity)
  if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")

  if (cs%debug) then
    call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, haloshift=0)
    call mom_forcing_chksum("after set_diffusivity ", fluxes, g, haloshift=0)
    call mom_thermovar_chksum("after set_diffusivity ", tv, g)
    call hchksum(kd, "after set_diffusivity Kd",g%HI,haloshift=0)
    call hchksum(kd_int, "after set_diffusivity Kd_Int",g%HI,haloshift=0)
  endif


  if (cs%useKPP) then
    call cpu_clock_begin(id_clock_kpp)
    ! KPP needs the surface buoyancy flux but does not update state variables.
    ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
    ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux
    ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second)
    ! unlike other instances where the fluxes are integrated in time over a time-step.
    call calculatebuoyancyflux2d(g, gv, fluxes, cs%optics, h, tv%T, tv%S, tv, &
                                 cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
    ! The KPP scheme calculates boundary layer diffusivities and non-local transport.
    ! MOM6 implementation of KPP matches the boundary layer to zero interior diffusivity,
    ! since the matching to nonzero interior diffusivity can be problematic.
    ! Changes: Kd_int. Sets: KPP_NLTheat, KPP_NLTscalar

!$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat)
!$OMP do
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_salt(i,j,k) = kd_int(i,j,k)
        kd_heat(i,j,k) = kd_int(i,j,k)
      enddo ; enddo ; enddo
    if (associated(visc%Kd_extra_S)) then
!$OMP do
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
      enddo ; enddo ; enddo
    endif
    if (associated(visc%Kd_extra_T)) then
!$OMP do
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
      enddo ; enddo ; enddo
    endif
!$OMP end parallel

    call kpp_calculate(cs%KPP_CSp, g, gv, h, tv%T, tv%S, u, v, tv%eqn_of_state, &
      fluxes%ustar, cs%KPP_buoy_flux, kd_heat, kd_salt, visc%Kv_turb, cs%KPP_NLTheat, cs%KPP_NLTscalar)
!$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat)

    if (.not. cs%KPPisPassive) then
!$OMP do
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_int(i,j,k) = min( kd_salt(i,j,k),  kd_heat(i,j,k) )
      enddo ; enddo ; enddo
      if (associated(visc%Kd_extra_S)) then
!$OMP do
        do k=1,nz+1 ; do j=js,je ; do i=is,ie
          visc%Kd_extra_S(i,j,k) = kd_salt(i,j,k) - kd_int(i,j,k)
        enddo ; enddo ; enddo
      endif
      if (associated(visc%Kd_extra_T)) then
!$OMP do
        do k=1,nz+1 ; do j=js,je ; do i=is,ie
          visc%Kd_extra_T(i,j,k) = kd_heat(i,j,k) - kd_int(i,j,k)
        enddo ; enddo ; enddo
      endif
    endif ! not passive
!$OMP end parallel
    call cpu_clock_end(id_clock_kpp)
    if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")
    if (cs%debug) then
      call mom_state_chksum("after KPP", u, v, h, g, gv, haloshift=0)
      call mom_forcing_chksum("after KPP", fluxes, g, haloshift=0)
      call mom_thermovar_chksum("after KPP", tv, g)
      call hchksum(kd, "after KPP Kd",g%HI,haloshift=0)
      call hchksum(kd_int, "after KPP Kd_Int",g%HI,haloshift=0)
    endif

  endif  ! endif for KPP

  ! Check for static instabilities and increase Kd_int where unstable
  if (cs%useConvection) call diffconvection_calculate(cs%Conv_CSp, &
         g, gv, h, tv%T, tv%S, tv%eqn_of_state, kd_int)

  if (cs%useKPP) then

    call cpu_clock_begin(id_clock_kpp)
    if (cs%debug) then
      call hchksum(cs%KPP_temp_flux, "before KPP_applyNLT netHeat",g%HI,haloshift=0, scale=gv%H_to_m)
      call hchksum(cs%KPP_salt_flux, "before KPP_applyNLT netSalt",g%HI,haloshift=0, scale=gv%H_to_m)
      call hchksum(cs%KPP_NLTheat, "before KPP_applyNLT NLTheat",g%HI,haloshift=0)
      call hchksum(cs%KPP_NLTscalar, "before KPP_applyNLT NLTscalar",g%HI,haloshift=0)
    endif
    ! Apply non-local transport of heat and salt
    ! Changes: tv%T, tv%S
    call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat,   cs%KPP_temp_flux, dt, tv%T, tv%C_p)
    call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, dt, tv%S)
    call cpu_clock_end(id_clock_kpp)
    if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
    if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g)

    if (cs%debug) then
      call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, haloshift=0)
      call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, haloshift=0)
      call mom_thermovar_chksum("after KPP_applyNLT ", tv, g)
    endif

  endif ! endif for KPP

  ! Differential diffusion done here.
  ! Changes: tv%T, tv%S
  ! If using matching within the KPP scheme, then this step needs to provide
  ! a diffusivity and happen before KPP.  But generally in MOM, we do not match
  ! KPP boundary layer to interior, so this diffusivity can be computed when convenient.
  if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T)) then
    call cpu_clock_begin(id_clock_differential_diff)

    call differential_diffuse_t_s(h, tv, visc, dt, g, gv)
    call cpu_clock_end(id_clock_differential_diff)
    if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")
    if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g)

    ! increment heat and salt diffusivity.
    ! CS%useKPP==.true. already has extra_T and extra_S included
    if(.not. cs%useKPP) then
      do k=2,nz ; do j=js,je ; do i=is,ie
        kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
        kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
      enddo ; enddo ; enddo
    endif


  endif


  ! This block sets ea, eb from Kd or Kd_int.
  ! If using ALE algorithm, set ea=eb=Kd_int on interfaces for
  ! use in the tri-diagonal solver.
  ! Otherwise, call entrainment_diffusive() which sets ea and eb
  ! based on KD and target densities (ie. does remapping as well).
  if (cs%useALEalgorithm) then

    do j=js,je ; do i=is,ie
      ea(i,j,1) = 0.
    enddo ; enddo
!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_neglect,h,ea,GV,dt,Kd_int,eb) &
!$OMP                          private(hval)
    do k=2,nz ; do j=js,je ; do i=is,ie
      hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
      ea(i,j,k) = (gv%m_to_H**2) * dt * hval * kd_int(i,j,k)
      eb(i,j,k-1) = ea(i,j,k)
    enddo ; enddo ; enddo
    do j=js,je ; do i=is,ie
      eb(i,j,nz) = 0.
    enddo ; enddo
    if (showcalltree) call calltree_waypoint("done setting ea,eb from Kd_int (diabatic)")

  else ! .not. CS%useALEalgorithm
    ! When not using ALE, calculate layer entrainments/detrainments from
    ! diffusivities and differences between layer and target densities
    call cpu_clock_begin(id_clock_entrain)
    ! Calculate appropriately limited diapycnal mass fluxes to account
    ! for diapycnal diffusion and advection.  Sets: ea, eb. Changes: kb
    call entrainment_diffusive(u, v, h, tv, fluxes, dt, g, gv, cs%entrain_diffusive_CSp, &
                               ea, eb, kb, kd_lay=kd, kd_int=kd_int)
    call cpu_clock_end(id_clock_entrain)
    if (showcalltree) call calltree_waypoint("done with Entrainment_diffusive (diabatic)")

  endif ! endif for (CS%useALEalgorithm)

  if (cs%debug) then
    call mom_forcing_chksum("after calc_entrain ", fluxes, g, haloshift=0)
    call mom_thermovar_chksum("after calc_entrain ", tv, g)
    call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, haloshift=0)
    call hchksum(ea, "after calc_entrain ea", g%HI, haloshift=0, scale=gv%H_to_m)
    call hchksum(eb, "after calc_entrain eb", g%HI, haloshift=0, scale=gv%H_to_m)
  endif

  ! Apply forcing when using the ALE algorithm
  if (cs%useALEalgorithm) then
    call cpu_clock_begin(id_clock_remap)

    ! Changes made to following fields:  h, tv%T and tv%S.

    ! save prior values for diagnostics
    if(cs%boundary_forcing_tendency_diag) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        h_diag(i,j,k)    = h(i,j,k)
        temp_diag(i,j,k) = tv%T(i,j,k)
        saln_diag(i,j,k) = tv%S(i,j,k)
      enddo ; enddo ; enddo
    endif

    do k=1,nz ; do j=js,je ; do i=is,ie
        h_prebound(i,j,k) = h(i,j,k)
    enddo ; enddo ; enddo
    if (cs%use_energetic_PBL) then

      skinbuoyflux(:,:) = 0.0
      call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, dt, fluxes, cs%optics, &
              h, tv, cs%aggregate_FW_forcing, cs%evap_CFL_limit,                         &
              cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux)

      if (cs%debug) then
        call hchksum(ea, "after applyBoundaryFluxes ea",g%HI,haloshift=0, scale=gv%H_to_m)
        call hchksum(eb, "after applyBoundaryFluxes eb",g%HI,haloshift=0, scale=gv%H_to_m)
        call hchksum(ctke, "after applyBoundaryFluxes cTKE",g%HI,haloshift=0)
        call hchksum(dsv_dt, "after applyBoundaryFluxes dSV_dT",g%HI,haloshift=0)
        call hchksum(dsv_ds, "after applyBoundaryFluxes dSV_dS",g%HI,haloshift=0)
      endif

      call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
      call energetic_pbl(h, u_h, v_h, tv, fluxes, dt, kd_epbl, g, gv, &
           cs%energetic_PBL_CSp, dsv_dt, dsv_ds, ctke, skinbuoyflux)

      ! If visc%MLD exists, copy the ePBL's MLD into it
      if (associated(visc%MLD)) then
        call energetic_pbl_get_mld(cs%energetic_PBL_CSp, visc%MLD, g)
        call pass_var(visc%MLD, g%domain, halo=1)
        hml(:,:) = visc%MLD(:,:)
      endif

      ! Augment the diffusivities due to those diagnosed in energetic_PBL.
      do k=2,nz ; do j=js,je ; do i=is,ie

        if (cs%ePBL_is_additive) then
          kd_add_here = kd_epbl(i,j,k)
          visc%Kv_turb(i,j,k) = visc%Kv_turb(i,j,k) + kd_epbl(i,j,k)
        else
          kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_turb(i,j,k), 0.0)
          visc%Kv_turb(i,j,k) = max(visc%Kv_turb(i,j,k), kd_epbl(i,j,k))
        endif
        ent_int = kd_add_here * (gv%m_to_H**2 * dt) / &
                    (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect)
        eb(i,j,k-1) = eb(i,j,k-1) + ent_int
        ea(i,j,k) = ea(i,j,k) + ent_int
        kd_int(i,j,k)  = kd_int(i,j,k) + kd_add_here

        ! for diagnostics
        kd_heat(i,j,k) = kd_heat(i,j,k) + kd_int(i,j,k)
        kd_salt(i,j,k) = kd_salt(i,j,k) + kd_int(i,j,k)

      enddo ; enddo ; enddo

      if (cs%debug) then
        call hchksum(ea, "after ePBL ea",g%HI,haloshift=0, scale=gv%H_to_m)
        call hchksum(eb, "after ePBL eb",g%HI,haloshift=0, scale=gv%H_to_m)
        call hchksum(kd_epbl, "after ePBL Kd_ePBL",g%HI,haloshift=0)
      endif

    else
      call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, dt, fluxes, cs%optics, &
                                    h, tv, cs%aggregate_FW_forcing, &
                                    cs%evap_CFL_limit, cs%minimum_forcing_depth)

    endif   ! endif for CS%use_energetic_PBL

    ! diagnose the tendencies due to boundary forcing
    if(cs%boundary_forcing_tendency_diag) then
      call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, g, gv, cs)
    endif

    call cpu_clock_end(id_clock_remap)
    if (cs%debug) then
      call mom_forcing_chksum("after applyBoundaryFluxes ", fluxes, g, haloshift=0)
      call mom_thermovar_chksum("after applyBoundaryFluxes ", tv, g)
      call mom_state_chksum("after applyBoundaryFluxes ", u, v, h, g, gv, haloshift=0)
    endif
    if (showcalltree) call calltree_waypoint("done with applyBoundaryFluxes (diabatic)")
    if (cs%debugConservation)  call mom_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, g)

  endif   ! endif for (CS%useALEalgorithm)

  ! Update h according to divergence of the difference between
  ! ea and eb. We keep a record of the original h in hold.
  ! In the following, the checks for negative values are to guard
  ! against instances where entrainment drives a layer to
  ! negative thickness.  This situation will never happen if
  ! enough iterations are permitted in Calculate_Entrainment.
  ! Even if too few iterations are allowed, it is still guarded
  ! against.  In other words the checks are probably unnecessary.
!$OMP parallel do default(none) shared(is,ie,js,je,nz,hold,h,eb,ea,GV)
  do j=js,je
    do i=is,ie
      hold(i,j,1) = h(i,j,1)
      h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2))
      hold(i,j,nz) = h(i,j,nz)
      h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1))
      if (h(i,j,1) <= 0.0) then
        h(i,j,1) = gv%Angstrom
      endif
      if (h(i,j,nz) <= 0.0) then
        h(i,j,nz) = gv%Angstrom
      endif
    enddo
    do k=2,nz-1 ; do i=is,ie
      hold(i,j,k) = h(i,j,k)
      h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
                    (eb(i,j,k) - ea(i,j,k+1)))
      if (h(i,j,k) <= 0.0) then
        h(i,j,k) = gv%Angstrom
      endif
    enddo ; enddo
  enddo
  if (cs%debug) then
    call mom_state_chksum("after negative check ", u, v, h, g, gv, haloshift=0)
    call mom_forcing_chksum("after negative check ", fluxes, g, haloshift=0)
    call mom_thermovar_chksum("after negative check ", tv, g)
  endif
  if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")
  if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g)


  ! Here, T and S are updated according to ea and eb.
  ! If using the bulk mixed layer, T and S are also updated
  ! by surface fluxes (in fluxes%*).
  ! This is a very long block.
  if (cs%bulkmixedlayer) then

    if (ASSOCIATED(tv%T)) then
      call cpu_clock_begin(id_clock_tridiag)
      ! Temperature and salinity (as state variables) are treated
      ! differently from other tracers to insure massless layers that
      ! are lighter than the mixed layer have temperatures and salinities
      ! that correspond to their prescribed densities.
      if (cs%massless_match_targets) then
!$OMP parallel do default (none) shared(is,ie,js,je,nkmb,hold,h_neglect,eb,ea,nz,kb,tv) &
!$OMP                           private(h_tr,b1,d1,c1,b_denom_1)
        do j=js,je
          do i=is,ie
            h_tr = hold(i,j,1) + h_neglect
            b1(i) = 1.0 / (h_tr + eb(i,j,1))
            d1(i) = h_tr * b1(i)
            tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1))
            tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1))
          enddo
          do k=2,nkmb ; do i=is,ie
            c1(i,k) = eb(i,j,k-1) * b1(i)
            h_tr = hold(i,j,k) + h_neglect
            b_denom_1 = h_tr + d1(i)*ea(i,j,k)
            b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
            if (k<nkmb) d1(i) = b_denom_1 * b1(i)
            tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
            tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
          enddo ; enddo

          do k=nkmb+1,nz ; do i=is,ie
            if (k == kb(i,j)) then
              c1(i,k) = eb(i,j,k-1) * b1(i)
              d1(i) = (((eb(i,j,nkmb)-eb(i,j,k-1)) + hold(i,j,nkmb) + h_neglect) + &
                       d1(i)*ea(i,j,nkmb)) * b1(i)
              h_tr = hold(i,j,k) + h_neglect
              b_denom_1 = h_tr + d1(i)*ea(i,j,k)
              b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
              d1(i) = b_denom_1 * b1(i)
              tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,nkmb))
              tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,nkmb))
            elseif (k > kb(i,j)) then
              c1(i,k) = eb(i,j,k-1) * b1(i)
              h_tr = hold(i,j,k) + h_neglect
              b_denom_1 = h_tr + d1(i)*ea(i,j,k)
              b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
              d1(i) = b_denom_1 * b1(i)
              tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
              tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
            elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j))
              !   The bottommost buffer layer might entrain all the mass from some
              ! of the interior layers that are thin and lighter in the coordinate
              ! density than that buffer layer.  The T and S of these newly
              ! massless interior layers are unchanged.
              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)
              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)
            endif
          enddo ; enddo

          do k=nz-1,nkmb,-1 ; do i=is,ie
            if (k >= kb(i,j)) then
              tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
              tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
            endif
          enddo ; enddo
          do i=is,ie ; if (kb(i,j) <= nz) then
            tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j))
            tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j))
          endif ; enddo
          do k=nkmb-1,1,-1 ; do i=is,ie
            tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
            tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
          enddo ; enddo
        enddo ! end of j loop
      else ! .not. massless_match_targets
        ! This simpler form allows T & S to be too dense for the layers
        ! between the buffer layers and the interior.
        ! Changes: T, S
        if (cs%tracer_tridiag) then
          call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
          call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
        else
          call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
        endif
      endif ! massless_match_targets
      call cpu_clock_end(id_clock_tridiag)

    endif ! endif for ASSOCIATED(T)
    if (cs%debugConservation) call mom_state_stats('BML tridiag', u, v, h, tv%T, tv%S, g)

    if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
      ! The mixed layer code has already been called, but there is some needed
      ! bookkeeping.
!$OMP parallel do default(none) shared(is,ie,js,je,nz,hold,h_orig,ea,eaml,eb,ebml)
      do k=1,nz ; do j=js,je ; do i=is,ie
        hold(i,j,k) = h_orig(i,j,k)
        ea(i,j,k) = ea(i,j,k) + eaml(i,j,k)
        eb(i,j,k) = eb(i,j,k) + ebml(i,j,k)
      enddo ; enddo ; enddo
      if (cs%debug) then
        call hchksum(ea, "after ea = ea + eaml",g%HI,haloshift=0, scale=gv%H_to_m)
        call hchksum(eb, "after eb = eb + ebml",g%HI,haloshift=0, scale=gv%H_to_m)
      endif
    endif

    if (cs%ML_mix_first < 1.0) then
    !  Call the mixed layer code now, perhaps for a second time.
    !  This subroutine (1)  Cools the mixed layer.
    !    (2) Performs convective adjustment by mixed layer entrainment.
    !    (3) Heats the mixed layer and causes it to detrain to
    !        Monin-Obukhov depth or minimum mixed layer depth.
    !    (4) Uses any remaining TKE to drive mixed layer entrainment.
    !    (5) Possibly splits the buffer layer into two isopycnal layers.

      call find_uv_at_h(u, v, hold, u_h, v_h, g, gv, ea, eb)
      if (cs%debug) call mom_state_chksum("find_uv_at_h1 ", u, v, h, g, gv, haloshift=0)

      dt_mix = min(dt,dt*(1.0 - cs%ML_mix_first))
      call cpu_clock_begin(id_clock_mixedlayer)
      ! Changes: h, tv%T, tv%S, ea and eb  (G is also inout???)
      call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, &
                      g, gv, cs%bulkmixedlayer_CSp, cs%optics, &
                      hml, cs%aggregate_FW_forcing, dt, last_call=.true.)

      if (cs%salt_reject_below_ML) &
        call insert_brine(h, tv, g, gv, fluxes, nkmb, cs%diabatic_aux_CSp, dt_mix, &
                          cs%id_brine_lay)

      !  Keep salinity from falling below a small but positive threshold.
      !  This constraint is needed for SIS1 ice model, which can extract
      !  more salt than is present in the ocean. SIS2 does not suffer
      !  from this limitation, in which case we can let salinity=0 and still
      !  have salt conserved with SIS2 ice. So for SIS2, we can run with
      !  BOUND_SALINITY=False in MOM.F90.
      if (ASSOCIATED(tv%S) .and. ASSOCIATED(tv%salt_deficit)) &
        call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

      call cpu_clock_end(id_clock_mixedlayer)
      if (showcalltree) call calltree_waypoint("done with 2nd bulkmixedlayer (diabatic)")
      if (cs%debugConservation) call mom_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, g)
    endif

  else  ! following block for when NOT using BULKMIXEDLAYER


    ! calculate change in temperature & salinity due to dia-coordinate surface diffusion
    if (ASSOCIATED(tv%T)) then

      if (cs%debug) then
        call hchksum(ea, "before triDiagTS ea ",g%HI,haloshift=0, scale=gv%H_to_m)
        call hchksum(eb, "before triDiagTS eb ",g%HI,haloshift=0, scale=gv%H_to_m)
      endif
      call cpu_clock_begin(id_clock_tridiag)

      if(cs%diabatic_diff_tendency_diag) then
        do k=1,nz ; do j=js,je ; do i=is,ie
          temp_diag(i,j,k) = tv%T(i,j,k)
          saln_diag(i,j,k) = tv%S(i,j,k)
        enddo ; enddo ; enddo
      endif

      ! Changes T and S via the tridiagonal solver; no change to h
      if(cs%tracer_tridiag) then
          call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
          call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
      else
        call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
      endif

      ! diagnose temperature, salinity, heat, and salt tendencies
      if(cs%diabatic_diff_tendency_diag) then
         call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, cs)
      endif

      call cpu_clock_end(id_clock_tridiag)
      if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")

    endif  ! endif corresponding to if (ASSOCIATED(tv%T))
    if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g)


  endif  ! endif for the BULKMIXEDLAYER block


  if (cs%debug) then
    call mom_state_chksum("after mixed layer ", u, v, h, g, gv, haloshift=0)
    call mom_thermovar_chksum("after mixed layer ", tv, g)
    call hchksum(ea, "after mixed layer ea", g%HI, scale=gv%H_to_m)
    call hchksum(eb, "after mixed layer eb", g%HI, scale=gv%H_to_m)
  endif

  if (.not. cs%useALEalgorithm) then
    call cpu_clock_begin(id_clock_remap)
    call regularize_layers(h, tv, dt, ea, eb, g, gv, cs%regularize_layers_CSp)
    call cpu_clock_end(id_clock_remap)
    if (showcalltree) call calltree_waypoint("done with regularize_layers (diabatic)")
    if (cs%debugConservation) call mom_state_stats('regularize_layers', u, v, h, tv%T, tv%S, g)
  endif

  ! Whenever thickness changes let the diag manager know, as the
  ! target grids for vertical remapping may need to be regenerated.
  call diag_update_remap_grids(cs%diag)

  ! diagnostics
  if ((cs%id_Tdif > 0) .or. (cs%id_Tdif_z > 0) .or. &
      (cs%id_Tadv > 0) .or. (cs%id_Tadv_z > 0)) then
    do j=js,je ; do i=is,ie
      tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
      tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
    enddo ; enddo
!$OMP parallel do default(none) shared(is,ie,js,je,nz,Tdif_flx,Idt,ea,eb,Tadv_flx,tv)
    do k=2,nz ; do j=js,je ; do i=is,ie
      tdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
                        (tv%T(i,j,k-1) - tv%T(i,j,k))
      tadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
                    0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
    enddo ; enddo ; enddo
  endif
  if ((cs%id_Sdif > 0) .or. (cs%id_Sdif_z > 0) .or. &
      (cs%id_Sadv > 0) .or. (cs%id_Sadv_z > 0)) then
    do j=js,je ; do i=is,ie
      sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
      sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
    enddo ; enddo
!$OMP parallel do default(none) shared(is,ie,js,je,nz,Sdif_flx,Idt,ea,eb,Sadv_flx,tv)
    do k=2,nz ; do j=js,je ; do i=is,ie
      sdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
                        (tv%S(i,j,k-1) - tv%S(i,j,k))
      sadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
                    0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
    enddo ; enddo ; enddo
  endif

  ! mixing of passive tracers from massless boundary layers to interior
  call cpu_clock_begin(id_clock_tracers)
  if (cs%mix_boundary_tracers) then
    tr_ea_bbl = sqrt(dt*cs%Kd_BBL_tr)
!$OMP parallel do default(none) shared(is,ie,js,je,ebtr,nz,G,GV,h,dt,CS,h_neglect,  &
!$OMP                                  ea,eb,Tr_ea_BBL,eatr,visc,hold,h_neglect2 )  &
!$OMP                          private(htot,in_boundary,add_ent)
    do j=js,je
      do i=is,ie
        ebtr(i,j,nz) = eb(i,j,nz)
        htot(i) = 0.0
        in_boundary(i) = (g%mask2dT(i,j) > 0.0)
      enddo
      do k=nz,2,-1 ; do i=is,ie
        if (in_boundary(i)) then
          htot(i) = htot(i) + h(i,j,k)
          !   If diapycnal mixing has been suppressed because this is a massless
          ! layer near the bottom, add some mixing of tracers between these
          ! layers.  This flux is based on the harmonic mean of the two
          ! thicknesses, as this corresponds pretty closely (to within
          ! differences in the density jumps between layers) with what is done
          ! in the calculation of the fluxes in the first place.  Kd_min_tr
          ! should be much less than the values that have been set in Kd,
          ! perhaps a molecular diffusivity.
          add_ent = ((dt * cs%Kd_min_tr) * gv%m_to_H**2) * &
                    ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
                     (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
                    0.5*(ea(i,j,k) + eb(i,j,k-1))
          if (htot(i) < tr_ea_bbl) then
            add_ent = max(0.0, add_ent, &
                          (tr_ea_bbl - htot(i)) - min(ea(i,j,k),eb(i,j,k-1)))
          elseif (add_ent < 0.0) then
            add_ent = 0.0 ; in_boundary(i) = .false.
          endif

          ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
          eatr(i,j,k) = ea(i,j,k) + add_ent
        else
          ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k)
        endif
        if (associated(visc%Kd_extra_S)) then ; if (visc%Kd_extra_S(i,j,k) > 0.0) then
          add_ent = ((dt * visc%Kd_extra_S(i,j,k)) * gv%m_to_H**2) / &
             (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
              h_neglect)
          ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
          eatr(i,j,k) = eatr(i,j,k) + add_ent
        endif ; endif
      enddo ; enddo
      do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo

    enddo

    if (cs%useALEalgorithm) then
    ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
    ! so hold should be h_orig
      call call_tracer_column_fns(h_prebound, h, ea, eb, fluxes, hml, dt, g, gv, tv, &
                                cs%optics, cs%tracer_flow_CSp, cs%debug, &
                                evap_cfl_limit = cs%evap_CFL_limit, &
                                minimum_forcing_depth = cs%minimum_forcing_depth)
    else
      call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
                                cs%optics, cs%tracer_flow_CSp, cs%debug)
    endif

  elseif (associated(visc%Kd_extra_S)) then  ! extra diffusivity for passive tracers

    do j=js,je ; do i=is,ie
      ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1)
    enddo ; enddo
!$OMP parallel do default(none) shared(nz,is,ie,js,je,visc,dt,GV,h,hold,h_neglect,&
!$OMP                                  ebtr,eb,eatr,ea )                          &
!$OMP                          private(add_ent)
    do k=nz,2,-1 ; do j=js,je ; do i=is,ie
      if (visc%Kd_extra_S(i,j,k) > 0.0) then
        add_ent = ((dt * visc%Kd_extra_S(i,j,k)) * gv%m_to_H**2) / &
           (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
            h_neglect)
      else
        add_ent = 0.0
      endif
      ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
      eatr(i,j,k) = ea(i,j,k) + add_ent
    enddo ; enddo ; enddo

    if (cs%useALEalgorithm) then
    ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
      call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
                                  cs%optics, cs%tracer_flow_CSp, cs%debug,&
                                  evap_cfl_limit = cs%evap_CFL_limit, &
                                  minimum_forcing_depth = cs%minimum_forcing_depth)
    else
      call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
                                  cs%optics, cs%tracer_flow_CSp, cs%debug)
    endif

  else
    if (cs%useALEalgorithm) then
    ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
      call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
                                  cs%optics, cs%tracer_flow_CSp, cs%debug, &
                                  evap_cfl_limit = cs%evap_CFL_limit, &
                                  minimum_forcing_depth = cs%minimum_forcing_depth)
    else
      call call_tracer_column_fns(hold, h, ea, eb, fluxes, hml, dt, g, gv, tv, &
                                  cs%optics, cs%tracer_flow_CSp, cs%debug)
    endif

  endif  ! (CS%mix_boundary_tracers)



  call cpu_clock_end(id_clock_tracers)


  ! sponges
  if (cs%use_sponge) then
    call cpu_clock_begin(id_clock_sponge)
    if (associated(cs%ALE_sponge_CSp)) then
      ! ALE sponge
      call apply_ale_sponge(h, dt, g, cs%ALE_sponge_CSp)
    else
      ! Layer mode sponge
      if (cs%bulkmixedlayer .and. ASSOCIATED(tv%eqn_of_state)) then
        do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo
!$OMP parallel do default(none) shared(js,je,p_ref_cv,Rcv_ml,is,ie,tv)
        do j=js,je
           call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, rcv_ml(:,j), &
                               is, ie-is+1, tv%eqn_of_state)
        enddo
        call apply_sponge(h, dt, g, gv, ea, eb, cs%sponge_CSp, rcv_ml)
      else
        call apply_sponge(h, dt, g, gv, ea, eb, cs%sponge_CSp)
      endif
    endif
    call cpu_clock_end(id_clock_sponge)
    if (cs%debug) then
      call mom_state_chksum("apply_sponge ", u, v, h, g, gv, haloshift=0)
      call mom_thermovar_chksum("apply_sponge ", tv, g)
    endif
  endif ! CS%use_sponge


!$OMP parallel default(none) shared(is,ie,js,je,nz,CDp,Idt,G,GV,ea,eb,CS,hold) private(net_ent)
!   Save the diapycnal mass fluxes as a diagnostic field.
  if (ASSOCIATED(cdp%diapyc_vel)) then
!$OMP do
    do j=js,je
      do k=2,nz ; do i=is,ie
        cdp%diapyc_vel(i,j,k) = idt * (gv%H_to_m * (ea(i,j,k) - eb(i,j,k-1)))
      enddo ; enddo
      do i=is,ie
        cdp%diapyc_vel(i,j,1) = 0.0
        cdp%diapyc_vel(i,j,nz+1) = 0.0
      enddo
    enddo
  endif

! For momentum, it is only the net flux that homogenizes within
! the mixed layer.  Vertical viscosity that is proportional to the
! mixed layer turbulence is applied elsewhere.
  if (cs%bulkmixedlayer) then
    if (cs%debug) then
      call hchksum(ea, "before net flux rearrangement ea",g%HI, scale=gv%H_to_m)
      call hchksum(eb, "before net flux rearrangement eb",g%HI, scale=gv%H_to_m)
    endif
!$OMP do
    do j=js,je
      do k=2,gv%nkml ; do i=is,ie
        net_ent = ea(i,j,k) - eb(i,j,k-1)
        ea(i,j,k) = max(net_ent, 0.0)
        eb(i,j,k-1) = max(-net_ent, 0.0)
      enddo ; enddo
    enddo
    if (cs%debug) then
      call hchksum(ea, "after net flux rearrangement ea",g%HI, scale=gv%H_to_m)
      call hchksum(eb, "after net flux rearrangement eb",g%HI, scale=gv%H_to_m)
    endif
  endif


! Initialize halo regions of ea, eb, and hold to default values.
!$OMP do
  do k=1,nz
    do i=is-1,ie+1
      hold(i,js-1,k) = gv%Angstrom ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0
      hold(i,je+1,k) = gv%Angstrom ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0
    enddo
    do j=js,je
      hold(is-1,j,k) = gv%Angstrom ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0
      hold(ie+1,j,k) = gv%Angstrom ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0
    enddo
  enddo
!$OMP end parallel

  call cpu_clock_begin(id_clock_pass)
  if (g%symmetric) then
    call create_group_pass(cs%pass_hold_eb_ea,hold,g%Domain,to_all+omit_corners,halo=1)
    call create_group_pass(cs%pass_hold_eb_ea,eb,g%Domain,to_all+omit_corners,halo=1)
    call create_group_pass(cs%pass_hold_eb_ea,ea,g%Domain,to_all+omit_corners,halo=1)
  else
    call create_group_pass(cs%pass_hold_eb_ea,hold,g%Domain,to_west+to_south+omit_corners,halo=1)
    call create_group_pass(cs%pass_hold_eb_ea,eb,g%Domain,to_west+to_south+omit_corners,halo=1)
    call create_group_pass(cs%pass_hold_eb_ea,ea,g%Domain,to_west+to_south+omit_corners,halo=1)
  endif
  call do_group_pass(cs%pass_hold_eb_ea,g%Domain)
  call cpu_clock_end(id_clock_pass)

  if (.not. cs%useALEalgorithm) then
    !  Use a tridiagonal solver to determine effect of the diapycnal
    !  advection on velocity field. It is assumed that water leaves
    !  or enters the ocean with the surface velocity.
    if (cs%debug) then
      call mom_state_chksum("before u/v tridiag ", u, v, h, g, gv, haloshift=0)
      call hchksum(ea, "before u/v tridiag ea",g%HI, scale=gv%H_to_m)
      call hchksum(eb, "before u/v tridiag eb",g%HI, scale=gv%H_to_m)
      call hchksum(hold, "before u/v tridiag hold",g%HI, scale=gv%H_to_m)
    endif
    call cpu_clock_begin(id_clock_tridiag)
!$OMP parallel do default(none) shared(js,je,Isq,Ieq,ADp,u,hold,ea,h_neglect,eb,nz,Idt) &
!$OMP                          private(hval,b1,d1,c1,eaval)
    do j=js,je
      do i=isq,ieq
        if (ASSOCIATED(adp%du_dt_dia)) adp%du_dt_dia(i,j,1) = u(i,j,1)
        hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect
        b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))
        d1(i) = hval * b1(i)
        u(i,j,1) = b1(i) * (hval * u(i,j,1))
      enddo
      do k=2,nz ; do i=isq,ieq
        if (ASSOCIATED(adp%du_dt_dia)) adp%du_dt_dia(i,j,k) = u(i,j,k)
        c1(i,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(i)
        eaval = ea(i,j,k) + ea(i+1,j,k)
        hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect
        b1(i) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(i)*eaval))
        d1(i) = (hval + d1(i)*eaval) * b1(i)
        u(i,j,k) = (hval*u(i,j,k) + eaval*u(i,j,k-1))*b1(i)
      enddo ; enddo
      do k=nz-1,1,-1 ; do i=isq,ieq
        u(i,j,k) = u(i,j,k) + c1(i,k+1)*u(i,j,k+1)
        if (ASSOCIATED(adp%du_dt_dia)) &
          adp%du_dt_dia(i,j,k) = (u(i,j,k) - adp%du_dt_dia(i,j,k)) * idt
      enddo ; enddo
      if (ASSOCIATED(adp%du_dt_dia)) then
        do i=isq,ieq
          adp%du_dt_dia(i,j,nz) = (u(i,j,nz)-adp%du_dt_dia(i,j,nz)) * idt
        enddo
      endif
    enddo
    if (cs%debug) then
      call mom_state_chksum("aft 1st loop tridiag ", u, v, h, g, gv, haloshift=0)
    endif
!$OMP parallel do default(none) shared(Jsq,Jeq,is,ie,ADp,v,hold,ea,h_neglect,eb,nz,Idt) &
!$OMP                          private(hval,b1,d1,c1,eaval)
    do j=jsq,jeq
      do i=is,ie
        if (ASSOCIATED(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,1) = v(i,j,1)
        hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect
        b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))
        d1(i) = hval * b1(i)
        v(i,j,1) = b1(i) * (hval * v(i,j,1))
      enddo
      do k=2,nz ; do i=is,ie
        if (ASSOCIATED(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,k) = v(i,j,k)
        c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)
        eaval = ea(i,j,k) + ea(i,j+1,k)
        hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect
        b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))
        d1(i) = (hval + d1(i)*eaval) * b1(i)
        v(i,j,k) = (hval*v(i,j,k) + eaval*v(i,j,k-1))*b1(i)
      enddo ; enddo
      do k=nz-1,1,-1 ; do i=is,ie
        v(i,j,k) = v(i,j,k) + c1(i,k+1)*v(i,j,k+1)
        if (ASSOCIATED(adp%dv_dt_dia)) &
          adp%dv_dt_dia(i,j,k) = (v(i,j,k) - adp%dv_dt_dia(i,j,k)) * idt
      enddo ; enddo
      if (ASSOCIATED(adp%dv_dt_dia)) then
        do i=is,ie
          adp%dv_dt_dia(i,j,nz) = (v(i,j,nz)-adp%dv_dt_dia(i,j,nz)) * idt
        enddo
      endif
    enddo
    call cpu_clock_end(id_clock_tridiag)
    if (cs%debug) then
      call mom_state_chksum("after u/v tridiag ", u, v, h, g, gv, haloshift=0)
    endif
  endif ! useALEalgorithm

  ! Frazil formation keeps temperature above the freezing point.
  ! make_frazil is deliberately called at both the beginning and at
  ! the end of the diabatic processes.
  if (ASSOCIATED(tv%T) .AND. ASSOCIATED(tv%frazil)) then

    if(cs%frazil_tendency_diag) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        temp_diag(i,j,k) = tv%T(i,j,k)
      enddo ; enddo ; enddo
    endif

    if (ASSOCIATED(fluxes%p_surf_full)) then
      call make_frazil(h, tv, g, gv, cs%diabatic_aux_CSp, fluxes%p_surf_full)
    else
      call make_frazil(h, tv, g, gv, cs%diabatic_aux_CSp)
    endif

    if (cs%frazil_tendency_diag) then
      call diagnose_frazil_tendency(tv, h, temp_diag, dt, g, gv, cs, 2)
    endif

    if (showcalltree) call calltree_waypoint("done with 2nd make_frazil (diabatic)")
    if (cs%debugConservation) call mom_state_stats('2nd make_frazil', u, v, h, tv%T, tv%S, g)

  endif  ! endif for frazil


  ! Diagnose the diapycnal diffusivities and other related quantities.
  if (cs%id_Kd_interface > 0) call post_data(cs%id_Kd_interface, kd_int,  cs%diag)
  if (cs%id_Kd_heat      > 0) call post_data(cs%id_Kd_heat,      kd_heat, cs%diag)
  if (cs%id_Kd_salt      > 0) call post_data(cs%id_Kd_salt,      kd_salt, cs%diag)
  if (cs%id_Kd_ePBL      > 0) call post_data(cs%id_Kd_ePBL,      kd_epbl, cs%diag)

  if (cs%id_ea       > 0) call post_data(cs%id_ea,       ea, cs%diag)
  if (cs%id_eb       > 0) call post_data(cs%id_eb,       eb, cs%diag)

  if (cs%id_dudt_dia > 0) call post_data(cs%id_dudt_dia, adp%du_dt_dia,  cs%diag)
  if (cs%id_dvdt_dia > 0) call post_data(cs%id_dvdt_dia, adp%dv_dt_dia,  cs%diag)
  if (cs%id_wd       > 0) call post_data(cs%id_wd,       cdp%diapyc_vel, cs%diag)

  if (cs%id_MLD_003 > 0 .or. cs%id_subMLN2 > 0 .or. cs%id_mlotstsq > 0) then
    call diagnosemldbydensitydifference(cs%id_MLD_003, h, tv, 0.03, g, gv, cs%diag, &
                                        id_n2subml=cs%id_subMLN2, id_mldsq=cs%id_mlotstsq)
  endif
  if (cs%id_MLD_0125 > 0) then
    call diagnosemldbydensitydifference(cs%id_MLD_0125, h, tv, 0.125, g, gv, cs%diag)
  endif
  if (cs%id_MLD_user > 0) then
    call diagnosemldbydensitydifference(cs%id_MLD_user, h, tv, cs%MLDdensityDifference, g, gv, cs%diag)
  endif

  if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)
  if (cs%id_Tadv > 0) call post_data(cs%id_Tadv, tadv_flx, cs%diag)
  if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)
  if (cs%id_Sadv > 0) call post_data(cs%id_Sadv, sadv_flx, cs%diag)
  if (cs%use_int_tides) then
    if (cs%id_cg1 > 0) call post_data(cs%id_cg1, cn(:,:,1),cs%diag)
    do m=1,cs%nMode
      if (cs%id_cn(m) > 0) call post_data(cs%id_cn(m),cn(:,:,m),cs%diag)
    enddo
  endif

  num_z_diags = 0
  if (cs%id_Kd_z > 0) then
    num_z_diags = num_z_diags + 1
    z_ids(num_z_diags) = cs%id_Kd_z ; z_ptrs(num_z_diags)%p => kd_int
  endif
  if (cs%id_Tdif_z > 0) then
    num_z_diags = num_z_diags + 1
    z_ids(num_z_diags) = cs%id_Tdif_z ; z_ptrs(num_z_diags)%p => tdif_flx
  endif
  if (cs%id_Tadv_z > 0) then
    num_z_diags = num_z_diags + 1
    z_ids(num_z_diags) = cs%id_Tadv_z ; z_ptrs(num_z_diags)%p => tadv_flx
  endif
  if (cs%id_Sdif_z > 0) then
    num_z_diags = num_z_diags + 1
    z_ids(num_z_diags) = cs%id_Sdif_z ; z_ptrs(num_z_diags)%p => sdif_flx
  endif
  if (cs%id_Sadv_z > 0) then
    num_z_diags = num_z_diags + 1
    z_ids(num_z_diags) = cs%id_Sadv_z ; z_ptrs(num_z_diags)%p => sadv_flx
  endif

  if (num_z_diags > 0) &
    call calc_zint_diags(h, z_ptrs, z_ids, num_z_diags, g, gv, cs%diag_to_Z_CSp)

  if (cs%debugConservation) call mom_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, g)
  if (showcalltree) call calltree_leave("diabatic()")

Here is the call graph for this function:

diabatic_driver_end()

subroutine, public mom_diabatic_driver::diabatic_driver_end

(

type(diabatic_cs), pointer

CS

)

Routine to close the diabatic driver module.

Parameters

cs	module control structure

Definition at line 2314 of file MOM_diabatic_driver.F90.

References mom_diffconvection::diffconvection_end(), mom_entrain_diffusive::entrain_diffusive_end(), and mom_opacity::opacity_end().

  type(diabatic_cs), pointer :: cs    !< module control structure

  if (.not.associated(cs)) return

  call diabatic_aux_end(cs%diabatic_aux_CSp)

  call entrain_diffusive_end(cs%entrain_diffusive_CSp)
  call set_diffusivity_end(cs%set_diff_CSp)
  if (cs%useKPP) then
    deallocate( cs%KPP_buoy_flux )
    deallocate( cs%KPP_temp_flux )
    deallocate( cs%KPP_salt_flux )
  endif
  if (cs%useKPP) then
    deallocate( cs%KPP_NLTheat )
    deallocate( cs%KPP_NLTscalar )
    call kpp_end(cs%KPP_CSp)
  endif
  if (cs%useConvection) call diffconvection_end(cs%Conv_CSp)
  if (cs%use_energetic_PBL) &
    call energetic_pbl_end(cs%energetic_PBL_CSp)
  if (cs%debug_energy_req) &
    call diapyc_energy_req_end(cs%diapyc_en_rec_CSp)

  if (associated(cs%optics)) then
    call opacity_end(cs%opacity_CSp, cs%optics)
    deallocate(cs%optics)
  endif
  if (associated(cs)) deallocate(cs)

Here is the call graph for this function:

diabatic_driver_init()

subroutine, public mom_diabatic_driver::diabatic_driver_init	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		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,
		type(diag_to_z_cs), pointer	diag_to_Z_CSp
	)

This routine initializes the diabatic driver module.

Parameters

[in]	time	model time
[in,out]	g	model grid structure
[in]	gv	model vertical grid structure
[in]	param_file	file to parse for parameter values
[in]	usealealgorithm	logical for whether to use ALE remapping
[in,out]	diag	structure to regulate diagnostic output
[in,out]	adp	pointers to accelerations in momentum equations, to enable diagnostics, like energy budgets
[in,out]	cdp	pointers to terms in continuity equations
	cs	module control structure
	tracer_flow_csp	pointer to control structure of tracer flow control module
	sponge_csp	pointer to the sponge module control structure
	ale_sponge_csp	pointer to the ALE sponge module control structure
	diag_to_z_csp	pointer to the Z-diagnostics control structure

Definition at line 1789 of file MOM_diabatic_driver.F90.

References mom_bulk_mixed_layer::bulkmixedlayer_init(), mom_cvmix_shear::cvmix_shear_is_used(), mom_geothermal::geothermal_init(), id_clock_differential_diff, id_clock_entrain, id_clock_geothermal, id_clock_kpp, id_clock_mixedlayer, id_clock_pass, id_clock_remap, id_clock_set_diffusivity, id_clock_sponge, id_clock_tracers, id_clock_tridiag, mom_internal_tides::internal_tides_init(), mom_kappa_shear::kappa_shear_is_used(), mom_opacity::opacity_init(), mom_diag_to_z::register_zint_diag(), mom_regularize_layers::regularize_layers_init(), and mom_io::var_desc().

  type(time_type),         intent(in)    :: time             !< model time
  type(ocean_grid_type),   intent(inout) :: g                !< model grid structure
  type(verticalgrid_type), intent(in)    :: gv               !< model vertical grid structure
  type(param_file_type),   intent(in)    :: param_file       !< file to parse for parameter values
  logical,                 intent(in)    :: usealealgorithm  !< logical for whether to use ALE remapping
  type(diag_ctrl), target, intent(inout) :: diag             !< structure to regulate diagnostic output
  type(accel_diag_ptrs),   intent(inout) :: adp              !< pointers to accelerations in momentum equations,
                                                             !! to enable diagnostics, like energy budgets
  type(cont_diag_ptrs),    intent(inout) :: cdp              !< pointers to terms in continuity equations
  type(diabatic_cs),       pointer       :: cs               !< module control structure
  type(tracer_flow_control_cs), pointer  :: tracer_flow_csp  !< pointer to control structure of tracer flow control module
  type(sponge_cs),         pointer       :: sponge_csp       !< pointer to the sponge module control structure
  type(ale_sponge_cs),     pointer       :: ale_sponge_csp   !< pointer to the ALE sponge module control structure
  type(diag_to_z_cs),      pointer       :: diag_to_z_csp    !< pointer to the Z-diagnostics control structure

  real    :: kd
  integer :: num_mode
  logical :: use_temperature, differentialdiffusion
  type(vardesc) :: vd

! This "include" declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mod  = "MOM_diabatic_driver" ! This module's name.
  character(len=48)  :: thickness_units
  character(len=40)  :: var_name
  character(len=160) :: var_descript
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nbands, m
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (associated(cs)) then
    call mom_error(warning, "diabatic_driver_init called with an "// &
                            "associated control structure.")
    return
  else
    allocate(cs)
  endif

  cs%diag => diag
  if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
  if (associated(sponge_csp))      cs%sponge_CSp      => sponge_csp
  if (associated(ale_sponge_csp))  cs%ALE_sponge_CSp  => ale_sponge_csp
  if (associated(diag_to_z_csp))   cs%diag_to_Z_CSp   => diag_to_z_csp

  cs%useALEalgorithm = usealealgorithm
  cs%bulkmixedlayer = (gv%nkml > 0)

  ! Set default, read and log parameters
  call log_version(param_file, mod, version, &
                   "The following parameters are used for diabatic processes.")

  call get_param(param_file, mod, "SPONGE", cs%use_sponge, &
                 "If true, sponges may be applied anywhere in the domain. \n"//&
                 "The exact location and properties of those sponges are \n"//&
                 "specified via calls to initialize_sponge and possibly \n"//&
                 "set_up_sponge_field.", default=.false.)
  call get_param(param_file, mod, "ENABLE_THERMODYNAMICS", use_temperature, &
                 "If true, temperature and salinity are used as state \n"//&
                 "variables.", default=.true.)
  call get_param(param_file, mod, "ENERGETICS_SFC_PBL", cs%use_energetic_PBL, &
                 "If true, use an implied energetics planetary boundary \n"//&
                 "layer scheme to determine the diffusivity and viscosity \n"//&
                 "in the surface boundary layer.", default=.false.)
  call get_param(param_file, mod, "EPBL_IS_ADDITIVE", cs%ePBL_is_additive, &
                 "If true, the diffusivity from ePBL is added to all\n"//&
                 "other diffusivities. Otherwise, the larger of kappa-\n"//&
                 "shear and ePBL diffusivities are used.", default=.true.)
  call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialdiffusion, &
                 "If true, apply parameterization of double-diffusion.", &
                 default=.false. )
  cs%use_kappa_shear = kappa_shear_is_used(param_file)
  cs%use_CVMix_shear = cvmix_shear_is_used(param_file)
  if (cs%bulkmixedlayer) then
    call get_param(param_file, mod, "ML_MIX_FIRST", cs%ML_mix_first, &
                 "The fraction of the mixed layer mixing that is applied \n"//&
                 "before interior diapycnal mixing.  0 by default.", &
                 units="nondim", default=0.0)
    call get_param(param_file, mod, "NKBL", cs%nkbl, default=2, do_not_log=.true.)
  else
    cs%ML_mix_first = 0.0
  endif
  if (use_temperature) then
    call get_param(param_file, mod, "DO_GEOTHERMAL", cs%use_geothermal, &
                 "If true, apply geothermal heating.", default=.false.)
  else
    cs%use_geothermal = .false.
  endif
  call get_param(param_file, mod, "INTERNAL_TIDES", cs%use_int_tides, &
                 "If true, use the code that advances a separate set of \n"//&
                 "equations for the internal tide energy density.", default=.false.)
  cs%nMode = 1
  if (cs%use_int_tides) then
    ! SET NUMBER OF MODES TO CONSIDER
    call get_param(param_file, mod, "INTERNAL_TIDE_MODES", cs%nMode, &
                 "The number of distinct internal tide modes \n"//&
                 "that will be calculated.", default=1, do_not_log=.true.)

    ! The following parameters are used in testing the internal tide code.
    ! GET LOCATION AND DURATION OF ENERGY POINT SOURCE FOR TESTING (BDM)
    call get_param(param_file, mod, "INTERNAL_TIDE_SOURCE_TEST", cs%int_tide_source_test, &
                 "If true, apply an arbitrary generation site for internal tide testing", &
                 default=.false.)
    if(cs%int_tide_source_test)then
      call get_param(param_file, mod, "INTERNAL_TIDE_SOURCE_X", cs%int_tide_source_x, &
                 "X Location of generation site for internal tide", default=1.)
      call get_param(param_file, mod, "INTERNAL_TIDE_SOURCE_Y", cs%int_tide_source_y, &
                 "Y Location of generation site for internal tide", default=1.)
      call get_param(param_file, mod, "INTERNAL_TIDE_SOURCE_TLEN_DAYS", cs%tlen_days, &
                 "Time interval from start of experiment for adding wave source", &
                 units="days", default=0)
      cs%time_max_source = increment_time(time,0,days=cs%tlen_days)
    endif
    ! GET UNIFORM MODE VELOCITY FOR TESTING (BDM)
    call get_param(param_file, mod, "UNIFORM_CG", cs%uniform_cg, &
                 "If true, set cg = cg_test everywhere for test case", default=.false.)
    if(cs%uniform_cg)then
      call get_param(param_file, mod, "CG_TEST", cs%cg_test, &
                 "Uniform group velocity of internal tide for test case", default=1.)
    endif
  endif

  call get_param(param_file, mod, "MASSLESS_MATCH_TARGETS", &
                                cs%massless_match_targets, &
                 "If true, the temperature and salinity of massless layers \n"//&
                 "are kept consistent with their target densities. \n"//&
                 "Otherwise the properties of massless layers evolve \n"//&
                 "diffusively to match massive neighboring layers.", &
                 default=.true.)

  call get_param(param_file, mod, "AGGREGATE_FW_FORCING", cs%aggregate_FW_forcing, &
                 "If true, the net incoming and outgoing fresh water fluxes are combined\n"//&
                 "and applied as either incoming or outgoing depending on the sign of the net.\n"//&
                 "If false, the net incoming fresh water flux is added to the model and\n"//&
                 "thereafter the net outgoing is removed from the updated state."//&
                 "into the first non-vanished layer for which the column remains stable", &
                 default=.true.)

  call get_param(param_file, mod, "DEBUG", cs%debug, &
                 "If true, write out verbose debugging data.", default=.false.)
  call get_param(param_file, mod, "DEBUG_CONSERVATION", cs%debugConservation, &
                 "If true, monitor conservation and extrema.", default=.false.)

  call get_param(param_file, mod, "DEBUG_ENERGY_REQ", cs%debug_energy_req, &
                 "If true, debug the energy requirements.", default=.false., do_not_log=.true.)
  call get_param(param_file, mod, "MIX_BOUNDARY_TRACERS", cs%mix_boundary_tracers, &
                 "If true, mix the passive tracers in massless layers at \n"//&
                 "the bottom into the interior as though a diffusivity of \n"//&
                 "KD_MIN_TR were operating.", default=.true.)

  if (cs%mix_boundary_tracers) then
    call get_param(param_file, mod, "KD", kd, fail_if_missing=.true.)
    call get_param(param_file, mod, "KD_MIN_TR", cs%Kd_min_tr, &
                 "A minimal diffusivity that should always be applied to \n"//&
                 "tracers, especially in massless layers near the bottom. \n"//&
                 "The default is 0.1*KD.", units="m2 s-1", default=0.1*kd)
    call get_param(param_file, mod, "KD_BBL_TR", cs%Kd_BBL_tr, &
                 "A bottom boundary layer tracer diffusivity that will \n"//&
                 "allow for explicitly specified bottom fluxes. The \n"//&
                 "entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) \n"//&
                 "over the same distance.", units="m2 s-1", default=0.)
  endif

  call get_param(param_file, mod, "TRACER_TRIDIAG", cs%tracer_tridiag, &
                 "If true, use the passive tracer tridiagonal solver for T and S\n", &
                 default=.false.)

  call get_param(param_file, mod, "MINIMUM_FORCING_DEPTH", cs%minimum_forcing_depth, &
                 "The smallest depth over which forcing can be applied. This\n"//&
                 "only takes effect when near-surface layers become thin\n"//&
                 "relative to this scale, in which case the forcing tendencies\n"//&
                 "scaled down by distributing the forcing over this depth scale.", &
                 units="m", default=0.001)
  call get_param(param_file, mod, "EVAP_CFL_LIMIT", cs%evap_CFL_limit, &
                 "The largest fraction of a layer than can be lost to forcing\n"//&
                 "(e.g. evaporation, sea-ice formation) in one time-step. The unused\n"//&
                 "mass loss is passed down through the column.", &
                 units="nondim", default=0.8)


  ! Register all available diagnostics for this module.
  if (gv%Boussinesq) then ; thickness_units = "meter"
  else ; thickness_units = "kilogram meter-2" ; endif

  cs%id_ea = register_diag_field('ocean_model','ea',diag%axesTL,time, &
      'Layer entrainment from above per timestep','meter')
  cs%id_eb = register_diag_field('ocean_model','eb',diag%axesTL,time, &
      'Layer entrainment from below per timestep', 'meter')
  cs%id_dudt_dia = register_diag_field('ocean_model','dudt_dia',diag%axesCuL,time, &
      'Zonal Acceleration from Diapycnal Mixing', 'meter second-2')
  cs%id_dvdt_dia = register_diag_field('ocean_model','dvdt_dia',diag%axesCvL,time, &
      'Meridional Acceleration from Diapycnal Mixing', 'meter second-2')
  cs%id_wd = register_diag_field('ocean_model','wd',diag%axesTi,time, &
      'Diapycnal Velocity', 'meter second-1')
  if (cs%use_int_tides) then
    cs%id_cg1 = register_diag_field('ocean_model','cn1', diag%axesT1, &
                 time, 'First baroclinic mode (eigen) speed', 'm s-1')
    allocate(cs%id_cn(cs%nMode)) ; cs%id_cn(:) = -1
    do m=1,cs%nMode
      write(var_name, '("cn_mode",i1)') m
      write(var_descript, '("Baroclinic (eigen) speed of mode ",i1)') m
      cs%id_cn(m) = register_diag_field('ocean_model',var_name, diag%axesT1, &
                   time, var_descript, 'm s-1')
      call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
    enddo
  endif

  cs%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff",diag%axesTi, &
      time, "Diffusive diapycnal temperature flux across interfaces", &
      "degC meter second-1")
  cs%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv",diag%axesTi, &
      time, "Advective diapycnal temperature flux across interfaces", &
      "degC meter second-1")
  cs%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff",diag%axesTi, &
      time, "Diffusive diapycnal salnity flux across interfaces", &
      "PSU meter second-1")
  cs%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv",diag%axesTi, &
      time, "Advective diapycnal salnity flux across interfaces", &
      "PSU meter second-1")
  cs%id_MLD_003 = register_diag_field('ocean_model','MLD_003',diag%axesT1,time,        &
      'Mixed layer depth (delta rho = 0.03)', 'meter', cmor_field_name='mlotst',       &
      cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', cmor_units='m', &
      cmor_standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t')
  cs%id_mlotstsq = register_diag_field('ocean_model','mlotstsq',diag%axesT1,time,      &
      long_name='Square of Ocean Mixed Layer Thickness Defined by Sigma T',            &
      standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t',units='m2')
  cs%id_MLD_0125 = register_diag_field('ocean_model','MLD_0125',diag%axesT1,time, &
      'Mixed layer depth (delta rho = 0.125)', 'meter')
  cs%id_subMLN2  = register_diag_field('ocean_model','subML_N2',diag%axesT1,time, &
      'Squared buoyancy frequency below mixed layer', 's-2')
  cs%id_MLD_user = register_diag_field('ocean_model','MLD_user',diag%axesT1,time, &
      'Mixed layer depth (used defined)', 'meter')
  call get_param(param_file, mod, "DIAG_MLD_DENSITY_DIFF", cs%MLDdensityDifference, &
                 "The density difference used to determine a diagnostic mixed\n"//&
                 "layer depth, MLD_user, following the definition of Levitus 1982. \n"//&
                 "The MLD is the depth at which the density is larger than the\n"//&
                 "surface density by the specified amount.", units='kg/m3', default=0.1)

  ! diagnostics making use of the z-gridding code
  if (associated(diag_to_z_csp)) then
    vd = var_desc("Kd_interface", "meter2 second-1", &
                  "Diapycnal diffusivity at interfaces, interpolated to z", z_grid='z')
    cs%id_Kd_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
    vd = var_desc("Tflx_dia_diff", "degC meter second-1", &
                  "Diffusive diapycnal temperature flux across interfaces, interpolated to z", &
                  z_grid='z')
    cs%id_Tdif_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
    vd = var_desc("Tflx_dia_adv", "degC meter second-1", &
                  "Advective diapycnal temperature flux across interfaces, interpolated to z",&
                  z_grid='z')
    cs%id_Tadv_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
    vd = var_desc("Sflx_dia_diff", "PSU meter second-1", &
                  "Diffusive diapycnal salinity flux across interfaces, interpolated to z",&
                  z_grid='z')
    cs%id_Sdif_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
    vd = var_desc("Sflx_dia_adv", "PSU meter second-1", &
                  "Advective diapycnal salinity flux across interfaces, interpolated to z",&
                  z_grid='z')
    cs%id_Sadv_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
  endif

  if (cs%id_dudt_dia > 0) call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
  if (cs%id_dvdt_dia > 0) call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
  if (cs%id_wd > 0)       call safe_alloc_ptr(cdp%diapyc_vel,isd,ied,jsd,jed,nz+1)

  !call set_diffusivity_init(Time, G, param_file, diag, CS%set_diff_CSp, diag_to_Z_CSp, CS%int_tide_CSp)
  cs%id_Kd_interface = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, time, &
      'Total diapycnal diffusivity at interfaces', 'meter2 second-1')
  if (cs%use_energetic_PBL) then
      cs%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, time, &
          'ePBL diapycnal diffusivity at interfaces', 'meter2 second-1')
  endif

  cs%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, time, &
      'Total diapycnal diffusivity for heat at interfaces', 'meter2 second-1',     &
       cmor_field_name='difvho', cmor_units='m2 s-1',                              &
       cmor_standard_name='ocean_vertical_heat_diffusivity',                       &
       cmor_long_name='Ocean vertical heat diffusivity')
  cs%id_Kd_salt = register_diag_field('ocean_model', 'Kd_salt', diag%axesTi, time, &
      'Total diapycnal diffusivity for salt at interfaces', 'meter2 second-1',     &
       cmor_field_name='difvso', cmor_units='m2 s-1',                              &
       cmor_standard_name='ocean_vertical_salt_diffusivity',                       &
       cmor_long_name='Ocean vertical salt diffusivity')

  ! CS%useKPP is set to True if KPP-scheme is to be used, False otherwise.
  ! KPP_init() allocated CS%KPP_Csp and also sets CS%KPPisPassive
  cs%useKPP = kpp_init(param_file, g, diag, time, cs%KPP_CSp, passive=cs%KPPisPassive)
  if (cs%useKPP) then
    allocate( cs%KPP_NLTheat(isd:ied,jsd:jed,nz+1) )   ; cs%KPP_NLTheat(:,:,:)   = 0.
    allocate( cs%KPP_NLTscalar(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_NLTscalar(:,:,:) = 0.
  endif
  if (cs%useKPP) then
    allocate( cs%KPP_buoy_flux(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_buoy_flux(:,:,:) = 0.
    allocate( cs%KPP_temp_flux(isd:ied,jsd:jed) )      ; cs%KPP_temp_flux(:,:)   = 0.
    allocate( cs%KPP_salt_flux(isd:ied,jsd:jed) )      ; cs%KPP_salt_flux(:,:)   = 0.
  endif

  call get_param(param_file, mod, "SALT_REJECT_BELOW_ML", cs%salt_reject_below_ML, &
                 "If true, place salt from brine rejection below the mixed layer,\n"// &
                 "into the first non-vanished layer for which the column remains stable", &
                 default=.false.)

  if (cs%salt_reject_below_ML) then
      cs%id_brine_lay = register_diag_field('ocean_model','brine_layer',diag%axesT1,time, &
      'Brine insertion layer','none')
  endif


  ! diagnostics for tendencies of temp and saln due to diabatic processes;
  ! available only for ALE algorithm.
  if (cs%useALEalgorithm) then
    cs%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', &
        'diabatic_diff_temp_tendency', diag%axesTL, time,              &
        'Diabatic diffusion temperature tendency', 'Degree C per second')
    if (cs%id_diabatic_diff_temp_tend > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    cs%id_diabatic_diff_saln_tend = register_diag_field('ocean_model',&
        'diabatic_diff_saln_tendency', diag%axesTL, time,             &
        'Diabatic diffusion salinity tendency', 'PPT per second')
    if (cs%id_diabatic_diff_saln_tend > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    cs%id_diabatic_diff_heat_tend = register_diag_field('ocean_model',                                                 &
        'diabatic_heat_tendency', diag%axesTL, time,                                                                   &
        'Diabatic diffusion heat tendency',                                                                            &
        'Watts/m2',cmor_field_name='opottempdiff', cmor_units='W m-2',                                                 &
        cmor_standard_name=                                                                                            &
        'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_due_to_parameterized_dianeutral_mixing',&
        cmor_long_name =                                                                                               &
        'Tendency of sea water potential temperature expressed as heat content due to parameterized dianeutral mixing',&
        v_extensive=.true.)
    if (cs%id_diabatic_diff_heat_tend > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    cs%id_diabatic_diff_salt_tend = register_diag_field('ocean_model',                                     &
        'diabatic_salt_tendency', diag%axesTL, time,                                                       &
        'Diabatic diffusion of salt tendency',                                                             &
        'kg/(m2 * s)',cmor_field_name='osaltdiff', cmor_units='kg m-2 s-1',                                &
        cmor_standard_name=                                                                                &
        'tendency_of_sea_water_salinity_expressed_as_salt_content_due_to_parameterized_dianeutral_mixing', &
        cmor_long_name =                                                                                   &
        'Tendency of sea water salinity expressed as salt content due to parameterized dianeutral mixing', &
        v_extensive=.true.)
    if (cs%id_diabatic_diff_salt_tend > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    ! This diagnostic should equal to roundoff if all is working well.
    cs%id_diabatic_diff_heat_tend_2d = register_diag_field('ocean_model',                                                               &
        'diabatic_heat_tendency_2d', diag%axesT1, time,                                                                                 &
        'Depth integrated diabatic diffusion heat tendency',                                                                            &
        'Watts/m2',cmor_field_name='opottempdiff_2d', cmor_units='W m-2',                                                               &
        cmor_standard_name=                                                                                                             &
        'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_due_to_parameterized_dianeutral_mixing_depth_integrated',&
        cmor_long_name =                                                                                                                &
        'Tendency of sea water potential temperature expressed as heat content due to parameterized dianeutral mixing depth integrated')
    if (cs%id_diabatic_diff_heat_tend_2d > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    ! This diagnostic should equal to roundoff if all is working well.
    cs%id_diabatic_diff_salt_tend_2d = register_diag_field('ocean_model',                                                  &
        'diabatic_salt_tendency_2d', diag%axesT1, time,                                                                    &
        'Depth integrated diabatic diffusion salt tendency',                                                               &
        'kg/(m2 * s)',cmor_field_name='osaltdiff_2d', cmor_units='kg m-2 s-1',                                             &
        cmor_standard_name=                                                                                                &
        'tendency_of_sea_water_salinity_expressed_as_salt_content_due_to_parameterized_dianeutral_mixing_depth_integrated',&
        cmor_long_name =                                                                                                   &
        'Tendency of sea water salinity expressed as salt content due to parameterized dianeutral mixing depth integrated')
    if (cs%id_diabatic_diff_salt_tend_2d > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    ! diagnostics for tendencies of temp and saln due to boundary forcing;
    ! available only for ALE algorithm.
    cs%id_boundary_forcing_temp_tend = register_diag_field('ocean_model',&
        'boundary_forcing_temp_tendency', diag%axesTL, time,             &
        'Boundary forcing temperature tendency', 'Degree C per second')
    if (cs%id_boundary_forcing_temp_tend > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    cs%id_boundary_forcing_saln_tend = register_diag_field('ocean_model',&
        'boundary_forcing_saln_tendency', diag%axesTL, time,             &
        'Boundary forcing saln tendency', 'PPT per second')
    if (cs%id_boundary_forcing_saln_tend > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    cs%id_boundary_forcing_heat_tend = register_diag_field('ocean_model',&
        'boundary_forcing_heat_tendency', diag%axesTL, time,             &
        'Boundary forcing heat tendency','Watts/m2')
    if (cs%id_boundary_forcing_heat_tend > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    cs%id_boundary_forcing_salt_tend = register_diag_field('ocean_model',&
        'boundary_forcing_salt_tendency', diag%axesTL, time,             &
        'Boundary forcing salt tendency','kg m-2 s-1')
    if (cs%id_boundary_forcing_salt_tend > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    ! This diagnostic should equal to surface heat flux if all is working well.
    cs%id_boundary_forcing_heat_tend_2d = register_diag_field('ocean_model',&
        'boundary_forcing_heat_tendency_2d', diag%axesT1, time,             &
        'Depth integrated boundary forcing of ocean heat','Watts/m2')
    if (cs%id_boundary_forcing_heat_tend_2d > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    ! This diagnostic should equal to surface salt flux if all is working well.
    cs%id_boundary_forcing_salt_tend_2d = register_diag_field('ocean_model',&
        'boundary_forcing_salt_tendency_2d', diag%axesT1, time,             &
        'Depth integrated boundary forcing of ocean salt','kg m-2 s-1')
    if (cs%id_boundary_forcing_salt_tend_2d > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif
  endif

  ! diagnostics for tendencies of temp and heat due to frazil

  ! diagnostic for tendency of temp due to frazil
  cs%id_frazil_temp_tend = register_diag_field('ocean_model',&
      'frazil_temp_tendency', diag%axesTL, time,             &
      'Temperature tendency due to frazil formation', 'Degree C per second')
  if (cs%id_frazil_temp_tend > 0) then
    cs%frazil_tendency_diag = .true.
  endif

  ! diagnostic for tendency of heat due to frazil
  cs%id_frazil_heat_tend = register_diag_field('ocean_model',&
      'frazil_heat_tendency', diag%axesTL, time,             &
      'Heat tendency due to frazil formation','Watts/m2')
  if (cs%id_frazil_heat_tend > 0) then
    cs%frazil_tendency_diag = .true.
  endif

  ! if all is working propertly, this diagnostic should equal to hfsifrazil
  cs%id_frazil_heat_tend_2d = register_diag_field('ocean_model',&
      'frazil_heat_tendency_2d', diag%axesT1, time,             &
      'Depth integrated heat tendency due to frazil formation','Watts/m2')
  if (cs%id_frazil_heat_tend_2d > 0) then
    cs%frazil_tendency_diag = .true.
  endif

  if (cs%frazil_tendency_diag) then
    allocate(cs%frazil_temp_diag(isd:ied,jsd:jed,nz) ) ; cs%frazil_temp_diag(:,:,:) = 0.
    allocate(cs%frazil_heat_diag(isd:ied,jsd:jed,nz) ) ; cs%frazil_heat_diag(:,:,:) = 0.
  endif


  ! CS%useConvection is set to True IF convection will be used, otherwise False.
  ! CS%Conv_CSp is allocated by diffConvection_init()
  cs%useConvection = diffconvection_init(param_file, g, diag, time, cs%Conv_CSp)

  call entrain_diffusive_init(time, g, gv, param_file, diag, cs%entrain_diffusive_CSp)

  ! initialize the geothermal heating module
  if (cs%use_geothermal) &
    call geothermal_init(time, g, param_file, diag, cs%geothermal_CSp)

  ! initialize module for internal tide induced mixing
  if (cs%use_int_tides) then
    call int_tide_input_init(time, g, gv, param_file, diag, cs%int_tide_input_CSp, &
                             cs%int_tide_input)
    call internal_tides_init(time, g, gv, param_file, diag, cs%int_tide_CSp)
  endif

  ! initialize module for setting diffusivities
  call set_diffusivity_init(time, g, gv, param_file, diag, cs%set_diff_CSp, diag_to_z_csp, cs%int_tide_CSp)


  ! set up the clocks for this module
  id_clock_entrain = cpu_clock_id('(Ocean diabatic entrain)', grain=clock_module)
  if (cs%bulkmixedlayer) &
    id_clock_mixedlayer = cpu_clock_id('(Ocean mixed layer)', grain=clock_module)
  id_clock_remap = cpu_clock_id('(Ocean vert remap)', grain=clock_module)
  if (cs%use_geothermal) &
    id_clock_geothermal = cpu_clock_id('(Ocean geothermal)', grain=clock_routine)
  id_clock_set_diffusivity = cpu_clock_id('(Ocean set_diffusivity)', grain=clock_module)
  id_clock_kpp = cpu_clock_id('(Ocean KPP)', grain=clock_module)
  id_clock_tracers = cpu_clock_id('(Ocean tracer_columns)', grain=clock_module_driver+5)
  if (cs%use_sponge) &
    id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=clock_module)
  id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=clock_routine)
  id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=clock_routine)
  id_clock_differential_diff = -1 ; if (differentialdiffusion) &
    id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=clock_routine)

  ! initialize the auxiliary diabatic driver module
  call diabatic_aux_init(time, g, gv, param_file, diag, cs%diabatic_aux_CSp, &
                         cs%useALEalgorithm, cs%use_energetic_PBL)

  ! initialize the boundary layer modules
  if (cs%bulkmixedlayer) &
    call bulkmixedlayer_init(time, g, gv, param_file, diag, cs%bulkmixedlayer_CSp)
  if (cs%use_energetic_PBL) &
    call energetic_pbl_init(time, g, gv, param_file, diag, cs%energetic_PBL_CSp)

  call regularize_layers_init(time, g, param_file, diag, cs%regularize_layers_CSp)

  if (cs%debug_energy_req) &
    call diapyc_energy_req_init(time, g, param_file, diag, cs%diapyc_en_rec_CSp)

  ! obtain information about the number of bands for penetrative shortwave
  if (use_temperature) then
    call get_param(param_file, mod, "PEN_SW_NBANDS", nbands, default=1)
    if (nbands > 0) then
      allocate(cs%optics)
      call opacity_init(time, g, param_file, diag, cs%tracer_flow_CSp, cs%opacity_CSp, cs%optics)
    endif
  endif
  cs%nsw = 0
  if (ASSOCIATED(cs%optics)) cs%nsw = cs%optics%nbands

Here is the call graph for this function:

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(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]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	tv	points to updated thermodynamic fields
[in]	h	thickness after boundary flux application (m or kg/m2)
[in]	temp_old	temperature prior to boundary flux application
[in]	saln_old	salinity prior to boundary flux application (PPT)
[in]	h_old	thickness prior to boundary flux application (m or kg/m2)
[in]	dt	time step (sec)
	cs	module control structure

Definition at line 1608 of file MOM_diabatic_driver.F90.

Referenced by diabatic().

  type(ocean_grid_type),                    intent(in) :: g        !< ocean grid structure
  type(verticalgrid_type),                  intent(in) :: gv       !< ocean vertical grid structure
  type(thermo_var_ptrs),                    intent(in) :: tv       !< points to updated thermodynamic fields
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h        !< thickness after boundary flux application (m or kg/m2)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: temp_old !< temperature prior to boundary flux application
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: saln_old !< salinity prior to boundary flux application (PPT)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old    !< thickness prior to boundary flux application (m or kg/m2)
  real,                                     intent(in) :: dt       !< time step (sec)
  type(diabatic_cs),                        pointer    :: cs       !< module control structure

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
  real, dimension(SZI_(G),SZJ_(G))         :: work_2d
  real    :: idt
  integer :: i, j, k, is, ie, js, je, nz

  is  = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  idt = 1/dt
  work_3d(:,:,:) = 0.0
  work_2d(:,:)   = 0.0

  ! temperature tendency
  if(cs%id_boundary_forcing_temp_tend > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt
    enddo ; enddo ; enddo
    call post_data(cs%id_boundary_forcing_temp_tend, work_3d, cs%diag)
  endif

  ! heat tendency
  if(cs%id_boundary_forcing_heat_tend > 0 .or. cs%id_boundary_forcing_heat_tend_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = gv%H_to_kg_m2 * tv%C_p * idt * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * temp_old(i,j,k))
    enddo ; enddo ; enddo
    if(cs%id_boundary_forcing_heat_tend > 0) then
      call post_data(cs%id_boundary_forcing_heat_tend, work_3d, cs%diag)
    endif
    if(cs%id_boundary_forcing_heat_tend_2d > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = 0.0
        do k=1,nz
          work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_boundary_forcing_heat_tend_2d, work_2d, cs%diag)
    endif
  endif


  ! salinity tendency
  if(cs%id_boundary_forcing_saln_tend > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*idt
    enddo ; enddo ; enddo
    call post_data(cs%id_boundary_forcing_saln_tend, work_3d, cs%diag)
  endif

  ! salt tendency
  if(cs%id_boundary_forcing_salt_tend > 0 .or. cs%id_boundary_forcing_salt_tend_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = gv%H_to_kg_m2 * cs%ppt2mks * idt * (h(i,j,k) * tv%S(i,j,k) - h_old(i,j,k) * saln_old(i,j,k))
    enddo ; enddo ; enddo
    if(cs%id_boundary_forcing_salt_tend > 0) then
      call post_data(cs%id_boundary_forcing_salt_tend, work_3d, cs%diag)
    endif
    if(cs%id_boundary_forcing_salt_tend_2d > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = 0.0
        do k=1,nz
          work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_boundary_forcing_salt_tend_2d, work_2d, cs%diag)
    endif
  endif

Here is the caller graph for this function:

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(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]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	tv	points to updated thermodynamic fields
[in]	h	thickness (m or kg/m2)
[in]	temp_old	temperature prior to diabatic physics
[in]	saln_old	salinity prior to diabatic physics (PPT)
[in]	dt	time step (sec)
	cs	module control structure

Definition at line 1525 of file MOM_diabatic_driver.F90.

Referenced by diabatic().

  type(ocean_grid_type),                     intent(in) :: g        !< ocean grid structure
  type(verticalgrid_type),                   intent(in) :: gv       !< ocean vertical grid structure
  type(thermo_var_ptrs),                     intent(in) :: tv       !< points to updated thermodynamic fields
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h        !< thickness (m or kg/m2)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: temp_old !< temperature prior to diabatic physics
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: saln_old !< salinity prior to diabatic physics (PPT)
  real,                                      intent(in) :: dt       !< time step (sec)
  type(diabatic_cs),                         pointer    :: cs       !< module control structure

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
  real, dimension(SZI_(G),SZJ_(G))         :: work_2d
  real    :: idt
  integer :: i, j, k, is, ie, js, je, nz

  is  = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  idt = 1/dt
  work_3d(:,:,:) = 0.0
  work_2d(:,:)   = 0.0


  ! temperature tendency
  do k=1,nz ; do j=js,je ; do i=is,ie
    work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt
  enddo ; enddo ; enddo
  if(cs%id_diabatic_diff_temp_tend > 0) then
    call post_data(cs%id_diabatic_diff_temp_tend, work_3d, cs%diag)
  endif

  ! heat tendency
  if(cs%id_diabatic_diff_heat_tend > 0 .or. cs%id_diabatic_diff_heat_tend_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = h(i,j,k) * gv%H_to_kg_m2 * tv%C_p * work_3d(i,j,k)
    enddo ; enddo ; enddo
    if(cs%id_diabatic_diff_heat_tend > 0) then
      call post_data(cs%id_diabatic_diff_heat_tend, work_3d, cs%diag)
    endif
    if(cs%id_diabatic_diff_heat_tend_2d > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = 0.0
        do k=1,nz
          work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_diabatic_diff_heat_tend_2d, work_2d, cs%diag)
    endif
  endif

  ! salinity tendency
  do k=1,nz ; do j=js,je ; do i=is,ie
    work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*idt
  enddo ; enddo ; enddo
  if(cs%id_diabatic_diff_saln_tend > 0) then
    call post_data(cs%id_diabatic_diff_saln_tend, work_3d, cs%diag)
  endif

  ! salt tendency
  if(cs%id_diabatic_diff_salt_tend > 0 .or. cs%id_diabatic_diff_salt_tend_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = h(i,j,k) * gv%H_to_kg_m2 * cs%ppt2mks * work_3d(i,j,k)
    enddo ; enddo ; enddo
    if(cs%id_diabatic_diff_salt_tend > 0) then
      call post_data(cs%id_diabatic_diff_salt_tend, work_3d, cs%diag)
    endif
    if(cs%id_diabatic_diff_salt_tend_2d > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = 0.0
        do k=1,nz
          work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_diabatic_diff_salt_tend_2d, work_2d, cs%diag)
    endif
  endif

Here is the caller graph for this function:

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(diabatic_cs), pointer  CS, integer, intent(in)  ncall  )

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]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	tv	points to updated thermodynamic fields
[in]	h	thickness (m or kg/m2)
[in]	temp_old	temperature prior to frazil formation
[in]	dt	time step (sec)
[in]	ncall	the first or second call of this routine
	cs	module control structure

Definition at line 1691 of file MOM_diabatic_driver.F90.

Referenced by diabatic().

  type(ocean_grid_type),                    intent(in) :: g        !< ocean grid structure
  type(verticalgrid_type),                  intent(in) :: gv       !< ocean vertical grid structure
  type(thermo_var_ptrs),                    intent(in) :: tv       !< points to updated thermodynamic fields
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h        !< thickness (m or kg/m2)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: temp_old !< temperature prior to frazil formation
  real,                                     intent(in) :: dt       !< time step (sec)
  integer,                                  intent(in) :: ncall    !< the first or second call of this routine
  type(diabatic_cs),                        pointer    :: cs       !< module control structure

  real, dimension(SZI_(G),SZJ_(G))         :: work_2d
  real    :: idt
  integer :: i, j, k, is, ie, js, je, nz

  is  = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  idt = 1/dt
  work_2d(:,:) = 0.0

  ! zero the tendencies at start of first call
  if(ncall == 1) then
    cs%frazil_heat_diag(:,:,:) = 0.0
    cs%frazil_temp_diag(:,:,:) = 0.0
  endif

  ! temperature tendency
  if(cs%id_frazil_temp_tend > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%frazil_temp_diag(i,j,k) = cs%frazil_temp_diag(i,j,k) + idt * (tv%T(i,j,k)-temp_old(i,j,k))
    enddo ; enddo ; enddo
    if(ncall == 2) then
      call post_data(cs%id_frazil_temp_tend, cs%frazil_temp_diag(:,:,:), cs%diag)
    endif
  endif

  ! heat tendency
  if(cs%id_frazil_heat_tend > 0 .or. cs%id_frazil_heat_tend_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%frazil_heat_diag(i,j,k) = cs%frazil_heat_diag(i,j,k) + &
                                   gv%H_to_kg_m2 * tv%C_p * h(i,j,k) * idt * (tv%T(i,j,k)-temp_old(i,j,k))
    enddo ; enddo ; enddo
    if(cs%id_frazil_heat_tend  > 0 .and. ncall == 2) then
      call post_data(cs%id_frazil_heat_tend, cs%frazil_heat_diag(:,:,:), cs%diag)
    endif

    ! As a consistency check, we must have
    ! FRAZIL_HEAT_TENDENCY_2d = HFSIFRAZIL
    if(cs%id_frazil_heat_tend_2d > 0 .and. ncall == 2) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = 0.0
        do k=1,nz
          work_2d(i,j) = work_2d(i,j) + cs%frazil_heat_diag(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_frazil_heat_tend_2d, work_2d, cs%diag)
    endif
  endif

Here is the caller graph for this function:

extract_diabatic_member()

subroutine, public mom_diabatic_driver::extract_diabatic_member	(	type(diabatic_cs), intent(in)	CS,
		type(opacity_cs), intent(out), optional, pointer	opacity_CSp,
		type(optics_type), intent(out), optional, pointer	optics_CSp,
		real, intent(out), optional	evap_CFL_limit,
		real, intent(out), optional	minimum_forcing_depth
	)

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

Definition at line 1484 of file MOM_diabatic_driver.F90.

  type(diabatic_cs),  intent(in   ) :: cs
  ! All output arguments are optional
  type(opacity_cs),   pointer, optional, intent(  out) :: opacity_csp
  type(optics_type),  pointer, optional, intent(  out) :: optics_csp
  real,                        optional, intent(  out) :: evap_cfl_limit
  real,                        optional, intent(  out) :: minimum_forcing_depth

  ! Pointers to control structures
  if (present(opacity_csp)) opacity_csp => cs%opacity_CSp
  if (present(optics_csp))  optics_csp  => cs%optics

  ! Constants within diabatic_CS
  if (present(evap_cfl_limit))        evap_cfl_limit = cs%evap_CFL_limit
  if (present(minimum_forcing_depth)) minimum_forcing_depth = cs%minimum_forcing_depth

Variable Documentation

id_clock_differential_diff

integer mom_diabatic_driver::id_clock_differential_diff

private

Definition at line 227 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

id_clock_entrain

integer mom_diabatic_driver::id_clock_entrain

Definition at line 225 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity

id_clock_geothermal

integer mom_diabatic_driver::id_clock_geothermal

private

Definition at line 227 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap

id_clock_kpp

integer mom_diabatic_driver::id_clock_kpp

private

Definition at line 228 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

integer :: id_clock_kpp

id_clock_mixedlayer

integer mom_diabatic_driver::id_clock_mixedlayer

private

Definition at line 225 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

id_clock_pass

integer mom_diabatic_driver::id_clock_pass

private

Definition at line 226 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

id_clock_remap

integer mom_diabatic_driver::id_clock_remap

private

Definition at line 227 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

id_clock_set_diffusivity

integer mom_diabatic_driver::id_clock_set_diffusivity

private

Definition at line 225 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

id_clock_sponge

integer mom_diabatic_driver::id_clock_sponge

private

Definition at line 226 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

id_clock_tracers

integer mom_diabatic_driver::id_clock_tracers

private

Definition at line 226 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().

integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge

id_clock_tridiag

integer mom_diabatic_driver::id_clock_tridiag

private

Definition at line 226 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), and diabatic_driver_init().