MOM6
Data Types | Functions/Subroutines
mom_diabatic_aux Module Reference

Detailed Description

Provides functions for some diabatic processes such as fraxil, brine rejection, tendency due to surface flux divergence.

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.

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.

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_aux_cs
 Control structure for diabatic_aux. More...
 

Functions/Subroutines

subroutine, public make_frazil (h, tv, G, GV, US, CS, p_surf, halo)
 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 [Q R Z ~> J m-2]) in tvfrazil. More...
 
subroutine, public differential_diffuse_t_s (h, T, S, Kd_T, Kd_S, dt, G, GV)
 This subroutine applies double diffusion to T & S, assuming no diapycal mass fluxes, using a simple triadiagonal solver. More...
 
subroutine, public adjust_salt (h, tv, G, GV, CS, halo)
 This subroutine keeps salinity from falling below a small but positive threshold. This usually occurs when the ice model attempts to extract more salt then is actually available to it from the ocean. More...
 
subroutine, public tridiagts (G, GV, is, ie, js, je, hold, ea, eb, T, S)
 This is a simple tri-diagonal solver for T and S. "Simple" means it only uses arrays hold, ea and eb. More...
 
subroutine, public find_uv_at_h (u, v, h, u_h, v_h, G, GV, US, ea, eb)
 This subroutine calculates u_h and v_h (velocities at thickness points), optionally using the entrainment amounts passed in as arguments. More...
 
subroutine, public set_pen_shortwave (optics, fluxes, G, GV, US, CS, opacity_CSp, tracer_flow_CSp)
 
subroutine, public diagnosemldbydensitydifference (id_MLD, h, tv, densityDiff, G, GV, US, diagPtr, id_N2subML, id_MLDsq, dz_subML)
 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 diagnosemldbyenergy (id_MLD, h, tv, G, GV, US, Mixing_Energy, diagPtr)
 Diagnose a mixed layer depth (MLD) determined by the depth a given energy value would mix. This routine is appropriate in MOM_diabatic_driver due to its position within the time stepping. More...
 
subroutine, public applyboundaryfluxesinout (CS, G, GV, US, dt, fluxes, optics, nsw, 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, US, param_file, diag, CS, useALEalgorithm, use_ePBL)
 This subroutine initializes the parameters and control structure of the diabatic_aux module. More...
 
subroutine, public diabatic_aux_end (CS)
 This subroutine initializes the control structure and any related memory for the diabatic_aux module. More...
 

Function/Subroutine Documentation

◆ adjust_salt()

subroutine, public mom_diabatic_aux::adjust_salt ( 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,
integer, intent(in), optional  halo 
)

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

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]hLayer thicknesses [H ~> m or kg m-2]
[in,out]tvStructure containing pointers to any available thermodynamic fields.
[in]csThe control structure returned by a previous call to diabatic_aux_init.
[in]haloHalo width over which to work

Definition at line 322 of file MOM_diabatic_aux.F90.

323  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
324  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
325  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
326  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
327  type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any
328  !! available thermodynamic fields.
329  type(diabatic_aux_CS), intent(in) :: CS !< The control structure returned by a previous
330  !! call to diabatic_aux_init.
331  integer, optional, intent(in) :: halo !< Halo width over which to work
332 
333  ! local variables
334  real :: salt_add_col(SZI_(G),SZJ_(G)) !< The accumulated salt requirement [ppt R Z ~> gSalt m-2]
335  real :: S_min !< The minimum salinity [ppt].
336  real :: mc !< A layer's mass [R Z ~> kg m-2].
337  integer :: i, j, k, is, ie, js, je, nz
338 
339  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
340  if (present(halo)) then
341  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
342  endif
343 
344 ! call cpu_clock_begin(id_clock_adjust_salt)
345 
346  s_min = tv%min_salinity
347 
348  salt_add_col(:,:) = 0.0
349 
350  !$OMP parallel do default(shared) private(mc)
351  do j=js,je
352  do k=nz,1,-1 ; do i=is,ie
353  if ( (g%mask2dT(i,j) > 0.0) .and. &
354  ((tv%S(i,j,k) < s_min) .or. (salt_add_col(i,j) > 0.0)) ) then
355  mc = gv%H_to_RZ * h(i,j,k)
356  if (h(i,j,k) <= 10.0*gv%Angstrom_H) then
357  ! Very thin layers should not be adjusted by the salt flux
358  if (tv%S(i,j,k) < s_min) then
359  salt_add_col(i,j) = salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k))
360  tv%S(i,j,k) = s_min
361  endif
362  elseif (salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k)) <= 0.0) then
363  tv%S(i,j,k) = tv%S(i,j,k) - salt_add_col(i,j) / mc
364  salt_add_col(i,j) = 0.0
365  else
366  salt_add_col(i,j) = salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k))
367  tv%S(i,j,k) = s_min
368  endif
369  endif
370  enddo ; enddo
371  do i=is,ie
372  tv%salt_deficit(i,j) = tv%salt_deficit(i,j) + salt_add_col(i,j)
373  enddo
374  enddo
375 ! call cpu_clock_end(id_clock_adjust_salt)
376 

◆ 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,
type(unit_scale_type), intent(in)  US,
real, intent(in)  dt,
type(forcing), intent(inout)  fluxes,
type(optics_type), pointer  optics,
integer, intent(in)  nsw,
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
csControl structure for diabatic_aux
[in]gGrid structure
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in]dtTime-step over which forcing is applied [T ~> s]
[in,out]fluxesSurface fluxes container
opticsOptical properties container
[in]nswThe number of frequency bands of penetrating shortwave radiation
[in,out]hLayer thickness [H ~> m or kg m-2]
[in,out]tvStructure containing pointers to any available thermodynamic fields.
[in]aggregate_fw_forcingIf False, treat in/out fluxes separately.
[in]evap_cfl_limitThe largest fraction of a layer that can be evaporated in one time-step [nondim].
[in]minimum_forcing_depthThe smallest depth over which heat and freshwater fluxes is applied [H ~> m or kg m-2].
[out]ctkeTurbulent kinetic energy requirement to mix
[out]dsv_dtPartial derivative of specific volume with
[out]dsv_dsPartial derivative of specific volume with
[out]skinbuoyfluxBuoyancy flux at surface [Z2 T-3 ~> m2 s-3].

Definition at line 928 of file MOM_diabatic_aux.F90.

