\hypertarget{namespacemom__diabatic__aux}{}\section{mom\+\_\+diabatic\+\_\+aux Module Reference}
\label{namespacemom__diabatic__aux}\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}


\subsection{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 Mac\+Dougall and Dewar (J\+PO, 1997) is used for combined diapycnal advection and diffusion, calculated implicitly and potentially with the Richardson number dependent mixing, as described by Hallberg (M\+WR, 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. \subsection*{Data Types}
\begin{DoxyCompactItemize}
\item 
type \hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}
\begin{DoxyCompactList}\small\item\em Control structure for diabatic\+\_\+aux. \end{DoxyCompactList}\end{DoxyCompactItemize}
\subsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_adb99fd4e2092c17cd69143945733a6d9}{make\+\_\+frazil} (h, tv, G, GV, US, CS, p\+\_\+surf, halo)
\begin{DoxyCompactList}\small\item\em 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 \mbox{[}Q R Z $\sim$$>$ J m-\/2\mbox{]}) in tvfrazil. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_aa9248acf0b1cf246a79092077bcb0b46}{differential\+\_\+diffuse\+\_\+t\+\_\+s} (h, tv, visc, dt, G, GV)
\begin{DoxyCompactList}\small\item\em This subroutine applies double diffusion to T \& S, assuming no diapycal mass fluxes, using a simple triadiagonal solver. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_a737df97b7e0205b0680c3b901ff6e2ff}{adjust\+\_\+salt} (h, tv, G, GV, CS, halo)
\begin{DoxyCompactList}\small\item\em 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. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_acde4c76c9a489b82cba0bac5ec1b1198}{tridiagts} (G, GV, is, ie, js, je, hold, ea, eb, T, S)
\begin{DoxyCompactList}\small\item\em This is a simple tri-\/diagonal solver for T and S. \char`\"{}\+Simple\char`\"{} means it only uses arrays hold, ea and eb. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_a3b3416d0ee87cd1a0d54a910fb681e93}{find\+\_\+uv\+\_\+at\+\_\+h} (u, v, h, u\+\_\+h, v\+\_\+h, G, GV, US, ea, eb)
\begin{DoxyCompactList}\small\item\em This subroutine calculates u\+\_\+h and v\+\_\+h (velocities at thickness points), optionally using the entrainment amounts passed in as arguments. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_adce369573ff355b820320d9c3b048693}{set\+\_\+pen\+\_\+shortwave} (optics, fluxes, G, GV, US, CS, opacity\+\_\+\+C\+Sp, tracer\+\_\+flow\+\_\+\+C\+Sp)
\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_a43ac98433f0f2f0e6cf4f8a3c9622ec4}{diagnosemldbydensitydifference} (id\+\_\+\+M\+LD, h, tv, density\+Diff, G, GV, US, diag\+Ptr, id\+\_\+\+N2sub\+ML, id\+\_\+\+M\+L\+Dsq, dz\+\_\+sub\+ML)
\begin{DoxyCompactList}\small\item\em Diagnose a mixed layer depth (M\+LD) determined by a given density difference with the surface. This routine is appropriate in M\+O\+M\+\_\+diabatic\+\_\+driver due to its position within the time stepping. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_a695474893b92e3e3418a2bed871a5d7e}{diagnosemldbyenergy} (id\+\_\+\+M\+LD, h, tv, G, GV, US, Mixing\+\_\+\+Energy, diag\+Ptr)
\begin{DoxyCompactList}\small\item\em Diagnose a mixed layer depth (M\+LD) determined by the depth a given energy value would mix. This routine is appropriate in M\+O\+M\+\_\+diabatic\+\_\+driver due to its position within the time stepping. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_ae7a279c765fa370c302532c13a1adaca}{applyboundaryfluxesinout} (CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, aggregate\+\_\+\+F\+W\+\_\+forcing, evap\+\_\+\+C\+F\+L\+\_\+limit, minimum\+\_\+forcing\+\_\+depth, c\+T\+KE, d\+S\+V\+\_\+dT, d\+S\+V\+\_\+dS, Skin\+Buoy\+Flux)
\begin{DoxyCompactList}\small\item\em Update the thickness, temperature, and salinity due to thermodynamic boundary forcing (contained in fluxes type) applied to h, tvT and tvS, and calculate the T\+KE implications of this heating. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_a967697feb7ea22e1430d2b673ebab544}{diabatic\+\_\+aux\+\_\+init} (Time, G, GV, US, param\+\_\+file, diag, CS, use\+A\+L\+Ealgorithm, use\+\_\+e\+P\+BL)
\begin{DoxyCompactList}\small\item\em This subroutine initializes the parameters and control structure of the diabatic\+\_\+aux module. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__diabatic__aux_a071e066ca5ed6619c6a9515fb07f5567}{diabatic\+\_\+aux\+\_\+end} (CS)
\begin{DoxyCompactList}\small\item\em This subroutine initializes the control structure and any related memory for the diabatic\+\_\+aux module. \end{DoxyCompactList}\end{DoxyCompactItemize}
\subsection*{Variables}
\textbf{ }\par
\begin{DoxyCompactItemize}
\item 
\mbox{\Hypertarget{namespacemom__diabatic__aux_a1fffddaa09c8650a74e045c2fe29e256}\label{namespacemom__diabatic__aux_a1fffddaa09c8650a74e045c2fe29e256}} 
integer \hyperlink{namespacemom__diabatic__aux_a1fffddaa09c8650a74e045c2fe29e256}{id\+\_\+clock\+\_\+uv\+\_\+at\+\_\+h}
\begin{DoxyCompactList}\small\item\em C\+PU time clock I\+Ds. \end{DoxyCompactList}\item 
\mbox{\Hypertarget{namespacemom__diabatic__aux_a87671b375428f3a27c18c61d6a93d035}\label{namespacemom__diabatic__aux_a87671b375428f3a27c18c61d6a93d035}} 
integer \hyperlink{namespacemom__diabatic__aux_a87671b375428f3a27c18c61d6a93d035}{id\+\_\+clock\+\_\+frazil}
\begin{DoxyCompactList}\small\item\em C\+PU time clock I\+Ds. \end{DoxyCompactList}\end{DoxyCompactItemize}



\subsection{Function/\+Subroutine Documentation}
\mbox{\Hypertarget{namespacemom__diabatic__aux_a737df97b7e0205b0680c3b901ff6e2ff}\label{namespacemom__diabatic__aux_a737df97b7e0205b0680c3b901ff6e2ff}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!adjust\+\_\+salt@{adjust\+\_\+salt}}
\index{adjust\+\_\+salt@{adjust\+\_\+salt}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{adjust\+\_\+salt()}{adjust\_salt()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::adjust\+\_\+salt (\begin{DoxyParamCaption}\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{h,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(inout)}]{tv,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), intent(in)}]{CS,  }\item[{integer, intent(in), optional}]{halo }\end{DoxyParamCaption})}



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. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\
\hline
\mbox{\tt in}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\
\hline
\mbox{\tt in}  & {\em cs} & The control structure returned by a previous call to diabatic\+\_\+aux\+\_\+init.\\
\hline
\mbox{\tt in}  & {\em halo} & Halo width over which to work \\
\hline
\end{DoxyParams}


Definition at line 330 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
330   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: g\textcolor{comment}{    !< The ocean's grid structure}
331   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{   !< The ocean's vertical grid structure}
332   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
333                            \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{    !< Layer thicknesses [H ~> m or kg m-2]}
334   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),   \textcolor{keywordtype}{intent(inout)} :: tv\textcolor{comment}{   !< Structure containing pointers to any}
335 \textcolor{comment}{                                                 !! available thermodynamic fields.}
336   \textcolor{keywordtype}{type}(diabatic\_aux\_cs),   \textcolor{keywordtype}{intent(in)}    :: cs\textcolor{comment}{   !< The control structure returned by a previous}
337 \textcolor{comment}{                                                 !! call to diabatic\_aux\_init.}
338   \textcolor{keywordtype}{integer},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: halo\textcolor{comment}{ !< Halo width over which to work}
339 
340   \textcolor{comment}{! local variables}
341   \textcolor{keywordtype}{real} :: salt\_add\_col(szi\_(g),szj\_(g))\textcolor{comment}{ !< The accumulated salt requirement [ppt R Z ~> gSalt m-2]}
342   \textcolor{keywordtype}{real} :: s\_min\textcolor{comment}{      !< The minimum salinity [ppt].}
343   \textcolor{keywordtype}{real} :: mc\textcolor{comment}{         !< A layer's mass [R Z ~> kg m-2].}
344   \textcolor{keywordtype}{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   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(halo)) \textcolor{keywordflow}{then}
348     is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
349 \textcolor{keywordflow}{  endif}
350 
351 \textcolor{comment}{!  call cpu\_clock\_begin(id\_clock\_adjust\_salt)}
352 
353   s\_min = tv%min\_salinity
354 
355   salt\_add\_col(:,:) = 0.0
356 
357   \textcolor{comment}{!$OMP parallel do default(shared) private(mc)}
358   \textcolor{keywordflow}{do} j=js,je
359     \textcolor{keywordflow}{do} k=nz,1,-1 ; \textcolor{keywordflow}{do} i=is,ie
360       \textcolor{keywordflow}{if} ( (g%mask2dT(i,j) > 0.0) .and. &
361            ((tv%S(i,j,k) < s\_min) .or. (salt\_add\_col(i,j) > 0.0)) ) \textcolor{keywordflow}{then}
362         mc = gv%H\_to\_RZ * h(i,j,k)
363         \textcolor{keywordflow}{if} (h(i,j,k) <= 10.0*gv%Angstrom\_H) \textcolor{keywordflow}{then}
364           \textcolor{comment}{! Very thin layers should not be adjusted by the salt flux}
365           \textcolor{keywordflow}{if} (tv%S(i,j,k) < s\_min) \textcolor{keywordflow}{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 \textcolor{keywordflow}{          endif}
369         \textcolor{keywordflow}{elseif} (salt\_add\_col(i,j) + mc * (s\_min - tv%S(i,j,k)) <= 0.0) \textcolor{keywordflow}{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         \textcolor{keywordflow}{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 \textcolor{keywordflow}{        endif}
376 \textcolor{keywordflow}{      endif}
377 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
378     \textcolor{keywordflow}{do} i=is,ie
379       tv%salt\_deficit(i,j) = tv%salt\_deficit(i,j) + salt\_add\_col(i,j)
380 \textcolor{keywordflow}{    enddo}
381 \textcolor{keywordflow}{  enddo}
382 \textcolor{comment}{!  call cpu\_clock\_end(id\_clock\_adjust\_salt)}
383 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_ae7a279c765fa370c302532c13a1adaca}\label{namespacemom__diabatic__aux_ae7a279c765fa370c302532c13a1adaca}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!applyboundaryfluxesinout@{applyboundaryfluxesinout}}
\index{applyboundaryfluxesinout@{applyboundaryfluxesinout}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{applyboundaryfluxesinout()}{applyboundaryfluxesinout()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::applyboundaryfluxesinout (\begin{DoxyParamCaption}\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), pointer}]{CS,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{real, intent(in)}]{dt,  }\item[{type(forcing), intent(inout)}]{fluxes,  }\item[{type(optics\+\_\+type), pointer}]{optics,  }\item[{integer, intent(in)}]{nsw,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{h,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(inout)}]{tv,  }\item[{logical, intent(in)}]{aggregate\+\_\+\+F\+W\+\_\+forcing,  }\item[{real, intent(in)}]{evap\+\_\+\+C\+F\+L\+\_\+limit,  }\item[{real, intent(in)}]{minimum\+\_\+forcing\+\_\+depth,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out), optional}]{c\+T\+KE,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out), optional}]{d\+S\+V\+\_\+dT,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out), optional}]{d\+S\+V\+\_\+dS,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g)), intent(out), optional}]{Skin\+Buoy\+Flux }\end{DoxyParamCaption})}



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


