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...

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

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

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

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

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

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

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

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

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

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

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

subroutine, public	diabatic (u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, GV, US, CS, WAVES)
	This subroutine imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers. More...

subroutine	diabatic_ale_legacy (u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, GV, US, CS, WAVES)
	Applies diabatic forcing and diapycnal mixing of temperature, salinity and other tracers for use with an ALE algorithm. This version uses an older set of algorithms compared with diabatic_ALE. More...

subroutine	diabatic_ale (u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, GV, US, CS, Waves)
	This subroutine imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers. More...

subroutine	layered_diabatic (u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, GV, US, CS, WAVES)
	Imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers using the original MOM6 algorithms. More...

subroutine, public	extract_diabatic_member (CS, opacity_CSp, optics_CSp, evap_CFL_limit, minimum_forcing_depth, KPP_CSp, energetic_PBL_CSp, diabatic_aux_CSp)
	Returns pointers or values of members within the diabatic_CS type. For extensibility, each returned argument is an optional argument. More...

subroutine, public	adiabatic (h, tv, fluxes, dt, G, GV, US, CS)
	Routine called for adiabatic physics. More...

subroutine	diagnose_diabatic_diff_tendency (tv, h, temp_old, saln_old, dt, G, GV, 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)
	This routine diagnoses tendencies for temperature and heat from frazil formation. This routine is called twice from within subroutine diabatic; at start and at end of the diabatic processes. The impacts from frazil are generally a function of depth. Hence, when checking heat budget, be sure to remove HFSIFRAZIL from HFDS in k=1. More...

subroutine, public	adiabatic_driver_init (Time, G, param_file, diag, CS, tracer_flow_CSp)
	A simplified version of diabatic_driver_init that will allow tracer column functions to be called without allowing any of the diabatic processes to be used. More...

subroutine, public	diabatic_driver_init (Time, G, GV, US, param_file, useALEalgorithm, diag, ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, ALE_sponge_CSp)
	This routine initializes the diabatic driver module. More...

subroutine, public	diabatic_driver_end (CS)
	Routine to close the diabatic driver module. More...

Function/Subroutine Documentation

◆ adiabatic()

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

Routine called for adiabatic physics.

Parameters

[in,out]	g	ocean grid structure
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields
[in,out]	fluxes	boundary fluxes
[in]	dt	time step [T ~> s]
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
	cs	module control structure

Definition at line 2892 of file MOM_diabatic_driver.F90.

   type(ocean_grid_type),   intent(inout) :: G      !< ocean grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                            intent(inout) :: h      !< thickness [H ~> m or kg m-2]
   type(thermo_var_ptrs),   intent(inout) :: tv     !< points to thermodynamic fields
   type(forcing),           intent(inout) :: fluxes !< boundary fluxes
   real,                    intent(in)    :: dt     !< time step [T ~> s]
   type(verticalGrid_type), intent(in)    :: GV     !< ocean vertical grid structure
   type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
   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, us, tv, &
                               cs%optics, cs%tracer_flow_CSp, cs%debug)
  

References mom_tracer_flow_control::call_tracer_column_fns().

Referenced by mom::step_mom_thermo().

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
	)

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	pointer to control structure of the tracer flow control module

Definition at line 3155 of file MOM_diabatic_driver.F90.

   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  !< pointer to control structure of the
                                                              !! tracer flow control module
  
 ! This "include" declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mdl = "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
  
 ! Set default, read and log parameters
   call log_version(param_file, mdl, version, &
                    "The following parameters are used for diabatic processes.")
  

Referenced by mom::initialize_mom().

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(time_type), intent(in)	Time_end,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(diabatic_cs), pointer	CS,
		type(wave_parameters_cs), optional, pointer	WAVES
	)

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 259 of file MOM_diabatic_driver.F90.

   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 [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v         !< meridional velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h         !< thickness [H ~> m or kg m-2]
   type(thermo_var_ptrs),                     intent(inout) :: tv        !< points to thermodynamic fields
                                                                         !! unused have NULL ptrs
   real, dimension(:,:),                      pointer       :: Hml       !< Active mixed layer depth [m]
   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
   type(accel_diag_ptrs),                     intent(inout) :: ADp       !< related points to accelerations in momentum
                                                                         !! equations, to enable the later derived
                                                                         !! diagnostics, like energy budgets
   type(cont_diag_ptrs),                      intent(inout) :: CDp       !< points to terms in continuity equations
   real,                                      intent(in)    :: dt        !< time increment [T ~> s]
   type(time_type),                           intent(in)    :: Time_end  !< Time at the end of the interval
   type(unit_scale_type),                     intent(in)    :: US        !< A dimensional unit scaling type
   type(diabatic_CS),                         pointer       :: CS        !< module control structure
   type(Wave_parameters_CS),        optional, pointer       :: Waves     !< Surface gravity waves
  
   ! local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
     eta      ! Interface heights before diapycnal mixing [m].
   real, dimension(SZI_(G),SZJ_(G),CS%nMode) :: &
     cn_IGW   ! baroclinic internal gravity wave speeds
   real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag             ! diagnostic array for temp
   integer :: i, j, k, m, is, ie, js, je, nz
   logical :: showCallTree ! If true, show the call tree
  
   if (g%ke == 1) return
  
   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  
   if (.not. associated(cs)) call mom_error(fatal, "MOM_diabatic_driver: "// &
          "Module must be initialized before it is used.")
   if (dt == 0.0) call mom_error(fatal, "MOM_diabatic_driver: "// &
         "diabatic was called with a zero length timestep.")
   if (dt < 0.0) call mom_error(fatal, "MOM_diabatic_driver: "// &
         "diabatic was called with a negative timestep.")
  
   showcalltree = calltree_showquery()
  
   ! Offer diagnostics of various state varables at the start of diabatic
   ! these are mostly for debugging purposes.
   if (cs%id_u_predia > 0) call post_data(cs%id_u_predia, u, cs%diag)
   if (cs%id_v_predia > 0) call post_data(cs%id_v_predia, v, cs%diag)
   if (cs%id_h_predia > 0) call post_data(cs%id_h_predia, h, cs%diag)
   if (cs%id_T_predia > 0) call post_data(cs%id_T_predia, tv%T, cs%diag)
   if (cs%id_S_predia > 0) call post_data(cs%id_S_predia, tv%S, cs%diag)
   if (cs%id_e_predia > 0) then
     call find_eta(h, tv, g, gv, us, eta, eta_to_m=1.0)
     call post_data(cs%id_e_predia, eta, cs%diag)
   endif
  
   if (cs%debug) then
     call mom_state_chksum("Start of diabatic ", u, v, h, g, gv, us, haloshift=0)
     call mom_forcing_chksum("Start of diabatic", fluxes, g, us, haloshift=0)
   endif
   if (cs%debugConservation) call mom_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, g, gv, us)
  
   if (cs%debug_energy_req) &
     call diapyc_energy_req_test(h, dt, tv, g, gv, us, cs%diapyc_en_rec_CSp)
  
   call cpu_clock_begin(id_clock_set_diffusivity)
   call set_bbl_tke(u, v, h, fluxes, visc, g, gv, us, 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
     ! For frazil diagnostic, the first call covers the first half of the time step
     call enable_averages(0.5*dt, time_end - real_to_time(0.5*us%T_to_s*dt), cs%diag)
     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, halo=cs%halo_TS_diff)
     else
       call make_frazil(h, tv, g, gv, cs%diabatic_aux_CSp, halo=cs%halo_TS_diff)
     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, 0.5*dt, g, gv, cs)
       if (cs%id_frazil_h > 0) call post_data(cs%id_frazil_h, h, cs%diag)
     endif
     call disable_averaging(cs%diag)
   endif ! associated(tv%T) .AND. associated(tv%frazil)
   if (cs%debugConservation) call mom_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, g, gv, us)
  
  
   if (cs%use_int_tides) then
     ! This block provides an interface for the unresolved low-mode internal tide module (BDM).
     call set_int_tide_input(u, v, h, tv, fluxes, cs%int_tide_input, dt, g, gv, us, &
                             cs%int_tide_input_CSp)
     cn_igw(:,:,:) = 0.0
     if (cs%uniform_test_cg > 0.0) then
       do m=1,cs%nMode ; cn_igw(:,:,m) = cs%uniform_test_cg ; enddo
     else
       call wave_speeds(h, tv, g, gv, us, cs%nMode, cn_igw, full_halos=.true.)
     endif
  
     call propagate_int_tide(h, tv, cn_igw, cs%int_tide_input%TKE_itidal_input, cs%int_tide_input%tideamp, &
                             cs%int_tide_input%Nb, dt, g, gv, us, cs%int_tide_CSp)
     if (showcalltree) call calltree_waypoint("done with propagate_int_tide (diabatic)")
   endif ! end CS%use_int_tides
  
   if (cs%useALEalgorithm .and. cs%use_legacy_diabatic) then
     call diabatic_ale_legacy(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
                       g, gv, us, cs, waves)
   elseif (cs%useALEalgorithm) then
     call diabatic_ale(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
                       g, gv, us, cs, waves)
   else
     call layered_diabatic(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
                           g, gv, us, cs, waves)
   endif
  
  
   call cpu_clock_begin(id_clock_pass)
   if (associated(visc%Kv_shear)) &
     call pass_var(visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)
   call cpu_clock_end(id_clock_pass)
  
   call disable_averaging(cs%diag)
   ! 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
     call enable_averages(0.5*dt, time_end, cs%diag)
     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, 0.5*dt, g, gv, cs)
       if (cs%id_frazil_h > 0 ) call post_data(cs%id_frazil_h, h, cs%diag)
     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, gv, us)
     call disable_averaging(cs%diag)
  
   endif  ! endif for frazil
  
  
   ! Diagnose mixed layer depths.
   call enable_averages(dt, time_end, 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*us%kg_m3_to_R, g, gv, us, cs%diag, &
                                         id_n2subml=cs%id_subMLN2, id_mldsq=cs%id_mlotstsq, dz_subml=cs%dz_subML_N2)
   endif
   if (cs%id_MLD_0125 > 0) then
     call diagnosemldbydensitydifference(cs%id_MLD_0125, h, tv, 0.125*us%kg_m3_to_R, g, gv, us, cs%diag)
   endif
   if (cs%id_MLD_user > 0) then
     call diagnosemldbydensitydifference(cs%id_MLD_user, h, tv, cs%MLDdensityDifference, g, gv, us, cs%diag)
   endif
   if (cs%use_int_tides) then
     if (cs%id_cg1 > 0) call post_data(cs%id_cg1, cn_igw(:,:,1),cs%diag)
     do m=1,cs%nMode ; if (cs%id_cn(m) > 0) call post_data(cs%id_cn(m), cn_igw(:,:,m), cs%diag) ; enddo
   endif
   call disable_averaging(cs%diag)
  
   if (cs%debugConservation) call mom_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, g, gv, us)
  

References mom_error_handler::calltree_waypoint(), diabatic_ale(), diabatic_ale_legacy(), diagnose_frazil_tendency(), mom_diapyc_energy_req::diapyc_energy_req_test(), id_clock_pass, id_clock_set_diffusivity, layered_diabatic(), mom_checksum_packages::mom_state_stats(), mom_time_manager::real_to_time(), and mom_wave_speed::wave_speeds().

Here is the call graph for this function:

◆ diabatic_ale()

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

private

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

Parameters