932  type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
933  type(ocean_grid_type), intent(in) :: G !< Grid structure
934  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
935  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
936  real, intent(in) :: dt !< Time-step over which forcing is applied [T ~> s]
937  type(forcing), intent(inout) :: fluxes !< Surface fluxes container
938  type(optics_type), pointer :: optics !< Optical properties container
939  integer, intent(in) :: nsw !< The number of frequency bands of penetrating
940  !! shortwave radiation
941  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
942  intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
943  type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any
944  !! available thermodynamic fields.
945  logical, intent(in) :: aggregate_FW_forcing !< If False, treat in/out fluxes separately.
946  real, intent(in) :: evap_CFL_limit !< The largest fraction of a layer that
947  !! can be evaporated in one time-step [nondim].
948  real, intent(in) :: minimum_forcing_depth !< The smallest depth over which
949  !! heat and freshwater fluxes is applied [H ~> m or kg m-2].
950  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
951  optional, intent(out) :: cTKE !< Turbulent kinetic energy requirement to mix
952  !! forcing through each layer [R Z3 T-2 ~> J m-2]
953  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
954  optional, intent(out) :: dSV_dT !< Partial derivative of specific volume with
955  !! potential temperature [R-1 degC-1 ~> m3 kg-1 degC-1].
956  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
957  optional, intent(out) :: dSV_dS !< Partial derivative of specific volume with
958  !! salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
959  real, dimension(SZI_(G),SZJ_(G)), &
960  optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3].
961 
962  ! Local variables
963  integer, parameter :: maxGroundings = 5
964  integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings)
965  real :: H_limit_fluxes
966  real :: IforcingDepthScale
967  real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
968  real :: dThickness, dTemp, dSalt
969  real :: fractionOfForcing, hOld, Ithickness
970  real :: RivermixConst ! A constant used in implementing river mixing [R Z2 T-1 ~> Pa s].
971 
972  real, dimension(SZI_(G)) :: &
973  d_pres, & ! pressure change across a layer [R L2 T-2 ~> Pa]
974  p_lay, & ! average pressure in a layer [R L2 T-2 ~> Pa]
975  pres, & ! pressure at an interface [R L2 T-2 ~> Pa]
976  netMassInOut, & ! surface water fluxes [H ~> m or kg m-2] over time step
977  netMassIn, & ! mass entering ocean surface [H ~> m or kg m-2] over a time step
978  netMassOut, & ! mass leaving ocean surface [H ~> m or kg m-2] over a time step
979  netHeat, & ! heat via surface fluxes excluding Pen_SW_bnd and netMassOut
980  ! [degC H ~> degC m or degC kg m-2]
981  netsalt, & ! surface salt flux ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
982  ! [ppt H ~> ppt m or ppt kg m-2]
983  nonpensw, & ! non-downwelling SW, which is absorbed at ocean surface
984  ! [degC H ~> degC m or degC kg m-2]
985  surfpressure, & ! Surface pressure (approximated as 0.0) [R L2 T-2 ~> Pa]
986  drhodt, & ! change in density per change in temperature [R degC-1 ~> kg m-3 degC-1]
987  drhods, & ! change in density per change in salinity [R ppt-1 ~> kg m-3 ppt-1]
988  netheat_rate, & ! netheat but for dt=1 [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
989  netsalt_rate, & ! netsalt but for dt=1 (e.g. returns a rate)
990  ! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
991  netmassinout_rate! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1]
992  real, dimension(SZI_(G), SZK_(G)) :: &
993  h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2]
994  T2d, & ! A 2-d copy of the layer temperatures [degC]
995  pen_TKE_2d, & ! The TKE required to homogenize the heating by shortwave radiation within
996  ! a layer [R Z3 T-2 ~> J m-2]
997  dsv_dt_2d ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
998  real, dimension(SZI_(G)) :: &
999  netPen_rate ! The surface penetrative shortwave heating rate summed over all bands
1000  ! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1001  real, dimension(max(nsw,1),SZI_(G)) :: &
1002  Pen_SW_bnd, & ! The penetrative shortwave heating integrated over a timestep by band
1003  ! [degC H ~> degC m or degC kg m-2]
1004  pen_sw_bnd_rate ! The penetrative shortwave heating rate by band
1005  ! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1006  real, dimension(max(nsw,1),SZI_(G),SZK_(G)) :: &
1007  opacityBand ! The opacity (inverse of the exponential absorption length) of each frequency
1008  ! band of shortwave radation in each layer [H-1 ~> m-1 or m2 kg-1]
1009  real, dimension(maxGroundings) :: hGrounding ! Thickness added by each grounding event [H ~> m or kg m-2]
1010  real :: Temp_in, Salin_in
1011  real :: g_Hconv2 ! A conversion factor for use in the TKE calculation
1012  ! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
1013  real :: GoRho ! g_Earth times a unit conversion factor divided by density
1014  ! [Z T-2 R-1 ~> m4 s-2 kg-1]
1015  logical :: calculate_energetics
1016  logical :: calculate_buoyancy
1017  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1018  integer :: i, j, is, ie, js, je, k, nz, n, nb
1019  character(len=45) :: mesg
1020 
1021  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1022 
1023  idt = 1.0 / dt
1024 
1025  calculate_energetics = (present(ctke) .and. present(dsv_dt) .and. present(dsv_ds))
1026  calculate_buoyancy = present(skinbuoyflux)
1027  if (calculate_buoyancy) skinbuoyflux(:,:) = 0.0
1028  if (present(ctke)) ctke(:,:,:) = 0.0
1029  g_hconv2 = (us%L_to_Z**2*gv%g_Earth * gv%H_to_RZ) * gv%H_to_RZ
1030  eosdom(:) = eos_domain(g%HI)
1031 
1032  ! Only apply forcing if fluxes%sw is associated.
1033  if (.not.associated(fluxes%sw) .and. .not.calculate_energetics) return
1034 
1035  if (calculate_buoyancy) then
1036  surfpressure(:) = 0.0
1037  gorho = us%L_to_Z**2*gv%g_Earth / gv%Rho0
1038  endif
1039 
1040  ! H_limit_fluxes is used by extractFluxes1d to scale down fluxes if the total
1041  ! depth of the ocean is vanishing. It does not (yet) handle a value of zero.
1042  ! To accommodate vanishing upper layers, we need to allow for an instantaneous
1043  ! distribution of forcing over some finite vertical extent. The bulk mixed layer
1044  ! code handles this issue properly.
1045  h_limit_fluxes = max(gv%Angstrom_H, 1.e-30*gv%m_to_H)
1046 
1047  ! diagnostic to see if need to create mass to avoid grounding
1048  if (cs%id_createdH>0) cs%createdH(:,:) = 0.
1049  numberofgroundings = 0
1050 
1051  !$OMP parallel do default(none) shared(is,ie,js,je,nz,h,tv,nsw,G,GV,US,optics,fluxes, &
1052  !$OMP H_limit_fluxes,numberOfGroundings,iGround,jGround,&
1053  !$OMP nonPenSW,hGrounding,CS,Idt,aggregate_FW_forcing, &
1054  !$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, &
1055  !$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho, &
1056  !$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2) &
1057  !$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, &
1058  !$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, &
1059  !$OMP IforcingDepthScale, &
1060  !$OMP dThickness,dTemp,dSalt,hOld,Ithickness, &
1061  !$OMP netMassIn,pres,d_pres,p_lay,dSV_dT_2d, &
1062  !$OMP netmassinout_rate,netheat_rate,netsalt_rate, &
1063  !$OMP drhodt,drhods,pen_sw_bnd_rate, &
1064  !$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst) &
1065  !$OMP firstprivate(SurfPressure)
1066  do j=js,je
1067  ! Work in vertical slices for efficiency
1068 
1069  ! Copy state into 2D-slice arrays
1070  do k=1,nz ; do i=is,ie
1071  h2d(i,k) = h(i,j,k)
1072  t2d(i,k) = tv%T(i,j,k)
1073  enddo ; enddo
1074 
1075  if (calculate_energetics) then
1076  ! The partial derivatives of specific volume with temperature and
1077  ! salinity need to be precalculated to avoid having heating of
1078  ! tiny layers give nonsensical values.
1079  if (associated(tv%p_surf)) then
1080  do i=is,ie ; pres(i) = tv%p_surf(i,j) ; enddo
1081  else
1082  do i=is,ie ; pres(i) = 0.0 ; enddo
1083  endif
1084  do k=1,nz
1085  do i=is,ie
1086  d_pres(i) = (gv%g_Earth * gv%H_to_RZ) * h2d(i,k)
1087  p_lay(i) = pres(i) + 0.5*d_pres(i)
1088  pres(i) = pres(i) + d_pres(i)
1089  enddo
1090  call calculate_specific_vol_derivs(t2d(:,k), tv%S(:,j,k), p_lay(:), &
1091  dsv_dt(:,j,k), dsv_ds(:,j,k), tv%eqn_of_state, eosdom)
1092  do i=is,ie ; dsv_dt_2d(i,k) = dsv_dt(i,j,k) ; enddo
1093  enddo
1094  pen_tke_2d(:,:) = 0.0
1095  endif
1096 
1097  ! Nothing more is done on this j-slice if there is no buoyancy forcing.
1098  if (.not.associated(fluxes%sw)) cycle
1099 
1100  if (nsw>0) call extract_optics_slice(optics, j, g, gv, opacity=opacityband, opacity_scale=(1.0/gv%m_to_H))
1101 
1102  ! The surface forcing is contained in the fluxes type.
1103  ! We aggregate the thermodynamic forcing for a time step into the following:
1104  ! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step
1105  ! = lprec + fprec + vprec + evap + lrunoff + frunoff
1106  ! note that lprec generally has sea ice melt/form included.
1107  ! netMassOut = net mass leaving ocean surface [H ~> m or kg m-2] over a time step.
1108  ! netMassOut < 0 means mass leaves ocean.
1109  ! netHeat = heat via surface fluxes [degC H ~> degC m or degC kg m-2], excluding the part
1110  ! contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0.
1111  ! netSalt = surface salt fluxes [ppt H ~> dppt m or gSalt m-2]
1112  ! Pen_SW_bnd = components to penetrative shortwave radiation split according to bands.
1113  ! This field provides that portion of SW from atmosphere that in fact
1114  ! enters to the ocean and participates in pentrative SW heating.
1115  ! nonpenSW = non-downwelling SW flux, which is absorbed in ocean surface
1116  ! (in tandem w/ LW,SENS,LAT); saved only for diagnostic purposes.
1117 
1118  !----------------------------------------------------------------------------------------
1119  !BGR-June 26, 2017{
1120  !Temporary action to preserve answers while fixing a bug.
1121  ! To fix a bug in a diagnostic calculation, applyboundaryfluxesinout now returns
1122  ! the surface buoyancy flux. Previously, extractbuoyancyflux2d was called, meaning
1123  ! a second call to extractfluxes1d (causing the diagnostic net_heat to be incorrect).
1124  ! Note that this call to extract buoyancyflux2d was AFTER applyboundaryfluxesinout,
1125  ! which means it used the T/S fields after this routine. Therefore, the surface
1126  ! buoyancy flux is computed here at the very end of this routine for legacy reasons.
1127  ! A few specific notes follow:
1128  ! 1) The old method did not included river/calving contributions to heat flux. This
1129  ! is kept consistent here via commenting code in the present extractFluxes1d <_rate>
1130  ! outputs, but we may reconsider this approach.
1131  ! 2) The old method computed the buoyancy flux rate directly (by setting dt=1), instead
1132  ! of computing the integrated value (and dividing by dt). Hence the required
1133  ! additional outputs from extractFluxes1d.
1134  ! *** This is because: A*dt/dt =/= A due to round off.
1135  ! 3) The old method computed buoyancy flux after this routine, meaning the returned
1136  ! surface fluxes (from extractfluxes1d) must be recorded for use later in the code.
1137  ! We could (and maybe should) move that loop up to before the surface fluxes are
1138  ! applied, but this will change answers.
1139  ! For all these reasons we compute additional values of <_rate> which are preserved
1140  ! for the buoyancy flux calculation and reproduce the old answers.
1141  ! In the future this needs more detailed investigation to make sure everything is
1142  ! consistent and correct. These details shouldnt significantly effect climate,
1143  ! but do change answers.
1144  !-----------------------------------------------------------------------------------------
1145  if (calculate_buoyancy) then
1146  call extractfluxes1d(g, gv, us, fluxes, optics, nsw, j, dt, &
1147  h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
1148  h2d, t2d, netmassinout, netmassout, netheat, netsalt, &
1149  pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw, &
1150  net_heat_rate=netheat_rate, net_salt_rate=netsalt_rate, &
1151  netmassinout_rate=netmassinout_rate, pen_sw_bnd_rate=pen_sw_bnd_rate)
1152  else
1153  call extractfluxes1d(g, gv, us, fluxes, optics, nsw, j, dt, &
1154  h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
1155  h2d, t2d, netmassinout, netmassout, netheat, netsalt, &
1156  pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw)
1157  endif
1158  ! ea is for passive tracers
1159  do i=is,ie
1160  ! ea(i,j,1) = netMassInOut(i)
1161  if (aggregate_fw_forcing) then
1162  netmassout(i) = netmassinout(i)
1163  netmassin(i) = 0.
1164  else
1165  netmassin(i) = netmassinout(i) - netmassout(i)
1166  endif
1167  if (g%mask2dT(i,j)>0.0) then
1168  fluxes%netMassOut(i,j) = netmassout(i)
1169  fluxes%netMassIn(i,j) = netmassin(i)
1170  else
1171  fluxes%netMassOut(i,j) = 0.0
1172  fluxes%netMassIn(i,j) = 0.0
1173  endif
1174  enddo
1175 
1176  ! Apply the surface boundary fluxes in three steps:
1177  ! A/ update mass, temp, and salinity due to all terms except mass leaving
1178  ! ocean (and corresponding outward heat content), and ignoring penetrative SW.
1179  ! B/ update mass, salt, temp from mass leaving ocean.
1180  ! C/ update temp due to penetrative SW
1181  do i=is,ie
1182  if (g%mask2dT(i,j)>0.) then
1183 
1184  ! A/ Update mass, temp, and salinity due to incoming mass flux.
1185  do k=1,1
1186 
1187  ! Change in state due to forcing
1188  dthickness = netmassin(i) ! Since we are adding mass, we can use all of it
1189  dtemp = 0.
1190  dsalt = 0.
1191 
1192  ! Update the forcing by the part to be consumed within the present k-layer.
1193  ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
1194  netmassin(i) = netmassin(i) - dthickness
1195  ! This line accounts for the temperature of the mass exchange
1196  temp_in = t2d(i,k)
1197  salin_in = 0.0
1198  dtemp = dtemp + dthickness*temp_in
1199 
1200  ! Diagnostics of heat content associated with mass fluxes
1201  if (associated(fluxes%heat_content_massin)) &
1202  fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1203  t2d(i,k) * max(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
1204  if (associated(fluxes%heat_content_massout)) &
1205  fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
1206  t2d(i,k) * min(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
1207  if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1208  t2d(i,k) * dthickness * gv%H_to_RZ
1209 
1210  ! Determine the energetics of river mixing before updating the state.
1211  if (calculate_energetics .and. associated(fluxes%lrunoff) .and. cs%do_rivermix) then
1212  ! Here we add an additional source of TKE to the mixed layer where river
1213  ! is present to simulate unresolved estuaries. The TKE input, TKE_river in
1214  ! [Z3 T-3 ~> m3 s-3], is diagnosed as follows:
1215  ! TKE_river = 0.5*rivermix_depth*g*(1/rho)*drho_ds*
1216  ! River*(Samb - Sriver) = CS%mstar*U_star^3
1217  ! where River is in units of [Z T-1 ~> m s-1].
1218  ! Samb = Ambient salinity at the mouth of the estuary
1219  ! rivermix_depth = The prescribed depth over which to mix river inflow
1220  ! drho_ds = The gradient of density wrt salt at the ambient surface salinity.
1221  ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
1222  if (gv%Boussinesq) then
1223  rivermixconst = -0.5*(cs%rivermix_depth*dt) * ( us%L_to_Z**2*gv%g_Earth ) * gv%Rho0
1224  else
1225  rivermixconst = -0.5*(cs%rivermix_depth*dt) * gv%Rho0 * ( us%L_to_Z**2*gv%g_Earth )
1226  endif
1227  ctke(i,j,k) = ctke(i,j,k) + max(0.0, rivermixconst*dsv_ds(i,j,1) * &
1228  (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * tv%S(i,j,1))
1229  endif
1230 
1231  ! Update state
1232  hold = h2d(i,k) ! Keep original thickness in hand
1233  h2d(i,k) = h2d(i,k) + dthickness ! New thickness
1234  if (h2d(i,k) > 0.0) then
1235  if (calculate_energetics .and. (dthickness > 0.)) then
1236  ! Calculate the energy required to mix the newly added water over
1237  ! the topmost grid cell.
1238  ctke(i,j,k) = ctke(i,j,k) + 0.5*g_hconv2*(hold*dthickness) * &
1239  ((t2d(i,k) - temp_in) * dsv_dt(i,j,k) + (tv%S(i,j,k) - salin_in) * dsv_ds(i,j,k))
1240  endif
1241  ithickness = 1.0/h2d(i,k) ! Inverse new thickness
1242  ! The "if"s below avoid changing T/S by roundoff unnecessarily
1243  if (dthickness /= 0. .or. dtemp /= 0.) t2d(i,k) = (hold*t2d(i,k) + dtemp)*ithickness
1244  if (dthickness /= 0. .or. dsalt /= 0.) tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness
1245 
1246  endif
1247 
1248  enddo ! k=1,1
1249 
1250  ! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt.
1251  do k=1,nz
1252 
1253  ! Place forcing into this layer if this layer has nontrivial thickness.
1254  ! For layers thin relative to 1/IforcingDepthScale, then distribute
1255  ! forcing into deeper layers.
1256  iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth - netmassout(i) )
1257  ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
1258  fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale)
1259 
1260  ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we
1261  ! limit the forcing applied to this cell, leaving the remaining forcing to
1262  ! be distributed downwards.
1263  if (-fractionofforcing*netmassout(i) > evap_cfl_limit*h2d(i,k)) then
1264  fractionofforcing = -evap_cfl_limit*h2d(i,k)/netmassout(i)
1265  endif
1266 
1267  ! Change in state due to forcing
1268 
1269  dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) )
1270  dtemp = fractionofforcing*netheat(i)
1271  ! ### The 0.9999 here should become a run-time parameter?
1272  dsalt = max( fractionofforcing*netsalt(i), -0.9999*h2d(i,k)*tv%S(i,j,k))
1273 
1274  ! Update the forcing by the part to be consumed within the present k-layer.
1275  ! If fractionOfForcing = 1, then new netMassOut vanishes.
1276  netmassout(i) = netmassout(i) - dthickness
1277  netheat(i) = netheat(i) - dtemp
1278  netsalt(i) = netsalt(i) - dsalt
1279 
1280  ! This line accounts for the temperature of the mass exchange
1281  dtemp = dtemp + dthickness*t2d(i,k)
1282 
1283  ! Diagnostics of heat content associated with mass fluxes
1284  if (associated(fluxes%heat_content_massin)) &
1285  fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1286  t2d(i,k) * max(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
1287  if (associated(fluxes%heat_content_massout)) &
1288  fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
1289  t2d(i,k) * min(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
1290  if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1291  t2d(i,k) * dthickness * gv%H_to_RZ
1292 
1293  ! Update state by the appropriate increment.
1294  hold = h2d(i,k) ! Keep original thickness in hand
1295  h2d(i,k) = h2d(i,k) + dthickness ! New thickness
1296 
1297  if (h2d(i,k) > 0.) then
1298  if (calculate_energetics) then
1299  ! Calculate the energy required to mix the newly added water over the topmost grid
1300  ! cell, assuming that the fluxes of heat and salt and rejected brine are initially
1301  ! applied in vanishingly thin layers at the top of the layer before being mixed
1302  ! throughout the layer. Note that dThickness is always <= 0 here, and that
1303  ! negative cTKE is a deficit that will need to be filled later.
1304  ctke(i,j,k) = ctke(i,j,k) - (0.5*h2d(i,k)*g_hconv2) * &
1305  ((dtemp - dthickness*t2d(i,k)) * dsv_dt(i,j,k) + &
1306  (dsalt - dthickness*tv%S(i,j,k)) * dsv_ds(i,j,k))
1307  endif
1308  ithickness = 1.0/h2d(i,k) ! Inverse of new thickness
1309  t2d(i,k) = (hold*t2d(i,k) + dtemp)*ithickness
1310  tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness
1311  elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling
1312  call forcing_singlepointprint(fluxes,g,i,j,'applyBoundaryFluxesInOut (h<0)')
1313  write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
1314  write(0,*) 'applyBoundaryFluxesInOut(): netT,netS,netH=',netheat(i),netsalt(i),netmassinout(i)
1315  write(0,*) 'applyBoundaryFluxesInOut(): dT,dS,dH=',dtemp,dsalt,dthickness
1316  write(0,*) 'applyBoundaryFluxesInOut(): h(n),h(n+1),k=',hold,h2d(i,k),k
1317  call mom_error(fatal, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
1318  "Complete mass loss in column!")
1319  endif
1320 
1321  enddo ! k
1322 
1323  ! Check if trying to apply fluxes over land points
1324  elseif ((abs(netheat(i))+abs(netsalt(i))+abs(netmassin(i))+abs(netmassout(i)))>0.) then
1325 
1326  if (.not. cs%ignore_fluxes_over_land) then
1327  call forcing_singlepointprint(fluxes,g,i,j,'applyBoundaryFluxesInOut (land)')
1328  write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
1329  write(0,*) 'applyBoundaryFluxesInOut(): netHeat,netSalt,netMassIn,netMassOut=',&
1330  netheat(i),netsalt(i),netmassin(i),netmassout(i)
1331 
1332  call mom_error(fatal, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
1333  "Mass loss over land?")
1334  endif
1335 
1336  endif
1337 
1338  ! If anything remains after the k-loop, then we have grounded out, which is a problem.
1339  if (netmassin(i)+netmassout(i) /= 0.0) then
1340 !$OMP critical
1341  numberofgroundings = numberofgroundings +1
1342  if (numberofgroundings<=maxgroundings) then
1343  iground(numberofgroundings) = i ! Record i,j location of event for
1344  jground(numberofgroundings) = j ! warning message
1345  hgrounding(numberofgroundings) = netmassin(i)+netmassout(i)
1346  endif
1347 !$OMP end critical
1348  if (cs%id_createdH>0) cs%createdH(i,j) = cs%createdH(i,j) - (netmassin(i)+netmassout(i))/dt
1349  endif
1350 
1351  enddo ! i
1352 
1353  ! Step C/ in the application of fluxes
1354  ! Heat by the convergence of penetrating SW.
1355  ! SW penetrative heating uses the updated thickness from above.
1356 
1357  ! Save temperature before increment with SW heating
1358  ! and initialize CS%penSWflux_diag to zero.
1359  if (cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0) then
1360  do k=1,nz ; do i=is,ie
1361  cs%penSW_diag(i,j,k) = t2d(i,k)
1362  cs%penSWflux_diag(i,j,k) = 0.0
1363  enddo ; enddo
1364  k=nz+1 ; do i=is,ie
1365  cs%penSWflux_diag(i,j,k) = 0.0
1366  enddo
1367  endif
1368 
1369  if (calculate_energetics) then
1370  call absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h_limit_fluxes, &
1371  .false., .true., t2d, pen_sw_bnd, tke=pen_tke_2d, dsv_dt=dsv_dt_2d)
1372  k = 1 ! For setting break-points.
1373  do k=1,nz ; do i=is,ie
1374  ctke(i,j,k) = ctke(i,j,k) + pen_tke_2d(i,k)
1375  enddo ; enddo
1376  else
1377  call absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h_limit_fluxes, &
1378  .false., .true., t2d, pen_sw_bnd)
1379  endif
1380 
1381 
1382  ! Step D/ copy updated thickness and temperature
1383  ! 2d slice now back into model state.
1384  do k=1,nz ; do i=is,ie
1385  h(i,j,k) = h2d(i,k)
1386  tv%T(i,j,k) = t2d(i,k)
1387  enddo ; enddo
1388 
1389  ! Diagnose heating [Q R Z T-1 ~> W m-2] applied to a grid cell from SW penetration
1390  ! Also diagnose the penetrative SW heat flux at base of layer.
1391  if (cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0) then
1392 
1393  ! convergence of SW into a layer
1394  do k=1,nz ; do i=is,ie
1395  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_RZ
1396  enddo ; enddo
1397 
1398  ! Perform a cumulative sum upwards from bottom to
1399  ! diagnose penetrative SW flux at base of tracer cell.
1400  ! CS%penSWflux_diag(i,j,k=1) is penetrative shortwave at top of ocean.
1401  ! CS%penSWflux_diag(i,j,k=kbot+1) is zero, since assume no SW penetrates rock.
1402  ! CS%penSWflux_diag = rsdo and CS%penSW_diag = rsdoabsorb
1403  ! rsdoabsorb(k) = rsdo(k) - rsdo(k+1), so that rsdo(k) = rsdo(k+1) + rsdoabsorb(k)
1404  if (cs%id_penSWflux_diag > 0) then
1405  do k=nz,1,-1 ; do i=is,ie
1406  cs%penSWflux_diag(i,j,k) = cs%penSW_diag(i,j,k) + cs%penSWflux_diag(i,j,k+1)
1407  enddo ; enddo
1408  endif
1409 
1410  endif
1411 
1412  ! Fill CS%nonpenSW_diag
1413  if (cs%id_nonpenSW_diag > 0) then
1414  do i=is,ie
1415  cs%nonpenSW_diag(i,j) = nonpensw(i) * idt * tv%C_p * gv%H_to_RZ
1416  enddo
1417  endif
1418 
1419  ! BGR: Get buoyancy flux to return for ePBL
1420  ! We want the rate, so we use the rate values returned from extractfluxes1d.
1421  ! Note that the *dt values could be divided by dt here, but
1422  ! 1) Answers will change due to round-off
1423  ! 2) Be sure to save their values BEFORE fluxes are used.
1424  if (calculate_buoyancy) then
1425  drhodt(:) = 0.0
1426  drhods(:) = 0.0
1427  netpen_rate(:) = 0.0
1428  ! Sum over bands and attenuate as a function of depth.
1429  ! netPen_rate is the netSW as a function of depth, but only the surface value is used here,
1430  ! in which case the values of dt, h, optics and H_limit_fluxes are irrelevant. Consider
1431  ! writing a shorter and simpler variant to handle this very limited case.
1432  ! call sumSWoverBands(G, GV, US, h2d(:,:), optics_nbands(optics), optics, j, dt, &
1433  ! H_limit_fluxes, .true., pen_SW_bnd_rate, netPen)
1434  do i=is,ie ; do nb=1,nsw ; netpen_rate(i) = netpen_rate(i) + pen_sw_bnd_rate(nb,i) ; enddo ; enddo
1435 
1436  ! Density derivatives
1437  if (associated(tv%p_surf)) then ; do i=is,ie ; surfpressure(i) = tv%p_surf(i,j) ; enddo ; endif
1438  call calculate_density_derivs(t2d(:,1), tv%S(:,j,1), surfpressure, drhodt, drhods, &
1439  tv%eqn_of_state, eosdom)
1440  ! 1. Adjust netSalt to reflect dilution effect of FW flux
1441  ! 2. Add in the SW heating for purposes of calculating the net
1442  ! surface buoyancy flux affecting the top layer.
1443  ! 3. Convert to a buoyancy flux, excluding penetrating SW heating
1444  ! BGR-Jul 5, 2017: The contribution of SW heating here needs investigated for ePBL.
1445  do i=is,ie
1446  skinbuoyflux(i,j) = - gorho * gv%H_to_Z * &
1447  (drhods(i) * (netsalt_rate(i) - tv%S(i,j,1)*netmassinout_rate(i)) + &
1448  drhodt(i) * ( netheat_rate(i) + netpen_rate(i)) ) ! [Z2 T-3 ~> m2 s-3]
1449  enddo
1450  endif
1451 
1452  enddo ! j-loop finish
1453 
1454  ! Post the diagnostics
1455  if (cs%id_createdH > 0) call post_data(cs%id_createdH , cs%createdH , cs%diag)
1456  if (cs%id_penSW_diag > 0) call post_data(cs%id_penSW_diag , cs%penSW_diag , cs%diag)
1457  if (cs%id_penSWflux_diag > 0) call post_data(cs%id_penSWflux_diag, cs%penSWflux_diag, cs%diag)
1458  if (cs%id_nonpenSW_diag > 0) call post_data(cs%id_nonpenSW_diag , cs%nonpenSW_diag , cs%diag)
1459 
1460 ! The following check will be ignored if ignore_fluxes_over_land = true
1461  if (numberofgroundings>0 .and. .not. cs%ignore_fluxes_over_land) then
1462  do i = 1, min(numberofgroundings, maxgroundings)
1463  call forcing_singlepointprint(fluxes,g,iground(i),jground(i),'applyBoundaryFluxesInOut (grounding)')
1464  write(mesg(1:45),'(3es15.3)') g%geoLonT( iground(i), jground(i) ), &
1465  g%geoLatT( iground(i), jground(i)), hgrounding(i)*gv%H_to_m
1466  call mom_error(warning, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
1467  "Mass created. x,y,dh= "//trim(mesg), all_print=.true.)
1468  enddo
1469 
1470  if (numberofgroundings - maxgroundings > 0) then
1471  write(mesg, '(i4)') numberofgroundings - maxgroundings
1472  call mom_error(warning, "MOM_diabatic_driver:F90, applyBoundaryFluxesInOut(): "//&
1473  trim(mesg) // " groundings remaining")
1474  endif
1475  endif
1476 

◆ diabatic_aux_end()

subroutine, public mom_diabatic_aux::diabatic_aux_end ( type(diabatic_aux_cs), pointer  CS)

This subroutine initializes the control structure and any related memory for the diabatic_aux module.

Parameters
csThe control structure returned by a previous call to diabatic_aux_init; it is deallocated here.

Definition at line 1644 of file MOM_diabatic_aux.F90.

1645  type(diabatic_aux_CS), pointer :: CS !< The control structure returned by a previous
1646  !! call to diabatic_aux_init; it is deallocated here.
1647 
1648  if (.not.associated(cs)) return
1649 
1650  if (cs%id_createdH >0) deallocate(cs%createdH)
1651  if (cs%id_penSW_diag >0) deallocate(cs%penSW_diag)
1652  if (cs%id_penSWflux_diag >0) deallocate(cs%penSWflux_diag)
1653  if (cs%id_nonpenSW_diag >0) deallocate(cs%nonpenSW_diag)
1654 
1655  if (associated(cs)) deallocate(cs)
1656 

◆ diabatic_aux_init()

subroutine, public mom_diabatic_aux::diabatic_aux_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
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 
)

This subroutine initializes the parameters and control structure of the diabatic_aux module.

Parameters
[in]timeThe current model time.
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in]param_fileA structure to parse for run-time parameters
[in,out]diagA structure used to regulate diagnostic output
csA pointer to the control structure for the diabatic_aux module, which is initialized here.
[in]usealealgorithmIf true, use the ALE algorithm rather than layered mode.
[in]use_epblIf true, use the implicit energetics planetary boundary layer scheme to determine the diffusivity in the surface boundary layer.

