MOM6
Data Types
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...
 
integer id_clock_uv_at_h
 CPU time clock IDs. More...
 
integer id_clock_frazil
 CPU time clock IDs. More...
 
subroutine, public make_frazil (h, tv, G, GV, 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 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 insert_brine (h, tv, G, GV, US, fluxes, nkmb, CS, dt, id_brine_lay)
 Insert salt from brine rejection into the first layer below the mixed layer which both contains mass and in which the change in layer density remains stable after the addition of salt via brine rejection. 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, 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 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 326 of file MOM_diabatic_aux.F90.

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

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

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

References mom_opacity::absorbremainingsw(), mom_eos::calculate_specific_vol_derivs(), mom_forcing_type::extractfluxes1d(), and mom_forcing_type::forcing_singlepointprint().

Here is the call graph for this function:

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

1560  type(diabatic_aux_CS), pointer :: CS !< The control structure returned by a previous
1561  !! call to diabatic_aux_init; it is deallocated here.
1562 
1563  if (.not.associated(cs)) return
1564 
1565  if (cs%id_createdH >0) deallocate(cs%createdH)
1566  if (cs%id_penSW_diag >0) deallocate(cs%penSW_diag)
1567  if (cs%id_penSWflux_diag >0) deallocate(cs%penSWflux_diag)
1568  if (cs%id_nonpenSW_diag >0) deallocate(cs%nonpenSW_diag)
1569 
1570  if (associated(cs)) deallocate(cs)
1571 

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

1395  type(time_type), target, intent(in) :: Time !< The current model time.
1396  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1397  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1398  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1399  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1400  type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output
1401  type(diabatic_aux_CS), pointer :: CS !< A pointer to the control structure for the
1402  !! diabatic_aux module, which is initialized here.
1403  logical, intent(in) :: useALEalgorithm !< If true, use the ALE algorithm rather
1404  !! than layered mode.
1405  logical, intent(in) :: use_ePBL !< If true, use the implicit energetics planetary
1406  !! boundary layer scheme to determine the diffusivity
1407  !! in the surface boundary layer.
1408 
1409 ! This "include" declares and sets the variable "version".
1410 #include "version_variable.h"
1411  character(len=40) :: mdl = "MOM_diabatic_aux" ! This module's name.
1412  character(len=48) :: thickness_units
1413  character(len=200) :: inputdir ! The directory where NetCDF input files
1414  character(len=240) :: chl_filename ! A file from which chl_a concentrations are to be read.
1415  character(len=128) :: chl_file ! Data containing chl_a concentrations. Used
1416  ! when var_pen_sw is defined and reading from file.
1417  character(len=32) :: chl_varname ! Name of chl_a variable in chl_file.
1418  logical :: use_temperature ! True if thermodynamics are enabled.
1419  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nbands
1420  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1421  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1422 
1423  if (associated(cs)) then
1424  call mom_error(warning, "diabatic_aux_init called with an "// &
1425  "associated control structure.")
1426  return
1427  else
1428  allocate(cs)
1429  endif
1430 
1431  cs%diag => diag
1432  cs%Time => time
1433 
1434 ! Set default, read and log parameters
1435  call log_version(param_file, mdl, version, &
1436  "The following parameters are used for auxiliary diabatic processes.")
1437 
1438  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
1439  "If true, temperature and salinity are used as state "//&
1440  "variables.", default=.true.)
1441 
1442  call get_param(param_file, mdl, "RECLAIM_FRAZIL", cs%reclaim_frazil, &
1443  "If true, try to use any frazil heat deficit to cool any "//&
1444  "overlying layers down to the freezing point, thereby "//&
1445  "avoiding the creation of thin ice when the SST is above "//&
1446  "the freezing point.", default=.true.)
1447  call get_param(param_file, mdl, "PRESSURE_DEPENDENT_FRAZIL", &
1448  cs%pressure_dependent_frazil, &
1449  "If true, use a pressure dependent freezing temperature "//&
1450  "when making frazil. The default is false, which will be "//&
1451  "faster but is inappropriate with ice-shelf cavities.", &
1452  default=.false.)
1453 
1454  if (use_epbl) then
1455  call get_param(param_file, mdl, "IGNORE_FLUXES_OVER_LAND", cs%ignore_fluxes_over_land,&
1456  "If true, the model does not check if fluxes are being applied "//&
1457  "over land points. This is needed when the ocean is coupled "//&
1458  "with ice shelves and sea ice, since the sea ice mask needs to "//&
1459  "be different than the ocean mask to avoid sea ice formation "//&
1460  "under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
1461  call get_param(param_file, mdl, "DO_RIVERMIX", cs%do_rivermix, &
1462  "If true, apply additional mixing wherever there is "//&
1463  "runoff, so that it is mixed down to RIVERMIX_DEPTH "//&
1464  "if the ocean is that deep.", default=.false.)
1465  if (cs%do_rivermix) &
1466  call get_param(param_file, mdl, "RIVERMIX_DEPTH", cs%rivermix_depth, &
1467  "The depth to which rivers are mixed if DO_RIVERMIX is "//&
1468  "defined.", units="m", default=0.0, scale=us%m_to_Z)
1469  else
1470  cs%do_rivermix = .false. ; cs%rivermix_depth = 0.0 ; cs%ignore_fluxes_over_land = .false.
1471  endif
1472 
1473  if (gv%nkml == 0) then
1474  call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
1475  "If true, use the fluxes%runoff_Hflx field to set the "//&
1476  "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
1477  default=.false.)
1478  call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
1479  "If true, use the fluxes%calving_Hflx field to set the "//&
1480  "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
1481  default=.false.)
1482  else
1483  cs%use_river_heat_content = .false.
1484  cs%use_calving_heat_content = .false.
1485  endif
1486 
1487  if (usealealgorithm) then
1488  cs%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, &
1489  time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", &
1490  "m s-1", conversion=gv%H_to_m*us%s_to_T)
1491  if (cs%id_createdH>0) allocate(cs%createdH(isd:ied,jsd:jed))
1492 
1493  ! diagnostic for heating of a grid cell from convergence of SW heat into the cell
1494  cs%id_penSW_diag = register_diag_field('ocean_model', 'rsdoabsorb', &
1495  diag%axesTL, time, 'Convergence of Penetrative Shortwave Flux in Sea Water Layer',&
1496  'W m-2', standard_name='net_rate_of_absorption_of_shortwave_energy_in_ocean_layer',v_extensive=.true.)
1497 
1498  ! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces)
1499  ! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock).
1500  cs%id_penSWflux_diag = register_diag_field('ocean_model', 'rsdo', &
1501  diag%axesTi, time, 'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',&
1502  'W m-2', standard_name='downwelling_shortwave_flux_in_sea_water')
1503 
1504  ! need both arrays for the SW diagnostics (one for flux, one for convergence)
1505  if (cs%id_penSW_diag>0 .or. cs%id_penSWflux_diag>0) then
1506  allocate(cs%penSW_diag(isd:ied,jsd:jed,nz))
1507  cs%penSW_diag(:,:,:) = 0.0
1508  allocate(cs%penSWflux_diag(isd:ied,jsd:jed,nz+1))
1509  cs%penSWflux_diag(:,:,:) = 0.0
1510  endif
1511 
1512  ! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface)
1513  cs%id_nonpenSW_diag = register_diag_field('ocean_model', 'nonpenSW', &
1514  diag%axesT1, time, &
1515  'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',&
1516  'W m-2', standard_name='nondownwelling_shortwave_flux_in_sea_water')
1517  if (cs%id_nonpenSW_diag > 0) then
1518  allocate(cs%nonpenSW_diag(isd:ied,jsd:jed))
1519  cs%nonpenSW_diag(:,:) = 0.0
1520  endif
1521  endif
1522 
1523  if (use_temperature) then
1524  call get_param(param_file, mdl, "VAR_PEN_SW", cs%var_pen_sw, &
1525  "If true, use one of the CHL_A schemes specified by "//&
1526  "OPACITY_SCHEME to determine the e-folding depth of "//&
1527  "incoming short wave radiation.", default=.false.)
1528  if (cs%var_pen_sw) then
1529 
1530  call get_param(param_file, mdl, "CHL_FROM_FILE", cs%chl_from_file, &
1531  "If true, chl_a is read from a file.", default=.true.)
1532  if (cs%chl_from_file) then
1533  call time_interp_external_init()
1534 
1535  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1536  call get_param(param_file, mdl, "CHL_FILE", chl_file, &
1537  "CHL_FILE is the file containing chl_a concentrations in "//&
1538  "the variable CHL_A. It is used when VAR_PEN_SW and "//&
1539  "CHL_FROM_FILE are true.", fail_if_missing=.true.)
1540  chl_filename = trim(slasher(inputdir))//trim(chl_file)
1541  call log_param(param_file, mdl, "INPUTDIR/CHL_FILE", chl_filename)
1542  call get_param(param_file, mdl, "CHL_VARNAME", chl_varname, &
1543  "Name of CHL_A variable in CHL_FILE.", default='CHL_A')
1544  cs%sbc_chl = init_external_field(chl_filename, trim(chl_varname), domain=g%Domain%mpp_domain)
1545  endif
1546 
1547  cs%id_chl = register_diag_field('ocean_model', 'Chl_opac', diag%axesT1, time, &
1548  'Surface chlorophyll A concentration used to find opacity', 'mg m-3')
1549  endif
1550  endif
1551 
1552  id_clock_uv_at_h = cpu_clock_id('(Ocean find_uv_at_h)', grain=clock_routine)
1553  id_clock_frazil = cpu_clock_id('(Ocean frazil)', grain=clock_routine)
1554 

