Skip to content

Commit

Permalink
+Add runtime parameter SALT_EXTRACTION_LIMIT
Browse files Browse the repository at this point in the history 
  Added the new runtime parameter SALT_EXTRACTION_LIMIT to specify a previously
hard-coded limit on the fraction of the salt in a layer that can be extracted by
the surface fluxes within a timestep.  Also added do_not_log flags based on the
value of ENABLE_THERMODYNAMICS to avoid logging parameters that are only used if
thermodynamics are active in cases when it is not.  In addition, the doxygen
trailer describing this module had been inappropriately borrowed from another
module, so it was completely rewritten to describe what the routines in this
module actually do.  The missing doxygen description of set_pen_shortwave was
also added, and a number of spelling errors in comments were corrected.  By
default all answers are bitwise identical, but there is a new entry in the
MOM_parameter_doc.all files for some configurations and meaningless entries have
been removed from others.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Jan 3, 2023
1 parent 994c29c commit f388abb
Showing 1 changed file with 76 additions and 56 deletions.
132 changes: 76 additions & 56 deletions src/parameterizations/vertical/MOM_diabatic_aux.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,8 @@ module MOM_diabatic_aux
logical :: do_rivermix = .false. !< Provide additional TKE to mix river runoff at the
!! river mouths to a depth of "rivermix_depth"
real :: rivermix_depth = 0.0 !< The depth to which rivers are mixed if do_rivermix = T [Z ~> m].
real :: dSalt_frac_max !< An upper limit on the fraction of the salt in a layer that can be
!! lost to the net surface salt fluxes within a timestep [nondim]
logical :: reclaim_frazil !< If true, try to use any frazil heat deficit to
!! to cool the topmost layer down to the freezing
!! point. The default is true.
Expand Down Expand Up @@ -220,7 +222,7 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo)
end subroutine make_frazil

!> This subroutine applies double diffusion to T & S, assuming no diapycnal mass
!! fluxes, using a simple triadiagonal solver.
!! fluxes, using a simple tridiagonal solver.
subroutine differential_diffuse_T_S(h, T, S, Kd_T, Kd_S, dt, G, GV)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
Expand Down Expand Up @@ -507,7 +509,7 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb, zero_mix)
real :: a_n(SZI_(G)), a_s(SZI_(G)) ! Fractional weights of the neighboring velocity points [nondim]
real :: a_e(SZI_(G)), a_w(SZI_(G)) ! Fractional weights of the neighboring velocity points [nondim]
real :: sum_area ! A sum of adjacent areas [L2 ~> m2]
real :: Idenom ! The inverse of the denomninator in a weighted average [L-2 ~> m-2]
real :: Idenom ! The inverse of the denominator in a weighted average [L-2 ~> m-2]
logical :: mix_vertically, zero_mixing
integer :: i, j, k, is, ie, js, je, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
Expand Down Expand Up @@ -594,12 +596,15 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb, zero_mix)
call cpu_clock_end(id_clock_uv_at_h)
end subroutine find_uv_at_h