Definition at line 1480 of file MOM_diabatic_aux.F90.

1481  type(time_type), target, intent(in) :: Time !< The current model time.
1482  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1483  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1484  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1485  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1486  type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output
1487  type(diabatic_aux_CS), pointer :: CS !< A pointer to the control structure for the
1488  !! diabatic_aux module, which is initialized here.
1489  logical, intent(in) :: useALEalgorithm !< If true, use the ALE algorithm rather
1490  !! than layered mode.
1491  logical, intent(in) :: use_ePBL !< If true, use the implicit energetics planetary
1492  !! boundary layer scheme to determine the diffusivity
1493  !! in the surface boundary layer.
1494 
1495 ! This "include" declares and sets the variable "version".
1496 #include "version_variable.h"
1497  character(len=40) :: mdl = "MOM_diabatic_aux" ! This module's name.
1498  character(len=48) :: thickness_units
1499  character(len=200) :: inputdir ! The directory where NetCDF input files
1500  character(len=240) :: chl_filename ! A file from which chl_a concentrations are to be read.
1501  character(len=128) :: chl_file ! Data containing chl_a concentrations. Used
1502  ! when var_pen_sw is defined and reading from file.
1503  character(len=32) :: chl_varname ! Name of chl_a variable in chl_file.
1504  logical :: use_temperature ! True if thermodynamics are enabled.
1505  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nbands
1506  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1507  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1508 
1509  if (associated(cs)) then
1510  call mom_error(warning, "diabatic_aux_init called with an "// &
1511  "associated control structure.")
1512  return
1513  else
1514  allocate(cs)
1515  endif
1516 
1517  cs%diag => diag
1518  cs%Time => time
1519 
1520 ! Set default, read and log parameters
1521  call log_version(param_file, mdl, version, &
1522  "The following parameters are used for auxiliary diabatic processes.")
1523 
1524  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
1525  "If true, temperature and salinity are used as state "//&
1526  "variables.", default=.true.)
1527 
1528  call get_param(param_file, mdl, "RECLAIM_FRAZIL", cs%reclaim_frazil, &
1529  "If true, try to use any frazil heat deficit to cool any "//&
1530  "overlying layers down to the freezing point, thereby "//&
1531  "avoiding the creation of thin ice when the SST is above "//&
1532  "the freezing point.", default=.true.)
1533  call get_param(param_file, mdl, "PRESSURE_DEPENDENT_FRAZIL", &
1534  cs%pressure_dependent_frazil, &
1535  "If true, use a pressure dependent freezing temperature "//&
1536  "when making frazil. The default is false, which will be "//&
1537  "faster but is inappropriate with ice-shelf cavities.", &
1538  default=.false.)
1539 
1540  if (use_epbl) then
1541  call get_param(param_file, mdl, "IGNORE_FLUXES_OVER_LAND", cs%ignore_fluxes_over_land,&
1542  "If true, the model does not check if fluxes are being applied "//&
1543  "over land points. This is needed when the ocean is coupled "//&
1544  "with ice shelves and sea ice, since the sea ice mask needs to "//&
1545  "be different than the ocean mask to avoid sea ice formation "//&
1546  "under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
1547  call get_param(param_file, mdl, "DO_RIVERMIX", cs%do_rivermix, &
1548  "If true, apply additional mixing wherever there is "//&
1549  "runoff, so that it is mixed down to RIVERMIX_DEPTH "//&
1550  "if the ocean is that deep.", default=.false.)
1551  if (cs%do_rivermix) &
1552  call get_param(param_file, mdl, "RIVERMIX_DEPTH", cs%rivermix_depth, &
1553  "The depth to which rivers are mixed if DO_RIVERMIX is "//&
1554  "defined.", units="m", default=0.0, scale=us%m_to_Z)
1555  else
1556  cs%do_rivermix = .false. ; cs%rivermix_depth = 0.0 ; cs%ignore_fluxes_over_land = .false.
1557  endif
1558 
1559  if (gv%nkml == 0) then
1560  call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
1561  "If true, use the fluxes%runoff_Hflx field to set the "//&
1562  "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
1563  default=.false.)
1564  call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
1565  "If true, use the fluxes%calving_Hflx field to set the "//&
1566  "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
1567  default=.false.)
1568  else
1569  cs%use_river_heat_content = .false.
1570  cs%use_calving_heat_content = .false.
1571  endif
1572 
1573  if (usealealgorithm) then
1574  cs%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, &
1575  time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", &
1576  "m s-1", conversion=gv%H_to_m*us%s_to_T)
1577  if (cs%id_createdH>0) allocate(cs%createdH(isd:ied,jsd:jed))
1578 
1579  ! diagnostic for heating of a grid cell from convergence of SW heat into the cell
1580  cs%id_penSW_diag = register_diag_field('ocean_model', 'rsdoabsorb', &
1581  diag%axesTL, time, 'Convergence of Penetrative Shortwave Flux in Sea Water Layer',&
1582  'W m-2', conversion=us%QRZ_T_to_W_m2, &
1583  standard_name='net_rate_of_absorption_of_shortwave_energy_in_ocean_layer', v_extensive=.true.)
1584 
1585  ! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces)
1586  ! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock).
1587  cs%id_penSWflux_diag = register_diag_field('ocean_model', 'rsdo', &
1588  diag%axesTi, time, 'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',&
1589  'W m-2', conversion=us%QRZ_T_to_W_m2, standard_name='downwelling_shortwave_flux_in_sea_water')
1590 
1591  ! need both arrays for the SW diagnostics (one for flux, one for convergence)
1592  if (cs%id_penSW_diag>0 .or. cs%id_penSWflux_diag>0) then
1593  allocate(cs%penSW_diag(isd:ied,jsd:jed,nz)) ; cs%penSW_diag(:,:,:) = 0.0
1594  allocate(cs%penSWflux_diag(isd:ied,jsd:jed,nz+1)) ; cs%penSWflux_diag(:,:,:) = 0.0
1595  endif
1596 
1597  ! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface)
1598  cs%id_nonpenSW_diag = register_diag_field('ocean_model', 'nonpenSW', &
1599  diag%axesT1, time, &
1600  'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',&
1601  'W m-2', conversion=us%QRZ_T_to_W_m2, &
1602  standard_name='nondownwelling_shortwave_flux_in_sea_water')
1603  if (cs%id_nonpenSW_diag > 0) then
1604  allocate(cs%nonpenSW_diag(isd:ied,jsd:jed)) ; cs%nonpenSW_diag(:,:) = 0.0
1605  endif
1606  endif
1607 
1608  if (use_temperature) then
1609  call get_param(param_file, mdl, "VAR_PEN_SW", cs%var_pen_sw, &
1610  "If true, use one of the CHL_A schemes specified by "//&
1611  "OPACITY_SCHEME to determine the e-folding depth of "//&
1612  "incoming short wave radiation.", default=.false.)
1613  if (cs%var_pen_sw) then
1614 
1615  call get_param(param_file, mdl, "CHL_FROM_FILE", cs%chl_from_file, &
1616  "If true, chl_a is read from a file.", default=.true.)
1617  if (cs%chl_from_file) then
1618  call time_interp_external_init()
1619 
1620  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1621  call get_param(param_file, mdl, "CHL_FILE", chl_file, &
1622  "CHL_FILE is the file containing chl_a concentrations in "//&
1623  "the variable CHL_A. It is used when VAR_PEN_SW and "//&
1624  "CHL_FROM_FILE are true.", fail_if_missing=.true.)
1625  chl_filename = trim(slasher(inputdir))//trim(chl_file)
1626  call log_param(param_file, mdl, "INPUTDIR/CHL_FILE", chl_filename)
1627  call get_param(param_file, mdl, "CHL_VARNAME", chl_varname, &
1628  "Name of CHL_A variable in CHL_FILE.", default='CHL_A')
1629  cs%sbc_chl = init_external_field(chl_filename, trim(chl_varname), domain=g%Domain%mpp_domain)
1630  endif
1631 
1632  cs%id_chl = register_diag_field('ocean_model', 'Chl_opac', diag%axesT1, time, &
1633  'Surface chlorophyll A concentration used to find opacity', 'mg m-3')
1634  endif
1635  endif
1636 
1637  id_clock_uv_at_h = cpu_clock_id('(Ocean find_uv_at_h)', grain=clock_routine)
1638  id_clock_frazil = cpu_clock_id('(Ocean frazil)', grain=clock_routine)
1639 