References id_clock_frazil, and id_clock_uv_at_h.

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

717  type(ocean_grid_type), intent(in) :: G !< Grid type
718  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
719  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
720  integer, intent(in) :: id_MLD !< Handle (ID) of MLD diagnostic
721  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
722  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
723  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
724  !! available thermodynamic fields.
725  real, intent(in) :: densityDiff !< Density difference to determine MLD [R ~> kg m-3]
726  type(diag_ctrl), pointer :: diagPtr !< Diagnostics structure
727  integer, optional, intent(in) :: id_N2subML !< Optional handle (ID) of subML stratification
728  integer, optional, intent(in) :: id_MLDsq !< Optional handle (ID) of squared MLD
729  real, optional, intent(in) :: dz_subML !< The distance over which to calculate N2subML
730  !! or 50 m if missing [Z ~> m]
731 
732  ! Local variables
733  real, dimension(SZI_(G)) :: deltaRhoAtKm1, deltaRhoAtK ! Density differences [R ~> kg m-3].
734  real, dimension(SZI_(G)) :: pRef_MLD, pRef_N2 ! Reference pressures [Pa].
735  real, dimension(SZI_(G)) :: H_subML, dH_N2 ! Summed thicknesses used in N2 calculation [H ~> m].
736  real, dimension(SZI_(G)) :: T_subML, T_deeper ! Temperatures used in the N2 calculation [degC].
737  real, dimension(SZI_(G)) :: S_subML, S_deeper ! Salinities used in the N2 calculation [ppt].
738  real, dimension(SZI_(G)) :: rho_subML, rho_deeper ! Densities used in the N2 calculation [R ~> kg m-3].
739  real, dimension(SZI_(G)) :: dK, dKm1 ! Depths [Z ~> m].
740  real, dimension(SZI_(G)) :: rhoSurf ! Density used in finding the mixedlayer depth [R ~> kg m-3].
741  real, dimension(SZI_(G), SZJ_(G)) :: MLD ! Diagnosed mixed layer depth [Z ~> m].
742  real, dimension(SZI_(G), SZJ_(G)) :: subMLN2 ! Diagnosed stratification below ML [T-2 ~> s-2].
743  real, dimension(SZI_(G), SZJ_(G)) :: MLD2 ! Diagnosed MLD^2 [Z2 ~> m2].
744  logical, dimension(SZI_(G)) :: N2_region_set ! If true, all necessary values for calculating N2
745  ! have been stored already.
746  real :: gE_Rho0 ! The gravitational acceleration divided by a mean density [Z T-2 R-1 ~> m4 s-2 kg-1].
747  real :: dH_subML ! Depth below ML over which to diagnose stratification [H ~> m].
748  real :: aFac ! A nondimensional factor [nondim]
749  real :: ddRho ! A density difference [R ~> kg m-3]
750  integer :: i, j, is, ie, js, je, k, nz, id_N2, id_SQ
751 
752  id_n2 = -1 ; if (PRESENT(id_n2subml)) id_n2 = id_n2subml
753 
754  id_sq = -1 ; if (PRESENT(id_mldsq)) id_sq = id_mldsq
755 
756  ge_rho0 = us%L_to_Z**2*gv%g_Earth / (gv%Rho0)
757  dh_subml = 50.*gv%m_to_H ; if (present(dz_subml)) dh_subml = gv%Z_to_H*dz_subml
758 
759  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
760 
761  pref_mld(:) = 0.0
762  do j=js,je
763  do i=is,ie ; dk(i) = 0.5 * h(i,j,1) * gv%H_to_Z ; enddo ! Depth of center of surface layer
764  call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pref_mld, rhosurf, is, ie-is+1, &
765  tv%eqn_of_state, scale=us%kg_m3_to_R)
766  do i=is,ie
767  deltarhoatk(i) = 0.
768  mld(i,j) = 0.
769  if (id_n2>0) then
770  submln2(i,j) = 0.0
771  h_subml(i) = h(i,j,1) ; dh_n2(i) = 0.0
772  t_subml(i) = 0.0 ; s_subml(i) = 0.0 ; t_deeper(i) = 0.0 ; s_deeper(i) = 0.0
773  n2_region_set(i) = (g%mask2dT(i,j)<0.5) ! Only need to work on ocean points.
774  endif
775  enddo
776  do k=2,nz
777  do i=is,ie
778  dkm1(i) = dk(i) ! Depth of center of layer K-1
779  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
780  enddo
781 
782  ! Prepare to calculate stratification, N2, immediately below the mixed layer by finding
783  ! the cells that extend over at least dz_subML.
784  if (id_n2>0) then
785  do i=is,ie
786  if (mld(i,j)==0.0) then ! Still in the mixed layer.
787  h_subml(i) = h_subml(i) + h(i,j,k)
788  elseif (.not.n2_region_set(i)) then ! This block is below the mixed layer, but N2 has not been found yet.
789  if (dh_n2(i)==0.0) then ! Record the temperature, salinity, pressure, immediately below the ML
790  t_subml(i) = tv%T(i,j,k) ; s_subml(i) = tv%S(i,j,k)
791  h_subml(i) = h_subml(i) + 0.5 * h(i,j,k) ! Start midway through this layer.
792  dh_n2(i) = 0.5 * h(i,j,k)
793  elseif (dh_n2(i) + h(i,j,k) < dh_subml) then
794  dh_n2(i) = dh_n2(i) + h(i,j,k)
795  else ! This layer includes the base of the region where N2 is calculated.
796  t_deeper(i) = tv%T(i,j,k) ; s_deeper(i) = tv%S(i,j,k)
797  dh_n2(i) = dh_n2(i) + 0.5 * h(i,j,k)
798  n2_region_set(i) = .true.
799  endif
800  endif
801  enddo ! i-loop
802  endif ! id_N2>0
803 
804  ! Mixed-layer depth, using sigma-0 (surface reference pressure)
805  do i=is,ie ; deltarhoatkm1(i) = deltarhoatk(i) ; enddo ! Store value from previous iteration of K
806  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_mld, deltarhoatk, is, ie-is+1, &
807  tv%eqn_of_state, scale=us%kg_m3_to_R)
808  do i = is, ie
809  deltarhoatk(i) = deltarhoatk(i) - rhosurf(i) ! Density difference between layer K and surface
810  ddrho = deltarhoatk(i) - deltarhoatkm1(i)
811  if ((mld(i,j)==0.) .and. (ddrho>0.) .and. &
812  (deltarhoatkm1(i)<densitydiff) .and. (deltarhoatk(i)>=densitydiff)) then
813  afac = ( densitydiff - deltarhoatkm1(i) ) / ddrho
814  mld(i,j) = dk(i) * afac + dkm1(i) * (1. - afac)
815  endif
816  if (id_sq > 0) mld2(i,j) = mld(i,j)**2
817  enddo ! i-loop
818  enddo ! k-loop
819  do i=is,ie
820  if ((mld(i,j)==0.) .and. (deltarhoatk(i)<densitydiff)) mld(i,j) = dk(i) ! Assume mixing to the bottom
821  enddo
822 
823  if (id_n2>0) then ! Now actually calculate stratification, N2, below the mixed layer.
824  do i=is,ie ; pref_n2(i) = gv%H_to_Pa * (h_subml(i) + 0.5*dh_n2(i)) ; enddo
825  ! if ((.not.N2_region_set(i)) .and. (dH_N2(i) > 0.5*dH_subML)) then
826  ! ! Use whatever stratification we can, measured over whatever distance is available?
827  ! T_deeper(i) = tv%T(i,j,nz) ; S_deeper(i) = tv%S(i,j,nz)
828  ! N2_region_set(i) = .true.
829  ! endif
830  call calculate_density(t_subml, s_subml, pref_n2, rho_subml, is, ie-is+1, &
831  tv%eqn_of_state, scale=us%kg_m3_to_R)
832  call calculate_density(t_deeper, s_deeper, pref_n2, rho_deeper, is, ie-is+1, &
833  tv%eqn_of_state, scale=us%kg_m3_to_R)
834  do i=is,ie ; if ((g%mask2dT(i,j)>0.5) .and. n2_region_set(i)) then
835  submln2(i,j) = ge_rho0 * (rho_deeper(i) - rho_subml(i)) / (gv%H_to_z * dh_n2(i))
836  endif ; enddo
837  endif
838  enddo ! j-loop
839 
840  if (id_mld > 0) call post_data(id_mld, mld, diagptr)
841  if (id_n2 > 0) call post_data(id_n2, submln2, diagptr)
842  if (id_sq > 0) call post_data(id_sq, mld2, diagptr)
843 

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

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

