\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_a0b81bac82b18bd9f1e3e5a3beda78966}{differential\+\_\+diffuse\+\_\+t\+\_\+s} (h, T, S, Kd\+\_\+T, Kd\+\_\+S, 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 323 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
323   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: g\textcolor{comment}{    !< The ocean's grid structure}
324   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{   !< The ocean's vertical grid structure}
325   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
326                            \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{    !< Layer thicknesses [H ~> m or kg m-2]}
327   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),   \textcolor{keywordtype}{intent(inout)} :: tv\textcolor{comment}{   !< Structure containing pointers to any}
328 \textcolor{comment}{                                                 !! available thermodynamic fields.}
329   \textcolor{keywordtype}{type}(diabatic\_aux\_cs),   \textcolor{keywordtype}{intent(in)}    :: cs\textcolor{comment}{   !< The control structure returned by a previous}
330 \textcolor{comment}{                                                 !! call to diabatic\_aux\_init.}
331   \textcolor{keywordtype}{integer},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: halo\textcolor{comment}{ !< Halo width over which to work}
332 
333   \textcolor{comment}{! local variables}
334   \textcolor{keywordtype}{real} :: salt\_add\_col(szi\_(g),szj\_(g))\textcolor{comment}{ !< The accumulated salt requirement [ppt R Z ~> gSalt m-2]}
335   \textcolor{keywordtype}{real} :: s\_min\textcolor{comment}{      !< The minimum salinity [ppt].}
336   \textcolor{keywordtype}{real} :: mc\textcolor{comment}{         !< A layer's mass [R Z ~> kg m-2].}
337   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, nz
338 
339   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
340   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(halo)) \textcolor{keywordflow}{then}
341     is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
342 \textcolor{keywordflow}{  endif}
343 
344 \textcolor{comment}{!  call cpu\_clock\_begin(id\_clock\_adjust\_salt)}
345 
346   s\_min = tv%min\_salinity
347 
348   salt\_add\_col(:,:) = 0.0
349 
350   \textcolor{comment}{!$OMP parallel do default(shared) private(mc)}
351   \textcolor{keywordflow}{do} j=js,je
352     \textcolor{keywordflow}{do} k=nz,1,-1 ; \textcolor{keywordflow}{do} i=is,ie
353       \textcolor{keywordflow}{if} ( (g%mask2dT(i,j) > 0.0) .and. &
354            ((tv%S(i,j,k) < s\_min) .or. (salt\_add\_col(i,j) > 0.0)) ) \textcolor{keywordflow}{then}
355         mc = gv%H\_to\_RZ * h(i,j,k)
356         \textcolor{keywordflow}{if} (h(i,j,k) <= 10.0*gv%Angstrom\_H) \textcolor{keywordflow}{then}
357           \textcolor{comment}{! Very thin layers should not be adjusted by the salt flux}
358           \textcolor{keywordflow}{if} (tv%S(i,j,k) < s\_min) \textcolor{keywordflow}{then}
359             salt\_add\_col(i,j) = salt\_add\_col(i,j) +  mc * (s\_min - tv%S(i,j,k))
360             tv%S(i,j,k) = s\_min
361 \textcolor{keywordflow}{          endif}
362         \textcolor{keywordflow}{elseif} (salt\_add\_col(i,j) + mc * (s\_min - tv%S(i,j,k)) <= 0.0) \textcolor{keywordflow}{then}
363           tv%S(i,j,k) = tv%S(i,j,k) - salt\_add\_col(i,j) / mc
364           salt\_add\_col(i,j) = 0.0
365         \textcolor{keywordflow}{else}
366           salt\_add\_col(i,j) = salt\_add\_col(i,j) + mc * (s\_min - tv%S(i,j,k))
367           tv%S(i,j,k) = s\_min
368 \textcolor{keywordflow}{        endif}
369 \textcolor{keywordflow}{      endif}
370 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
371     \textcolor{keywordflow}{do} i=is,ie
372       tv%salt\_deficit(i,j) = tv%salt\_deficit(i,j) + salt\_add\_col(i,j)
373 \textcolor{keywordflow}{    enddo}
374 \textcolor{keywordflow}{  enddo}
375 \textcolor{comment}{!  call cpu\_clock\_end(id\_clock\_adjust\_salt)}
376 
\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 932 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


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


\begin{DoxyCode}
1645   \textcolor{keywordtype}{type}(diabatic\_aux\_cs), \textcolor{keywordtype}{pointer} :: cs\textcolor{comment}{ !< The control structure returned by a previous}
1646 \textcolor{comment}{                                       !! call to diabatic\_aux\_init; it is deallocated here.}
1647 
1648   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{return}
1649 
1650   \textcolor{keywordflow}{if} (cs%id\_createdH       >0) \textcolor{keyword}{deallocate}(cs%createdH)
1651   \textcolor{keywordflow}{if} (cs%id\_penSW\_diag     >0) \textcolor{keyword}{deallocate}(cs%penSW\_diag)
1652   \textcolor{keywordflow}{if} (cs%id\_penSWflux\_diag >0) \textcolor{keyword}{deallocate}(cs%penSWflux\_diag)
1653   \textcolor{keywordflow}{if} (cs%id\_nonpenSW\_diag  >0) \textcolor{keyword}{deallocate}(cs%nonpenSW\_diag)
1654 
1655   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{deallocate}(cs)
1656 
\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 1481 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


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


\begin{DoxyCode}
593   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{           !< Grid type}
594   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{          !< ocean vertical grid structure}
595   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{          !< A dimensional unit scaling type}
596   \textcolor{keywordtype}{integer},                 \textcolor{keywordtype}{intent(in)} :: id\_mld\textcolor{comment}{      !< Handle (ID) of MLD diagnostic}
597   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
598                            \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{           !< Layer thickness [H ~> m or kg m-2]}
599   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),   \textcolor{keywordtype}{intent(in)} :: tv\textcolor{comment}{          !< Structure containing pointers to any}
600 \textcolor{comment}{                                                     !! available thermodynamic fields.}
601   \textcolor{keywordtype}{real},                    \textcolor{keywordtype}{intent(in)} :: densitydiff\textcolor{comment}{ !< Density difference to determine MLD [R ~> kg m-3]}
602   \textcolor{keywordtype}{type}(diag\_ctrl),         \textcolor{keywordtype}{pointer}    :: diagptr\textcolor{comment}{     !< Diagnostics structure}
603   \textcolor{keywordtype}{integer},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: id\_n2subml\textcolor{comment}{  !< Optional handle (ID) of subML stratification}
604   \textcolor{keywordtype}{integer},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: id\_mldsq\textcolor{comment}{    !< Optional handle (ID) of squared MLD}
605   \textcolor{keywordtype}{real},          \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: dz\_subml\textcolor{comment}{    !< The distance over which to calculate N2subML}
606 \textcolor{comment}{                                                     !! or 50 m if missing [Z ~> m]}
607 
608   \textcolor{comment}{! Local variables}
609   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: deltarhoatkm1, deltarhoatk \textcolor{comment}{! Density differences [R ~> kg m-3].}
610   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: pref\_mld, pref\_n2 \textcolor{comment}{! Reference pressures [R L2 T-2 ~> Pa].}
611   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: h\_subml, dh\_n2    \textcolor{comment}{! Summed thicknesses used in N2 calculation [H ~> m].}
612   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: t\_subml, t\_deeper \textcolor{comment}{! Temperatures used in the N2 calculation [degC].}
613   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: s\_subml, s\_deeper \textcolor{comment}{! Salinities used in the N2 calculation [ppt].}
614   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: rho\_subml, rho\_deeper \textcolor{comment}{! Densities used in the N2 calculation [R ~> kg m-3].}
615   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: dk, dkm1          \textcolor{comment}{! Depths [Z ~> m].}
616   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: rhosurf          \textcolor{comment}{! Density used in finding the mixedlayer depth [R ~> kg
       m-3].}
617   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G))} :: mld     \textcolor{comment}{! Diagnosed mixed layer depth [Z ~> m].}
618   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G))} :: submln2 \textcolor{comment}{! Diagnosed stratification below ML [T-2 ~> s-2].}
619   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G))} :: mld2    \textcolor{comment}{! Diagnosed MLD^2 [Z2 ~> m2].}
620   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: n2\_region\_set \textcolor{comment}{! If true, all necessary values for calculating N2}
621                                                \textcolor{comment}{! have been stored already.}
622   \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].}
623   \textcolor{keywordtype}{real} :: dh\_subml         \textcolor{comment}{! Depth below ML over which to diagnose stratification [H ~> m].}
624   \textcolor{keywordtype}{real} :: afac             \textcolor{comment}{! A nondimensional factor [nondim]}
625   \textcolor{keywordtype}{real} :: ddrho            \textcolor{comment}{! A density difference [R ~> kg m-3]}
626   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(2)} :: eosdom \textcolor{comment}{! The i-computational domain for the equation of state}
627   \textcolor{keywordtype}{integer} :: i, j, is, ie, js, je, k, nz, id\_n2, id\_sq
628 
629   id\_n2 = -1 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{PRESENT}(id\_n2subml)) id\_n2 = id\_n2subml
630 
631   id\_sq = -1 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{PRESENT}(id\_mldsq)) id\_sq = id\_mldsq
632 
633   ge\_rho0 = us%L\_to\_Z**2*gv%g\_Earth / (gv%Rho0)
634   dh\_subml = 50.*gv%m\_to\_H  ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(dz\_subml)) dh\_subml = gv%Z\_to\_H*dz\_subml
635 
636   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
637 
638   pref\_mld(:) = 0.0
639   eosdom(:) = eos\_domain(g%HI)
640   \textcolor{keywordflow}{do} j=js,je
641     \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}
642     \textcolor{keyword}{call }calculate\_density(tv%T(:,j,1), tv%S(:,j,1), pref\_mld, rhosurf, tv%eqn\_of\_state, eosdom)
643     \textcolor{keywordflow}{do} i=is,ie
644       deltarhoatk(i) = 0.
645       mld(i,j) = 0.
646       \textcolor{keywordflow}{if} (id\_n2>0) \textcolor{keywordflow}{then}
647         submln2(i,j) = 0.0
648         h\_subml(i) = h(i,j,1) ; dh\_n2(i) = 0.0
649         t\_subml(i) = 0.0  ; s\_subml(i) = 0.0 ; t\_deeper(i) = 0.0 ; s\_deeper(i) = 0.0
650         n2\_region\_set(i) = (g%mask2dT(i,j)<0.5) \textcolor{comment}{! Only need to work on ocean points.}
651 \textcolor{keywordflow}{      endif}
652 \textcolor{keywordflow}{    enddo}
653     \textcolor{keywordflow}{do} k=2,nz
654       \textcolor{keywordflow}{do} i=is,ie
655         dkm1(i) = dk(i) \textcolor{comment}{! Depth of center of layer K-1}
656         dk(i) = dk(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * gv%H\_to\_Z \textcolor{comment}{! Depth of center of layer K}
657 \textcolor{keywordflow}{      enddo}
658 
659       \textcolor{comment}{! Prepare to calculate stratification, N2, immediately below the mixed layer by finding}
660       \textcolor{comment}{! the cells that extend over at least dz\_subML.}
661       \textcolor{keywordflow}{if} (id\_n2>0) \textcolor{keywordflow}{then}
662          \textcolor{keywordflow}{do} i=is,ie
663           \textcolor{keywordflow}{if} (mld(i,j)==0.0) \textcolor{keywordflow}{then}  \textcolor{comment}{! Still in the mixed layer.}
664             h\_subml(i) = h\_subml(i) + h(i,j,k)
665           \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.}
666             \textcolor{keywordflow}{if} (dh\_n2(i)==0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! Record the temperature, salinity, pressure, immediately below the ML}
667               t\_subml(i) = tv%T(i,j,k) ; s\_subml(i) = tv%S(i,j,k)
668               h\_subml(i) = h\_subml(i) + 0.5 * h(i,j,k) \textcolor{comment}{! Start midway through this layer.}
669               dh\_n2(i) = 0.5 * h(i,j,k)
670             \textcolor{keywordflow}{elseif} (dh\_n2(i) + h(i,j,k) < dh\_subml) \textcolor{keywordflow}{then}
671               dh\_n2(i) = dh\_n2(i) + h(i,j,k)
672             \textcolor{keywordflow}{else}  \textcolor{comment}{! This layer includes the base of the region where N2 is calculated.}
673               t\_deeper(i) = tv%T(i,j,k) ; s\_deeper(i) = tv%S(i,j,k)
674               dh\_n2(i) = dh\_n2(i) + 0.5 * h(i,j,k)
675               n2\_region\_set(i) = .true.
676 \textcolor{keywordflow}{            endif}
677 \textcolor{keywordflow}{          endif}
678 \textcolor{keywordflow}{        enddo} \textcolor{comment}{! i-loop}
679 \textcolor{keywordflow}{      endif} \textcolor{comment}{! id\_N2>0}
680 
681       \textcolor{comment}{! Mixed-layer depth, using sigma-0 (surface reference pressure)}
682       \textcolor{keywordflow}{do} i=is,ie ; deltarhoatkm1(i) = deltarhoatk(i) ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! Store value from previous iteration of K}
683       \textcolor{keyword}{call }calculate\_density(tv%T(:,j,k), tv%S(:,j,k), pref\_mld, deltarhoatk, tv%eqn\_of\_state, eosdom)
684       \textcolor{keywordflow}{do} i = is, ie
685         deltarhoatk(i) = deltarhoatk(i) - rhosurf(i) \textcolor{comment}{! Density difference between layer K and surface}
686         ddrho = deltarhoatk(i) - deltarhoatkm1(i)
687         \textcolor{keywordflow}{if} ((mld(i,j)==0.) .and. (ddrho>0.) .and. &
688             (deltarhoatkm1(i)<densitydiff) .and. (deltarhoatk(i)>=densitydiff)) \textcolor{keywordflow}{then}
689           afac = ( densitydiff - deltarhoatkm1(i) ) / ddrho
690           mld(i,j) = dk(i) * afac + dkm1(i) * (1. - afac)
691 \textcolor{keywordflow}{        endif}
692         \textcolor{keywordflow}{if} (id\_sq > 0) mld2(i,j) = mld(i,j)**2
693 \textcolor{keywordflow}{      enddo} \textcolor{comment}{! i-loop}
694 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! k-loop}
695     \textcolor{keywordflow}{do} i=is,ie
696       \textcolor{keywordflow}{if} ((mld(i,j)==0.) .and. (deltarhoatk(i)<densitydiff)) mld(i,j) = dk(i) \textcolor{comment}{! Assume mixing to the bottom}
697 \textcolor{keywordflow}{    enddo}
698 
699     \textcolor{keywordflow}{if} (id\_n2>0) \textcolor{keywordflow}{then}  \textcolor{comment}{! Now actually calculate stratification, N2, below the mixed layer.}
700       \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}
701       \textcolor{comment}{! if ((.not.N2\_region\_set(i)) .and. (dH\_N2(i) > 0.5*dH\_subML)) then}
702       \textcolor{comment}{!    ! Use whatever stratification we can, measured over whatever distance is available?}
703       \textcolor{comment}{!    T\_deeper(i) = tv%T(i,j,nz) ; S\_deeper(i) = tv%S(i,j,nz)}
704       \textcolor{comment}{!    N2\_region\_set(i) = .true.}
705       \textcolor{comment}{! endif}
706       \textcolor{keyword}{call }calculate\_density(t\_subml, s\_subml, pref\_n2, rho\_subml, tv%eqn\_of\_state, eosdom)
707       \textcolor{keyword}{call }calculate\_density(t\_deeper, s\_deeper, pref\_n2, rho\_deeper, tv%eqn\_of\_state, eosdom)
708       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} ((g%mask2dT(i,j)>0.5) .and. n2\_region\_set(i)) \textcolor{keywordflow}{then}
709         submln2(i,j) =  ge\_rho0 * (rho\_deeper(i) - rho\_subml(i)) / (gv%H\_to\_z * dh\_n2(i))
710 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
711 \textcolor{keywordflow}{    endif}
712 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! j-loop}
713 
714   \textcolor{keywordflow}{if} (id\_mld > 0) \textcolor{keyword}{call }post\_data(id\_mld, mld, diagptr)
715   \textcolor{keywordflow}{if} (id\_n2 > 0)  \textcolor{keyword}{call }post\_data(id\_n2, submln2, diagptr)
716   \textcolor{keywordflow}{if} (id\_sq > 0)  \textcolor{keyword}{call }post\_data(id\_sq, mld2, diagptr)
717 
\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 723 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
723   \textcolor{comment}{! Author: Brandon Reichl}
724   \textcolor{comment}{! Date: October 2, 2020}
725   \textcolor{comment}{! //}
726   \textcolor{comment}{! *Note that gravity is assumed constant everywhere and divided out of all calculations.}
727   \textcolor{comment}{!}
728   \textcolor{comment}{! This code has been written to step through the columns layer by layer, summing the PE}
729   \textcolor{comment}{! change inferred by mixing the layer with all layers above.  When the change exceeds a}
730   \textcolor{comment}{! threshold (determined by input array Mixing\_Energy), the code needs to solve for how far}
731   \textcolor{comment}{! into this layer the threshold PE change occurs (assuming constant density layers).}
732   \textcolor{comment}{! This is expressed here via solving the function F(X) = 0 where:}
733   \textcolor{comment}{! F(X) = 0.5 * ( Ca*X^3/(D1+X) + Cb*X^2/(D1+X) + Cc*X/(D1+X) + Dc/(D1+X)}
734   \textcolor{comment}{!                + Ca2*X^2 + Cb2*X + Cc2)}
735   \textcolor{comment}{! where all coefficients are determined by the previous mixed layer depth, the}
736   \textcolor{comment}{! density of the previous mixed layer, the present layer thickness, and the present}
737   \textcolor{comment}{! layer density.  This equation is worked out by computing the total PE assuming constant}
738   \textcolor{comment}{! density in the mixed layer as well as in the remaining part of the present layer that is}
739   \textcolor{comment}{! not mixed.}
740   \textcolor{comment}{! To solve for X in this equation a Newton's method iteration is employed, which}
741   \textcolor{comment}{! converges extremely quickly (usually 1 guess) since this equation turns out being rather}
742   \textcolor{comment}{! lienar for PE change with increasing X.}
743   \textcolor{comment}{! Input parameters:}
744   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(3)},   \textcolor{keywordtype}{intent(in)} :: id\_mld\textcolor{comment}{      !< Energy output diag IDs}
745   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{           !< Grid type}
746   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{          !< ocean vertical grid structure}
747   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{          !< A dimensional unit scaling type}
748   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(3)},      \textcolor{keywordtype}{intent(in)} :: mixing\_energy\textcolor{comment}{ !< Energy values for up to 3 MLDs}
749   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
750                            \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{           !< Layer thickness [H ~> m or kg m-2]}
751   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),   \textcolor{keywordtype}{intent(in)} :: tv\textcolor{comment}{          !< Structure containing pointers to any}
752 \textcolor{comment}{                                                     !! available thermodynamic fields.}
753   \textcolor{keywordtype}{type}(diag\_ctrl),         \textcolor{keywordtype}{pointer}    :: diagptr\textcolor{comment}{     !< Diagnostics structure}
754 
755   \textcolor{comment}{! Local variables}
756   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G),3)} :: mld     \textcolor{comment}{! Diagnosed mixed layer depth [Z ~> m].}
757   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: z\_l, z\_u, dz, rho\_c, pref\_mld
758   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(3)} :: pe\_threshold
759 
760   \textcolor{keywordtype}{real} :: ig, e\_g
761   \textcolor{keywordtype}{real} :: pe\_threshold\_fraction, pe, pe\_mixed, pe\_mixed\_tst
762   \textcolor{keywordtype}{real} :: rhodz\_ml, h\_ml, rhodz\_ml\_tst, h\_ml\_tst
763   \textcolor{keywordtype}{real} :: rho\_ml
764 
765   \textcolor{keywordtype}{real} :: r1, d1, r2, d2
766   \textcolor{keywordtype}{real} :: ca, cb,d ,cc, cd, ca2, cb2, c, cc2
767   \textcolor{keywordtype}{real} :: gx, gpx, hx, ihx, hpx, ix, ipx, fgx, fpx, x, x2
768 
769   \textcolor{keywordtype}{integer} :: it, im
770   \textcolor{keywordtype}{integer} :: i, j, is, ie, js, je, k, nz
771 
772   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
773 
774   pref\_mld(:) = 0.0
775   mld(:,:,:) = 0.0
776   pe\_threshold\_fraction = 1.e-4 \textcolor{comment}{!Fixed threshold of 0.01%, could be runtime.}
777 
778   \textcolor{keywordflow}{do} im=1,3
779     pe\_threshold(im) = mixing\_energy(im)/gv%g\_earth
780 \textcolor{keywordflow}{  enddo}
781 
782   \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is,ie
783     \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.0) \textcolor{keywordflow}{then}
784 
785       \textcolor{keyword}{call }calculate\_density(tv%T(i,j,:), tv%S(i,j,:), pref\_mld, rho\_c, 1, nz, &
786                              tv%eqn\_of\_state, scale=us%kg\_m3\_to\_R)
787 
788       \textcolor{keywordflow}{do} k=1,nz
789         dz(k) = h(i,j,k) * gv%H\_to\_Z
790 \textcolor{keywordflow}{      enddo}
791       z\_u(1) = 0.0
792       z\_l(1) = -dz(1)
793       \textcolor{keywordflow}{do} k=2,nz
794         z\_u(k) = z\_l(k-1)
795         z\_l(k) = z\_l(k-1)-dz(k)
796 \textcolor{keywordflow}{      enddo}
797 
798       \textcolor{keywordflow}{do} im=1,3
799 
800         \textcolor{comment}{! Initialize these for each column-wise calculation}
801         pe = 0.0
802         rhodz\_ml = 0.0
803         h\_ml = 0.0
804         rhodz\_ml\_tst = 0.0
805         h\_ml\_tst = 0.0
806         pe\_mixed = 0.0
807 
808         \textcolor{keywordflow}{do} k=1,nz
809 
810           \textcolor{comment}{! This is the unmixed PE cumulative sum from top down}
811           pe = pe + 0.5 * rho\_c(k) * (z\_u(k)**2 - z\_l(k)**2)
812 
813           \textcolor{comment}{! This is the depth and integral of density}
814           h\_ml\_tst = h\_ml + dz(k)
815           rhodz\_ml\_tst = rhodz\_ml + rho\_c(k) * dz(k)
816 
817           \textcolor{comment}{! The average density assuming all layers including this were mixed}
818           rho\_ml = rhodz\_ml\_tst/h\_ml\_tst
819 
820           \textcolor{comment}{! The PE assuming all layers including this were mixed}
821           \textcolor{comment}{! Note that 0. could be replaced with "Surface", which doesn't have to be 0}
822           \textcolor{comment}{! but 0 is a good reference value.}
823           pe\_mixed\_tst = 0.5 * rho\_ml * (0.**2 - (0. - h\_ml\_tst)**2)
824 
825           \textcolor{comment}{! Check if we supplied enough energy to mix to this layer}
826           \textcolor{keywordflow}{if} (pe\_mixed\_tst - pe <= pe\_threshold(im)) \textcolor{keywordflow}{then}
827             h\_ml = h\_ml\_tst
828             rhodz\_ml = rhodz\_ml\_tst
829 
830           \textcolor{keywordflow}{else} \textcolor{comment}{! If not, we need to solve where the energy ran out}
831             \textcolor{comment}{! This will be done with a Newton's method iteration:}
832 
833             r1 = rhodz\_ml / h\_ml \textcolor{comment}{! The density of the mixed layer (not including this layer)}
834             d1 = h\_ml \textcolor{comment}{! The thickness of the mixed layer (not including this layer)}
835             r2 = rho\_c(k) \textcolor{comment}{! The density of this layer}
836             d2 = dz(k) \textcolor{comment}{! The thickness of this layer}
837 
838             \textcolor{comment}{! This block could be used to calculate the function coefficients if}
839             \textcolor{comment}{! we don't reference all values to a surface designated as z=0}
840             \textcolor{comment}{! S = Surface}
841             \textcolor{comment}{! Ca  = -(R2)}
842             \textcolor{comment}{! Cb  = -( (R1*D1) + R2*(2.*D1-2.*S) )}
843             \textcolor{comment}{! D   = D1**2. - 2.*D1*S}
844             \textcolor{comment}{! Cc  = -( R1*D1*(2.*D1-2.*S) + R2*D )}
845             \textcolor{comment}{! Cd  = -(R1*D1*D)}
846             \textcolor{comment}{! Ca2 = R2}
847             \textcolor{comment}{! Cb2 = R2*(2*D1-2*S)}
848             \textcolor{comment}{! C   = S**2 + D2**2 + D1**2 - 2*D1*S - 2.*D2*S +2.*D1*D2}
849             \textcolor{comment}{! Cc2 = R2*(D+S**2-C)}
850             \textcolor{comment}{!}
851             \textcolor{comment}{! If the surface is S = 0, it simplifies to:}
852             ca  = -r2
853             cb  = -(r1 * d1 + r2 * (2. * d1))
854             d   = d1**2
855             cc  = -(r1 * d1 * (2. * d1) + (r2 * d))
856             cd  = -r1 * (d1 * d)
857             ca2 = r2
858             cb2 = r2 * (2. * d1)
859             c   = d2**2 + d1**2 + 2. * (d1 * d2)
860             cc2 = r2 * (d - c)
861 
862             \textcolor{comment}{! First guess for an iteration using Newton's method}
863             x = dz(k) * 0.5
864 
865             it=0
866             \textcolor{keywordflow}{do} \textcolor{keywordflow}{while}(it<10)\textcolor{comment}{!We can iterate up to 10 times}
867               \textcolor{comment}{! We are trying to solve the function:}
868               \textcolor{comment}{! F(x) = G(x)/H(x)+I(x)}
869               \textcolor{comment}{! for where F(x) = PE+PE\_threshold, or equivalently for where}
870               \textcolor{comment}{! F(x) = G(x)/H(x)+I(x) - (PE+PE\_threshold) = 0}
871               \textcolor{comment}{! We also need the derivative of this function for the Newton's method iteration}
872               \textcolor{comment}{! F'(x) = (G'(x)H(x)-G(x)H'(x))/H(x)^2 + I'(x)}
873               \textcolor{comment}{! G and its derivative}
874               gx = 0.5 * (ca * (x*x*x) + cb * x**2 + cc * x + cd)
875               gpx = 0.5 * (3. * (ca * x**2) + 2. * (cb * x) + cc)
876               \textcolor{comment}{! H, its inverse, and its derivative}
877               hx = d1 + x
878               ihx = 1. / hx
879               hpx = 1.
880               \textcolor{comment}{! I and its derivative}
881               ix = 0.5 * (ca2 * x**2 + cb2 * x + cc2)
882               ipx = 0.5 * (2. * ca2 * x + cb2)
883 
884               \textcolor{comment}{! The Function and its derivative:}
885               pe\_mixed = gx * ihx + ix
886               fgx = pe\_mixed - (pe + pe\_threshold(im))
887               fpx = (gpx * hx - hpx * gx) * ihx**2 + ipx
888 
889               \textcolor{comment}{! Check if our solution is within the threshold bounds, if not update}
890               \textcolor{comment}{! using Newton's method.  This appears to converge almost always in}
891               \textcolor{comment}{! one step because the function is very close to linear in most applications.}
892               \textcolor{keywordflow}{if} (abs(fgx) > pe\_threshold(im) * pe\_threshold\_fraction) \textcolor{keywordflow}{then}
893                 x2 = x - fgx / fpx
894                 it = it + 1
895                 \textcolor{keywordflow}{if} (x2 < 0. .or. x2 > dz(k)) \textcolor{keywordflow}{then}
896                   \textcolor{comment}{! The iteration seems to be robust, but we need to do something *if*}
897                   \textcolor{comment}{! things go wrong... How should we treat failed iteration?}
898                   \textcolor{comment}{! Present solution: Stop trying to compute and just say we can't mix this layer.}
899                   x=0
900                   \textcolor{keywordflow}{exit}
901                 \textcolor{keywordflow}{else}
902                   x = x2
903 \textcolor{keywordflow}{                endif}
904               \textcolor{keywordflow}{else}
905                 exit\textcolor{comment}{! Quit the iteration}
906 \textcolor{keywordflow}{              endif}
907 \textcolor{keywordflow}{            enddo}
908             h\_ml = h\_ml + x
909             exit\textcolor{comment}{! Quit looping through the column}
910 \textcolor{keywordflow}{          endif}
911 \textcolor{keywordflow}{        enddo}
912         mld(i,j,im) = h\_ml
913 \textcolor{keywordflow}{      enddo}
914     \textcolor{keywordflow}{else}
915       mld(i,j,:) = 0.0
916 \textcolor{keywordflow}{    endif}
917 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
918 
919   \textcolor{keywordflow}{if} (id\_mld(1) > 0) \textcolor{keyword}{call }post\_data(id\_mld(1), mld(:,:,1), diagptr)
920   \textcolor{keywordflow}{if} (id\_mld(2) > 0) \textcolor{keyword}{call }post\_data(id\_mld(2), mld(:,:,2), diagptr)
921   \textcolor{keywordflow}{if} (id\_mld(3) > 0) \textcolor{keyword}{call }post\_data(id\_mld(3), mld(:,:,3), diagptr)
922 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__diabatic__aux_a0b81bac82b18bd9f1e3e5a3beda78966}\label{namespacemom__diabatic__aux_a0b81bac82b18bd9f1e3e5a3beda78966}} 
\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[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(inout)}]{T,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(inout)}]{S,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke+1), intent(inout)}]{Kd\+\_\+T,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke+1), intent(in)}]{Kd\+\_\+S,  }\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 t} & Potential temperature \mbox{[}degC\mbox{]}.\\
\hline
\mbox{\tt in,out}  & {\em s} & Salinity \mbox{[}P\+SU\mbox{]} or \mbox{[}g\+Salt/kg\mbox{]}, generically \mbox{[}ppt\mbox{]}.\\
\hline
\mbox{\tt in,out}  & {\em kd\+\_\+t} & The extra diffusivity of temperature due to\\
\hline
\mbox{\tt in}  & {\em kd\+\_\+s} & The extra diffusivity of salinity due to\\
\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}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
232                            \textcolor{keywordtype}{intent(inout)} :: t\textcolor{comment}{    !< Potential temperature [degC].}
233   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
234                            \textcolor{keywordtype}{intent(inout)} :: s\textcolor{comment}{    !< Salinity [PSU] or [gSalt/kg], generically [ppt].}
235   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G)+1)}, &
236                            \textcolor{keywordtype}{intent(inout)}    :: kd\_t\textcolor{comment}{ !< The extra diffusivity of temperature due to}
237 \textcolor{comment}{                                                 !! double diffusion relative to the diffusivity of}
238 \textcolor{comment}{                                                 !! diffusivity of density [Z2 T-1 ~> m2 s-1].}
239   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G)+1)}, &
240                            \textcolor{keywordtype}{intent(in)}    :: kd\_s\textcolor{comment}{ !< The extra diffusivity of salinity due to}
241 \textcolor{comment}{                                                 !! double diffusion relative to the diffusivity of}
242 \textcolor{comment}{                                                 !! diffusivity of density [Z2 T-1 ~> m2 s-1].}
243   \textcolor{keywordtype}{real},                    \textcolor{keywordtype}{intent(in)}    :: dt\textcolor{comment}{   !<  Time increment [T ~> s].}
244 
245   \textcolor{comment}{! local variables}
246   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: &
247     b1\_t, b1\_s, &  \textcolor{comment}{!  Variables used by the tridiagonal solvers of T & S [H ~> m or kg m-2].}
248     d1\_t, d1\_s     \textcolor{comment}{!  Variables used by the tridiagonal solvers [nondim].}
249   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZK\_(G))} :: &
250     c1\_t, c1\_s     \textcolor{comment}{!  Variables used by the tridiagonal solvers [H ~> m or kg m-2].}
251   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZK\_(G)+1)} :: &
252     mix\_t, mix\_s   \textcolor{comment}{!  Mixing distances in both directions across each interface [H ~> m or kg m-2].}
253   \textcolor{keywordtype}{real} :: h\_tr         \textcolor{comment}{! h\_tr is h at tracer points with a tiny thickness}
254                        \textcolor{comment}{! added to ensure positive definiteness [H ~> m or kg m-2].}
255   \textcolor{keywordtype}{real} :: h\_neglect    \textcolor{comment}{! A thickness that is so small it is usually lost}
256                        \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
257   \textcolor{keywordtype}{real} :: i\_h\_int      \textcolor{comment}{! The inverse of the thickness associated with an}
258                        \textcolor{comment}{! interface [H-1 ~> m-1 or m2 kg-1].}
259   \textcolor{keywordtype}{real} :: b\_denom\_t    \textcolor{comment}{! The first term in the denominators for the expressions}
260   \textcolor{keywordtype}{real} :: b\_denom\_s    \textcolor{comment}{! for b1\_T and b1\_S, both [H ~> m or kg m-2].}
261   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, nz
262 
263   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
264   h\_neglect = gv%H\_subroundoff
265 
266   \textcolor{comment}{!$OMP parallel do default(private) shared(is,ie,js,je,h,h\_neglect,dt,Kd\_T,Kd\_S,G,GV,T,S,nz)}
267   \textcolor{keywordflow}{do} j=js,je
268     \textcolor{keywordflow}{do} i=is,ie
269       i\_h\_int = 1.0 / (0.5 * (h(i,j,1) + h(i,j,2)) + h\_neglect)
270       mix\_t(i,2) = ((dt * kd\_t(i,j,2)) * gv%Z\_to\_H**2) * i\_h\_int
271       mix\_s(i,2) = ((dt * kd\_s(i,j,2)) * gv%Z\_to\_H**2) * i\_h\_int
272 
273       h\_tr = h(i,j,1) + h\_neglect
274       b1\_t(i) = 1.0 / (h\_tr + mix\_t(i,2))
275       b1\_s(i) = 1.0 / (h\_tr + mix\_s(i,2))
276       d1\_t(i) = h\_tr * b1\_t(i)
277       d1\_s(i) = h\_tr * b1\_s(i)
278       t(i,j,1) = (b1\_t(i)*h\_tr)*t(i,j,1)
279       s(i,j,1) = (b1\_s(i)*h\_tr)*s(i,j,1)
280 \textcolor{keywordflow}{    enddo}
281     \textcolor{keywordflow}{do} k=2,nz-1 ; \textcolor{keywordflow}{do} i=is,ie
282       \textcolor{comment}{! Calculate the mixing across the interface below this layer.}
283       i\_h\_int = 1.0 / (0.5 * (h(i,j,k) + h(i,j,k+1)) + h\_neglect)
284       mix\_t(i,k+1) = ((dt * kd\_t(i,j,k+1)) * gv%Z\_to\_H**2) * i\_h\_int
285       mix\_s(i,k+1) = ((dt * kd\_s(i,j,k+1)) * gv%Z\_to\_H**2) * i\_h\_int
286 
287       c1\_t(i,k) = mix\_t(i,k) * b1\_t(i)
288       c1\_s(i,k) = mix\_s(i,k) * b1\_s(i)
289 
290       h\_tr = h(i,j,k) + h\_neglect
291       b\_denom\_t = h\_tr + d1\_t(i)*mix\_t(i,k)
292       b\_denom\_s = h\_tr + d1\_s(i)*mix\_s(i,k)
293       b1\_t(i) = 1.0 / (b\_denom\_t + mix\_t(i,k+1))
294       b1\_s(i) = 1.0 / (b\_denom\_s + mix\_s(i,k+1))
295       d1\_t(i) = b\_denom\_t * b1\_t(i)
296       d1\_s(i) = b\_denom\_s * b1\_s(i)
297 
298       t(i,j,k) = b1\_t(i) * (h\_tr*t(i,j,k) + mix\_t(i,k)*t(i,j,k-1))
299       s(i,j,k) = b1\_s(i) * (h\_tr*s(i,j,k) + mix\_s(i,k)*s(i,j,k-1))
300 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
301     \textcolor{keywordflow}{do} i=is,ie
302       c1\_t(i,nz) = mix\_t(i,nz) * b1\_t(i)
303       c1\_s(i,nz) = mix\_s(i,nz) * b1\_s(i)
304 
305       h\_tr = h(i,j,nz) + h\_neglect
306       b1\_t(i) = 1.0 / (h\_tr + d1\_t(i)*mix\_t(i,nz))
307       b1\_s(i) = 1.0 / (h\_tr + d1\_s(i)*mix\_s(i,nz))
308 
309       t(i,j,nz) = b1\_t(i) * (h\_tr*t(i,j,nz) + mix\_t(i,nz)*t(i,j,nz-1))
310       s(i,j,nz) = b1\_s(i) * (h\_tr*s(i,j,nz) + mix\_s(i,nz)*s(i,j,nz-1))
311 \textcolor{keywordflow}{    enddo}
312     \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=is,ie
313       t(i,j,k) = t(i,j,k) + c1\_t(i,k+1)*t(i,j,k+1)
314       s(i,j,k) = s(i,j,k) + c1\_s(i,k+1)*s(i,j,k+1)
315 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
316 \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 431 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