◆ 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(unit_scale_type), intent(in)  US,
type(diag_ctrl), pointer  diagPtr,
integer, intent(in), optional  id_N2subML,
integer, intent(in), optional  id_MLDsq,
real, intent(in), optional  dz_subML 
)

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]gGrid type
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in]id_mldHandle (ID) of MLD diagnostic
[in]hLayer thickness [H ~> m or kg m-2]
[in]tvStructure containing pointers to any available thermodynamic fields.
[in]densitydiffDensity difference to determine MLD [R ~> kg m-3]
diagptrDiagnostics structure
[in]id_n2submlOptional handle (ID) of subML stratification
[in]id_mldsqOptional handle (ID) of squared MLD
[in]dz_submlThe distance over which to calculate N2subML or 50 m if missing [Z ~> m]

Definition at line 591 of file MOM_diabatic_aux.F90.

593  type(ocean_grid_type), intent(in) :: G !< Grid type
594  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
595  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
596  integer, intent(in) :: id_MLD !< Handle (ID) of MLD diagnostic
597  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
598  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
599  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
600  !! available thermodynamic fields.
601  real, intent(in) :: densityDiff !< Density difference to determine MLD [R ~> kg m-3]
602  type(diag_ctrl), pointer :: diagPtr !< Diagnostics structure
603  integer, optional, intent(in) :: id_N2subML !< Optional handle (ID) of subML stratification
604  integer, optional, intent(in) :: id_MLDsq !< Optional handle (ID) of squared MLD
605  real, optional, intent(in) :: dz_subML !< The distance over which to calculate N2subML
606  !! or 50 m if missing [Z ~> m]
607 
608  ! Local variables
609  real, dimension(SZI_(G)) :: deltaRhoAtKm1, deltaRhoAtK ! Density differences [R ~> kg m-3].
610  real, dimension(SZI_(G)) :: pRef_MLD, pRef_N2 ! Reference pressures [R L2 T-2 ~> Pa].
611  real, dimension(SZI_(G)) :: H_subML, dH_N2 ! Summed thicknesses used in N2 calculation [H ~> m].
612  real, dimension(SZI_(G)) :: T_subML, T_deeper ! Temperatures used in the N2 calculation [degC].
613  real, dimension(SZI_(G)) :: S_subML, S_deeper ! Salinities used in the N2 calculation [ppt].
614  real, dimension(SZI_(G)) :: rho_subML, rho_deeper ! Densities used in the N2 calculation [R ~> kg m-3].
615  real, dimension(SZI_(G)) :: dK, dKm1 ! Depths [Z ~> m].
616  real, dimension(SZI_(G)) :: rhoSurf ! Density used in finding the mixedlayer depth [R ~> kg m-3].
617  real, dimension(SZI_(G), SZJ_(G)) :: MLD ! Diagnosed mixed layer depth [Z ~> m].
618  real, dimension(SZI_(G), SZJ_(G)) :: subMLN2 ! Diagnosed stratification below ML [T-2 ~> s-2].
619  real, dimension(SZI_(G), SZJ_(G)) :: MLD2 ! Diagnosed MLD^2 [Z2 ~> m2].
620  logical, dimension(SZI_(G)) :: N2_region_set ! If true, all necessary values for calculating N2
621  ! have been stored already.
622  real :: gE_Rho0 ! The gravitational acceleration divided by a mean density [Z T-2 R-1 ~> m4 s-2 kg-1].
623  real :: dH_subML ! Depth below ML over which to diagnose stratification [H ~> m].
624  real :: aFac ! A nondimensional factor [nondim]
625  real :: ddRho ! A density difference [R ~> kg m-3]
626  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
627  integer :: i, j, is, ie, js, je, k, nz, id_N2, id_SQ
628 
629  id_n2 = -1 ; if (PRESENT(id_n2subml)) id_n2 = id_n2subml
630 
631  id_sq = -1 ; if (PRESENT(id_mldsq)) id_sq = id_mldsq
632 
633  ge_rho0 = us%L_to_Z**2*gv%g_Earth / (gv%Rho0)
634  dh_subml = 50.*gv%m_to_H ; if (present(dz_subml)) dh_subml = gv%Z_to_H*dz_subml
635 
636  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
637 
638  pref_mld(:) = 0.0
639  eosdom(:) = eos_domain(g%HI)
640  do j=js,je
641  do i=is,ie ; dk(i) = 0.5 * h(i,j,1) * gv%H_to_Z ; enddo ! Depth of center of surface layer
642  call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pref_mld, rhosurf, tv%eqn_of_state, eosdom)
643  do i=is,ie
644  deltarhoatk(i) = 0.
645  mld(i,j) = 0.
646  if (id_n2>0) then
647  submln2(i,j) = 0.0
648  h_subml(i) = h(i,j,1) ; dh_n2(i) = 0.0
649  t_subml(i) = 0.0 ; s_subml(i) = 0.0 ; t_deeper(i) = 0.0 ; s_deeper(i) = 0.0
650  n2_region_set(i) = (g%mask2dT(i,j)<0.5) ! Only need to work on ocean points.
651  endif
652  enddo
653  do k=2,nz
654  do i=is,ie
655  dkm1(i) = dk(i) ! Depth of center of layer K-1
656  dk(i) = dk(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * gv%H_to_Z ! Depth of center of layer K
657  enddo
658 
659  ! Prepare to calculate stratification, N2, immediately below the mixed layer by finding
660  ! the cells that extend over at least dz_subML.
661  if (id_n2>0) then
662  do i=is,ie
663  if (mld(i,j)==0.0) then ! Still in the mixed layer.
664  h_subml(i) = h_subml(i) + h(i,j,k)
665  elseif (.not.n2_region_set(i)) then ! This block is below the mixed layer, but N2 has not been found yet.
666  if (dh_n2(i)==0.0) then ! Record the temperature, salinity, pressure, immediately below the ML
667  t_subml(i) = tv%T(i,j,k) ; s_subml(i) = tv%S(i,j,k)
668  h_subml(i) = h_subml(i) + 0.5 * h(i,j,k) ! Start midway through this layer.
669  dh_n2(i) = 0.5 * h(i,j,k)
670  elseif (dh_n2(i) + h(i,j,k) < dh_subml) then
671  dh_n2(i) = dh_n2(i) + h(i,j,k)
672  else ! This layer includes the base of the region where N2 is calculated.
673  t_deeper(i) = tv%T(i,j,k) ; s_deeper(i) = tv%S(i,j,k)
674  dh_n2(i) = dh_n2(i) + 0.5 * h(i,j,k)
675  n2_region_set(i) = .true.
676  endif
677  endif
678  enddo ! i-loop
679  endif ! id_N2>0
680 
681  ! Mixed-layer depth, using sigma-0 (surface reference pressure)
682  do i=is,ie ; deltarhoatkm1(i) = deltarhoatk(i) ; enddo ! Store value from previous iteration of K
683  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_mld, deltarhoatk, tv%eqn_of_state, eosdom)
684  do i = is, ie
685  deltarhoatk(i) = deltarhoatk(i) - rhosurf(i) ! Density difference between layer K and surface
686  ddrho = deltarhoatk(i) - deltarhoatkm1(i)
687  if ((mld(i,j)==0.) .and. (ddrho>0.) .and. &
688  (deltarhoatkm1(i)<densitydiff) .and. (deltarhoatk(i)>=densitydiff)) then
689  afac = ( densitydiff - deltarhoatkm1(i) ) / ddrho
690  mld(i,j) = dk(i) * afac + dkm1(i) * (1. - afac)
691  endif
692  if (id_sq > 0) mld2(i,j) = mld(i,j)**2
693  enddo ! i-loop
694  enddo ! k-loop
695  do i=is,ie
696  if ((mld(i,j)==0.) .and. (deltarhoatk(i)<densitydiff)) mld(i,j) = dk(i) ! Assume mixing to the bottom
697  enddo
698 
699  if (id_n2>0) then ! Now actually calculate stratification, N2, below the mixed layer.
700  do i=is,ie ; pref_n2(i) = (gv%g_Earth * gv%H_to_RZ) * (h_subml(i) + 0.5*dh_n2(i)) ; enddo
701  ! if ((.not.N2_region_set(i)) .and. (dH_N2(i) > 0.5*dH_subML)) then
702  ! ! Use whatever stratification we can, measured over whatever distance is available?
703  ! T_deeper(i) = tv%T(i,j,nz) ; S_deeper(i) = tv%S(i,j,nz)
704  ! N2_region_set(i) = .true.
705  ! endif
706  call calculate_density(t_subml, s_subml, pref_n2, rho_subml, tv%eqn_of_state, eosdom)
707  call calculate_density(t_deeper, s_deeper, pref_n2, rho_deeper, tv%eqn_of_state, eosdom)
708  do i=is,ie ; if ((g%mask2dT(i,j)>0.5) .and. n2_region_set(i)) then
709  submln2(i,j) = ge_rho0 * (rho_deeper(i) - rho_subml(i)) / (gv%H_to_z * dh_n2(i))
710  endif ; enddo
711  endif
712  enddo ! j-loop
713 
714  if (id_mld > 0) call post_data(id_mld, mld, diagptr)
715  if (id_n2 > 0) call post_data(id_n2, submln2, diagptr)
716  if (id_sq > 0) call post_data(id_sq, mld2, diagptr)
717 