\begin{DoxyParams}[1]{Parameters}
 & {\em cs} & Control structure for diabatic\+\_\+aux\\
\hline
\mbox{\tt in}  & {\em g} & Grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & ocean vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em dt} & Time-\/step over which forcing is applied \mbox{[}T $\sim$$>$ s\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em fluxes} & Surface fluxes container\\
\hline
 & {\em optics} & Optical properties container\\
\hline
\mbox{\tt in}  & {\em nsw} & The number of frequency bands of penetrating shortwave radiation\\
\hline
\mbox{\tt in,out}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\
\hline
\mbox{\tt in}  & {\em aggregate\+\_\+fw\+\_\+forcing} & If False, treat in/out fluxes separately.\\
\hline
\mbox{\tt in}  & {\em evap\+\_\+cfl\+\_\+limit} & The largest fraction of a layer that can be evaporated in one time-\/step \mbox{[}nondim\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em minimum\+\_\+forcing\+\_\+depth} & The smallest depth over which heat and freshwater fluxes is applied \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}.\\
\hline
\mbox{\tt out}  & {\em ctke} & Turbulent kinetic energy requirement to mix\\
\hline
\mbox{\tt out}  & {\em dsv\+\_\+dt} & Partial derivative of specific volume with\\
\hline
\mbox{\tt out}  & {\em dsv\+\_\+ds} & Partial derivative of specific volume with\\
\hline
\mbox{\tt out}  & {\em skinbuoyflux} & Buoyancy flux at surface \mbox{[}Z2 T-\/3 $\sim$$>$ m2 s-\/3\mbox{]}. \\
\hline
\end{DoxyParams}


Definition at line 939 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


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



This subroutine initializes the control structure and any related memory for the diabatic\+\_\+aux module. 


\begin{DoxyParams}{Parameters}
{\em cs} & The control structure returned by a previous call to diabatic\+\_\+aux\+\_\+init; it is deallocated here. \\
\hline
\end{DoxyParams}


Definition at line 1652 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
1652   \textcolor{keywordtype}{type}(diabatic\_aux\_cs), \textcolor{keywordtype}{pointer} :: cs\textcolor{comment}{ !< The control structure returned by a previous}
1653 \textcolor{comment}{                                       !! call to diabatic\_aux\_init; it is deallocated here.}
1654 
1655   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{return}
1656 
1657   \textcolor{keywordflow}{if} (cs%id\_createdH       >0) \textcolor{keyword}{deallocate}(cs%createdH)
1658   \textcolor{keywordflow}{if} (cs%id\_penSW\_diag     >0) \textcolor{keyword}{deallocate}(cs%penSW\_diag)
1659   \textcolor{keywordflow}{if} (cs%id\_penSWflux\_diag >0) \textcolor{keyword}{deallocate}(cs%penSWflux\_diag)
1660   \textcolor{keywordflow}{if} (cs%id\_nonpenSW\_diag  >0) \textcolor{keyword}{deallocate}(cs%nonpenSW\_diag)
1661 
1662   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{deallocate}(cs)
1663 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_a967697feb7ea22e1430d2b673ebab544}\label{namespacemom__diabatic__aux_a967697feb7ea22e1430d2b673ebab544}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!diabatic\+\_\+aux\+\_\+init@{diabatic\+\_\+aux\+\_\+init}}
\index{diabatic\+\_\+aux\+\_\+init@{diabatic\+\_\+aux\+\_\+init}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{diabatic\+\_\+aux\+\_\+init()}{diabatic\_aux\_init()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::diabatic\+\_\+aux\+\_\+init (\begin{DoxyParamCaption}\item[{type(time\+\_\+type), intent(in), target}]{Time,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(param\+\_\+file\+\_\+type), intent(in)}]{param\+\_\+file,  }\item[{type(diag\+\_\+ctrl), intent(inout), target}]{diag,  }\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), pointer}]{CS,  }\item[{logical, intent(in)}]{use\+A\+L\+Ealgorithm,  }\item[{logical, intent(in)}]{use\+\_\+e\+P\+BL }\end{DoxyParamCaption})}



This subroutine initializes the parameters and control structure of the diabatic\+\_\+aux module. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em time} & The current model time.\\
\hline
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em param\+\_\+file} & A structure to parse for run-\/time parameters\\
\hline
\mbox{\tt in,out}  & {\em diag} & A structure used to regulate diagnostic output\\
\hline
 & {\em cs} & A pointer to the control structure for the diabatic\+\_\+aux module, which is initialized here.\\
\hline
\mbox{\tt in}  & {\em usealealgorithm} & If true, use the A\+LE algorithm rather than layered mode.\\
\hline
\mbox{\tt in}  & {\em use\+\_\+epbl} & If true, use the implicit energetics planetary boundary layer scheme to determine the diffusivity in the surface boundary layer. \\
\hline
\end{DoxyParams}