\begin{DoxyCode}
431   \textcolor{keywordtype}{type}(ocean\_grid\_type),     \textcolor{keywordtype}{intent(in)}  :: g\textcolor{comment}{    !< The ocean's grid structure}
432   \textcolor{keywordtype}{type}(verticalgrid\_type),   \textcolor{keywordtype}{intent(in)}  :: gv\textcolor{comment}{   !< The ocean's vertical grid structure}
433   \textcolor{keywordtype}{type}(unit\_scale\_type),     \textcolor{keywordtype}{intent(in)}  :: us\textcolor{comment}{   !< A dimensional unit scaling type}
434   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G))}, &
435                              \textcolor{keywordtype}{intent(in)}  :: u\textcolor{comment}{    !< The zonal velocity [L T-1 ~> m s-1]}
436   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G))}, &
437                              \textcolor{keywordtype}{intent(in)}  :: v\textcolor{comment}{    !< The meridional velocity [L T-1 ~> m s-1]}
438   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
439                              \textcolor{keywordtype}{intent(in)}  :: h\textcolor{comment}{    !< Layer thicknesses [H ~> m or kg m-2]}
440   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
441                              \textcolor{keywordtype}{intent(out)}   :: u\_h\textcolor{comment}{ !< Zonal velocity interpolated to h points [L T-1 ~> m
       s-1].}