◆ diagnosemldbyenergy()

subroutine, public mom_diabatic_aux::diagnosemldbyenergy ( integer, dimension(3), intent(in)  id_MLD,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension(3), intent(in)  Mixing_Energy,
type(diag_ctrl), pointer  diagPtr 
)

Diagnose a mixed layer depth (MLD) determined by the depth a given energy value would mix. This routine is appropriate in MOM_diabatic_driver due to its position within the time stepping.

Parameters
[in]id_mldEnergy output diag IDs
[in]gGrid type
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in]mixing_energyEnergy values for up to 3 MLDs
[in]hLayer thickness [H ~> m or kg m-2]
[in]tvStructure containing pointers to any available thermodynamic fields.
diagptrDiagnostics structure

Definition at line 722 of file MOM_diabatic_aux.F90.

723  ! Author: Brandon Reichl
724  ! Date: October 2, 2020
725  ! //
726  ! *Note that gravity is assumed constant everywhere and divided out of all calculations.
727  !
728  ! This code has been written to step through the columns layer by layer, summing the PE
729  ! change inferred by mixing the layer with all layers above. When the change exceeds a
730  ! threshold (determined by input array Mixing_Energy), the code needs to solve for how far
731  ! into this layer the threshold PE change occurs (assuming constant density layers).
732  ! This is expressed here via solving the function F(X) = 0 where:
733  ! F(X) = 0.5 * ( Ca*X^3/(D1+X) + Cb*X^2/(D1+X) + Cc*X/(D1+X) + Dc/(D1+X)
734  ! + Ca2*X^2 + Cb2*X + Cc2)
735  ! where all coefficients are determined by the previous mixed layer depth, the
736  ! density of the previous mixed layer, the present layer thickness, and the present
737  ! layer density. This equation is worked out by computing the total PE assuming constant
738  ! density in the mixed layer as well as in the remaining part of the present layer that is
739  ! not mixed.
740  ! To solve for X in this equation a Newton's method iteration is employed, which
741  ! converges extremely quickly (usually 1 guess) since this equation turns out being rather
742  ! lienar for PE change with increasing X.
743  ! Input parameters:
744  integer, dimension(3), intent(in) :: id_MLD !< Energy output diag IDs
745  type(ocean_grid_type), intent(in) :: G !< Grid type
746  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
747  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
748  real, dimension(3), intent(in) :: Mixing_Energy !< Energy values for up to 3 MLDs
749  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
750  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
751  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
752  !! available thermodynamic fields.
753  type(diag_ctrl), pointer :: diagPtr !< Diagnostics structure
754 
755  ! Local variables
756  real, dimension(SZI_(G), SZJ_(G),3) :: MLD ! Diagnosed mixed layer depth [Z ~> m].
757  real, dimension(SZK_(G)) :: Z_L, Z_U, dZ, Rho_c, pRef_MLD
758  real, dimension(3) :: PE_threshold
759 
760  real :: ig, E_g
761  real :: PE_Threshold_fraction, PE, PE_Mixed, PE_Mixed_TST
762  real :: RhoDZ_ML, H_ML, RhoDZ_ML_TST, H_ML_TST
763  real :: Rho_ML
764 
765  real :: R1, D1, R2, D2
766  real :: Ca, Cb,D ,Cc, Cd, Ca2, Cb2, C, Cc2
767  real :: Gx, Gpx, Hx, iHx, Hpx, Ix, Ipx, Fgx, Fpx, X, X2
768 
769  integer :: IT, iM
770  integer :: i, j, is, ie, js, je, k, nz
771 
772  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
773 
774  pref_mld(:) = 0.0
775  mld(:,:,:) = 0.0
776  pe_threshold_fraction = 1.e-4 !Fixed threshold of 0.01%, could be runtime.
777 
778  do im=1,3
779  pe_threshold(im) = mixing_energy(im)/gv%g_earth
780  enddo
781 
782  do j=js,je; do i=is,ie
783  if (g%mask2dT(i,j) > 0.0) then
784 
785  call calculate_density(tv%T(i,j,:), tv%S(i,j,:), pref_mld, rho_c, 1, nz, &
786  tv%eqn_of_state, scale=us%kg_m3_to_R)
787 
788  do k=1,nz
789  dz(k) = h(i,j,k) * gv%H_to_Z
790  enddo
791  z_u(1) = 0.0
792  z_l(1) = -dz(1)
793  do k=2,nz
794  z_u(k) = z_l(k-1)
795  z_l(k) = z_l(k-1)-dz(k)
796  enddo
797 
798  do im=1,3
799 
800  ! Initialize these for each column-wise calculation
801  pe = 0.0
802  rhodz_ml = 0.0
803  h_ml = 0.0
804  rhodz_ml_tst = 0.0
805  h_ml_tst = 0.0
806  pe_mixed = 0.0
807 
808  do k=1,nz
809 
810  ! This is the unmixed PE cumulative sum from top down
811  pe = pe + 0.5 * rho_c(k) * (z_u(k)**2 - z_l(k)**2)
812 
813  ! This is the depth and integral of density
814  h_ml_tst = h_ml + dz(k)
815  rhodz_ml_tst = rhodz_ml + rho_c(k) * dz(k)
816 
817  ! The average density assuming all layers including this were mixed
818  rho_ml = rhodz_ml_tst/h_ml_tst
819 
820  ! The PE assuming all layers including this were mixed
821  ! Note that 0. could be replaced with "Surface", which doesn't have to be 0
822  ! but 0 is a good reference value.
823  pe_mixed_tst = 0.5 * rho_ml * (0.**2 - (0. - h_ml_tst)**2)
824 
825  ! Check if we supplied enough energy to mix to this layer
826  if (pe_mixed_tst - pe <= pe_threshold(im)) then
827  h_ml = h_ml_tst
828  rhodz_ml = rhodz_ml_tst
829 
830  else ! If not, we need to solve where the energy ran out
831  ! This will be done with a Newton's method iteration:
832 
833  r1 = rhodz_ml / h_ml ! The density of the mixed layer (not including this layer)
834  d1 = h_ml ! The thickness of the mixed layer (not including this layer)
835  r2 = rho_c(k) ! The density of this layer
836  d2 = dz(k) ! The thickness of this layer
837 
838  ! This block could be used to calculate the function coefficients if
839  ! we don't reference all values to a surface designated as z=0
840  ! S = Surface
841  ! Ca = -(R2)
842  ! Cb = -( (R1*D1) + R2*(2.*D1-2.*S) )
843  ! D = D1**2. - 2.*D1*S
844  ! Cc = -( R1*D1*(2.*D1-2.*S) + R2*D )
845  ! Cd = -(R1*D1*D)
846  ! Ca2 = R2
847  ! Cb2 = R2*(2*D1-2*S)
848  ! C = S**2 + D2**2 + D1**2 - 2*D1*S - 2.*D2*S +2.*D1*D2
849  ! Cc2 = R2*(D+S**2-C)
850  !
851  ! If the surface is S = 0, it simplifies to:
852  ca = -r2
853  cb = -(r1 * d1 + r2 * (2. * d1))
854  d = d1**2
855  cc = -(r1 * d1 * (2. * d1) + (r2 * d))
856  cd = -r1 * (d1 * d)
857  ca2 = r2
858  cb2 = r2 * (2. * d1)
859  c = d2**2 + d1**2 + 2. * (d1 * d2)
860  cc2 = r2 * (d - c)
861 
862  ! First guess for an iteration using Newton's method
863  x = dz(k) * 0.5
864 
865  it=0
866  do while(it<10)!We can iterate up to 10 times
867  ! We are trying to solve the function:
868  ! F(x) = G(x)/H(x)+I(x)
869  ! for where F(x) = PE+PE_threshold, or equivalently for where
870  ! F(x) = G(x)/H(x)+I(x) - (PE+PE_threshold) = 0
871  ! We also need the derivative of this function for the Newton's method iteration
872  ! F'(x) = (G'(x)H(x)-G(x)H'(x))/H(x)^2 + I'(x)
873  ! G and its derivative
874  gx = 0.5 * (ca * (x*x*x) + cb * x**2 + cc * x + cd)
875  gpx = 0.5 * (3. * (ca * x**2) + 2. * (cb * x) + cc)
876  ! H, its inverse, and its derivative
877  hx = d1 + x
878  ihx = 1. / hx
879  hpx = 1.
880  ! I and its derivative
881  ix = 0.5 * (ca2 * x**2 + cb2 * x + cc2)
882  ipx = 0.5 * (2. * ca2 * x + cb2)
883 
884  ! The Function and its derivative:
885  pe_mixed = gx * ihx + ix
886  fgx = pe_mixed - (pe + pe_threshold(im))
887  fpx = (gpx * hx - hpx * gx) * ihx**2 + ipx
888 
889  ! Check if our solution is within the threshold bounds, if not update
890  ! using Newton's method. This appears to converge almost always in
891  ! one step because the function is very close to linear in most applications.
892  if (abs(fgx) > pe_threshold(im) * pe_threshold_fraction) then
893  x2 = x - fgx / fpx
894  it = it + 1
895  if (x2 < 0. .or. x2 > dz(k)) then
896  ! The iteration seems to be robust, but we need to do something *if*
897  ! things go wrong... How should we treat failed iteration?
898  ! Present solution: Stop trying to compute and just say we can't mix this layer.
899  x=0
900  exit
901  else
902  x = x2
903  endif
904  else
905  exit! Quit the iteration
906  endif
907  enddo
908  h_ml = h_ml + x
909  exit! Quit looping through the column
910  endif
911  enddo
912  mld(i,j,im) = h_ml
913  enddo
914  else
915  mld(i,j,:) = 0.0
916  endif
917  enddo ; enddo
918 
919  if (id_mld(1) > 0) call post_data(id_mld(1), mld(:,:,1), diagptr)
920  if (id_mld(2) > 0) call post_data(id_mld(2), mld(:,:,2), diagptr)
921  if (id_mld(3) > 0) call post_data(id_mld(3), mld(:,:,3), diagptr)
922 

