\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 \mbox{\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 \mbox{\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 \mbox{\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 \mbox{\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 \mbox{\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 \mbox{\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 \mbox{\hyperlink{namespacemom__diabatic__aux_adce369573ff355b820320d9c3b048693}{set\+\_\+pen\+\_\+shortwave}} (optics, fluxes, G, GV, US, CS, opacity\+\_\+\+C\+Sp, tracer\+\_\+flow\+\_\+\+C\+Sp) \item subroutine, public \mbox{\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 \mbox{\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 \mbox{\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 \mbox{\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 \mbox{\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 \mbox{\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 \mbox{\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(\mbox{\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(\mbox{\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(\mbox{\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(\mbox{\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)) \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)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(\mbox{\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(\mbox{\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}