mom_diabatic_aux Module Reference

Data Types
type	diabatic_aux_cs
	Control structure for diabatic_aux. More...

Functions/Subroutines
subroutine, public	make_frazil (h, tv, G, GV, CS, p_surf)
subroutine, public	differential_diffuse_t_s (h, tv, visc, dt, G, GV)
subroutine, public	adjust_salt (h, tv, G, GV, CS)
subroutine, public	insert_brine (h, tv, G, GV, fluxes, nkmb, CS, dt, id_brine_lay)
subroutine, public	tridiagts (G, GV, is, ie, js, je, hold, ea, eb, T, S)
subroutine, public	find_uv_at_h (u, v, h, u_h, v_h, G, GV, ea, eb)
subroutine, public	diagnosemldbydensitydifference (id_MLD, h, tv, densityDiff, G, GV, diagPtr, id_N2subML, id_MLDsq)
	Diagnose a mixed layer depth (MLD) determined by a given density difference with the surface. This routine is appropriate in MOM_diabatic_driver due to its position within the time stepping. More...
subroutine, public	applyboundaryfluxesinout (CS, G, GV, dt, fluxes, optics, h, tv, aggregate_FW_forcing, evap_CFL_limit, minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux)
	Update the thickness, temperature, and salinity due to thermodynamic boundary forcing (contained in fluxes type) applied to h, tvT and tvS, and calculate the TKE implications of this heating. More...
subroutine, public	diabatic_aux_init (Time, G, GV, param_file, diag, CS, useALEalgorithm, use_ePBL)
subroutine, public	diabatic_aux_end (CS)

Variables
integer	id_clock_uv_at_h
integer	id_clock_frazil

Function/Subroutine Documentation

adjust_salt()

subroutine, public mom_diabatic_aux::adjust_salt	(	real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		type(thermo_var_ptrs), intent(inout)	tv,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(diabatic_aux_cs), intent(in)	CS
	)

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	h	Layer thicknesses, in H (usually m or kg m-2)

Definition at line 368 of file MOM_diabatic_aux.F90.

  type(ocean_grid_type),                 intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),               intent(in)    :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                 intent(inout) :: tv
  type(diabatic_aux_cs),                 intent(in)    :: cs

!  Keep salinity from falling below a small but positive threshold
!  This occurs when the ice model attempts to extract more salt then
!  is actually available to it from the ocean.

! Arguments: h - Layer thickness, in m.
!  (in/out)  tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 diabatic_driver_init.
  real :: salt_add_col(szi_(g),szj_(g)) ! The accumulated salt requirement
  real :: s_min      ! The minimum salinity
  real :: mc         ! A layer's mass kg  m-2 .
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

!  call cpu_clock_begin(id_clock_adjust_salt)

  s_min = 0.01

  salt_add_col(:,:) = 0.0

!$OMP parallel do default(none) shared(is,ie,js,je,nz,G,GV,tv,h,salt_add_col, S_min) &
!$OMP                          private(mc)
  do j=js,je
    do k=nz,1,-1 ; do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .and. &
           ((tv%S(i,j,k) < s_min) .or. (salt_add_col(i,j) > 0.0))) then
        mc = gv%H_to_kg_m2 * h(i,j,k)
        if (h(i,j,k) <= 10.0*gv%Angstrom) then
          ! Very thin layers should not be adjusted by the salt flux
          if (tv%S(i,j,k) < s_min) then
            salt_add_col(i,j) = salt_add_col(i,j) +  mc * (s_min - tv%S(i,j,k))
            tv%S(i,j,k) = s_min
          endif
        else
          if (salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k)) <= 0.0) then
            tv%S(i,j,k) = tv%S(i,j,k) - salt_add_col(i,j)/mc
            salt_add_col(i,j) = 0.0
          else
            salt_add_col(i,j) = salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k))
            tv%S(i,j,k) = s_min
          endif
        endif
      endif
    enddo ; enddo
    do i=is,ie
      tv%salt_deficit(i,j) = tv%salt_deficit(i,j) + salt_add_col(i,j)
    enddo
  enddo
!  call cpu_clock_end(id_clock_adjust_salt)

applyboundaryfluxesinout()

subroutine, public mom_diabatic_aux::applyboundaryfluxesinout	(	type(diabatic_aux_cs), pointer	CS,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, intent(in)	dt,
		type(forcing), intent(inout)	fluxes,
		type(optics_type), pointer	optics,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	h,
		type(thermo_var_ptrs), intent(inout)	tv,
		logical, intent(in)	aggregate_FW_forcing,
		real, intent(in)	evap_CFL_limit,
		real, intent(in)	minimum_forcing_depth,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(out), optional	cTKE,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(out), optional	dSV_dT,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(out), optional	dSV_dS,
		real, dimension(szi_(g),szj_(g)), intent(out), optional	SkinBuoyFlux
	)

Update the thickness, temperature, and salinity due to thermodynamic boundary forcing (contained in fluxes type) applied to h, tvT and tvS, and calculate the TKE implications of this heating.

Parameters

	cs	Control structure for diabatic_aux
[in]	g	Grid structure
[in]	gv	ocean vertical grid structure
[in]	dt	Time-step over which forcing is applied (s)
[in,out]	fluxes	Surface fluxes container
	optics	Optical properties container
[in,out]	h	Layer thickness in H units
[in,out]	tv	Thermodynamics container
[in]	aggregate_fw_forcing	If False, treat in/out fluxes separately.
[in]	evap_cfl_limit	The largest fraction of a layer that can be evaporated in one time-step (non-dim).
[in]	minimum_forcing_depth	The smallest depth over which heat and freshwater fluxes is applied, in m.
[out]	ctke	Turbulent kinetic energy requirement to mix forcing through each layer, in W m-2
[out]	dsv_dt	Partial derivative of specific volume with potential temperature, in m3 kg-1 K-1.
[out]	dsv_ds	Partial derivative of specific a volume with potential salinity, in m3 kg-1 / (g kg-1).
[out]	skinbuoyflux	Buoyancy flux at surface in m2 s-3

Definition at line 814 of file MOM_diabatic_aux.F90.

References mom_shortwave_abs::absorbremainingsw(), mom_eos::calculate_density_derivs(), mom_eos::calculate_specific_vol_derivs(), mom_forcing_type::extractfluxes1d(), mom_forcing_type::forcing_singlepointprint(), and mom_shortwave_abs::sumswoverbands().

Referenced by mom_diabatic_driver::diabatic().

  type(diabatic_aux_cs),                 pointer       :: cs !< Control structure for diabatic_aux
  type(ocean_grid_type),                 intent(in)    :: g  !< Grid structure
  type(verticalgrid_type),               intent(in)    :: gv !< ocean vertical grid structure
  real,                                  intent(in)    :: dt !< Time-step over which forcing is applied (s)
  type(forcing),                         intent(inout) :: fluxes !< Surface fluxes container
  type(optics_type),                     pointer       :: optics !< Optical properties container
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h  !< Layer thickness in H units
  type(thermo_var_ptrs),                 intent(inout) :: tv !< Thermodynamics container
  !> If False, treat in/out fluxes separately.
  logical,                               intent(in)    :: aggregate_fw_forcing
  !> The largest fraction of a layer that can be evaporated in one time-step (non-dim).
  real,                                  intent(in)   :: evap_cfl_limit
  !> The smallest depth over which heat and freshwater fluxes is applied, in m.
  real,                                  intent(in)   :: minimum_forcing_depth
  !> Turbulent kinetic energy requirement to mix forcing through each layer, in W m-2
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: ctke
  !> Partial derivative of specific volume with potential temperature, in m3 kg-1 K-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: dsv_dt
  !> Partial derivative of specific a volume with potential salinity, in m3 kg-1 / (g kg-1).
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: dsv_ds
  !> Buoyancy flux at surface in m2 s-3
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: skinbuoyflux

  ! Local variables
  integer, parameter :: maxgroundings = 5
  integer :: numberofgroundings, iground(maxgroundings), jground(maxgroundings)
  real :: h_limit_fluxes, iforcingdepthscale, idt
  real :: dthickness, dtemp, dsalt
  real :: fractionofforcing, hold, ithickness
  real :: rivermixconst  ! A constant used in implementing river mixing, in Pa s.
  real, dimension(SZI_(G)) :: &
    d_pres,       &  ! pressure change across a layer (Pa)
    p_lay,        &  ! average pressure in a layer (Pa)
    pres,         &  ! pressure at an interface (Pa)
    netmassinout, &  ! surface water fluxes (H units) over time step
    netmassin,    &  ! mass entering ocean surface (H units) over a time step
    netmassout,   &  ! mass leaving ocean surface (H units) over a time step
    netheat,      &  ! heat (degC * H) via surface fluxes, excluding
                     ! Pen_SW_bnd and netMassOut
    netsalt,      &  ! surface salt flux ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
    nonpensw,     &  ! non-downwelling SW, which is absorbed at ocean surface
    surfpressure, &  ! Surface pressure (approximated as 0.0)
    drhodt,       &  ! change in density per change in temperature
    drhods,       &  ! change in density per change in salinity
    netheat_rate, &  ! netheat but for dt=1 (e.g. returns a rate)
    netsalt_rate, &  ! netsalt but for dt=1 (e.g. returns a rate)
    netmassinout_rate! netmassinout but for dt=1 (e.g. returns a rate)
  real, dimension(SZI_(G), SZK_(G))                     :: h2d, t2d
  real, dimension(SZI_(G), SZK_(G))                     :: pen_tke_2d, dsv_dt_2d
  real, dimension(SZI_(G),SZK_(G)+1)                    :: netpen
  real, dimension(max(optics%nbands,1),SZI_(G))         :: pen_sw_bnd, pen_sw_bnd_rate
                                                           !^ _rate is w/ dt=1
  real, dimension(max(optics%nbands,1),SZI_(G),SZK_(G)) :: opacityband
  real                                                  :: hgrounding(maxgroundings)
  real    :: temp_in, salin_in
  real    :: i_g_earth, g_hconv2
  real    :: gorho
  logical :: calculate_energetics
  logical :: calculate_buoyancy
  integer :: i, j, is, ie, js, je, k, nz, n, nsw
  integer :: start, npts
  character(len=45) :: mesg

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  ! Only apply forcing if fluxes%sw is associated.
  if (.not.ASSOCIATED(fluxes%sw)) return