!> Estimate the optical properties of the water column and determine the penetrating shortwave
!! radiation by band, extracting the relevant information from the fluxes type and storing it
!! in the optics type for later application. This routine is effectively a wrapper for
!! set_opacity with added error handling and diagnostics.
subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity, tracer_flow_CSp)
type(optics_type), pointer :: optics !< An optics structure that has will contain
!! information about shortwave fluxes and absorption.
type(forcing), intent(inout) :: fluxes !< points to forcing fields
!! unused fields have NULL ptrs
!! unused fields have NULL pointers
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
Expand Down Expand Up @@ -702,7 +707,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
id_N2 = id_N2subML
dH_subML = GV%Z_to_H*dz_subML
else
call MOM_error(FATAL, "When the a diagnostic of the subML stratification is "//&
call MOM_error(FATAL, "When the diagnostic of the subML stratification is "//&
"requested by providing id_N2_subML to diagnoseMLDbyDensityDifference, "//&
"the distance over which to calculate that distance must also be provided.")
endif
Expand Down Expand Up @@ -737,10 +742,10 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
! the cells that extend over at least dz_subML.
if (id_N2>0) then
do i=is,ie
if (MLD(i,j)==0.0) then ! Still in the mixed layer.
if (MLD(i,j) == 0.0) then ! Still in the mixed layer.
H_subML(i) = H_subML(i) + h(i,j,k)
elseif (.not.N2_region_set(i)) then ! This block is below the mixed layer, but N2 has not been found yet.
if (dH_N2(i)==0.0) then ! Record the temperature, salinity, pressure, immediately below the ML
if (dH_N2(i) == 0.0) then ! Record the temperature, salinity, pressure, immediately below the ML
T_subML(i) = tv%T(i,j,k) ; S_subML(i) = tv%S(i,j,k)
H_subML(i) = H_subML(i) + 0.5 * h(i,j,k) ! Start midway through this layer.
dH_N2(i) = 0.5 * h(i,j,k)
Expand All @@ -761,16 +766,16 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
do i = is, ie
deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface
ddRho = deltaRhoAtK(i) - deltaRhoAtKm1(i)
if ((MLD(i,j)==0.) .and. (ddRho>0.) .and. &
(deltaRhoAtKm1(i)<densityDiff) .and. (deltaRhoAtK(i)>=densityDiff)) then
if ((MLD(i,j) == 0.) .and. (ddRho > 0.) .and. &
(deltaRhoAtKm1(i) < densityDiff) .and. (deltaRhoAtK(i) >= densityDiff)) then
aFac = ( densityDiff - deltaRhoAtKm1(i) ) / ddRho
MLD(i,j) = dK(i) * aFac + dKm1(i) * (1. - aFac)
endif
if (id_SQ > 0) MLD2(i,j) = MLD(i,j)**2
enddo ! i-loop
enddo ! k-loop
do i=is,ie
if ((MLD(i,j)==0.) .and. (deltaRhoAtK(i)<densityDiff)) MLD(i,j) = dK(i) ! Assume mixing to the bottom
if ((MLD(i,j) == 0.) .and. (deltaRhoAtK(i) < densityDiff)) MLD(i,j) = dK(i) ! Assume mixing to the bottom
enddo

if (id_N2>0) then ! Now actually calculate stratification, N2, below the mixed layer.
Expand Down Expand Up @@ -1293,7 +1298,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
! B/ update mass, salt, temp from mass leaving ocean.
! C/ update temp due to penetrative SW
do i=is,ie
if (G%mask2dT(i,j)>0.) then
if (G%mask2dT(i,j) > 0.) then

! A/ Update mass, temp, and salinity due to incoming mass flux.
do k=1,1
Expand Down Expand Up @@ -1333,7 +1338,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
! where River is in units of [Z T-1 ~> m s-1].
! Samb = Ambient salinity at the mouth of the estuary
! rivermix_depth = The prescribed depth over which to mix river inflow
! drho_ds = The gradient of density wrt salt at the ambient surface salinity.
! drho_ds = The derivative of density with salt at the ambient surface salinity.
! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
if (GV%Boussinesq) then
RivermixConst = -0.5*(CS%rivermix_depth*dt) * ( US%L_to_Z**2*GV%g_Earth ) * GV%Rho0
Expand Down Expand Up @@ -1384,8 +1389,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t

dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) )
dTemp = fractionOfForcing*netHeat(i)
! ### The 0.9999 here should become a run-time parameter?
dSalt = max( fractionOfForcing*netSalt(i), -0.9999*h2d(i,k)*tv%S(i,j,k))
dSalt = max( fractionOfForcing*netSalt(i), -CS%dSalt_frac_max * h2d(i,k) * tv%S(i,j,k))