[in,out]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	u	zonal velocity [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
	cs	module control structure
	waves	Surface gravity waves

Definition at line 1229 of file MOM_diabatic_driver.F90.

   type(ocean_grid_type),                     intent(inout) :: G         !< ocean grid structure
   type(verticalGrid_type),                   intent(in)    :: GV        !< ocean vertical grid structure
   type(unit_scale_type),                     intent(in)    :: US        !< A dimensional unit scaling type
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u         !< zonal velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v         !< meridional velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h         !< thickness [H ~> m or kg m-2]
   type(thermo_var_ptrs),                     intent(inout) :: tv        !< points to thermodynamic fields
                                                                         !! unused have NULL ptrs
   real, dimension(:,:),                      pointer       :: Hml       !< Active mixed layer depth [m]
   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
   type(accel_diag_ptrs),                     intent(inout) :: ADp       !< related points to accelerations in momentum
                                                                         !! equations, to enable the later derived
                                                                         !! diagnostics, like energy budgets
   type(cont_diag_ptrs),                      intent(inout) :: CDp       !< points to terms in continuity equations
   real,                                      intent(in)    :: dt        !< time increment [T ~> s]
   type(time_type),                           intent(in)    :: Time_end  !< Time at the end of the interval
   type(diabatic_CS),                         pointer       :: CS        !< module control structure
   type(Wave_parameters_CS),        optional, pointer       :: Waves     !< Surface gravity waves
  
   ! local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
     ea_s,     &  ! amount of fluid entrained from the layer above within
                  ! one time step  [H ~> m or kg m-2]
     eb_s,     &  ! amount of fluid entrained from the layer below within
                  ! one time step  [H ~> m or kg m-2]
     ea_t,     &  ! amount of fluid entrained from the layer above within
                  ! one time step  [H ~> m or kg m-2]
     eb_t,     &  ! amount of fluid entrained from the layer below within
                  ! one time step  [H ~> m or kg m-2]
     kd_lay, &    ! diapycnal diffusivity of layers [Z2 T-1 ~> m2 s-1]
     h_orig, &    ! initial layer thicknesses [H ~> m or kg m-2]
     h_prebound, & ! initial layer thicknesses [H ~> m or kg m-2]
 !    hold,   &    ! layer thickness before diapycnal entrainment, and later
                  ! the initial layer thicknesses (if a mixed layer is used),
                  ! [H ~> m or kg m-2]
     dsv_dt, &    ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
     dsv_ds, &    ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
     ctke,   &    ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
     u_h,    &    ! zonal and meridional velocities at thickness points after
     v_h          ! entrainment [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJ_(G)) :: &
     Rcv_ml, &   ! coordinate density of mixed layer, used for applying sponges
     SkinBuoyFlux! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], 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 :: net_ent  ! The net of ea-eb at an interface.
  
   real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: &
     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 [H ~> m or kg m-2] (for Bous or non-Bouss)
  
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), target :: &
     Kd_int,   & ! diapycnal diffusivity of interfaces [Z2 T-1 ~> m2 s-1]
     Kd_heat,  & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1]
     Kd_salt,  & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1]
     Kd_ePBL,  & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]
     Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
     Tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
     Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
     Sadv_flx    ! advective diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
  
   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
                        ! [H ~> m or kg m-2]
   real :: h_neglect    ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected
                        ! [H ~> m or kg m-2]
   real :: h_neglect2   ! h_neglect^2 [H2 ~> m2 or kg2 m-4]
   real :: add_ent      ! Entrainment that needs to be added when mixing tracers
                        ! [H ~> m or kg m-2]
   real :: eaval        ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
   real :: hval         ! hval is 2*h at velocity grid points [H ~> m or kg m-2]
   real :: h_tr         ! h_tr is h at tracer points with a tiny thickness
                        ! added to ensure positive definiteness [H ~> m or kg m-2]
   real :: Tr_ea_BBL    ! The diffusive tracer thickness in the BBL that is
                        ! coupled to the bottom within a timestep [H ~> m or kg m-2]
  
   real :: htot(SZIB_(G))             ! The summed thickness from the bottom [H ~> m or kg m-2].
   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 ~> m or kg m-2]
   real :: Idt     ! The inverse time step [T-1 ~> s-1]
  
   integer :: dir_flag     ! An integer encoding the directions in which to do halo updates.
   logical :: showCallTree ! If true, show the call tree
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo
  
   integer :: ig, jg      ! global indices for testing testing itide point source (BDM)
   real :: Kd_add_here    ! An added diffusivity [Z2 T-1 ~> m2 s-1].
  
   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
  
   showcalltree = calltree_showquery()
   if (showcalltree) call calltree_enter("diabatic_ALE(), MOM_diabatic_driver.F90")
  
   if (.not. (cs%useALEalgorithm)) call mom_error(fatal, "MOM_diabatic_driver: "// &
          "The ALE algorithm must be enabled when using MOM_diabatic_driver.")
  
   ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
   call enable_averages(dt, time_end, cs%diag)
  
   if (cs%use_geothermal) then
     halo = cs%halo_TS_diff
     !$OMP parallel do default(shared)
     do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo
       h_orig(i,j,k) = h(i,j,k) ; eatr(i,j,k) = 0.0 ; ebtr(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, eatr, ebtr, g, gv, us, cs%geothermal_CSp, halo=cs%halo_TS_diff)
     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, gv, us)
   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_pen_shortwave 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_pen_shortwave(cs%optics, fluxes, g, gv, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
  
   if (cs%debug) call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)
  
   if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then
     if (cs%use_geothermal) then
       call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, us, eatr, ebtr)
       if (cs%debug) then
         call hchksum(eatr, "after find_uv_at_h eatr",g%HI, scale=gv%H_to_m)
         call hchksum(ebtr, "after find_uv_at_h ebtr",g%HI, scale=gv%H_to_m)
       endif
     else
       call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
     endif
     if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
   endif
  
   call cpu_clock_begin(id_clock_set_diffusivity)
   ! Sets: Kd_lay, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S and visc%TKE_turb
   ! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow
   call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics,  &
                        visc, dt, g, gv, us,cs%set_diff_CSp, kd_lay, 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, us, haloshift=0)
     call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)
     call mom_thermovar_chksum("after set_diffusivity ", tv, g)
     call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
   endif
  
   ! Set diffusivities for heat and salt separately
  
   !$OMP parallel do default(shared)
   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
   ! Add contribution from double diffusion
   if (associated(visc%Kd_extra_S)) then
     !$OMP parallel do default(shared)
     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 parallel do default(shared)
     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
  
   if (cs%debug) then
     call hchksum(kd_heat, "after set_diffusivity Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     call hchksum(kd_salt, "after set_diffusivity Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
   endif
  
   if (cs%useKPP) then
     call cpu_clock_begin(id_clock_kpp)
     ! total vertical viscosity in the interior is represented via visc%Kv_shear
     do k=1,nz+1 ; do j=js,je ; do i=is,ie
       visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k)
     enddo ; enddo ; enddo
  
     ! 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, us, 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.
  
     call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv%eqn_of_state, &
                          fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
  
     call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
                        kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
  
     if (associated(hml)) then
       !$OMP parallel default(shared)
       call kpp_get_bld(cs%KPP_CSp, hml(:,:), g)
       !$OMP end parallel
       call pass_var(hml, g%domain, halo=1)
       ! If visc%MLD exists, copy KPP's BLD into it
       if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
     endif
  
     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, us, haloshift=0)
       call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)
       call mom_thermovar_chksum("after KPP", tv, g)
       call hchksum(kd_heat, "after KPP Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
       call hchksum(kd_salt, "after KPP Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     endif
  
   endif  ! endif for KPP
  
  
   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, &
                                     us%T_to_s*dt, tv%T, tv%C_p)
     call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, &
                                     us%T_to_s*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, gv, us)
  
     if (cs%debug) then
       call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)
       call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)
       call mom_thermovar_chksum("after KPP_applyNLT ", tv, g)
     endif
   endif ! endif for KPP
  
   ! This is the "old" method for applying differential diffusion.
   ! Changes: tv%T, tv%S
   if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T) .and. &
       (.not.cs%use_CVMix_ddiff)) 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, gv, us)
  
     ! increment heat and salt diffusivity.
     ! CS%useKPP==.true. already has extra_T and extra_S included
     if (.not. cs%useKPP) then
       !$OMP parallel do default(shared)
       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
  
   ! Calculate vertical mixing due to convection (computed via CVMix)
   if (cs%use_CVMix_conv) then
     call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml)
     ! Increment vertical diffusion and viscosity due to convection
     !$OMP parallel do default(shared)
     do k=1,nz+1 ; do j=js,je ; do i=is,ie
       kd_heat(i,j,k) = kd_heat(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
       kd_salt(i,j,k) = kd_salt(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
       if (cs%useKPP) then
         visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
       else
         visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
       endif
     enddo ; enddo ; enddo
   endif
  
   ! Save fields before boundary forcing is applied for tendency 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
  
   ! Apply forcing
   call cpu_clock_begin(id_clock_remap)
  
   ! Changes made to following fields:  h, tv%T and tv%S.
   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, us, dt, fluxes, cs%optics, &
             optics_nbands(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_t, "after applyBoundaryFluxes ea_t",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_t, "after applyBoundaryFluxes eb_t",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(ea_s, "after applyBoundaryFluxes ea_s",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_s, "after applyBoundaryFluxes eb_s",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(ctke, "after applyBoundaryFluxes cTKE",g%HI,haloshift=0, &
                    scale=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**2)
       call hchksum(dsv_dt, "after applyBoundaryFluxes dSV_dT",g%HI,haloshift=0, scale=us%kg_m3_to_R)
       call hchksum(dsv_ds, "after applyBoundaryFluxes dSV_dS",g%HI,haloshift=0, scale=us%kg_m3_to_R)
     endif
  
     call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
     call energetic_pbl(h, u_h, v_h, tv, fluxes, dt, kd_epbl, g, gv, us, &
                        cs%energetic_PBL_CSp, dsv_dt, dsv_ds, ctke, skinbuoyflux, waves=waves)
  
     if (associated(hml)) then
       call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hml(:,:), g, us)
       call pass_var(hml, g%domain, halo=1)
       ! If visc%MLD exists, copy ePBL's MLD into it
       if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
     elseif (associated(visc%MLD)) then
       call energetic_pbl_get_mld(cs%energetic_PBL_CSp, visc%MLD, g, us)
       call pass_var(visc%MLD, g%domain, halo=1)
     endif
  
     ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
     do k=2,nz ; do j=js,je ; do i=is,ie
       !### These expressions assume a Prandtl number of 1.
       if (cs%ePBL_is_additive) then
         kd_add_here = kd_epbl(i,j,k)
         visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + kd_epbl(i,j,k)
       else
         kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)
         visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), kd_epbl(i,j,k))
       endif
  
       kd_heat(i,j,k) = kd_heat(i,j,k) + kd_add_here
       kd_salt(i,j,k) = kd_salt(i,j,k) + kd_add_here
  
     enddo ; enddo ; enddo
  
     if (cs%debug) then
       call hchksum(ea_t, "after ePBL ea_t",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_t, "after ePBL eb_t",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(ea_s, "after ePBL ea_s",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_s, "after ePBL eb_s",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(kd_epbl, "after ePBL Kd_ePBL", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     endif
  
   else
     call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
                                   optics_nbands(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
   ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme
   !  so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards
   if (cs%boundary_forcing_tendency_diag) then
     call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, g, gv, cs)
     if (cs%id_boundary_forcing_h > 0) call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h = h_diag)
   endif
   ! Boundary fluxes may have changed T, S, and h
   call diag_update_remap_grids(cs%diag)
   call cpu_clock_end(id_clock_remap)
   if (cs%debug) then
     call mom_forcing_chksum("after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)
     call mom_thermovar_chksum("after applyBoundaryFluxes ", tv, g)
     call mom_state_chksum("after applyBoundaryFluxes ", u, v, h, g, gv, us, 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, gv, us)
  
   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, gv, us)
  
   ! calculate change in temperature & salinity due to dia-coordinate surface diffusion
   if (associated(tv%T)) then
  
     if (cs%debug) then
       call hchksum(ea_t, "before triDiagTS ea_t ",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_t, "before triDiagTS eb_t ",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(ea_s, "before triDiagTS ea_s ",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_s, "before triDiagTS eb_s ",g%HI,haloshift=0, scale=gv%H_to_m)
     endif
  
     call cpu_clock_begin(id_clock_tridiag)
     !  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)
  
     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
  
     ! set ea_t=eb_t=Kd_heat and ea_s=eb_s=Kd_salt on interfaces for use in the
     ! tri-diagonal solver.
  
     do j=js,je ; do i=is,ie
       ea_t(i,j,1) = 0.; ea_s(i,j,1) = 0.
     enddo ; enddo
  
     !$OMP parallel do default(shared) 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_t(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_heat(i,j,k)
       eb_t(i,j,k-1) = ea_t(i,j,k)
       ea_s(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_salt(i,j,k)
       eb_s(i,j,k-1) = ea_s(i,j,k)
     enddo ; enddo ; enddo
     do j=js,je ; do i=is,ie
       eb_t(i,j,nz) = 0. ; eb_s(i,j,nz) = 0.
     enddo ; enddo
     if (showcalltree) call calltree_waypoint("done setting ea_t,ea_s,eb_t,eb_s from Kd_heat" //&
        "and Kd_salt (diabatic)")
  
   ! Initialize halo regions of ea, eb, and hold to default values.
   !$OMP parallel do default(shared)
     do k=1,nz
       do i=is-1,ie+1
         ea_t(i,js-1,k) = 0.0 ; eb_t(i,js-1,k) = 0.0
         ea_s(i,js-1,k) = 0.0 ; eb_s(i,js-1,k) = 0.0
         ea_t(i,je+1,k) = 0.0 ; eb_t(i,je+1,k) = 0.0
         ea_s(i,je+1,k) = 0.0 ; eb_s(i,je+1,k) = 0.0
       enddo
       do j=js,je
         ea_t(is-1,j,k) = 0.0 ; eb_t(is-1,j,k) = 0.0
         ea_s(is-1,j,k) = 0.0 ; eb_s(is-1,j,k) = 0.0
         ea_t(ie+1,j,k) = 0.0 ; eb_t(ie+1,j,k) = 0.0
         ea_s(ie+1,j,k) = 0.0 ; eb_s(ie+1,j,k) = 0.0
       enddo
     enddo
  
   ! Changes T and S via the tridiagonal solver; no change to h
     call tracer_vertdiff(h, ea_t, eb_t, dt, tv%T, g, gv)
     call tracer_vertdiff(h, ea_s, eb_s, dt, tv%S, g, gv)
  
  
     ! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below
     if (cs%diabatic_diff_tendency_diag) then
       call diagnose_diabatic_diff_tendency(tv, h, 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, gv, us)
  
   if (cs%debug) then
     call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)
     call mom_thermovar_chksum("after mixed layer ", tv, 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
   idt = 1.0 / dt
   if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 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(shared)
     do k=2,nz ; do j=js,je ; do i=is,ie
       tdif_flx(i,j,k) = (idt * 0.5*(ea_t(i,j,k) + eb_t(i,j,k-1))) * &
                         (tv%T(i,j,k-1) - tv%T(i,j,k))
       tadv_flx(i,j,k) = (idt * (ea_t(i,j,k) - eb_t(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_Sadv > 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(shared)
     do k=2,nz ; do j=js,je ; do i=is,ie
       sdif_flx(i,j,k) = (idt * 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))) * &
                         (tv%S(i,j,k-1) - tv%S(i,j,k))
       sadv_flx(i,j,k) = (idt * (ea_s(i,j,k) - eb_s(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 = gv%Z_to_H * sqrt(dt*cs%Kd_BBL_tr)
     !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
     do j=js,je
       do i=is,ie
         ebtr(i,j,nz) = eb_s(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_int,
           ! perhaps a molecular diffusivity.
           add_ent = ((dt * cs%Kd_min_tr) * gv%Z_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_s(i,j,k) + eb_s(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_s(i,j,k),eb_s(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_s(i,j,k-1) + add_ent
           eatr(i,j,k) = ea_s(i,j,k) + add_ent
         else
           ebtr(i,j,k-1) = eb_s(i,j,k-1) ; eatr(i,j,k) = ea_s(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%Z_to_H**2) / &
              (0.5 * (h(i,j,k-1) + h(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_s(i,j,1) ; enddo
  
     enddo
  
     ! 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_s, eb_s, fluxes, hml, dt, g, gv, us, tv, &
                               cs%optics, cs%tracer_flow_CSp, cs%debug, &
                               evap_cfl_limit = cs%evap_CFL_limit, &
                               minimum_forcing_depth=cs%minimum_forcing_depth)
  
   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_s(i,j,nz) ; eatr(i,j,1) = ea_s(i,j,1)
     enddo ; enddo
     !$OMP parallel do default(shared) 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%Z_to_H**2) / &
            (0.5 * (h(i,j,k-1) + h(i,j,k))  + &
             h_neglect)
       else
         add_ent = 0.0
       endif
       ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
       eatr(i,j,k) = ea_s(i,j,k) + add_ent
     enddo ; enddo ; enddo
  
     ! 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, us, tv, &
                                 cs%optics, cs%tracer_flow_CSp, cs%debug,&
                                 evap_cfl_limit = cs%evap_CFL_limit, &
                                 minimum_forcing_depth=cs%minimum_forcing_depth)
   else
     ! 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, us, tv, &
                                 cs%optics, cs%tracer_flow_CSp, cs%debug, &
                                 evap_cfl_limit = cs%evap_CFL_limit, &
                                 minimum_forcing_depth=cs%minimum_forcing_depth)
   endif  ! (CS%mix_boundary_tracers)
  
   call cpu_clock_end(id_clock_tracers)
  
   ! Apply ALE sponge
   if (cs%use_sponge) then
     call cpu_clock_begin(id_clock_sponge)
     if (associated(cs%ALE_sponge_CSp)) then
       call apply_ale_sponge(h, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)
     endif
  
     call cpu_clock_end(id_clock_sponge)
     if (cs%debug) then
       call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)
       call mom_thermovar_chksum("apply_sponge ", tv, g)
     endif
   endif ! CS%use_sponge
  
   call cpu_clock_begin(id_clock_pass)
   if (g%symmetric) then ; dir_flag = to_all+omit_corners
   else ; dir_flag = to_west+to_south+omit_corners ; endif
   call create_group_pass(cs%pass_hold_eb_ea, eb_t, g%Domain, dir_flag, halo=1)
   call create_group_pass(cs%pass_hold_eb_ea, eb_s, g%Domain, dir_flag, halo=1)
   call create_group_pass(cs%pass_hold_eb_ea, ea_t, g%Domain, dir_flag, halo=1)
   call create_group_pass(cs%pass_hold_eb_ea, ea_s, g%Domain, dir_flag, halo=1)
   call do_group_pass(cs%pass_hold_eb_ea, g%Domain)
   ! visc%Kv_slow is not in the group pass because it has larger vertical extent.
   if (associated(visc%Kv_slow)) &
     call pass_var(visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)
   call cpu_clock_end(id_clock_pass)
  
   call disable_averaging(cs%diag)
  
   ! Diagnose the diapycnal diffusivities and other related quantities.
   call enable_averages(dt, time_end, cs%diag)
  
   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_t       > 0) call post_data(cs%id_ea_t,       ea_t, cs%diag)
   if (cs%id_eb_t       > 0) call post_data(cs%id_eb_t,       eb_t, cs%diag)
   if (cs%id_ea_s       > 0) call post_data(cs%id_ea_s,       ea_s, cs%diag)
   if (cs%id_eb_s       > 0) call post_data(cs%id_eb_s,       eb_s, 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_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)
  
   call disable_averaging(cs%diag)
  
   if (showcalltree) call calltree_leave("diabatic_ALE()")
  

References mom_ale_sponge::apply_ale_sponge(), mom_cvmix_conv::calculate_cvmix_conv(), 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(), mom_domains::do_group_pass(), mom_energetic_pbl::energetic_pbl_get_mld(), mom_geothermal::geothermal(), id_clock_differential_diff, id_clock_geothermal, id_clock_kpp, id_clock_pass, id_clock_remap, id_clock_set_diffusivity, id_clock_sponge, id_clock_tracers, id_clock_tridiag, mom_cvmix_kpp::kpp_nonlocaltransport_saln(), mom_cvmix_kpp::kpp_nonlocaltransport_temp(), mom_checksum_packages::mom_state_stats(), mom_variables::mom_thermovar_chksum(), mom_opacity::optics_nbands(), mom_diabatic_aux::set_pen_shortwave(), and mom_tracer_diabatic::tracer_vertdiff().

Referenced by diabatic().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ diabatic_ale_legacy()

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

private

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

Parameters

[in,out]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	u	zonal velocity [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
	cs	module control structure
	waves	Surface gravity waves

Definition at line 446 of file MOM_diabatic_driver.F90.

   type(ocean_grid_type),                     intent(inout) :: G         !< ocean grid structure
   type(verticalGrid_type),                   intent(in)    :: GV        !< ocean vertical grid structure
   type(unit_scale_type),                     intent(in)    :: US        !< A dimensional unit scaling type
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u         !< zonal velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v         !< meridional velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h         !< thickness [H ~> m or kg m-2]
   type(thermo_var_ptrs),                     intent(inout) :: tv        !< points to thermodynamic fields
                                                                         !! unused have NULL ptrs
   real, dimension(:,:),                      pointer       :: Hml       !< Active mixed layer depth [m]
   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
   type(accel_diag_ptrs),                     intent(inout) :: ADp       !< related points to accelerations in momentum
                                                                         !! equations, to enable the later derived
                                                                         !! diagnostics, like energy budgets
   type(cont_diag_ptrs),                      intent(inout) :: CDp       !< points to terms in continuity equations
   real,                                      intent(in)    :: dt        !< time increment [T ~> s]
   type(time_type),                           intent(in)    :: Time_end  !< Time at the end of the interval
   type(diabatic_CS),                         pointer       :: CS        !< module control structure
   type(Wave_parameters_CS),        optional, pointer       :: Waves     !< Surface gravity waves
  
   ! local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
     ea_s,     &  ! amount of fluid entrained from the layer above within
                  ! one time step  [H ~> m or kg m-2]
     eb_s,     &  ! amount of fluid entrained from the layer below within
                  ! one time step  [H ~> m or kg m-2]
     ea_t,     &  ! amount of fluid entrained from the layer above within
                  ! one time step  [H ~> m or kg m-2]
     eb_t,     &  ! amount of fluid entrained from the layer below within
                  ! one time step  [H ~> m or kg m-2]
     kd_lay, &    ! diapycnal diffusivity of layers [Z2 T-1 ~> m2 s-1]
     h_orig, &    ! initial layer thicknesses [H ~> m or kg m-2]
     h_prebound, & ! initial layer thicknesses [H ~> m or kg m-2]
     hold,   &    ! layer thickness before diapycnal entrainment, and later
                  ! the initial layer thicknesses (if a mixed layer is used),
                  ! [H ~> m or kg m-2]
     dsv_dt, &    ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
     dsv_ds, &    ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
     ctke,   &    ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
     u_h,    &    ! zonal and meridional velocities at thickness points after
     v_h          ! entrainment [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJ_(G)) :: &
     Rcv_ml, &   ! coordinate density of mixed layer, used for applying sponges
     SkinBuoyFlux! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], 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 :: net_ent  ! The net of ea-eb at an interface.
  
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
     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 [H ~> m or kg m-2] (for Bous or non-Bouss)
  
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), target :: &
     Kd_int,   & ! diapycnal diffusivity of interfaces [Z2 T-1 ~> m2 s-1]
     Kd_heat,  & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1]
     Kd_salt,  & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1]
     Kd_ePBL,  & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]
     Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
     Tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
     Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
     Sadv_flx    ! advective diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
  
   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
                        ! [H ~> m or kg m-2]
   real :: h_neglect    ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected
                        ! [H ~> m or kg m-2]
   real :: h_neglect2   ! h_neglect^2 [H2 ~> m2 or kg2 m-4]
   real :: add_ent      ! Entrainment that needs to be added when mixing tracers
                        ! [H ~> m or kg m-2]
   real :: eaval        ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
   real :: hval         ! hval is 2*h at velocity grid points [H ~> m or kg m-2]
   real :: h_tr         ! h_tr is h at tracer points with a tiny thickness
                        ! added to ensure positive definiteness [H ~> m or kg m-2]
   real :: Tr_ea_BBL    ! The diffusive tracer thickness in the BBL that is
                        ! coupled to the bottom within a timestep [H ~> m or kg m-2]
  
   real :: htot(SZIB_(G))             ! The summed thickness from the bottom [H ~> m or kg m-2].
   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 ~> m or kg m-2]
   real :: Idt     ! The inverse time step [T-1 ~> s-1]
  
   integer :: dir_flag     ! An integer encoding the directions in which to do halo updates.
   logical :: showCallTree ! If true, show the call tree
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo
  
   integer :: ig, jg      ! global indices for testing testing itide point source (BDM)
   real :: Kd_add_here    ! An added diffusivity [Z2 T-1 ~> m2 s-1].
  
   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
  
   showcalltree = calltree_showquery()
   if (showcalltree) call calltree_enter("diabatic_ALE_legacy(), MOM_diabatic_driver.F90")
 !   if (showCallTree) call callTree_enter("diabatic_ALE(), MOM_diabatic_driver.F90")
  
   ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
   call enable_averages(dt, time_end, cs%diag)
  
   if (cs%use_geothermal) then
     halo = cs%halo_TS_diff
     !$OMP parallel do default(shared)
     do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo
       h_orig(i,j,k) = h(i,j,k) ; eatr(i,j,k) = 0.0 ; ebtr(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, eatr, ebtr, g, gv, us, cs%geothermal_CSp, halo=cs%halo_TS_diff)
     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, gv, us)
   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_pen_shortwave 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_pen_shortwave(cs%optics, fluxes, g, gv, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
  
   if (cs%debug) call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)
  
   if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then
     if (cs%use_geothermal) then
       call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, us, eatr, ebtr)
       if (cs%debug) then
         call hchksum(eatr, "after find_uv_at_h eatr",g%HI, scale=gv%H_to_m)
         call hchksum(ebtr, "after find_uv_at_h ebtr",g%HI, scale=gv%H_to_m)
       endif
     else
       call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
     endif
     if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
   endif
  
   call cpu_clock_begin(id_clock_set_diffusivity)
   ! Sets: Kd_lay, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S and visc%TKE_turb
   ! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow
   call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, &
                        visc, dt, g, gv, us, cs%set_diff_CSp, kd_lay, kd_int)
   call cpu_clock_end(id_clock_set_diffusivity)
   if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")
  
   ! Set diffusivities for heat and salt separately
  
   if (.not.cs%use_legacy_diabatic .or. cs%useKPP) then
     !$OMP parallel do default(shared)
     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
     ! Add contribution from double diffusion
     if (associated(visc%Kd_extra_S)) then
       !$OMP parallel do default(shared)
       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 parallel do default(shared)
       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
   endif
  
   if (cs%debug) then
     call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)
     call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)
     call mom_thermovar_chksum("after set_diffusivity ", tv, g)
     call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     call hchksum(kd_heat, "after set_diffusivity Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     call hchksum(kd_salt, "after set_diffusivity Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
   endif
  
   if (cs%useKPP) then
     call cpu_clock_begin(id_clock_kpp)
     ! total vertical viscosity in the interior is represented via visc%Kv_shear
     if (.not.cs%use_legacy_diabatic) then
       do k=1,nz+1 ; do j=js,je ; do i=is,ie
         visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k)
       enddo ; enddo ; enddo
     endif
  
     ! 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, us, 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.
  
     call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv%eqn_of_state, &
                          fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
  
     call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
                        kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
  
     if (associated(hml)) then
       !$OMP parallel default(shared)
       call kpp_get_bld(cs%KPP_CSp, hml(:,:), g)
       !$OMP end parallel
       call pass_var(hml, g%domain, halo=1)
       ! If visc%MLD exists, copy KPP's BLD into it
       if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
     endif
  
     if (cs%use_legacy_diabatic .and. .not.cs%KPPisPassive) then
       !$OMP parallel do default(shared)
       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 parallel do default(shared)
         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 parallel do default(shared)
         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
  
     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, us, haloshift=0)
       call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)
       call mom_thermovar_chksum("after KPP", tv, g)
       call hchksum(kd_heat, "after KPP Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
       call hchksum(kd_salt, "after KPP Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     endif
  
   endif  ! endif for KPP
  
  
   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, &
                                     us%T_to_s*dt, tv%T, tv%C_p)
     call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, &
                                     us%T_to_s*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, gv, us)
  
     if (cs%debug) then
       call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)
       call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)
       call mom_thermovar_chksum("after KPP_applyNLT ", tv, g)
     endif
   endif ! endif for KPP
  
   ! This is the "old" method for applying differential diffusion.
   ! Changes: tv%T, tv%S
   if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T) .and. &
       (cs%use_legacy_diabatic .or. .not.cs%use_CVMix_ddiff)) 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, gv, us)
  
     ! increment heat and salt diffusivity.
     ! CS%useKPP==.true. already has extra_T and extra_S included
     if (.not. cs%useKPP) then
       !$OMP parallel do default(shared)
       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
  
   ! Calculate vertical mixing due to convection (computed via CVMix)
   if (cs%use_CVMix_conv) then
     call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml)
     ! Increment vertical diffusion and viscosity due to convection
     if (cs%use_legacy_diabatic) then
       !$OMP parallel do default(shared)
       do k=1,nz+1 ; do j=js,je ; do i=is,ie
         kd_int(i,j,k) = kd_int(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
         visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
       enddo ; enddo ; enddo
     else
       !$OMP parallel do default(shared)
       do k=1,nz+1 ; do j=js,je ; do i=is,ie
         kd_heat(i,j,k) = kd_heat(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
         kd_salt(i,j,k) = kd_salt(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
         if (cs%useKPP) then
           visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
         else
           visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
         endif
       enddo ; enddo ; enddo
     endif
   endif
  
   ! This block sets ea, eb from h and Kd_int.
   if (cs%use_legacy_diabatic) then
     do j=js,je ; do i=is,ie
       ea_s(i,j,1) = 0.0
     enddo ; enddo
     !$OMP parallel do default(shared)  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_s(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_int(i,j,k)
       eb_s(i,j,k-1) = ea_s(i,j,k)
       ea_t(i,j,k-1) = ea_s(i,j,k-1) ; eb_t(i,j,k-1) = eb_s(i,j,k-1)
     enddo ; enddo ; enddo
     do j=js,je ; do i=is,ie
       eb_s(i,j,nz) = 0.0
       ea_t(i,j,nz) = ea_s(i,j,nz) ; eb_t(i,j,nz) = eb_s(i,j,nz)
     enddo ; enddo
     if (showcalltree) call calltree_waypoint("done setting ea,eb from Kd_int (diabatic)")
   endif
  
   if (cs%debug) then
     call mom_forcing_chksum("after calc_entrain ", fluxes, g, us, haloshift=0)
     call mom_thermovar_chksum("after calc_entrain ", tv, g)
     call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, us, haloshift=0)
     call hchksum(ea_s, "after calc_entrain ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
     call hchksum(eb_s, "after calc_entrain eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
   endif
  
   ! Save fields before boundary forcing is applied for tendency 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
  
   ! Apply forcing
   call cpu_clock_begin(id_clock_remap)
  
   ! Changes made to following fields:  h, tv%T and tv%S.
   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, us, dt, fluxes, cs%optics, &
             optics_nbands(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_t, "after applyBoundaryFluxes ea_t",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_t, "after applyBoundaryFluxes eb_t",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(ea_s, "after applyBoundaryFluxes ea_s",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_s, "after applyBoundaryFluxes eb_s",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(ctke, "after applyBoundaryFluxes cTKE", g%HI, haloshift=0, &
                    scale=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**2)
       call hchksum(dsv_dt, "after applyBoundaryFluxes dSV_dT", g%HI, haloshift=0, scale=us%kg_m3_to_R)
       call hchksum(dsv_ds, "after applyBoundaryFluxes dSV_dS", g%HI, haloshift=0, scale=us%kg_m3_to_R)
     endif
  
     call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
     call energetic_pbl(h, u_h, v_h, tv, fluxes, dt, kd_epbl, g, gv, us, &
                        cs%energetic_PBL_CSp, dsv_dt, dsv_ds, ctke, skinbuoyflux, waves=waves)
  
     if (associated(hml)) then
       call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hml(:,:), g, us)
       call pass_var(hml, g%domain, halo=1)
       ! If visc%MLD exists, copy ePBL's MLD into it
       if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
     elseif (associated(visc%MLD)) then
       call energetic_pbl_get_mld(cs%energetic_PBL_CSp, visc%MLD, g, us)
       call pass_var(visc%MLD, g%domain, halo=1)
     endif
  
     ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
     do k=2,nz ; do j=js,je ; do i=is,ie
       !### These expressions assume a Prandtl number of 1.
       if (cs%ePBL_is_additive) then
         kd_add_here = kd_epbl(i,j,k)
         visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + kd_epbl(i,j,k)
       else
         kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)
         visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), kd_epbl(i,j,k))
       endif
  
       if (cs%use_legacy_diabatic) then
         ent_int = kd_add_here * (gv%Z_to_H**2 * dt) / &
                     (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect)
         eb_s(i,j,k-1) = eb_s(i,j,k-1) + ent_int
         ea_s(i,j,k) = ea_s(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)
       else
         kd_heat(i,j,k) = kd_heat(i,j,k) + kd_add_here
         kd_salt(i,j,k) = kd_salt(i,j,k) + kd_add_here
       endif
  
     enddo ; enddo ; enddo
  
     if (cs%debug) then
       call hchksum(ea_t, "after ePBL ea_t",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_t, "after ePBL eb_t",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(ea_s, "after ePBL ea_s",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_s, "after ePBL eb_s",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(kd_epbl, "after ePBL Kd_ePBL", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     endif
  
   else
     call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
                                   optics_nbands(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
   ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme
   !  so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards
   if (cs%boundary_forcing_tendency_diag) then
     call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, g, gv, cs)
     if (cs%id_boundary_forcing_h > 0) call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h = h_diag)
   endif
   ! Boundary fluxes may have changed T, S, and h
   call diag_update_remap_grids(cs%diag)
   call cpu_clock_end(id_clock_remap)
   if (cs%debug) then
     call mom_forcing_chksum("after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)
     call mom_thermovar_chksum("after applyBoundaryFluxes ", tv, g)
     call mom_state_chksum("after applyBoundaryFluxes ", u, v, h, g, gv, us, 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, gv, us)
  
   ! 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 code may be unnecessary, but the negative-thickness checks do appear to change
   !    answers slightly in some cases.
   if (cs%use_legacy_diabatic) then
     !$OMP parallel do default(shared)
     do j=js,je
       do i=is,ie
         hold(i,j,1) = h(i,j,1)
    ! Does nothing with ALE:  h(i,j,1) = h(i,j,1) + (eb_s(i,j,1) - ea_s(i,j,2))
         hold(i,j,nz) = h(i,j,nz)
    ! Does nothing with ALE:  h(i,j,nz) = h(i,j,nz) + (ea_s(i,j,nz) - eb_s(i,j,nz-1))
         if (h(i,j,1) <= 0.0) h(i,j,1) = gv%Angstrom_H
         if (h(i,j,nz) <= 0.0) h(i,j,nz) = gv%Angstrom_H
       enddo
       do k=2,nz-1 ; do i=is,ie
         hold(i,j,k) = h(i,j,k)
    ! Does nothing with ALE:  h(i,j,k) = h(i,j,k) + ((ea_s(i,j,k) - eb_s(i,j,k-1)) + &
    !                                                (eb_s(i,j,k) - ea_s(i,j,k+1)))
         if (h(i,j,k) <= 0.0) h(i,j,k) = gv%Angstrom_H
       enddo ; enddo
     enddo
     ! Checks for negative thickness may have changed layer thicknesses
     call diag_update_remap_grids(cs%diag)
   endif
  
   if (cs%debug) then
     call mom_state_chksum("after negative check ", u, v, h, g, gv, us, haloshift=0)
     call mom_forcing_chksum("after negative check ", fluxes, g, us, 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, gv, us)
  
   ! calculate change in temperature & salinity due to dia-coordinate surface diffusion
   if (associated(tv%T)) then
  
     if (cs%debug) then
       call hchksum(ea_t, "before triDiagTS ea_t ",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_t, "before triDiagTS eb_t ",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(ea_s, "before triDiagTS ea_s ",g%HI,haloshift=0, scale=gv%H_to_m)
       call hchksum(eb_s, "before triDiagTS eb_s ",g%HI,haloshift=0, scale=gv%H_to_m)
     endif
  
     call cpu_clock_begin(id_clock_tridiag)
     !  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)
  
     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
  
     if (cs%use_legacy_diabatic) then
       ! Changes T and S via the tridiagonal solver; no change to h
       do k=1,nz ; do j=js,je ; do i=is,ie
         ea_t(i,j,k) = ea_s(i,j,k) ; eb_t(i,j,k) = eb_s(i,j,k)
       enddo ; enddo ; enddo
       if (cs%tracer_tridiag) then
         call tracer_vertdiff(hold, ea_t, eb_t, dt, tv%T, g, gv)
         call tracer_vertdiff(hold, ea_s, eb_s, dt, tv%S, g, gv)
       else
         call tridiagts(g, gv, is, ie, js, je, hold, ea_s, eb_s, tv%T, tv%S)
       endif
  
       ! diagnose temperature, salinity, heat, and salt tendencies
       ! Note: hold here refers to the thicknesses from before the dual-entraintment when using
       ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will (not?) have changed
       ! In either case, tendencies should be posted on hold
       if (cs%diabatic_diff_tendency_diag) then
         call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, cs)
         if (cs%id_diabatic_diff_h > 0) call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h = hold)
       endif
     else
       ! Set ea_t=eb_t based on Kd_heat and ea_s=eb_s based on Kd_salt on interfaces for use in the tri-diagonal solver.
  
       do j=js,je ; do i=is,ie
         ea_t(i,j,1) = 0.; ea_s(i,j,1) = 0.
       enddo ; enddo
  
       !$OMP parallel do default(shared) 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_t(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_heat(i,j,k)
         eb_t(i,j,k-1) = ea_t(i,j,k)
         ea_s(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_salt(i,j,k)
         eb_s(i,j,k-1) = ea_s(i,j,k)
       enddo ; enddo ; enddo
       do j=js,je ; do i=is,ie
         eb_t(i,j,nz) = 0. ; eb_s(i,j,nz) = 0.
       enddo ; enddo
       if (showcalltree) call calltree_waypoint("done setting ea_t,ea_s,eb_t,eb_s from Kd_heat" //&
          "and Kd_salt (diabatic)")
  
     ! Changes T and S via the tridiagonal solver; no change to h
       call tracer_vertdiff(h, ea_t, eb_t, dt, tv%T, g, gv)
       call tracer_vertdiff(h, ea_s, eb_s, dt, tv%S, g, gv)
  
       ! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below
       if (cs%diabatic_diff_tendency_diag) &
         call diagnose_diabatic_diff_tendency(tv, h, 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, gv, us)
  
   if (cs%debug) then
     call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)
     call mom_thermovar_chksum("after mixed layer ", tv, g)
   endif
  
   ! Whenever thickness changes let the diag manager know, as the
   ! target grids for vertical remapping may need to be regenerated.
   if (cs%id_dudt_dia > 0 .or. cs%id_dvdt_dia > 0) &
     ! Remapped d[uv]dt_dia require east/north halo updates of h
     call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
   call diag_update_remap_grids(cs%diag)
  
   ! diagnostics
   idt = 1.0 / dt
   if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 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(shared)
     do k=2,nz ; do j=js,je ; do i=is,ie
       tdif_flx(i,j,k) = (idt * 0.5*(ea_t(i,j,k) + eb_t(i,j,k-1))) * &
                         (tv%T(i,j,k-1) - tv%T(i,j,k))
       tadv_flx(i,j,k) = (idt * (ea_t(i,j,k) - eb_t(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_Sadv > 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(shared)
     do k=2,nz ; do j=js,je ; do i=is,ie
       sdif_flx(i,j,k) = (idt * 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))) * &
                         (tv%S(i,j,k-1) - tv%S(i,j,k))
       sadv_flx(i,j,k) = (idt * (ea_s(i,j,k) - eb_s(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 = gv%Z_to_H * sqrt(dt*cs%Kd_BBL_tr)
     !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
     do j=js,je
       do i=is,ie
         ebtr(i,j,nz) = eb_s(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_int,
           ! perhaps a molecular diffusivity.
           add_ent = ((dt * cs%Kd_min_tr) * gv%Z_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_s(i,j,k) + eb_s(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_s(i,j,k),eb_s(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_s(i,j,k-1) + add_ent
           eatr(i,j,k) = ea_s(i,j,k) + add_ent
         else
           ebtr(i,j,k-1) = eb_s(i,j,k-1) ; eatr(i,j,k) = ea_s(i,j,k)
         endif
  
         if (associated(visc%Kd_extra_S)) then ; if (visc%Kd_extra_S(i,j,k) > 0.0) then
           if (cs%use_legacy_diabatic) then
             add_ent = ((dt * visc%Kd_extra_S(i,j,k)) * gv%Z_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 = ((dt * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
                (0.5 * (h(i,j,k-1) + h(i,j,k)) + &
                 h_neglect)
           endif
           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_s(i,j,1) ; enddo
  
     enddo
  
     ! 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_s, eb_s, fluxes, hml, dt, g, gv, us, tv, &
                               cs%optics, cs%tracer_flow_CSp, cs%debug, &
                               evap_cfl_limit = cs%evap_CFL_limit, &
                               minimum_forcing_depth=cs%minimum_forcing_depth)
  
   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_s(i,j,nz) ; eatr(i,j,1) = ea_s(i,j,1)
     enddo ; enddo
     !$OMP parallel do default(shared) 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
         if (cs%use_legacy_diabatic) then
           add_ent = ((dt * visc%Kd_extra_S(i,j,k)) * gv%Z_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 = ((dt * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
              (0.5 * (h(i,j,k-1) + h(i,j,k)) + &
               h_neglect)
         endif
       else
         add_ent = 0.0
       endif
       ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
       eatr(i,j,k) = ea_s(i,j,k) + add_ent
     enddo ; enddo ; enddo
  
     ! 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, us, tv, &
                                 cs%optics, cs%tracer_flow_CSp, cs%debug,&
                                 evap_cfl_limit = cs%evap_CFL_limit, &
                                 minimum_forcing_depth=cs%minimum_forcing_depth)
   else
     ! 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, us, tv, &
                                 cs%optics, cs%tracer_flow_CSp, cs%debug, &
                                 evap_cfl_limit = cs%evap_CFL_limit, &
                                 minimum_forcing_depth=cs%minimum_forcing_depth)
   endif  ! (CS%mix_boundary_tracers)
  
   call cpu_clock_end(id_clock_tracers)
  
   ! Apply ALE sponge
   if (cs%use_sponge) then
     call cpu_clock_begin(id_clock_sponge)
     if (associated(cs%ALE_sponge_CSp)) then
       call apply_ale_sponge(h, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)
     endif
  
     call cpu_clock_end(id_clock_sponge)
     if (cs%debug) then
       call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)
       call mom_thermovar_chksum("apply_sponge ", tv, g)
     endif
   endif ! CS%use_sponge
  
   call disable_averaging(cs%diag)
   ! Diagnose the diapycnal diffusivities and other related quantities.
   call enable_averages(dt, time_end, cs%diag)
  
   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_s, cs%diag)
   if (cs%id_eb         > 0) call post_data(cs%id_eb,         eb_s, cs%diag)
   if (cs%id_ea_t       > 0) call post_data(cs%id_ea_t,       ea_t, cs%diag)
   if (cs%id_eb_t       > 0) call post_data(cs%id_eb_t,       eb_t, cs%diag)
   if (cs%id_ea_s       > 0) call post_data(cs%id_ea_s,       ea_s, cs%diag)
   if (cs%id_eb_s       > 0) call post_data(cs%id_eb_s,       eb_s, 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_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)
  
   call disable_averaging(cs%diag)
  
   if (showcalltree) call calltree_leave("diabatic_ALE_legacy()")
  

References mom_ale_sponge::apply_ale_sponge(), mom_cvmix_conv::calculate_cvmix_conv(), 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(), mom_energetic_pbl::energetic_pbl_get_mld(), mom_geothermal::geothermal(), id_clock_differential_diff, id_clock_geothermal, id_clock_kpp, id_clock_remap, id_clock_set_diffusivity, id_clock_sponge, id_clock_tracers, id_clock_tridiag, mom_cvmix_kpp::kpp_nonlocaltransport_saln(), mom_cvmix_kpp::kpp_nonlocaltransport_temp(), mom_checksum_packages::mom_state_stats(), mom_variables::mom_thermovar_chksum(), mom_opacity::optics_nbands(), mom_diabatic_aux::set_pen_shortwave(), and mom_tracer_diabatic::tracer_vertdiff().

Referenced by diabatic().

Here is the call graph for this function:

Here is the caller 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 3758 of file MOM_diabatic_driver.F90.

   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%use_tidal_mixing) call tidal_mixing_end(cs%tidal_mixing_CSp)
  
   if (cs%use_CVMix_conv) call cvmix_conv_end(cs%CVMix_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
  
   ! GMM, the following is commented out because arrays in
   ! CS%diag_grids_prev are neither pointers or allocatables
   ! and, therefore, cannot be deallocated.
  
   !call diag_grid_storage_end(CS%diag_grids_prev)
  
   deallocate(cs)
  

References mom_cvmix_conv::cvmix_conv_end(), mom_entrain_diffusive::entrain_diffusive_end(), and mom_tidal_mixing::tidal_mixing_end().

Referenced by mom::mom_end().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ diabatic_driver_init()

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

This routine initializes the diabatic driver module.

Parameters

	time	model time
[in,out]	g	model grid structure
[in]	gv	model vertical grid structure
[in]	us	A dimensional unit scaling type
[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 the tracer flow control module
	sponge_csp	pointer to the sponge module control structure
	ale_sponge_csp	pointer to the ALE sponge module control structure

Definition at line 3187 of file MOM_diabatic_driver.F90.

   type(time_type), target                :: Time             !< model time
   type(ocean_grid_type),   intent(inout) :: G                !< model grid structure
   type(verticalGrid_type), intent(in)    :: GV               !< model vertical grid structure
   type(unit_scale_type),   intent(in)    :: US               !< A dimensional unit scaling type
   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 the
                                                              !! 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
  
   real    :: Kd
   integer :: num_mode
   logical :: use_temperature, differentialDiffusion
  
 ! This "include" declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mdl = "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
   cs%Time => time
  
   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
  
   cs%useALEalgorithm = usealealgorithm
   cs%bulkmixedlayer = (gv%nkml > 0)
  
   ! Set default, read and log parameters
   call log_version(param_file, mdl, version, &
                    "The following parameters are used for diabatic processes.")
   call get_param(param_file, mdl, "USE_LEGACY_DIABATIC_DRIVER", cs%use_legacy_diabatic, &
                  "If true, use a legacy version of the diabatic subroutine. "//&
                  "This is temporary and is needed to avoid change in answers.", &
                  default=.true.)
   call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &
                  "If true, sponges may be applied anywhere in the domain. "//&
                  "The exact location and properties of those sponges are "//&
                  "specified via calls to initialize_sponge and possibly "//&
                  "set_up_sponge_field.", default=.false.)
   call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
                  "If true, temperature and salinity are used as state "//&
                  "variables.", default=.true.)
   call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", cs%use_energetic_PBL, &
                  "If true, use an implied energetics planetary boundary "//&
                  "layer scheme to determine the diffusivity and viscosity "//&
                  "in the surface boundary layer.", default=.false.)
   call get_param(param_file, mdl, "EPBL_IS_ADDITIVE", cs%ePBL_is_additive, &
                  "If true, the diffusivity from ePBL is added to all "//&
                  "other diffusivities. Otherwise, the larger of kappa-shear "//&
                  "and ePBL diffusivities are used.", default=.true.)
   call get_param(param_file, mdl, "DOUBLE_DIFFUSION", differentialdiffusion, &
                  "If true, apply parameterization of double-diffusion.", &
                  default=.false. )
   call get_param(param_file, mdl, "USE_KPP", cs%use_KPP, &
                  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
                  "to calculate diffusivities and non-local transport in the OBL.", &
                  default=.false., do_not_log=.true.)
   cs%use_CVMix_ddiff = cvmix_ddiff_is_used(param_file)
  
   if (cs%use_CVMix_ddiff .and. differentialdiffusion) then
     call mom_error(fatal, 'diabatic_driver_init: '// &
            'Multiple double-diffusion options selected (DOUBLE_DIFFUSION and'//&
            'USE_CVMIX_DDIFF), please disable all but one option to proceed.')
   endif
  
   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, mdl, "ML_MIX_FIRST", cs%ML_mix_first, &
                  "The fraction of the mixed layer mixing that is applied "//&
                  "before interior diapycnal mixing.  0 by default.", &
                  units="nondim", default=0.0)
     call get_param(param_file, mdl, "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, mdl, "DO_GEOTHERMAL", cs%use_geothermal, &
                  "If true, apply geothermal heating.", default=.false.)
   else
     cs%use_geothermal = .false.
   endif
   call get_param(param_file, mdl, "INTERNAL_TIDES", cs%use_int_tides, &
                  "If true, use the code that advances a separate set of "//&
                  "equations for the internal tide energy density.", default=.false.)
   cs%nMode = 1
   if (cs%use_int_tides) then
     call get_param(param_file, mdl, "INTERNAL_TIDE_MODES", cs%nMode, &
                  "The number of distinct internal tide modes "//&
                  "that will be calculated.", default=1, do_not_log=.true.)
     call get_param(param_file, mdl, "UNIFORM_TEST_CG", cs%uniform_test_cg, &
                  "If positive, a uniform group velocity of internal tide for test case", &
                  default=-1., units="m s-1", scale=us%m_s_to_L_T)
   endif
  
   call get_param(param_file, mdl, "MASSLESS_MATCH_TARGETS", &
                                 cs%massless_match_targets, &
                  "If true, the temperature and salinity of massless layers "//&
                  "are kept consistent with their target densities. "//&
                  "Otherwise the properties of massless layers evolve "//&
                  "diffusively to match massive neighboring layers.", &
                  default=.true.)
  
   call get_param(param_file, mdl, "AGGREGATE_FW_FORCING", cs%aggregate_FW_forcing, &
                  "If true, the net incoming and outgoing fresh water fluxes are combined "//&
                  "and applied as either incoming or outgoing depending on the sign of the net. "//&
                  "If false, the net incoming fresh water flux is added to the model and "//&
                  "thereafter the net outgoing is removed from the topmost non-vanished "//&
                  "layers of the updated state.", default=.true.)
  
   call get_param(param_file, mdl, "DEBUG", cs%debug, &
                  "If true, write out verbose debugging data.", &
                  default=.false., debuggingparam=.true.)
   call get_param(param_file, mdl, "DEBUG_CONSERVATION", cs%debugConservation, &
                  "If true, monitor conservation and extrema.", &
                  default=.false., debuggingparam=.true.)
  
   call get_param(param_file, mdl, "DEBUG_ENERGY_REQ", cs%debug_energy_req, &
                  "If true, debug the energy requirements.", default=.false., do_not_log=.true.)
   call get_param(param_file, mdl, "MIX_BOUNDARY_TRACERS", cs%mix_boundary_tracers, &
                  "If true, mix the passive tracers in massless layers at "//&
                  "the bottom into the interior as though a diffusivity of "//&
                  "KD_MIN_TR were operating.", default=.true.)
  
   if (cs%mix_boundary_tracers) then
     call get_param(param_file, mdl, "KD", kd, fail_if_missing=.true.)
     call get_param(param_file, mdl, "KD_MIN_TR", cs%Kd_min_tr, &
                  "A minimal diffusivity that should always be applied to "//&
                  "tracers, especially in massless layers near the bottom. "//&
                  "The default is 0.1*KD.", units="m2 s-1", default=0.1*kd, scale=us%m2_s_to_Z2_T)
     call get_param(param_file, mdl, "KD_BBL_TR", cs%Kd_BBL_tr, &
                  "A bottom boundary layer tracer diffusivity that will "//&
                  "allow for explicitly specified bottom fluxes. The "//&
                  "entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) "//&
                  "over the same distance.", units="m2 s-1", default=0., scale=us%m2_s_to_Z2_T)
   endif
  
   call get_param(param_file, mdl, "TRACER_TRIDIAG", cs%tracer_tridiag, &
                  "If true, use the passive tracer tridiagonal solver for T and S", &
                  default=.false.)
  
   call get_param(param_file, mdl, "MINIMUM_FORCING_DEPTH", cs%minimum_forcing_depth, &
                  "The smallest depth over which forcing can be applied. This "//&
                  "only takes effect when near-surface layers become thin "//&
                  "relative to this scale, in which case the forcing tendencies "//&
                  "scaled down by distributing the forcing over this depth scale.", &
                  units="m", default=0.001, scale=gv%m_to_H)
   call get_param(param_file, mdl, "EVAP_CFL_LIMIT", cs%evap_CFL_limit, &
                  "The largest fraction of a layer than can be lost to forcing "//&
                  "(e.g. evaporation, sea-ice formation) in one time-step. The unused "//&
                  "mass loss is passed down through the column.", &
                  units="nondim", default=0.8)
  
  
   ! Register all available diagnostics for this module.
   thickness_units = get_thickness_units(gv)
  
   cs%id_ea_t = register_diag_field('ocean_model','ea_t',diag%axesTL,time, &
       'Layer (heat) entrainment from above per timestep','m', &
       conversion=gv%H_to_m)
   cs%id_eb_t = register_diag_field('ocean_model','eb_t',diag%axesTL,time, &
       'Layer (heat) entrainment from below per timestep', 'm', &
       conversion=gv%H_to_m)
   cs%id_ea_s = register_diag_field('ocean_model','ea_s',diag%axesTL,time, &
       'Layer (salt) entrainment from above per timestep','m', &
       conversion=gv%H_to_m)
   cs%id_eb_s = register_diag_field('ocean_model','eb_s',diag%axesTL,time, &
       'Layer (salt) entrainment from below per timestep', 'm', &
       conversion=gv%H_to_m)
   ! used by layer diabatic
   cs%id_ea = register_diag_field('ocean_model','ea',diag%axesTL,time, &
       'Layer entrainment from above per timestep','m', &
       conversion=gv%H_to_m)
   cs%id_eb = register_diag_field('ocean_model','eb',diag%axesTL,time, &
       'Layer entrainment from below per timestep', 'm', &
       conversion=gv%H_to_m)
   cs%id_wd = register_diag_field('ocean_model','wd',diag%axesTi,time, &
     'Diapycnal velocity', 'm s-1', conversion=gv%H_to_m)
   if (cs%id_wd > 0) call safe_alloc_ptr(cdp%diapyc_vel,isd,ied,jsd,jed,nz+1)
  
   cs%id_dudt_dia = register_diag_field('ocean_model','dudt_dia',diag%axesCuL,time, &
       'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)
   cs%id_dvdt_dia = register_diag_field('ocean_model','dvdt_dia',diag%axesCvL,time, &
       'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)
  
   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
  
   if (use_temperature) then
     cs%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff",diag%axesTi, &
         time, "Diffusive diapycnal temperature flux across interfaces", &
         "degC m s-1", conversion=gv%H_to_m*us%s_to_T)
     cs%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv",diag%axesTi, &
         time, "Advective diapycnal temperature flux across interfaces", &
         "degC m s-1", conversion=gv%H_to_m*us%s_to_T)
     cs%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff",diag%axesTi, &
         time, "Diffusive diapycnal salnity flux across interfaces", &
         "psu m s-1", conversion=gv%H_to_m*us%s_to_T)
     cs%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv",diag%axesTi, &
         time, "Advective diapycnal salnity flux across interfaces", &
         "psu m s-1", conversion=gv%H_to_m*us%s_to_T)
     cs%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, time, &
         'Mixed layer depth (delta rho = 0.03)', 'm', conversion=us%Z_to_m, &
         cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', &
         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', conversion=us%Z_to_m**2)
     cs%id_MLD_0125 = register_diag_field('ocean_model','MLD_0125',diag%axesT1,time, &
         'Mixed layer depth (delta rho = 0.125)', 'm', conversion=us%Z_to_m)
     cs%id_subMLN2  = register_diag_field('ocean_model','subML_N2',diag%axesT1,time, &
         'Squared buoyancy frequency below mixed layer', 's-2', conversion=us%s_to_T**2)
     cs%id_MLD_user = register_diag_field('ocean_model','MLD_user',diag%axesT1,time, &
         'Mixed layer depth (used defined)', 'm', conversion=us%Z_to_m)
   endif
   call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", cs%MLDdensityDifference, &
                  "The density difference used to determine a diagnostic mixed "//&
                  "layer depth, MLD_user, following the definition of Levitus 1982. "//&
                  "The MLD is the depth at which the density is larger than the "//&
                  "surface density by the specified amount.", &
                  units='kg/m3', default=0.1, scale=us%kg_m3_to_R)
   call get_param(param_file, mdl, "DIAG_DEPTH_SUBML_N2", cs%dz_subML_N2, &
                  "The distance over which to calculate a diagnostic of the "//&
                  "stratification at the base of the mixed layer.", &
                  units='m', default=50.0, scale=us%m_to_Z)
  
   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)
  
   ! diagnostics for values prior to diabatic and prior to ALE
   cs%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, time, &
       'Zonal velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)
   cs%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, time, &
       'Meridional velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)
   cs%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, time, &
       'Layer Thickness before diabatic forcing', trim(thickness_units), &
       conversion=gv%H_to_MKS, v_extensive=.true.)
   cs%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, time, &
       'Interface Heights before diabatic forcing', 'm')
   if (use_temperature) then
     cs%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, time, &
         'Potential Temperature', 'degC')
     cs%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, time, &
         'Salinity', 'PSU')
   endif
  
  
   !call set_diffusivity_init(Time, G, param_file, diag, CS%set_diff_CSp, CS%int_tide_CSp)
   cs%id_Kd_interface = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, time, &
       'Total diapycnal diffusivity at interfaces', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
   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', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
   endif
  
   cs%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, time, &
       'Total diapycnal diffusivity for heat at interfaces', 'm2 s-1', conversion=us%Z2_T_to_m2_s, &
        cmor_field_name='difvho',                                                   &
        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', 'm2 s-1',  conversion=us%Z2_T_to_m2_s, &
        cmor_field_name='difvso',                                                   &
        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, gv, us, 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
  
   if (cs%useKPP .and. differentialdiffusion) then
     call mom_error(fatal, 'diabatic_driver_init: '// &
            'DOUBLE_DIFFUSION (old method) does not work with KPP. Please'//&
            'set DOUBLE_DIFFUSION=False and USE_CVMIX_DDIFF=True.')
   endif
  
   call get_param(param_file, mdl, "SALT_REJECT_BELOW_ML", cs%salt_reject_below_ML, &
                  "If true, place salt from brine rejection below the mixed layer, "// &
                  "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.
   ! diagnostics for tendencies of temp and heat due to frazil
   cs%id_diabatic_diff_h = register_diag_field('ocean_model', 'diabatic_diff_h', diag%axesTL, time, &
       long_name='Cell thickness used during diabatic diffusion', units='m', &
       conversion=gv%H_to_m, v_extensive=.true.)
   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', 'degC s-1', conversion=us%s_to_T)
     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', 'psu s-1', conversion=us%s_to_T)
     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',                                                        &
         'W m-2', conversion=us%s_to_T, cmor_field_name='opottempdiff',                             &
         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 m-2 s-1', conversion=us%s_to_T, cmor_field_name='osaltdiff',                 &
         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',                                     &
         'W m-2', conversion=us%s_to_T, cmor_field_name='opottempdiff_2d',                        &
         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 m-2 s-1', conversion=us%s_to_T, cmor_field_name='osaltdiff_2d',              &
         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 thickness temp and saln due to boundary forcing
     ! available only for ALE algorithm.
   ! diagnostics for tendencies of temp and heat due to frazil
     cs%id_boundary_forcing_h = register_diag_field('ocean_model', 'boundary_forcing_h', diag%axesTL, time, &
       long_name='Cell thickness after applying boundary forcing', units='m', &
       conversion=gv%H_to_m, v_extensive=.true.)
     cs%id_boundary_forcing_h_tendency = register_diag_field('ocean_model',   &
         'boundary_forcing_h_tendency', diag%axesTL, time,                &
         'Cell thickness tendency due to boundary forcing', 'm s-1', conversion=us%s_to_T, &
          v_extensive = .true.)
     if (cs%id_boundary_forcing_h_tendency > 0) then
       cs%boundary_forcing_tendency_diag = .true.
     endif
  
     cs%id_boundary_forcing_temp_tend = register_diag_field('ocean_model',&
         'boundary_forcing_temp_tendency', diag%axesTL, time,             &
         'Boundary forcing temperature tendency', 'degC s-1', conversion=us%s_to_T)
     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', 'psu s-1', conversion=us%s_to_T)
     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', 'W m-2', conversion=us%s_to_T,  &
         v_extensive = .true.)
     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', conversion=us%s_to_T, &
         v_extensive = .true.)
     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', 'W m-2', conversion=us%s_to_T)
     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', conversion=us%s_to_T)
     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
   cs%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, time, &
       long_name='Cell Thickness', standard_name='cell_thickness', units='m', &
       conversion=gv%H_to_m, v_extensive=.true.)
  
   ! 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', 'degC s-1', conversion=us%s_to_T)
   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', 'W m-2', conversion=us%s_to_T, v_extensive=.true.)
   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', 'W m-2', conversion=us%s_to_T)
   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%use_tidal_mixing is set to True if an internal tidal dissipation scheme is to be used.
   cs%use_tidal_mixing = tidal_mixing_init(time, g, gv, us, param_file, diag, &
                                           cs%tidal_mixing_CSp)
  
   ! CS%use_CVMix_conv is set to True if CVMix convection will be used, otherwise
   ! False.
   cs%use_CVMix_conv = cvmix_conv_init(time, g, gv, us, param_file, diag, cs%CVMix_conv_csp)
  
   call entrain_diffusive_init(time, g, gv, us, param_file, diag, cs%entrain_diffusive_CSp)
  
   ! initialize the geothermal heating module
   if (cs%use_geothermal) &
     call geothermal_init(time, g, gv, us, 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, us, param_file, diag, cs%int_tide_input_CSp, &
                              cs%int_tide_input)
     call internal_tides_init(time, g, gv, us, param_file, diag, cs%int_tide_CSp)
   endif
  
   ! initialize module for setting diffusivities
   call set_diffusivity_init(time, g, gv, us, param_file, diag, cs%set_diff_CSp, &
                             cs%int_tide_CSp, cs%tidal_mixing_CSp, cs%halo_TS_diff)
  
   ! 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, us, 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, us, param_file, diag, cs%bulkmixedlayer_CSp)
   if (cs%use_energetic_PBL) &
     call energetic_pbl_init(time, g, gv, us, param_file, diag, cs%energetic_PBL_CSp)
  
   call regularize_layers_init(time, g, gv, param_file, diag, cs%regularize_layers_CSp)
  
   if (cs%debug_energy_req) &
     call diapyc_energy_req_init(time, g, gv, us, 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, mdl, "PEN_SW_NBANDS", nbands, default=1)
     if (nbands > 0) then
       allocate(cs%optics)
       call opacity_init(time, g, gv, us, param_file, diag, cs%opacity_CSp, cs%optics)
     endif
   endif
  
   ! Initialize the diagnostic grid storage
   call diag_grid_storage_init(cs%diag_grids_prev, g, diag)
  