#define _OLD_ALG_
  nsw = optics%nbands
  idt = 1.0/dt

  calculate_energetics = (present(ctke) .and. present(dsv_dt) .and. present(dsv_ds))
  calculate_buoyancy = present(skinbuoyflux)
  if (calculate_buoyancy) skinbuoyflux(:,:) = 0.0
  i_g_earth = 1.0 / gv%g_Earth
  g_hconv2 = gv%g_Earth * gv%H_to_kg_m2**2

  if (present(ctke)) ctke(:,:,:) = 0.0
  if (calculate_buoyancy) then
    surfpressure(:) = 0.0
    gorho       = gv%g_Earth / gv%Rho0
    start       = 1 + g%isc - g%isd
    npts        = 1 + g%iec - g%isc
  endif

  ! H_limit_fluxes is used by extractFluxes1d to scale down fluxes if the total
  ! depth of the ocean is vanishing. It does not (yet) handle a value of zero.
  ! To accommodate vanishing upper layers, we need to allow for an instantaneous
  ! distribution of forcing over some finite vertical extent. The bulk mixed layer
  ! code handles this issue properly.
  h_limit_fluxes = max(gv%Angstrom, 1.e-30*gv%m_to_H)

  ! diagnostic to see if need to create mass to avoid grounding
  if (cs%id_createdH>0) cs%createdH(:,:) = 0.
  numberofgroundings = 0

!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,tv,nsw,G,GV,optics,fluxes,dt,    &
!$OMP                                  H_limit_fluxes,IforcingDepthScale,                &
!$OMP                                  numberOfGroundings,iGround,jGround,nonPenSW,      &
!$OMP                                  hGrounding,CS,Idt,aggregate_FW_forcing,           &
!$OMP                                  minimum_forcing_depth,evap_CFL_limit,             &
!$OMP                                  calculate_buoyancy,netPen,SkinBuoyFlux,GoRho,     &
!$OMP                                  calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2) &
!$OMP                          private(opacityBand,h2d,T2d,netMassInOut,netMassOut,      &
!$OMP                                  netHeat,netSalt,Pen_SW_bnd,fractionOfForcing,     &
!$OMP                                  dThickness,dTemp,dSalt,hOld,Ithickness,           &
!$OMP                                  netMassIn,pres,d_pres,p_lay,dSV_dT_2d,            &
!$OMP                                  netmassinout_rate,netheat_rate,netsalt_rate,      &
!$OMP                                  drhodt,drhods,pen_sw_bnd_rate,SurfPressure,       &
!$OMP                                  start,npts,                                       &
!$OMP                                  pen_TKE_2d,Temp_in,Salin_in,RivermixConst)

  ! Work in vertical slices for efficiency
  do j=js,je

    ! Copy state into 2D-slice arrays
    do k=1,nz ; do i=is,ie
      h2d(i,k) = h(i,j,k)
      t2d(i,k) = tv%T(i,j,k)
      do n=1,nsw
        opacityband(n,i,k) = (1.0 / gv%m_to_H)*optics%opacity_band(n,i,j,k)
      enddo
    enddo ; enddo

    if (calculate_energetics) then
      ! The partial derivatives of specific volume with temperature and
      ! salinity need to be precalculated to avoid having heating of
      ! tiny layers give nonsensical values.
      do i=is,ie ; pres(i) = 0.0 ; enddo ! Add surface pressure?
      do k=1,nz
        do i=is,ie
          d_pres(i) = gv%g_Earth * gv%H_to_kg_m2 * h2d(i,k)
          p_lay(i) = pres(i) + 0.5*d_pres(i)
          pres(i) = pres(i) + d_pres(i)
        enddo
        call calculate_specific_vol_derivs(t2d(:,k), tv%S(:,j,k), p_lay(:),&
                 dsv_dt(:,j,k), dsv_ds(:,j,k), is, ie-is+1, tv%eqn_of_state)
        do i=is,ie ; dsv_dt_2d(i,k) = dsv_dt(i,j,k) ; enddo