557  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
558  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
559  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
560  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
561  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
562  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
563  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
564  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
565  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
566  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
567  intent(out) :: u_h !< Zonal velocity interpolated to h points [L T-1 ~> m s-1].
568  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
569  intent(out) :: v_h !< Meridional velocity interpolated to h points [L T-1 ~> m s-1].
570  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
571  optional, intent(in) :: ea !< The amount of fluid entrained from the layer
572  !! above within this time step [H ~> m or kg m-2].
573  !! Omitting ea is the same as setting it to 0.
574  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
575  optional, intent(in) :: eb !< The amount of fluid entrained from the layer
576  !! below within this time step [H ~> m or kg m-2].
577  !! Omitting eb is the same as setting it to 0.
578 
579  ! local variables
580  real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
581  real :: h_neglect ! A thickness that is so small it is usually lost
582  ! in roundoff and can be neglected [H ~> m or kg m-2].
583  real :: b1(SZI_(G)), d1(SZI_(G)), c1(SZI_(G),SZK_(G))
584  real :: a_n(SZI_(G)), a_s(SZI_(G)) ! Fractional weights of the neighboring
585  real :: a_e(SZI_(G)), a_w(SZI_(G)) ! velocity points, ~1/2 in the open
586  ! ocean, nondimensional.
587  real :: sum_area, Idenom
588  logical :: mix_vertically
589  integer :: i, j, k, is, ie, js, je, nz
590  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
591  call cpu_clock_begin(id_clock_uv_at_h)
592  h_neglect = gv%H_subroundoff
593 
594  mix_vertically = present(ea)
595  if (present(ea) .neqv. present(eb)) call mom_error(fatal, &
596  "find_uv_at_h: Either both ea and eb or neither one must be present "// &
597  "in call to find_uv_at_h.")
598 !$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,h,h_neglect, &
599 !$OMP eb,u_h,u,v_h,v,nz,ea) &
600 !$OMP private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1)
601  do j=js,je
602  do i=is,ie
603  sum_area = g%areaCu(i-1,j) + g%areaCu(i,j)
604  if (sum_area>0.0) then
605  idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
606  a_w(i) = g%areaCu(i-1,j) * idenom
607  a_e(i) = g%areaCu(i,j) * idenom
608  else
609  a_w(i) = 0.0 ; a_e(i) = 0.0
610  endif
611 
612  sum_area = g%areaCv(i,j-1) + g%areaCv(i,j)
613  if (sum_area>0.0) then
614  idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
615  a_s(i) = g%areaCv(i,j-1) * idenom
616  a_n(i) = g%areaCv(i,j) * idenom
617  else
618  a_s(i) = 0.0 ; a_n(i) = 0.0
619  endif
620  enddo
621 
622  if (mix_vertically) then
623  do i=is,ie
624  b_denom_1 = h(i,j,1) + h_neglect
625  b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
626  d1(i) = b_denom_1 * b1(i)
627  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))
628  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))
629  enddo
630  do k=2,nz ; do i=is,ie
631  c1(i,k) = eb(i,j,k-1) * b1(i)
632  b_denom_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h_neglect
633  b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
634  d1(i) = b_denom_1 * b1(i)
635  u_h(i,j,k) = (h(i,j,k) * (a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)) + &
636  ea(i,j,k)*u_h(i,j,k-1))*b1(i)
637  v_h(i,j,k) = (h(i,j,k) * (a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)) + &
638  ea(i,j,k)*v_h(i,j,k-1))*b1(i)
639  enddo ; enddo
640  do k=nz-1,1,-1 ; do i=is,ie
641  u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1)
642  v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1)
643  enddo ; enddo
644  else
645  do k=1,nz ; do i=is,ie
646  u_h(i,j,k) = a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)
647  v_h(i,j,k) = a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)
648  enddo ; enddo
649  endif
650  enddo
651 
652  call cpu_clock_end(id_clock_uv_at_h)

