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, tv, visc, 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...
 

Variables

integer id_clock_uv_at_h
 CPU time clock IDs.
 
integer id_clock_frazil
 CPU time clock IDs.
 

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 330 of file MOM_diabatic_aux.F90.

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

◆ 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 939 of file MOM_diabatic_aux.F90.

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

◆ 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 1652 of file MOM_diabatic_aux.F90.

1652  type(diabatic_aux_cs), pointer :: cs !< The control structure returned by a previous
1653  !! call to diabatic_aux_init; it is deallocated here.
1654 
1655  if (.not.associated(cs)) return
1656 
1657  if (cs%id_createdH >0) deallocate(cs%createdH)
1658  if (cs%id_penSW_diag >0) deallocate(cs%penSW_diag)
1659  if (cs%id_penSWflux_diag >0) deallocate(cs%penSWflux_diag)
1660  if (cs%id_nonpenSW_diag >0) deallocate(cs%nonpenSW_diag)
1661 
1662  if (associated(cs)) deallocate(cs)
1663 

◆ 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 1488 of file MOM_diabatic_aux.F90.

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

◆ 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 600 of file MOM_diabatic_aux.F90.

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

◆ 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 730 of file MOM_diabatic_aux.F90.

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

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

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]tvStructure containing pointers to any available thermodynamic fields.
[in]viscStructure containing vertical viscosities, bottom boundary layer properies, and related fields.
[in]dtTime increment [T ~> s].