!        do i=is,ie
!          dT_to_dPE(i,k) = I_G_Earth * d_pres(i) * p_lay(i) * dSV_dT(i,j,k)
!          dS_to_dPE(i,k) = I_G_Earth * d_pres(i) * p_lay(i) * dSV_dS(i,j,k)
!        enddo
      enddo
      pen_tke_2d(:,:) = 0.0
    endif

    ! The surface forcing is contained in the fluxes type.
    ! We aggregate the thermodynamic forcing for a time step into the following:
    ! netMassInOut = surface water fluxes (H units) over time step
    !              = lprec + fprec + vprec + evap + lrunoff + frunoff
    !                note that lprec generally has sea ice melt/form included.
    ! netMassOut   = net mass leaving ocean surface (H units) over a time step.
    !                netMassOut < 0 means mass leaves ocean.
    ! netHeat      = heat (degC * H) via surface fluxes, excluding the part
    !                contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0.
    ! netSalt      = surface salt fluxes ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
    ! Pen_SW_bnd   = components to penetrative shortwave radiation split according to bands.
    !                This field provides that portion of SW from atmosphere that in fact
    !                enters to the ocean and participates in pentrative SW heating.
    ! nonpenSW     = non-downwelling SW flux, which is absorbed in ocean surface
    !                (in tandem w/ LW,SENS,LAT); saved only for diagnostic purposes.

    !----------------------------------------------------------------------------------------
    !BGR-June 26, 2017{
    !Temporary action to preserve answers while fixing a bug.
    ! To fix a bug in a diagnostic calculation, applyboundaryfluxesinout now returns
    !  the surface buoyancy flux. Previously, extractbuoyancyflux2d was called, meaning
    !  a second call to extractfluxes1d (causing the diagnostic net_heat to be incorrect).
    !  Note that this call to extract buoyancyflux2d was AFTER applyboundaryfluxesinout,
    !  which means it used the T/S fields after this routine.  Therefore, the surface
    !  buoyancy flux is computed here at the very end of this routine for legacy reasons.
    !  A few specific notes follow:
    !     1) The old method did not included river/calving contributions to heat flux.  This
    !        is kept consistent here via commenting code in the present extractFluxes1d <_rate>
    !        outputs, but we may reconsider this approach.
    !     2) The old method computed the buoyancy flux rate directly (by setting dt=1), instead
    !        of computing the integrated value (and dividing by dt). Hence the required
    !        additional outputs from extractFluxes1d.
    !          *** This is because: A*dt/dt =/=  A due to round off.
    !     3) The old method computed buoyancy flux after this routine, meaning the returned
    !        surface fluxes (from extractfluxes1d) must be recorded for use later in the code.
    !        We could (and maybe should) move that loop up to before the surface fluxes are
    !        applied, but this will change answers.
    !     For all these reasons we compute additional values of <_rate> which are preserved
    !     for the buoyancy flux calculation and reproduce the old answers.
    !   In the future this needs more detailed investigation to make sure everything is
    !   consistent and correct. These details shouldnt significantly effect climate,
    !   but do change answers.
    !-----------------------------------------------------------------------------------------
    if (calculate_buoyancy) then
      call extractfluxes1d(g, gv, fluxes, optics, nsw, j, dt,                        &
                  h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
                  h2d, t2d, netmassinout, netmassout, netheat, netsalt,                   &
                  pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw,                &
                  net_heat_rate=netheat_rate,net_salt_rate=netsalt_rate,                  &
                  netmassinout_rate=netmassinout_rate,pen_sw_bnd_rate=pen_sw_bnd_rate)
    else
      call extractfluxes1d(g, gv, fluxes, optics, nsw, j, dt,                        &
                  h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
                  h2d, t2d, netmassinout, netmassout, netheat, netsalt,                   &
                  pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw)
   endif
    ! ea is for passive tracers
    do i=is,ie
    !  ea(i,j,1) = netMassInOut(i)
      if (aggregate_fw_forcing) then
        netmassout(i) = netmassinout(i)
        netmassin(i) = 0.
      else
        netmassin(i) = netmassinout(i) - netmassout(i)
      endif
      if (g%mask2dT(i,j)>0.0) then
        fluxes%netMassOut(i,j) = netmassout(i)
        fluxes%netMassIn(i,j) = netmassin(i)
      else
        fluxes%netMassOut(i,j) = 0.0
        fluxes%netMassIn(i,j) = 0.0
      endif
    enddo

    ! Apply the surface boundary fluxes in three steps:
    ! A/ update mass, temp, and salinity due to all terms except mass leaving
    !    ocean (and corresponding outward heat content), and ignoring penetrative SW.
    ! 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

        ! A/ Update mass, temp, and salinity due to incoming mass flux.
        do k=1,1

          ! Change in state due to forcing
          dthickness = netmassin(i) ! Since we are adding mass, we can use all of it
          dtemp = 0.
          dsalt = 0.

          ! Update the forcing by the part to be consumed within the present k-layer.
          ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
          netmassin(i) = netmassin(i) - dthickness
          ! This line accounts for the temperature of the mass exchange
          temp_in = t2d(i,k)
          salin_in = 0.0
          dtemp = dtemp + dthickness*temp_in

          ! Diagnostics of heat content associated with mass fluxes
          if (ASSOCIATED(fluxes%heat_content_massin))                             &
            fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) +   &
                         t2d(i,k) * max(0.,dthickness) * gv%H_to_kg_m2 * fluxes%C_p * idt
          if (ASSOCIATED(fluxes%heat_content_massout))                            &
            fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
                         t2d(i,k) * min(0.,dthickness) * gv%H_to_kg_m2 * fluxes%C_p * idt
          if (ASSOCIATED(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
                         t2d(i,k) * dthickness * gv%H_to_kg_m2

          ! Determine the energetics of river mixing before updating the state.
          if (calculate_energetics .and. associated(fluxes%lrunoff) .and. cs%do_rivermix) then
            ! Here we add an additional source of TKE to the mixed layer where river
            ! is present to simulate unresolved estuaries. The TKE input is diagnosed
            ! as follows:
            !   TKE_river[m3 s-3] = 0.5*rivermix_depth*g*(1/rho)*drho_ds*
            !                       River*(Samb - Sriver) = CS%mstar*U_star^3
            ! where River is in units of 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.
            ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
            rivermixconst = -0.5*(cs%rivermix_depth*dt)*gv%m_to_H*gv%H_to_Pa

            ctke(i,j,k) = ctke(i,j,k) + max(0.0, rivermixconst*dsv_ds(i,j,1) * &
                  (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * tv%S(i,j,1))
          endif

          ! Update state
          hold     = h2d(i,k)               ! Keep original thickness in hand
          h2d(i,k) = h2d(i,k) + dthickness  ! New thickness
          if (h2d(i,k) > 0.0) then
            if (calculate_energetics .and. (dthickness > 0.)) then
              ! Calculate the energy required to mix the newly added water over
              ! the topmost grid cell.  ###CHECK THE SIGNS!!!
              ctke(i,j,k) = ctke(i,j,k) + 0.5*g_hconv2*(hold*dthickness) * &
                 ((t2d(i,k) - temp_in) * dsv_dt(i,j,k) + (tv%S(i,j,k) - salin_in) * dsv_ds(i,j,k))
            endif
            ithickness  = 1.0/h2d(i,k)      ! Inverse new thickness
            ! The "if"s below avoid changing T/S by roundoff unnecessarily
            if (dthickness /= 0. .or. dtemp /= 0.) t2d(i,k)    = (hold*t2d(i,k)    + dtemp)*ithickness
            if (dthickness /= 0. .or. dsalt /= 0.) tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness

          endif

        enddo ! k=1,1

        ! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt.
        do k=1,nz

          ! Place forcing into this layer if this layer has nontrivial thickness.
          ! For layers thin relative to 1/IforcingDepthScale, then distribute
          ! forcing into deeper layers.
          iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth*gv%m_to_H - netmassout(i) )
          ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
          fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale)

          ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we
          ! limit the forcing applied to this cell, leaving the remaining forcing to
          ! be distributed downwards.
          if (-fractionofforcing*netmassout(i) > evap_cfl_limit*h2d(i,k)) then
            fractionofforcing = -evap_cfl_limit*h2d(i,k)/netmassout(i)
          endif

          ! Change in state due to forcing

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

          ! Update the forcing by the part to be consumed within the present k-layer.
          ! If fractionOfForcing = 1, then new netMassOut vanishes.
          netmassout(i) = netmassout(i) - dthickness
          netheat(i) = netheat(i) - dtemp
          netsalt(i) = netsalt(i) - dsalt

          ! This line accounts for the temperature of the mass exchange
          dtemp = dtemp + dthickness*t2d(i,k)

          ! Diagnostics of heat content associated with mass fluxes
          if (ASSOCIATED(fluxes%heat_content_massin))                             &
            fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) +   &
                         tv%T(i,j,k) * max(0.,dthickness) * gv%H_to_kg_m2 * fluxes%C_p * idt
          if (ASSOCIATED(fluxes%heat_content_massout))                            &
            fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
                         tv%T(i,j,k) * min(0.,dthickness) * gv%H_to_kg_m2 * fluxes%C_p * idt
          if (ASSOCIATED(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
                         tv%T(i,j,k) * dthickness * gv%H_to_kg_m2
!NOTE tv%T should be T2d

          ! Update state by the appropriate increment.
          hold     = h2d(i,k)               ! Keep original thickness in hand
          h2d(i,k) = h2d(i,k) + dthickness  ! New thickness

          if (h2d(i,k) > 0.) then
            if (calculate_energetics) then
              ! Calculate the energy required to mix the newly added water over
              ! the topmost grid cell, assuming that the fluxes of heat and salt
              ! and rejected brine are initially applied in vanishingly thin
              ! layers at the top of the layer before being mixed throughout
              ! the layer.  Note that dThickness is always <= 0. ###CHECK THE SIGNS!!!
              ctke(i,j,k) = ctke(i,j,k) - (0.5*h2d(i,k)*g_hconv2) * &
                 ((dtemp - dthickness*t2d(i,k)) * dsv_dt(i,j,k) + &
                  (dsalt - dthickness*tv%S(i,j,k)) * dsv_ds(i,j,k))
            endif
            ithickness  = 1.0/h2d(i,k) ! Inverse of new thickness
            t2d(i,k)    = (hold*t2d(i,k) + dtemp)*ithickness
            tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness
          elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling
            call forcing_singlepointprint(fluxes,g,i,j,'applyBoundaryFluxesInOut (h<0)')
            write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
            write(0,*) 'applyBoundaryFluxesInOut(): netT,netS,netH=',netheat(i),netsalt(i),netmassinout(i)
            write(0,*) 'applyBoundaryFluxesInOut(): dT,dS,dH=',dtemp,dsalt,dthickness
            write(0,*) 'applyBoundaryFluxesInOut(): h(n),h(n+1),k=',hold,h2d(i,k),k
            call mom_error(fatal, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
                           "Complete mass loss in column!")
          endif

        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

        if (.not. cs%ignore_fluxes_over_land) then
           call forcing_singlepointprint(fluxes,g,i,j,'applyBoundaryFluxesInOut (land)')
           write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
           write(0,*) 'applyBoundaryFluxesInOut(): netHeat,netSalt,netMassIn,netMassOut=',&
                   netheat(i),netsalt(i),netmassin(i),netmassout(i)

           call mom_error(fatal, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
                                 "Mass loss over land?")
        endif

      endif

      ! If anything remains after the k-loop, then we have grounded out, which is a problem.
      if (netmassin(i)+netmassout(i) /= 0.0) then
!$OMP critical
        numberofgroundings = numberofgroundings +1
        if (numberofgroundings<=maxgroundings) then
          iground(numberofgroundings) = i ! Record i,j location of event for
          jground(numberofgroundings) = j ! warning message
          hgrounding(numberofgroundings) = netmassin(i)+netmassout(i)
        endif
!$OMP end critical
        if (cs%id_createdH>0) cs%createdH(i,j) = cs%createdH(i,j) - (netmassin(i)+netmassout(i))/dt
      endif

    enddo ! i

    ! Step C/ in the application of fluxes
    ! Heat by the convergence of penetrating SW.
    ! SW penetrative heating uses the updated thickness from above.

    ! Save temperature before increment with SW heating
    ! and initialize CS%penSWflux_diag to zero.
    if(cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0) then
      do k=1,nz ; do i=is,ie
        cs%penSW_diag(i,j,k)     = t2d(i,k)
        cs%penSWflux_diag(i,j,k) = 0.0
      enddo ; enddo
      k=nz+1 ; do i=is,ie
         cs%penSWflux_diag(i,j,k) = 0.0
      enddo
    endif

    if (calculate_energetics) then
      call absorbremainingsw(g, gv, h2d, opacityband, nsw, j, dt, h_limit_fluxes, &
                             .false., .true., t2d, pen_sw_bnd, tke=pen_tke_2d, dsv_dt=dsv_dt_2d)
      k = 1 ! For setting break-points.
      do k=1,nz ; do i=is,ie
        ctke(i,j,k) = ctke(i,j,k) + pen_tke_2d(i,k)
      enddo ; enddo
    else
      call absorbremainingsw(g, gv, h2d, opacityband, nsw, j, dt, h_limit_fluxes, &
                             .false., .true., t2d, pen_sw_bnd)
    endif


    ! Step D/ copy updated thickness and temperature
    ! 2d slice now back into model state.
    do k=1,nz ; do i=is,ie
      h(i,j,k)    = h2d(i,k)
      tv%T(i,j,k) = t2d(i,k)
    enddo ; enddo

    ! Diagnose heating (W/m2) applied to a grid cell from SW penetration
    ! Also diagnose the penetrative SW heat flux at base of layer.
    if(cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0) then

      ! convergence of SW into a layer
      do k=1,nz ; do i=is,ie
        cs%penSW_diag(i,j,k) = (t2d(i,k)-cs%penSW_diag(i,j,k))*h(i,j,k) * idt * tv%C_p * gv%H_to_kg_m2
      enddo ; enddo

      ! Perform a cumulative sum upwards from bottom to
      ! diagnose penetrative SW flux at base of tracer cell.
      ! CS%penSWflux_diag(i,j,k=1)    is penetrative shortwave at top of ocean.
      ! CS%penSWflux_diag(i,j,k=kbot+1) is zero, since assume no SW penetrates rock.
      ! CS%penSWflux_diag = rsdo  and CS%penSW_diag = rsdoabsorb
      ! rsdoabsorb(k) = rsdo(k) - rsdo(k+1), so that rsdo(k) = rsdo(k+1) + rsdoabsorb(k)
      if(cs%id_penSWflux_diag > 0) then
        do k=nz,1,-1 ; do i=is,ie
          cs%penSWflux_diag(i,j,k) = cs%penSW_diag(i,j,k) + cs%penSWflux_diag(i,j,k+1)
        enddo ; enddo
      endif

    endif

    ! Fill CS%nonpenSW_diag
    if(cs%id_nonpenSW_diag > 0) then
      do i=is,ie
        cs%nonpenSW_diag(i,j) = nonpensw(i)
      enddo
    endif

    ! BGR: Get buoyancy flux to return for ePBL
    !  We want the rate, so we use the rate values returned from extractfluxes1d.
    !  Note that the *dt values could be divided by dt here, but
    !  1) Answers will change due to round-off
    !  2) Be sure to save their values BEFORE fluxes are used.
    if (calculate_buoyancy) then
       drhodt(:) = 0.0
       drhods(:) = 0.0
       netpen(:,:) = 0.0
       ! Sum over bands and attenuate as a function of depth
       ! netPen is the netSW as a function of depth
       call sumswoverbands(g, gv, h2d(:,:), optics%opacity_band(:,:,j,:), nsw, j, dt, &
            h_limit_fluxes, .true., pen_sw_bnd_rate, netpen)
       ! Density derivatives
       call calculate_density_derivs(t2d(:,1), tv%S(:,j,1), surfpressure, &
            drhodt, drhods, start, npts, tv%eqn_of_state)
       ! 1. Adjust netSalt to reflect dilution effect of FW flux
       ! 2. Add in the SW heating for purposes of calculating the net
       ! surface buoyancy flux affecting the top layer.
       ! 3. Convert to a buoyancy flux, excluding penetrating SW heating
       !    BGR-Jul 5, 2017: The contribution of SW heating here needs investigated for ePBL.
       skinbuoyflux(g%isc:g%iec,j) = - gorho * ( drhods(g%isc:g%iec) * (netsalt_rate(g%isc:g%iec) &
            - tv%S(g%isc:g%iec,j,1) * netmassinout_rate(g%isc:g%iec)* gv%H_to_m )&
            + drhodt(g%isc:g%iec) * ( netheat_rate(g%isc:g%iec) +        &
            netpen(g%isc:g%iec,1))) * gv%H_to_m ! m^2/s^3
    endif

  enddo ! j-loop finish

  ! Post the diagnostics
  if (cs%id_createdH       > 0) call post_data(cs%id_createdH      , cs%createdH      , cs%diag)
  if (cs%id_penSW_diag     > 0) call post_data(cs%id_penSW_diag    , cs%penSW_diag    , cs%diag)
  if (cs%id_penSWflux_diag > 0) call post_data(cs%id_penSWflux_diag, cs%penSWflux_diag, cs%diag)
  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
    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) ), &
                             g%geoLatT( iground(i), jground(i)) , hgrounding(i)
      call mom_error(warning, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
                              "Mass created. x,y,dh= "//trim(mesg), all_print=.true.)
    enddo

    if (numberofgroundings - maxgroundings > 0) then
      write(mesg, '(i4)') numberofgroundings - maxgroundings
      call mom_error(warning, "MOM_diabatic_driver:F90, applyBoundaryFluxesInOut(): "//&
                              trim(mesg) // " groundings remaining")
    endif
  endif

Here is the call graph for this function:

Here is the caller graph for this function:

diabatic_aux_end()

subroutine, public mom_diabatic_aux::diabatic_aux_end

(

type(diabatic_aux_cs), pointer

CS

)

Definition at line 1459 of file MOM_diabatic_aux.F90.

  type(diabatic_aux_cs), pointer :: cs

  if (.not.associated(cs)) return

  if (cs%id_createdH       >0) deallocate(cs%createdH)
  if (cs%id_penSW_diag     >0) deallocate(cs%penSW_diag)
  if (cs%id_penSWflux_diag >0) deallocate(cs%penSWflux_diag)
  if (cs%id_nonpenSW_diag  >0) deallocate(cs%nonpenSW_diag)

  if (associated(cs)) deallocate(cs)

diabatic_aux_init()

subroutine, public mom_diabatic_aux::diabatic_aux_init	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(diabatic_aux_cs), pointer	CS,
		logical, intent(in)	useALEalgorithm,
		logical, intent(in)	use_ePBL
	)

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	param_file	A structure to parse for run-time parameters

Definition at line 1331 of file MOM_diabatic_aux.F90.

References id_clock_frazil, and id_clock_uv_at_h.

  type(time_type),         intent(in)    :: time
  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(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(diag_ctrl), target, intent(inout) :: diag
  type(diabatic_aux_cs),   pointer       :: cs
  logical,                 intent(in)    :: usealealgorithm
  logical,                 intent(in)    :: use_epbl

! Arguments:
!  (in)     Time       = current model time
!  (in)     G          = ocean grid structure
!  (in)     GV - The ocean's vertical grid structure.
!  (in)     param_file = structure indicating the open file to parse for parameter values
!  (in)     diag       = structure used to regulate diagnostic output
!  (in/out) CS         = pointer set to point to the control structure for this module
!  (in)     use_ePBL   = If true, use the implicit energetics planetary boundary
!                        layer scheme to determine the diffusivity in the
!                        surface boundary layer.
  type(vardesc) :: vd

! This "include" declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mod  = "MOM_diabatic_aux" ! This module's name.
  character(len=48)  :: thickness_units
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nbands
  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_aux_init called with an "// &
                            "associated control structure.")
    return
  else
    allocate(cs)
   endif

  cs%diag => diag

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

  call get_param(param_file, mod, "RECLAIM_FRAZIL", cs%reclaim_frazil, &
                 "If true, try to use any frazil heat deficit to cool any\n"//&
                 "overlying layers down to the freezing point, thereby \n"//&
                 "avoiding the creation of thin ice when the SST is above \n"//&
                 "the freezing point.", default=.true.)
  call get_param(param_file, mod, "PRESSURE_DEPENDENT_FRAZIL", &
                                cs%pressure_dependent_frazil, &
                 "If true, use a pressure dependent freezing temperature \n"//&
                 "when making frazil. The default is false, which will be \n"//&
                 "faster but is inappropriate with ice-shelf cavities.", &
                 default=.false.)

  if (use_epbl) then
    call get_param(param_file, mod, "IGNORE_FLUXES_OVER_LAND", cs%ignore_fluxes_over_land,&
         "If true, the model does not check if fluxes are being applied\n"//&
         "over land points. This is needed when the ocean is coupled \n"//&
         "with ice shelves and sea ice, since the sea ice mask needs to \n"//&
         "be different than the ocean mask to avoid sea ice formation \n"//&
         "under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
    call get_param(param_file, mod, "DO_RIVERMIX", cs%do_rivermix, &
                 "If true, apply additional mixing whereever there is \n"//&
                 "runoff, so that it is mixed down to RIVERMIX_DEPTH \n"//&
                 "if the ocean is that deep.", default=.false.)
    if (cs%do_rivermix) &
      call get_param(param_file, mod, "RIVERMIX_DEPTH", cs%rivermix_depth, &
                 "The depth to which rivers are mixed if DO_RIVERMIX is \n"//&
                 "defined.", units="m", default=0.0)
  else ; cs%do_rivermix = .false. ; cs%rivermix_depth = 0.0 ; endif
  if (gv%nkml == 0) then
    call get_param(param_file, mod, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
                   "If true, use the fluxes%runoff_Hflx field to set the \n"//&
                   "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
                   default=.false.)
    call get_param(param_file, mod, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
                   "If true, use the fluxes%calving_Hflx field to set the \n"//&
                   "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
                   default=.false.)
  else
    cs%use_river_heat_content = .false.
    cs%use_calving_heat_content = .false.
  endif

  if (usealealgorithm) then
    cs%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, &
        time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", &
        "meter second-1")
    if (cs%id_createdH>0) allocate(cs%createdH(isd:ied,jsd:jed))

    ! diagnostic for heating of a grid cell from convergence of SW heat into the cell
    cs%id_penSW_diag = register_diag_field('ocean_model', 'rsdoabsorb',                     &
          diag%axesTL, time, 'Convergence of Penetrative Shortwave Flux in Sea Water Layer',&
          'Watt meter-2', standard_name='net_rate_of_absorption_of_shortwave_energy_in_ocean_layer',v_extensive=.true.)

    ! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces)
    ! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock).
    cs%id_penSWflux_diag = register_diag_field('ocean_model', 'rsdo',                               &
          diag%axesTi, time, 'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',&
          'Watt meter-2', standard_name='downwelling_shortwave_flux_in_sea_water')

    ! need both arrays for the SW diagnostics (one for flux, one for convergence)
    if (cs%id_penSW_diag>0 .or. cs%id_penSWflux_diag>0) then
       allocate(cs%penSW_diag(isd:ied,jsd:jed,nz))
       cs%penSW_diag(:,:,:) = 0.0
       allocate(cs%penSWflux_diag(isd:ied,jsd:jed,nz+1))
       cs%penSWflux_diag(:,:,:) = 0.0
    endif

    ! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface)
    cs%id_nonpenSW_diag = register_diag_field('ocean_model', 'nonpenSW',                       &
          diag%axesT1, time,                                                                   &
          'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',&
          'Watt meter-2', standard_name='nondownwelling_shortwave_flux_in_sea_water')
    if (cs%id_nonpenSW_diag > 0) then
       allocate(cs%nonpenSW_diag(isd:ied,jsd:jed))
       cs%nonpenSW_diag(:,:) = 0.0
    endif
  endif

  id_clock_uv_at_h = cpu_clock_id('(Ocean find_uv_at_h)', grain=clock_routine)
  id_clock_frazil  = cpu_clock_id('(Ocean frazil)', grain=clock_routine)