References mom_bulk_mixed_layer::bulkmixedlayer_init(), mom_cvmix_ddiff::cvmix_ddiff_is_used(), mom_cvmix_shear::cvmix_shear_is_used(), mom_geothermal::geothermal_init(), mom_verticalgrid::get_thickness_units(), 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(), and mom_regularize_layers::regularize_layers_init().

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 [H ~> m or kg m-2]
[in]	temp_old	temperature prior to boundary flux application [degC]
[in]	saln_old	salinity prior to boundary flux application [ppt]
[in]	h_old	thickness prior to boundary flux application [H ~> m or kg m-2]
[in]	dt	time step [T ~> s]
	cs	module control structure

Definition at line 3008 of file MOM_diabatic_driver.F90.

   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 [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                            intent(in) :: temp_old !< temperature prior to boundary flux application [degC]
   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 [H ~> m or kg m-2]
   real,                    intent(in) :: dt       !< time step [T ~> s]
   type(diabatic_CS),       pointer    :: CS       !< module control structure
  
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
   real, dimension(SZI_(G),SZJ_(G))         :: work_2d
   real :: Idt  ! The inverse of the timestep [T-1 ~> s-1]
   real :: ppt2mks = 0.001  ! Conversion factor from g/kg to kg/kg.
   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 = 0.0 ; if (dt > 0.0) idt = 1. / dt
   work_3d(:,:,:) = 0.0
   work_2d(:,:)   = 0.0
  
   ! Thickness tendency
   if (cs%id_boundary_forcing_h_tendency > 0) then
     do k=1,nz ; do j=js,je ; do i=is,ie
       work_3d(i,j,k) = (h(i,j,k) - h_old(i,j,k))*idt
     enddo ; enddo ; enddo
     call post_data(cs%id_boundary_forcing_h_tendency, work_3d, cs%diag, alt_h=h_old)
   endif
  
   ! 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, alt_h=h_old)
   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, alt_h = h_old)
     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, alt_h = h_old)
   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 * 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, alt_h = h_old)
     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
  

