\hypertarget{namespacemom__tracer__diabatic}{}\section{mom\+\_\+tracer\+\_\+diabatic Module Reference} \label{namespacemom__tracer__diabatic}\index{mom\+\_\+tracer\+\_\+diabatic@{mom\+\_\+tracer\+\_\+diabatic}} \subsection{Detailed Description} This module contains routines that implement physical fluxes of tracers (e.\+g. due to surface fluxes or mixing). These are intended to be called from call\+\_\+tracer\+\_\+column\+\_\+fns in the \mbox{\hyperlink{namespaceMOM__tracer__flow__control}{M\+O\+M\+\_\+tracer\+\_\+flow\+\_\+control}} module. \subsection*{Functions/\+Subroutines} \begin{DoxyCompactItemize} \item subroutine, public \mbox{\hyperlink{namespacemom__tracer__diabatic_ac5d57973547cc4ed3a89808d3910943e}{tracer\+\_\+vertdiff}} (h\+\_\+old, ea, eb, dt, tr, G, GV, sfc\+\_\+flux, btm\+\_\+flux, btm\+\_\+reservoir, sink\+\_\+rate, convert\+\_\+flux\+\_\+in) \begin{DoxyCompactList}\small\item\em This subroutine solves a tridiagonal equation for the final tracer concentrations after the dual-\/entrainments, and possibly sinking or surface and bottom sources, are applied. The sinking is implemented with an fully implicit upwind advection scheme. Alternate time units can be used for the timestep, surface and bottom fluxes and sink\+\_\+rate provided they are all consistent. \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespacemom__tracer__diabatic_ad4d3d4de0f2b84c15bccc5eb2f767df3}{applytracerboundaryfluxesinout}} (G, GV, Tr, dt, fluxes, h, evap\+\_\+\+C\+F\+L\+\_\+limit, minimum\+\_\+forcing\+\_\+depth, in\+\_\+flux\+\_\+optional, out\+\_\+flux\+\_\+optional, update\+\_\+h\+\_\+opt) \begin{DoxyCompactList}\small\item\em This routine is modeled after apply\+Boundary\+Fluxes\+In\+Out in \mbox{\hyperlink{MOM__diabatic__aux_8F90_source}{M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90}} N\+O\+TE\+: Please note that in this routine sfc\+\_\+flux gets set to zero to ensure that the surface flux of the tracer does not get applied again during a subsequent call to tracer\+\_\+vertdif. \end{DoxyCompactList}\end{DoxyCompactItemize} \subsection{Function/\+Subroutine Documentation} \mbox{\Hypertarget{namespacemom__tracer__diabatic_ad4d3d4de0f2b84c15bccc5eb2f767df3}\label{namespacemom__tracer__diabatic_ad4d3d4de0f2b84c15bccc5eb2f767df3}} \index{mom\+\_\+tracer\+\_\+diabatic@{mom\+\_\+tracer\+\_\+diabatic}!applytracerboundaryfluxesinout@{applytracerboundaryfluxesinout}} \index{applytracerboundaryfluxesinout@{applytracerboundaryfluxesinout}!mom\+\_\+tracer\+\_\+diabatic@{mom\+\_\+tracer\+\_\+diabatic}} \subsubsection{\texorpdfstring{applytracerboundaryfluxesinout()}{applytracerboundaryfluxesinout()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+tracer\+\_\+diabatic\+::applytracerboundaryfluxesinout (\begin{DoxyParamCaption}\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{Tr, }\item[{real, intent(in)}]{dt, }\item[{type(forcing), intent(in)}]{fluxes, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{h, }\item[{real, intent(in)}]{evap\+\_\+\+C\+F\+L\+\_\+limit, }\item[{real, intent(in)}]{minimum\+\_\+forcing\+\_\+depth, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g)), intent(in), optional}]{in\+\_\+flux\+\_\+optional, }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g)), intent(in), optional}]{out\+\_\+flux\+\_\+optional, }\item[{logical, intent(in), optional}]{update\+\_\+h\+\_\+opt }\end{DoxyParamCaption})} This routine is modeled after apply\+Boundary\+Fluxes\+In\+Out in \mbox{\hyperlink{MOM__diabatic__aux_8F90_source}{M\+O\+M\+\_\+diabatic\+\_\+aux.\+F90}} N\+O\+TE\+: Please note that in this routine sfc\+\_\+flux gets set to zero to ensure that the surface flux of the tracer does not get applied again during a subsequent call to tracer\+\_\+vertdif. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em g} & Grid structure\\ \hline \mbox{\tt in} & {\em gv} & ocean vertical grid structure\\ \hline \mbox{\tt in,out} & {\em tr} & Tracer concentration on T-\/cell\\ \hline \mbox{\tt in} & {\em dt} & Time-\/step over which forcing is applied \mbox{[}T $\sim$$>$ s\mbox{]}\\ \hline \mbox{\tt in} & {\em fluxes} & Surface fluxes container\\ \hline \mbox{\tt in,out} & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em evap\+\_\+cfl\+\_\+limit} & Limit on the fraction of the water that can be fluxed out of the top layer in a timestep \mbox{[}nondim\mbox{]}\\ \hline \mbox{\tt in} & {\em minimum\+\_\+forcing\+\_\+depth} & The smallest depth over which fluxes can be applied \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em in\+\_\+flux\+\_\+optional} & The total time-\/integrated amount of tracer that enters with freshwater\\ \hline \mbox{\tt in} & {\em out\+\_\+flux\+\_\+optional} & The total time-\/integrated amount of tracer that leaves with freshwater\\ \hline \mbox{\tt in} & {\em update\+\_\+h\+\_\+opt} & Optional flag to determine whether h should be updated \\ \hline \end{DoxyParams} Definition at line 230 of file M\+O\+M\+\_\+tracer\+\_\+diabatic.\+F90. \begin{DoxyCode} 230 231 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in )} :: G\textcolor{comment}{ !< Grid structure} 232 \textcolor{keywordtype}{type}(verticalGrid\_type), \textcolor{keywordtype}{intent(in )} :: GV\textcolor{comment}{ !< ocean vertical grid structure} 233 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(inout)} :: Tr\textcolor{comment}{ !< Tracer concentration on T-cell} 234 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in )} :: dt\textcolor{comment}{ !< Time-step over which forcing is applied [T ~> s]} 235 \textcolor{keywordtype}{type}(forcing), \textcolor{keywordtype}{intent(in )} :: fluxes\textcolor{comment}{ !< Surface fluxes container} 236 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(inout)} :: h\textcolor{comment}{ !< Layer thickness [H ~> m or kg m-2]} 237 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in )} :: evap\_CFL\_limit\textcolor{comment}{ !< Limit on the fraction of the} 238 \textcolor{comment}{ !! water that can be fluxed out of the top} 239 \textcolor{comment}{ !! layer in a timestep [nondim]} 240 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in )} :: minimum\_forcing\_depth\textcolor{comment}{ !< The smallest depth over} 241 \textcolor{comment}{ !! which fluxes can be applied [H ~> m or kg m-2]} 242 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in )} :: in\_flux\_optional\textcolor{comment}{ !< The total time-integrated} 243 \textcolor{comment}{ !! amount of tracer that enters with freshwater} 244 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: out\_flux\_optional\textcolor{comment}{ !< The total time-integrated} 245 \textcolor{comment}{ !! amount of tracer that leaves with freshwater} 246 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: update\_h\_opt\textcolor{comment}{ !< Optional flag to determine whether} 247 \textcolor{comment}{ !! h should be updated} 248 249 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: maxGroundings = 5 250 \textcolor{keywordtype}{integer} :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings) 251 \textcolor{keywordtype}{real} :: H\_limit\_fluxes, IforcingDepthScale 252 \textcolor{keywordtype}{real} :: dThickness, dTracer 253 \textcolor{keywordtype}{real} :: fractionOfForcing, hOld, Ithickness 254 \textcolor{keywordtype}{real} :: RivermixConst \textcolor{comment}{! A constant used in implementing river mixing [Pa s].} 255 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: & 256 netMassInOut, & \textcolor{comment}{! surface water fluxes [H ~> m or kg m-2] over time step} 257 netMassIn, & \textcolor{comment}{! mass entering ocean surface [H ~> m or kg m-2] over a time step} 258 netMassOut \textcolor{comment}{! mass leaving ocean surface [H ~> m or kg m-2] over a time step} 259 260 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZK\_(G))} :: h2d, Tr2d 261 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))} :: in\_flux \textcolor{comment}{! The total time-integrated amount of tracer!} 262 \textcolor{comment}{! that enters with freshwater} 263 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))} :: out\_flux \textcolor{comment}{! The total time-integrated amount of tracer!} 264 \textcolor{comment}{! that leaves with freshwater} 265 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: in\_flux\_1d, out\_flux\_1d 266 \textcolor{keywordtype}{real} :: hGrounding(maxGroundings) 267 \textcolor{keywordtype}{real} :: Tr\_in 268 \textcolor{keywordtype}{logical} :: update\_h 269 \textcolor{keywordtype}{integer} :: i, j, is, ie, js, je, k, nz, n, nsw 270 \textcolor{keywordtype}{character(len=45)} :: mesg 271 272 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke 273 274 \textcolor{comment}{! If no freshwater fluxes, nothing needs to be done in this routine} 275 \textcolor{keywordflow}{if} ( (.not. \textcolor{keyword}{associated}(fluxes%netMassIn)) .or. (.not. \textcolor{keyword}{associated}(fluxes%netMassOut)) ) \textcolor{keywordflow}{return} 276 277 in\_flux(:,:) = 0.0 ; out\_flux(:,:) = 0.0 278 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(in\_flux\_optional)) \textcolor{keywordflow}{then} 279 \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie 280 in\_flux(i,j) = in\_flux\_optional(i,j) 281 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 282 \textcolor{keywordflow}{ endif} 283 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(out\_flux\_optional)) \textcolor{keywordflow}{then} 284 \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie 285 out\_flux(i,j) = out\_flux\_optional(i,j) 286 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 287 \textcolor{keywordflow}{ endif} 288 289 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(update\_h\_opt)) \textcolor{keywordflow}{then} 290 update\_h = update\_h\_opt 291 \textcolor{keywordflow}{else} 292 update\_h = .true. 293 \textcolor{keywordflow}{ endif} 294 295 numberofgroundings = 0 296 297 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,Tr,G,GV,fluxes,dt, &} 298 \textcolor{comment}{!$OMP IforcingDepthScale,minimum\_forcing\_depth, &} 299 \textcolor{comment}{!$OMP numberOfGroundings,iGround,jGround,update\_h, &} 300 \textcolor{comment}{!$OMP in\_flux,out\_flux,hGrounding,evap\_CFL\_limit) &} 301 \textcolor{comment}{!$OMP private(h2d,Tr2d,netMassInOut,netMassOut, &} 302 \textcolor{comment}{!$OMP in\_flux\_1d,out\_flux\_1d,fractionOfForcing, &} 303 \textcolor{comment}{!$OMP dThickness,dTracer,hOld,Ithickness, &} 304 \textcolor{comment}{!$OMP netMassIn, Tr\_in)} 305 306 \textcolor{comment}{! Work in vertical slices for efficiency} 307 \textcolor{keywordflow}{do} j=js,je 308 309 \textcolor{comment}{! Copy state into 2D-slice arrays} 310 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie 311 h2d(i,k) = h(i,j,k) 312 tr2d(i,k) = tr(i,j,k) 313 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 314 315 \textcolor{keywordflow}{do} i = is,ie 316 in\_flux\_1d(i) = in\_flux(i,j) 317 out\_flux\_1d(i) = out\_flux(i,j) 318 \textcolor{keywordflow}{ enddo} 319 \textcolor{comment}{! The surface forcing is contained in the fluxes type.} 320 \textcolor{comment}{! We aggregate the thermodynamic forcing for a time step into the following:} 321 \textcolor{comment}{! These should have been set and stored during a call to applyBoundaryFluxesInOut} 322 \textcolor{comment}{! netMassIn = net mass entering at ocean surface over a timestep} 323 \textcolor{comment}{! netMassOut = net mass leaving ocean surface [H ~> m or kg m-2] over a time step.} 324 \textcolor{comment}{! netMassOut < 0 means mass leaves ocean.} 325 326 \textcolor{comment}{! Note here that the aggregateFW flag has already been taken care of in the call to} 327 \textcolor{comment}{! applyBoundaryFluxesInOut} 328 \textcolor{keywordflow}{do} i=is,ie 329 netmassout(i) = fluxes%netMassOut(i,j) 330 netmassin(i) = fluxes%netMassIn(i,j) 331 \textcolor{keywordflow}{ enddo} 332 333 \textcolor{comment}{! Apply the surface boundary fluxes in three steps:} 334 \textcolor{comment}{! A/ update concentration from mass entering the ocean} 335 \textcolor{comment}{! B/ update concentration from mass leaving ocean.} 336 \textcolor{keywordflow}{do} i=is,ie 337 \textcolor{keywordflow}{if} (g%mask2dT(i,j)>0.) \textcolor{keywordflow}{then} 338 339 \textcolor{comment}{! A/ Update tracer due to incoming mass flux.} 340 \textcolor{keywordflow}{do} k=1,1 341 342 \textcolor{comment}{! Change in state due to forcing} 343 dthickness = netmassin(i) \textcolor{comment}{! Since we are adding mass, we can use all of it} 344 dtracer = 0. 345 346 \textcolor{comment}{! Update the forcing by the part to be consumed within the present k-layer.} 347 \textcolor{comment}{! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.} 348 netmassin(i) = netmassin(i) - dthickness 349 dtracer = dtracer + in\_flux\_1d(i) 350 in\_flux\_1d(i) = 0.0 351 352 \textcolor{comment}{! Update state} 353 hold = h2d(i,k) \textcolor{comment}{! Keep original thickness in hand} 354 h2d(i,k) = h2d(i,k) + dthickness \textcolor{comment}{! New thickness} 355 \textcolor{keywordflow}{if} (h2d(i,k) > 0.0) \textcolor{keywordflow}{then} 356 ithickness = 1.0/h2d(i,k) \textcolor{comment}{! Inverse new thickness} 357 \textcolor{comment}{! The "if"s below avoid changing T/S by roundoff unnecessarily} 358 \textcolor{keywordflow}{if} (dthickness /= 0. .or. dtracer /= 0.) tr2d(i,k) = (hold*tr2d(i,k)+ dtracer)*ithickness 359 \textcolor{keywordflow}{ endif} 360 361 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! k=1,1} 362 363 \textcolor{comment}{! B/ Update tracer from mass leaving ocean} 364 \textcolor{keywordflow}{do} k=1,nz 365 366 \textcolor{comment}{! Place forcing into this layer if this layer has nontrivial thickness.} 367 \textcolor{comment}{! For layers thin relative to 1/IforcingDepthScale, then distribute} 368 \textcolor{comment}{! forcing into deeper layers.} 369 iforcingdepthscale = 1. / max(gv%H\_subroundoff, minimum\_forcing\_depth - netmassout(i) ) 370 \textcolor{comment}{! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.} 371 fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale) 372 373 \textcolor{comment}{! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we} 374 \textcolor{comment}{! limit the forcing applied to this cell, leaving the remaining forcing to} 375 \textcolor{comment}{! be distributed downwards.} 376 \textcolor{keywordflow}{if} (-fractionofforcing*netmassout(i) > evap\_cfl\_limit*h2d(i,k)) \textcolor{keywordflow}{then} 377 fractionofforcing = -evap\_cfl\_limit*h2d(i,k)/netmassout(i) 378 \textcolor{keywordflow}{ endif} 379 380 \textcolor{comment}{! Change in state due to forcing} 381 dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) ) 382 \textcolor{comment}{! Note this is slightly different to how salt is currently treated} 383 dtracer = fractionofforcing*out\_flux\_1d(i) 384 385 \textcolor{comment}{! Update the forcing by the part to be consumed within the present k-layer.} 386 \textcolor{comment}{! If fractionOfForcing = 1, then new netMassOut vanishes.} 387 netmassout(i) = netmassout(i) - dthickness 388 out\_flux\_1d(i) = out\_flux\_1d(i) - dtracer 389 390 \textcolor{comment}{! Update state by the appropriate increment.} 391 hold = h2d(i,k) \textcolor{comment}{! Keep original thickness in hand} 392 h2d(i,k) = h2d(i,k) + dthickness \textcolor{comment}{! New thickness} 393 \textcolor{keywordflow}{if} (h2d(i,k) > 0.) \textcolor{keywordflow}{then} 394 ithickness = 1.0/h2d(i,k) \textcolor{comment}{! Inverse of new thickness} 395 tr2d(i,k) = (hold*tr2d(i,k) + dtracer)*ithickness 396 \textcolor{keywordflow}{ endif} 397 398 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! k} 399 400 \textcolor{keywordflow}{ endif} 401 402 \textcolor{comment}{! If anything remains after the k-loop, then we have grounded out, which is a problem.} 403 \textcolor{keywordflow}{if} (abs(in\_flux\_1d(i))+abs(out\_flux\_1d(i)) /= 0.0) \textcolor{keywordflow}{then} 404 \textcolor{comment}{!$OMP critical} 405 numberofgroundings = numberofgroundings +1 406 \textcolor{keywordflow}{if} (numberofgroundings<=maxgroundings) \textcolor{keywordflow}{then} 407 iground(numberofgroundings) = i \textcolor{comment}{! Record i,j location of event for} 408 jground(numberofgroundings) = j \textcolor{comment}{! warning message} 409 hgrounding(numberofgroundings) = abs(in\_flux\_1d(i))+abs(out\_flux\_1d(i)) 410 \textcolor{keywordflow}{ endif} 411 \textcolor{comment}{!$OMP end critical} 412 \textcolor{keywordflow}{ endif} 413 414 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! i} 415 416 \textcolor{comment}{! Step C/ copy updated tracer concentration from the 2d slice now back into model state.} 417 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie 418 tr(i,j,k) = tr2d(i,k) 419 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 420 421 \textcolor{keywordflow}{if} (update\_h) \textcolor{keywordflow}{then} 422 \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie 423 h(i,j,k) = h2d(i,k) 424 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 425 \textcolor{keywordflow}{ endif} 426 427 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! j-loop finish} 428 429 \textcolor{keywordflow}{if} (numberofgroundings>0) \textcolor{keywordflow}{then} 430 \textcolor{keywordflow}{do} i = 1, min(numberofgroundings, maxgroundings) 431 \textcolor{keyword}{write}(mesg(1:45),\textcolor{stringliteral}{'(3es15.3)'}) g%geoLonT( iground(i), jground(i) ), & 432 g%geoLatT( iground(i), jground(i)) , hgrounding(i) 433 \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"MOM\_tracer\_diabatic.F90, applyTracerBoundaryFluxesInOut(): "}//& 434 \textcolor{stringliteral}{"Tracer created. x,y,dh= "}//trim(mesg), all\_print=.true.) 435 \textcolor{keywordflow}{ enddo} 436 437 \textcolor{keywordflow}{if} (numberofgroundings - maxgroundings > 0) \textcolor{keywordflow}{then} 438 \textcolor{keyword}{write}(mesg, \textcolor{stringliteral}{'(i4)'}) numberofgroundings - maxgroundings 439 \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"MOM\_tracer\_vertical.F90, applyTracerBoundaryFluxesInOut(): "}//& 440 trim(mesg) // \textcolor{stringliteral}{" groundings remaining"}) 441 \textcolor{keywordflow}{ endif} 442 \textcolor{keywordflow}{ endif} 443 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__tracer__diabatic_ac5d57973547cc4ed3a89808d3910943e}\label{namespacemom__tracer__diabatic_ac5d57973547cc4ed3a89808d3910943e}} \index{mom\+\_\+tracer\+\_\+diabatic@{mom\+\_\+tracer\+\_\+diabatic}!tracer\+\_\+vertdiff@{tracer\+\_\+vertdiff}} \index{tracer\+\_\+vertdiff@{tracer\+\_\+vertdiff}!mom\+\_\+tracer\+\_\+diabatic@{mom\+\_\+tracer\+\_\+diabatic}} \subsubsection{\texorpdfstring{tracer\+\_\+vertdiff()}{tracer\_vertdiff()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+tracer\+\_\+diabatic\+::tracer\+\_\+vertdiff (\begin{DoxyParamCaption}\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, gv \%ke), intent(in)}]{h\+\_\+old, }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, gv \%ke), intent(in)}]{ea, }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, gv \%ke), intent(in)}]{eb, }\item[{real, intent(in)}]{dt, }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, gv \%ke), intent(inout)}]{tr, }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed), intent(in), optional}]{sfc\+\_\+flux, }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed), intent(in), optional}]{btm\+\_\+flux, }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed), intent(inout), optional}]{btm\+\_\+reservoir, }\item[{real, intent(in), optional}]{sink\+\_\+rate, }\item[{logical, intent(in), optional}]{convert\+\_\+flux\+\_\+in }\end{DoxyParamCaption})} This subroutine solves a tridiagonal equation for the final tracer concentrations after the dual-\/entrainments, and possibly sinking or surface and bottom sources, are applied. The sinking is implemented with an fully implicit upwind advection scheme. Alternate time units can be used for the timestep, surface and bottom fluxes and sink\+\_\+rate provided they are all consistent. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em g} & ocean grid structure\\ \hline \mbox{\tt in} & {\em gv} & ocean vertical grid structure\\ \hline \mbox{\tt in} & {\em h\+\_\+old} & layer thickness before entrainment \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em ea} & amount of fluid entrained from the layer above \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em eb} & amount of fluid entrained from the layer below \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\ \hline \mbox{\tt in,out} & {\em tr} & tracer concentration in concentration units \mbox{[}CU\mbox{]}\\ \hline \mbox{\tt in} & {\em dt} & amount of time covered by this call \mbox{[}T $\sim$$>$ s\mbox{]}\\ \hline \mbox{\tt in} & {\em sfc\+\_\+flux} & surface flux of the tracer in units of \mbox{[}CU kg m-\/2 T-\/1 $\sim$$>$ CU kg m-\/2 s-\/1\mbox{]} or \mbox{[}CU H $\sim$$>$ CU m or CU kg m-\/2\mbox{]} if convert\+\_\+flux\+\_\+in is .false.\\ \hline \mbox{\tt in} & {\em btm\+\_\+flux} & The (negative upward) bottom flux of the tracer in \mbox{[}CU kg m-\/2 T-\/1 $\sim$$>$ CU kg m-\/2 s-\/1\mbox{]} or \mbox{[}CU H $\sim$$>$ CU m or CU kg m-\/2\mbox{]} if\\ \hline \mbox{\tt in,out} & {\em btm\+\_\+reservoir} & amount of tracer in a bottom reservoir \mbox{[}CU kg m-\/2\mbox{]}; formerly \mbox{[}CU m\mbox{]}\\ \hline \mbox{\tt in} & {\em sink\+\_\+rate} & rate at which the tracer sinks \mbox{[}m T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\ \hline \mbox{\tt in} & {\em convert\+\_\+flux\+\_\+in} & True if the specified sfc\+\_\+flux needs to be integrated in time \\ \hline \end{DoxyParams} Definition at line 27 of file M\+O\+M\+\_\+tracer\+\_\+diabatic.\+F90. \begin{DoxyCode} 27 \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: G\textcolor{comment}{ !< ocean grid structure} 28 \textcolor{keywordtype}{type}(verticalGrid\_type), \textcolor{keywordtype}{intent(in)} :: GV\textcolor{comment}{ !< ocean vertical grid structure} 29 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(GV))}, \textcolor{keywordtype}{intent(in)} :: h\_old\textcolor{comment}{ !< layer thickness before entrainment} 30 \textcolor{comment}{ !! [H ~> m or kg m-2]} 31 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(GV))}, \textcolor{keywordtype}{intent(in)} :: ea\textcolor{comment}{ !< amount of fluid entrained from the layer} 32 \textcolor{comment}{ !! above [H ~> m or kg m-2]} 33 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(GV))}, \textcolor{keywordtype}{intent(in)} :: eb\textcolor{comment}{ !< amount of fluid entrained from the layer} 34 \textcolor{comment}{ !! below [H ~> m or kg m-2]} 35 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(GV))}, \textcolor{keywordtype}{intent(inout)} :: tr\textcolor{comment}{ !< tracer concentration in concentration units [CU]} 36 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: dt\textcolor{comment}{ !< amount of time covered by this call [T ~> s]} 37 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, \textcolor{keywordtype}{optional},\textcolor{keywordtype}{intent(in)} :: sfc\_flux\textcolor{comment}{ !< surface flux of the tracer in units of} 38 \textcolor{comment}{ !! [CU kg m-2 T-1 ~> CU kg m-2 s-1] or} 39 \textcolor{comment}{ !! [CU H ~> CU m or CU kg m-2] if} 40 \textcolor{comment}{ !! convert\_flux\_in is .false.} 41 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, \textcolor{keywordtype}{optional},\textcolor{keywordtype}{intent(in)} :: btm\_flux\textcolor{comment}{ !< The (negative upward) bottom flux of the} 42 \textcolor{comment}{ !! tracer in [CU kg m-2 T-1 ~> CU kg m-2 s-1] or} 43 \textcolor{comment}{ !! [CU H ~> CU m or CU kg m-2] if} 44 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, \textcolor{keywordtype}{optional},\textcolor{keywordtype}{intent(inout)} :: btm\_reservoir\textcolor{comment}{ !< amount of tracer in a bottom reservoir} 45 \textcolor{comment}{ !! [CU kg m-2]; formerly [CU m]} 46 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional},\textcolor{keywordtype}{intent(in)} :: sink\_rate\textcolor{comment}{ !< rate at which the tracer sinks} 47 \textcolor{comment}{ !! [m T-1 ~> m s-1]} 48 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional},\textcolor{keywordtype}{intent(in)} :: convert\_flux\_in\textcolor{comment}{ !< True if the specified sfc\_flux needs} 49 \textcolor{comment}{ !! to be integrated in time} 50 51 \textcolor{comment}{! local variables} 52 \textcolor{keywordtype}{real} :: sink\_dist\textcolor{comment}{ !< The distance the tracer sinks in a time step [H ~> m or kg m-2].} 53 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))} :: & 54 sfc\_src, & !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2]. 55 btm\_src\textcolor{comment}{ !< The time-integrated bottom source of the tracer [CU H ~> CU m or CU kg m-2].} 56 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: & 57 b1, & !< b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. 58 d1 \textcolor{comment}{!! d1=1-c1 is used by the tridiagonal solver, nondimensional.} 59 \textcolor{keywordtype}{real} :: c1(SZI\_(G),SZK\_(GV))\textcolor{comment}{ !< c1 is used by the tridiagonal solver [nondim].} 60 \textcolor{keywordtype}{real} :: h\_minus\_dsink(SZI\_(G),SZK\_(GV))\textcolor{comment}{ !< The layer thickness minus the} 61 \textcolor{comment}{ !! difference in sinking rates across the layer [H ~> m or kg m-2].} 62 \textcolor{comment}{ !! By construction, 0 <= h\_minus\_dsink < h\_work.} 63 \textcolor{keywordtype}{real} :: sink(SZI\_(G),SZK\_(GV)+1)\textcolor{comment}{ !< The tracer's sinking distances at the} 64 \textcolor{comment}{ !! interfaces, limited to prevent characteristics from} 65 \textcolor{comment}{ !! crossing within a single timestep [H ~> m or kg m-2].} 66 \textcolor{keywordtype}{real} :: b\_denom\_1\textcolor{comment}{ !< The first term in the denominator of b1 [H ~> m or kg m-2].} 67 \textcolor{keywordtype}{real} :: h\_tr\textcolor{comment}{ !< h\_tr is h at tracer points with a h\_neglect added to} 68 \textcolor{comment}{ !! ensure positive definiteness [H ~> m or kg m-2].} 69 \textcolor{keywordtype}{real} :: h\_neglect\textcolor{comment}{ !< A thickness that is so small it is usually lost} 70 \textcolor{comment}{ !! in roundoff and can be neglected [H ~> m or kg m-2].} 71 \textcolor{keywordtype}{logical} :: convert\_flux = .true. 72 73 74 \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, nz 75 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke 76 77 \textcolor{keywordflow}{if} (nz == 1) \textcolor{keywordflow}{then} 78 \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"MOM\_tracer\_diabatic.F90, tracer\_vertdiff called "}//& 79 \textcolor{stringliteral}{"with only one vertical level"}) 80 \textcolor{keywordflow}{return} 81 \textcolor{keywordflow}{ endif} 82 83 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(convert\_flux\_in)) convert\_flux = convert\_flux\_in 84 h\_neglect = gv%H\_subroundoff 85 sink\_dist = 0.0 86 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(sink\_rate)) sink\_dist = (dt*sink\_rate) * gv%m\_to\_H 87 \textcolor{comment}{!$OMP parallel default(none) shared(is,ie,js,je,sfc\_src,btm\_src,sfc\_flux,dt,G,GV,btm\_flux, &} 88 \textcolor{comment}{!$OMP sink\_rate,btm\_reservoir,nz,sink\_dist,ea, &} 89 \textcolor{comment}{!$OMP h\_old,convert\_flux,h\_neglect,eb,tr) &} 90 \textcolor{comment}{!$OMP private(sink,h\_minus\_dsink,b\_denom\_1,b1,d1,h\_tr,c1)} 91 \textcolor{comment}{!$OMP do} 92 \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is,ie ; sfc\_src(i,j) = 0.0 ; btm\_src(i,j) = 0.0 ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 93 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(sfc\_flux)) \textcolor{keywordflow}{then} 94 \textcolor{keywordflow}{if} (convert\_flux) \textcolor{keywordflow}{then} 95 \textcolor{comment}{!$OMP do} 96 \textcolor{keywordflow}{do} j = js, je; \textcolor{keywordflow}{do} i = is,ie 97 sfc\_src(i,j) = (sfc\_flux(i,j)*dt) * gv%kg\_m2\_to\_H 98 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 99 \textcolor{keywordflow}{else} 100 \textcolor{comment}{!$OMP do} 101 \textcolor{keywordflow}{do} j = js, je; \textcolor{keywordflow}{do} i = is,ie 102 sfc\_src(i,j) = sfc\_flux(i,j) 103 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 104 \textcolor{keywordflow}{ endif} 105 \textcolor{keywordflow}{ endif} 106 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(btm\_flux)) \textcolor{keywordflow}{then} 107 \textcolor{keywordflow}{if} (convert\_flux) \textcolor{keywordflow}{then} 108 \textcolor{comment}{!$OMP do} 109 \textcolor{keywordflow}{do} j = js, je; \textcolor{keywordflow}{do} i = is,ie 110 btm\_src(i,j) = (btm\_flux(i,j)*dt) * gv%kg\_m2\_to\_H 111 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 112 \textcolor{keywordflow}{else} 113 \textcolor{comment}{!$OMP do} 114 \textcolor{keywordflow}{do} j = js, je; \textcolor{keywordflow}{do} i = is,ie 115 btm\_src(i,j) = btm\_flux(i,j) 116 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 117 \textcolor{keywordflow}{ endif} 118 \textcolor{keywordflow}{ endif} 119 120 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(sink\_rate)) \textcolor{keywordflow}{then} 121 \textcolor{comment}{!$OMP do} 122 \textcolor{keywordflow}{do} j=js,je 123 \textcolor{comment}{! Find the sinking rates at all interfaces, limiting them if necesary} 124 \textcolor{comment}{! so that the characteristics do not cross within a timestep.} 125 \textcolor{comment}{! If a non-constant sinking rate were used, that would be incorprated} 126 \textcolor{comment}{! here.} 127 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(btm\_reservoir)) \textcolor{keywordflow}{then} 128 \textcolor{keywordflow}{do} i=is,ie ; sink(i,nz+1) = sink\_dist ;\textcolor{keywordflow}{ enddo} 129 \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie 130 sink(i,k) = sink\_dist ; h\_minus\_dsink(i,k) = h\_old(i,j,k) 131 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 132 \textcolor{keywordflow}{else} 133 \textcolor{keywordflow}{do} i=is,ie ; sink(i,nz+1) = 0.0 ;\textcolor{keywordflow}{ enddo} 134 \textcolor{comment}{! Find the limited sinking distance at the interfaces.} 135 \textcolor{keywordflow}{do} k=nz,2,-1 ; \textcolor{keywordflow}{do} i=is,ie 136 \textcolor{keywordflow}{if} (sink(i,k+1) >= sink\_dist) \textcolor{keywordflow}{then} 137 sink(i,k) = sink\_dist 138 h\_minus\_dsink(i,k) = h\_old(i,j,k) + (sink(i,k+1) - sink(i,k)) 139 \textcolor{keywordflow}{elseif} (sink(i,k+1) + h\_old(i,j,k) < sink\_dist) \textcolor{keywordflow}{then} 140 sink(i,k) = sink(i,k+1) + h\_old(i,j,k) 141 h\_minus\_dsink(i,k) = 0.0 142 \textcolor{keywordflow}{else} 143 sink(i,k) = sink\_dist 144 h\_minus\_dsink(i,k) = (h\_old(i,j,k) + sink(i,k+1)) - sink(i,k) 145 \textcolor{keywordflow}{ endif} 146 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 147 \textcolor{keywordflow}{ endif} 148 \textcolor{keywordflow}{do} i=is,ie 149 sink(i,1) = 0.0 ; h\_minus\_dsink(i,1) = (h\_old(i,j,1) + sink(i,2)) 150 \textcolor{keywordflow}{ enddo} 151 152 \textcolor{comment}{! Now solve the tridiagonal equation for the tracer concentrations.} 153 \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then} 154 b\_denom\_1 = h\_minus\_dsink(i,1) + ea(i,j,1) + h\_neglect 155 b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,1)) 156 d1(i) = b\_denom\_1 * b1(i) 157 h\_tr = h\_old(i,j,1) + h\_neglect 158 tr(i,j,1) = (b1(i)*h\_tr)*tr(i,j,1) + sfc\_src(i,j) 159 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} 160 \textcolor{keywordflow}{do} k=2,nz-1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then} 161 c1(i,k) = eb(i,j,k-1) * b1(i) 162 b\_denom\_1 = h\_minus\_dsink(i,k) + d1(i) * (ea(i,j,k) + sink(i,k)) + & 163 h\_neglect 164 b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,k)) 165 d1(i) = b\_denom\_1 * b1(i) 166 h\_tr = h\_old(i,j,k) + h\_neglect 167 tr(i,j,k) = b1(i) * (h\_tr * tr(i,j,k) + & 168 (ea(i,j,k) + sink(i,k)) * tr(i,j,k-1)) 169 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 170 \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then} 171 c1(i,nz) = eb(i,j,nz-1) * b1(i) 172 b\_denom\_1 = h\_minus\_dsink(i,nz) + d1(i) * (ea(i,j,nz) + sink(i,nz)) + & 173 h\_neglect 174 b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,nz)) 175 h\_tr = h\_old(i,j,nz) + h\_neglect 176 tr(i,j,nz) = b1(i) * ((h\_tr * tr(i,j,nz) + btm\_src(i,j)) + & 177 (ea(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1)) 178 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} 179 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(btm\_reservoir)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j)>0.5) \textcolor{keywordflow}{then} 180 btm\_reservoir(i,j) = btm\_reservoir(i,j) + & 181 (sink(i,nz+1)*tr(i,j,nz)) * gv%H\_to\_kg\_m2 182 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 183 184 \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then} 185 tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) 186 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 187 \textcolor{keywordflow}{ enddo} 188 \textcolor{keywordflow}{else} 189 \textcolor{comment}{!$OMP do} 190 \textcolor{keywordflow}{do} j=js,je 191 \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j) > -0.5) \textcolor{keywordflow}{then} 192 h\_tr = h\_old(i,j,1) + h\_neglect 193 b\_denom\_1 = h\_tr + ea(i,j,1) 194 b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,1)) 195 d1(i) = h\_tr * b1(i) 196 tr(i,j,1) = (b1(i)*h\_tr)*tr(i,j,1) + sfc\_src(i,j) 197 \textcolor{keywordflow}{ endif} 198 \textcolor{keywordflow}{ enddo} 199 \textcolor{keywordflow}{do} k=2,nz-1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j) > -0.5) \textcolor{keywordflow}{then} 200 c1(i,k) = eb(i,j,k-1) * b1(i) 201 h\_tr = h\_old(i,j,k) + h\_neglect 202 b\_denom\_1 = h\_tr + d1(i) * ea(i,j,k) 203 b1(i) = 1.0 / (b\_denom\_1 + eb(i,j,k)) 204 d1(i) = b\_denom\_1 * b1(i) 205 tr(i,j,k) = b1(i) * (h\_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1)) 206 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 207 \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j) > -0.5) \textcolor{keywordflow}{then} 208 c1(i,nz) = eb(i,j,nz-1) * b1(i) 209 h\_tr = h\_old(i,j,nz) + h\_neglect 210 b\_denom\_1 = h\_tr + d1(i)*ea(i,j,nz) 211 b1(i) = 1.0 / ( b\_denom\_1 + eb(i,j,nz)) 212 tr(i,j,nz) = b1(i) * (( h\_tr * tr(i,j,nz) + btm\_src(i,j)) + & 213 ea(i,j,nz) * tr(i,j,nz-1)) 214 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} 215 \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j) > -0.5) \textcolor{keywordflow}{then} 216 tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) 217 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 218 \textcolor{keywordflow}{ enddo} 219 \textcolor{keywordflow}{ endif} 220 221 \textcolor{comment}{!$OMP end parallel} 222 \end{DoxyCode}