diagnosemldbydensitydifference()

subroutine, public mom_diabatic_aux::diagnosemldbydensitydifference	(	integer, intent(in)	id_MLD,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		type(thermo_var_ptrs), intent(in)	tv,
		real, intent(in)	densityDiff,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(diag_ctrl), pointer	diagPtr,
		integer, intent(in), optional	id_N2subML,
		integer, intent(in), optional	id_MLDsq
	)

Diagnose a mixed layer depth (MLD) determined by a given density difference with the surface. This routine is appropriate in MOM_diabatic_driver due to its position within the time stepping.

Parameters

[in]	g	Grid type
[in]	gv	ocean vertical grid structure
[in]	id_mld	Handle (ID) of MLD diagnostic
[in]	h	Layer thickness
[in]	tv	Thermodynamics type
[in]	densitydiff	Density difference to determine MLD (kg/m3)
	diagptr	Diagnostics structure
[in]	id_n2subml	Optional handle (ID) of subML stratification
[in]	id_mldsq	Optional handle (ID) of squared MLD

Definition at line 706 of file MOM_diabatic_aux.F90.

Referenced by mom_diabatic_driver::diabatic().

  type(ocean_grid_type),                 intent(in) :: g           !< Grid type
  type(verticalgrid_type),               intent(in) :: gv          !< ocean vertical grid structure
  integer,                               intent(in) :: id_mld      !< Handle (ID) of MLD diagnostic
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h        !< Layer thickness
  type(thermo_var_ptrs),                 intent(in) :: tv          !< Thermodynamics type
  real,                                  intent(in) :: densitydiff !< Density difference to determine MLD (kg/m3)
  type(diag_ctrl),                       pointer    :: diagptr     !< Diagnostics structure
  integer,                     optional, intent(in) :: id_n2subml  !< Optional handle (ID) of subML stratification
  integer,                     optional, intent(in) :: id_mldsq    !< Optional handle (ID) of squared MLD

  ! Local variables
  real, dimension(SZI_(G))          :: rhosurf, deltarhoatkm1, deltarhoatk, dk, dkm1, pref_mld
  real, dimension(SZI_(G))          :: rhoatk, rho1, d1, pref_n2 ! Used for N2
  real, dimension(SZI_(G), SZJ_(G)) :: mld ! Diagnosed mixed layer depth
  real, dimension(SZI_(G), SZJ_(G)) :: submln2 ! Diagnosed stratification below ML
  real, dimension(SZI_(G), SZJ_(G)) :: mld2 ! Diagnosed MLD^2
  real, parameter                   :: dz_subml = 50. ! Depth below ML over which to diagnose stratification (m)
  integer :: i, j, is, ie, js, je, k, nz, id_n2, id_sq
  real :: afac, ddrho

  id_n2 = -1
  if (PRESENT(id_n2subml)) id_n2 = id_n2subml

  id_sq = -1
  if (PRESENT(id_n2subml)) id_sq = id_mldsq

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  pref_mld(:) = 0. ; pref_n2(:) = 0.
  do j=js,je
    do i=is,ie ; dk(i) = 0.5 * h(i,j,1) * gv%H_to_m ; enddo ! Depth of center of surface layer
    call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pref_mld, rhosurf, is, ie-is+1, tv%eqn_of_state)
    do i=is,ie
      deltarhoatk(i) = 0.
      mld(i,j) = 0.
      if (id_n2>0) then
        submln2(i,j) = 0.
        rho1(i) = 0.
        d1(i) = 0.
        pref_n2(i) = gv%g_Earth * gv%Rho0 * h(i,j,1) * gv%H_to_m ! Boussinesq approximation!!!! ?????
        !### This should be: pRef_N2(i) = GV%g_Earth * GV%H_to_kg_m2 * h(i,j,1) ! This might change answers at roundoff.
      endif
    enddo
    do k=2,nz
      do i=is,ie
        dkm1(i) = dk(i) ! Depth of center of layer K-1
        dk(i) = dk(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * gv%H_to_m ! Depth of center of layer K
      enddo

      ! Stratification, N2, immediately below the mixed layer, averaged over at least 50 m.
      if (id_n2>0) then
        do i=is,ie
          pref_n2(i) = pref_n2(i) + gv%g_Earth * gv%Rho0 * h(i,j,k) * gv%H_to_m ! Boussinesq approximation!!!! ?????
          !### This should be: pRef_N2(i) = pRev_N2(i) + GV%g_Earth * GV%H_to_kg_m2 * h(i,j,k) ! This might change answers at roundoff.
        enddo
        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_n2, rhoatk, is, ie-is+1, tv%eqn_of_state)
        do i=is,ie
          if (mld(i,j)>0. .and. submln2(i,j)==0.) then ! This block is below the mixed layer
            if (d1(i)==0.) then ! Record the density, depth and pressure, immediately below the ML
              rho1(i) = rhoatk(i)
              d1(i) = dk(i)
              !### It looks to me like there is bad logic here. - RWH
              ! Use pressure at the bottom of the upper layer used in calculating d/dz rho
              pref_n2(i) = pref_n2(i) + gv%g_Earth * gv%Rho0 * h(i,j,k) * gv%H_to_m ! Boussinesq approximation!!!! ?????
              !### This line should be: pRef_N2(i) = pRev_N2(i) + GV%g_Earth * GV%H_to_kg_m2 * h(i,j,k) ! This might change answers at roundoff.
            endif
            if (d1(i)>0. .and. dk(i)-d1(i)>=dz_subml) then
              submln2(i,j) = gv%g_Earth/ gv%Rho0 * (rho1(i)-rhoatk(i)) / (d1(i) - dk(i))
            endif
          endif
        enddo ! i-loop
      endif ! id_N2>0

      ! Mixed-layer depth, using sigma-0 (surface reference pressure)
      do i=is,ie ; deltarhoatkm1(i) = deltarhoatk(i) ; enddo ! Store value from previous iteration of K
      call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_mld, deltarhoatk, is, ie-is+1, tv%eqn_of_state)
      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
          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 (id_N2>0 .and. subMLN2(i,j)==0. .and. d1(i)>0. .and. dK(i)-d1(i)>0.) then
   !    ! Use what ever stratification we can, measured over what ever distance is available
   !    subMLN2(i,j) = GV%g_Earth/ GV%Rho0 * (rho1(i)-rhoAtK(i)) / (d1(i) - dK(i))
   !  endif
    enddo
  enddo ! j-loop

  if (id_mld > 0) call post_data(id_mld, mld, diagptr)
  if (id_n2 > 0)  call post_data(id_n2, submln2 , diagptr)
  if (id_sq > 0)  call post_data(id_sq, mld2, diagptr)