Definition at line 1488 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
1488   \textcolor{keywordtype}{type}(time\_type), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(in)}    :: time\textcolor{comment}{ !< The current model time.}
1489   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: g\textcolor{comment}{    !< The ocean's grid structure}
1490   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{   !< The ocean's vertical grid structure}
1491   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)}    :: us\textcolor{comment}{   !< A dimensional unit scaling type}
1492   \textcolor{keywordtype}{type}(param\_file\_type),   \textcolor{keywordtype}{intent(in)}    :: param\_file\textcolor{comment}{ !< A structure to parse for run-time parameters}
1493   \textcolor{keywordtype}{type}(diag\_ctrl), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(inout)} :: diag\textcolor{comment}{ !< A structure used to regulate diagnostic output}
1494   \textcolor{keywordtype}{type}(diabatic\_aux\_cs),   \textcolor{keywordtype}{pointer}       :: cs\textcolor{comment}{   !< A pointer to the control structure for the}
1495 \textcolor{comment}{                                                 !! diabatic\_aux module, which is initialized here.}
1496   \textcolor{keywordtype}{logical},                 \textcolor{keywordtype}{intent(in)}    :: usealealgorithm\textcolor{comment}{ !< If true, use the ALE algorithm rather}
1497 \textcolor{comment}{                                                 !! than layered mode.}
1498   \textcolor{keywordtype}{logical},                 \textcolor{keywordtype}{intent(in)}    :: use\_epbl\textcolor{comment}{ !< If true, use the implicit energetics planetary}
1499 \textcolor{comment}{                                                 !! boundary layer scheme to determine the diffusivity}
1500 \textcolor{comment}{                                                 !! in the surface boundary layer.}
1501 
1502 \textcolor{comment}{! This "include" declares and sets the variable "version".}
1503 \textcolor{preprocessor}{#include "version\_variable.h"}
1504 \textcolor{preprocessor}{}  \textcolor{keywordtype}{character(len=40)}  :: mdl  = \textcolor{stringliteral}{"MOM\_diabatic\_aux"} \textcolor{comment}{! This module's name.}
1505   \textcolor{keywordtype}{character(len=48)}  :: thickness\_units
1506   \textcolor{keywordtype}{character(len=200)} :: inputdir   \textcolor{comment}{! The directory where NetCDF input files}
1507   \textcolor{keywordtype}{character(len=240)} :: chl\_filename \textcolor{comment}{! A file from which chl\_a concentrations are to be read.}
1508   \textcolor{keywordtype}{character(len=128)} :: chl\_file \textcolor{comment}{! Data containing chl\_a concentrations. Used}
1509                                  \textcolor{comment}{! when var\_pen\_sw is defined and reading from file.}
1510   \textcolor{keywordtype}{character(len=32)}  :: chl\_varname \textcolor{comment}{! Name of chl\_a variable in chl\_file.}
1511   \textcolor{keywordtype}{logical} :: use\_temperature     \textcolor{comment}{! True if thermodynamics are enabled.}
1512   \textcolor{keywordtype}{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   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{then}
1517     \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"diabatic\_aux\_init called with an "}// &
1518                             \textcolor{stringliteral}{"associated control structure."})
1519     \textcolor{keywordflow}{return}
1520   \textcolor{keywordflow}{else}
1521     \textcolor{keyword}{allocate}(cs)
1522 \textcolor{keywordflow}{   endif}
1523 
1524   cs%diag => diag
1525   cs%Time => time
1526 
1527 \textcolor{comment}{! Set default, read and log parameters}
1528   \textcolor{keyword}{call }log\_version(param\_file, mdl, version, &
1529                    \textcolor{stringliteral}{"The following parameters are used for auxiliary diabatic processes."})
1530 
1531   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"ENABLE\_THERMODYNAMICS"}, use\_temperature, &
1532                  \textcolor{stringliteral}{"If true, temperature and salinity are used as state "}//&
1533                  \textcolor{stringliteral}{"variables."}, default=.true.)
1534 
1535   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"RECLAIM\_FRAZIL"}, cs%reclaim\_frazil, &
1536                  \textcolor{stringliteral}{"If true, try to use any frazil heat deficit to cool any "}//&
1537                  \textcolor{stringliteral}{"overlying layers down to the freezing point, thereby "}//&
1538                  \textcolor{stringliteral}{"avoiding the creation of thin ice when the SST is above "}//&
1539                  \textcolor{stringliteral}{"the freezing point."}, default=.true.)
1540   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"PRESSURE\_DEPENDENT\_FRAZIL"}, &
1541                                 cs%pressure\_dependent\_frazil, &
1542                  \textcolor{stringliteral}{"If true, use a pressure dependent freezing temperature "}//&
1543                  \textcolor{stringliteral}{"when making frazil. The default is false, which will be "}//&
1544                  \textcolor{stringliteral}{"faster but is inappropriate with ice-shelf cavities."}, &
1545                  default=.false.)
1546 
1547   \textcolor{keywordflow}{if} (use\_epbl) \textcolor{keywordflow}{then}
1548     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"IGNORE\_FLUXES\_OVER\_LAND"}, cs%ignore\_fluxes\_over\_land,&
1549          \textcolor{stringliteral}{"If true, the model does not check if fluxes are being applied "}//&
1550          \textcolor{stringliteral}{"over land points. This is needed when the ocean is coupled "}//&
1551          \textcolor{stringliteral}{"with ice shelves and sea ice, since the sea ice mask needs to "}//&
1552          \textcolor{stringliteral}{"be different than the ocean mask to avoid sea ice formation "}//&
1553          \textcolor{stringliteral}{"under ice shelves. This flag only works when use\_ePBL = True."}, default=.false.)
1554     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DO\_RIVERMIX"}, cs%do\_rivermix, &
1555                  \textcolor{stringliteral}{"If true, apply additional mixing wherever there is "}//&
1556                  \textcolor{stringliteral}{"runoff, so that it is mixed down to RIVERMIX\_DEPTH "}//&
1557                  \textcolor{stringliteral}{"if the ocean is that deep."}, default=.false.)
1558     \textcolor{keywordflow}{if} (cs%do\_rivermix) &
1559       \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"RIVERMIX\_DEPTH"}, cs%rivermix\_depth, &
1560                  \textcolor{stringliteral}{"The depth to which rivers are mixed if DO\_RIVERMIX is "}//&
1561                  \textcolor{stringliteral}{"defined."}, units=\textcolor{stringliteral}{"m"}, default=0.0, scale=us%m\_to\_Z)
1562   \textcolor{keywordflow}{else}
1563     cs%do\_rivermix = .false. ; cs%rivermix\_depth = 0.0 ; cs%ignore\_fluxes\_over\_land = .false.
1564 \textcolor{keywordflow}{  endif}
1565 
1566   \textcolor{keywordflow}{if} (gv%nkml == 0) \textcolor{keywordflow}{then}
1567     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_RIVER\_HEAT\_CONTENT"}, cs%use\_river\_heat\_content, &
1568                    \textcolor{stringliteral}{"If true, use the fluxes%runoff\_Hflx field to set the "}//&
1569                    \textcolor{stringliteral}{"heat carried by runoff, instead of using SST*CP*liq\_runoff."}, &
1570                    default=.false.)
1571     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_CALVING\_HEAT\_CONTENT"}, cs%use\_calving\_heat\_content, &
1572                    \textcolor{stringliteral}{"If true, use the fluxes%calving\_Hflx field to set the "}//&
1573                    \textcolor{stringliteral}{"heat carried by runoff, instead of using SST*CP*froz\_runoff."}, &
1574                    default=.false.)
1575   \textcolor{keywordflow}{else}
1576     cs%use\_river\_heat\_content = .false.
1577     cs%use\_calving\_heat\_content = .false.
1578 \textcolor{keywordflow}{  endif}
1579 
1580   \textcolor{keywordflow}{if} (usealealgorithm) \textcolor{keywordflow}{then}
1581     cs%id\_createdH = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'},\textcolor{stringliteral}{"created\_H"},diag%axesT1, &
1582         time, \textcolor{stringliteral}{"The volume flux added to stop the ocean from drying out and becoming negative in depth"}, &
1583         \textcolor{stringliteral}{"m s-1"}, conversion=gv%H\_to\_m*us%s\_to\_T)
1584     \textcolor{keywordflow}{if} (cs%id\_createdH>0) \textcolor{keyword}{allocate}(cs%createdH(isd:ied,jsd:jed))
1585 
1586     \textcolor{comment}{! diagnostic for heating of a grid cell from convergence of SW heat into the cell}
1587     cs%id\_penSW\_diag = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'rsdoabsorb'},                     &
1588           diag%axesTL, time, \textcolor{stringliteral}{'Convergence of Penetrative Shortwave Flux in Sea Water Layer'},&
1589           \textcolor{stringliteral}{'W m-2'}, conversion=us%QRZ\_T\_to\_W\_m2, &
1590           standard\_name=\textcolor{stringliteral}{'net\_rate\_of\_absorption\_of\_shortwave\_energy\_in\_ocean\_layer'}, v\_extensive=.true.)
1591 
1592     \textcolor{comment}{! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces)}
1593     \textcolor{comment}{! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock).}
1594     cs%id\_penSWflux\_diag = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'rsdo'},                               &
1595           diag%axesTi, time, \textcolor{stringliteral}{'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface'},&
1596           \textcolor{stringliteral}{'W m-2'}, conversion=us%QRZ\_T\_to\_W\_m2, standard\_name=\textcolor{stringliteral}{'downwelling\_shortwave\_flux\_in\_sea\_water'})
1597 
1598     \textcolor{comment}{! need both arrays for the SW diagnostics (one for flux, one for convergence)}
1599     \textcolor{keywordflow}{if} (cs%id\_penSW\_diag>0 .or. cs%id\_penSWflux\_diag>0) \textcolor{keywordflow}{then}
1600       \textcolor{keyword}{allocate}(cs%penSW\_diag(isd:ied,jsd:jed,nz)) ; cs%penSW\_diag(:,:,:) = 0.0
1601       \textcolor{keyword}{allocate}(cs%penSWflux\_diag(isd:ied,jsd:jed,nz+1)) ; cs%penSWflux\_diag(:,:,:) = 0.0
1602 \textcolor{keywordflow}{    endif}
1603 
1604     \textcolor{comment}{! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface)}
1605     cs%id\_nonpenSW\_diag = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'nonpenSW'},                       &
1606           diag%axesT1, time,                                                                   &
1607           \textcolor{stringliteral}{'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)'},&
1608           \textcolor{stringliteral}{'W m-2'}, conversion=us%QRZ\_T\_to\_W\_m2, &
1609           standard\_name=\textcolor{stringliteral}{'nondownwelling\_shortwave\_flux\_in\_sea\_water'})
1610     \textcolor{keywordflow}{if} (cs%id\_nonpenSW\_diag > 0) \textcolor{keywordflow}{then}
1611       \textcolor{keyword}{allocate}(cs%nonpenSW\_diag(isd:ied,jsd:jed)) ; cs%nonpenSW\_diag(:,:) = 0.0
1612 \textcolor{keywordflow}{    endif}
1613 \textcolor{keywordflow}{  endif}
1614 
1615   \textcolor{keywordflow}{if} (use\_temperature) \textcolor{keywordflow}{then}
1616     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"VAR\_PEN\_SW"}, cs%var\_pen\_sw, &
1617                    \textcolor{stringliteral}{"If true, use one of the CHL\_A schemes specified by "}//&
1618                    \textcolor{stringliteral}{"OPACITY\_SCHEME to determine the e-folding depth of "}//&
1619                    \textcolor{stringliteral}{"incoming short wave radiation."}, default=.false.)
1620     \textcolor{keywordflow}{if} (cs%var\_pen\_sw) \textcolor{keywordflow}{then}
1621 
1622       \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CHL\_FROM\_FILE"}, cs%chl\_from\_file, &
1623                    \textcolor{stringliteral}{"If true, chl\_a is read from a file."}, default=.true.)
1624       \textcolor{keywordflow}{if} (cs%chl\_from\_file) \textcolor{keywordflow}{then}
1625         \textcolor{keyword}{call }time\_interp\_external\_init()
1626 
1627         \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"INPUTDIR"}, inputdir, default=\textcolor{stringliteral}{"."})
1628         \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CHL\_FILE"}, chl\_file, &
1629                    \textcolor{stringliteral}{"CHL\_FILE is the file containing chl\_a concentrations in "}//&
1630                    \textcolor{stringliteral}{"the variable CHL\_A. It is used when VAR\_PEN\_SW and "}//&
1631                    \textcolor{stringliteral}{"CHL\_FROM\_FILE are true."}, fail\_if\_missing=.true.)
1632         chl\_filename = trim(slasher(inputdir))//trim(chl\_file)
1633         \textcolor{keyword}{call }log\_param(param\_file, mdl, \textcolor{stringliteral}{"INPUTDIR/CHL\_FILE"}, chl\_filename)
1634         \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CHL\_VARNAME"}, chl\_varname, &
1635                    \textcolor{stringliteral}{"Name of CHL\_A variable in CHL\_FILE."}, default=\textcolor{stringliteral}{'CHL\_A'})
1636         cs%sbc\_chl = init\_external\_field(chl\_filename, trim(chl\_varname), domain=g%Domain%mpp\_domain)
1637 \textcolor{keywordflow}{      endif}
1638 
1639       cs%id\_chl = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'Chl\_opac'}, diag%axesT1, time, &
1640           \textcolor{stringliteral}{'Surface chlorophyll A concentration used to find opacity'}, \textcolor{stringliteral}{'mg m-3'})
1641 \textcolor{keywordflow}{    endif}
1642 \textcolor{keywordflow}{  endif}
1643 
1644   id\_clock\_uv\_at\_h = cpu\_clock\_id(\textcolor{stringliteral}{'(Ocean find\_uv\_at\_h)'}, grain=clock\_routine)
1645   id\_clock\_frazil  = cpu\_clock\_id(\textcolor{stringliteral}{'(Ocean frazil)'}, grain=clock\_routine)
1646 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_a43ac98433f0f2f0e6cf4f8a3c9622ec4}\label{namespacemom__diabatic__aux_a43ac98433f0f2f0e6cf4f8a3c9622ec4}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!diagnosemldbydensitydifference@{diagnosemldbydensitydifference}}
\index{diagnosemldbydensitydifference@{diagnosemldbydensitydifference}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{diagnosemldbydensitydifference()}{diagnosemldbydensitydifference()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::diagnosemldbydensitydifference (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{id\+\_\+\+M\+LD,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv,  }\item[{real, intent(in)}]{density\+Diff,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(diag\+\_\+ctrl), pointer}]{diag\+Ptr,  }\item[{integer, intent(in), optional}]{id\+\_\+\+N2sub\+ML,  }\item[{integer, intent(in), optional}]{id\+\_\+\+M\+L\+Dsq,  }\item[{real, intent(in), optional}]{dz\+\_\+sub\+ML }\end{DoxyParamCaption})}



Diagnose a mixed layer depth (M\+LD) determined by a given density difference with the surface. This routine is appropriate in M\+O\+M\+\_\+diabatic\+\_\+driver due to its position within the time stepping. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & Grid type\\
\hline
\mbox{\tt in}  & {\em gv} & ocean vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em id\+\_\+mld} & Handle (ID) of M\+LD diagnostic\\
\hline
\mbox{\tt in}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\
\hline
\mbox{\tt in}  & {\em densitydiff} & Density difference to determine M\+LD \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]}\\
\hline
 & {\em diagptr} & Diagnostics structure\\