Referenced by diabatic_ale(), and diabatic_ale_legacy().

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 [H ~> m or kg m-2]
[in]	temp_old	temperature prior to diabatic physics
[in]	saln_old	salinity prior to diabatic physics [ppt]
[in]	dt	time step [T ~> s]
	cs	module control structure

Definition at line 2916 of file MOM_diabatic_driver.F90.

   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 [H ~> m or kg m-2]
   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 [T ~> s]
   type(diabatic_CS),                         pointer    :: CS       !< module control structure
  
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
   real, dimension(SZI_(G),SZJ_(G))         :: work_2d
   real :: Idt  ! The inverse of the timestep [T-1 ~> s-1]
   real :: ppt2mks = 0.001  ! Conversion factor from g/kg to kg/kg.
   integer :: i, j, k, is, ie, js, je, nz
   logical :: do_saln_tend   ! Calculate salinity-based tendency diagnosics
  
   is  = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
   idt = 0.0 ; if (dt > 0.0) 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, alt_h=h)
   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, alt_h=h)
     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_saln_tend = cs%id_diabatic_diff_saln_tend > 0 &
     .or. cs%id_diabatic_diff_salt_tend > 0 &
     .or. cs%id_diabatic_diff_salt_tend_2d > 0
  
   if (do_saln_tend) 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
  
     if (cs%id_diabatic_diff_saln_tend > 0) &
       call post_data(cs%id_diabatic_diff_saln_tend, work_3d, cs%diag, alt_h=h)
  
     ! 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 * 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, alt_h=h)
       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
   endif
  