◆ differential_diffuse_t_s()

subroutine, public mom_diabatic_aux::differential_diffuse_t_s ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  T,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  S,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(inout)  Kd_T,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(in)  Kd_S,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV 
)

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

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]hLayer thicknesses [H ~> m or kg m-2]
[in,out]tPotential temperature [degC].
[in,out]sSalinity [PSU] or [gSalt/kg], generically [ppt].
[in,out]kd_tThe extra diffusivity of temperature due to
[in]kd_sThe extra diffusivity of salinity due to
[in]dtTime increment [T ~> s].

Definition at line 226 of file MOM_diabatic_aux.F90.

227  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
228  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
229  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
230  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
231  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
232  intent(inout) :: T !< Potential temperature [degC].
233  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
234  intent(inout) :: S !< Salinity [PSU] or [gSalt/kg], generically [ppt].
235  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
236  intent(inout) :: Kd_T !< The extra diffusivity of temperature due to
237  !! double diffusion relative to the diffusivity of
238  !! diffusivity of density [Z2 T-1 ~> m2 s-1].
239  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
240  intent(in) :: Kd_S !< The extra diffusivity of salinity due to
241  !! double diffusion relative to the diffusivity of
242  !! diffusivity of density [Z2 T-1 ~> m2 s-1].
243  real, intent(in) :: dt !< Time increment [T ~> s].
244 
245  ! local variables
246  real, dimension(SZI_(G)) :: &
247  b1_T, b1_S, & ! Variables used by the tridiagonal solvers of T & S [H ~> m or kg m-2].
248  d1_T, d1_S ! Variables used by the tridiagonal solvers [nondim].
249  real, dimension(SZI_(G),SZK_(G)) :: &
250  c1_T, c1_S ! Variables used by the tridiagonal solvers [H ~> m or kg m-2].
251  real, dimension(SZI_(G),SZK_(G)+1) :: &
252  mix_T, mix_S ! Mixing distances in both directions across each interface [H ~> m or kg m-2].
253  real :: h_tr ! h_tr is h at tracer points with a tiny thickness
254  ! added to ensure positive definiteness [H ~> m or kg m-2].
255  real :: h_neglect ! A thickness that is so small it is usually lost
256  ! in roundoff and can be neglected [H ~> m or kg m-2].
257  real :: I_h_int ! The inverse of the thickness associated with an
258  ! interface [H-1 ~> m-1 or m2 kg-1].
259  real :: b_denom_T ! The first term in the denominators for the expressions
260  real :: b_denom_S ! for b1_T and b1_S, both [H ~> m or kg m-2].
261  integer :: i, j, k, is, ie, js, je, nz
262 
263  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
264  h_neglect = gv%H_subroundoff
265 
266  !$OMP parallel do default(private) shared(is,ie,js,je,h,h_neglect,dt,Kd_T,Kd_S,G,GV,T,S,nz)
267  do j=js,je
268  do i=is,ie
269  i_h_int = 1.0 / (0.5 * (h(i,j,1) + h(i,j,2)) + h_neglect)
270  mix_t(i,2) = ((dt * kd_t(i,j,2)) * gv%Z_to_H**2) * i_h_int
271  mix_s(i,2) = ((dt * kd_s(i,j,2)) * gv%Z_to_H**2) * i_h_int
272 
273  h_tr = h(i,j,1) + h_neglect
274  b1_t(i) = 1.0 / (h_tr + mix_t(i,2))
275  b1_s(i) = 1.0 / (h_tr + mix_s(i,2))
276  d1_t(i) = h_tr * b1_t(i)
277  d1_s(i) = h_tr * b1_s(i)
278  t(i,j,1) = (b1_t(i)*h_tr)*t(i,j,1)
279  s(i,j,1) = (b1_s(i)*h_tr)*s(i,j,1)
280  enddo
281  do k=2,nz-1 ; do i=is,ie
282  ! Calculate the mixing across the interface below this layer.
283  i_h_int = 1.0 / (0.5 * (h(i,j,k) + h(i,j,k+1)) + h_neglect)
284  mix_t(i,k+1) = ((dt * kd_t(i,j,k+1)) * gv%Z_to_H**2) * i_h_int
285  mix_s(i,k+1) = ((dt * kd_s(i,j,k+1)) * gv%Z_to_H**2) * i_h_int
286 
287  c1_t(i,k) = mix_t(i,k) * b1_t(i)
288  c1_s(i,k) = mix_s(i,k) * b1_s(i)
289 
290  h_tr = h(i,j,k) + h_neglect
291  b_denom_t = h_tr + d1_t(i)*mix_t(i,k)
292  b_denom_s = h_tr + d1_s(i)*mix_s(i,k)
293  b1_t(i) = 1.0 / (b_denom_t + mix_t(i,k+1))
294  b1_s(i) = 1.0 / (b_denom_s + mix_s(i,k+1))
295  d1_t(i) = b_denom_t * b1_t(i)
296  d1_s(i) = b_denom_s * b1_s(i)
297 
298  t(i,j,k) = b1_t(i) * (h_tr*t(i,j,k) + mix_t(i,k)*t(i,j,k-1))
299  s(i,j,k) = b1_s(i) * (h_tr*s(i,j,k) + mix_s(i,k)*s(i,j,k-1))
300  enddo ; enddo
301  do i=is,ie
302  c1_t(i,nz) = mix_t(i,nz) * b1_t(i)
303  c1_s(i,nz) = mix_s(i,nz) * b1_s(i)
304 
305  h_tr = h(i,j,nz) + h_neglect
306  b1_t(i) = 1.0 / (h_tr + d1_t(i)*mix_t(i,nz))
307  b1_s(i) = 1.0 / (h_tr + d1_s(i)*mix_s(i,nz))
308 
309  t(i,j,nz) = b1_t(i) * (h_tr*t(i,j,nz) + mix_t(i,nz)*t(i,j,nz-1))
310  s(i,j,nz) = b1_s(i) * (h_tr*s(i,j,nz) + mix_s(i,nz)*s(i,j,nz-1))
311  enddo
312  do k=nz-1,1,-1 ; do i=is,ie
313  t(i,j,k) = t(i,j,k) + c1_t(i,k+1)*t(i,j,k+1)
314  s(i,j,k) = s(i,j,k) + c1_s(i,k+1)*s(i,j,k+1)
315  enddo ; enddo
316  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,
type(unit_scale_type), intent(in)  US,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in), optional  ea,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in), optional  eb 
)