\hline
\mbox{\tt in}  & {\em id\+\_\+n2subml} & Optional handle (ID) of sub\+ML stratification\\
\hline
\mbox{\tt in}  & {\em id\+\_\+mldsq} & Optional handle (ID) of squared M\+LD\\
\hline
\mbox{\tt in}  & {\em dz\+\_\+subml} & The distance over which to calculate N2sub\+ML or 50 m if missing \mbox{[}Z $\sim$$>$ m\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 600 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
600   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{           !< Grid type}
601   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{          !< ocean vertical grid structure}
602   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{          !< A dimensional unit scaling type}
603   \textcolor{keywordtype}{integer},                 \textcolor{keywordtype}{intent(in)} :: id\_mld\textcolor{comment}{      !< Handle (ID) of MLD diagnostic}
604   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
605                            \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{           !< Layer thickness [H ~> m or kg m-2]}
606   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),   \textcolor{keywordtype}{intent(in)} :: tv\textcolor{comment}{          !< Structure containing pointers to any}
607 \textcolor{comment}{                                                     !! available thermodynamic fields.}
608   \textcolor{keywordtype}{real},                    \textcolor{keywordtype}{intent(in)} :: densitydiff\textcolor{comment}{ !< Density difference to determine MLD [R ~> kg m-3]}
609   \textcolor{keywordtype}{type}(diag\_ctrl),         \textcolor{keywordtype}{pointer}    :: diagptr\textcolor{comment}{     !< Diagnostics structure}
610   \textcolor{keywordtype}{integer},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: id\_n2subml\textcolor{comment}{  !< Optional handle (ID) of subML stratification}
611   \textcolor{keywordtype}{integer},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: id\_mldsq\textcolor{comment}{    !< Optional handle (ID) of squared MLD}
612   \textcolor{keywordtype}{real},          \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: dz\_subml\textcolor{comment}{    !< The distance over which to calculate N2subML}
613 \textcolor{comment}{                                                     !! or 50 m if missing [Z ~> m]}
614 
615   \textcolor{comment}{! Local variables}
616   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: deltarhoatkm1, deltarhoatk \textcolor{comment}{! Density differences [R ~> kg m-3].}
617   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: pref\_mld, pref\_n2 \textcolor{comment}{! Reference pressures [R L2 T-2 ~> Pa].}
618   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: h\_subml, dh\_n2    \textcolor{comment}{! Summed thicknesses used in N2 calculation [H ~> m].}
619   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: t\_subml, t\_deeper \textcolor{comment}{! Temperatures used in the N2 calculation [degC].}
620   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: s\_subml, s\_deeper \textcolor{comment}{! Salinities used in the N2 calculation [ppt].}
621   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: rho\_subml, rho\_deeper \textcolor{comment}{! Densities used in the N2 calculation [R ~> kg m-3].}
622   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: dk, dkm1          \textcolor{comment}{! Depths [Z ~> m].}
623   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: rhosurf          \textcolor{comment}{! Density used in finding the mixedlayer depth [R ~> kg
       m-3].}
624   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G))} :: mld     \textcolor{comment}{! Diagnosed mixed layer depth [Z ~> m].}
625   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G))} :: submln2 \textcolor{comment}{! Diagnosed stratification below ML [T-2 ~> s-2].}
626   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G))} :: mld2    \textcolor{comment}{! Diagnosed MLD^2 [Z2 ~> m2].}
627   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: n2\_region\_set \textcolor{comment}{! If true, all necessary values for calculating N2}
628                                                \textcolor{comment}{! have been stored already.}
629   \textcolor{keywordtype}{real} :: ge\_rho0          \textcolor{comment}{! The gravitational acceleration divided by a mean density [Z T-2 R-1 ~> m4 s-2
       kg-1].}
630   \textcolor{keywordtype}{real} :: dh\_subml         \textcolor{comment}{! Depth below ML over which to diagnose stratification [H ~> m].}
631   \textcolor{keywordtype}{real} :: afac             \textcolor{comment}{! A nondimensional factor [nondim]}
632   \textcolor{keywordtype}{real} :: ddrho            \textcolor{comment}{! A density difference [R ~> kg m-3]}
633   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(2)} :: eosdom \textcolor{comment}{! The i-computational domain for the equation of state}
634   \textcolor{keywordtype}{integer} :: i, j, is, ie, js, je, k, nz, id\_n2, id\_sq
635 
636   id\_n2 = -1 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{PRESENT}(id\_n2subml)) id\_n2 = id\_n2subml
637 
638   id\_sq = -1 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{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  ; \textcolor{keywordflow}{if} (\textcolor{keyword}{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   \textcolor{keywordflow}{do} j=js,je
648     \textcolor{keywordflow}{do} i=is,ie ; dk(i) = 0.5 * h(i,j,1) * gv%H\_to\_Z ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! Depth of center of surface layer}
649     \textcolor{keyword}{call }calculate\_density(tv%T(:,j,1), tv%S(:,j,1), pref\_mld, rhosurf, tv%eqn\_of\_state, eosdom)
650     \textcolor{keywordflow}{do} i=is,ie
651       deltarhoatk(i) = 0.
652       mld(i,j) = 0.
653       \textcolor{keywordflow}{if} (id\_n2>0) \textcolor{keywordflow}{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) \textcolor{comment}{! Only need to work on ocean points.}
658 \textcolor{keywordflow}{      endif}
659 \textcolor{keywordflow}{    enddo}
660     \textcolor{keywordflow}{do} k=2,nz
661       \textcolor{keywordflow}{do} i=is,ie
662         dkm1(i) = dk(i) \textcolor{comment}{! 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 \textcolor{comment}{! Depth of center of layer K}
664 \textcolor{keywordflow}{      enddo}
665 
666       \textcolor{comment}{! Prepare to calculate stratification, N2, immediately below the mixed layer by finding}
667       \textcolor{comment}{! the cells that extend over at least dz\_subML.}
668       \textcolor{keywordflow}{if} (id\_n2>0) \textcolor{keywordflow}{then}
669          \textcolor{keywordflow}{do} i=is,ie
670           \textcolor{keywordflow}{if} (mld(i,j)==0.0) \textcolor{keywordflow}{then}  \textcolor{comment}{! Still in the mixed layer.}
671             h\_subml(i) = h\_subml(i) + h(i,j,k)
672           \textcolor{keywordflow}{elseif} (.not.n2\_region\_set(i)) \textcolor{keywordflow}{then} \textcolor{comment}{! This block is below the mixed layer, but N2 has not been
       found yet.}
673             \textcolor{keywordflow}{if} (dh\_n2(i)==0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! 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) \textcolor{comment}{! Start midway through this layer.}
676               dh\_n2(i) = 0.5 * h(i,j,k)
677             \textcolor{keywordflow}{elseif} (dh\_n2(i) + h(i,j,k) < dh\_subml) \textcolor{keywordflow}{then}
678               dh\_n2(i) = dh\_n2(i) + h(i,j,k)
679             \textcolor{keywordflow}{else}  \textcolor{comment}{! 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 \textcolor{keywordflow}{            endif}
684 \textcolor{keywordflow}{          endif}
685 \textcolor{keywordflow}{        enddo} \textcolor{comment}{! i-loop}
686 \textcolor{keywordflow}{      endif} \textcolor{comment}{! id\_N2>0}
687 
688       \textcolor{comment}{! Mixed-layer depth, using sigma-0 (surface reference pressure)}
689       \textcolor{keywordflow}{do} i=is,ie ; deltarhoatkm1(i) = deltarhoatk(i) ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! Store value from previous iteration of K}
690       \textcolor{keyword}{call }calculate\_density(tv%T(:,j,k), tv%S(:,j,k), pref\_mld, deltarhoatk, tv%eqn\_of\_state, eosdom)
691       \textcolor{keywordflow}{do} i = is, ie
692         deltarhoatk(i) = deltarhoatk(i) - rhosurf(i) \textcolor{comment}{! Density difference between layer K and surface}
693         ddrho = deltarhoatk(i) - deltarhoatkm1(i)
694         \textcolor{keywordflow}{if} ((mld(i,j)==0.) .and. (ddrho>0.) .and. &
695             (deltarhoatkm1(i)<densitydiff) .and. (deltarhoatk(i)>=densitydiff)) \textcolor{keywordflow}{then}
696           afac = ( densitydiff - deltarhoatkm1(i) ) / ddrho
697           mld(i,j) = dk(i) * afac + dkm1(i) * (1. - afac)
698 \textcolor{keywordflow}{        endif}
699         \textcolor{keywordflow}{if} (id\_sq > 0) mld2(i,j) = mld(i,j)**2
700 \textcolor{keywordflow}{      enddo} \textcolor{comment}{! i-loop}
701 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! k-loop}
702     \textcolor{keywordflow}{do} i=is,ie
703       \textcolor{keywordflow}{if} ((mld(i,j)==0.) .and. (deltarhoatk(i)<densitydiff)) mld(i,j) = dk(i) \textcolor{comment}{! Assume mixing to the bottom}
704 \textcolor{keywordflow}{    enddo}
705 
706     \textcolor{keywordflow}{if} (id\_n2>0) \textcolor{keywordflow}{then}  \textcolor{comment}{! Now actually calculate stratification, N2, below the mixed layer.}
707       \textcolor{keywordflow}{do} i=is,ie ; pref\_n2(i) = (gv%g\_Earth * gv%H\_to\_RZ) * (h\_subml(i) + 0.5*dh\_n2(i)) ;\textcolor{keywordflow}{ enddo}
708       \textcolor{comment}{! if ((.not.N2\_region\_set(i)) .and. (dH\_N2(i) > 0.5*dH\_subML)) then}
709       \textcolor{comment}{!    ! Use whatever stratification we can, measured over whatever distance is available?}
710       \textcolor{comment}{!    T\_deeper(i) = tv%T(i,j,nz) ; S\_deeper(i) = tv%S(i,j,nz)}
711       \textcolor{comment}{!    N2\_region\_set(i) = .true.}
712       \textcolor{comment}{! endif}
713       \textcolor{keyword}{call }calculate\_density(t\_subml, s\_subml, pref\_n2, rho\_subml, tv%eqn\_of\_state, eosdom)
714       \textcolor{keyword}{call }calculate\_density(t\_deeper, s\_deeper, pref\_n2, rho\_deeper, tv%eqn\_of\_state, eosdom)
715       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} ((g%mask2dT(i,j)>0.5) .and. n2\_region\_set(i)) \textcolor{keywordflow}{then}
716         submln2(i,j) =  ge\_rho0 * (rho\_deeper(i) - rho\_subml(i)) / (gv%H\_to\_z * dh\_n2(i))
717 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
718 \textcolor{keywordflow}{    endif}
719 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! j-loop}
720 
721   \textcolor{keywordflow}{if} (id\_mld > 0) \textcolor{keyword}{call }post\_data(id\_mld, mld, diagptr)
722   \textcolor{keywordflow}{if} (id\_n2 > 0)  \textcolor{keyword}{call }post\_data(id\_n2, submln2, diagptr)
723   \textcolor{keywordflow}{if} (id\_sq > 0)  \textcolor{keyword}{call }post\_data(id\_sq, mld2, diagptr)
724 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_a695474893b92e3e3418a2bed871a5d7e}\label{namespacemom__diabatic__aux_a695474893b92e3e3418a2bed871a5d7e}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!diagnosemldbyenergy@{diagnosemldbyenergy}}
\index{diagnosemldbyenergy@{diagnosemldbyenergy}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{diagnosemldbyenergy()}{diagnosemldbyenergy()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::diagnosemldbyenergy (\begin{DoxyParamCaption}\item[{integer, dimension(3), intent(in)}]{id\+\_\+\+M\+LD,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{h,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{real, dimension(3), intent(in)}]{Mixing\+\_\+\+Energy,  }\item[{type(diag\+\_\+ctrl), pointer}]{diag\+Ptr }\end{DoxyParamCaption})}



Diagnose a mixed layer depth (M\+LD) determined by the depth a given energy value would mix. This routine is appropriate in M\+O\+M\+\_\+diabatic\+\_\+driver due to its position within the time stepping. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em id\+\_\+mld} & Energy output diag I\+Ds\\
\hline
\mbox{\tt in}  & {\em g} & Grid type\\
\hline
\mbox{\tt in}  & {\em gv} & ocean vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em mixing\+\_\+energy} & Energy values for up to 3 M\+L\+Ds\\
\hline
\mbox{\tt in}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\
\hline
 & {\em diagptr} & Diagnostics structure \\