Referenced by diabatic_ale(), diabatic_ale_legacy(), and layered_diabatic().

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
	)

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
	cs	module control structure
[in]	tv	points to updated thermodynamic fields
[in]	h	thickness [H ~> m or kg m-2]
[in]	temp_old	temperature prior to frazil formation [degC]
[in]	dt	time step [T ~> s]

Definition at line 3104 of file MOM_diabatic_driver.F90.

   type(ocean_grid_type),                    intent(in) :: G        !< ocean grid structure
   type(verticalGrid_type),                  intent(in) :: GV       !< ocean vertical grid structure
   type(diabatic_CS),                        pointer    :: CS       !< module control 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 [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: temp_old !< temperature prior to frazil formation [degC]
   real,                                     intent(in) :: dt       !< time step [T ~> s]
  
   real, dimension(SZI_(G),SZJ_(G))         :: work_2d
   real    :: Idt ! The inverse of the timestep [T-1 ~> s-1]
   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 = 0.0 ; if (dt > 0.0) idt = 1. / dt
  
   ! 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) = idt * (tv%T(i,j,k)-temp_old(i,j,k))
     enddo ; enddo ; enddo
     call post_data(cs%id_frazil_temp_tend, cs%frazil_temp_diag(:,:,:), cs%diag)
   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) = 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) call post_data(cs%id_frazil_heat_tend, cs%frazil_heat_diag(:,:,:), cs%diag)
  
     ! As a consistency check, we must have
     ! FRAZIL_HEAT_TENDENCY_2d = HFSIFRAZIL
     if (cs%id_frazil_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) + 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
  