This subroutine calculates u_h and v_h (velocities at thickness points), optionally using the entrainment amounts passed in as arguments.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in]uThe zonal velocity [L T-1 ~> m s-1]
[in]vThe meridional velocity [L T-1 ~> m s-1]
[in]hLayer thicknesses [H ~> m or kg m-2]
[out]u_hZonal velocity interpolated to h points [L T-1 ~> m s-1].
[out]v_hMeridional velocity interpolated to h points [L T-1 ~> m s-1].
[in]eaThe amount of fluid entrained from the layer
[in]ebThe amount of fluid entrained from the layer

Definition at line 430 of file MOM_diabatic_aux.F90.

431  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
432  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
433  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
434  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
435  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
436  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
437  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
438  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
439  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
440  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
441  intent(out) :: u_h !< Zonal velocity interpolated to h points [L T-1 ~> m s-1].
442  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
443  intent(out) :: v_h !< Meridional velocity interpolated to h points [L T-1 ~> m s-1].
444  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
445  optional, intent(in) :: ea !< The amount of fluid entrained from the layer
446  !! above within this time step [H ~> m or kg m-2].
447  !! Omitting ea is the same as setting it to 0.
448  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
449  optional, intent(in) :: eb !< The amount of fluid entrained from the layer
450  !! below within this time step [H ~> m or kg m-2].
451  !! Omitting eb is the same as setting it to 0.
452 
453  ! local variables
454  real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
455  real :: h_neglect ! A thickness that is so small it is usually lost
456  ! in roundoff and can be neglected [H ~> m or kg m-2].
457  real :: b1(SZI_(G)), d1(SZI_(G)), c1(SZI_(G),SZK_(G))
458  real :: a_n(SZI_(G)), a_s(SZI_(G)) ! Fractional weights of the neighboring
459  real :: a_e(SZI_(G)), a_w(SZI_(G)) ! velocity points, ~1/2 in the open
460  ! ocean, nondimensional.
461  real :: sum_area, Idenom
462  logical :: mix_vertically
463  integer :: i, j, k, is, ie, js, je, nz
464  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
465  call cpu_clock_begin(id_clock_uv_at_h)
466  h_neglect = gv%H_subroundoff
467 
468  mix_vertically = present(ea)
469  if (present(ea) .neqv. present(eb)) call mom_error(fatal, &
470  "find_uv_at_h: Either both ea and eb or neither one must be present "// &
471  "in call to find_uv_at_h.")
472 !$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,h,h_neglect, &
473 !$OMP eb,u_h,u,v_h,v,nz,ea) &
474 !$OMP private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1)
475  do j=js,je
476  do i=is,ie
477  sum_area = g%areaCu(i-1,j) + g%areaCu(i,j)
478  if (sum_area>0.0) then
479  idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
480  a_w(i) = g%areaCu(i-1,j) * idenom
481  a_e(i) = g%areaCu(i,j) * idenom
482  else
483  a_w(i) = 0.0 ; a_e(i) = 0.0
484  endif
485 
486  sum_area = g%areaCv(i,j-1) + g%areaCv(i,j)
487  if (sum_area>0.0) then
488  idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
489  a_s(i) = g%areaCv(i,j-1) * idenom
490  a_n(i) = g%areaCv(i,j) * idenom
491  else
492  a_s(i) = 0.0 ; a_n(i) = 0.0
493  endif
494  enddo
495 
496  if (mix_vertically) then
497  do i=is,ie
498  b_denom_1 = h(i,j,1) + h_neglect
499  b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
500  d1(i) = b_denom_1 * b1(i)
501  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))
502  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))
503  enddo
504  do k=2,nz ; do i=is,ie
505  c1(i,k) = eb(i,j,k-1) * b1(i)
506  b_denom_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h_neglect
507  b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
508  d1(i) = b_denom_1 * b1(i)
509  u_h(i,j,k) = (h(i,j,k) * (a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)) + &
510  ea(i,j,k)*u_h(i,j,k-1))*b1(i)
511  v_h(i,j,k) = (h(i,j,k) * (a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)) + &
512  ea(i,j,k)*v_h(i,j,k-1))*b1(i)
513  enddo ; enddo
514  do k=nz-1,1,-1 ; do i=is,ie
515  u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1)
516  v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1)
517  enddo ; enddo
518  else
519  do k=1,nz ; do i=is,ie
520  u_h(i,j,k) = a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)
521  v_h(i,j,k) = a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)
522  enddo ; enddo
523  endif
524  enddo
525 
526  call cpu_clock_end(id_clock_uv_at_h)

◆ make_frazil()

subroutine, public mom_diabatic_aux::make_frazil ( 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(unit_scale_type), intent(in)  US,
type(diabatic_aux_cs), intent(in)  CS,
real, dimension(szi_(g),szj_(g)), intent(in), optional  p_surf,
integer, intent(in), optional  halo 
)

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 [Q R Z ~> J m-2]) in tvfrazil.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]hLayer thicknesses [H ~> m or kg m-2]
[in,out]tvStructure containing pointers to any available thermodynamic fields.
[in]usA dimensional unit scaling type
[in]csThe control structure returned by a previous call to diabatic_aux_init.
[in]p_surfThe pressure at the ocean surface [R L2 T-2 ~> Pa].
[in]haloHalo width over which to calculate frazil

Definition at line 103 of file MOM_diabatic_aux.F90.