! Update the forcing by the part to be consumed within the present k-layer.
! If fractionOfForcing = 1, then new netMassOut vanishes.
Expand Down Expand Up @@ -1441,7 +1445,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
enddo ! k

! Check if trying to apply fluxes over land points
elseif ((abs(netHeat(i))+abs(netSalt(i))+abs(netMassIn(i))+abs(netMassOut(i)))>0.) then
elseif ((abs(netHeat(i)) + abs(netSalt(i)) + abs(netMassIn(i)) + abs(netMassOut(i))) > 0.) then

if (.not. CS%ignore_fluxes_over_land) then
call forcing_SinglePointPrint(fluxes,G,i,j,'applyBoundaryFluxesInOut (land)')
Expand Down Expand Up @@ -1579,7 +1583,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
if (CS%id_nonpenSW_diag > 0) call post_data(CS%id_nonpenSW_diag , CS%nonpenSW_diag , CS%diag)

! The following check will be ignored if ignore_fluxes_over_land = true
if (numberOfGroundings>0 .and. .not. CS%ignore_fluxes_over_land) then
if ((numberOfGroundings > 0) .and. .not.CS%ignore_fluxes_over_land) then
do i = 1, min(numberOfGroundings, maxGroundings)
call forcing_SinglePointPrint(fluxes,G,iGround(i),jGround(i),'applyBoundaryFluxesInOut (grounding)')
write(mesg(1:45),'(3es15.3)') G%geoLonT( iGround(i), jGround(i) ), &
Expand Down Expand Up @@ -1613,8 +1617,8 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori
!! boundary layer scheme to determine the diffusivity
!! in the surface boundary layer.

! This "include" declares and sets the variable "version".
#include "version_variable.h"
! This "include" declares and sets the variable "version".
# include "version_variable.h"
character(len=40) :: mdl = "MOM_diabatic_aux" ! This module's name.
character(len=200) :: inputdir ! The directory where NetCDF input files
character(len=240) :: chl_filename ! A file from which chl_a concentrations are to be read.
Expand Down Expand Up @@ -1642,28 +1646,31 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori
"The following parameters are used for auxiliary diabatic processes.")

call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
"If true, temperature and salinity are used as state "//&
"variables.", default=.true.)
"If true, temperature and salinity are used as state variables.", default=.true.)

call get_param(param_file, mdl, "RECLAIM_FRAZIL", CS%reclaim_frazil, &
"If true, try to use any frazil heat deficit to cool any "//&
"overlying layers down to the freezing point, thereby "//&
"avoiding the creation of thin ice when the SST is above "//&
"the freezing point.", default=.true.)
call get_param(param_file, mdl, "PRESSURE_DEPENDENT_FRAZIL", &
CS%pressure_dependent_frazil, &
"the freezing point.", default=.true., do_not_log=.not.use_temperature)
call get_param(param_file, mdl, "SALT_EXTRACTION_LIMIT", CS%dSalt_frac_max, &
"An upper limit on the fraction of the salt in a layer that can be lost to the "//&
"net surface salt fluxes within a timestep.", &
units="nondim", default=0.9999, do_not_log=.not.use_temperature)
CS%dSalt_frac_max = max(min(CS%dSalt_frac_max, 1.0), 0.0)
call get_param(param_file, mdl, "PRESSURE_DEPENDENT_FRAZIL", CS%pressure_dependent_frazil, &
"If true, use a pressure dependent freezing temperature "//&
"when making frazil. The default is false, which will be "//&
"faster but is inappropriate with ice-shelf cavities.", &
default=.false.)
default=.false., do_not_log=.not.use_temperature)