Referenced by diabatic().

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), optional, pointer	opacity_CSp,
		type(optics_type), optional, pointer	optics_CSp,
		real, intent(out), optional	evap_CFL_limit,
		real, intent(out), optional	minimum_forcing_depth,
		type(kpp_cs), optional, pointer	KPP_CSp,
		type(energetic_pbl_cs), optional, pointer	energetic_PBL_CSp,
		type(diabatic_aux_cs), optional, pointer	diabatic_aux_CSp
	)

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

Parameters

[in]	cs	module control structure
	opacity_csp	A pointer to be set to the opacity control structure
	optics_csp	A pointer to be set to the optics control structure
	kpp_csp	A pointer to be set to the KPP CS
	energetic_pbl_csp	A pointer to be set to the ePBL CS
[out]	evap_cfl_limit	The largest fraction of a layer that can be evaporated in one time-step [nondim].
[out]	minimum_forcing_depth	The smallest depth over which heat and freshwater fluxes are applied [H ~> m or kg m-2].
	diabatic_aux_csp	A pointer to be set to the diabatic_aux control structure

Definition at line 2865 of file MOM_diabatic_driver.F90.

   type(diabatic_CS), intent(in   )           :: CS !< module control structure
   ! All output arguments are optional
   type(opacity_CS),  optional, pointer       :: opacity_CSp !< A pointer to be set to the opacity control structure
   type(optics_type), optional, pointer       :: optics_CSp  !< A pointer to be set to the optics control structure
   type(KPP_CS),      optional, pointer       :: KPP_CSp     !< A pointer to be set to the KPP CS
   type(energetic_PBL_CS), optional, pointer  :: energetic_PBL_CSp !< A pointer to be set to the ePBL CS
   real,              optional, intent(  out) :: evap_CFL_limit !<The largest fraction of a layer that can be
                                                             !! evaporated in one time-step [nondim].
   real,              optional, intent(  out) :: minimum_forcing_depth !< The smallest depth over which heat
                                                             !! and freshwater fluxes are applied [H ~> m or kg m-2].
   type(diabatic_aux_CS), optional, pointer   :: diabatic_aux_CSp !< A pointer to be set to the diabatic_aux
                                                             !! control structure
  
   ! Pointers to control structures
   if (present(opacity_csp))       opacity_csp => cs%opacity_CSp
   if (present(optics_csp))        optics_csp  => cs%optics
   if (present(kpp_csp))           kpp_csp     => cs%KPP_CSp
   if (present(energetic_pbl_csp)) energetic_pbl_csp => cs%energetic_PBL_CSp
  
   ! 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
  