Definition at line 227 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  type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any
232  !! available thermodynamic fields.
233  type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom
234  !! boundary layer properies, and related fields.
235  real, intent(in) :: dt !< Time increment [T ~> s].
236 
237  ! local variables
238  real, dimension(SZI_(G)) :: &
239  b1_t, b1_s, & ! Variables used by the tridiagonal solvers of T & S [H ~> m or kg m-2].
240  d1_t, d1_s ! Variables used by the tridiagonal solvers [nondim].
241  real, dimension(SZI_(G),SZK_(G)) :: &
242  c1_t, c1_s ! Variables used by the tridiagonal solvers [H ~> m or kg m-2].
243  real, dimension(SZI_(G),SZK_(G)+1) :: &
244  mix_t, mix_s ! Mixing distances in both directions across each interface [H ~> m or kg m-2].
245  real :: h_tr ! h_tr is h at tracer points with a tiny thickness
246  ! added to ensure positive definiteness [H ~> m or kg m-2].
247  real :: h_neglect ! A thickness that is so small it is usually lost
248  ! in roundoff and can be neglected [H ~> m or kg m-2].
249  real :: i_h_int ! The inverse of the thickness associated with an
250  ! interface [H-1 ~> m-1 or m2 kg-1].
251  real :: b_denom_t ! The first term in the denominators for the expressions
252  real :: b_denom_s ! for b1_T and b1_S, both [H ~> m or kg m-2].
253  real, dimension(:,:,:), pointer :: t=>null(), s=>null()
254  real, dimension(:,:,:), pointer :: kd_t=>null(), kd_s=>null() ! Diffusivities [Z2 T-1 ~> m2 s-1].
255  integer :: i, j, k, is, ie, js, je, nz
256 
257  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
258  h_neglect = gv%H_subroundoff
259 
260  if (.not.associated(tv%T)) call mom_error(fatal, &
261  "differential_diffuse_T_S: Called with an unassociated tv%T")
262  if (.not.associated(tv%S)) call mom_error(fatal, &
263  "differential_diffuse_T_S: Called with an unassociated tv%S")
264  if (.not.associated(visc%Kd_extra_T)) call mom_error(fatal, &
265  "differential_diffuse_T_S: Called with an unassociated visc%Kd_extra_T")
266  if (.not.associated(visc%Kd_extra_S)) call mom_error(fatal, &
267  "differential_diffuse_T_S: Called with an unassociated visc%Kd_extra_S")
268 
269  t => tv%T ; s => tv%S
270  kd_t => visc%Kd_extra_T ; kd_s => visc%Kd_extra_S
271 !$OMP parallel do default(none) shared(is,ie,js,je,h,h_neglect,dt,Kd_T,Kd_S,G,GV,T,S,nz) &
272 !$OMP private(I_h_int,mix_T,mix_S,h_tr,b1_T,b1_S, &
273 !$OMP d1_T,d1_S,c1_T,c1_S,b_denom_T,b_denom_S)
274  do j=js,je
275  do i=is,ie
276  i_h_int = 1.0 / (0.5 * (h(i,j,1) + h(i,j,2)) + h_neglect)
277  mix_t(i,2) = ((dt * kd_t(i,j,2)) * gv%Z_to_H**2) * i_h_int
278  mix_s(i,2) = ((dt * kd_s(i,j,2)) * gv%Z_to_H**2) * i_h_int
279 
280  h_tr = h(i,j,1) + h_neglect
281  b1_t(i) = 1.0 / (h_tr + mix_t(i,2))
282  b1_s(i) = 1.0 / (h_tr + mix_s(i,2))
283  d1_t(i) = h_tr * b1_t(i)
284  d1_s(i) = h_tr * b1_s(i)
285  t(i,j,1) = (b1_t(i)*h_tr)*t(i,j,1)
286  s(i,j,1) = (b1_s(i)*h_tr)*s(i,j,1)
287  enddo
288  do k=2,nz-1 ; do i=is,ie
289  ! Calculate the mixing across the interface below this layer.
290  i_h_int = 1.0 / (0.5 * (h(i,j,k) + h(i,j,k+1)) + h_neglect)
291  mix_t(i,k+1) = ((dt * kd_t(i,j,k+1)) * gv%Z_to_H**2) * i_h_int
292  mix_s(i,k+1) = ((dt * kd_s(i,j,k+1)) * gv%Z_to_H**2) * i_h_int
293 
294  c1_t(i,k) = mix_t(i,k) * b1_t(i)
295  c1_s(i,k) = mix_s(i,k) * b1_s(i)
296 
297  h_tr = h(i,j,k) + h_neglect
298  b_denom_t = h_tr + d1_t(i)*mix_t(i,k)
299  b_denom_s = h_tr + d1_s(i)*mix_s(i,k)
300  b1_t(i) = 1.0 / (b_denom_t + mix_t(i,k+1))
301  b1_s(i) = 1.0 / (b_denom_s + mix_s(i,k+1))
302  d1_t(i) = b_denom_t * b1_t(i)
303  d1_s(i) = b_denom_s * b1_s(i)
304 
305  t(i,j,k) = b1_t(i) * (h_tr*t(i,j,k) + mix_t(i,k)*t(i,j,k-1))
306  s(i,j,k) = b1_s(i) * (h_tr*s(i,j,k) + mix_s(i,k)*s(i,j,k-1))
307  enddo ; enddo
308  do i=is,ie
309  c1_t(i,nz) = mix_t(i,nz) * b1_t(i)
310  c1_s(i,nz) = mix_s(i,nz) * b1_s(i)
311 
312  h_tr = h(i,j,nz) + h_neglect
313  b1_t(i) = 1.0 / (h_tr + d1_t(i)*mix_t(i,nz))
314  b1_s(i) = 1.0 / (h_tr + d1_s(i)*mix_s(i,nz))
315 
316  t(i,j,nz) = b1_t(i) * (h_tr*t(i,j,nz) + mix_t(i,nz)*t(i,j,nz-1))
317  s(i,j,nz) = b1_s(i) * (h_tr*s(i,j,nz) + mix_s(i,nz)*s(i,j,nz-1))
318  enddo
319  do k=nz-1,1,-1 ; do i=is,ie
320  t(i,j,k) = t(i,j,k) + c1_t(i,k+1)*t(i,j,k+1)
321  s(i,j,k) = s(i,j,k) + c1_s(i,k+1)*s(i,j,k+1)
322  enddo ; enddo
323  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 438 of file MOM_diabatic_aux.F90.

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

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

◆ 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 389 of file MOM_diabatic_aux.F90.

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