\hline
\end{DoxyParams}


Definition at line 730 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
730   \textcolor{comment}{! Author: Brandon Reichl}
731   \textcolor{comment}{! Date: October 2, 2020}
732   \textcolor{comment}{! //}
733   \textcolor{comment}{! *Note that gravity is assumed constant everywhere and divided out of all calculations.}
734   \textcolor{comment}{!}
735   \textcolor{comment}{! This code has been written to step through the columns layer by layer, summing the PE}
736   \textcolor{comment}{! change inferred by mixing the layer with all layers above.  When the change exceeds a}
737   \textcolor{comment}{! threshold (determined by input array Mixing\_Energy), the code needs to solve for how far}
738   \textcolor{comment}{! into this layer the threshold PE change occurs (assuming constant density layers).}
739   \textcolor{comment}{! This is expressed here via solving the function F(X) = 0 where:}
740   \textcolor{comment}{! F(X) = 0.5 * ( Ca*X^3/(D1+X) + Cb*X^2/(D1+X) + Cc*X/(D1+X) + Dc/(D1+X)}
741   \textcolor{comment}{!                + Ca2*X^2 + Cb2*X + Cc2)}
742   \textcolor{comment}{! where all coefficients are determined by the previous mixed layer depth, the}
743   \textcolor{comment}{! density of the previous mixed layer, the present layer thickness, and the present}
744   \textcolor{comment}{! layer density.  This equation is worked out by computing the total PE assuming constant}
745   \textcolor{comment}{! density in the mixed layer as well as in the remaining part of the present layer that is}
746   \textcolor{comment}{! not mixed.}
747   \textcolor{comment}{! To solve for X in this equation a Newton's method iteration is employed, which}
748   \textcolor{comment}{! converges extremely quickly (usually 1 guess) since this equation turns out being rather}
749   \textcolor{comment}{! lienar for PE change with increasing X.}
750   \textcolor{comment}{! Input parameters:}
751   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(3)},   \textcolor{keywordtype}{intent(in)} :: id\_mld\textcolor{comment}{      !< Energy output diag IDs}
752   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{           !< Grid type}
753   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{          !< ocean vertical grid structure}
754   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{          !< A dimensional unit scaling type}
755   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(3)},      \textcolor{keywordtype}{intent(in)} :: mixing\_energy\textcolor{comment}{ !< Energy values for up to 3 MLDs}
756   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
757                            \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{           !< Layer thickness [H ~> m or kg m-2]}
758   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),   \textcolor{keywordtype}{intent(in)} :: tv\textcolor{comment}{          !< Structure containing pointers to any}
759 \textcolor{comment}{                                                     !! available thermodynamic fields.}
760   \textcolor{keywordtype}{type}(diag\_ctrl),         \textcolor{keywordtype}{pointer}    :: diagptr\textcolor{comment}{     !< Diagnostics structure}
761 
762   \textcolor{comment}{! Local variables}
763   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G),3)} :: mld     \textcolor{comment}{! Diagnosed mixed layer depth [Z ~> m].}
764   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: z\_l, z\_u, dz, rho\_c, pref\_mld
765   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(3)} :: pe\_threshold
766 
767   \textcolor{keywordtype}{real} :: ig, e\_g
768   \textcolor{keywordtype}{real} :: pe\_threshold\_fraction, pe, pe\_mixed, pe\_mixed\_tst
769   \textcolor{keywordtype}{real} :: rhodz\_ml, h\_ml, rhodz\_ml\_tst, h\_ml\_tst
770   \textcolor{keywordtype}{real} :: rho\_ml
771 
772   \textcolor{keywordtype}{real} :: r1, d1, r2, d2
773   \textcolor{keywordtype}{real} :: ca, cb,d ,cc, cd, ca2, cb2, c, cc2
774   \textcolor{keywordtype}{real} :: gx, gpx, hx, ihx, hpx, ix, ipx, fgx, fpx, x, x2
775 
776   \textcolor{keywordtype}{integer} :: it, im
777   \textcolor{keywordtype}{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 \textcolor{comment}{!Fixed threshold of 0.01%, could be runtime.}
784 
785   \textcolor{keywordflow}{do} im=1,3
786     pe\_threshold(im) = mixing\_energy(im)/gv%g\_earth
787 \textcolor{keywordflow}{  enddo}
788 
789   \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is,ie
790     \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.0) \textcolor{keywordflow}{then}
791 
792       \textcolor{keyword}{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       \textcolor{keywordflow}{do} k=1,nz
796         dz(k) = h(i,j,k) * gv%H\_to\_Z
797 \textcolor{keywordflow}{      enddo}
798       z\_u(1) = 0.0
799       z\_l(1) = -dz(1)
800       \textcolor{keywordflow}{do} k=2,nz
801         z\_u(k) = z\_l(k-1)
802         z\_l(k) = z\_l(k-1)-dz(k)
803 \textcolor{keywordflow}{      enddo}
804 
805       \textcolor{keywordflow}{do} im=1,3
806 
807         \textcolor{comment}{! 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         \textcolor{keywordflow}{do} k=1,nz
816 
817           \textcolor{comment}{! 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           \textcolor{comment}{! 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           \textcolor{comment}{! The average density assuming all layers including this were mixed}
825           rho\_ml = rhodz\_ml\_tst/h\_ml\_tst
826 
827           \textcolor{comment}{! The PE assuming all layers including this were mixed}
828           \textcolor{comment}{! Note that 0. could be replaced with "Surface", which doesn't have to be 0}
829           \textcolor{comment}{! but 0 is a good reference value.}
830           pe\_mixed\_tst = 0.5*rho\_ml*(0.**2-(0.-h\_ml\_tst)**2)
831 
832           \textcolor{comment}{! Check if we supplied enough energy to mix to this layer}
833           \textcolor{keywordflow}{if} (pe\_mixed\_tst-pe<=pe\_threshold(im)) \textcolor{keywordflow}{then}
834             h\_ml = h\_ml\_tst
835             rhodz\_ml = rhodz\_ml\_tst
836 
837           \textcolor{keywordflow}{else} \textcolor{comment}{! If not, we need to solve where the energy ran out}
838             \textcolor{comment}{! This will be done with a Newton's method iteration:}
839 
840             r1 = rhodz\_ml/h\_ml \textcolor{comment}{! The density of the mixed layer (not including this layer)}
841             d1 = h\_ml \textcolor{comment}{! The thickness of the mixed layer (not including this layer)}
842             r2 = rho\_c(k) \textcolor{comment}{! The density of this layer}
843             d2 = dz(k) \textcolor{comment}{! The thickness of this layer}
844 
845             \textcolor{comment}{! This block could be used to calculate the function coefficients if}
846             \textcolor{comment}{! we don't reference all values to a surface designated as z=0}
847             \textcolor{comment}{! S = Surface}
848             \textcolor{comment}{! Ca  = -(R2)}
849             \textcolor{comment}{! Cb  = -( (R1*D1) + R2*(2.*D1-2.*S) )}
850             \textcolor{comment}{! D   = D1**2. - 2.*D1*S}
851             \textcolor{comment}{! Cc  = -( R1*D1*(2.*D1-2.*S) + R2*D )}
852             \textcolor{comment}{! Cd  = -(R1*D1*D)}
853             \textcolor{comment}{! Ca2 = R2}
854             \textcolor{comment}{! Cb2 = R2*(2*D1-2*S)}
855             \textcolor{comment}{! C   = S**2 + D2**2 + D1**2 - 2*D1*S - 2.*D2*S +2.*D1*D2}
856             \textcolor{comment}{! Cc2 = R2*(D+S**2-C)}
857             \textcolor{comment}{!}
858             \textcolor{comment}{! 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             \textcolor{comment}{! First guess for an iteration using Newton's method}
870             x = dz(k)*0.5
871 
872             it=0
873             \textcolor{keywordflow}{do} \textcolor{keywordflow}{while}(it<10)\textcolor{comment}{!We can iterate up to 10 times}
874               \textcolor{comment}{! We are trying to solve the function:}
875               \textcolor{comment}{! F(x) = G(x)/H(x)+I(x)}
876               \textcolor{comment}{! for where F(x) = PE+PE\_threshold, or equivalently for where}
877               \textcolor{comment}{! F(x) = G(x)/H(x)+I(x) - (PE+PE\_threshold) = 0}
878               \textcolor{comment}{! We also need the derivative of this function for the Newton's method iteration}
879               \textcolor{comment}{! F'(x) = (G'(x)H(x)-G(x)H'(x))/H(x)^2 + I'(x)}
880               \textcolor{comment}{! 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               \textcolor{comment}{! H, its inverse, and its derivative}
884               hx = (d1+x)
885               ihx = 1./hx
886               hpx = 1.
887               \textcolor{comment}{! I and its derivative}
888               ix = 0.5*(ca2*x**2. + cb2*x + cc2)
889               ipx = 0.5*(2.*ca2*x+cb2)
890 
891               \textcolor{comment}{! 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               \textcolor{comment}{! Check if our solution is within the threshold bounds, if not update}
897               \textcolor{comment}{! using Newton's method.  This appears to converge almost always in}
898               \textcolor{comment}{! one step because the function is very close to linear in most applications.}
899               \textcolor{keywordflow}{if} (abs(fgx)>pe\_threshold(im)*pe\_threshold\_fraction) \textcolor{keywordflow}{then}
900                 x2 = x - fgx/fpx
901                 it = it + 1
902                 \textcolor{keywordflow}{if} (x2<0. .or. x2>dz(k)) \textcolor{keywordflow}{then}
903                   \textcolor{comment}{! The iteration seems to be robust, but we need to do something *if*}
904                   \textcolor{comment}{! things go wrong... How should we treat failed iteration?}
905                   \textcolor{comment}{! Present solution: Stop trying to compute and just say we can't mix this layer.}
906                   x=0
907                   \textcolor{keywordflow}{exit}
908                 \textcolor{keywordflow}{else}
909                   x = x2
910 \textcolor{keywordflow}{                endif}
911               \textcolor{keywordflow}{else}
912                 exit\textcolor{comment}{! Quit the iteration}
913 \textcolor{keywordflow}{              endif}
914 \textcolor{keywordflow}{            enddo}
915             h\_ml = h\_ml+x
916             exit\textcolor{comment}{! Quit looping through the column}
917 \textcolor{keywordflow}{          endif}
918 \textcolor{keywordflow}{        enddo}
919         mld(i,j,im) = h\_ml
920 \textcolor{keywordflow}{      enddo}
921     \textcolor{keywordflow}{else}
922       mld(i,j,:) = 0.0
923 \textcolor{keywordflow}{    endif}
924 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
925 
926   \textcolor{keywordflow}{if} (id\_mld(1) > 0) \textcolor{keyword}{call }post\_data(id\_mld(1), mld(:,:,1), diagptr)
927   \textcolor{keywordflow}{if} (id\_mld(2) > 0) \textcolor{keyword}{call }post\_data(id\_mld(2), mld(:,:,2), diagptr)
928   \textcolor{keywordflow}{if} (id\_mld(3) > 0) \textcolor{keyword}{call }post\_data(id\_mld(3), mld(:,:,3), diagptr)
929 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_aa9248acf0b1cf246a79092077bcb0b46}\label{namespacemom__diabatic__aux_aa9248acf0b1cf246a79092077bcb0b46}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!differential\+\_\+diffuse\+\_\+t\+\_\+s@{differential\+\_\+diffuse\+\_\+t\+\_\+s}}
\index{differential\+\_\+diffuse\+\_\+t\+\_\+s@{differential\+\_\+diffuse\+\_\+t\+\_\+s}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{differential\+\_\+diffuse\+\_\+t\+\_\+s()}{differential\_diffuse\_t\_s()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::differential\+\_\+diffuse\+\_\+t\+\_\+s (\begin{DoxyParamCaption}\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{h,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(inout)}]{tv,  }\item[{type(vertvisc\+\_\+type), intent(in)}]{visc,  }\item[{real, intent(in)}]{dt,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV }\end{DoxyParamCaption})}



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


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\
\hline
\mbox{\tt in}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\
\hline
\mbox{\tt in}  & {\em visc} & Structure containing vertical viscosities, bottom boundary layer properies, and related fields.\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}. \\
\hline
\end{DoxyParams}


Definition at line 227 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
227   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: g\textcolor{comment}{    !< The ocean's grid structure}
228   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{   !< The ocean's vertical grid structure}
229   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
230                            \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{    !< Layer thicknesses [H ~> m or kg m-2]}
231   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),   \textcolor{keywordtype}{intent(inout)} :: tv\textcolor{comment}{   !< Structure containing pointers to any}
232 \textcolor{comment}{                                                 !! available thermodynamic fields.}
233   \textcolor{keywordtype}{type}(vertvisc\_type),     \textcolor{keywordtype}{intent(in)}    :: visc\textcolor{comment}{ !< Structure containing vertical viscosities, bottom}
234 \textcolor{comment}{                                                 !! boundary layer properies, and related fields.}
235   \textcolor{keywordtype}{real},                    \textcolor{keywordtype}{intent(in)}    :: dt\textcolor{comment}{   !<  Time increment [T ~> s].}
236 
237   \textcolor{comment}{! local variables}
238   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: &
239     b1\_t, b1\_s, &  \textcolor{comment}{!  Variables used by the tridiagonal solvers of T & S [H ~> m or kg m-2].}
240     d1\_t, d1\_s     \textcolor{comment}{!  Variables used by the tridiagonal solvers [nondim].}
241   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZK\_(G))} :: &
242     c1\_t, c1\_s     \textcolor{comment}{!  Variables used by the tridiagonal solvers [H ~> m or kg m-2].}
243   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZK\_(G)+1)} :: &
244     mix\_t, mix\_s   \textcolor{comment}{!  Mixing distances in both directions across each interface [H ~> m or kg m-2].}
245   \textcolor{keywordtype}{real} :: h\_tr         \textcolor{comment}{! h\_tr is h at tracer points with a tiny thickness}
246                        \textcolor{comment}{! added to ensure positive definiteness [H ~> m or kg m-2].}
247   \textcolor{keywordtype}{real} :: h\_neglect    \textcolor{comment}{! A thickness that is so small it is usually lost}
248                        \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
249   \textcolor{keywordtype}{real} :: i\_h\_int      \textcolor{comment}{! The inverse of the thickness associated with an}
250                        \textcolor{comment}{! interface [H-1 ~> m-1 or m2 kg-1].}
251   \textcolor{keywordtype}{real} :: b\_denom\_t    \textcolor{comment}{! The first term in the denominators for the expressions}
252   \textcolor{keywordtype}{real} :: b\_denom\_s    \textcolor{comment}{! for b1\_T and b1\_S, both [H ~> m or kg m-2].}
253   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:,:)}, \textcolor{keywordtype}{pointer} :: t=>null(), s=>null()
254   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:,:)}, \textcolor{keywordtype}{pointer} :: kd\_t=>null(), kd\_s=>null() \textcolor{comment}{! Diffusivities [Z2 T-1 ~> m2 s-1].}
255   \textcolor{keywordtype}{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   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(tv%T)) \textcolor{keyword}{call }mom\_error(fatal, &
261       \textcolor{stringliteral}{"differential\_diffuse\_T\_S: Called with an unassociated tv%T"})
262   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(tv%S)) \textcolor{keyword}{call }mom\_error(fatal, &
263       \textcolor{stringliteral}{"differential\_diffuse\_T\_S: Called with an unassociated tv%S"})
264   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(visc%Kd\_extra\_T)) \textcolor{keyword}{call }mom\_error(fatal, &
265       \textcolor{stringliteral}{"differential\_diffuse\_T\_S: Called with an unassociated visc%Kd\_extra\_T"})
266   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(visc%Kd\_extra\_S)) \textcolor{keyword}{call }mom\_error(fatal, &
267       \textcolor{stringliteral}{"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 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,h,h\_neglect,dt,Kd\_T,Kd\_S,G,GV,T,S,nz) &}
272 \textcolor{comment}{!$OMP                          private(I\_h\_int,mix\_T,mix\_S,h\_tr,b1\_T,b1\_S, &}
273 \textcolor{comment}{!$OMP                                  d1\_T,d1\_S,c1\_T,c1\_S,b\_denom\_T,b\_denom\_S)}
274   \textcolor{keywordflow}{do} j=js,je
275     \textcolor{keywordflow}{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 \textcolor{keywordflow}{    enddo}
288     \textcolor{keywordflow}{do} k=2,nz-1 ; \textcolor{keywordflow}{do} i=is,ie
289       \textcolor{comment}{! 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 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
308     \textcolor{keywordflow}{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 \textcolor{keywordflow}{    enddo}
319     \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{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 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
323 \textcolor{keywordflow}{  enddo}
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_a3b3416d0ee87cd1a0d54a910fb681e93}\label{namespacemom__diabatic__aux_a3b3416d0ee87cd1a0d54a910fb681e93}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!find\+\_\+uv\+\_\+at\+\_\+h@{find\+\_\+uv\+\_\+at\+\_\+h}}
\index{find\+\_\+uv\+\_\+at\+\_\+h@{find\+\_\+uv\+\_\+at\+\_\+h}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{find\+\_\+uv\+\_\+at\+\_\+h()}{find\_uv\_at\_h()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::find\+\_\+uv\+\_\+at\+\_\+h (\begin{DoxyParamCaption}\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{u,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)), intent(in)}]{v,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out)}]{u\+\_\+h,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out)}]{v\+\_\+h,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in), optional}]{ea,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in), optional}]{eb }\end{DoxyParamCaption})}