Here is the caller graph for this function:

differential_diffuse_t_s()

subroutine, public mom_diabatic_aux::differential_diffuse_t_s	(	real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		type(thermo_var_ptrs), intent(inout)	tv,
		type(vertvisc_type), intent(in)	visc,
		real, intent(in)	dt,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV
	)

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	h	Layer thicknesses, in H (usually m or kg m-2)

Definition at line 259 of file MOM_diabatic_aux.F90.

  type(ocean_grid_type),                 intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),               intent(in)    :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                 intent(inout) :: tv
  type(vertvisc_type),                   intent(in)    :: visc
  real,                                  intent(in)    :: dt

! This subroutine applies double diffusion to T & S, assuming no diapycal mass
! fluxes, using a simple triadiagonal solver.

! Arguments: h - Layer thickness, in m or kg m-2.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      visc - A structure containing vertical viscosities, bottom boundary
!                   layer properies, and related fields.
!  (in)      dt - Time increment, in s.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.

  real, dimension(SZI_(G)) :: &
    b1_t, b1_s, &  !  Variables used by the tridiagonal solvers of T & S, in H.
    d1_t, d1_s     !  Variables used by the tridiagonal solvers, nondim.
  real, dimension(SZI_(G),SZK_(G)) :: &
    c1_t, c1_s     !  Variables used by the tridiagonal solvers, in m or kg m-2.
  real, dimension(SZI_(G),SZK_(G)+1) :: &
    mix_t, mix_s   !  Mixing distances in both directions across each
                   !  interface, in m or kg m-2.
  real :: h_tr         ! h_tr is h at tracer points with a tiny thickness
                       ! added to ensure positive definiteness, in m or kg m-2.
  real :: h_neglect    ! A thickness that is so small it is usually lost
                       ! in roundoff and can be neglected, in m or kg m-2.
  real :: i_h_int      ! The inverse of the thickness associated with an
                       ! interface, in m-1 or m2 kg-1.
  real :: b_denom_t    ! The first term in the denominators for the expressions
  real :: b_denom_s    ! for b1_T and b1_S, both in m or kg m-2.

  integer :: i, j, k, is, ie, js, je, nz
  real, pointer :: t(:,:,:), s(:,:,:), kd_t(:,:,:), kd_s(:,:,:)
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  h_neglect = gv%H_subroundoff

  if (.not.associated(tv%T)) call mom_error(fatal, &
      "differential_diffuse_T_S: Called with an unassociated tv%T")
  if (.not.associated(tv%S)) call mom_error(fatal, &
      "differential_diffuse_T_S: Called with an unassociated tv%S")
  if (.not.associated(visc%Kd_extra_T)) call mom_error(fatal, &
      "differential_diffuse_T_S: Called with an unassociated visc%Kd_extra_T")
  if (.not.associated(visc%Kd_extra_S)) call mom_error(fatal, &
      "differential_diffuse_T_S: Called with an unassociated visc%Kd_extra_S")

  t => tv%T ; s => tv%S
  kd_t => visc%Kd_extra_T ; kd_s => visc%Kd_extra_S