Referenced by mom_lateral_boundary_diffusion::lateral_boundary_diffusion_init(), mom_neutral_diffusion::neutral_diffusion_init(), and mom_offline_main::offline_transport_init().

Here is the caller graph for this function:

◆ layered_diabatic()

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

private

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

Parameters

[in,out]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	u	zonal velocity [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
	cs	module control structure
	waves	Surface gravity waves

Definition at line 1908 of file MOM_diabatic_driver.F90.

   type(ocean_grid_type),                     intent(inout) :: G         !< ocean grid structure
   type(verticalGrid_type),                   intent(in)    :: GV        !< ocean vertical grid structure
   type(unit_scale_type),                     intent(in)    :: US        !< A dimensional unit scaling type
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u         !< zonal velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v         !< meridional velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h         !< thickness [H ~> m or kg m-2]
   type(thermo_var_ptrs),                     intent(inout) :: tv        !< points to thermodynamic fields
                                                                         !! unused have NULL ptrs
   real, dimension(:,:),                      pointer       :: Hml       !< Active mixed layer depth [m]
   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
   type(accel_diag_ptrs),                     intent(inout) :: ADp       !< related points to accelerations in momentum
                                                                         !! equations, to enable the later derived
                                                                         !! diagnostics, like energy budgets
   type(cont_diag_ptrs),                      intent(inout) :: CDp       !< points to terms in continuity equations
   real,                                      intent(in)    :: dt        !< time increment [T ~> s]
   type(time_type),                           intent(in)    :: Time_end  !< Time at the end of the interval
   type(diabatic_CS),                         pointer       :: CS        !< module control structure
   type(Wave_parameters_CS),        optional, pointer       :: Waves     !< Surface gravity waves
  
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
     ea,     &    ! amount of fluid entrained from the layer above within
                  ! one time step  [H ~> m or kg m-2]
     eb,     &    ! amount of fluid entrained from the layer below within
                  ! one time step  [H ~> m or kg m-2]
     kd_lay, &    ! diapycnal diffusivity of layers [Z2 T-1 ~> m2 s-1]
     h_orig, &    ! initial layer thicknesses [H ~> m or kg m-2]
     h_prebound, & ! initial layer thicknesses [H ~> m or kg m-2]
     hold,   &    ! layer thickness before diapycnal entrainment, and later
                  ! the initial layer thicknesses (if a mixed layer is used),
                  ! [H ~> m or kg m-2]
     dsv_dt, &    ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
     dsv_ds, &    ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
     u_h,    &    ! zonal and meridional velocities at thickness points after
     v_h          ! entrainment [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJ_(G)) :: &
     Rcv_ml, &   ! coordinate density of mixed layer, used for applying sponges
     SkinBuoyFlux! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], 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 :: 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 [H ~> m or kg m-2] (for Bous or non-Bouss)
  
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), target :: &
     Kd_int,   & ! diapycnal diffusivity of interfaces [Z2 T-1 ~> m2 s-1]
     Kd_heat,  & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1]
     Kd_salt,  & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1]
     Kd_ePBL,  & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]
     Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
     Tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
     Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
     Sadv_flx    ! advective diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
  
   ! The following 5 variables are only used with a bulk mixed layer.
   real, pointer, dimension(:,:,:) :: &
     eaml, &  ! The equivalent of ea due to mixed layer processes [H ~> m or kg m-2].
     ebml     ! The equivalent of eb due to mixed layer processes [H ~> m or kg m-2].
              ! eaml and ebml are 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 layer [nondim]
  
   real :: p_ref_cv(SZI_(G))      ! Reference pressure for the potential
                                  ! density which defines the coordinate
                                  ! variable, set to P_Ref [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
                        ! [H ~> m or kg m-2]
   real :: h_neglect    ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected
                        ! [H ~> m or kg m-2]
   real :: h_neglect2   ! h_neglect^2 [H2 ~> m2 or kg2 m-4]
   real :: add_ent      ! Entrainment that needs to be added when mixing tracers
                        ! [H ~> m or kg m-2]
   real :: eaval        ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
   real :: hval         ! hval is 2*h at velocity grid points [H ~> m or kg m-2]
   real :: h_tr         ! h_tr is h at tracer points with a tiny thickness
                        ! added to ensure positive definiteness [H ~> m or kg m-2]
   real :: Tr_ea_BBL    ! The diffusive tracer thickness in the BBL that is
                        ! coupled to the bottom within a timestep [H ~> m or kg m-2]
  
   real :: htot(SZIB_(G))             ! The summed thickness from the bottom [H ~> m or kg m-2].
   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 ~> m or kg m-2]
   real :: dt_mix  ! The amount of time over which to apply mixing [T ~> s]
   real :: Idt     ! The inverse time step [T-1 ~> s-1]
   real :: Idt_accel  ! The inverse time step times rescaling factors [T-1 ~> s-1]
  
   integer :: dir_flag     ! An integer encoding the directions in which to do halo updates.
   logical :: showCallTree ! If true, show the call tree
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo
  
   integer :: ig, jg      ! global indices for testing testing itide point source (BDM)
   real :: Kd_add_here    ! An added diffusivity [Z2 T-1 ~> m2 s-1].
  
   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
  
  
   showcalltree = calltree_showquery()
   if (showcalltree) call calltree_enter("layered_diabatic(), MOM_diabatic_driver.F90")
  
   ! set equivalence between the same bits of memory for these arrays
   eaml => eatr ; ebml => ebtr
  
   ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
   call enable_averages(dt, time_end, cs%diag)
  
   if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
     halo = cs%halo_TS_diff
     !$OMP parallel do default(shared)
     do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo
       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, us, cs%geothermal_CSp, halo=cs%halo_TS_diff)
     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, gv, us)
   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_pen_shortwave 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_pen_shortwave(cs%optics, fluxes, g, gv, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
  
   if (cs%bulkmixedlayer) then
     if (cs%debug) call mom_forcing_chksum("Before mixedlayer", fluxes, g, us, haloshift=0)
  
     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, us)
  
       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, us, 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, us, 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, us, 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, us, haloshift=0)
         call mom_forcing_chksum("After mixedlayer", fluxes, g, us, 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, gv, us)
     endif
   endif
  
   if (cs%debug) &
     call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)
   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, us, 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, us)
     endif
     if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
   endif
  
   call cpu_clock_begin(id_clock_set_diffusivity)
   ! Sets: Kd_lay, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S and visc%TKE_turb
   ! Also changes: visc%Kd_shear and visc%Kv_shear
   if ((cs%halo_TS_diff > 0) .and. (cs%ML_mix_first > 0.0)) then
     if (associated(tv%T)) call pass_var(tv%T, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
     if (associated(tv%T)) call pass_var(tv%S, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
     call pass_var(h, g%domain, halo=cs%halo_TS_diff, complete=.true.)
   endif
   call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, &
                        visc, dt, g, gv, us, cs%set_diff_CSp, kd_lay, 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, us, haloshift=0)
     call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)
     call mom_thermovar_chksum("after set_diffusivity ", tv, g)
     call hchksum(kd_lay, "after set_diffusivity Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
   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, us, 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.
  
     ! Set diffusivities for heat and salt separately
  
     !$OMP parallel do default(shared)
     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
     ! Add contribution from double diffusion
     if (associated(visc%Kd_extra_S)) then
       !$OMP parallel do default(shared)
       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 parallel do default(shared)
       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
  
     call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv%eqn_of_state, &
                          fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
  
     call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
                        kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
  
     if (associated(hml)) then
       call kpp_get_bld(cs%KPP_CSp, hml(:,:), g)
       call pass_var(hml, g%domain, halo=1)
       ! If visc%MLD exists, copy KPP's BLD into it
       if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
     endif
  
     if (.not. cs%KPPisPassive) then
       !$OMP parallel do default(shared)
       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 parallel do default(shared)
         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 parallel do default(shared)
         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
  
     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, us, haloshift=0)
       call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)
       call mom_thermovar_chksum("after KPP", tv, g)
       call hchksum(kd_lay, "after KPP Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
       call hchksum(kd_int, "after KPP Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     endif
  
   endif  ! endif for KPP
  
   ! Add vertical diff./visc. due to convection (computed via CVMix)
   if (cs%use_CVMix_conv) then
     call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml)
  
     do k=1,nz+1 ; do j=js,je ; do i=is,ie
       kd_int(i,j,k) = kd_int(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
       visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
     enddo ; enddo ; enddo
   endif
  
   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, &
                                     us%T_to_s*dt, tv%T, tv%C_p)
     call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, &
                                     us%T_to_s*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, gv, us)
  
     if (cs%debug) then
       call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)
       call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, 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 (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, gv, us)
  
     ! increment heat and salt diffusivity.
     ! CS%useKPP==.true. already has extra_T and extra_S included
     if (.not. cs%useKPP) then
       !$OMP parallel do default(shared)
       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.
   ! Otherwise, call entrainment_diffusive() which sets ea and eb
   ! based on KD and target densities (ie. does remapping as well).
   ! 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(h, tv, fluxes, dt, g, gv, us, cs%entrain_diffusive_CSp, &
                              ea, eb, kb, kd_lay=kd_lay, kd_int=kd_int)
   call cpu_clock_end(id_clock_entrain)
   if (showcalltree) call calltree_waypoint("done with Entrainment_diffusive (diabatic)")
  
   if (cs%debug) then
     call mom_forcing_chksum("after calc_entrain ", fluxes, g, us, haloshift=0)
     call mom_thermovar_chksum("after calc_entrain ", tv, g)
     call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, us, 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
  
   ! Save fields before boundary forcing is applied for tendency 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
  
   ! 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(shared)
   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_H
       endif
       if (h(i,j,nz) <= 0.0) then
         h(i,j,nz) = gv%Angstrom_H
       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_H
       endif
     enddo ; enddo
   enddo
   ! Checks for negative thickness may have changed layer thicknesses
   call diag_update_remap_grids(cs%diag)
  
   if (cs%debug) then
     call mom_state_chksum("after negative check ", u, v, h, g, gv, us, haloshift=0)
     call mom_forcing_chksum("after negative check ", fluxes, g, us, 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, gv, us)
  
   ! 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 (shared) 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, gv, us)
  
     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(shared)
       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, us, ea, eb)
       if (cs%debug) call mom_state_chksum("find_uv_at_h1 ", u, v, h, g, gv, us, 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, us, 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, us, 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, gv, us)
     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)
  
       !  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)
  
       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
       ! Note: hold here refers to the thicknesses from before the dual-entraintment when using
       ! the bulk mixed layer scheme, so tendencies should be posted on hold.
       if (cs%diabatic_diff_tendency_diag) then
         call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, cs)
         if (cs%id_diabatic_diff_h > 0) call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h=hold)
       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, gv, us)
  
   endif  ! endif for the BULKMIXEDLAYER block
  
   if (cs%debug) then
     call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, 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
  
   call cpu_clock_begin(id_clock_remap)
   call regularize_layers(h, tv, dt, ea, eb, g, gv, us, 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, gv, us)
  
   ! Whenever thickness changes let the diag manager know, as the
   ! target grids for vertical remapping may need to be regenerated.
   if (cs%id_dudt_dia > 0 .or. cs%id_dvdt_dia > 0) &
     ! Remapped d[uv]dt_dia require east/north halo updates of h
     call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
   call diag_update_remap_grids(cs%diag)
  
   ! diagnostics
   idt = 1.0 / dt
   if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 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(shared)
     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_Sadv > 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(shared)
     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 = gv%Z_to_H * sqrt(dt*cs%Kd_BBL_tr)
     !$OMP parallel do default(shared) 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_lay,
           ! perhaps a molecular diffusivity.
           add_ent = ((dt * cs%Kd_min_tr) * gv%Z_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%Z_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
  
     call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
                               cs%optics, cs%tracer_flow_CSp, cs%debug)
  
   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(shared) 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%Z_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
  
     call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
                                 cs%optics, cs%tracer_flow_CSp, cs%debug)
  
   else
     call call_tracer_column_fns(hold, h, ea, eb, fluxes, hml, dt, g, gv, us, tv, &
                                 cs%optics, cs%tracer_flow_CSp, cs%debug)
  
   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)
     ! 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(shared)
       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, scale=us%kg_m3_to_R)
       enddo
       call apply_sponge(h, dt, g, gv, us, ea, eb, cs%sponge_CSp, rcv_ml)
     else
       call apply_sponge(h, dt, g, gv, us, ea, eb, cs%sponge_CSp)
     endif
     call cpu_clock_end(id_clock_sponge)
     if (cs%debug) then
       call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)
       call mom_thermovar_chksum("apply_sponge ", tv, g)
     endif
   endif ! CS%use_sponge
  
 !   Save the diapycnal mass fluxes as a diagnostic field.
   if (associated(cdp%diapyc_vel)) then
     !$OMP parallel do default(shared)
     do j=js,je
       do k=2,nz ; do i=is,ie
         cdp%diapyc_vel(i,j,k) = us%s_to_T*idt * (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 parallel do default(shared) private(net_ent)
     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 parallel do default(shared)
   do k=1,nz
     do i=is-1,ie+1
       hold(i,js-1,k) = gv%Angstrom_H ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0
       hold(i,je+1,k) = gv%Angstrom_H ; 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_H ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0
       hold(ie+1,j,k) = gv%Angstrom_H ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0
     enddo
   enddo
  
   call cpu_clock_begin(id_clock_pass)
   if (g%symmetric) then ; dir_flag = to_all+omit_corners
   else ; dir_flag = to_west+to_south+omit_corners ; endif
   call create_group_pass(cs%pass_hold_eb_ea, hold, g%Domain, dir_flag, halo=1)
   call create_group_pass(cs%pass_hold_eb_ea, eb, g%Domain, dir_flag, halo=1)
   call create_group_pass(cs%pass_hold_eb_ea, ea, g%Domain, dir_flag, halo=1)
   call do_group_pass(cs%pass_hold_eb_ea, g%Domain)
   call cpu_clock_end(id_clock_pass)
  
   !  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, us, 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)
   idt_accel = 1.0 / dt
   !$OMP parallel do default(shared) 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_accel
     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_accel
       enddo
     endif
   enddo
   if (cs%debug) then
     call mom_state_chksum("aft 1st loop tridiag ", u, v, h, g, gv, us, haloshift=0)
   endif
   !$OMP parallel do default(shared) 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_accel
     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_accel
       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, us, haloshift=0)
   endif
  
   call disable_averaging(cs%diag)
   ! Diagnose the diapycnal diffusivities and other related quantities.
   call enable_averages(dt, time_end, cs%diag)
  
   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_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)
  
   call disable_averaging(cs%diag)
  
   if (showcalltree) call calltree_leave("layered_diabatic()")
  