References id_clock_uv_at_h.

◆ insert_brine()

subroutine, public mom_diabatic_aux::insert_brine ( 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(forcing), intent(in)  fluxes,
integer, intent(in)  nkmb,
type(diabatic_aux_cs), intent(in)  CS,
real, intent(in)  dt,
integer, intent(in)  id_brine_lay 
)

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

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]fluxesA structure of thermodynamic surface fluxes
[in]nkmbThe number of layers in the mixed and buffer layers
[in]csThe control structure returned by a previous call to diabatic_aux_init
[in]dtThe thermodynamic time step [T ~> s].
[in]id_brine_layThe handle for a diagnostic of which layer receivees the brine.

Definition at line 386 of file MOM_diabatic_aux.F90.

386  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
387  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
388  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
389  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
390  type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any
391  !! available thermodynamic fields
392  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
393  type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
394  integer, intent(in) :: nkmb !< The number of layers in the mixed and buffer layers
395  type(diabatic_aux_CS), intent(in) :: CS !< The control structure returned by a previous
396  !! call to diabatic_aux_init
397  real, intent(in) :: dt !< The thermodynamic time step [T ~> s].
398  integer, intent(in) :: id_brine_lay !< The handle for a diagnostic of
399  !! which layer receivees the brine.
400 
401  ! local variables
402  real :: salt(SZI_(G)) ! The amount of salt rejected from sea ice [ppt R Z ~> gramSalt m-2]
403  real :: dzbr(SZI_(G)) ! Cumulative depth over which brine is distributed [H ~> m to kg m-2]
404  real :: inject_layer(SZI_(G),SZJ_(G)) ! diagnostic
405 
406  real :: p_ref_cv(SZI_(G))
407  real :: T(SZI_(G),SZK_(G))
408  real :: S(SZI_(G),SZK_(G))
409  real :: h_2d(SZI_(G),SZK_(G)) ! A 2-d slice of h with a minimum thickness [H ~> m to kg m-2]
410  real :: Rcv(SZI_(G),SZK_(G))
411  real :: s_new,R_new,t0,scale, cdz
412  integer :: i, j, k, is, ie, js, je, nz, ks
413 
414  real :: brine_dz ! minumum thickness over which to distribute brine [H ~> m or kg m-2]
415  real, parameter :: s_max = 45.0 ! salinity bound
416 
417  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
418 
419  if (.not.associated(fluxes%salt_flux)) return
420 
421  !### Injecting the brine into a single layer with a prescribed thickness seems problematic,
422  ! because it is not convergent when resolution becomes very fine. I think that this whole
423  ! subroutine needs to be revisited.- RWH
424 
425  p_ref_cv(:) = tv%P_ref
426  brine_dz = 1.0*gv%m_to_H
427 
428  inject_layer(:,:) = nz
429 
430  do j=js,je
431 
432  salt(:)=0.0 ; dzbr(:)=0.0
433 
434  do i=is,ie ; if (g%mask2dT(i,j) > 0.) then
435  salt(i) = dt * (1000. * fluxes%salt_flux(i,j))
436  endif ; enddo
437 
438  do k=1,nz
439  do i=is,ie
440  t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
441  ! avoid very small thickness
442  h_2d(i,k) = max(h(i,j,k), gv%Angstrom_H)
443  enddo
444 
445  call calculate_density(t(:,k), s(:,k), p_ref_cv, rcv(:,k), is, &
446  ie-is+1, tv%eqn_of_state)
447  enddo
448 
449  ! First, try to find an interior layer where inserting all the salt
450  ! will not cause the layer to become statically unstable.
451  ! Bias towards deeper layers.
452 
453  do k=nkmb+1,nz-1 ; do i=is,ie
454  if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
455  s_new = s(i,k) + salt(i) / (gv%H_to_RZ * h_2d(i,k))
456  t0 = t(i,k)
457  call calculate_density(t0, s_new, tv%P_Ref, r_new, tv%eqn_of_state)
458  if (r_new < 0.5*(rcv(i,k)+rcv(i,k+1)) .and. s_new<s_max) then
459  dzbr(i) = dzbr(i)+h_2d(i,k)
460  inject_layer(i,j) = min(inject_layer(i,j),real(k))
461  endif
462  endif
463  enddo ; enddo
464 
465  ! Then try to insert into buffer layers if they exist
466  do k=nkmb,gv%nkml+1,-1 ; do i=is,ie
467  if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
468  dzbr(i) = dzbr(i) + h_2d(i,k)
469  inject_layer(i,j) = min(inject_layer(i,j), real(k))
470  endif
471  enddo ; enddo
472 
473  ! finally if unable to find a layer to insert, then place in mixed layer
474 
475  do k=1,gv%nkml ; do i=is,ie
476  if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
477  dzbr(i) = dzbr(i) + h_2d(i,k)
478  inject_layer(i,j) = min(inject_layer(i,j), real(k))
479  endif
480  enddo ; enddo
481 
482 
483  do i=is,ie
484  if ((g%mask2dT(i,j) > 0.0) .and. salt(i) > 0.) then
485  ! if (dzbr(i)< brine_dz) call MOM_error(FATAL,"insert_brine: failed")
486  ks = inject_layer(i,j)
487  cdz = 0.0
488  do k=ks,nz
489  scale = h_2d(i,k) / dzbr(i)
490  cdz = cdz + h_2d(i,k)
491  !### I think that the logic of this line is wrong. Moving it down a line
492  ! would seem to make more sense. - RWH
493  if (cdz > brine_dz) exit
494  tv%S(i,j,k) = tv%S(i,j,k) + scale*salt(i) / (gv%H_to_RZ * h_2d(i,k))
495  enddo
496  endif
497  enddo
498 
499  enddo
500 
501  if (cs%id_brine_lay > 0) call post_data(cs%id_brine_lay, inject_layer, cs%diag)
502 