!$OMP parallel do default(none) shared(is,ie,js,je,h,h_neglect,dt,Kd_T,Kd_S,G,GV,T,S,nz) &
!$OMP                          private(I_h_int,mix_T,mix_S,h_tr,b1_T,b1_S, &
!$OMP                                  d1_T,d1_S,c1_T,c1_S,b_denom_T,b_denom_S)
  do j=js,je
    do i=is,ie
      i_h_int = 1.0 / (0.5 * (h(i,j,1) + h(i,j,2)) + h_neglect)
      mix_t(i,2) = ((dt * kd_t(i,j,2)) * gv%m_to_H**2) * i_h_int
      mix_s(i,2) = ((dt * kd_s(i,j,2)) * gv%m_to_H**2) * i_h_int

      h_tr = h(i,j,1) + h_neglect
      b1_t(i) = 1.0 / (h_tr + mix_t(i,2))
      b1_s(i) = 1.0 / (h_tr + mix_s(i,2))
      d1_t(i) = h_tr * b1_t(i)
      d1_s(i) = h_tr * b1_s(i)
      t(i,j,1) = (b1_t(i)*h_tr)*t(i,j,1)
      s(i,j,1) = (b1_s(i)*h_tr)*s(i,j,1)
    enddo
    do k=2,nz-1 ; do i=is,ie
      ! Calculate the mixing across the interface below this layer.
      i_h_int = 1.0 / (0.5 * (h(i,j,k) + h(i,j,k+1)) + h_neglect)
      mix_t(i,k+1) = ((dt * kd_t(i,j,k+1)) * gv%m_to_H**2) * i_h_int
      mix_s(i,k+1) = ((dt * kd_s(i,j,k+1)) * gv%m_to_H**2) * i_h_int

      c1_t(i,k) = mix_t(i,k) * b1_t(i)
      c1_s(i,k) = mix_s(i,k) * b1_s(i)

      h_tr = h(i,j,k) + h_neglect
      b_denom_t = h_tr + d1_t(i)*mix_t(i,k)
      b_denom_s = h_tr + d1_s(i)*mix_s(i,k)
      b1_t(i) = 1.0 / (b_denom_t + mix_t(i,k+1))
      b1_s(i) = 1.0 / (b_denom_s + mix_s(i,k+1))
      d1_t(i) = b_denom_t * b1_t(i)
      d1_s(i) = b_denom_s * b1_s(i)

      t(i,j,k) = b1_t(i) * (h_tr*t(i,j,k) + mix_t(i,k)*t(i,j,k-1))
      s(i,j,k) = b1_s(i) * (h_tr*s(i,j,k) + mix_s(i,k)*s(i,j,k-1))
    enddo ; enddo
    do i=is,ie
      c1_t(i,nz) = mix_t(i,nz) * b1_t(i)
      c1_s(i,nz) = mix_s(i,nz) * b1_s(i)

      h_tr = h(i,j,nz) + h_neglect
      b1_t(i) = 1.0 / (h_tr + d1_t(i)*mix_t(i,nz))
      b1_s(i) = 1.0 / (h_tr + d1_s(i)*mix_s(i,nz))

      t(i,j,nz) = b1_t(i) * (h_tr*t(i,j,nz) + mix_t(i,nz)*t(i,j,nz-1))
      s(i,j,nz) = b1_s(i) * (h_tr*s(i,j,nz) + mix_s(i,nz)*s(i,j,nz-1))
    enddo
    do k=nz-1,1,-1 ; do i=is,ie
      t(i,j,k) = t(i,j,k) + c1_t(i,k+1)*t(i,j,k+1)
      s(i,j,k) = s(i,j,k) + c1_s(i,k+1)*s(i,j,k+1)
    enddo ; enddo
  enddo