References mom_sponge::apply_sponge(), mom_bulk_mixed_layer::bulkmixedlayer(), mom_cvmix_conv::calculate_cvmix_conv(), 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_diabatic_diff_tendency(), mom_domains::do_group_pass(), 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_cvmix_kpp::kpp_nonlocaltransport_saln(), mom_cvmix_kpp::kpp_nonlocaltransport_temp(), mom_checksum_packages::mom_state_stats(), mom_variables::mom_thermovar_chksum(), mom_regularize_layers::regularize_layers(), mom_diabatic_aux::set_pen_shortwave(), and mom_tracer_diabatic::tracer_vertdiff().

Referenced by diabatic().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

◆ id_clock_differential_diff

integer mom_diabatic_driver::id_clock_differential_diff

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 250 of file MOM_diabatic_driver.F90.

Referenced by diabatic_ale(), diabatic_ale_legacy(), diabatic_driver_init(), and layered_diabatic().

◆ id_clock_entrain

integer mom_diabatic_driver::id_clock_entrain

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 248 of file MOM_diabatic_driver.F90.

248 integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity

Referenced by diabatic_driver_init(), and layered_diabatic().

◆ id_clock_geothermal

integer mom_diabatic_driver::id_clock_geothermal

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 250 of file MOM_diabatic_driver.F90.

250 integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap

Referenced by diabatic_ale(), diabatic_ale_legacy(), diabatic_driver_init(), and layered_diabatic().

◆ id_clock_kpp

integer mom_diabatic_driver::id_clock_kpp

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 251 of file MOM_diabatic_driver.F90.

251 integer :: id_clock_kpp

Referenced by diabatic_ale(), diabatic_ale_legacy(), diabatic_driver_init(), and layered_diabatic().

◆ id_clock_mixedlayer

integer mom_diabatic_driver::id_clock_mixedlayer

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 248 of file MOM_diabatic_driver.F90.

Referenced by diabatic_driver_init(), and layered_diabatic().

◆ id_clock_pass

integer mom_diabatic_driver::id_clock_pass

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 249 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), diabatic_ale(), diabatic_driver_init(), and layered_diabatic().

◆ id_clock_remap

integer mom_diabatic_driver::id_clock_remap

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 250 of file MOM_diabatic_driver.F90.

Referenced by diabatic_ale(), diabatic_ale_legacy(), diabatic_driver_init(), and layered_diabatic().

◆ id_clock_set_diffusivity

integer mom_diabatic_driver::id_clock_set_diffusivity

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 248 of file MOM_diabatic_driver.F90.

Referenced by diabatic(), diabatic_ale(), diabatic_ale_legacy(), diabatic_driver_init(), and layered_diabatic().

◆ id_clock_sponge

integer mom_diabatic_driver::id_clock_sponge

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 249 of file MOM_diabatic_driver.F90.

Referenced by diabatic_ale(), diabatic_ale_legacy(), diabatic_driver_init(), and layered_diabatic().

◆ id_clock_tracers

integer mom_diabatic_driver::id_clock_tracers

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 249 of file MOM_diabatic_driver.F90.

249 integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge

Referenced by diabatic_ale(), diabatic_ale_legacy(), diabatic_driver_init(), and layered_diabatic().

◆ id_clock_tridiag

integer mom_diabatic_driver::id_clock_tridiag

private

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 [L T-1 ~> m s-1]
[in,out]	v	meridional velocity [L T-1 ~> m s-1]
[in,out]	h	thickness [H ~> m or kg m-2]
[in,out]	tv	points to thermodynamic fields unused have NULL ptrs
	hml	Active mixed layer depth [m]
[in,out]	fluxes	points to forcing fields unused fields have NULL ptrs
[in,out]	visc	vertical viscosities, BBL properies, and
[in,out]	adp	related points to accelerations in momentum equations, to enable the later derived diagnostics, like energy budgets
[in,out]	cdp	points to terms in continuity equations
[in]	dt	time increment [T ~> s]
[in]	time_end	Time at the end of the interval
[in]	us	A dimensional unit scaling type
	cs	module control structure
	waves	Surface gravity waves

Definition at line 249 of file MOM_diabatic_driver.F90.

Referenced by diabatic_ale(), diabatic_ale_legacy(), diabatic_driver_init(), and layered_diabatic().

Detailed Description

Outline of MOM diabatic

Data Types

Function/Subroutine Documentation

◆ adiabatic()

◆ adiabatic_driver_init()

◆ diabatic()

◆ diabatic_ale()

◆ diabatic_ale_legacy()

◆ diabatic_driver_end()

◆ diabatic_driver_init()

◆ diagnose_boundary_forcing_tendency()

◆ diagnose_diabatic_diff_tendency()

◆ diagnose_frazil_tendency()

◆ extract_diabatic_member()

◆ layered_diabatic()

Variable Documentation

◆ 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