◆ 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(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 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]csThe control structure returned by a previous call to diabatic_aux_init.
[in]p_surfThe pressure at the ocean surface [Pa].
[in]haloHalo width over which to calculate frazil

Definition at line 103 of file MOM_diabatic_aux.F90.

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

References id_clock_frazil.

◆ 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(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.
csControl structure for diabatic_aux
opacity_cspThe control structure for the opacity module.
tracer_flow_cspA pointer to the control structure of the tracer modules.

Definition at line 657 of file MOM_diabatic_aux.F90.

657  type(optics_type), pointer :: optics !< An optics structure that has will contain
658  !! information about shortwave fluxes and absorption.
659  type(forcing), intent(inout) :: fluxes !< points to forcing fields
660  !! unused fields have NULL ptrs
661  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
662  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
663  type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
664  type(opacity_CS), pointer :: opacity_CSp !< The control structure for the opacity module.
665  type(tracer_flow_control_CS), pointer :: tracer_flow_CSp !< A pointer to the control structure of the tracer modules.
666 
667  ! Local variables
668  real, dimension(SZI_(G),SZJ_(G)) :: chl_2d !< Vertically uniform chlorophyll-A concentractions [mg m-3]
669  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: chl_3d !< The chlorophyll-A concentractions of each layer [mg m-3]
670  character(len=128) :: mesg
671  integer :: i, j, k, is, ie, js, je
672  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
673 
674  if (.not.associated(optics)) return
675 
676  if (cs%var_pen_sw) then
677  if (cs%chl_from_file) then
678  ! Only the 2-d surface chlorophyll can be read in from a file. The
679  ! same value is assumed for all layers.
680  call time_interp_external(cs%sbc_chl, cs%Time, chl_2d)
681  do j=js,je ; do i=is,ie
682  if ((g%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0)) then
683  write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&
684  & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
685  chl_2d(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j)
686  call mom_error(fatal, "MOM_diabatic_aux set_pen_shortwave: "//trim(mesg))
687  endif
688  enddo ; enddo
689 
690  if (cs%id_chl > 0) call post_data(cs%id_chl, chl_2d, cs%diag)
691 
692  call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
693  fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, opacity_csp, chl_2d=chl_2d)
694  else
695  if (.not.associated(tracer_flow_csp)) call mom_error(fatal, &
696  "The tracer flow control structure must be associated when the model sets "//&
697  "the chlorophyll internally in set_pen_shortwave.")
698  call get_chl_from_model(chl_3d, g, tracer_flow_csp)
699 
700  if (cs%id_chl > 0) call post_data(cs%id_chl, chl_3d(:,:,1), cs%diag)
701 
702  call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
703  fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, opacity_csp, chl_3d=chl_3d)
704  endif
705  else
706  call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
707  fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, opacity_csp)
708  endif
709 