find_uv_at_h()

subroutine, public mom_diabatic_aux::find_uv_at_h	(	real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)	u,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)	v,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)	u_h,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)	v_h,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(in), optional	ea,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(in), optional	eb
	)

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	u	The zonal velocity, in m s-1
[in]	v	The meridional velocity, in m s-1
[in]	h	Layer thicknesses, in H (usually m or kg m-2)

Definition at line 600 of file MOM_diabatic_aux.F90.

References id_clock_uv_at_h.

Referenced by mom_diabatic_driver::diabatic().

  type(ocean_grid_type),                     intent(in)  :: g    !< The ocean's grid structure
  type(verticalgrid_type),                   intent(in)  :: gv   !< The ocean's vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)  :: u    !< The zonal velocity, in m s-1
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)  :: v    !< The meridional velocity, in m s-1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)  :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(out) :: u_h, v_h
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in), optional  :: ea, eb
!   This subroutine calculates u_h and v_h (velocities at thickness
! points), optionally using the entrainments (in m) passed in as arguments.

! Arguments: u - Zonal velocity, in m s-1.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m or kg m-2.
!  (out)     u_h - The zonal velocity at thickness points after
!                  entrainment, in m s-1.
!  (out)     v_h - The meridional velocity at thickness points after
!                  entrainment, in m s-1.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in, opt) ea - The amount of fluid entrained from the layer above within
!                 this time step, in units of m or kg m-2.  Omitting ea is the
!                 same as setting it to 0.
!  (in, opt) eb - The amount of fluid entrained from the layer below within
!                 this time step, in units of m or kg m-2.  Omitting eb is the
!                 same as setting it to 0.  ea and eb must either be both
!                 present or both absent.

  real :: b_denom_1    ! The first term in the denominator of b1 in m or kg m-2.
  real :: h_neglect    ! A thickness that is so small it is usually lost
                       ! in roundoff and can be neglected, in m or kg m-2.
  real :: b1(szi_(g)), d1(szi_(g)), c1(szi_(g),szk_(g))
  real :: a_n(szi_(g)), a_s(szi_(g))  ! Fractional weights of the neighboring
  real :: a_e(szi_(g)), a_w(szi_(g))  ! velocity points, ~1/2 in the open
                                      ! ocean, nondimensional.
  real :: s, idenom
  logical :: mix_vertically
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  call cpu_clock_begin(id_clock_uv_at_h)
  h_neglect = gv%H_subroundoff

  mix_vertically = present(ea)
  if (present(ea) .neqv. present(eb)) call mom_error(fatal, &
      "find_uv_at_h: Either both ea and eb or neither one must be present "// &
      "in call to find_uv_at_h.")
!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,h,h_neglect, &
!$OMP                                  eb,u_h,u,v_h,v,nz,ea)                     &
!$OMP                          private(s,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1)
  do j=js,je
    do i=is,ie
      s = g%areaCu(i-1,j)+g%areaCu(i,j)
      if (s>0.0) then
        idenom = sqrt(0.5*g%IareaT(i,j)/s)
        a_w(i) = g%areaCu(i-1,j)*idenom
        a_e(i) = g%areaCu(i,j)*idenom
      else
        a_w(i) = 0.0 ; a_e(i) = 0.0
      endif

      s = g%areaCv(i,j-1)+g%areaCv(i,j)
      if (s>0.0) then
        idenom = sqrt(0.5*g%IareaT(i,j)/s)
        a_s(i) = g%areaCv(i,j-1)*idenom
        a_n(i) = g%areaCv(i,j)*idenom
      else
        a_s(i) = 0.0 ; a_n(i) = 0.0
      endif
    enddo

    if (mix_vertically) then
      do i=is,ie
        b_denom_1 = h(i,j,1) + h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
        d1(i) = b_denom_1 * b1(i)
        u_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_e(i)*u(i,j,1) + a_w(i)*u(i-1,j,1))
        v_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_n(i)*v(i,j,1) + a_s(i)*v(i,j-1,1))
      enddo
      do k=2,nz ; do i=is,ie
        c1(i,k) = eb(i,j,k-1) * b1(i)
        b_denom_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
        d1(i) = b_denom_1 * b1(i)
        u_h(i,j,k) = (h(i,j,k) * (a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)) + &
                      ea(i,j,k)*u_h(i,j,k-1))*b1(i)
        v_h(i,j,k) = (h(i,j,k) * (a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)) + &
                      ea(i,j,k)*v_h(i,j,k-1))*b1(i)
      enddo ; enddo
      do k=nz-1,1,-1 ; do i=is,ie
        u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1)
        v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1)
      enddo ; enddo
    else
      do k=1,nz ; do i=is,ie
        u_h(i,j,k) = a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)
        v_h(i,j,k) = a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)
      enddo ; enddo
    endif
  enddo

  call cpu_clock_end(id_clock_uv_at_h)

Here is the caller graph for this function:

insert_brine()

subroutine, public mom_diabatic_aux::insert_brine	(	real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		type(thermo_var_ptrs), intent(inout)	tv,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(forcing), intent(in)	fluxes,
		integer, intent(in)	nkmb,
		type(diabatic_aux_cs), intent(in)	CS,
		real, intent(in)	dt,
		integer, intent(in)	id_brine_lay
	)

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	h	Layer thicknesses, in H (usually m or kg m-2)

Definition at line 430 of file MOM_diabatic_aux.F90.

  type(ocean_grid_type),                 intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),               intent(in)    :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                 intent(inout) :: tv
  type(forcing),                         intent(in)    :: fluxes
  integer,                               intent(in)    :: nkmb
  type(diabatic_aux_cs),                 intent(in)    :: cs
  real,                                  intent(in)    :: dt
  integer,                               intent(in)    :: id_brine_lay

! Insert salt from brine rejection into the first layer below
! the mixed layer which both contains mass and in which the
! change in layer density remains stable after the addition
! of salt via brine rejection.

! Arguments: h - Layer thickness, in m.
!  (in/out)  tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      fluxes = A structure containing pointers to any possible
!                     forcing fields; unused fields have NULL ptrs.
!  (in)      nkmb - The number of layers in the mixed and buffer layers.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 diabatic_driver_init.

  real :: salt(szi_(g)) ! The amount of salt rejected from
                        !  sea ice. [grams]
  real :: dzbr(szi_(g)) ! cumulative depth over which brine is distributed
  real :: inject_layer(szi_(g),szj_(g)) ! diagnostic

  real :: p_ref_cv(szi_(g))
  real :: t(szi_(g),szk_(g))
  real :: s(szi_(g),szk_(g))
  real :: h_2d(szi_(g),szk_(g))
  real :: rcv(szi_(g),szk_(g))
  real :: mc  ! A layer's mass in kg m-2 .
  real :: s_new,r_new,t0,scale, cdz
  integer :: i, j, k, is, ie, js, je, nz, ks

  real, parameter :: brine_dz = 1.0  ! minumum thickness over which to distribute brine
  real, parameter :: s_max = 45.0    ! salinity bound

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  if (.not.ASSOCIATED(fluxes%salt_flux)) return

  p_ref_cv(:)  = tv%P_ref

  inject_layer = nz

  do j=js,je

    salt(:)=0.0 ; dzbr(:)=0.0

    do i=is,ie ; if (g%mask2dT(i,j) > 0.) then
      salt(i) = dt * (1000. * fluxes%salt_flux(i,j))
    endif ; enddo

    do k=1,nz
      do i=is,ie
        t(i,k)=tv%T(i,j,k); s(i,k)=tv%S(i,j,k)
        ! avoid very small thickness
        h_2d(i,k)=max(h(i,j,k), gv%Angstrom)
      enddo

      call calculate_density(t(:,k), s(:,k), p_ref_cv, rcv(:,k), is, &
                             ie-is+1, tv%eqn_of_state)
    enddo

    ! First, try to find an interior layer where inserting all the salt
    ! will not cause the layer to become statically unstable.
    ! Bias towards deeper layers.

    do k=nkmb+1,nz-1 ; do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
        mc = gv%H_to_kg_m2 * h_2d(i,k)
        s_new = s(i,k) + salt(i)/mc
        t0 = t(i,k)
        call calculate_density(t0,s_new,tv%P_Ref,r_new,tv%eqn_of_state)
        if (r_new < 0.5*(rcv(i,k)+rcv(i,k+1)) .and. s_new<s_max) then
          dzbr(i)=dzbr(i)+h_2d(i,k)
          inject_layer(i,j) = min(inject_layer(i,j),real(k))
        endif
      endif
    enddo ; enddo

    ! Then try to insert into buffer layers if they exist
    do k=nkmb,gv%nkml+1,-1 ; do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
        mc = gv%H_to_kg_m2 * h_2d(i,k)
        dzbr(i)=dzbr(i)+h_2d(i,k)
        inject_layer(i,j) = min(inject_layer(i,j),real(k))
      endif
    enddo ; enddo

    ! finally if unable to find a layer to insert, then place in mixed layer

    do k=1,gv%nkml ; do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
        mc = gv%H_to_kg_m2 * h_2d(i,k)
        dzbr(i)=dzbr(i)+h_2d(i,k)
        inject_layer(i,j) = min(inject_layer(i,j),real(k))
      endif
    enddo ; enddo


    do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .and. salt(i) > 0.) then
    !   if (dzbr(i)< brine_dz) call MOM_error(FATAL,"insert_brine: failed")
        ks=inject_layer(i,j)
        cdz=0.0
        do k=ks,nz
          mc = gv%H_to_kg_m2 * h_2d(i,k)
          scale = h_2d(i,k)/dzbr(i)
          cdz=cdz+h_2d(i,k)
          if (cdz > 1.0) exit
          tv%S(i,j,k) = tv%S(i,j,k) + scale*salt(i)/mc
        enddo
      endif
    enddo

  enddo

  if (cs%id_brine_lay > 0) call post_data(cs%id_brine_lay,inject_layer,cs%diag)