442   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
443                              \textcolor{keywordtype}{intent(out)}   :: v\_h\textcolor{comment}{ !< Meridional velocity interpolated to h points [L T-1 ~>
       m s-1].}
444   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
445                      \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: ea\textcolor{comment}{ !< The amount of fluid entrained from the layer}
446 \textcolor{comment}{                                                 !! above within this time step [H ~> m or kg m-2].}
447 \textcolor{comment}{                                                 !! Omitting ea is the same as setting it to 0.}
448   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
449                      \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: eb\textcolor{comment}{ !< The amount of fluid entrained from the layer}
450 \textcolor{comment}{                                                 !! below within this time step [H ~> m or kg m-2].}
451 \textcolor{comment}{                                                 !! Omitting eb is the same as setting it to 0.}
452 
453   \textcolor{comment}{! local variables}
454   \textcolor{keywordtype}{real} :: b\_denom\_1    \textcolor{comment}{! The first term in the denominator of b1 [H ~> m or kg m-2].}
455   \textcolor{keywordtype}{real} :: h\_neglect    \textcolor{comment}{! A thickness that is so small it is usually lost}
456                        \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
457   \textcolor{keywordtype}{real} :: b1(szi\_(g)), d1(szi\_(g)), c1(szi\_(g),szk\_(g))
458   \textcolor{keywordtype}{real} :: a\_n(szi\_(g)), a\_s(szi\_(g))  \textcolor{comment}{! Fractional weights of the neighboring}
459   \textcolor{keywordtype}{real} :: a\_e(szi\_(g)), a\_w(szi\_(g))  \textcolor{comment}{! velocity points, ~1/2 in the open}
460                                       \textcolor{comment}{! ocean, nondimensional.}
461   \textcolor{keywordtype}{real} :: sum\_area, idenom
462   \textcolor{keywordtype}{logical} :: mix\_vertically
463   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, nz
464   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
465   \textcolor{keyword}{call }cpu\_clock\_begin(id\_clock\_uv\_at\_h)
466   h\_neglect = gv%H\_subroundoff
467 
468   mix\_vertically = \textcolor{keyword}{present}(ea)
469   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(ea) .neqv. \textcolor{keyword}{present}(eb)) \textcolor{keyword}{call }mom\_error(fatal, &
470       \textcolor{stringliteral}{"find\_uv\_at\_h: Either both ea and eb or neither one must be present "}// &
471       \textcolor{stringliteral}{"in call to find\_uv\_at\_h."})
472 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix\_vertically,h,h\_neglect, &}
473 \textcolor{comment}{!$OMP                                  eb,u\_h,u,v\_h,v,nz,ea)                     &}
474 \textcolor{comment}{!$OMP                          private(sum\_area,Idenom,a\_w,a\_e,a\_s,a\_n,b\_denom\_1,b1,d1,c1)}
475   \textcolor{keywordflow}{do} j=js,je
476     \textcolor{keywordflow}{do} i=is,ie
477       sum\_area = g%areaCu(i-1,j) + g%areaCu(i,j)
478       \textcolor{keywordflow}{if} (sum\_area>0.0) \textcolor{keywordflow}{then}
479         idenom = sqrt(0.5*g%IareaT(i,j) / sum\_area)
480         a\_w(i) = g%areaCu(i-1,j) * idenom
481         a\_e(i) = g%areaCu(i,j) * idenom
482       \textcolor{keywordflow}{else}
483         a\_w(i) = 0.0 ; a\_e(i) = 0.0
484 \textcolor{keywordflow}{      endif}
485 
486       sum\_area = g%areaCv(i,j-1) + g%areaCv(i,j)
487       \textcolor{keywordflow}{if} (sum\_area>0.0) \textcolor{keywordflow}{then}
488         idenom = sqrt(0.5*g%IareaT(i,j) / sum\_area)
489         a\_s(i) = g%areaCv(i,j-1) * idenom
490         a\_n(i) = g%areaCv(i,j) * idenom
491       \textcolor{keywordflow}{else}
492         a\_s(i) = 0.0 ; a\_n(i) = 0.0
493 \textcolor{keywordflow}{      endif}
494 \textcolor{keywordflow}{    enddo}
495 
496     \textcolor{keywordflow}{if} (mix\_vertically) \textcolor{keywordflow}{then}
497       \textcolor{keywordflow}{do} i=is,ie
498         b\_denom\_1 = h(i,j,1) + h\_neglect
499         b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,1))
500         d1(i) = b\_denom\_1 * b1(i)
501         u\_h(i,j,1) = (h(i,j,1)*b1(i)) * (a\_e(i)*u(i,j,1) + a\_w(i)*u(i-1,j,1))
502         v\_h(i,j,1) = (h(i,j,1)*b1(i)) * (a\_n(i)*v(i,j,1) + a\_s(i)*v(i,j-1,1))
503 \textcolor{keywordflow}{      enddo}
504       \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie
505         c1(i,k) = eb(i,j,k-1) * b1(i)
506         b\_denom\_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h\_neglect
507         b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,k))
508         d1(i) = b\_denom\_1 * b1(i)
509         u\_h(i,j,k) = (h(i,j,k) * (a\_e(i)*u(i,j,k) + a\_w(i)*u(i-1,j,k)) + &
510                       ea(i,j,k)*u\_h(i,j,k-1))*b1(i)
511         v\_h(i,j,k) = (h(i,j,k) * (a\_n(i)*v(i,j,k) + a\_s(i)*v(i,j-1,k)) + &
512                       ea(i,j,k)*v\_h(i,j,k-1))*b1(i)
513 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
514       \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=is,ie
515         u\_h(i,j,k) = u\_h(i,j,k) + c1(i,k+1)*u\_h(i,j,k+1)
516         v\_h(i,j,k) = v\_h(i,j,k) + c1(i,k+1)*v\_h(i,j,k+1)
517 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
518     \textcolor{keywordflow}{else}
519       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
520         u\_h(i,j,k) = a\_e(i)*u(i,j,k) + a\_w(i)*u(i-1,j,k)
521         v\_h(i,j,k) = a\_n(i)*v(i,j,k) + a\_s(i)*v(i,j-1,k)
522 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
523 \textcolor{keywordflow}{    endif}
524 \textcolor{keywordflow}{  enddo}
525 
526   \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 531 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90.


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


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