\hypertarget{namespacemom__diabatic__aux}{}\section{mom\+\_\+diabatic\+\_\+aux Module Reference} \label{namespacemom__diabatic__aux}\index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsection{Detailed Description} Provides functions for some diabatic processes such as fraxil, brine rejection, tendency due to surface flux divergence. This module contains the subroutines that, along with the subroutines that it calls, implements diapycnal mass and momentum fluxes and a bulk mixed layer. The diapycnal diffusion can be used without the bulk mixed layer. diabatic first determines the (diffusive) diapycnal mass fluxes based on the convergence of the buoyancy fluxes within each layer. The dual-\/stream entrainment scheme of Mac\+Dougall and Dewar (J\+PO, 1997) is used for combined diapycnal advection and diffusion, calculated implicitly and potentially with the Richardson number dependent mixing, as described by Hallberg (M\+WR, 2000). Diapycnal advection is fundamentally the residual of diapycnal diffusion, so the fully implicit upwind differencing scheme that is used is entirely appropriate. The downward buoyancy flux in each layer is determined from an implicit calculation based on the previously calculated flux of the layer above and an estimated flux in the layer below. This flux is subject to the following conditions\+: (1) the flux in the top and bottom layers are set by the boundary conditions, and (2) no layer may be driven below an Angstrom thick-\/ ness. If there is a bulk mixed layer, the buffer layer is treat-\/ ed as a fixed density layer with vanishingly small diffusivity. diabatic takes 5 arguments\+: the two velocities (u and v), the thicknesses (h), a structure containing the forcing fields, and the length of time over which to act (dt). The velocities and thickness are taken as inputs and modified within the subroutine. There is no limit on the time step. \subsection*{Data Types} \begin{DoxyCompactItemize} \item type \hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs} \begin{DoxyCompactList}\small\item\em Control structure for diabatic\+\_\+aux. \end{DoxyCompactList}\end{DoxyCompactItemize} \subsection*{Functions/\+Subroutines} \begin{DoxyCompactItemize} \item subroutine, public \hyperlink{namespacemom__diabatic__aux_adb99fd4e2092c17cd69143945733a6d9}{make\+\_\+frazil} (h, tv, G, GV, US, CS, p\+\_\+surf, halo) \begin{DoxyCompactList}\small\item\em Frazil formation keeps the temperature above the freezing point. This subroutine warms any water that is colder than the (currently surface) freezing point up to the freezing point and accumulates the required heat (in \mbox{[}Q R Z $\sim$$>$ J m-\/2\mbox{]}) in tvfrazil. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__diabatic__aux_aa9248acf0b1cf246a79092077bcb0b46}{differential\+\_\+diffuse\+\_\+t\+\_\+s} (h, tv, visc, dt, G, GV) \begin{DoxyCompactList}\small\item\em This subroutine applies double diffusion to T \& S, assuming no diapycal mass fluxes, using a simple triadiagonal solver. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__diabatic__aux_a737df97b7e0205b0680c3b901ff6e2ff}{adjust\+\_\+salt} (h, tv, G, GV, CS, halo) \begin{DoxyCompactList}\small\item\em This subroutine keeps salinity from falling below a small but positive threshold. This usually occurs when the ice model attempts to extract more salt then is actually available to it from the ocean. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__diabatic__aux_acde4c76c9a489b82cba0bac5ec1b1198}{tridiagts} (G, GV, is, ie, js, je, hold, ea, eb, T, S) \begin{DoxyCompactList}\small\item\em This is a simple tri-\/diagonal solver for T and S. \char`\"{}\+Simple\char`\"{} means it only uses arrays hold, ea and eb. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__diabatic__aux_a3b3416d0ee87cd1a0d54a910fb681e93}{find\+\_\+uv\+\_\+at\+\_\+h} (u, v, h, u\+\_\+h, v\+\_\+h, G, GV, US, ea, eb) \begin{DoxyCompactList}\small\item\em This subroutine calculates u\+\_\+h and v\+\_\+h (velocities at thickness points), optionally using the entrainment amounts passed in as arguments. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__diabatic__aux_adce369573ff355b820320d9c3b048693}{set\+\_\+pen\+\_\+shortwave} (optics, fluxes, G, GV, US, CS, opacity\+\_\+\+C\+Sp, tracer\+\_\+flow\+\_\+\+C\+Sp) \item subroutine, public \hyperlink{namespacemom__diabatic__aux_a43ac98433f0f2f0e6cf4f8a3c9622ec4}{diagnosemldbydensitydifference} (id\+\_\+\+M\+LD, h, tv, density\+Diff, G, GV, US, diag\+Ptr, id\+\_\+\+N2sub\+ML, id\+\_\+\+M\+L\+Dsq, dz\+\_\+sub\+ML) \begin{DoxyCompactList}\small\item\em Diagnose a mixed layer depth (M\+LD) determined by a given density difference with the surface. This routine is appropriate in M\+O\+M\+\_\+diabatic\+\_\+driver due to its position within the time stepping. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__diabatic__aux_a695474893b92e3e3418a2bed871a5d7e}{diagnosemldbyenergy} (id\+\_\+\+M\+LD, h, tv, G, GV, US, Mixing\+\_\+\+Energy, diag\+Ptr) \begin{DoxyCompactList}\small\item\em Diagnose a mixed layer depth (M\+LD) determined by the depth a given energy value would mix. This routine is appropriate in M\+O\+M\+\_\+diabatic\+\_\+driver due to its position within the time stepping. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__diabatic__aux_ae7a279c765fa370c302532c13a1adaca}{applyboundaryfluxesinout} (CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, aggregate\+\_\+\+F\+W\+\_\+forcing, evap\+\_\+\+C\+F\+L\+\_\+limit, minimum\+\_\+forcing\+\_\+depth, c\+T\+KE, d\+S\+V\+\_\+dT, d\+S\+V\+\_\+dS, Skin\+Buoy\+Flux) \begin{DoxyCompactList}\small\item\em Update the thickness, temperature, and salinity due to thermodynamic boundary forcing (contained in fluxes type) applied to h, tvT and tvS, and calculate the T\+KE implications of this heating. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__diabatic__aux_a967697feb7ea22e1430d2b673ebab544}{diabatic\+\_\+aux\+\_\+init} (Time, G, GV, US, param\+\_\+file, diag, CS, use\+A\+L\+Ealgorithm, use\+\_\+e\+P\+BL) \begin{DoxyCompactList}\small\item\em This subroutine initializes the parameters and control structure of the diabatic\+\_\+aux module. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__diabatic__aux_a071e066ca5ed6619c6a9515fb07f5567}{diabatic\+\_\+aux\+\_\+end} (CS) \begin{DoxyCompactList}\small\item\em This subroutine initializes the control structure and any related memory for the diabatic\+\_\+aux module. \end{DoxyCompactList}\end{DoxyCompactItemize} \subsection*{Variables} \textbf{ }\par \begin{DoxyCompactItemize} \item \mbox{\Hypertarget{namespacemom__diabatic__aux_a1fffddaa09c8650a74e045c2fe29e256}\label{namespacemom__diabatic__aux_a1fffddaa09c8650a74e045c2fe29e256}} integer \hyperlink{namespacemom__diabatic__aux_a1fffddaa09c8650a74e045c2fe29e256}{id\+\_\+clock\+\_\+uv\+\_\+at\+\_\+h} \begin{DoxyCompactList}\small\item\em C\+PU time clock I\+Ds. \end{DoxyCompactList}\item \mbox{\Hypertarget{namespacemom__diabatic__aux_a87671b375428f3a27c18c61d6a93d035}\label{namespacemom__diabatic__aux_a87671b375428f3a27c18c61d6a93d035}} integer \hyperlink{namespacemom__diabatic__aux_a87671b375428f3a27c18c61d6a93d035}{id\+\_\+clock\+\_\+frazil} \begin{DoxyCompactList}\small\item\em C\+PU time clock I\+Ds. \end{DoxyCompactList}\end{DoxyCompactItemize} \subsection{Function/\+Subroutine Documentation} \mbox{\Hypertarget{namespacemom__diabatic__aux_a737df97b7e0205b0680c3b901ff6e2ff}\label{namespacemom__diabatic__aux_a737df97b7e0205b0680c3b901ff6e2ff}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!adjust\+\_\+salt@{adjust\+\_\+salt}} \index{adjust\+\_\+salt@{adjust\+\_\+salt}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{adjust\+\_\+salt()}{adjust\_salt()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::adjust\+\_\+salt (\begin{DoxyParamCaption}\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{h, }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(inout)}]{tv, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), intent(in)}]{CS, }\item[{integer, intent(in), optional}]{halo }\end{DoxyParamCaption})} This subroutine keeps salinity from falling below a small but positive threshold. This usually occurs when the ice model attempts to extract more salt then is actually available to it from the ocean. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em g} & The ocean\textquotesingle{}s grid structure\\ \hline \mbox{\tt in} & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\ \hline \mbox{\tt in} & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in,out} & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\ \hline \mbox{\tt in} & {\em cs} & The control structure returned by a previous call to diabatic\+\_\+aux\+\_\+init.\\ \hline \mbox{\tt in} & {\em halo} & Halo width over which to work \\ \hline \end{DoxyParams} Definition at line 330 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 330 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< The ocean's grid structure} 331 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< The ocean's vertical grid structure} 332 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 333 \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< Layer thicknesses [H ~> m or kg m-2]} 334 \textcolor{keywordtype}{type}(thermo\_var\_ptrs), \textcolor{keywordtype}{intent(inout)} :: tv\textcolor{comment}{ !< Structure containing pointers to any} 335 \textcolor{comment}{ !! available thermodynamic fields.} 336 \textcolor{keywordtype}{type}(diabatic\_aux\_cs), \textcolor{keywordtype}{intent(in)} :: cs\textcolor{comment}{ !< The control structure returned by a previous} 337 \textcolor{comment}{ !! call to diabatic\_aux\_init.} 338 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: halo\textcolor{comment}{ !< Halo width over which to work} 339 340 \textcolor{comment}{! local variables} 341 \textcolor{keywordtype}{real} :: salt\_add\_col(szi\_(g),szj\_(g))\textcolor{comment}{ !< The accumulated salt requirement [ppt R Z ~> gSalt m-2]} 342 \textcolor{keywordtype}{real} :: s\_min\textcolor{comment}{ !< The minimum salinity [ppt].} 343 \textcolor{keywordtype}{real} :: mc\textcolor{comment}{ !< A layer's mass [R Z ~> kg m-2].} 344 \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, nz 345 346 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke 347 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(halo)) \textcolor{keywordflow}{then} 348 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo 349 \textcolor{keywordflow}{ endif} 350 351 \textcolor{comment}{! call cpu\_clock\_begin(id\_clock\_adjust\_salt)} 352 353 s\_min = tv%min\_salinity 354 355 salt\_add\_col(:,:) = 0.0 356 357 \textcolor{comment}{!$OMP parallel do default(shared) private(mc)} 358 \textcolor{keywordflow}{do} j=js,je 359 \textcolor{keywordflow}{do} k=nz,1,-1 ; \textcolor{keywordflow}{do} i=is,ie 360 \textcolor{keywordflow}{if} ( (g%mask2dT(i,j) > 0.0) .and. & 361 ((tv%S(i,j,k) < s\_min) .or. (salt\_add\_col(i,j) > 0.0)) ) \textcolor{keywordflow}{then} 362 mc = gv%H\_to\_RZ * h(i,j,k) 363 \textcolor{keywordflow}{if} (h(i,j,k) <= 10.0*gv%Angstrom\_H) \textcolor{keywordflow}{then} 364 \textcolor{comment}{! Very thin layers should not be adjusted by the salt flux} 365 \textcolor{keywordflow}{if} (tv%S(i,j,k) < s\_min) \textcolor{keywordflow}{then} 366 salt\_add\_col(i,j) = salt\_add\_col(i,j) + mc * (s\_min - tv%S(i,j,k)) 367 tv%S(i,j,k) = s\_min 368 \textcolor{keywordflow}{ endif} 369 \textcolor{keywordflow}{elseif} (salt\_add\_col(i,j) + mc * (s\_min - tv%S(i,j,k)) <= 0.0) \textcolor{keywordflow}{then} 370 tv%S(i,j,k) = tv%S(i,j,k) - salt\_add\_col(i,j) / mc 371 salt\_add\_col(i,j) = 0.0 372 \textcolor{keywordflow}{else} 373 salt\_add\_col(i,j) = salt\_add\_col(i,j) + mc * (s\_min - tv%S(i,j,k)) 374 tv%S(i,j,k) = s\_min 375 \textcolor{keywordflow}{ endif} 376 \textcolor{keywordflow}{ endif} 377 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 378 \textcolor{keywordflow}{do} i=is,ie 379 tv%salt\_deficit(i,j) = tv%salt\_deficit(i,j) + salt\_add\_col(i,j) 380 \textcolor{keywordflow}{ enddo} 381 \textcolor{keywordflow}{ enddo} 382 \textcolor{comment}{! call cpu\_clock\_end(id\_clock\_adjust\_salt)} 383 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_ae7a279c765fa370c302532c13a1adaca}\label{namespacemom__diabatic__aux_ae7a279c765fa370c302532c13a1adaca}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!applyboundaryfluxesinout@{applyboundaryfluxesinout}} \index{applyboundaryfluxesinout@{applyboundaryfluxesinout}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{applyboundaryfluxesinout()}{applyboundaryfluxesinout()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::applyboundaryfluxesinout (\begin{DoxyParamCaption}\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), pointer}]{CS, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, intent(in)}]{dt, }\item[{type(forcing), intent(inout)}]{fluxes, }\item[{type(optics\+\_\+type), pointer}]{optics, }\item[{integer, intent(in)}]{nsw, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{h, }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(inout)}]{tv, }\item[{logical, intent(in)}]{aggregate\+\_\+\+F\+W\+\_\+forcing, }\item[{real, intent(in)}]{evap\+\_\+\+C\+F\+L\+\_\+limit, }\item[{real, intent(in)}]{minimum\+\_\+forcing\+\_\+depth, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out), optional}]{c\+T\+KE, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out), optional}]{d\+S\+V\+\_\+dT, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out), optional}]{d\+S\+V\+\_\+dS, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g)), intent(out), optional}]{Skin\+Buoy\+Flux }\end{DoxyParamCaption})} Update the thickness, temperature, and salinity due to thermodynamic boundary forcing (contained in fluxes type) applied to h, tvT and tvS, and calculate the T\+KE implications of this heating. \begin{DoxyParams}[1]{Parameters} & {\em cs} & Control structure for diabatic\+\_\+aux\\ \hline \mbox{\tt in} & {\em g} & Grid structure\\ \hline \mbox{\tt in} & {\em gv} & ocean vertical grid structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in} & {\em dt} & Time-\/step over which forcing is applied \mbox{[}T $\sim$$>$ s\mbox{]}\\ \hline \mbox{\tt in,out} & {\em fluxes} & Surface fluxes container\\ \hline & {\em optics} & Optical properties container\\ \hline \mbox{\tt in} & {\em nsw} & The number of frequency bands of penetrating shortwave radiation\\ \hline \mbox{\tt in,out} & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in,out} & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\ \hline \mbox{\tt in} & {\em aggregate\+\_\+fw\+\_\+forcing} & If False, treat in/out fluxes separately.\\ \hline \mbox{\tt in} & {\em evap\+\_\+cfl\+\_\+limit} & The largest fraction of a layer that can be evaporated in one time-\/step \mbox{[}nondim\mbox{]}.\\ \hline \mbox{\tt in} & {\em minimum\+\_\+forcing\+\_\+depth} & The smallest depth over which heat and freshwater fluxes is applied \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}.\\ \hline \mbox{\tt out} & {\em ctke} & Turbulent kinetic energy requirement to mix\\ \hline \mbox{\tt out} & {\em dsv\+\_\+dt} & Partial derivative of specific volume with\\ \hline \mbox{\tt out} & {\em dsv\+\_\+ds} & Partial derivative of specific volume with\\ \hline \mbox{\tt out} & {\em skinbuoyflux} & Buoyancy flux at surface \mbox{[}Z2 T-\/3 $\sim$$>$ m2 s-\/3\mbox{]}. \\ \hline \end{DoxyParams} Definition at line 939 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 939 \textcolor{keywordtype}{type}(diabatic\_aux\_cs), \textcolor{keywordtype}{pointer} :: cs\textcolor{comment}{ !< Control structure for diabatic\_aux} 940 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< Grid structure} 941 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< ocean vertical grid structure} 942 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 943 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: dt\textcolor{comment}{ !< Time-step over which forcing is applied [T ~> s]} 944 \textcolor{keywordtype}{type}(forcing), \textcolor{keywordtype}{intent(inout)} :: fluxes\textcolor{comment}{ !< Surface fluxes container} 945 \textcolor{keywordtype}{type}(optics\_type), \textcolor{keywordtype}{pointer} :: optics\textcolor{comment}{ !< Optical properties container} 946 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: nsw\textcolor{comment}{ !< The number of frequency bands of penetrating} 947 \textcolor{comment}{ !! shortwave radiation} 948 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 949 \textcolor{keywordtype}{intent(inout)} :: h\textcolor{comment}{ !< Layer thickness [H ~> m or kg m-2]} 950 \textcolor{keywordtype}{type}(thermo\_var\_ptrs), \textcolor{keywordtype}{intent(inout)} :: tv\textcolor{comment}{ !< Structure containing pointers to any} 951 \textcolor{comment}{ !! available thermodynamic fields.} 952 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{intent(in)} :: aggregate\_fw\_forcing\textcolor{comment}{ !< If False, treat in/out fluxes separately.} 953 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: evap\_cfl\_limit\textcolor{comment}{ !< The largest fraction of a layer that} 954 \textcolor{comment}{ !! can be evaporated in one time-step [nondim].} 955 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: minimum\_forcing\_depth\textcolor{comment}{ !< The smallest depth over which} 956 \textcolor{comment}{ !! heat and freshwater fluxes is applied [H ~> m or kg m-2].} 957 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 958 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(out)} :: ctke\textcolor{comment}{ !< Turbulent kinetic energy requirement to mix} 959 \textcolor{comment}{ !! forcing through each layer [R Z3 T-2 ~> J m-2]} 960 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 961 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(out)} :: dsv\_dt\textcolor{comment}{ !< Partial derivative of specific volume with} 962 \textcolor{comment}{ !! potential temperature [R-1 degC-1 ~> m3 kg-1 degC-1].} 963 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 964 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(out)} :: dsv\_ds\textcolor{comment}{ !< Partial derivative of specific volume with} 965 \textcolor{comment}{ !! salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].} 966 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, & 967 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(out)} :: skinbuoyflux\textcolor{comment}{ !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3].} 968 969 \textcolor{comment}{! Local variables} 970 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: maxgroundings = 5 971 \textcolor{keywordtype}{integer} :: numberofgroundings, iground(maxgroundings), jground(maxgroundings) 972 \textcolor{keywordtype}{real} :: h\_limit\_fluxes 973 \textcolor{keywordtype}{real} :: iforcingdepthscale 974 \textcolor{keywordtype}{real} :: idt \textcolor{comment}{! The inverse of the timestep [T-1 ~> s-1]} 975 \textcolor{keywordtype}{real} :: dthickness, dtemp, dsalt 976 \textcolor{keywordtype}{real} :: fractionofforcing, hold, ithickness 977 \textcolor{keywordtype}{real} :: rivermixconst \textcolor{comment}{! A constant used in implementing river mixing [R Z2 T-1 ~> Pa s].} 978 979 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: & 980 d\_pres, & \textcolor{comment}{! pressure change across a layer [R L2 T-2 ~> Pa]} 981 p\_lay, & \textcolor{comment}{! average pressure in a layer [R L2 T-2 ~> Pa]} 982 pres, & \textcolor{comment}{! pressure at an interface [R L2 T-2 ~> Pa]} 983 netmassinout, & \textcolor{comment}{! surface water fluxes [H ~> m or kg m-2] over time step} 984 netmassin, & \textcolor{comment}{! mass entering ocean surface [H ~> m or kg m-2] over a time step} 985 netmassout, & \textcolor{comment}{! mass leaving ocean surface [H ~> m or kg m-2] over a time step} 986 netheat, & \textcolor{comment}{! heat via surface fluxes excluding Pen\_SW\_bnd and netMassOut} 987 \textcolor{comment}{! [degC H ~> degC m or degC kg m-2]} 988 netsalt, & \textcolor{comment}{! surface salt flux ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )} 989 \textcolor{comment}{! [ppt H ~> ppt m or ppt kg m-2]} 990 nonpensw, & \textcolor{comment}{! non-downwelling SW, which is absorbed at ocean surface} 991 \textcolor{comment}{! [degC H ~> degC m or degC kg m-2]} 992 surfpressure, & \textcolor{comment}{! Surface pressure (approximated as 0.0) [R L2 T-2 ~> Pa]} 993 drhodt, & \textcolor{comment}{! change in density per change in temperature [R degC-1 ~> kg m-3 degC-1]} 994 drhods, & \textcolor{comment}{! change in density per change in salinity [R ppt-1 ~> kg m-3 ppt-1]} 995 netheat\_rate, & \textcolor{comment}{! netheat but for dt=1 [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]} 996 netsalt\_rate, & \textcolor{comment}{! netsalt but for dt=1 (e.g. returns a rate)} 997 \textcolor{comment}{! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]} 998 netmassinout\_rate\textcolor{comment}{! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1]} 999 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZK\_(G))} :: & 1000 h2d, & \textcolor{comment}{! A 2-d copy of the thicknesses [H ~> m or kg m-2]} 1001 t2d, & \textcolor{comment}{! A 2-d copy of the layer temperatures [degC]} 1002 pen\_tke\_2d, & \textcolor{comment}{! The TKE required to homogenize the heating by shortwave radiation within} 1003 \textcolor{comment}{! a layer [R Z3 T-2 ~> J m-2]} 1004 dsv\_dt\_2d \textcolor{comment}{! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]} 1005 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: & 1006 netpen\_rate \textcolor{comment}{! The surface penetrative shortwave heating rate summed over all bands} 1007 \textcolor{comment}{! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]} 1008 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(max(nsw,1),SZI\_(G))} :: & 1009 pen\_sw\_bnd, & \textcolor{comment}{! The penetrative shortwave heating integrated over a timestep by band} 1010 \textcolor{comment}{! [degC H ~> degC m or degC kg m-2]} 1011 pen\_sw\_bnd\_rate \textcolor{comment}{! The penetrative shortwave heating rate by band} 1012 \textcolor{comment}{! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]} 1013 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(max(nsw,1),SZI\_(G),SZK\_(G))} :: & 1014 opacityband \textcolor{comment}{! The opacity (inverse of the exponential absorption length) of each frequency} 1015 \textcolor{comment}{! band of shortwave radation in each layer [H-1 ~> m-1 or m2 kg-1]} 1016 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(maxGroundings)} :: hgrounding \textcolor{comment}{! Thickness added by each grounding event [H ~> m or kg m-2]} 1017 \textcolor{keywordtype}{real} :: temp\_in, salin\_in 1018 \textcolor{keywordtype}{real} :: g\_hconv2 \textcolor{comment}{! A conversion factor for use in the TKE calculation} 1019 \textcolor{comment}{! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].} 1020 \textcolor{keywordtype}{real} :: gorho \textcolor{comment}{! g\_Earth times a unit conversion factor divided by density} 1021 \textcolor{comment}{! [Z T-2 R-1 ~> m4 s-2 kg-1]} 1022 \textcolor{keywordtype}{logical} :: calculate\_energetics 1023 \textcolor{keywordtype}{logical} :: calculate\_buoyancy 1024 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(2)} :: eosdom \textcolor{comment}{! The i-computational domain for the equation of state} 1025 \textcolor{keywordtype}{integer} :: i, j, is, ie, js, je, k, nz, n, nb 1026 \textcolor{keywordtype}{character(len=45)} :: mesg 1027 1028 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke 1029 1030 idt = 1.0 / dt 1031 1032 calculate\_energetics = (\textcolor{keyword}{present}(ctke) .and. \textcolor{keyword}{present}(dsv\_dt) .and. \textcolor{keyword}{present}(dsv\_ds)) 1033 calculate\_buoyancy = \textcolor{keyword}{present}(skinbuoyflux) 1034 \textcolor{keywordflow}{if} (calculate\_buoyancy) skinbuoyflux(:,:) = 0.0 1035 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(ctke)) ctke(:,:,:) = 0.0 1036 g\_hconv2 = (us%L\_to\_Z**2*gv%g\_Earth * gv%H\_to\_RZ) * gv%H\_to\_RZ 1037 eosdom(:) = eos\_domain(g%HI) 1038 1039 \textcolor{comment}{! Only apply forcing if fluxes%sw is associated.} 1040 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(fluxes%sw) .and. .not.calculate\_energetics) \textcolor{keywordflow}{return} 1041 1042 \textcolor{keywordflow}{if} (calculate\_buoyancy) \textcolor{keywordflow}{then} 1043 surfpressure(:) = 0.0 1044 gorho = us%L\_to\_Z**2*gv%g\_Earth / gv%Rho0 1045 \textcolor{keywordflow}{ endif} 1046 1047 \textcolor{comment}{! H\_limit\_fluxes is used by extractFluxes1d to scale down fluxes if the total} 1048 \textcolor{comment}{! depth of the ocean is vanishing. It does not (yet) handle a value of zero.} 1049 \textcolor{comment}{! To accommodate vanishing upper layers, we need to allow for an instantaneous} 1050 \textcolor{comment}{! distribution of forcing over some finite vertical extent. The bulk mixed layer} 1051 \textcolor{comment}{! code handles this issue properly.} 1052 h\_limit\_fluxes = max(gv%Angstrom\_H, 1.e-30*gv%m\_to\_H) 1053 1054 \textcolor{comment}{! diagnostic to see if need to create mass to avoid grounding} 1055 \textcolor{keywordflow}{if} (cs%id\_createdH>0) cs%createdH(:,:) = 0. 1056 numberofgroundings = 0 1057 1058 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,tv,nsw,G,GV,US,optics,fluxes, &} 1059 \textcolor{comment}{!$OMP H\_limit\_fluxes,numberOfGroundings,iGround,jGround,&} 1060 \textcolor{comment}{!$OMP nonPenSW,hGrounding,CS,Idt,aggregate\_FW\_forcing, &} 1061 \textcolor{comment}{!$OMP minimum\_forcing\_depth,evap\_CFL\_limit,dt,EOSdom, &} 1062 \textcolor{comment}{!$OMP calculate\_buoyancy,netPen\_rate,SkinBuoyFlux,GoRho, &} 1063 \textcolor{comment}{!$OMP calculate\_energetics,dSV\_dT,dSV\_dS,cTKE,g\_Hconv2) &} 1064 \textcolor{comment}{!$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, &} 1065 \textcolor{comment}{!$OMP netHeat,netSalt,Pen\_SW\_bnd,fractionOfForcing, &} 1066 \textcolor{comment}{!$OMP IforcingDepthScale, &} 1067 \textcolor{comment}{!$OMP dThickness,dTemp,dSalt,hOld,Ithickness, &} 1068 \textcolor{comment}{!$OMP netMassIn,pres,d\_pres,p\_lay,dSV\_dT\_2d, &} 1069 \textcolor{comment}{!$OMP netmassinout\_rate,netheat\_rate,netsalt\_rate, &} 1070 \textcolor{comment}{!$OMP drhodt,drhods,pen\_sw\_bnd\_rate, &} 1071 \textcolor{comment}{!$OMP pen\_TKE\_2d,Temp\_in,Salin\_in,RivermixConst) &} 1072 \textcolor{comment}{!$OMP firstprivate(SurfPressure)} 1073 \textcolor{keywordflow}{do} j=js,je 1074 \textcolor{comment}{! Work in vertical slices for efficiency} 1075 1076 \textcolor{comment}{! Copy state into 2D-slice arrays} 1077 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie 1078 h2d(i,k) = h(i,j,k) 1079 t2d(i,k) = tv%T(i,j,k) 1080 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 1081 1082 \textcolor{keywordflow}{if} (calculate\_energetics) \textcolor{keywordflow}{then} 1083 \textcolor{comment}{! The partial derivatives of specific volume with temperature and} 1084 \textcolor{comment}{! salinity need to be precalculated to avoid having heating of} 1085 \textcolor{comment}{! tiny layers give nonsensical values.} 1086 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%p\_surf)) \textcolor{keywordflow}{then} 1087 \textcolor{keywordflow}{do} i=is,ie ; pres(i) = tv%p\_surf(i,j) ;\textcolor{keywordflow}{ enddo} 1088 \textcolor{keywordflow}{else} 1089 \textcolor{keywordflow}{do} i=is,ie ; pres(i) = 0.0 ;\textcolor{keywordflow}{ enddo} 1090 \textcolor{keywordflow}{ endif} 1091 \textcolor{keywordflow}{do} k=1,nz 1092 \textcolor{keywordflow}{do} i=is,ie 1093 d\_pres(i) = (gv%g\_Earth * gv%H\_to\_RZ) * h2d(i,k) 1094 p\_lay(i) = pres(i) + 0.5*d\_pres(i) 1095 pres(i) = pres(i) + d\_pres(i) 1096 \textcolor{keywordflow}{ enddo} 1097 \textcolor{keyword}{call }calculate\_specific\_vol\_derivs(t2d(:,k), tv%S(:,j,k), p\_lay(:), & 1098 dsv\_dt(:,j,k), dsv\_ds(:,j,k), tv%eqn\_of\_state, eosdom) 1099 \textcolor{keywordflow}{do} i=is,ie ; dsv\_dt\_2d(i,k) = dsv\_dt(i,j,k) ;\textcolor{keywordflow}{ enddo} 1100 \textcolor{keywordflow}{ enddo} 1101 pen\_tke\_2d(:,:) = 0.0 1102 \textcolor{keywordflow}{ endif} 1103 1104 \textcolor{comment}{! Nothing more is done on this j-slice if there is no buoyancy forcing.} 1105 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(fluxes%sw)) cycle 1106 1107 \textcolor{keywordflow}{if} (nsw>0) \textcolor{keyword}{call }extract\_optics\_slice(optics, j, g, gv, opacity=opacityband, opacity\_scale=(1.0/gv %m\_to\_H)) 1108 1109 \textcolor{comment}{! The surface forcing is contained in the fluxes type.} 1110 \textcolor{comment}{! We aggregate the thermodynamic forcing for a time step into the following:} 1111 \textcolor{comment}{! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step} 1112 \textcolor{comment}{! = lprec + fprec + vprec + evap + lrunoff + frunoff} 1113 \textcolor{comment}{! note that lprec generally has sea ice melt/form included.} 1114 \textcolor{comment}{! netMassOut = net mass leaving ocean surface [H ~> m or kg m-2] over a time step.} 1115 \textcolor{comment}{! netMassOut < 0 means mass leaves ocean.} 1116 \textcolor{comment}{! netHeat = heat via surface fluxes [degC H ~> degC m or degC kg m-2], excluding the part} 1117 \textcolor{comment}{! contained in Pen\_SW\_bnd; and excluding heat\_content of netMassOut < 0.} 1118 \textcolor{comment}{! netSalt = surface salt fluxes [ppt H ~> dppt m or gSalt m-2]} 1119 \textcolor{comment}{! Pen\_SW\_bnd = components to penetrative shortwave radiation split according to bands.} 1120 \textcolor{comment}{! This field provides that portion of SW from atmosphere that in fact} 1121 \textcolor{comment}{! enters to the ocean and participates in pentrative SW heating.} 1122 \textcolor{comment}{! nonpenSW = non-downwelling SW flux, which is absorbed in ocean surface} 1123 \textcolor{comment}{! (in tandem w/ LW,SENS,LAT); saved only for diagnostic purposes.} 1124 1125 \textcolor{comment}{!----------------------------------------------------------------------------------------} 1126 \textcolor{comment}{!BGR-June 26, 2017\{} 1127 \textcolor{comment}{!Temporary action to preserve answers while fixing a bug.} 1128 \textcolor{comment}{! To fix a bug in a diagnostic calculation, applyboundaryfluxesinout now returns} 1129 \textcolor{comment}{! the surface buoyancy flux. Previously, extractbuoyancyflux2d was called, meaning} 1130 \textcolor{comment}{! a second call to extractfluxes1d (causing the diagnostic net\_heat to be incorrect).} 1131 \textcolor{comment}{! Note that this call to extract buoyancyflux2d was AFTER applyboundaryfluxesinout,} 1132 \textcolor{comment}{! which means it used the T/S fields after this routine. Therefore, the surface} 1133 \textcolor{comment}{! buoyancy flux is computed here at the very end of this routine for legacy reasons.} 1134 \textcolor{comment}{! A few specific notes follow:} 1135 \textcolor{comment}{! 1) The old method did not included river/calving contributions to heat flux. This} 1136 \textcolor{comment}{! is kept consistent here via commenting code in the present extractFluxes1d <\_rate>} 1137 \textcolor{comment}{! outputs, but we may reconsider this approach.} 1138 \textcolor{comment}{! 2) The old method computed the buoyancy flux rate directly (by setting dt=1), instead} 1139 \textcolor{comment}{! of computing the integrated value (and dividing by dt). Hence the required} 1140 \textcolor{comment}{! additional outputs from extractFluxes1d.} 1141 \textcolor{comment}{! *** This is because: A*dt/dt =/= A due to round off.} 1142 \textcolor{comment}{! 3) The old method computed buoyancy flux after this routine, meaning the returned} 1143 \textcolor{comment}{! surface fluxes (from extractfluxes1d) must be recorded for use later in the code.} 1144 \textcolor{comment}{! We could (and maybe should) move that loop up to before the surface fluxes are} 1145 \textcolor{comment}{! applied, but this will change answers.} 1146 \textcolor{comment}{! For all these reasons we compute additional values of <\_rate> which are preserved} 1147 \textcolor{comment}{! for the buoyancy flux calculation and reproduce the old answers.} 1148 \textcolor{comment}{! In the future this needs more detailed investigation to make sure everything is} 1149 \textcolor{comment}{! consistent and correct. These details shouldnt significantly effect climate,} 1150 \textcolor{comment}{! but do change answers.} 1151 \textcolor{comment}{!-----------------------------------------------------------------------------------------} 1152 \textcolor{keywordflow}{if} (calculate\_buoyancy) \textcolor{keywordflow}{then} 1153 \textcolor{keyword}{call }extractfluxes1d(g, gv, us, fluxes, optics, nsw, j, dt, & 1154 h\_limit\_fluxes, cs%use\_river\_heat\_content, cs%use\_calving\_heat\_content, & 1155 h2d, t2d, netmassinout, netmassout, netheat, netsalt, & 1156 pen\_sw\_bnd, tv, aggregate\_fw\_forcing, nonpensw=nonpensw, & 1157 net\_heat\_rate=netheat\_rate, net\_salt\_rate=netsalt\_rate, & 1158 netmassinout\_rate=netmassinout\_rate, pen\_sw\_bnd\_rate=pen\_sw\_bnd\_rate) 1159 \textcolor{keywordflow}{else} 1160 \textcolor{keyword}{call }extractfluxes1d(g, gv, us, fluxes, optics, nsw, j, dt, & 1161 h\_limit\_fluxes, cs%use\_river\_heat\_content, cs%use\_calving\_heat\_content, & 1162 h2d, t2d, netmassinout, netmassout, netheat, netsalt, & 1163 pen\_sw\_bnd, tv, aggregate\_fw\_forcing, nonpensw=nonpensw) 1164 \textcolor{keywordflow}{ endif} 1165 \textcolor{comment}{! ea is for passive tracers} 1166 \textcolor{keywordflow}{do} i=is,ie 1167 \textcolor{comment}{! ea(i,j,1) = netMassInOut(i)} 1168 \textcolor{keywordflow}{if} (aggregate\_fw\_forcing) \textcolor{keywordflow}{then} 1169 netmassout(i) = netmassinout(i) 1170 netmassin(i) = 0. 1171 \textcolor{keywordflow}{else} 1172 netmassin(i) = netmassinout(i) - netmassout(i) 1173 \textcolor{keywordflow}{ endif} 1174 \textcolor{keywordflow}{if} (g%mask2dT(i,j)>0.0) \textcolor{keywordflow}{then} 1175 fluxes%netMassOut(i,j) = netmassout(i) 1176 fluxes%netMassIn(i,j) = netmassin(i) 1177 \textcolor{keywordflow}{else} 1178 fluxes%netMassOut(i,j) = 0.0 1179 fluxes%netMassIn(i,j) = 0.0 1180 \textcolor{keywordflow}{ endif} 1181 \textcolor{keywordflow}{ enddo} 1182 1183 \textcolor{comment}{! Apply the surface boundary fluxes in three steps:} 1184 \textcolor{comment}{! A/ update mass, temp, and salinity due to all terms except mass leaving} 1185 \textcolor{comment}{! ocean (and corresponding outward heat content), and ignoring penetrative SW.} 1186 \textcolor{comment}{! B/ update mass, salt, temp from mass leaving ocean.} 1187 \textcolor{comment}{! C/ update temp due to penetrative SW} 1188 \textcolor{keywordflow}{do} i=is,ie 1189 \textcolor{keywordflow}{if} (g%mask2dT(i,j)>0.) \textcolor{keywordflow}{then} 1190 1191 \textcolor{comment}{! A/ Update mass, temp, and salinity due to incoming mass flux.} 1192 \textcolor{keywordflow}{do} k=1,1 1193 1194 \textcolor{comment}{! Change in state due to forcing} 1195 dthickness = netmassin(i) \textcolor{comment}{! Since we are adding mass, we can use all of it} 1196 dtemp = 0. 1197 dsalt = 0. 1198 1199 \textcolor{comment}{! Update the forcing by the part to be consumed within the present k-layer.} 1200 \textcolor{comment}{! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.} 1201 netmassin(i) = netmassin(i) - dthickness 1202 \textcolor{comment}{! This line accounts for the temperature of the mass exchange} 1203 temp\_in = t2d(i,k) 1204 salin\_in = 0.0 1205 dtemp = dtemp + dthickness*temp\_in 1206 1207 \textcolor{comment}{! Diagnostics of heat content associated with mass fluxes} 1208 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(fluxes%heat\_content\_massin)) & 1209 fluxes%heat\_content\_massin(i,j) = fluxes%heat\_content\_massin(i,j) + & 1210 t2d(i,k) * max(0.,dthickness) * gv%H\_to\_RZ * fluxes%C\_p * idt 1211 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(fluxes%heat\_content\_massout)) & 1212 fluxes%heat\_content\_massout(i,j) = fluxes%heat\_content\_massout(i,j) + & 1213 t2d(i,k) * min(0.,dthickness) * gv%H\_to\_RZ * fluxes%C\_p * idt 1214 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + & 1215 t2d(i,k) * dthickness * gv%H\_to\_RZ 1216 1217 \textcolor{comment}{! Determine the energetics of river mixing before updating the state.} 1218 \textcolor{keywordflow}{if} (calculate\_energetics .and. \textcolor{keyword}{associated}(fluxes%lrunoff) .and. cs%do\_rivermix) \textcolor{keywordflow}{then} 1219 \textcolor{comment}{! Here we add an additional source of TKE to the mixed layer where river} 1220 \textcolor{comment}{! is present to simulate unresolved estuaries. The TKE input, TKE\_river in} 1221 \textcolor{comment}{! [Z3 T-3 ~> m3 s-3], is diagnosed as follows:} 1222 \textcolor{comment}{! TKE\_river = 0.5*rivermix\_depth*g*(1/rho)*drho\_ds*} 1223 \textcolor{comment}{! River*(Samb - Sriver) = CS%mstar*U\_star^3} 1224 \textcolor{comment}{! where River is in units of [Z T-1 ~> m s-1].} 1225 \textcolor{comment}{! Samb = Ambient salinity at the mouth of the estuary} 1226 \textcolor{comment}{! rivermix\_depth = The prescribed depth over which to mix river inflow} 1227 \textcolor{comment}{! drho\_ds = The gradient of density wrt salt at the ambient surface salinity.} 1228 \textcolor{comment}{! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)} 1229 \textcolor{keywordflow}{if} (gv%Boussinesq) \textcolor{keywordflow}{then} 1230 rivermixconst = -0.5*(cs%rivermix\_depth*dt) * ( us%L\_to\_Z**2*gv%g\_Earth ) * gv%Rho0 1231 \textcolor{keywordflow}{else} 1232 rivermixconst = -0.5*(cs%rivermix\_depth*dt) * gv%Rho0 * ( us%L\_to\_Z**2*gv%g\_Earth ) 1233 \textcolor{keywordflow}{ endif} 1234 ctke(i,j,k) = ctke(i,j,k) + max(0.0, rivermixconst*dsv\_ds(i,j,1) * & 1235 (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * tv%S(i,j,1)) 1236 \textcolor{keywordflow}{ endif} 1237 1238 \textcolor{comment}{! Update state} 1239 hold = h2d(i,k) \textcolor{comment}{! Keep original thickness in hand} 1240 h2d(i,k) = h2d(i,k) + dthickness \textcolor{comment}{! New thickness} 1241 \textcolor{keywordflow}{if} (h2d(i,k) > 0.0) \textcolor{keywordflow}{then} 1242 \textcolor{keywordflow}{if} (calculate\_energetics .and. (dthickness > 0.)) \textcolor{keywordflow}{then} 1243 \textcolor{comment}{! Calculate the energy required to mix the newly added water over} 1244 \textcolor{comment}{! the topmost grid cell.} 1245 ctke(i,j,k) = ctke(i,j,k) + 0.5*g\_hconv2*(hold*dthickness) * & 1246 ((t2d(i,k) - temp\_in) * dsv\_dt(i,j,k) + (tv%S(i,j,k) - salin\_in) * dsv\_ds(i,j,k)) 1247 \textcolor{keywordflow}{ endif} 1248 ithickness = 1.0/h2d(i,k) \textcolor{comment}{! Inverse new thickness} 1249 \textcolor{comment}{! The "if"s below avoid changing T/S by roundoff unnecessarily} 1250 \textcolor{keywordflow}{if} (dthickness /= 0. .or. dtemp /= 0.) t2d(i,k) = (hold*t2d(i,k) + dtemp)*ithickness 1251 \textcolor{keywordflow}{if} (dthickness /= 0. .or. dsalt /= 0.) tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness 1252 1253 \textcolor{keywordflow}{ endif} 1254 1255 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! k=1,1} 1256 1257 \textcolor{comment}{! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt.} 1258 \textcolor{keywordflow}{do} k=1,nz 1259 1260 \textcolor{comment}{! Place forcing into this layer if this layer has nontrivial thickness.} 1261 \textcolor{comment}{! For layers thin relative to 1/IforcingDepthScale, then distribute} 1262 \textcolor{comment}{! forcing into deeper layers.} 1263 iforcingdepthscale = 1. / max(gv%H\_subroundoff, minimum\_forcing\_depth - netmassout(i) ) 1264 \textcolor{comment}{! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.} 1265 fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale) 1266 1267 \textcolor{comment}{! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we} 1268 \textcolor{comment}{! limit the forcing applied to this cell, leaving the remaining forcing to} 1269 \textcolor{comment}{! be distributed downwards.} 1270 \textcolor{keywordflow}{if} (-fractionofforcing*netmassout(i) > evap\_cfl\_limit*h2d(i,k)) \textcolor{keywordflow}{then} 1271 fractionofforcing = -evap\_cfl\_limit*h2d(i,k)/netmassout(i) 1272 \textcolor{keywordflow}{ endif} 1273 1274 \textcolor{comment}{! Change in state due to forcing} 1275 1276 dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) ) 1277 dtemp = fractionofforcing*netheat(i) 1278 \textcolor{comment}{! ### The 0.9999 here should become a run-time parameter?} 1279 dsalt = max( fractionofforcing*netsalt(i), -0.9999*h2d(i,k)*tv%S(i,j,k)) 1280 1281 \textcolor{comment}{! Update the forcing by the part to be consumed within the present k-layer.} 1282 \textcolor{comment}{! If fractionOfForcing = 1, then new netMassOut vanishes.} 1283 netmassout(i) = netmassout(i) - dthickness 1284 netheat(i) = netheat(i) - dtemp 1285 netsalt(i) = netsalt(i) - dsalt 1286 1287 \textcolor{comment}{! This line accounts for the temperature of the mass exchange} 1288 dtemp = dtemp + dthickness*t2d(i,k) 1289 1290 \textcolor{comment}{! Diagnostics of heat content associated with mass fluxes} 1291 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(fluxes%heat\_content\_massin)) & 1292 fluxes%heat\_content\_massin(i,j) = fluxes%heat\_content\_massin(i,j) + & 1293 t2d(i,k) * max(0.,dthickness) * gv%H\_to\_RZ * fluxes%C\_p * idt 1294 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(fluxes%heat\_content\_massout)) & 1295 fluxes%heat\_content\_massout(i,j) = fluxes%heat\_content\_massout(i,j) + & 1296 t2d(i,k) * min(0.,dthickness) * gv%H\_to\_RZ * fluxes%C\_p * idt 1297 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + & 1298 t2d(i,k) * dthickness * gv%H\_to\_RZ 1299 1300 \textcolor{comment}{! Update state by the appropriate increment.} 1301 hold = h2d(i,k) \textcolor{comment}{! Keep original thickness in hand} 1302 h2d(i,k) = h2d(i,k) + dthickness \textcolor{comment}{! New thickness} 1303 1304 \textcolor{keywordflow}{if} (h2d(i,k) > 0.) \textcolor{keywordflow}{then} 1305 \textcolor{keywordflow}{if} (calculate\_energetics) \textcolor{keywordflow}{then} 1306 \textcolor{comment}{! Calculate the energy required to mix the newly added water over the topmost grid} 1307 \textcolor{comment}{! cell, assuming that the fluxes of heat and salt and rejected brine are initially} 1308 \textcolor{comment}{! applied in vanishingly thin layers at the top of the layer before being mixed} 1309 \textcolor{comment}{! throughout the layer. Note that dThickness is always <= 0 here, and that} 1310 \textcolor{comment}{! negative cTKE is a deficit that will need to be filled later.} 1311 ctke(i,j,k) = ctke(i,j,k) - (0.5*h2d(i,k)*g\_hconv2) * & 1312 ((dtemp - dthickness*t2d(i,k)) * dsv\_dt(i,j,k) + & 1313 (dsalt - dthickness*tv%S(i,j,k)) * dsv\_ds(i,j,k)) 1314 \textcolor{keywordflow}{ endif} 1315 ithickness = 1.0/h2d(i,k) \textcolor{comment}{! Inverse of new thickness} 1316 t2d(i,k) = (hold*t2d(i,k) + dtemp)*ithickness 1317 tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness 1318 \textcolor{keywordflow}{elseif} (h2d(i,k) < 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! h2d==0 is a special limit that needs no extra handling} 1319 \textcolor{keyword}{call }forcing\_singlepointprint(fluxes,g,i,j,\textcolor{stringliteral}{'applyBoundaryFluxesInOut (h<0)'}) 1320 \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{'applyBoundaryFluxesInOut(): lon,lat='},g%geoLonT(i,j),g%geoLatT(i,j) 1321 \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{'applyBoundaryFluxesInOut(): netT,netS,netH='},netheat(i),netsalt(i),netmassinout(i) 1322 \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{'applyBoundaryFluxesInOut(): dT,dS,dH='},dtemp,dsalt,dthickness 1323 \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{'applyBoundaryFluxesInOut(): h(n),h(n+1),k='},hold,h2d(i,k),k 1324 \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_diabatic\_driver.F90, applyBoundaryFluxesInOut(): "}//& 1325 \textcolor{stringliteral}{"Complete mass loss in column!"}) 1326 \textcolor{keywordflow}{ endif} 1327 1328 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! k} 1329 1330 \textcolor{comment}{! Check if trying to apply fluxes over land points} 1331 \textcolor{keywordflow}{elseif} ((abs(netheat(i))+abs(netsalt(i))+abs(netmassin(i))+abs(netmassout(i)))>0.) \textcolor{keywordflow}{then} 1332 1333 \textcolor{keywordflow}{if} (.not. cs%ignore\_fluxes\_over\_land) \textcolor{keywordflow}{then} 1334 \textcolor{keyword}{call }forcing\_singlepointprint(fluxes,g,i,j,\textcolor{stringliteral}{'applyBoundaryFluxesInOut (land)'}) 1335 \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{'applyBoundaryFluxesInOut(): lon,lat='},g%geoLonT(i,j),g%geoLatT(i,j) 1336 \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{'applyBoundaryFluxesInOut(): netHeat,netSalt,netMassIn,netMassOut='},& 1337 netheat(i),netsalt(i),netmassin(i),netmassout(i) 1338 1339 \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_diabatic\_driver.F90, applyBoundaryFluxesInOut(): "}//& 1340 \textcolor{stringliteral}{"Mass loss over land?"}) 1341 \textcolor{keywordflow}{ endif} 1342 1343 \textcolor{keywordflow}{ endif} 1344 1345 \textcolor{comment}{! If anything remains after the k-loop, then we have grounded out, which is a problem.} 1346 \textcolor{keywordflow}{if} (netmassin(i)+netmassout(i) /= 0.0) \textcolor{keywordflow}{then} 1347 \textcolor{comment}{!$OMP critical} 1348 numberofgroundings = numberofgroundings +1 1349 \textcolor{keywordflow}{if} (numberofgroundings<=maxgroundings) \textcolor{keywordflow}{then} 1350 iground(numberofgroundings) = i \textcolor{comment}{! Record i,j location of event for} 1351 jground(numberofgroundings) = j \textcolor{comment}{! warning message} 1352 hgrounding(numberofgroundings) = netmassin(i)+netmassout(i) 1353 \textcolor{keywordflow}{ endif} 1354 \textcolor{comment}{!$OMP end critical} 1355 \textcolor{keywordflow}{if} (cs%id\_createdH>0) cs%createdH(i,j) = cs%createdH(i,j) - (netmassin(i)+netmassout(i))/dt 1356 \textcolor{keywordflow}{ endif} 1357 1358 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! i} 1359 1360 \textcolor{comment}{! Step C/ in the application of fluxes} 1361 \textcolor{comment}{! Heat by the convergence of penetrating SW.} 1362 \textcolor{comment}{! SW penetrative heating uses the updated thickness from above.} 1363 1364 \textcolor{comment}{! Save temperature before increment with SW heating} 1365 \textcolor{comment}{! and initialize CS%penSWflux\_diag to zero.} 1366 \textcolor{keywordflow}{if} (cs%id\_penSW\_diag > 0 .or. cs%id\_penSWflux\_diag > 0) \textcolor{keywordflow}{then} 1367 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie 1368 cs%penSW\_diag(i,j,k) = t2d(i,k) 1369 cs%penSWflux\_diag(i,j,k) = 0.0 1370 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 1371 k=nz+1 ; \textcolor{keywordflow}{do} i=is,ie 1372 cs%penSWflux\_diag(i,j,k) = 0.0 1373 \textcolor{keywordflow}{ enddo} 1374 \textcolor{keywordflow}{ endif} 1375 1376 \textcolor{keywordflow}{if} (calculate\_energetics) \textcolor{keywordflow}{then} 1377 \textcolor{keyword}{call }absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h\_limit\_fluxes, & 1378 .false., .true., t2d, pen\_sw\_bnd, tke=pen\_tke\_2d, dsv\_dt=dsv\_dt\_2d) 1379 k = 1 \textcolor{comment}{! For setting break-points.} 1380 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie 1381 ctke(i,j,k) = ctke(i,j,k) + pen\_tke\_2d(i,k) 1382 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 1383 \textcolor{keywordflow}{else} 1384 \textcolor{keyword}{call }absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h\_limit\_fluxes, & 1385 .false., .true., t2d, pen\_sw\_bnd) 1386 \textcolor{keywordflow}{ endif} 1387 1388 1389 \textcolor{comment}{! Step D/ copy updated thickness and temperature} 1390 \textcolor{comment}{! 2d slice now back into model state.} 1391 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie 1392 h(i,j,k) = h2d(i,k) 1393 tv%T(i,j,k) = t2d(i,k) 1394 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 1395 1396 \textcolor{comment}{! Diagnose heating [Q R Z T-1 ~> W m-2] applied to a grid cell from SW penetration} 1397 \textcolor{comment}{! Also diagnose the penetrative SW heat flux at base of layer.} 1398 \textcolor{keywordflow}{if} (cs%id\_penSW\_diag > 0 .or. cs%id\_penSWflux\_diag > 0) \textcolor{keywordflow}{then} 1399 1400 \textcolor{comment}{! convergence of SW into a layer} 1401 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie 1402 cs%penSW\_diag(i,j,k) = (t2d(i,k)-cs%penSW\_diag(i,j,k))*h(i,j,k) * idt * tv%C\_p * gv%H\_to\_RZ 1403 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 1404 1405 \textcolor{comment}{! Perform a cumulative sum upwards from bottom to} 1406 \textcolor{comment}{! diagnose penetrative SW flux at base of tracer cell.} 1407 \textcolor{comment}{! CS%penSWflux\_diag(i,j,k=1) is penetrative shortwave at top of ocean.} 1408 \textcolor{comment}{! CS%penSWflux\_diag(i,j,k=kbot+1) is zero, since assume no SW penetrates rock.} 1409 \textcolor{comment}{! CS%penSWflux\_diag = rsdo and CS%penSW\_diag = rsdoabsorb} 1410 \textcolor{comment}{! rsdoabsorb(k) = rsdo(k) - rsdo(k+1), so that rsdo(k) = rsdo(k+1) + rsdoabsorb(k)} 1411 \textcolor{keywordflow}{if} (cs%id\_penSWflux\_diag > 0) \textcolor{keywordflow}{then} 1412 \textcolor{keywordflow}{do} k=nz,1,-1 ; \textcolor{keywordflow}{do} i=is,ie 1413 cs%penSWflux\_diag(i,j,k) = cs%penSW\_diag(i,j,k) + cs%penSWflux\_diag(i,j,k+1) 1414 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 1415 \textcolor{keywordflow}{ endif} 1416 1417 \textcolor{keywordflow}{ endif} 1418 1419 \textcolor{comment}{! Fill CS%nonpenSW\_diag} 1420 \textcolor{keywordflow}{if} (cs%id\_nonpenSW\_diag > 0) \textcolor{keywordflow}{then} 1421 \textcolor{keywordflow}{do} i=is,ie 1422 cs%nonpenSW\_diag(i,j) = nonpensw(i) * idt * tv%C\_p * gv%H\_to\_RZ 1423 \textcolor{keywordflow}{ enddo} 1424 \textcolor{keywordflow}{ endif} 1425 1426 \textcolor{comment}{! BGR: Get buoyancy flux to return for ePBL} 1427 \textcolor{comment}{! We want the rate, so we use the rate values returned from extractfluxes1d.} 1428 \textcolor{comment}{! Note that the *dt values could be divided by dt here, but} 1429 \textcolor{comment}{! 1) Answers will change due to round-off} 1430 \textcolor{comment}{! 2) Be sure to save their values BEFORE fluxes are used.} 1431 \textcolor{keywordflow}{if} (calculate\_buoyancy) \textcolor{keywordflow}{then} 1432 drhodt(:) = 0.0 1433 drhods(:) = 0.0 1434 netpen\_rate(:) = 0.0 1435 \textcolor{comment}{! Sum over bands and attenuate as a function of depth.} 1436 \textcolor{comment}{! netPen\_rate is the netSW as a function of depth, but only the surface value is used here,} 1437 \textcolor{comment}{! in which case the values of dt, h, optics and H\_limit\_fluxes are irrelevant. Consider} 1438 \textcolor{comment}{! writing a shorter and simpler variant to handle this very limited case.} 1439 \textcolor{comment}{! call sumSWoverBands(G, GV, US, h2d(:,:), optics\_nbands(optics), optics, j, dt, &} 1440 \textcolor{comment}{! H\_limit\_fluxes, .true., pen\_SW\_bnd\_rate, netPen)} 1441 \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{do} nb=1,nsw ; netpen\_rate(i) = netpen\_rate(i) + pen\_sw\_bnd\_rate(nb,i) ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 1442 1443 \textcolor{comment}{! Density derivatives} 1444 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%p\_surf)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=is,ie ; surfpressure(i) = tv%p\_surf(i,j) ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 1445 \textcolor{keyword}{call }calculate\_density\_derivs(t2d(:,1), tv%S(:,j,1), surfpressure, drhodt, drhods, & 1446 tv%eqn\_of\_state, eosdom) 1447 \textcolor{comment}{! 1. Adjust netSalt to reflect dilution effect of FW flux} 1448 \textcolor{comment}{! 2. Add in the SW heating for purposes of calculating the net} 1449 \textcolor{comment}{! surface buoyancy flux affecting the top layer.} 1450 \textcolor{comment}{! 3. Convert to a buoyancy flux, excluding penetrating SW heating} 1451 \textcolor{comment}{! BGR-Jul 5, 2017: The contribution of SW heating here needs investigated for ePBL.} 1452 \textcolor{keywordflow}{do} i=is,ie 1453 skinbuoyflux(i,j) = - gorho * gv%H\_to\_Z * & 1454 (drhods(i) * (netsalt\_rate(i) - tv%S(i,j,1)*netmassinout\_rate(i)) + & 1455 drhodt(i) * ( netheat\_rate(i) + netpen\_rate(i)) ) \textcolor{comment}{! [Z2 T-3 ~> m2 s-3]} 1456 \textcolor{keywordflow}{ enddo} 1457 \textcolor{keywordflow}{ endif} 1458 1459 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! j-loop finish} 1460 1461 \textcolor{comment}{! Post the diagnostics} 1462 \textcolor{keywordflow}{if} (cs%id\_createdH > 0) \textcolor{keyword}{call }post\_data(cs%id\_createdH , cs%createdH , cs%diag) 1463 \textcolor{keywordflow}{if} (cs%id\_penSW\_diag > 0) \textcolor{keyword}{call }post\_data(cs%id\_penSW\_diag , cs%penSW\_diag , cs%diag) 1464 \textcolor{keywordflow}{if} (cs%id\_penSWflux\_diag > 0) \textcolor{keyword}{call }post\_data(cs%id\_penSWflux\_diag, cs%penSWflux\_diag, cs%diag) 1465 \textcolor{keywordflow}{if} (cs%id\_nonpenSW\_diag > 0) \textcolor{keyword}{call }post\_data(cs%id\_nonpenSW\_diag , cs%nonpenSW\_diag , cs%diag) 1466 1467 \textcolor{comment}{! The following check will be ignored if ignore\_fluxes\_over\_land = true} 1468 \textcolor{keywordflow}{if} (numberofgroundings>0 .and. .not. cs%ignore\_fluxes\_over\_land) \textcolor{keywordflow}{then} 1469 \textcolor{keywordflow}{do} i = 1, min(numberofgroundings, maxgroundings) 1470 \textcolor{keyword}{call }forcing\_singlepointprint(fluxes,g,iground(i),jground(i),\textcolor{stringliteral}{'applyBoundaryFluxesInOut (grounding)'}) 1471 \textcolor{keyword}{write}(mesg(1:45),\textcolor{stringliteral}{'(3es15.3)'}) g%geoLonT( iground(i), jground(i) ), & 1472 g%geoLatT( iground(i), jground(i)), hgrounding(i)*gv%H\_to\_m 1473 \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"MOM\_diabatic\_driver.F90, applyBoundaryFluxesInOut(): "}//& 1474 \textcolor{stringliteral}{"Mass created. x,y,dh= "}//trim(mesg), all\_print=.true.) 1475 \textcolor{keywordflow}{ enddo} 1476 1477 \textcolor{keywordflow}{if} (numberofgroundings - maxgroundings > 0) \textcolor{keywordflow}{then} 1478 \textcolor{keyword}{write}(mesg, \textcolor{stringliteral}{'(i4)'}) numberofgroundings - maxgroundings 1479 \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"MOM\_diabatic\_driver:F90, applyBoundaryFluxesInOut(): "}//& 1480 trim(mesg) // \textcolor{stringliteral}{" groundings remaining"}) 1481 \textcolor{keywordflow}{ endif} 1482 \textcolor{keywordflow}{ endif} 1483 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_a071e066ca5ed6619c6a9515fb07f5567}\label{namespacemom__diabatic__aux_a071e066ca5ed6619c6a9515fb07f5567}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!diabatic\+\_\+aux\+\_\+end@{diabatic\+\_\+aux\+\_\+end}} \index{diabatic\+\_\+aux\+\_\+end@{diabatic\+\_\+aux\+\_\+end}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{diabatic\+\_\+aux\+\_\+end()}{diabatic\_aux\_end()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::diabatic\+\_\+aux\+\_\+end (\begin{DoxyParamCaption}\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), pointer}]{CS }\end{DoxyParamCaption})} This subroutine initializes the control structure and any related memory for the diabatic\+\_\+aux module. \begin{DoxyParams}{Parameters} {\em cs} & The control structure returned by a previous call to diabatic\+\_\+aux\+\_\+init; it is deallocated here. \\ \hline \end{DoxyParams} Definition at line 1652 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 1652 \textcolor{keywordtype}{type}(diabatic\_aux\_cs), \textcolor{keywordtype}{pointer} :: cs\textcolor{comment}{ !< The control structure returned by a previous} 1653 \textcolor{comment}{ !! call to diabatic\_aux\_init; it is deallocated here.} 1654 1655 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{return} 1656 1657 \textcolor{keywordflow}{if} (cs%id\_createdH >0) \textcolor{keyword}{deallocate}(cs%createdH) 1658 \textcolor{keywordflow}{if} (cs%id\_penSW\_diag >0) \textcolor{keyword}{deallocate}(cs%penSW\_diag) 1659 \textcolor{keywordflow}{if} (cs%id\_penSWflux\_diag >0) \textcolor{keyword}{deallocate}(cs%penSWflux\_diag) 1660 \textcolor{keywordflow}{if} (cs%id\_nonpenSW\_diag >0) \textcolor{keyword}{deallocate}(cs%nonpenSW\_diag) 1661 1662 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{deallocate}(cs) 1663 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_a967697feb7ea22e1430d2b673ebab544}\label{namespacemom__diabatic__aux_a967697feb7ea22e1430d2b673ebab544}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!diabatic\+\_\+aux\+\_\+init@{diabatic\+\_\+aux\+\_\+init}} \index{diabatic\+\_\+aux\+\_\+init@{diabatic\+\_\+aux\+\_\+init}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{diabatic\+\_\+aux\+\_\+init()}{diabatic\_aux\_init()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::diabatic\+\_\+aux\+\_\+init (\begin{DoxyParamCaption}\item[{type(time\+\_\+type), intent(in), target}]{Time, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{type(param\+\_\+file\+\_\+type), intent(in)}]{param\+\_\+file, }\item[{type(diag\+\_\+ctrl), intent(inout), target}]{diag, }\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), pointer}]{CS, }\item[{logical, intent(in)}]{use\+A\+L\+Ealgorithm, }\item[{logical, intent(in)}]{use\+\_\+e\+P\+BL }\end{DoxyParamCaption})} This subroutine initializes the parameters and control structure of the diabatic\+\_\+aux module. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em time} & The current model time.\\ \hline \mbox{\tt in} & {\em g} & The ocean\textquotesingle{}s grid structure\\ \hline \mbox{\tt in} & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in} & {\em param\+\_\+file} & A structure to parse for run-\/time parameters\\ \hline \mbox{\tt in,out} & {\em diag} & A structure used to regulate diagnostic output\\ \hline & {\em cs} & A pointer to the control structure for the diabatic\+\_\+aux module, which is initialized here.\\ \hline \mbox{\tt in} & {\em usealealgorithm} & If true, use the A\+LE algorithm rather than layered mode.\\ \hline \mbox{\tt in} & {\em use\+\_\+epbl} & If true, use the implicit energetics planetary boundary layer scheme to determine the diffusivity in the surface boundary layer. \\ \hline \end{DoxyParams} Definition at line 1488 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 1488 \textcolor{keywordtype}{type}(time\_type), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(in)} :: time\textcolor{comment}{ !< The current model time.} 1489 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< The ocean's grid structure} 1490 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< The ocean's vertical grid structure} 1491 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 1492 \textcolor{keywordtype}{type}(param\_file\_type), \textcolor{keywordtype}{intent(in)} :: param\_file\textcolor{comment}{ !< A structure to parse for run-time parameters} 1493 \textcolor{keywordtype}{type}(diag\_ctrl), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(inout)} :: diag\textcolor{comment}{ !< A structure used to regulate diagnostic output} 1494 \textcolor{keywordtype}{type}(diabatic\_aux\_cs), \textcolor{keywordtype}{pointer} :: cs\textcolor{comment}{ !< A pointer to the control structure for the} 1495 \textcolor{comment}{ !! diabatic\_aux module, which is initialized here.} 1496 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{intent(in)} :: usealealgorithm\textcolor{comment}{ !< If true, use the ALE algorithm rather} 1497 \textcolor{comment}{ !! than layered mode.} 1498 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{intent(in)} :: use\_epbl\textcolor{comment}{ !< If true, use the implicit energetics planetary} 1499 \textcolor{comment}{ !! boundary layer scheme to determine the diffusivity} 1500 \textcolor{comment}{ !! in the surface boundary layer.} 1501 1502 \textcolor{comment}{! This "include" declares and sets the variable "version".} 1503 \textcolor{preprocessor}{#include "version\_variable.h"} 1504 \textcolor{preprocessor}{} \textcolor{keywordtype}{character(len=40)} :: mdl = \textcolor{stringliteral}{"MOM\_diabatic\_aux"} \textcolor{comment}{! This module's name.} 1505 \textcolor{keywordtype}{character(len=48)} :: thickness\_units 1506 \textcolor{keywordtype}{character(len=200)} :: inputdir \textcolor{comment}{! The directory where NetCDF input files} 1507 \textcolor{keywordtype}{character(len=240)} :: chl\_filename \textcolor{comment}{! A file from which chl\_a concentrations are to be read.} 1508 \textcolor{keywordtype}{character(len=128)} :: chl\_file \textcolor{comment}{! Data containing chl\_a concentrations. Used} 1509 \textcolor{comment}{! when var\_pen\_sw is defined and reading from file.} 1510 \textcolor{keywordtype}{character(len=32)} :: chl\_varname \textcolor{comment}{! Name of chl\_a variable in chl\_file.} 1511 \textcolor{keywordtype}{logical} :: use\_temperature \textcolor{comment}{! True if thermodynamics are enabled.} 1512 \textcolor{keywordtype}{integer} :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nbands 1513 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke 1514 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB 1515 1516 \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{then} 1517 \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"diabatic\_aux\_init called with an "}// & 1518 \textcolor{stringliteral}{"associated control structure."}) 1519 \textcolor{keywordflow}{return} 1520 \textcolor{keywordflow}{else} 1521 \textcolor{keyword}{allocate}(cs) 1522 \textcolor{keywordflow}{ endif} 1523 1524 cs%diag => diag 1525 cs%Time => time 1526 1527 \textcolor{comment}{! Set default, read and log parameters} 1528 \textcolor{keyword}{call }log\_version(param\_file, mdl, version, & 1529 \textcolor{stringliteral}{"The following parameters are used for auxiliary diabatic processes."}) 1530 1531 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"ENABLE\_THERMODYNAMICS"}, use\_temperature, & 1532 \textcolor{stringliteral}{"If true, temperature and salinity are used as state "}//& 1533 \textcolor{stringliteral}{"variables."}, default=.true.) 1534 1535 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"RECLAIM\_FRAZIL"}, cs%reclaim\_frazil, & 1536 \textcolor{stringliteral}{"If true, try to use any frazil heat deficit to cool any "}//& 1537 \textcolor{stringliteral}{"overlying layers down to the freezing point, thereby "}//& 1538 \textcolor{stringliteral}{"avoiding the creation of thin ice when the SST is above "}//& 1539 \textcolor{stringliteral}{"the freezing point."}, default=.true.) 1540 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"PRESSURE\_DEPENDENT\_FRAZIL"}, & 1541 cs%pressure\_dependent\_frazil, & 1542 \textcolor{stringliteral}{"If true, use a pressure dependent freezing temperature "}//& 1543 \textcolor{stringliteral}{"when making frazil. The default is false, which will be "}//& 1544 \textcolor{stringliteral}{"faster but is inappropriate with ice-shelf cavities."}, & 1545 default=.false.) 1546 1547 \textcolor{keywordflow}{if} (use\_epbl) \textcolor{keywordflow}{then} 1548 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"IGNORE\_FLUXES\_OVER\_LAND"}, cs%ignore\_fluxes\_over\_land,& 1549 \textcolor{stringliteral}{"If true, the model does not check if fluxes are being applied "}//& 1550 \textcolor{stringliteral}{"over land points. This is needed when the ocean is coupled "}//& 1551 \textcolor{stringliteral}{"with ice shelves and sea ice, since the sea ice mask needs to "}//& 1552 \textcolor{stringliteral}{"be different than the ocean mask to avoid sea ice formation "}//& 1553 \textcolor{stringliteral}{"under ice shelves. This flag only works when use\_ePBL = True."}, default=.false.) 1554 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DO\_RIVERMIX"}, cs%do\_rivermix, & 1555 \textcolor{stringliteral}{"If true, apply additional mixing wherever there is "}//& 1556 \textcolor{stringliteral}{"runoff, so that it is mixed down to RIVERMIX\_DEPTH "}//& 1557 \textcolor{stringliteral}{"if the ocean is that deep."}, default=.false.) 1558 \textcolor{keywordflow}{if} (cs%do\_rivermix) & 1559 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"RIVERMIX\_DEPTH"}, cs%rivermix\_depth, & 1560 \textcolor{stringliteral}{"The depth to which rivers are mixed if DO\_RIVERMIX is "}//& 1561 \textcolor{stringliteral}{"defined."}, units=\textcolor{stringliteral}{"m"}, default=0.0, scale=us%m\_to\_Z) 1562 \textcolor{keywordflow}{else} 1563 cs%do\_rivermix = .false. ; cs%rivermix\_depth = 0.0 ; cs%ignore\_fluxes\_over\_land = .false. 1564 \textcolor{keywordflow}{ endif} 1565 1566 \textcolor{keywordflow}{if} (gv%nkml == 0) \textcolor{keywordflow}{then} 1567 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_RIVER\_HEAT\_CONTENT"}, cs%use\_river\_heat\_content, & 1568 \textcolor{stringliteral}{"If true, use the fluxes%runoff\_Hflx field to set the "}//& 1569 \textcolor{stringliteral}{"heat carried by runoff, instead of using SST*CP*liq\_runoff."}, & 1570 default=.false.) 1571 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_CALVING\_HEAT\_CONTENT"}, cs%use\_calving\_heat\_content, & 1572 \textcolor{stringliteral}{"If true, use the fluxes%calving\_Hflx field to set the "}//& 1573 \textcolor{stringliteral}{"heat carried by runoff, instead of using SST*CP*froz\_runoff."}, & 1574 default=.false.) 1575 \textcolor{keywordflow}{else} 1576 cs%use\_river\_heat\_content = .false. 1577 cs%use\_calving\_heat\_content = .false. 1578 \textcolor{keywordflow}{ endif} 1579 1580 \textcolor{keywordflow}{if} (usealealgorithm) \textcolor{keywordflow}{then} 1581 cs%id\_createdH = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'},\textcolor{stringliteral}{"created\_H"},diag%axesT1, & 1582 time, \textcolor{stringliteral}{"The volume flux added to stop the ocean from drying out and becoming negative in depth"}, & 1583 \textcolor{stringliteral}{"m s-1"}, conversion=gv%H\_to\_m*us%s\_to\_T) 1584 \textcolor{keywordflow}{if} (cs%id\_createdH>0) \textcolor{keyword}{allocate}(cs%createdH(isd:ied,jsd:jed)) 1585 1586 \textcolor{comment}{! diagnostic for heating of a grid cell from convergence of SW heat into the cell} 1587 cs%id\_penSW\_diag = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'rsdoabsorb'}, & 1588 diag%axesTL, time, \textcolor{stringliteral}{'Convergence of Penetrative Shortwave Flux in Sea Water Layer'},& 1589 \textcolor{stringliteral}{'W m-2'}, conversion=us%QRZ\_T\_to\_W\_m2, & 1590 standard\_name=\textcolor{stringliteral}{'net\_rate\_of\_absorption\_of\_shortwave\_energy\_in\_ocean\_layer'}, v\_extensive=.true.) 1591 1592 \textcolor{comment}{! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces)} 1593 \textcolor{comment}{! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock).} 1594 cs%id\_penSWflux\_diag = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'rsdo'}, & 1595 diag%axesTi, time, \textcolor{stringliteral}{'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface'},& 1596 \textcolor{stringliteral}{'W m-2'}, conversion=us%QRZ\_T\_to\_W\_m2, standard\_name=\textcolor{stringliteral}{'downwelling\_shortwave\_flux\_in\_sea\_water'}) 1597 1598 \textcolor{comment}{! need both arrays for the SW diagnostics (one for flux, one for convergence)} 1599 \textcolor{keywordflow}{if} (cs%id\_penSW\_diag>0 .or. cs%id\_penSWflux\_diag>0) \textcolor{keywordflow}{then} 1600 \textcolor{keyword}{allocate}(cs%penSW\_diag(isd:ied,jsd:jed,nz)) ; cs%penSW\_diag(:,:,:) = 0.0 1601 \textcolor{keyword}{allocate}(cs%penSWflux\_diag(isd:ied,jsd:jed,nz+1)) ; cs%penSWflux\_diag(:,:,:) = 0.0 1602 \textcolor{keywordflow}{ endif} 1603 1604 \textcolor{comment}{! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface)} 1605 cs%id\_nonpenSW\_diag = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'nonpenSW'}, & 1606 diag%axesT1, time, & 1607 \textcolor{stringliteral}{'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)'},& 1608 \textcolor{stringliteral}{'W m-2'}, conversion=us%QRZ\_T\_to\_W\_m2, & 1609 standard\_name=\textcolor{stringliteral}{'nondownwelling\_shortwave\_flux\_in\_sea\_water'}) 1610 \textcolor{keywordflow}{if} (cs%id\_nonpenSW\_diag > 0) \textcolor{keywordflow}{then} 1611 \textcolor{keyword}{allocate}(cs%nonpenSW\_diag(isd:ied,jsd:jed)) ; cs%nonpenSW\_diag(:,:) = 0.0 1612 \textcolor{keywordflow}{ endif} 1613 \textcolor{keywordflow}{ endif} 1614 1615 \textcolor{keywordflow}{if} (use\_temperature) \textcolor{keywordflow}{then} 1616 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"VAR\_PEN\_SW"}, cs%var\_pen\_sw, & 1617 \textcolor{stringliteral}{"If true, use one of the CHL\_A schemes specified by "}//& 1618 \textcolor{stringliteral}{"OPACITY\_SCHEME to determine the e-folding depth of "}//& 1619 \textcolor{stringliteral}{"incoming short wave radiation."}, default=.false.) 1620 \textcolor{keywordflow}{if} (cs%var\_pen\_sw) \textcolor{keywordflow}{then} 1621 1622 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CHL\_FROM\_FILE"}, cs%chl\_from\_file, & 1623 \textcolor{stringliteral}{"If true, chl\_a is read from a file."}, default=.true.) 1624 \textcolor{keywordflow}{if} (cs%chl\_from\_file) \textcolor{keywordflow}{then} 1625 \textcolor{keyword}{call }time\_interp\_external\_init() 1626 1627 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"INPUTDIR"}, inputdir, default=\textcolor{stringliteral}{"."}) 1628 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CHL\_FILE"}, chl\_file, & 1629 \textcolor{stringliteral}{"CHL\_FILE is the file containing chl\_a concentrations in "}//& 1630 \textcolor{stringliteral}{"the variable CHL\_A. It is used when VAR\_PEN\_SW and "}//& 1631 \textcolor{stringliteral}{"CHL\_FROM\_FILE are true."}, fail\_if\_missing=.true.) 1632 chl\_filename = trim(slasher(inputdir))//trim(chl\_file) 1633 \textcolor{keyword}{call }log\_param(param\_file, mdl, \textcolor{stringliteral}{"INPUTDIR/CHL\_FILE"}, chl\_filename) 1634 \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CHL\_VARNAME"}, chl\_varname, & 1635 \textcolor{stringliteral}{"Name of CHL\_A variable in CHL\_FILE."}, default=\textcolor{stringliteral}{'CHL\_A'}) 1636 cs%sbc\_chl = init\_external\_field(chl\_filename, trim(chl\_varname), domain=g%Domain%mpp\_domain) 1637 \textcolor{keywordflow}{ endif} 1638 1639 cs%id\_chl = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'Chl\_opac'}, diag%axesT1, time, & 1640 \textcolor{stringliteral}{'Surface chlorophyll A concentration used to find opacity'}, \textcolor{stringliteral}{'mg m-3'}) 1641 \textcolor{keywordflow}{ endif} 1642 \textcolor{keywordflow}{ endif} 1643 1644 id\_clock\_uv\_at\_h = cpu\_clock\_id(\textcolor{stringliteral}{'(Ocean find\_uv\_at\_h)'}, grain=clock\_routine) 1645 id\_clock\_frazil = cpu\_clock\_id(\textcolor{stringliteral}{'(Ocean frazil)'}, grain=clock\_routine) 1646 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_a43ac98433f0f2f0e6cf4f8a3c9622ec4}\label{namespacemom__diabatic__aux_a43ac98433f0f2f0e6cf4f8a3c9622ec4}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!diagnosemldbydensitydifference@{diagnosemldbydensitydifference}} \index{diagnosemldbydensitydifference@{diagnosemldbydensitydifference}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{diagnosemldbydensitydifference()}{diagnosemldbydensitydifference()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::diagnosemldbydensitydifference (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{id\+\_\+\+M\+LD, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h, }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv, }\item[{real, intent(in)}]{density\+Diff, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{type(diag\+\_\+ctrl), pointer}]{diag\+Ptr, }\item[{integer, intent(in), optional}]{id\+\_\+\+N2sub\+ML, }\item[{integer, intent(in), optional}]{id\+\_\+\+M\+L\+Dsq, }\item[{real, intent(in), optional}]{dz\+\_\+sub\+ML }\end{DoxyParamCaption})} Diagnose a mixed layer depth (M\+LD) determined by a given density difference with the surface. This routine is appropriate in M\+O\+M\+\_\+diabatic\+\_\+driver due to its position within the time stepping. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em g} & Grid type\\ \hline \mbox{\tt in} & {\em gv} & ocean vertical grid structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in} & {\em id\+\_\+mld} & Handle (ID) of M\+LD diagnostic\\ \hline \mbox{\tt in} & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\ \hline \mbox{\tt in} & {\em densitydiff} & Density difference to determine M\+LD \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]}\\ \hline & {\em diagptr} & Diagnostics structure\\ \hline \mbox{\tt in} & {\em id\+\_\+n2subml} & Optional handle (ID) of sub\+ML stratification\\ \hline \mbox{\tt in} & {\em id\+\_\+mldsq} & Optional handle (ID) of squared M\+LD\\ \hline \mbox{\tt in} & {\em dz\+\_\+subml} & The distance over which to calculate N2sub\+ML or 50 m if missing \mbox{[}Z $\sim$$>$ m\mbox{]} \\ \hline \end{DoxyParams} Definition at line 600 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 600 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< Grid type} 601 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< ocean vertical grid structure} 602 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 603 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: id\_mld\textcolor{comment}{ !< Handle (ID) of MLD diagnostic} 604 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 605 \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< Layer thickness [H ~> m or kg m-2]} 606 \textcolor{keywordtype}{type}(thermo\_var\_ptrs), \textcolor{keywordtype}{intent(in)} :: tv\textcolor{comment}{ !< Structure containing pointers to any} 607 \textcolor{comment}{ !! available thermodynamic fields.} 608 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: densitydiff\textcolor{comment}{ !< Density difference to determine MLD [R ~> kg m-3]} 609 \textcolor{keywordtype}{type}(diag\_ctrl), \textcolor{keywordtype}{pointer} :: diagptr\textcolor{comment}{ !< Diagnostics structure} 610 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: id\_n2subml\textcolor{comment}{ !< Optional handle (ID) of subML stratification} 611 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: id\_mldsq\textcolor{comment}{ !< Optional handle (ID) of squared MLD} 612 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: dz\_subml\textcolor{comment}{ !< The distance over which to calculate N2subML} 613 \textcolor{comment}{ !! or 50 m if missing [Z ~> m]} 614 615 \textcolor{comment}{! Local variables} 616 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: deltarhoatkm1, deltarhoatk \textcolor{comment}{! Density differences [R ~> kg m-3].} 617 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: pref\_mld, pref\_n2 \textcolor{comment}{! Reference pressures [R L2 T-2 ~> Pa].} 618 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: h\_subml, dh\_n2 \textcolor{comment}{! Summed thicknesses used in N2 calculation [H ~> m].} 619 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: t\_subml, t\_deeper \textcolor{comment}{! Temperatures used in the N2 calculation [degC].} 620 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: s\_subml, s\_deeper \textcolor{comment}{! Salinities used in the N2 calculation [ppt].} 621 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: rho\_subml, rho\_deeper \textcolor{comment}{! Densities used in the N2 calculation [R ~> kg m-3].} 622 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: dk, dkm1 \textcolor{comment}{! Depths [Z ~> m].} 623 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: rhosurf \textcolor{comment}{! Density used in finding the mixedlayer depth [R ~> kg m-3].} 624 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G))} :: mld \textcolor{comment}{! Diagnosed mixed layer depth [Z ~> m].} 625 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G))} :: submln2 \textcolor{comment}{! Diagnosed stratification below ML [T-2 ~> s-2].} 626 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G))} :: mld2 \textcolor{comment}{! Diagnosed MLD^2 [Z2 ~> m2].} 627 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: n2\_region\_set \textcolor{comment}{! If true, all necessary values for calculating N2} 628 \textcolor{comment}{! have been stored already.} 629 \textcolor{keywordtype}{real} :: ge\_rho0 \textcolor{comment}{! The gravitational acceleration divided by a mean density [Z T-2 R-1 ~> m4 s-2 kg-1].} 630 \textcolor{keywordtype}{real} :: dh\_subml \textcolor{comment}{! Depth below ML over which to diagnose stratification [H ~> m].} 631 \textcolor{keywordtype}{real} :: afac \textcolor{comment}{! A nondimensional factor [nondim]} 632 \textcolor{keywordtype}{real} :: ddrho \textcolor{comment}{! A density difference [R ~> kg m-3]} 633 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(2)} :: eosdom \textcolor{comment}{! The i-computational domain for the equation of state} 634 \textcolor{keywordtype}{integer} :: i, j, is, ie, js, je, k, nz, id\_n2, id\_sq 635 636 id\_n2 = -1 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{PRESENT}(id\_n2subml)) id\_n2 = id\_n2subml 637 638 id\_sq = -1 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{PRESENT}(id\_mldsq)) id\_sq = id\_mldsq 639 640 ge\_rho0 = us%L\_to\_Z**2*gv%g\_Earth / (gv%Rho0) 641 dh\_subml = 50.*gv%m\_to\_H ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(dz\_subml)) dh\_subml = gv%Z\_to\_H*dz\_subml 642 643 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke 644 645 pref\_mld(:) = 0.0 646 eosdom(:) = eos\_domain(g%HI) 647 \textcolor{keywordflow}{do} j=js,je 648 \textcolor{keywordflow}{do} i=is,ie ; dk(i) = 0.5 * h(i,j,1) * gv%H\_to\_Z ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! Depth of center of surface layer} 649 \textcolor{keyword}{call }calculate\_density(tv%T(:,j,1), tv%S(:,j,1), pref\_mld, rhosurf, tv%eqn\_of\_state, eosdom) 650 \textcolor{keywordflow}{do} i=is,ie 651 deltarhoatk(i) = 0. 652 mld(i,j) = 0. 653 \textcolor{keywordflow}{if} (id\_n2>0) \textcolor{keywordflow}{then} 654 submln2(i,j) = 0.0 655 h\_subml(i) = h(i,j,1) ; dh\_n2(i) = 0.0 656 t\_subml(i) = 0.0 ; s\_subml(i) = 0.0 ; t\_deeper(i) = 0.0 ; s\_deeper(i) = 0.0 657 n2\_region\_set(i) = (g%mask2dT(i,j)<0.5) \textcolor{comment}{! Only need to work on ocean points.} 658 \textcolor{keywordflow}{ endif} 659 \textcolor{keywordflow}{ enddo} 660 \textcolor{keywordflow}{do} k=2,nz 661 \textcolor{keywordflow}{do} i=is,ie 662 dkm1(i) = dk(i) \textcolor{comment}{! Depth of center of layer K-1} 663 dk(i) = dk(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * gv%H\_to\_Z \textcolor{comment}{! Depth of center of layer K} 664 \textcolor{keywordflow}{ enddo} 665 666 \textcolor{comment}{! Prepare to calculate stratification, N2, immediately below the mixed layer by finding} 667 \textcolor{comment}{! the cells that extend over at least dz\_subML.} 668 \textcolor{keywordflow}{if} (id\_n2>0) \textcolor{keywordflow}{then} 669 \textcolor{keywordflow}{do} i=is,ie 670 \textcolor{keywordflow}{if} (mld(i,j)==0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! Still in the mixed layer.} 671 h\_subml(i) = h\_subml(i) + h(i,j,k) 672 \textcolor{keywordflow}{elseif} (.not.n2\_region\_set(i)) \textcolor{keywordflow}{then} \textcolor{comment}{! This block is below the mixed layer, but N2 has not been found yet.} 673 \textcolor{keywordflow}{if} (dh\_n2(i)==0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! Record the temperature, salinity, pressure, immediately below the ML} 674 t\_subml(i) = tv%T(i,j,k) ; s\_subml(i) = tv%S(i,j,k) 675 h\_subml(i) = h\_subml(i) + 0.5 * h(i,j,k) \textcolor{comment}{! Start midway through this layer.} 676 dh\_n2(i) = 0.5 * h(i,j,k) 677 \textcolor{keywordflow}{elseif} (dh\_n2(i) + h(i,j,k) < dh\_subml) \textcolor{keywordflow}{then} 678 dh\_n2(i) = dh\_n2(i) + h(i,j,k) 679 \textcolor{keywordflow}{else} \textcolor{comment}{! This layer includes the base of the region where N2 is calculated.} 680 t\_deeper(i) = tv%T(i,j,k) ; s\_deeper(i) = tv%S(i,j,k) 681 dh\_n2(i) = dh\_n2(i) + 0.5 * h(i,j,k) 682 n2\_region\_set(i) = .true. 683 \textcolor{keywordflow}{ endif} 684 \textcolor{keywordflow}{ endif} 685 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! i-loop} 686 \textcolor{keywordflow}{ endif} \textcolor{comment}{! id\_N2>0} 687 688 \textcolor{comment}{! Mixed-layer depth, using sigma-0 (surface reference pressure)} 689 \textcolor{keywordflow}{do} i=is,ie ; deltarhoatkm1(i) = deltarhoatk(i) ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! Store value from previous iteration of K} 690 \textcolor{keyword}{call }calculate\_density(tv%T(:,j,k), tv%S(:,j,k), pref\_mld, deltarhoatk, tv%eqn\_of\_state, eosdom) 691 \textcolor{keywordflow}{do} i = is, ie 692 deltarhoatk(i) = deltarhoatk(i) - rhosurf(i) \textcolor{comment}{! Density difference between layer K and surface} 693 ddrho = deltarhoatk(i) - deltarhoatkm1(i) 694 \textcolor{keywordflow}{if} ((mld(i,j)==0.) .and. (ddrho>0.) .and. & 695 (deltarhoatkm1(i)=densitydiff)) \textcolor{keywordflow}{then} 696 afac = ( densitydiff - deltarhoatkm1(i) ) / ddrho 697 mld(i,j) = dk(i) * afac + dkm1(i) * (1. - afac) 698 \textcolor{keywordflow}{ endif} 699 \textcolor{keywordflow}{if} (id\_sq > 0) mld2(i,j) = mld(i,j)**2 700 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! i-loop} 701 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! k-loop} 702 \textcolor{keywordflow}{do} i=is,ie 703 \textcolor{keywordflow}{if} ((mld(i,j)==0.) .and. (deltarhoatk(i)0) \textcolor{keywordflow}{then} \textcolor{comment}{! Now actually calculate stratification, N2, below the mixed layer.} 707 \textcolor{keywordflow}{do} i=is,ie ; pref\_n2(i) = (gv%g\_Earth * gv%H\_to\_RZ) * (h\_subml(i) + 0.5*dh\_n2(i)) ;\textcolor{keywordflow}{ enddo} 708 \textcolor{comment}{! if ((.not.N2\_region\_set(i)) .and. (dH\_N2(i) > 0.5*dH\_subML)) then} 709 \textcolor{comment}{! ! Use whatever stratification we can, measured over whatever distance is available?} 710 \textcolor{comment}{! T\_deeper(i) = tv%T(i,j,nz) ; S\_deeper(i) = tv%S(i,j,nz)} 711 \textcolor{comment}{! N2\_region\_set(i) = .true.} 712 \textcolor{comment}{! endif} 713 \textcolor{keyword}{call }calculate\_density(t\_subml, s\_subml, pref\_n2, rho\_subml, tv%eqn\_of\_state, eosdom) 714 \textcolor{keyword}{call }calculate\_density(t\_deeper, s\_deeper, pref\_n2, rho\_deeper, tv%eqn\_of\_state, eosdom) 715 \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} ((g%mask2dT(i,j)>0.5) .and. n2\_region\_set(i)) \textcolor{keywordflow}{then} 716 submln2(i,j) = ge\_rho0 * (rho\_deeper(i) - rho\_subml(i)) / (gv%H\_to\_z * dh\_n2(i)) 717 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} 718 \textcolor{keywordflow}{ endif} 719 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! j-loop} 720 721 \textcolor{keywordflow}{if} (id\_mld > 0) \textcolor{keyword}{call }post\_data(id\_mld, mld, diagptr) 722 \textcolor{keywordflow}{if} (id\_n2 > 0) \textcolor{keyword}{call }post\_data(id\_n2, submln2, diagptr) 723 \textcolor{keywordflow}{if} (id\_sq > 0) \textcolor{keyword}{call }post\_data(id\_sq, mld2, diagptr) 724 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_a695474893b92e3e3418a2bed871a5d7e}\label{namespacemom__diabatic__aux_a695474893b92e3e3418a2bed871a5d7e}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!diagnosemldbyenergy@{diagnosemldbyenergy}} \index{diagnosemldbyenergy@{diagnosemldbyenergy}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{diagnosemldbyenergy()}{diagnosemldbyenergy()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::diagnosemldbyenergy (\begin{DoxyParamCaption}\item[{integer, dimension(3), intent(in)}]{id\+\_\+\+M\+LD, }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{h, }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, dimension(3), intent(in)}]{Mixing\+\_\+\+Energy, }\item[{type(diag\+\_\+ctrl), pointer}]{diag\+Ptr }\end{DoxyParamCaption})} Diagnose a mixed layer depth (M\+LD) determined by the depth a given energy value would mix. This routine is appropriate in M\+O\+M\+\_\+diabatic\+\_\+driver due to its position within the time stepping. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em id\+\_\+mld} & Energy output diag I\+Ds\\ \hline \mbox{\tt in} & {\em g} & Grid type\\ \hline \mbox{\tt in} & {\em gv} & ocean vertical grid structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in} & {\em mixing\+\_\+energy} & Energy values for up to 3 M\+L\+Ds\\ \hline \mbox{\tt in} & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\ \hline & {\em diagptr} & Diagnostics structure \\ \hline \end{DoxyParams} Definition at line 730 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 730 \textcolor{comment}{! Author: Brandon Reichl} 731 \textcolor{comment}{! Date: October 2, 2020} 732 \textcolor{comment}{! //} 733 \textcolor{comment}{! *Note that gravity is assumed constant everywhere and divided out of all calculations.} 734 \textcolor{comment}{!} 735 \textcolor{comment}{! This code has been written to step through the columns layer by layer, summing the PE} 736 \textcolor{comment}{! change inferred by mixing the layer with all layers above. When the change exceeds a} 737 \textcolor{comment}{! threshold (determined by input array Mixing\_Energy), the code needs to solve for how far} 738 \textcolor{comment}{! into this layer the threshold PE change occurs (assuming constant density layers).} 739 \textcolor{comment}{! This is expressed here via solving the function F(X) = 0 where:} 740 \textcolor{comment}{! F(X) = 0.5 * ( Ca*X^3/(D1+X) + Cb*X^2/(D1+X) + Cc*X/(D1+X) + Dc/(D1+X)} 741 \textcolor{comment}{! + Ca2*X^2 + Cb2*X + Cc2)} 742 \textcolor{comment}{! where all coefficients are determined by the previous mixed layer depth, the} 743 \textcolor{comment}{! density of the previous mixed layer, the present layer thickness, and the present} 744 \textcolor{comment}{! layer density. This equation is worked out by computing the total PE assuming constant} 745 \textcolor{comment}{! density in the mixed layer as well as in the remaining part of the present layer that is} 746 \textcolor{comment}{! not mixed.} 747 \textcolor{comment}{! To solve for X in this equation a Newton's method iteration is employed, which} 748 \textcolor{comment}{! converges extremely quickly (usually 1 guess) since this equation turns out being rather} 749 \textcolor{comment}{! lienar for PE change with increasing X.} 750 \textcolor{comment}{! Input parameters:} 751 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(3)}, \textcolor{keywordtype}{intent(in)} :: id\_mld\textcolor{comment}{ !< Energy output diag IDs} 752 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< Grid type} 753 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< ocean vertical grid structure} 754 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 755 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(3)}, \textcolor{keywordtype}{intent(in)} :: mixing\_energy\textcolor{comment}{ !< Energy values for up to 3 MLDs} 756 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 757 \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< Layer thickness [H ~> m or kg m-2]} 758 \textcolor{keywordtype}{type}(thermo\_var\_ptrs), \textcolor{keywordtype}{intent(in)} :: tv\textcolor{comment}{ !< Structure containing pointers to any} 759 \textcolor{comment}{ !! available thermodynamic fields.} 760 \textcolor{keywordtype}{type}(diag\_ctrl), \textcolor{keywordtype}{pointer} :: diagptr\textcolor{comment}{ !< Diagnostics structure} 761 762 \textcolor{comment}{! Local variables} 763 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G),3)} :: mld \textcolor{comment}{! Diagnosed mixed layer depth [Z ~> m].} 764 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: z\_l, z\_u, dz, rho\_c, pref\_mld 765 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(3)} :: pe\_threshold 766 767 \textcolor{keywordtype}{real} :: ig, e\_g 768 \textcolor{keywordtype}{real} :: pe\_threshold\_fraction, pe, pe\_mixed, pe\_mixed\_tst 769 \textcolor{keywordtype}{real} :: rhodz\_ml, h\_ml, rhodz\_ml\_tst, h\_ml\_tst 770 \textcolor{keywordtype}{real} :: rho\_ml 771 772 \textcolor{keywordtype}{real} :: r1, d1, r2, d2 773 \textcolor{keywordtype}{real} :: ca, cb,d ,cc, cd, ca2, cb2, c, cc2 774 \textcolor{keywordtype}{real} :: gx, gpx, hx, ihx, hpx, ix, ipx, fgx, fpx, x, x2 775 776 \textcolor{keywordtype}{integer} :: it, im 777 \textcolor{keywordtype}{integer} :: i, j, is, ie, js, je, k, nz 778 779 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke 780 781 pref\_mld(:) = 0.0 782 mld(:,:,:) = 0.0 783 pe\_threshold\_fraction = 1.e-4 \textcolor{comment}{!Fixed threshold of 0.01%, could be runtime.} 784 785 \textcolor{keywordflow}{do} im=1,3 786 pe\_threshold(im) = mixing\_energy(im)/gv%g\_earth 787 \textcolor{keywordflow}{ enddo} 788 789 \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is,ie 790 \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.0) \textcolor{keywordflow}{then} 791 792 \textcolor{keyword}{call }calculate\_density(tv%T(i,j,:), tv%S(i,j,:), pref\_mld, rho\_c, 1, nz, & 793 tv%eqn\_of\_state, scale=us%kg\_m3\_to\_R) 794 795 \textcolor{keywordflow}{do} k=1,nz 796 dz(k) = h(i,j,k) * gv%H\_to\_Z 797 \textcolor{keywordflow}{ enddo} 798 z\_u(1) = 0.0 799 z\_l(1) = -dz(1) 800 \textcolor{keywordflow}{do} k=2,nz 801 z\_u(k) = z\_l(k-1) 802 z\_l(k) = z\_l(k-1)-dz(k) 803 \textcolor{keywordflow}{ enddo} 804 805 \textcolor{keywordflow}{do} im=1,3 806 807 \textcolor{comment}{! Initialize these for each columnwise calculation} 808 pe = 0.0 809 rhodz\_ml = 0.0 810 h\_ml = 0.0 811 rhodz\_ml\_tst = 0.0 812 h\_ml\_tst = 0.0 813 pe\_mixed = 0.0 814 815 \textcolor{keywordflow}{do} k=1,nz 816 817 \textcolor{comment}{! This is the unmixed PE cummulative sum from top down} 818 pe = pe + 0.5*rho\_c(k)*(z\_u(k)**2-z\_l(k)**2) 819 820 \textcolor{comment}{! This is the depth and integral of density} 821 h\_ml\_tst = h\_ml + dz(k) 822 rhodz\_ml\_tst = rhodz\_ml + rho\_c(k)*dz(k) 823 824 \textcolor{comment}{! The average density assuming all layers including this were mixed} 825 rho\_ml = rhodz\_ml\_tst/h\_ml\_tst 826 827 \textcolor{comment}{! The PE assuming all layers including this were mixed} 828 \textcolor{comment}{! Note that 0. could be replaced with "Surface", which doesn't have to be 0} 829 \textcolor{comment}{! but 0 is a good reference value.} 830 pe\_mixed\_tst = 0.5*rho\_ml*(0.**2-(0.-h\_ml\_tst)**2) 831 832 \textcolor{comment}{! Check if we supplied enough energy to mix to this layer} 833 \textcolor{keywordflow}{if} (pe\_mixed\_tst-pe<=pe\_threshold(im)) \textcolor{keywordflow}{then} 834 h\_ml = h\_ml\_tst 835 rhodz\_ml = rhodz\_ml\_tst 836 837 \textcolor{keywordflow}{else} \textcolor{comment}{! If not, we need to solve where the energy ran out} 838 \textcolor{comment}{! This will be done with a Newton's method iteration:} 839 840 r1 = rhodz\_ml/h\_ml \textcolor{comment}{! The density of the mixed layer (not including this layer)} 841 d1 = h\_ml \textcolor{comment}{! The thickness of the mixed layer (not including this layer)} 842 r2 = rho\_c(k) \textcolor{comment}{! The density of this layer} 843 d2 = dz(k) \textcolor{comment}{! The thickness of this layer} 844 845 \textcolor{comment}{! This block could be used to calculate the function coefficients if} 846 \textcolor{comment}{! we don't reference all values to a surface designated as z=0} 847 \textcolor{comment}{! S = Surface} 848 \textcolor{comment}{! Ca = -(R2)} 849 \textcolor{comment}{! Cb = -( (R1*D1) + R2*(2.*D1-2.*S) )} 850 \textcolor{comment}{! D = D1**2. - 2.*D1*S} 851 \textcolor{comment}{! Cc = -( R1*D1*(2.*D1-2.*S) + R2*D )} 852 \textcolor{comment}{! Cd = -(R1*D1*D)} 853 \textcolor{comment}{! Ca2 = R2} 854 \textcolor{comment}{! Cb2 = R2*(2*D1-2*S)} 855 \textcolor{comment}{! C = S**2 + D2**2 + D1**2 - 2*D1*S - 2.*D2*S +2.*D1*D2} 856 \textcolor{comment}{! Cc2 = R2*(D+S**2-C)} 857 \textcolor{comment}{!} 858 \textcolor{comment}{! If the surface is S = 0, it simplifies to:} 859 ca = -(r2) 860 cb = -( r1*d1 + r2*(2.*d1) ) 861 d = d1**2. 862 cc = -( r1*d1*(2*d1) + r2*d ) 863 cd = -r1*d1*d 864 ca2 = r2 865 cb2 = r2*(2.*d1) 866 c = d2**2. + d1**2. + 2.*d1*d2 867 cc2 = r2*(d-c) 868 869 \textcolor{comment}{! First guess for an iteration using Newton's method} 870 x = dz(k)*0.5 871 872 it=0 873 \textcolor{keywordflow}{do} \textcolor{keywordflow}{while}(it<10)\textcolor{comment}{!We can iterate up to 10 times} 874 \textcolor{comment}{! We are trying to solve the function:} 875 \textcolor{comment}{! F(x) = G(x)/H(x)+I(x)} 876 \textcolor{comment}{! for where F(x) = PE+PE\_threshold, or equivalently for where} 877 \textcolor{comment}{! F(x) = G(x)/H(x)+I(x) - (PE+PE\_threshold) = 0} 878 \textcolor{comment}{! We also need the derivative of this function for the Newton's method iteration} 879 \textcolor{comment}{! F'(x) = (G'(x)H(x)-G(x)H'(x))/H(x)^2 + I'(x)} 880 \textcolor{comment}{! G and its derivative} 881 gx = 0.5*(ca*x**3+cb*x**2+cc*x+cd) 882 gpx = 0.5*(3*ca*x**2+2*cb*x+cc) 883 \textcolor{comment}{! H, its inverse, and its derivative} 884 hx = (d1+x) 885 ihx = 1./hx 886 hpx = 1. 887 \textcolor{comment}{! I and its derivative} 888 ix = 0.5*(ca2*x**2. + cb2*x + cc2) 889 ipx = 0.5*(2.*ca2*x+cb2) 890 891 \textcolor{comment}{! The Function and its derivative:} 892 pe\_mixed = gx*ihx+ix 893 fgx = (pe\_mixed-(pe+pe\_threshold(im))) 894 fpx = (gpx*hx-hpx*gx)*ihx**2+ipx 895 896 \textcolor{comment}{! Check if our solution is within the threshold bounds, if not update} 897 \textcolor{comment}{! using Newton's method. This appears to converge almost always in} 898 \textcolor{comment}{! one step because the function is very close to linear in most applications.} 899 \textcolor{keywordflow}{if} (abs(fgx)>pe\_threshold(im)*pe\_threshold\_fraction) \textcolor{keywordflow}{then} 900 x2 = x - fgx/fpx 901 it = it + 1 902 \textcolor{keywordflow}{if} (x2<0. .or. x2>dz(k)) \textcolor{keywordflow}{then} 903 \textcolor{comment}{! The iteration seems to be robust, but we need to do something *if*} 904 \textcolor{comment}{! things go wrong... How should we treat failed iteration?} 905 \textcolor{comment}{! Present solution: Stop trying to compute and just say we can't mix this layer.} 906 x=0 907 \textcolor{keywordflow}{exit} 908 \textcolor{keywordflow}{else} 909 x = x2 910 \textcolor{keywordflow}{ endif} 911 \textcolor{keywordflow}{else} 912 exit\textcolor{comment}{! Quit the iteration} 913 \textcolor{keywordflow}{ endif} 914 \textcolor{keywordflow}{ enddo} 915 h\_ml = h\_ml+x 916 exit\textcolor{comment}{! Quit looping through the column} 917 \textcolor{keywordflow}{ endif} 918 \textcolor{keywordflow}{ enddo} 919 mld(i,j,im) = h\_ml 920 \textcolor{keywordflow}{ enddo} 921 \textcolor{keywordflow}{else} 922 mld(i,j,:) = 0.0 923 \textcolor{keywordflow}{ endif} 924 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 925 926 \textcolor{keywordflow}{if} (id\_mld(1) > 0) \textcolor{keyword}{call }post\_data(id\_mld(1), mld(:,:,1), diagptr) 927 \textcolor{keywordflow}{if} (id\_mld(2) > 0) \textcolor{keyword}{call }post\_data(id\_mld(2), mld(:,:,2), diagptr) 928 \textcolor{keywordflow}{if} (id\_mld(3) > 0) \textcolor{keyword}{call }post\_data(id\_mld(3), mld(:,:,3), diagptr) 929 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_aa9248acf0b1cf246a79092077bcb0b46}\label{namespacemom__diabatic__aux_aa9248acf0b1cf246a79092077bcb0b46}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!differential\+\_\+diffuse\+\_\+t\+\_\+s@{differential\+\_\+diffuse\+\_\+t\+\_\+s}} \index{differential\+\_\+diffuse\+\_\+t\+\_\+s@{differential\+\_\+diffuse\+\_\+t\+\_\+s}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{differential\+\_\+diffuse\+\_\+t\+\_\+s()}{differential\_diffuse\_t\_s()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::differential\+\_\+diffuse\+\_\+t\+\_\+s (\begin{DoxyParamCaption}\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{h, }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(inout)}]{tv, }\item[{type(vertvisc\+\_\+type), intent(in)}]{visc, }\item[{real, intent(in)}]{dt, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV }\end{DoxyParamCaption})} This subroutine applies double diffusion to T \& S, assuming no diapycal mass fluxes, using a simple triadiagonal solver. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em g} & The ocean\textquotesingle{}s grid structure\\ \hline \mbox{\tt in} & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\ \hline \mbox{\tt in} & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in,out} & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\ \hline \mbox{\tt in} & {\em visc} & Structure containing vertical viscosities, bottom boundary layer properies, and related fields.\\ \hline \mbox{\tt in} & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}. \\ \hline \end{DoxyParams} Definition at line 227 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 227 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< The ocean's grid structure} 228 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< The ocean's vertical grid structure} 229 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 230 \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< Layer thicknesses [H ~> m or kg m-2]} 231 \textcolor{keywordtype}{type}(thermo\_var\_ptrs), \textcolor{keywordtype}{intent(inout)} :: tv\textcolor{comment}{ !< Structure containing pointers to any} 232 \textcolor{comment}{ !! available thermodynamic fields.} 233 \textcolor{keywordtype}{type}(vertvisc\_type), \textcolor{keywordtype}{intent(in)} :: visc\textcolor{comment}{ !< Structure containing vertical viscosities, bottom} 234 \textcolor{comment}{ !! boundary layer properies, and related fields.} 235 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: dt\textcolor{comment}{ !< Time increment [T ~> s].} 236 237 \textcolor{comment}{! local variables} 238 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: & 239 b1\_t, b1\_s, & \textcolor{comment}{! Variables used by the tridiagonal solvers of T & S [H ~> m or kg m-2].} 240 d1\_t, d1\_s \textcolor{comment}{! Variables used by the tridiagonal solvers [nondim].} 241 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZK\_(G))} :: & 242 c1\_t, c1\_s \textcolor{comment}{! Variables used by the tridiagonal solvers [H ~> m or kg m-2].} 243 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZK\_(G)+1)} :: & 244 mix\_t, mix\_s \textcolor{comment}{! Mixing distances in both directions across each interface [H ~> m or kg m-2].} 245 \textcolor{keywordtype}{real} :: h\_tr \textcolor{comment}{! h\_tr is h at tracer points with a tiny thickness} 246 \textcolor{comment}{! added to ensure positive definiteness [H ~> m or kg m-2].} 247 \textcolor{keywordtype}{real} :: h\_neglect \textcolor{comment}{! A thickness that is so small it is usually lost} 248 \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].} 249 \textcolor{keywordtype}{real} :: i\_h\_int \textcolor{comment}{! The inverse of the thickness associated with an} 250 \textcolor{comment}{! interface [H-1 ~> m-1 or m2 kg-1].} 251 \textcolor{keywordtype}{real} :: b\_denom\_t \textcolor{comment}{! The first term in the denominators for the expressions} 252 \textcolor{keywordtype}{real} :: b\_denom\_s \textcolor{comment}{! for b1\_T and b1\_S, both [H ~> m or kg m-2].} 253 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:,:)}, \textcolor{keywordtype}{pointer} :: t=>null(), s=>null() 254 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:,:)}, \textcolor{keywordtype}{pointer} :: kd\_t=>null(), kd\_s=>null() \textcolor{comment}{! Diffusivities [Z2 T-1 ~> m2 s-1].} 255 \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, nz 256 257 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke 258 h\_neglect = gv%H\_subroundoff 259 260 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(tv%T)) \textcolor{keyword}{call }mom\_error(fatal, & 261 \textcolor{stringliteral}{"differential\_diffuse\_T\_S: Called with an unassociated tv%T"}) 262 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(tv%S)) \textcolor{keyword}{call }mom\_error(fatal, & 263 \textcolor{stringliteral}{"differential\_diffuse\_T\_S: Called with an unassociated tv%S"}) 264 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(visc%Kd\_extra\_T)) \textcolor{keyword}{call }mom\_error(fatal, & 265 \textcolor{stringliteral}{"differential\_diffuse\_T\_S: Called with an unassociated visc%Kd\_extra\_T"}) 266 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(visc%Kd\_extra\_S)) \textcolor{keyword}{call }mom\_error(fatal, & 267 \textcolor{stringliteral}{"differential\_diffuse\_T\_S: Called with an unassociated visc%Kd\_extra\_S"}) 268 269 t => tv%T ; s => tv%S 270 kd\_t => visc%Kd\_extra\_T ; kd\_s => visc%Kd\_extra\_S 271 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,h,h\_neglect,dt,Kd\_T,Kd\_S,G,GV,T,S,nz) &} 272 \textcolor{comment}{!$OMP private(I\_h\_int,mix\_T,mix\_S,h\_tr,b1\_T,b1\_S, &} 273 \textcolor{comment}{!$OMP d1\_T,d1\_S,c1\_T,c1\_S,b\_denom\_T,b\_denom\_S)} 274 \textcolor{keywordflow}{do} j=js,je 275 \textcolor{keywordflow}{do} i=is,ie 276 i\_h\_int = 1.0 / (0.5 * (h(i,j,1) + h(i,j,2)) + h\_neglect) 277 mix\_t(i,2) = ((dt * kd\_t(i,j,2)) * gv%Z\_to\_H**2) * i\_h\_int 278 mix\_s(i,2) = ((dt * kd\_s(i,j,2)) * gv%Z\_to\_H**2) * i\_h\_int 279 280 h\_tr = h(i,j,1) + h\_neglect 281 b1\_t(i) = 1.0 / (h\_tr + mix\_t(i,2)) 282 b1\_s(i) = 1.0 / (h\_tr + mix\_s(i,2)) 283 d1\_t(i) = h\_tr * b1\_t(i) 284 d1\_s(i) = h\_tr * b1\_s(i) 285 t(i,j,1) = (b1\_t(i)*h\_tr)*t(i,j,1) 286 s(i,j,1) = (b1\_s(i)*h\_tr)*s(i,j,1) 287 \textcolor{keywordflow}{ enddo} 288 \textcolor{keywordflow}{do} k=2,nz-1 ; \textcolor{keywordflow}{do} i=is,ie 289 \textcolor{comment}{! Calculate the mixing across the interface below this layer.} 290 i\_h\_int = 1.0 / (0.5 * (h(i,j,k) + h(i,j,k+1)) + h\_neglect) 291 mix\_t(i,k+1) = ((dt * kd\_t(i,j,k+1)) * gv%Z\_to\_H**2) * i\_h\_int 292 mix\_s(i,k+1) = ((dt * kd\_s(i,j,k+1)) * gv%Z\_to\_H**2) * i\_h\_int 293 294 c1\_t(i,k) = mix\_t(i,k) * b1\_t(i) 295 c1\_s(i,k) = mix\_s(i,k) * b1\_s(i) 296 297 h\_tr = h(i,j,k) + h\_neglect 298 b\_denom\_t = h\_tr + d1\_t(i)*mix\_t(i,k) 299 b\_denom\_s = h\_tr + d1\_s(i)*mix\_s(i,k) 300 b1\_t(i) = 1.0 / (b\_denom\_t + mix\_t(i,k+1)) 301 b1\_s(i) = 1.0 / (b\_denom\_s + mix\_s(i,k+1)) 302 d1\_t(i) = b\_denom\_t * b1\_t(i) 303 d1\_s(i) = b\_denom\_s * b1\_s(i) 304 305 t(i,j,k) = b1\_t(i) * (h\_tr*t(i,j,k) + mix\_t(i,k)*t(i,j,k-1)) 306 s(i,j,k) = b1\_s(i) * (h\_tr*s(i,j,k) + mix\_s(i,k)*s(i,j,k-1)) 307 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 308 \textcolor{keywordflow}{do} i=is,ie 309 c1\_t(i,nz) = mix\_t(i,nz) * b1\_t(i) 310 c1\_s(i,nz) = mix\_s(i,nz) * b1\_s(i) 311 312 h\_tr = h(i,j,nz) + h\_neglect 313 b1\_t(i) = 1.0 / (h\_tr + d1\_t(i)*mix\_t(i,nz)) 314 b1\_s(i) = 1.0 / (h\_tr + d1\_s(i)*mix\_s(i,nz)) 315 316 t(i,j,nz) = b1\_t(i) * (h\_tr*t(i,j,nz) + mix\_t(i,nz)*t(i,j,nz-1)) 317 s(i,j,nz) = b1\_s(i) * (h\_tr*s(i,j,nz) + mix\_s(i,nz)*s(i,j,nz-1)) 318 \textcolor{keywordflow}{ enddo} 319 \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=is,ie 320 t(i,j,k) = t(i,j,k) + c1\_t(i,k+1)*t(i,j,k+1) 321 s(i,j,k) = s(i,j,k) + c1\_s(i,k+1)*s(i,j,k+1) 322 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 323 \textcolor{keywordflow}{ enddo} \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_a3b3416d0ee87cd1a0d54a910fb681e93}\label{namespacemom__diabatic__aux_a3b3416d0ee87cd1a0d54a910fb681e93}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!find\+\_\+uv\+\_\+at\+\_\+h@{find\+\_\+uv\+\_\+at\+\_\+h}} \index{find\+\_\+uv\+\_\+at\+\_\+h@{find\+\_\+uv\+\_\+at\+\_\+h}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{find\+\_\+uv\+\_\+at\+\_\+h()}{find\_uv\_at\_h()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::find\+\_\+uv\+\_\+at\+\_\+h (\begin{DoxyParamCaption}\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{u, }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)), intent(in)}]{v, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out)}]{u\+\_\+h, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out)}]{v\+\_\+h, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in), optional}]{ea, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in), optional}]{eb }\end{DoxyParamCaption})} This subroutine calculates u\+\_\+h and v\+\_\+h (velocities at thickness points), optionally using the entrainment amounts passed in as arguments. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em g} & The ocean\textquotesingle{}s grid structure\\ \hline \mbox{\tt in} & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in} & {\em u} & The zonal velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\ \hline \mbox{\tt in} & {\em v} & The meridional velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\ \hline \mbox{\tt in} & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt out} & {\em u\+\_\+h} & Zonal velocity interpolated to h points \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\ \hline \mbox{\tt out} & {\em v\+\_\+h} & Meridional velocity interpolated to h points \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\ \hline \mbox{\tt in} & {\em ea} & The amount of fluid entrained from the layer\\ \hline \mbox{\tt in} & {\em eb} & The amount of fluid entrained from the layer \\ \hline \end{DoxyParams} Definition at line 438 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 438 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< The ocean's grid structure} 439 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< The ocean's vertical grid structure} 440 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 441 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G))}, & 442 \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< The zonal velocity [L T-1 ~> m s-1]} 443 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G))}, & 444 \textcolor{keywordtype}{intent(in)} :: v\textcolor{comment}{ !< The meridional velocity [L T-1 ~> m s-1]} 445 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 446 \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< Layer thicknesses [H ~> m or kg m-2]} 447 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 448 \textcolor{keywordtype}{intent(out)} :: u\_h\textcolor{comment}{ !< Zonal velocity interpolated to h points [L T-1 ~> m s-1].} 449 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 450 \textcolor{keywordtype}{intent(out)} :: v\_h\textcolor{comment}{ !< Meridional velocity interpolated to h points [L T-1 ~> m s-1].} 451 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 452 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: ea\textcolor{comment}{ !< The amount of fluid entrained from the layer} 453 \textcolor{comment}{ !! above within this time step [H ~> m or kg m-2].} 454 \textcolor{comment}{ !! Omitting ea is the same as setting it to 0.} 455 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 456 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: eb\textcolor{comment}{ !< The amount of fluid entrained from the layer} 457 \textcolor{comment}{ !! below within this time step [H ~> m or kg m-2].} 458 \textcolor{comment}{ !! Omitting eb is the same as setting it to 0.} 459 460 \textcolor{comment}{! local variables} 461 \textcolor{keywordtype}{real} :: b\_denom\_1 \textcolor{comment}{! The first term in the denominator of b1 [H ~> m or kg m-2].} 462 \textcolor{keywordtype}{real} :: h\_neglect \textcolor{comment}{! A thickness that is so small it is usually lost} 463 \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].} 464 \textcolor{keywordtype}{real} :: b1(szi\_(g)), d1(szi\_(g)), c1(szi\_(g),szk\_(g)) 465 \textcolor{keywordtype}{real} :: a\_n(szi\_(g)), a\_s(szi\_(g)) \textcolor{comment}{! Fractional weights of the neighboring} 466 \textcolor{keywordtype}{real} :: a\_e(szi\_(g)), a\_w(szi\_(g)) \textcolor{comment}{! velocity points, ~1/2 in the open} 467 \textcolor{comment}{! ocean, nondimensional.} 468 \textcolor{keywordtype}{real} :: sum\_area, idenom 469 \textcolor{keywordtype}{logical} :: mix\_vertically 470 \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, nz 471 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke 472 \textcolor{keyword}{call }cpu\_clock\_begin(id\_clock\_uv\_at\_h) 473 h\_neglect = gv%H\_subroundoff 474 475 mix\_vertically = \textcolor{keyword}{present}(ea) 476 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(ea) .neqv. \textcolor{keyword}{present}(eb)) \textcolor{keyword}{call }mom\_error(fatal, & 477 \textcolor{stringliteral}{"find\_uv\_at\_h: Either both ea and eb or neither one must be present "}// & 478 \textcolor{stringliteral}{"in call to find\_uv\_at\_h."}) 479 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix\_vertically,h,h\_neglect, &} 480 \textcolor{comment}{!$OMP eb,u\_h,u,v\_h,v,nz,ea) &} 481 \textcolor{comment}{!$OMP private(sum\_area,Idenom,a\_w,a\_e,a\_s,a\_n,b\_denom\_1,b1,d1,c1)} 482 \textcolor{keywordflow}{do} j=js,je 483 \textcolor{keywordflow}{do} i=is,ie 484 sum\_area = g%areaCu(i-1,j) + g%areaCu(i,j) 485 \textcolor{keywordflow}{if} (sum\_area>0.0) \textcolor{keywordflow}{then} 486 idenom = sqrt(0.5*g%IareaT(i,j) / sum\_area) 487 a\_w(i) = g%areaCu(i-1,j) * idenom 488 a\_e(i) = g%areaCu(i,j) * idenom 489 \textcolor{keywordflow}{else} 490 a\_w(i) = 0.0 ; a\_e(i) = 0.0 491 \textcolor{keywordflow}{ endif} 492 493 sum\_area = g%areaCv(i,j-1) + g%areaCv(i,j) 494 \textcolor{keywordflow}{if} (sum\_area>0.0) \textcolor{keywordflow}{then} 495 idenom = sqrt(0.5*g%IareaT(i,j) / sum\_area) 496 a\_s(i) = g%areaCv(i,j-1) * idenom 497 a\_n(i) = g%areaCv(i,j) * idenom 498 \textcolor{keywordflow}{else} 499 a\_s(i) = 0.0 ; a\_n(i) = 0.0 500 \textcolor{keywordflow}{ endif} 501 \textcolor{keywordflow}{ enddo} 502 503 \textcolor{keywordflow}{if} (mix\_vertically) \textcolor{keywordflow}{then} 504 \textcolor{keywordflow}{do} i=is,ie 505 b\_denom\_1 = h(i,j,1) + h\_neglect 506 b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,1)) 507 d1(i) = b\_denom\_1 * b1(i) 508 u\_h(i,j,1) = (h(i,j,1)*b1(i)) * (a\_e(i)*u(i,j,1) + a\_w(i)*u(i-1,j,1)) 509 v\_h(i,j,1) = (h(i,j,1)*b1(i)) * (a\_n(i)*v(i,j,1) + a\_s(i)*v(i,j-1,1)) 510 \textcolor{keywordflow}{ enddo} 511 \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie 512 c1(i,k) = eb(i,j,k-1) * b1(i) 513 b\_denom\_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h\_neglect 514 b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,k)) 515 d1(i) = b\_denom\_1 * b1(i) 516 u\_h(i,j,k) = (h(i,j,k) * (a\_e(i)*u(i,j,k) + a\_w(i)*u(i-1,j,k)) + & 517 ea(i,j,k)*u\_h(i,j,k-1))*b1(i) 518 v\_h(i,j,k) = (h(i,j,k) * (a\_n(i)*v(i,j,k) + a\_s(i)*v(i,j-1,k)) + & 519 ea(i,j,k)*v\_h(i,j,k-1))*b1(i) 520 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 521 \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=is,ie 522 u\_h(i,j,k) = u\_h(i,j,k) + c1(i,k+1)*u\_h(i,j,k+1) 523 v\_h(i,j,k) = v\_h(i,j,k) + c1(i,k+1)*v\_h(i,j,k+1) 524 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 525 \textcolor{keywordflow}{else} 526 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie 527 u\_h(i,j,k) = a\_e(i)*u(i,j,k) + a\_w(i)*u(i-1,j,k) 528 v\_h(i,j,k) = a\_n(i)*v(i,j,k) + a\_s(i)*v(i,j-1,k) 529 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 530 \textcolor{keywordflow}{ endif} 531 \textcolor{keywordflow}{ enddo} 532 533 \textcolor{keyword}{call }cpu\_clock\_end(id\_clock\_uv\_at\_h) \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_adb99fd4e2092c17cd69143945733a6d9}\label{namespacemom__diabatic__aux_adb99fd4e2092c17cd69143945733a6d9}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!make\+\_\+frazil@{make\+\_\+frazil}} \index{make\+\_\+frazil@{make\+\_\+frazil}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{make\+\_\+frazil()}{make\_frazil()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::make\+\_\+frazil (\begin{DoxyParamCaption}\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h, }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(inout)}]{tv, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), intent(in)}]{CS, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g)), intent(in), optional}]{p\+\_\+surf, }\item[{integer, intent(in), optional}]{halo }\end{DoxyParamCaption})} Frazil formation keeps the temperature above the freezing point. This subroutine warms any water that is colder than the (currently surface) freezing point up to the freezing point and accumulates the required heat (in \mbox{[}Q R Z $\sim$$>$ J m-\/2\mbox{]}) in tvfrazil. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em g} & The ocean\textquotesingle{}s grid structure\\ \hline \mbox{\tt in} & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\ \hline \mbox{\tt in} & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in,out} & {\em tv} & Structure containing pointers to any available thermodynamic fields.\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in} & {\em cs} & The control structure returned by a previous call to diabatic\+\_\+aux\+\_\+init.\\ \hline \mbox{\tt in} & {\em p\+\_\+surf} & The pressure at the ocean surface \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]}.\\ \hline \mbox{\tt in} & {\em halo} & Halo width over which to calculate frazil \\ \hline \end{DoxyParams} Definition at line 104 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 104 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< The ocean's grid structure} 105 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< The ocean's vertical grid structure} 106 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, & 107 \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< Layer thicknesses [H ~> m or kg m-2]} 108 \textcolor{keywordtype}{type}(thermo\_var\_ptrs), \textcolor{keywordtype}{intent(inout)} :: tv\textcolor{comment}{ !< Structure containing pointers to any available} 109 \textcolor{comment}{ !! thermodynamic fields.} 110 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 111 \textcolor{keywordtype}{type}(diabatic\_aux\_cs), \textcolor{keywordtype}{intent(in)} :: cs\textcolor{comment}{ !< The control structure returned by a previous} 112 \textcolor{comment}{ !! call to diabatic\_aux\_init.} 113 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, & 114 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: p\_surf\textcolor{comment}{ !< The pressure at the ocean surface [R L2 T-2 ~> Pa].} 115 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: halo\textcolor{comment}{ !< Halo width over which to calculate frazil} 116 117 \textcolor{comment}{! Local variables} 118 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: & 119 fraz\_col, & \textcolor{comment}{! The accumulated heat requirement due to frazil [Q R Z ~> J m-2].} 120 t\_freeze, & \textcolor{comment}{! The freezing potential temperature at the current salinity [degC].} 121 ps \textcolor{comment}{! Surface pressure [R L2 T-2 ~> Pa]} 122 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZK\_(G))} :: & 123 pressure \textcolor{comment}{! The pressure at the middle of each layer [R L2 T-2 ~> Pa].} 124 \textcolor{keywordtype}{real} :: h\_to\_rl2\_t2 \textcolor{comment}{! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]} 125 \textcolor{keywordtype}{real} :: hc \textcolor{comment}{! A layer's heat capacity [Q R Z degC-1 ~> J m-2 degC-1].} 126 \textcolor{keywordtype}{logical} :: t\_fr\_set \textcolor{comment}{! True if the freezing point has been calculated for a} 127 \textcolor{comment}{! row of points.} 128 \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, nz 129 130 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke 131 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(halo)) \textcolor{keywordflow}{then} 132 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo 133 \textcolor{keywordflow}{ endif} 134 135 \textcolor{keyword}{call }cpu\_clock\_begin(id\_clock\_frazil) 136 137 \textcolor{keywordflow}{if} (.not.cs%pressure\_dependent\_frazil) \textcolor{keywordflow}{then} 138 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; pressure(i,k) = 0.0 ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 139 \textcolor{keywordflow}{else} 140 h\_to\_rl2\_t2 = gv%H\_to\_RZ * gv%g\_Earth 141 \textcolor{keywordflow}{ endif} 142 \textcolor{comment}{!$OMP parallel do default(shared) private(fraz\_col,T\_fr\_set,T\_freeze,hc,ps) &} 143 \textcolor{comment}{!$OMP firstprivate(pressure) ! pressure might be set above, so should be firstprivate} 144 \textcolor{keywordflow}{do} j=js,je 145 ps(:) = 0.0 146 \textcolor{keywordflow}{if} (\textcolor{keyword}{PRESENT}(p\_surf)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=is,ie 147 ps(i) = p\_surf(i,j) 148 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 149 150 \textcolor{keywordflow}{do} i=is,ie ; fraz\_col(i) = 0.0 ;\textcolor{keywordflow}{ enddo} 151 152 \textcolor{keywordflow}{if} (cs%pressure\_dependent\_frazil) \textcolor{keywordflow}{then} 153 \textcolor{keywordflow}{do} i=is,ie 154 pressure(i,1) = ps(i) + (0.5*h\_to\_rl2\_t2)*h(i,j,1) 155 \textcolor{keywordflow}{ enddo} 156 \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie 157 pressure(i,k) = pressure(i,k-1) + & 158 (0.5*h\_to\_rl2\_t2) * (h(i,j,k) + h(i,j,k-1)) 159 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 160 \textcolor{keywordflow}{ endif} 161 162 \textcolor{keywordflow}{if} (cs%reclaim\_frazil) \textcolor{keywordflow}{then} 163 t\_fr\_set = .false. 164 \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (tv%frazil(i,j) > 0.0) \textcolor{keywordflow}{then} 165 \textcolor{keywordflow}{if} (.not.t\_fr\_set) \textcolor{keywordflow}{then} 166 \textcolor{keyword}{call }calculate\_tfreeze(tv%S(i:,j,1), pressure(i:,1), t\_freeze(i:), & 167 1, ie-i+1, tv%eqn\_of\_state, pres\_scale=us%RL2\_T2\_to\_Pa) 168 t\_fr\_set = .true. 169 \textcolor{keywordflow}{ endif} 170 171 \textcolor{keywordflow}{if} (tv%T(i,j,1) > t\_freeze(i)) \textcolor{keywordflow}{then} 172 \textcolor{comment}{! If frazil had previously been formed, but the surface temperature is now} 173 \textcolor{comment}{! above freezing, cool the surface layer with the frazil heat deficit.} 174 hc = (tv%C\_p*gv%H\_to\_RZ) * h(i,j,1) 175 \textcolor{keywordflow}{if} (tv%frazil(i,j) - hc * (tv%T(i,j,1) - t\_freeze(i)) <= 0.0) \textcolor{keywordflow}{then} 176 tv%T(i,j,1) = tv%T(i,j,1) - tv%frazil(i,j) / hc 177 tv%frazil(i,j) = 0.0 178 \textcolor{keywordflow}{else} 179 tv%frazil(i,j) = tv%frazil(i,j) - hc * (tv%T(i,j,1) - t\_freeze(i)) 180 tv%T(i,j,1) = t\_freeze(i) 181 \textcolor{keywordflow}{ endif} 182 \textcolor{keywordflow}{ endif} 183 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} 184 \textcolor{keywordflow}{ endif} 185 186 \textcolor{keywordflow}{do} k=nz,1,-1 187 t\_fr\_set = .false. 188 \textcolor{keywordflow}{do} i=is,ie 189 \textcolor{keywordflow}{if} ((g%mask2dT(i,j) > 0.0) .and. & 190 ((tv%T(i,j,k) < 0.0) .or. (fraz\_col(i) > 0.0))) \textcolor{keywordflow}{then} 191 \textcolor{keywordflow}{if} (.not.t\_fr\_set) \textcolor{keywordflow}{then} 192 \textcolor{keyword}{call }calculate\_tfreeze(tv%S(i:,j,k), pressure(i:,k), t\_freeze(i:), & 193 1, ie-i+1, tv%eqn\_of\_state, pres\_scale=us%RL2\_T2\_to\_Pa) 194 t\_fr\_set = .true. 195 \textcolor{keywordflow}{ endif} 196 197 hc = (tv%C\_p*gv%H\_to\_RZ) * h(i,j,k) 198 \textcolor{keywordflow}{if} (h(i,j,k) <= 10.0*(gv%Angstrom\_H + gv%H\_subroundoff)) \textcolor{keywordflow}{then} 199 \textcolor{comment}{! Very thin layers should not be cooled by the frazil flux.} 200 \textcolor{keywordflow}{if} (tv%T(i,j,k) < t\_freeze(i)) \textcolor{keywordflow}{then} 201 fraz\_col(i) = fraz\_col(i) + hc * (t\_freeze(i) - tv%T(i,j,k)) 202 tv%T(i,j,k) = t\_freeze(i) 203 \textcolor{keywordflow}{ endif} 204 \textcolor{keywordflow}{elseif} ((fraz\_col(i) > 0.0) .or. (tv%T(i,j,k) < t\_freeze(i))) \textcolor{keywordflow}{then} 205 \textcolor{keywordflow}{if} (fraz\_col(i) + hc * (t\_freeze(i) - tv%T(i,j,k)) < 0.0) \textcolor{keywordflow}{then} 206 tv%T(i,j,k) = tv%T(i,j,k) - fraz\_col(i) / hc 207 fraz\_col(i) = 0.0 208 \textcolor{keywordflow}{else} 209 fraz\_col(i) = fraz\_col(i) + hc * (t\_freeze(i) - tv%T(i,j,k)) 210 tv%T(i,j,k) = t\_freeze(i) 211 \textcolor{keywordflow}{ endif} 212 \textcolor{keywordflow}{ endif} 213 \textcolor{keywordflow}{ endif} 214 \textcolor{keywordflow}{ enddo} 215 \textcolor{keywordflow}{ enddo} 216 \textcolor{keywordflow}{do} i=is,ie 217 tv%frazil(i,j) = tv%frazil(i,j) + fraz\_col(i) 218 \textcolor{keywordflow}{ enddo} 219 \textcolor{keywordflow}{ enddo} 220 \textcolor{keyword}{call }cpu\_clock\_end(id\_clock\_frazil) 221 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_adce369573ff355b820320d9c3b048693}\label{namespacemom__diabatic__aux_adce369573ff355b820320d9c3b048693}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!set\+\_\+pen\+\_\+shortwave@{set\+\_\+pen\+\_\+shortwave}} \index{set\+\_\+pen\+\_\+shortwave@{set\+\_\+pen\+\_\+shortwave}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{set\+\_\+pen\+\_\+shortwave()}{set\_pen\_shortwave()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::set\+\_\+pen\+\_\+shortwave (\begin{DoxyParamCaption}\item[{type(optics\+\_\+type), pointer}]{optics, }\item[{type(forcing), intent(inout)}]{fluxes, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{type(\hyperlink{structmom__diabatic__aux_1_1diabatic__aux__cs}{diabatic\+\_\+aux\+\_\+cs}), pointer}]{CS, }\item[{type(opacity\+\_\+cs), pointer}]{opacity\+\_\+\+C\+Sp, }\item[{type(tracer\+\_\+flow\+\_\+control\+\_\+cs), pointer}]{tracer\+\_\+flow\+\_\+\+C\+Sp }\end{DoxyParamCaption})} \begin{DoxyParams}[1]{Parameters} & {\em optics} & An optics structure that has will contain information about shortwave fluxes and absorption.\\ \hline \mbox{\tt in,out} & {\em fluxes} & points to forcing fields unused fields have N\+U\+LL ptrs\\ \hline \mbox{\tt in} & {\em g} & The ocean\textquotesingle{}s grid structure.\\ \hline \mbox{\tt in} & {\em gv} & The ocean\textquotesingle{}s vertical grid structure.\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline & {\em cs} & Control structure for diabatic\+\_\+aux\\ \hline & {\em opacity\+\_\+csp} & The control structure for the opacity module.\\ \hline & {\em tracer\+\_\+flow\+\_\+csp} & A pointer to the control structure organizing the tracer modules. \\ \hline \end{DoxyParams} Definition at line 538 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 538 \textcolor{keywordtype}{type}(optics\_type), \textcolor{keywordtype}{pointer} :: optics\textcolor{comment}{ !< An optics structure that has will contain} 539 \textcolor{comment}{ !! information about shortwave fluxes and absorption.} 540 \textcolor{keywordtype}{type}(forcing), \textcolor{keywordtype}{intent(inout)} :: fluxes\textcolor{comment}{ !< points to forcing fields} 541 \textcolor{comment}{ !! unused fields have NULL ptrs} 542 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< The ocean's grid structure.} 543 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< The ocean's vertical grid structure.} 544 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 545 \textcolor{keywordtype}{type}(diabatic\_aux\_cs), \textcolor{keywordtype}{pointer} :: cs\textcolor{comment}{ !< Control structure for diabatic\_aux} 546 \textcolor{keywordtype}{type}(opacity\_cs), \textcolor{keywordtype}{pointer} :: opacity\_csp\textcolor{comment}{ !< The control structure for the opacity module.} 547 \textcolor{keywordtype}{type}(tracer\_flow\_control\_cs), \textcolor{keywordtype}{pointer} :: tracer\_flow\_csp\textcolor{comment}{ !< A pointer to the control structure} 548 \textcolor{comment}{ !! organizing the tracer modules.} 549 550 \textcolor{comment}{! Local variables} 551 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))} :: chl\_2d\textcolor{comment}{ !< Vertically uniform chlorophyll-A concentractions [mg m-3]} 552 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(GV))} :: chl\_3d\textcolor{comment}{ !< The chlorophyll-A concentractions of each layer [mg m-3]} 553 \textcolor{keywordtype}{character(len=128)} :: mesg 554 \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je 555 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec 556 557 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(optics)) \textcolor{keywordflow}{return} 558 559 \textcolor{keywordflow}{if} (cs%var\_pen\_sw) \textcolor{keywordflow}{then} 560 \textcolor{keywordflow}{if} (cs%chl\_from\_file) \textcolor{keywordflow}{then} 561 \textcolor{comment}{! Only the 2-d surface chlorophyll can be read in from a file. The} 562 \textcolor{comment}{! same value is assumed for all layers.} 563 \textcolor{keyword}{call }time\_interp\_external(cs%sbc\_chl, cs%Time, chl\_2d) 564 \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie 565 \textcolor{keywordflow}{if} ((g%mask2dT(i,j) > 0.5) .and. (chl\_2d(i,j) < 0.0)) \textcolor{keywordflow}{then} 566 \textcolor{keyword}{write}(mesg,\textcolor{stringliteral}{'(" Time\_interp negative chl of ",(1pe12.4)," at i,j = ",&} 567 \textcolor{stringliteral}{}\textcolor{stringliteral}{ & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")'}) & 568 chl\_2d(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j) 569 \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_diabatic\_aux set\_pen\_shortwave: "}//trim(mesg)) 570 \textcolor{keywordflow}{ endif} 571 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 572 573 \textcolor{keywordflow}{if} (cs%id\_chl > 0) \textcolor{keyword}{call }post\_data(cs%id\_chl, chl\_2d, cs%diag) 574 575 \textcolor{keyword}{call }set\_opacity(optics, fluxes%sw, fluxes%sw\_vis\_dir, fluxes%sw\_vis\_dif, & 576 fluxes%sw\_nir\_dir, fluxes%sw\_nir\_dif, g, gv, us, opacity\_csp, chl\_2d=chl\_2d) 577 \textcolor{keywordflow}{else} 578 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(tracer\_flow\_csp)) \textcolor{keyword}{call }mom\_error(fatal, & 579 \textcolor{stringliteral}{"The tracer flow control structure must be associated when the model sets "}//& 580 \textcolor{stringliteral}{"the chlorophyll internally in set\_pen\_shortwave."}) 581 \textcolor{keyword}{call }get\_chl\_from\_model(chl\_3d, g, tracer\_flow\_csp) 582 583 \textcolor{keywordflow}{if} (cs%id\_chl > 0) \textcolor{keyword}{call }post\_data(cs%id\_chl, chl\_3d(:,:,1), cs%diag) 584 585 \textcolor{keyword}{call }set\_opacity(optics, fluxes%sw, fluxes%sw\_vis\_dir, fluxes%sw\_vis\_dif, & 586 fluxes%sw\_nir\_dir, fluxes%sw\_nir\_dif, g, gv, us, opacity\_csp, chl\_3d=chl\_3d) 587 \textcolor{keywordflow}{ endif} 588 \textcolor{keywordflow}{else} 589 \textcolor{keyword}{call }set\_opacity(optics, fluxes%sw, fluxes%sw\_vis\_dir, fluxes%sw\_vis\_dif, & 590 fluxes%sw\_nir\_dir, fluxes%sw\_nir\_dif, g, gv, us, opacity\_csp) 591 \textcolor{keywordflow}{ endif} 592 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__diabatic__aux_acde4c76c9a489b82cba0bac5ec1b1198}\label{namespacemom__diabatic__aux_acde4c76c9a489b82cba0bac5ec1b1198}} \index{mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}!tridiagts@{tridiagts}} \index{tridiagts@{tridiagts}!mom\+\_\+diabatic\+\_\+aux@{mom\+\_\+diabatic\+\_\+aux}} \subsubsection{\texorpdfstring{tridiagts()}{tridiagts()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+diabatic\+\_\+aux\+::tridiagts (\begin{DoxyParamCaption}\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{integer, intent(in)}]{is, }\item[{integer, intent(in)}]{ie, }\item[{integer, intent(in)}]{js, }\item[{integer, intent(in)}]{je, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{hold, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{ea, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{eb, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{T, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{S }\end{DoxyParamCaption})} This is a simple tri-\/diagonal solver for T and S. \char`\"{}\+Simple\char`\"{} means it only uses arrays hold, ea and eb. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em g} & The ocean\textquotesingle{}s grid structure\\ \hline \mbox{\tt in} & {\em gv} & The ocean\textquotesingle{}s vertical grid structure\\ \hline \mbox{\tt in} & {\em is} & The start i-\/index to work on.\\ \hline \mbox{\tt in} & {\em ie} & The end i-\/index to work on.\\ \hline \mbox{\tt in} & {\em js} & The start j-\/index to work on.\\ \hline \mbox{\tt in} & {\em je} & The end j-\/index to work on.\\ \hline \mbox{\tt in} & {\em hold} & The layer thicknesses before entrainment, \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}.\\ \hline \mbox{\tt in} & {\em ea} & The amount of fluid entrained from the layer above within this time step \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em eb} & The amount of fluid entrained from the layer below within this time step \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in,out} & {\em t} & Layer potential temperatures \mbox{[}degC\mbox{]}.\\ \hline \mbox{\tt in,out} & {\em s} & Layer salinities \mbox{[}ppt\mbox{]}. \\ \hline \end{DoxyParams} Definition at line 389 of file M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90. \begin{DoxyCode} 389 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{ !< The ocean's grid structure} 390 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< The ocean's vertical grid structure} 391 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: is\textcolor{comment}{ !< The start i-index to work on.} 392 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: ie\textcolor{comment}{ !< The end i-index to work on.} 393 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: js\textcolor{comment}{ !< The start j-index to work on.} 394 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: je\textcolor{comment}{ !< The end j-index to work on.} 395 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(in)} :: hold\textcolor{comment}{ !< The layer thicknesses before entrainment,} 396 \textcolor{comment}{ !! [H ~> m or kg m-2].} 397 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(in)} :: ea\textcolor{comment}{ !< The amount of fluid entrained from the layer} 398 \textcolor{comment}{ !! above within this time step [H ~> m or kg m-2]} 399 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(in)} :: eb\textcolor{comment}{ !< The amount of fluid entrained from the layer} 400 \textcolor{comment}{ !! below within this time step [H ~> m or kg m-2]} 401 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(inout)} :: t\textcolor{comment}{ !< Layer potential temperatures [degC].} 402 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(inout)} :: s\textcolor{comment}{ !< Layer salinities [ppt].} 403 404 \textcolor{comment}{! Local variables} 405 \textcolor{keywordtype}{real} :: b1(szib\_(g)), d1(szib\_(g)) \textcolor{comment}{! b1, c1, and d1 are variables used by the} 406 \textcolor{keywordtype}{real} :: c1(szib\_(g),szk\_(g)) \textcolor{comment}{! tridiagonal solver.} 407 \textcolor{keywordtype}{real} :: h\_tr, b\_denom\_1 408 \textcolor{keywordtype}{integer} :: i, j, k 409 410 \textcolor{comment}{!$OMP parallel do default(shared) private(h\_tr,b1,d1,c1,b\_denom\_1)} 411 \textcolor{keywordflow}{do} j=js,je 412 \textcolor{keywordflow}{do} i=is,ie 413 h\_tr = hold(i,j,1) + gv%H\_subroundoff 414 b1(i) = 1.0 / (h\_tr + eb(i,j,1)) 415 d1(i) = h\_tr * b1(i) 416 t(i,j,1) = (b1(i)*h\_tr)*t(i,j,1) 417 s(i,j,1) = (b1(i)*h\_tr)*s(i,j,1) 418 \textcolor{keywordflow}{ enddo} 419 \textcolor{keywordflow}{do} k=2,g%ke ; \textcolor{keywordflow}{do} i=is,ie 420 c1(i,k) = eb(i,j,k-1) * b1(i) 421 h\_tr = hold(i,j,k) + gv%H\_subroundoff 422 b\_denom\_1 = h\_tr + d1(i)*ea(i,j,k) 423 b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,k)) 424 d1(i) = b\_denom\_1 * b1(i) 425 t(i,j,k) = b1(i) * (h\_tr*t(i,j,k) + ea(i,j,k)*t(i,j,k-1)) 426 s(i,j,k) = b1(i) * (h\_tr*s(i,j,k) + ea(i,j,k)*s(i,j,k-1)) 427 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 428 \textcolor{keywordflow}{do} k=g%ke-1,1,-1 ; \textcolor{keywordflow}{do} i=is,ie 429 t(i,j,k) = t(i,j,k) + c1(i,k+1)*t(i,j,k+1) 430 s(i,j,k) = s(i,j,k) + c1(i,k+1)*s(i,j,k+1) 431 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 432 \textcolor{keywordflow}{ enddo} \end{DoxyCode}