This subroutine calculates u\+\_\+h and v\+\_\+h (velocities at thickness points), optionally using the entrainment amounts passed in as arguments. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em u} & The zonal velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em v} & The meridional velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt out}  & {\em u\+\_\+h} & Zonal velocity interpolated to h points \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\
\hline
\mbox{\tt out}  & {\em v\+\_\+h} & Meridional velocity interpolated to h points \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em ea} & The amount of fluid entrained from the layer\\
\hline
\mbox{\tt in}  & {\em eb} & The amount of fluid entrained from the layer \\
\hline
\end{DoxyParams}


Definition at line 438 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
438   \textcolor{keywordtype}{type}(ocean\_grid\_type),     \textcolor{keywordtype}{intent(in)}  :: g\textcolor{comment}{    !< The ocean's grid structure}
439   \textcolor{keywordtype}{type}(verticalgrid\_type),   \textcolor{keywordtype}{intent(in)}  :: gv\textcolor{comment}{   !< The ocean's vertical grid structure}
440   \textcolor{keywordtype}{type}(unit\_scale\_type),     \textcolor{keywordtype}{intent(in)}  :: us\textcolor{comment}{   !< A dimensional unit scaling type}
441   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G))}, &
442                              \textcolor{keywordtype}{intent(in)}  :: u\textcolor{comment}{    !< The zonal velocity [L T-1 ~> m s-1]}
443   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G))}, &
444                              \textcolor{keywordtype}{intent(in)}  :: v\textcolor{comment}{    !< The meridional velocity [L T-1 ~> m s-1]}
445   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
446                              \textcolor{keywordtype}{intent(in)}  :: h\textcolor{comment}{    !< Layer thicknesses [H ~> m or kg m-2]}
447   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
448                              \textcolor{keywordtype}{intent(out)}   :: u\_h\textcolor{comment}{ !< Zonal velocity interpolated to h points [L T-1 ~> m
       s-1].}
449   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
450                              \textcolor{keywordtype}{intent(out)}   :: v\_h\textcolor{comment}{ !< Meridional velocity interpolated to h points [L T-1 ~>
       m s-1].}