if (use_ePBL) then
call get_param(param_file, mdl, "IGNORE_FLUXES_OVER_LAND", CS%ignore_fluxes_over_land,&
"If true, the model does not check if fluxes are being applied "//&
"over land points. This is needed when the ocean is coupled "//&
"with ice shelves and sea ice, since the sea ice mask needs to "//&
"be different than the ocean mask to avoid sea ice formation "//&
"under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
"If true, the model does not check if fluxes are being applied "//&
"over land points. This is needed when the ocean is coupled "//&
"with ice shelves and sea ice, since the sea ice mask needs to "//&
"be different than the ocean mask to avoid sea ice formation "//&
"under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
call get_param(param_file, mdl, "DO_RIVERMIX", CS%do_rivermix, &
"If true, apply additional mixing wherever there is "//&
"runoff, so that it is mixed down to RIVERMIX_DEPTH "//&
Expand All @@ -1680,11 +1687,11 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori
call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", CS%use_river_heat_content, &
"If true, use the fluxes%runoff_Hflx field to set the "//&
"heat carried by runoff, instead of using SST*CP*liq_runoff.", &
default=.false.)
default=.false., do_not_log=.not.use_temperature)
call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", CS%use_calving_heat_content, &
"If true, use the fluxes%calving_Hflx field to set the "//&
"heat carried by runoff, instead of using SST*CP*froz_runoff.", &
default=.false.)
default=.false., do_not_log=.not.use_temperature)
else
CS%use_river_heat_content = .false.
CS%use_calving_heat_content = .false.
Expand Down Expand Up @@ -1778,32 +1785,45 @@ end subroutine diabatic_aux_end

!> \namespace mom_diabatic_aux
!!
!! This module contains the subroutines 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.
!! This module contains subroutines that apply various diabatic processes. Usually these
!! subroutines are called from the MOM_diabatic module. All of these routines use appropriate
!! limiters or logic to work properly with arbitrary layer thicknesses (including massless layers)
!! and an arbitrarily large timestep.
!!
!! The subroutine make_frazil facilitates the formation of frazil ice when the ocean water
!! drops below the in situ freezing point by heating the water to the freezing point and
!! accumulating the required heat for exchange with the sea-ice module.
!!
!! The subroutine adjust_salt adds salt as necessary to keep the salinity above a
!! specified minimum value, and keeps track of the cumulative additions. If the minimum
!! salinity is the natural value of 0, this routine should never do anything.
!!
!! The subroutine differential_diffuse_T_S solves a pair of tridiagonal equations for
!! the diffusion of temperatures and salinities with differing diffusivities.
!!
!! The subroutine triDiagTS solves a tridiagonal equations for the evolution of temperatures
!! and salinities due to net entrainment by layers and a diffusion with the same diffusivity.
!!
!! The subroutine triDiagTS_Eulerian solves a tridiagonal equations for the evolution of
!! temperatures and salinities due to diffusion with the same diffusivity, but no net entrainment.
!!
!! The subroutine find_uv_at_h interpolates velocities to thickness points, optionally also
!! using tridiagonal equations to solve for the impacts of net entrainment or mixing of
!! momentum between layers.
!!
!! The subroutine set_pen_shortwave determines the optical properties of the water column and
!! the net shortwave fluxes, and stores them in the optics type, working via calls to set_opacity.
!!
!! The subroutine diagnoseMLDbyDensityDifference diagnoses a mixed layer depth based on a
!! density difference criterion, and may also estimate the stratification of the water below
!! this diagnosed mixed layer.
!!
!! 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 fundamentally 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 treat-
!! ed as a fixed density layer with vanishingly small diffusivity.
!! The subroutine diagnoseMLDbyEnergy diagnoses a mixed layer depth based on a mixing-energy
!! criterion, as described by Reichl et al., 2022, JGR: Oceans, doi:10.1029/2021JC018140.
!!
!! 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.
!! The subroutine applyBoundaryFluxesInOut updates the layer thicknesses, temperatures and
!! salinities due to the application of the surface forcing. It may also calculate the implied
!! turbulent kinetic energy requirements for this forcing to be mixed over the model's finite
!! vertical resolution in the surface layers.

end module MOM_diabatic_aux

0 comments on commit f388abb

Please sign in to comment.