make_frazil()

subroutine, public mom_diabatic_aux::make_frazil	(	real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		type(thermo_var_ptrs), intent(inout)	tv,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(diabatic_aux_cs), intent(in)	CS,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional	p_surf
	)

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	h	Layer thicknesses, in H (usually m or kg m-2)

Definition at line 138 of file MOM_diabatic_aux.F90.

References id_clock_frazil.

  type(ocean_grid_type),                 intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),               intent(in)    :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                 intent(inout) :: tv
  type(diabatic_aux_cs),                 intent(in)    :: cs
  real, dimension(SZI_(G),SZJ_(G)), intent(in), optional :: p_surf

!   Frazil formation keeps the temperature above the freezing point.
! This subroutine warms any water that is colder than the (currently
! surface) freezing point up to the freezing point and accumulates
! the required heat (in J m-2) in tv%frazil.
!   The expression, below, for the freezing point of sea water comes
! from Millero (1978) via Appendix A of Gill, 1982.

! Arguments: h - Layer thickness, in m or kg m-2.
!  (in/out)  tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 diabatic_driver_init.
  real, dimension(SZI_(G)) :: &
    fraz_col, & ! The accumulated heat requirement due to frazil, in J.
    t_freeze, & ! The freezing potential temperature at the current salinity, C.
    ps          ! pressure
  real, dimension(SZI_(G),SZK_(G)) :: &
    pressure    ! The pressure at the middle of each layer in Pa.
  real :: hc    ! A layer's heat capacity in J m-2 K-1.
  logical :: t_fr_set  ! True if the freezing point has been calculated for a
                       ! row of points.
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  call cpu_clock_begin(id_clock_frazil)

  if (.not.cs%pressure_dependent_frazil) then
    do k=1,nz ; do i=is,ie ; pressure(i,k) = 0.0 ; enddo ; enddo
  endif
!$OMP parallel do default(none) shared(is,ie,js,je,CS,G,GV,h,nz,tv,p_surf) &
!$OMP                           private(fraz_col,T_fr_set,T_freeze,hc,ps,pressure)
  do j=js,je
    ps(:) = 0.0
    if (PRESENT(p_surf)) then ; do i=is,ie
      ps(i) = p_surf(i,j)
    enddo ; endif

    do i=is,ie ; fraz_col(i) = 0.0 ; enddo

    if (cs%pressure_dependent_frazil) then
      do i=is,ie
        pressure(i,1) = ps(i) + (0.5*gv%H_to_Pa)*h(i,j,1)
      enddo
      do k=2,nz ; do i=is,ie
        pressure(i,k) = pressure(i,k-1) + &
          (0.5*gv%H_to_Pa) * (h(i,j,k) + h(i,j,k-1))
      enddo ; enddo
    endif

    if (cs%reclaim_frazil) then
      t_fr_set = .false.
      do i=is,ie ; if (tv%frazil(i,j) > 0.0) then
        if (.not.t_fr_set) then
          call calculate_tfreeze(tv%S(i:,j,1), pressure(i:,1), t_freeze(i:), &
                                 1, ie-i+1, tv%eqn_of_state)
          t_fr_set = .true.
        endif

        if (tv%T(i,j,1) > t_freeze(i)) then
    ! If frazil had previously been formed, but the surface temperature is now
    ! above freezing, cool the surface layer with the frazil heat deficit.
          hc = (tv%C_p*gv%H_to_kg_m2) * h(i,j,1)
          if (tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i)) <= 0.0) then
            tv%T(i,j,1) = tv%T(i,j,1) - tv%frazil(i,j)/hc
            tv%frazil(i,j) = 0.0
          else
            tv%frazil(i,j) = tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i))
            tv%T(i,j,1) = t_freeze(i)
          endif
        endif
      endif ; enddo
    endif

    do k=nz,1,-1
      t_fr_set = .false.
      do i=is,ie
        if ((g%mask2dT(i,j) > 0.0) .and. &
            ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then
          if (.not.t_fr_set) then
            call calculate_tfreeze(tv%S(i:,j,k), pressure(i:,k), t_freeze(i:), &
                                   1, ie-i+1, tv%eqn_of_state)
            t_fr_set = .true.
          endif

          hc = (tv%C_p*gv%H_to_kg_m2) * h(i,j,k)
          if (h(i,j,k) <= 10.0*gv%Angstrom) then
            ! Very thin layers should not be cooled by the frazil flux.
            if (tv%T(i,j,k) < t_freeze(i)) then
              fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
              tv%T(i,j,k) = t_freeze(i)
            endif
          else
            if (fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k)) <= 0.0) then
              tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i)/hc
              fraz_col(i) = 0.0
            else
              fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
              tv%T(i,j,k) = t_freeze(i)
            endif
          endif
        endif
      enddo
    enddo
    do i=is,ie
      tv%frazil(i,j) = tv%frazil(i,j) + fraz_col(i)
    enddo
  enddo
  call cpu_clock_end(id_clock_frazil)

tridiagts()

subroutine, public mom_diabatic_aux::tridiagts	(	type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		integer, intent(in)	is,
		integer, intent(in)	ie,
		integer, intent(in)	js,
		integer, intent(in)	je,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	hold,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	ea,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	eb,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	T,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	S
	)

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure

Definition at line 559 of file MOM_diabatic_aux.F90.

! Simple tri-diagnonal solver for T and S
! "Simple" means it only uses arrays hold, ea and eb
  ! Arguments
  type(ocean_grid_type),                    intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),                  intent(in)    :: gv   !< The ocean's vertical grid structure
  integer,                                  intent(in)    :: is, ie, js, je
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: hold, ea, eb
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: t, s
  ! Local variables
  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 :: h_tr, b_denom_1
  integer :: i, j, k
!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,hold,eb,T,S,ea) &
!$OMP                          private(h_tr,b1,d1,c1,b_denom_1)
  do j=js,je
    do i=is,ie
      h_tr = hold(i,j,1) + gv%H_subroundoff
      b1(i) = 1.0 / (h_tr + eb(i,j,1))
      d1(i) = h_tr * b1(i)
      t(i,j,1) = (b1(i)*h_tr)*t(i,j,1)
      s(i,j,1) = (b1(i)*h_tr)*s(i,j,1)
    enddo
    do k=2,g%ke ; do i=is,ie
      c1(i,k) = eb(i,j,k-1) * b1(i)
      h_tr = hold(i,j,k) + gv%H_subroundoff
      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)
      t(i,j,k) = b1(i) * (h_tr*t(i,j,k) + ea(i,j,k)*t(i,j,k-1))
      s(i,j,k) = b1(i) * (h_tr*s(i,j,k) + ea(i,j,k)*s(i,j,k-1))
    enddo ; enddo
    do k=g%ke-1,1,-1 ; do i=is,ie
      t(i,j,k) = t(i,j,k) + c1(i,k+1)*t(i,j,k+1)
      s(i,j,k) = s(i,j,k) + c1(i,k+1)*s(i,j,k+1)
    enddo ; enddo
  enddo

Variable Documentation

id_clock_frazil

integer mom_diabatic_aux::id_clock_frazil

private

Definition at line 133 of file MOM_diabatic_aux.F90.

Referenced by diabatic_aux_init(), and make_frazil().

id_clock_uv_at_h

integer mom_diabatic_aux::id_clock_uv_at_h

Definition at line 133 of file MOM_diabatic_aux.F90.

Referenced by diabatic_aux_init(), and find_uv_at_h().

integer :: id_clock_uv_at_h, id_clock_frazil