451   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
452                      \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: ea\textcolor{comment}{ !< The amount of fluid entrained from the layer}
453 \textcolor{comment}{                                                 !! above within this time step [H ~> m or kg m-2].}
454 \textcolor{comment}{                                                 !! Omitting ea is the same as setting it to 0.}
455   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
456                      \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: eb\textcolor{comment}{ !< The amount of fluid entrained from the layer}
457 \textcolor{comment}{                                                 !! below within this time step [H ~> m or kg m-2].}
458 \textcolor{comment}{                                                 !! Omitting eb is the same as setting it to 0.}
459 
460   \textcolor{comment}{! local variables}
461   \textcolor{keywordtype}{real} :: b\_denom\_1    \textcolor{comment}{! The first term in the denominator of b1 [H ~> m or kg m-2].}
462   \textcolor{keywordtype}{real} :: h\_neglect    \textcolor{comment}{! A thickness that is so small it is usually lost}
463                        \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
464   \textcolor{keywordtype}{real} :: b1(szi\_(g)), d1(szi\_(g)), c1(szi\_(g),szk\_(g))
465   \textcolor{keywordtype}{real} :: a\_n(szi\_(g)), a\_s(szi\_(g))  \textcolor{comment}{! Fractional weights of the neighboring}
466   \textcolor{keywordtype}{real} :: a\_e(szi\_(g)), a\_w(szi\_(g))  \textcolor{comment}{! velocity points, ~1/2 in the open}
467                                       \textcolor{comment}{! ocean, nondimensional.}
468   \textcolor{keywordtype}{real} :: sum\_area, idenom
469   \textcolor{keywordtype}{logical} :: mix\_vertically
470   \textcolor{keywordtype}{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   \textcolor{keyword}{call }cpu\_clock\_begin(id\_clock\_uv\_at\_h)
473   h\_neglect = gv%H\_subroundoff
474 
475   mix\_vertically = \textcolor{keyword}{present}(ea)
476   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(ea) .neqv. \textcolor{keyword}{present}(eb)) \textcolor{keyword}{call }mom\_error(fatal, &
477       \textcolor{stringliteral}{"find\_uv\_at\_h: Either both ea and eb or neither one must be present "}// &
478       \textcolor{stringliteral}{"in call to find\_uv\_at\_h."})
479 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix\_vertically,h,h\_neglect, &}
480 \textcolor{comment}{!$OMP                                  eb,u\_h,u,v\_h,v,nz,ea)                     &}
481 \textcolor{comment}{!$OMP                          private(sum\_area,Idenom,a\_w,a\_e,a\_s,a\_n,b\_denom\_1,b1,d1,c1)}
482   \textcolor{keywordflow}{do} j=js,je
483     \textcolor{keywordflow}{do} i=is,ie
484       sum\_area = g%areaCu(i-1,j) + g%areaCu(i,j)
485       \textcolor{keywordflow}{if} (sum\_area>0.0) \textcolor{keywordflow}{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       \textcolor{keywordflow}{else}
490         a\_w(i) = 0.0 ; a\_e(i) = 0.0
491 \textcolor{keywordflow}{      endif}
492 
493       sum\_area = g%areaCv(i,j-1) + g%areaCv(i,j)
494       \textcolor{keywordflow}{if} (sum\_area>0.0) \textcolor{keywordflow}{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       \textcolor{keywordflow}{else}
499         a\_s(i) = 0.0 ; a\_n(i) = 0.0
500 \textcolor{keywordflow}{      endif}
501 \textcolor{keywordflow}{    enddo}
502 
503     \textcolor{keywordflow}{if} (mix\_vertically) \textcolor{keywordflow}{then}
504       \textcolor{keywordflow}{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 \textcolor{keywordflow}{      enddo}
511       \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{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 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
521       \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{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 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
525     \textcolor{keywordflow}{else}
526       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{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 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
530 \textcolor{keywordflow}{    endif}
531 \textcolor{keywordflow}{  enddo}
532 
533   \textcolor{keyword}{call }cpu\_clock\_end(id\_clock\_uv\_at\_h)
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_adb99fd4e2092c17cd69143945733a6d9}\label{namespacemom__diabatic__aux_adb99fd4e2092c17cd69143945733a6d9}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!make\+\_\+frazil@{make\+\_\+frazil}}
\index{make\+\_\+frazil@{make\+\_\+frazil}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{make\+\_\+frazil()}{make\_frazil()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::make\+\_\+frazil (\begin{DoxyParamCaption}\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(inout)}]{tv,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), intent(in)}]{CS,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g)), intent(in), optional}]{p\+\_\+surf,  }\item[{integer, intent(in), optional}]{halo }\end{DoxyParamCaption})}



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 \mbox{[}Q R Z $\sim$$>$ J m-\/2\mbox{]}) in tvfrazil. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\
\hline
\mbox{\tt in}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em cs} & The control structure returned by a previous call to diabatic\+\_\+aux\+\_\+init.\\
\hline
\mbox{\tt in}  & {\em p\+\_\+surf} & The pressure at the ocean surface \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em halo} & Halo width over which to calculate frazil \\
\hline
\end{DoxyParams}


Definition at line 104 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
104   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: g\textcolor{comment}{  !< The ocean's grid structure}
105   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{ !< The ocean's vertical grid structure}
106   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
107                            \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{  !< Layer thicknesses [H ~> m or kg m-2]}
108   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),   \textcolor{keywordtype}{intent(inout)} :: tv\textcolor{comment}{ !< Structure containing pointers to any available}
109 \textcolor{comment}{                                               !! thermodynamic fields.}
110   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)}    :: us\textcolor{comment}{   !< A dimensional unit scaling type}
111   \textcolor{keywordtype}{type}(diabatic\_aux\_cs),   \textcolor{keywordtype}{intent(in)}    :: cs\textcolor{comment}{ !< The control structure returned by a previous}
112 \textcolor{comment}{                                               !! call to diabatic\_aux\_init.}
113   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, &
114                  \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: p\_surf\textcolor{comment}{ !< The pressure at the ocean surface [R L2 T-2 ~> Pa].}
115   \textcolor{keywordtype}{integer},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: halo\textcolor{comment}{ !< Halo width over which to calculate frazil}
116 
117   \textcolor{comment}{! Local variables}
118   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: &
119     fraz\_col, & \textcolor{comment}{! The accumulated heat requirement due to frazil [Q R Z ~> J m-2].}
120     t\_freeze, & \textcolor{comment}{! The freezing potential temperature at the current salinity [degC].}
121     ps          \textcolor{comment}{! Surface pressure [R L2 T-2 ~> Pa]}
122   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZK\_(G))} :: &
123     pressure    \textcolor{comment}{! The pressure at the middle of each layer [R L2 T-2 ~> Pa].}
124   \textcolor{keywordtype}{real} :: h\_to\_rl2\_t2  \textcolor{comment}{! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or
       Pa m2 kg-1]}