References mom_tracer_flow_control::get_chl_from_model().

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

Here is the call graph for this function:
Here is the caller graph for this function:

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

508  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
509  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
510  integer, intent(in) :: is !< The start i-index to work on.
511  integer, intent(in) :: ie !< The end i-index to work on.
512  integer, intent(in) :: js !< The start j-index to work on.
513  integer, intent(in) :: je !< The end j-index to work on.
514  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hold !< The layer thicknesses before entrainment,
515  !! [H ~> m or kg m-2].
516  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: ea !< The amount of fluid entrained from the layer
517  !! above within this time step [H ~> m or kg m-2].
518  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: eb !< The amount of fluid entrained from the layer
519  !! below within this time step [H ~> m or kg m-2].
520  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: T !< Layer potential temperatures [degC].
521  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: S !< Layer salinities [ppt].
522 
523  ! Local variables
524  real :: b1(SZIB_(G)), d1(SZIB_(G)) ! b1, c1, and d1 are variables used by the
525  real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver.
526  real :: h_tr, b_denom_1
527  integer :: i, j, k
528 
529  !$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1)
530  do j=js,je
531  do i=is,ie
532  h_tr = hold(i,j,1) + gv%H_subroundoff
533  b1(i) = 1.0 / (h_tr + eb(i,j,1))
534  d1(i) = h_tr * b1(i)
535  t(i,j,1) = (b1(i)*h_tr)*t(i,j,1)
536  s(i,j,1) = (b1(i)*h_tr)*s(i,j,1)
537  enddo
538  do k=2,g%ke ; do i=is,ie
539  c1(i,k) = eb(i,j,k-1) * b1(i)
540  h_tr = hold(i,j,k) + gv%H_subroundoff
541  b_denom_1 = h_tr + d1(i)*ea(i,j,k)
542  b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
543  d1(i) = b_denom_1 * b1(i)
544  t(i,j,k) = b1(i) * (h_tr*t(i,j,k) + ea(i,j,k)*t(i,j,k-1))
545  s(i,j,k) = b1(i) * (h_tr*s(i,j,k) + ea(i,j,k)*s(i,j,k-1))
546  enddo ; enddo
547  do k=g%ke-1,1,-1 ; do i=is,ie
548  t(i,j,k) = t(i,j,k) + c1(i,k+1)*t(i,j,k+1)
549  s(i,j,k) = s(i,j,k) + c1(i,k+1)*s(i,j,k+1)
550  enddo ; enddo
551  enddo

Variable Documentation

◆ id_clock_frazil

integer mom_diabatic_aux::id_clock_frazil
private

CPU time clock IDs.

Definition at line 93 of file MOM_diabatic_aux.F90.

Referenced by diabatic_aux_init(), and make_frazil().

◆ id_clock_uv_at_h

integer mom_diabatic_aux::id_clock_uv_at_h

CPU time clock IDs.

Definition at line 93 of file MOM_diabatic_aux.F90.

93 integer :: id_clock_uv_at_h, id_clock_frazil

Referenced by diabatic_aux_init(), and find_uv_at_h().