104  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
105  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
106  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
107  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
108  type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any available
109  !! thermodynamic fields.
110  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
111  type(diabatic_aux_CS), intent(in) :: CS !< The control structure returned by a previous
112  !! call to diabatic_aux_init.
113  real, dimension(SZI_(G),SZJ_(G)), &
114  optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa].
115  integer, optional, intent(in) :: halo !< Halo width over which to calculate frazil
116 
117  ! Local variables
118  real, dimension(SZI_(G)) :: &
119  fraz_col, & ! The accumulated heat requirement due to frazil [Q R Z ~> J m-2].
120  T_freeze, & ! The freezing potential temperature at the current salinity [degC].
121  ps ! Surface pressure [R L2 T-2 ~> Pa]
122  real, dimension(SZI_(G),SZK_(G)) :: &
123  pressure ! The pressure at the middle of each layer [R L2 T-2 ~> Pa].
124  real :: H_to_RL2_T2 ! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
125  real :: hc ! A layer's heat capacity [Q R Z degC-1 ~> J m-2 degC-1].
126  logical :: T_fr_set ! True if the freezing point has been calculated for a
127  ! row of points.
128  integer :: i, j, k, is, ie, js, je, nz
129 
130  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
131  if (present(halo)) then
132  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
133  endif
134 
135  call cpu_clock_begin(id_clock_frazil)
136 
137  if (.not.cs%pressure_dependent_frazil) then
138  do k=1,nz ; do i=is,ie ; pressure(i,k) = 0.0 ; enddo ; enddo
139  else
140  h_to_rl2_t2 = gv%H_to_RZ * gv%g_Earth
141  endif
142  !$OMP parallel do default(shared) private(fraz_col,T_fr_set,T_freeze,hc,ps) &
143  !$OMP firstprivate(pressure) ! pressure might be set above, so should be firstprivate
144  do j=js,je
145  ps(:) = 0.0
146  if (PRESENT(p_surf)) then ; do i=is,ie
147  ps(i) = p_surf(i,j)
148  enddo ; endif
149 
150  do i=is,ie ; fraz_col(i) = 0.0 ; enddo
151 
152  if (cs%pressure_dependent_frazil) then
153  do i=is,ie
154  pressure(i,1) = ps(i) + (0.5*h_to_rl2_t2)*h(i,j,1)
155  enddo
156  do k=2,nz ; do i=is,ie
157  pressure(i,k) = pressure(i,k-1) + &
158  (0.5*h_to_rl2_t2) * (h(i,j,k) + h(i,j,k-1))
159  enddo ; enddo
160  endif
161 
162  if (cs%reclaim_frazil) then
163  t_fr_set = .false.
164  do i=is,ie ; if (tv%frazil(i,j) > 0.0) then
165  if (.not.t_fr_set) then
166  call calculate_tfreeze(tv%S(i:,j,1), pressure(i:,1), t_freeze(i:), &
167  1, ie-i+1, tv%eqn_of_state, pres_scale=us%RL2_T2_to_Pa)
168  t_fr_set = .true.
169  endif
170 
171  if (tv%T(i,j,1) > t_freeze(i)) then
172  ! If frazil had previously been formed, but the surface temperature is now
173  ! above freezing, cool the surface layer with the frazil heat deficit.
174  hc = (tv%C_p*gv%H_to_RZ) * h(i,j,1)
175  if (tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i)) <= 0.0) then
176  tv%T(i,j,1) = tv%T(i,j,1) - tv%frazil(i,j) / hc
177  tv%frazil(i,j) = 0.0
178  else
179  tv%frazil(i,j) = tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i))
180  tv%T(i,j,1) = t_freeze(i)
181  endif
182  endif
183  endif ; enddo
184  endif
185 
186  do k=nz,1,-1
187  t_fr_set = .false.
188  do i=is,ie
189  if ((g%mask2dT(i,j) > 0.0) .and. &
190  ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then
191  if (.not.t_fr_set) then
192  call calculate_tfreeze(tv%S(i:,j,k), pressure(i:,k), t_freeze(i:), &
193  1, ie-i+1, tv%eqn_of_state, pres_scale=us%RL2_T2_to_Pa)
194  t_fr_set = .true.
195  endif
196 
197  hc = (tv%C_p*gv%H_to_RZ) * h(i,j,k)
198  if (h(i,j,k) <= 10.0*(gv%Angstrom_H + gv%H_subroundoff)) then
199  ! Very thin layers should not be cooled by the frazil flux.
200  if (tv%T(i,j,k) < t_freeze(i)) then
201  fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
202  tv%T(i,j,k) = t_freeze(i)
203  endif
204  elseif ((fraz_col(i) > 0.0) .or. (tv%T(i,j,k) < t_freeze(i))) then
205  if (fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k)) < 0.0) then
206  tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i) / hc
207  fraz_col(i) = 0.0
208  else
209  fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
210  tv%T(i,j,k) = t_freeze(i)
211  endif
212  endif
213  endif
214  enddo
215  enddo
216  do i=is,ie
217  tv%frazil(i,j) = tv%frazil(i,j) + fraz_col(i)
218  enddo
219  enddo
220  call cpu_clock_end(id_clock_frazil)
221 

◆ set_pen_shortwave()

subroutine, public mom_diabatic_aux::set_pen_shortwave ( type(optics_type), pointer  optics,
type(forcing), intent(inout)  fluxes,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diabatic_aux_cs), pointer  CS,
type(opacity_cs), pointer  opacity_CSp,
type(tracer_flow_control_cs), pointer  tracer_flow_CSp 
)
Parameters
opticsAn optics structure that has will contain information about shortwave fluxes and absorption.
[in,out]fluxespoints to forcing fields unused fields have NULL ptrs
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
csControl structure for diabatic_aux
opacity_cspThe control structure for the opacity module.
tracer_flow_cspA pointer to the control structure organizing the tracer modules.

Definition at line 530 of file MOM_diabatic_aux.F90.

531  type(optics_type), pointer :: optics !< An optics structure that has will contain
532  !! information about shortwave fluxes and absorption.
533  type(forcing), intent(inout) :: fluxes !< points to forcing fields
534  !! unused fields have NULL ptrs
535  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
536  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
537  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
538  type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
539  type(opacity_CS), pointer :: opacity_CSp !< The control structure for the opacity module.
540  type(tracer_flow_control_CS), pointer :: tracer_flow_CSp !< A pointer to the control structure
541  !! organizing the tracer modules.
542 
543  ! Local variables
544  real, dimension(SZI_(G),SZJ_(G)) :: chl_2d !< Vertically uniform chlorophyll-A concentractions [mg m-3]
545  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: chl_3d !< The chlorophyll-A concentractions of each layer [mg m-3]
546  character(len=128) :: mesg
547  integer :: i, j, k, is, ie, js, je
548  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
549 
550  if (.not.associated(optics)) return
551 
552  if (cs%var_pen_sw) then
553  if (cs%chl_from_file) then
554  ! Only the 2-d surface chlorophyll can be read in from a file. The
555  ! same value is assumed for all layers.
556  call time_interp_external(cs%sbc_chl, cs%Time, chl_2d)
557  do j=js,je ; do i=is,ie
558  if ((g%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0)) then
559  write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&
560  & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
561  chl_2d(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j)
562  call mom_error(fatal, "MOM_diabatic_aux set_pen_shortwave: "//trim(mesg))
563  endif
564  enddo ; enddo
565 
566  if (cs%id_chl > 0) call post_data(cs%id_chl, chl_2d, cs%diag)
567 
568  call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
569  fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, us, opacity_csp, chl_2d=chl_2d)
570  else
571  if (.not.associated(tracer_flow_csp)) call mom_error(fatal, &
572  "The tracer flow control structure must be associated when the model sets "//&
573  "the chlorophyll internally in set_pen_shortwave.")
574  call get_chl_from_model(chl_3d, g, tracer_flow_csp)
575 
576  if (cs%id_chl > 0) call post_data(cs%id_chl, chl_3d(:,:,1), cs%diag)
577 
578  call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
579  fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, us, opacity_csp, chl_3d=chl_3d)
580  endif
581  else
582  call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
583  fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, us, opacity_csp)
584  endif
585 

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

This is a simple tri-diagonal solver for T and S. "Simple" means it only uses arrays hold, ea and eb.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]isThe start i-index to work on.
[in]ieThe end i-index to work on.
[in]jsThe start j-index to work on.
[in]jeThe end j-index to work on.
[in]holdThe layer thicknesses before entrainment, [H ~> m or kg m-2].
[in]eaThe amount of fluid entrained from the layer above within this time step [H ~> m or kg m-2]
[in]ebThe amount of fluid entrained from the layer below within this time step [H ~> m or kg m-2]
[in,out]tLayer potential temperatures [degC].
[in,out]sLayer salinities [ppt].

Definition at line 381 of file MOM_diabatic_aux.F90.

382  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
383  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
384  integer, intent(in) :: is !< The start i-index to work on.
385  integer, intent(in) :: ie !< The end i-index to work on.
386  integer, intent(in) :: js !< The start j-index to work on.
387  integer, intent(in) :: je !< The end j-index to work on.
388  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hold !< The layer thicknesses before entrainment,
389  !! [H ~> m or kg m-2].
390  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: ea !< The amount of fluid entrained from the layer
391  !! above within this time step [H ~> m or kg m-2]
392  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: eb !< The amount of fluid entrained from the layer
393  !! below within this time step [H ~> m or kg m-2]
394  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: T !< Layer potential temperatures [degC].
395  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: S !< Layer salinities [ppt].
396 
397  ! Local variables
398  real :: b1(SZIB_(G)), d1(SZIB_(G)) ! b1, c1, and d1 are variables used by the
399  real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver.
400  real :: h_tr, b_denom_1
401  integer :: i, j, k
402 
403  !$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1)
404  do j=js,je
405  do i=is,ie
406  h_tr = hold(i,j,1) + gv%H_subroundoff
407  b1(i) = 1.0 / (h_tr + eb(i,j,1))
408  d1(i) = h_tr * b1(i)
409  t(i,j,1) = (b1(i)*h_tr)*t(i,j,1)
410  s(i,j,1) = (b1(i)*h_tr)*s(i,j,1)
411  enddo
412  do k=2,g%ke ; do i=is,ie
413  c1(i,k) = eb(i,j,k-1) * b1(i)
414  h_tr = hold(i,j,k) + gv%H_subroundoff
415  b_denom_1 = h_tr + d1(i)*ea(i,j,k)
416  b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
417  d1(i) = b_denom_1 * b1(i)
418  t(i,j,k) = b1(i) * (h_tr*t(i,j,k) + ea(i,j,k)*t(i,j,k-1))
419  s(i,j,k) = b1(i) * (h_tr*s(i,j,k) + ea(i,j,k)*s(i,j,k-1))
420  enddo ; enddo
421  do k=g%ke-1,1,-1 ; do i=is,ie
422  t(i,j,k) = t(i,j,k) + c1(i,k+1)*t(i,j,k+1)
423  s(i,j,k) = s(i,j,k) + c1(i,k+1)*s(i,j,k+1)
424  enddo ; enddo
425  enddo