125   \textcolor{keywordtype}{real} :: hc    \textcolor{comment}{! A layer's heat capacity [Q R Z degC-1 ~> J m-2 degC-1].}
126   \textcolor{keywordtype}{logical} :: t\_fr\_set  \textcolor{comment}{! True if the freezing point has been calculated for a}
127                        \textcolor{comment}{! row of points.}
128   \textcolor{keywordtype}{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   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(halo)) \textcolor{keywordflow}{then}
132     is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
133 \textcolor{keywordflow}{  endif}
134 
135   \textcolor{keyword}{call }cpu\_clock\_begin(id\_clock\_frazil)
136 
137   \textcolor{keywordflow}{if} (.not.cs%pressure\_dependent\_frazil) \textcolor{keywordflow}{then}
138     \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; pressure(i,k) = 0.0 ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
139   \textcolor{keywordflow}{else}
140     h\_to\_rl2\_t2 = gv%H\_to\_RZ * gv%g\_Earth
141 \textcolor{keywordflow}{  endif}
142   \textcolor{comment}{!$OMP parallel do default(shared) private(fraz\_col,T\_fr\_set,T\_freeze,hc,ps)  &}
143   \textcolor{comment}{!$OMP                             firstprivate(pressure) ! pressure might be set above, so should be
       firstprivate}
144   \textcolor{keywordflow}{do} j=js,je
145     ps(:) = 0.0
146     \textcolor{keywordflow}{if} (\textcolor{keyword}{PRESENT}(p\_surf)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=is,ie
147       ps(i) = p\_surf(i,j)
148 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ endif}
149 
150     \textcolor{keywordflow}{do} i=is,ie ; fraz\_col(i) = 0.0 ;\textcolor{keywordflow}{ enddo}
151 
152     \textcolor{keywordflow}{if} (cs%pressure\_dependent\_frazil) \textcolor{keywordflow}{then}
153       \textcolor{keywordflow}{do} i=is,ie
154         pressure(i,1) = ps(i) + (0.5*h\_to\_rl2\_t2)*h(i,j,1)
155 \textcolor{keywordflow}{      enddo}
156       \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{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 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
160 \textcolor{keywordflow}{    endif}
161 
162     \textcolor{keywordflow}{if} (cs%reclaim\_frazil) \textcolor{keywordflow}{then}
163       t\_fr\_set = .false.
164       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (tv%frazil(i,j) > 0.0) \textcolor{keywordflow}{then}
165         \textcolor{keywordflow}{if} (.not.t\_fr\_set) \textcolor{keywordflow}{then}
166           \textcolor{keyword}{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 \textcolor{keywordflow}{        endif}
170 
171         \textcolor{keywordflow}{if} (tv%T(i,j,1) > t\_freeze(i)) \textcolor{keywordflow}{then}
172     \textcolor{comment}{! If frazil had previously been formed, but the surface temperature is now}
173     \textcolor{comment}{! 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           \textcolor{keywordflow}{if} (tv%frazil(i,j) - hc * (tv%T(i,j,1) - t\_freeze(i)) <= 0.0) \textcolor{keywordflow}{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           \textcolor{keywordflow}{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 \textcolor{keywordflow}{          endif}
182 \textcolor{keywordflow}{        endif}
183 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
184 \textcolor{keywordflow}{    endif}
185 
186     \textcolor{keywordflow}{do} k=nz,1,-1
187       t\_fr\_set = .false.
188       \textcolor{keywordflow}{do} i=is,ie
189         \textcolor{keywordflow}{if} ((g%mask2dT(i,j) > 0.0) .and. &
190             ((tv%T(i,j,k) < 0.0) .or. (fraz\_col(i) > 0.0))) \textcolor{keywordflow}{then}
191           \textcolor{keywordflow}{if} (.not.t\_fr\_set) \textcolor{keywordflow}{then}
192             \textcolor{keyword}{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 \textcolor{keywordflow}{          endif}
196 
197           hc = (tv%C\_p*gv%H\_to\_RZ) * h(i,j,k)
198           \textcolor{keywordflow}{if} (h(i,j,k) <= 10.0*(gv%Angstrom\_H + gv%H\_subroundoff)) \textcolor{keywordflow}{then}
199             \textcolor{comment}{! Very thin layers should not be cooled by the frazil flux.}
200             \textcolor{keywordflow}{if} (tv%T(i,j,k) < t\_freeze(i)) \textcolor{keywordflow}{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 \textcolor{keywordflow}{            endif}
204           \textcolor{keywordflow}{elseif} ((fraz\_col(i) > 0.0) .or. (tv%T(i,j,k) < t\_freeze(i))) \textcolor{keywordflow}{then}
205             \textcolor{keywordflow}{if} (fraz\_col(i) + hc * (t\_freeze(i) - tv%T(i,j,k)) < 0.0) \textcolor{keywordflow}{then}
206               tv%T(i,j,k) = tv%T(i,j,k) - fraz\_col(i) / hc
207               fraz\_col(i) = 0.0
208             \textcolor{keywordflow}{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 \textcolor{keywordflow}{            endif}
212 \textcolor{keywordflow}{          endif}
213 \textcolor{keywordflow}{        endif}
214 \textcolor{keywordflow}{      enddo}
215 \textcolor{keywordflow}{    enddo}
216     \textcolor{keywordflow}{do} i=is,ie
217       tv%frazil(i,j) = tv%frazil(i,j) + fraz\_col(i)
218 \textcolor{keywordflow}{    enddo}
219 \textcolor{keywordflow}{  enddo}
220   \textcolor{keyword}{call }cpu\_clock\_end(id\_clock\_frazil)
221 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_adce369573ff355b820320d9c3b048693}\label{namespacemom__diabatic__aux_adce369573ff355b820320d9c3b048693}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!set\+\_\+pen\+\_\+shortwave@{set\+\_\+pen\+\_\+shortwave}}
\index{set\+\_\+pen\+\_\+shortwave@{set\+\_\+pen\+\_\+shortwave}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{set\+\_\+pen\+\_\+shortwave()}{set\_pen\_shortwave()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::set\+\_\+pen\+\_\+shortwave (\begin{DoxyParamCaption}\item[{type(optics\+\_\+type), pointer}]{optics,  }\item[{type(forcing), intent(inout)}]{fluxes,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), pointer}]{CS,  }\item[{type(opacity\+\_\+cs), pointer}]{opacity\+\_\+\+C\+Sp,  }\item[{type(tracer\+\_\+flow\+\_\+control\+\_\+cs), pointer}]{tracer\+\_\+flow\+\_\+\+C\+Sp }\end{DoxyParamCaption})}


\begin{DoxyParams}[1]{Parameters}
 & {\em optics} & An optics structure that has will contain information about shortwave fluxes and absorption.\\
\hline
\mbox{\tt in,out}  & {\em fluxes} & points to forcing fields unused fields have N\+U\+LL ptrs\\
\hline
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure.\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure.\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
 & {\em cs} & Control structure for diabatic\+\_\+aux\\
\hline
 & {\em opacity\+\_\+csp} & The control structure for the opacity module.\\
\hline
 & {\em tracer\+\_\+flow\+\_\+csp} & A pointer to the control structure organizing the tracer modules. \\
\hline
\end{DoxyParams}


Definition at line 538 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
538   \textcolor{keywordtype}{type}(optics\_type),       \textcolor{keywordtype}{pointer}       :: optics\textcolor{comment}{ !< An optics structure that has will contain}
539 \textcolor{comment}{                                                   !! information about shortwave fluxes and absorption.}
540   \textcolor{keywordtype}{type}(forcing),           \textcolor{keywordtype}{intent(inout)} :: fluxes\textcolor{comment}{ !< points to forcing fields}
541 \textcolor{comment}{                                                   !! unused fields have NULL ptrs}
542   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: g\textcolor{comment}{      !< The ocean's grid structure.}
543   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{     !< The ocean's vertical grid structure.}
544   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)}    :: us\textcolor{comment}{     !< A dimensional unit scaling type}
545   \textcolor{keywordtype}{type}(diabatic\_aux\_cs),   \textcolor{keywordtype}{pointer}       :: cs\textcolor{comment}{     !< Control structure for diabatic\_aux}
546   \textcolor{keywordtype}{type}(opacity\_cs),        \textcolor{keywordtype}{pointer}       :: opacity\_csp\textcolor{comment}{ !< The control structure for the opacity module.}
547   \textcolor{keywordtype}{type}(tracer\_flow\_control\_cs), \textcolor{keywordtype}{pointer}  :: tracer\_flow\_csp\textcolor{comment}{ !< A pointer to the control structure}
548 \textcolor{comment}{                                                   !! organizing the tracer modules.}
549 
550   \textcolor{comment}{! Local variables}
551   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}          :: chl\_2d\textcolor{comment}{ !< Vertically uniform chlorophyll-A concentractions
       [mg m-3]}
552   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(GV))} :: chl\_3d\textcolor{comment}{ !< The chlorophyll-A concentractions of each layer
       [mg m-3]}
553   \textcolor{keywordtype}{character(len=128)} :: mesg
554   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je
555   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
556 
557   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(optics)) \textcolor{keywordflow}{return}
558 
559   \textcolor{keywordflow}{if} (cs%var\_pen\_sw) \textcolor{keywordflow}{then}
560     \textcolor{keywordflow}{if} (cs%chl\_from\_file) \textcolor{keywordflow}{then}
561       \textcolor{comment}{! Only the 2-d surface chlorophyll can be read in from a file.  The}
562       \textcolor{comment}{! same value is assumed for all layers.}
563       \textcolor{keyword}{call }time\_interp\_external(cs%sbc\_chl, cs%Time, chl\_2d)
564       \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
565         \textcolor{keywordflow}{if} ((g%mask2dT(i,j) > 0.5) .and. (chl\_2d(i,j) < 0.0)) \textcolor{keywordflow}{then}
566           \textcolor{keyword}{write}(mesg,\textcolor{stringliteral}{'(" Time\_interp negative chl of ",(1pe12.4)," at i,j = ",&}
567 \textcolor{stringliteral}{}\textcolor{stringliteral}{                    & 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           \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_diabatic\_aux set\_pen\_shortwave: "}//trim(mesg))
570 \textcolor{keywordflow}{        endif}
571 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
572 
573       \textcolor{keywordflow}{if} (cs%id\_chl > 0) \textcolor{keyword}{call }post\_data(cs%id\_chl, chl\_2d, cs%diag)
574 
575       \textcolor{keyword}{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     \textcolor{keywordflow}{else}
578       \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(tracer\_flow\_csp)) \textcolor{keyword}{call }mom\_error(fatal, &
579         \textcolor{stringliteral}{"The tracer flow control structure must be associated when the model sets "}//&
580         \textcolor{stringliteral}{"the chlorophyll internally in set\_pen\_shortwave."})
581       \textcolor{keyword}{call }get\_chl\_from\_model(chl\_3d, g, tracer\_flow\_csp)
582 
583       \textcolor{keywordflow}{if} (cs%id\_chl > 0) \textcolor{keyword}{call }post\_data(cs%id\_chl, chl\_3d(:,:,1), cs%diag)
584 
585       \textcolor{keyword}{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 \textcolor{keywordflow}{    endif}
588   \textcolor{keywordflow}{else}
589     \textcolor{keyword}{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 \textcolor{keywordflow}{  endif}
592 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_acde4c76c9a489b82cba0bac5ec1b1198}\label{namespacemom__diabatic__aux_acde4c76c9a489b82cba0bac5ec1b1198}} 
\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!tridiagts@{tridiagts}}
\index{tridiagts@{tridiagts}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}}
\subsubsection{\texorpdfstring{tridiagts()}{tridiagts()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::tridiagts (\begin{DoxyParamCaption}\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{integer, intent(in)}]{is,  }\item[{integer, intent(in)}]{ie,  }\item[{integer, intent(in)}]{js,  }\item[{integer, intent(in)}]{je,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{hold,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{ea,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{eb,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{T,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{S }\end{DoxyParamCaption})}



This is a simple tri-\/diagonal solver for T and S. \char`\"{}\+Simple\char`\"{} means it only uses arrays hold, ea and eb. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\
\hline
\mbox{\tt in}  & {\em is} & The start i-\/index to work on.\\
\hline
\mbox{\tt in}  & {\em ie} & The end i-\/index to work on.\\
\hline
\mbox{\tt in}  & {\em js} & The start j-\/index to work on.\\
\hline
\mbox{\tt in}  & {\em je} & The end j-\/index to work on.\\
\hline
\mbox{\tt in}  & {\em hold} & The layer thicknesses before entrainment, \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em ea} & The amount of fluid entrained from the layer above within this time step \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em eb} & The amount of fluid entrained from the layer below within this time step \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em t} & Layer potential temperatures \mbox{[}degC\mbox{]}.\\
\hline
\mbox{\tt in,out}  & {\em s} & Layer salinities \mbox{[}ppt\mbox{]}. \\
\hline
\end{DoxyParams}


Definition at line 389 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
389   \textcolor{keywordtype}{type}(ocean\_grid\_type),                    \textcolor{keywordtype}{intent(in)}    :: g\textcolor{comment}{    !< The ocean's grid structure}
390   \textcolor{keywordtype}{type}(verticalgrid\_type),                  \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{   !< The ocean's vertical grid structure}
391   \textcolor{keywordtype}{integer},                                  \textcolor{keywordtype}{intent(in)}    :: is\textcolor{comment}{   !< The start i-index to work on.}
392   \textcolor{keywordtype}{integer},                                  \textcolor{keywordtype}{intent(in)}    :: ie\textcolor{comment}{   !< The end i-index to work on.}
393   \textcolor{keywordtype}{integer},                                  \textcolor{keywordtype}{intent(in)}    :: js\textcolor{comment}{   !< The start j-index to work on.}
394   \textcolor{keywordtype}{integer},                                  \textcolor{keywordtype}{intent(in)}    :: je\textcolor{comment}{   !< The end j-index to work on.}
395   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(in)}    :: hold\textcolor{comment}{ !< The layer thicknesses before
       entrainment,}
396 \textcolor{comment}{                                                                  !! [H ~> m or kg m-2].}
397   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(in)}    :: ea\textcolor{comment}{ !< The amount of fluid entrained from the
       layer}
398 \textcolor{comment}{                                                                !! above within this time step [H ~> m or
       kg m-2]}
399   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(in)}    :: eb\textcolor{comment}{ !< The amount of fluid entrained from the
       layer}
400 \textcolor{comment}{                                                                !! below within this time step [H ~> m or
       kg m-2]}
401   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(inout)} :: t\textcolor{comment}{  !< Layer potential temperatures [degC].}
402   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(inout)} :: s\textcolor{comment}{  !< Layer salinities [ppt].}
403 
404   \textcolor{comment}{! Local variables}
405   \textcolor{keywordtype}{real} :: b1(szib\_(g)), d1(szib\_(g)) \textcolor{comment}{! b1, c1, and d1 are variables used by the}
406   \textcolor{keywordtype}{real} :: c1(szib\_(g),szk\_(g))       \textcolor{comment}{! tridiagonal solver.}
407   \textcolor{keywordtype}{real} :: h\_tr, b\_denom\_1
408   \textcolor{keywordtype}{integer} :: i, j, k
409 
410   \textcolor{comment}{!$OMP parallel do default(shared) private(h\_tr,b1,d1,c1,b\_denom\_1)}
411   \textcolor{keywordflow}{do} j=js,je
412     \textcolor{keywordflow}{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 \textcolor{keywordflow}{    enddo}
419     \textcolor{keywordflow}{do} k=2,g%ke ; \textcolor{keywordflow}{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 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
428     \textcolor{keywordflow}{do} k=g%ke-1,1,-1 ; \textcolor{keywordflow}{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 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
432 \textcolor{keywordflow}{  enddo}
\end{DoxyCode}
