\hypertarget{namespacep3m__functions}{}\section{p3m\+\_\+functions Module Reference} \label{namespacep3m__functions}\index{p3m\+\_\+functions@{p3m\+\_\+functions}} \subsection{Detailed Description} Cubic interpolation functions. Date of creation\+: 2008.\+06.\+09 L. White. This module contains p3m interpolation routines. p3m interpolation is performed by estimating the edge values and slopes and constructing a cubic polynomial. We then make sure that the edge values are bounded and continuous and we then modify the slopes to get a monotonic cubic curve. \subsection*{Functions/\+Subroutines} \begin{DoxyCompactItemize} \item subroutine, public \mbox{\hyperlink{namespacep3m__functions_a685d3ef292536dae810b2059ccaa5819}{p3m\+\_\+interpolation}} (N, h, u, edge\+\_\+values, ppoly\+\_\+S, ppoly\+\_\+coef, h\+\_\+neglect, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Set up a piecewise cubic interpolation from cell averages and estimated edge slopes and values. \end{DoxyCompactList}\item subroutine \mbox{\hyperlink{namespacep3m__functions_a44db59cb5df7f139e05b745746342fcf}{p3m\+\_\+limiter}} (N, h, u, edge\+\_\+values, ppoly\+\_\+S, ppoly\+\_\+coef, h\+\_\+neglect, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Adust a piecewise cubic reconstruction with a limiter that adjusts the edge values and slopes. \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespacep3m__functions_ab70eb9e69fc6586e1d6a371d2eeb44d1}{p3m\+\_\+boundary\+\_\+extrapolation}} (N, h, u, edge\+\_\+values, ppoly\+\_\+S, ppoly\+\_\+coef, h\+\_\+neglect, h\+\_\+neglect\+\_\+edge) \begin{DoxyCompactList}\small\item\em Calculate the edge values and slopes at boundary cells as part of building a piecewise cubic sub-\/grid scale profiles. \end{DoxyCompactList}\item subroutine \mbox{\hyperlink{namespacep3m__functions_a82c244f63ffd696c15d58a029e8e24f7}{build\+\_\+cubic\+\_\+interpolant}} (h, k, edge\+\_\+values, ppoly\+\_\+S, ppoly\+\_\+coef) \begin{DoxyCompactList}\small\item\em Build cubic interpolant in cell k. \end{DoxyCompactList}\item logical function \mbox{\hyperlink{namespacep3m__functions_a3c461688b7c3ae5b2b4a7cf3311b8a69}{is\+\_\+cubic\+\_\+monotonic}} (ppoly\+\_\+coef, k) \begin{DoxyCompactList}\small\item\em Check whether the cubic reconstruction in cell k is monotonic. \end{DoxyCompactList}\item subroutine \mbox{\hyperlink{namespacep3m__functions_adb96651fe725f11e90dec2b8509989b0}{monotonize\+\_\+cubic}} (h, u0\+\_\+l, u0\+\_\+r, sigma\+\_\+l, sigma\+\_\+r, slope, u1\+\_\+l, u1\+\_\+r) \begin{DoxyCompactList}\small\item\em Monotonize a cubic curve by modifying the edge slopes. \end{DoxyCompactList}\end{DoxyCompactItemize} \subsection*{Variables} \begin{DoxyCompactItemize} \item \mbox{\Hypertarget{namespacep3m__functions_aedfe0a0f65453ca61b90d4c601747c90}\label{namespacep3m__functions_aedfe0a0f65453ca61b90d4c601747c90}} real, parameter \mbox{\hyperlink{namespacep3m__functions_aedfe0a0f65453ca61b90d4c601747c90}{hneglect\+\_\+dflt}} = 1.E-\/30 \begin{DoxyCompactList}\small\item\em Default value of a negligible cell thickness. \end{DoxyCompactList}\item \mbox{\Hypertarget{namespacep3m__functions_abdb3a0dd6c75e81bbcb5f5c982c31e7b}\label{namespacep3m__functions_abdb3a0dd6c75e81bbcb5f5c982c31e7b}} real, parameter \mbox{\hyperlink{namespacep3m__functions_abdb3a0dd6c75e81bbcb5f5c982c31e7b}{hneglect\+\_\+edge\+\_\+dflt}} = 1.E-\/10 \begin{DoxyCompactList}\small\item\em Default value of a negligible edge thickness. \end{DoxyCompactList}\end{DoxyCompactItemize} \subsection{Function/\+Subroutine Documentation} \mbox{\Hypertarget{namespacep3m__functions_a82c244f63ffd696c15d58a029e8e24f7}\label{namespacep3m__functions_a82c244f63ffd696c15d58a029e8e24f7}} \index{p3m\+\_\+functions@{p3m\+\_\+functions}!build\+\_\+cubic\+\_\+interpolant@{build\+\_\+cubic\+\_\+interpolant}} \index{build\+\_\+cubic\+\_\+interpolant@{build\+\_\+cubic\+\_\+interpolant}!p3m\+\_\+functions@{p3m\+\_\+functions}} \subsubsection{\texorpdfstring{build\+\_\+cubic\+\_\+interpolant()}{build\_cubic\_interpolant()}} {\footnotesize\ttfamily subroutine p3m\+\_\+functions\+::build\+\_\+cubic\+\_\+interpolant (\begin{DoxyParamCaption}\item[{real, dimension(\+:), intent(in)}]{h, }\item[{integer, intent(in)}]{k, }\item[{real, dimension(\+:,\+:), intent(in)}]{edge\+\_\+values, }\item[{real, dimension(\+:,\+:), intent(in)}]{ppoly\+\_\+S, }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+coef }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}} Build cubic interpolant in cell k. Given edge values and edge slopes, compute coefficients of cubic in cell k. N\+O\+TE\+: edge values and slopes M\+U\+ST have been properly calculated prior to calling this routine. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em h} & cell widths (size N) \mbox{[}H\mbox{]}\\ \hline \mbox{\tt in} & {\em k} & The index of the cell to work on\\ \hline \mbox{\tt in} & {\em edge\+\_\+values} & Edge value of polynomial in arbitrary units \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in} & {\em ppoly\+\_\+s} & Edge slope of polynomial \mbox{[}A H-\/1\mbox{]}\\ \hline \mbox{\tt in,out} & {\em ppoly\+\_\+coef} & Coefficients of polynomial \mbox{[}A\mbox{]} \\ \hline \end{DoxyParams} Definition at line 344 of file P3\+M\+\_\+functions.\+F90. \begin{DoxyCode} 344 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths (size N) [H]} 345 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: k\textcolor{comment}{ !< The index of the cell to work on} 346 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(in)} :: edge\_values\textcolor{comment}{ !< Edge value of polynomial in arbitrary units [A]} 347 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(in)} :: ppoly\_S\textcolor{comment}{ !< Edge slope of polynomial [A H-1]} 348 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_coef\textcolor{comment}{ !< Coefficients of polynomial [A]} 349 350 \textcolor{comment}{! Local variables} 351 \textcolor{keywordtype}{real} :: u0\_l, u0\_r \textcolor{comment}{! edge values [A]} 352 \textcolor{keywordtype}{real} :: u1\_l, u1\_r \textcolor{comment}{! edge slopes times the cell width [A]} 353 \textcolor{keywordtype}{real} :: h\_c \textcolor{comment}{! cell width [H]} 354 \textcolor{keywordtype}{real} :: a0, a1, a2, a3 \textcolor{comment}{! cubic coefficients [A]} 355 356 h\_c = h(k) 357 358 u0\_l = edge\_values(k,1) 359 u0\_r = edge\_values(k,2) 360 361 u1\_l = ppoly\_s(k,1) * h\_c 362 u1\_r = ppoly\_s(k,2) * h\_c 363 364 a0 = u0\_l 365 a1 = u1\_l 366 a2 = 3.0 * ( u0\_r - u0\_l ) - u1\_r - 2.0 * u1\_l 367 a3 = u1\_r + u1\_l + 2.0 * ( u0\_l - u0\_r ) 368 369 ppoly\_coef(k,1) = a0 370 ppoly\_coef(k,2) = a1 371 ppoly\_coef(k,3) = a2 372 ppoly\_coef(k,4) = a3 373 \end{DoxyCode} \mbox{\Hypertarget{namespacep3m__functions_a3c461688b7c3ae5b2b4a7cf3311b8a69}\label{namespacep3m__functions_a3c461688b7c3ae5b2b4a7cf3311b8a69}} \index{p3m\+\_\+functions@{p3m\+\_\+functions}!is\+\_\+cubic\+\_\+monotonic@{is\+\_\+cubic\+\_\+monotonic}} \index{is\+\_\+cubic\+\_\+monotonic@{is\+\_\+cubic\+\_\+monotonic}!p3m\+\_\+functions@{p3m\+\_\+functions}} \subsubsection{\texorpdfstring{is\+\_\+cubic\+\_\+monotonic()}{is\_cubic\_monotonic()}} {\footnotesize\ttfamily logical function p3m\+\_\+functions\+::is\+\_\+cubic\+\_\+monotonic (\begin{DoxyParamCaption}\item[{real, dimension(\+:,\+:), intent(in)}]{ppoly\+\_\+coef, }\item[{integer, intent(in)}]{k }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}} Check whether the cubic reconstruction in cell k is monotonic. This function checks whether the cubic curve in cell k is monotonic. If so, returns 1. Otherwise, returns 0. The cubic is monotonic if the first derivative is single-\/signed in (0,1). Hence, we check whether the roots (if any) lie inside this interval. If there is no root or if both roots lie outside this interval, the cubic is monotonic. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em ppoly\+\_\+coef} & Coefficients of cubic polynomial in arbitary units \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in} & {\em k} & The index of the cell to work on \\ \hline \end{DoxyParams} Definition at line 386 of file P3\+M\+\_\+functions.\+F90. \begin{DoxyCode} 386 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(in)} :: ppoly\_coef\textcolor{comment}{ !< Coefficients of cubic polynomial in arbitary units [A]} 387 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: k\textcolor{comment}{ !< The index of the cell to work on} 388 \textcolor{comment}{! Local variables} 389 \textcolor{keywordtype}{real} :: a, b, c \textcolor{comment}{! Coefficients of the first derivative of the cubic [A]} 390 391 a = ppoly\_coef(k,2) 392 b = 2.0 * ppoly\_coef(k,3) 393 c = 3.0 * ppoly\_coef(k,4) 394 395 \textcolor{comment}{! Look for real roots of the quadratic derivative equation, c*x**2 + b*x + a = 0, in (0, 1)} 396 \textcolor{keywordflow}{if} (b*b - 4.0*a*c <= 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! The cubic is monotonic everywhere.} 397 is\_cubic\_monotonic = .true. 398 \textcolor{keywordflow}{elseif} (a * (a + (b + c)) < 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! The derivative changes sign between the endpoints of (0, 1)} 399 is\_cubic\_monotonic = .false. 400 \textcolor{keywordflow}{elseif} (b * (b + 2.0*c) < 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! The second derivative changes sign inside of (0, 1)} 401 is\_cubic\_monotonic = .false. 402 \textcolor{keywordflow}{else} 403 is\_cubic\_monotonic = .true. 404 \textcolor{keywordflow}{ endif} 405 \end{DoxyCode} \mbox{\Hypertarget{namespacep3m__functions_adb96651fe725f11e90dec2b8509989b0}\label{namespacep3m__functions_adb96651fe725f11e90dec2b8509989b0}} \index{p3m\+\_\+functions@{p3m\+\_\+functions}!monotonize\+\_\+cubic@{monotonize\+\_\+cubic}} \index{monotonize\+\_\+cubic@{monotonize\+\_\+cubic}!p3m\+\_\+functions@{p3m\+\_\+functions}} \subsubsection{\texorpdfstring{monotonize\+\_\+cubic()}{monotonize\_cubic()}} {\footnotesize\ttfamily subroutine p3m\+\_\+functions\+::monotonize\+\_\+cubic (\begin{DoxyParamCaption}\item[{real, intent(in)}]{h, }\item[{real, intent(in)}]{u0\+\_\+l, }\item[{real, intent(in)}]{u0\+\_\+r, }\item[{real, intent(in)}]{sigma\+\_\+l, }\item[{real, intent(in)}]{sigma\+\_\+r, }\item[{real, intent(in)}]{slope, }\item[{real, intent(inout)}]{u1\+\_\+l, }\item[{real, intent(inout)}]{u1\+\_\+r }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}} Monotonize a cubic curve by modifying the edge slopes. This routine takes care of monotonizing a cubic on \mbox{[}0,1\mbox{]} by modifying the edge slopes. The edge values are N\+OT modified. The cubic is entirely determined by the four degrees of freedom u0\+\_\+l, u0\+\_\+r, u1\+\_\+l and u1\+\_\+r. u1\+\_\+l and u1\+\_\+r are the edge slopes expressed in the G\+L\+O\+B\+AL coordinate system. The monotonization occurs as follows. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em h} & cell width \mbox{[}H\mbox{]}\\ \hline \mbox{\tt in} & {\em u0\+\_\+l} & left edge value in arbitrary units \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in} & {\em u0\+\_\+r} & right edge value \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in} & {\em sigma\+\_\+l} & left 2nd-\/order slopes \mbox{[}A H-\/1\mbox{]}\\ \hline \mbox{\tt in} & {\em sigma\+\_\+r} & right 2nd-\/order slopes \mbox{[}A H-\/1\mbox{]}\\ \hline \mbox{\tt in} & {\em slope} & limited P\+LM slope \mbox{[}A H-\/1\mbox{]}\\ \hline \mbox{\tt in,out} & {\em u1\+\_\+l} & left edge slopes \mbox{[}A H-\/1\mbox{]}\\ \hline \mbox{\tt in,out} & {\em u1\+\_\+r} & right edge slopes \mbox{[}A H-\/1\mbox{]} \\ \hline \end{DoxyParams} Definition at line 436 of file P3\+M\+\_\+functions.\+F90. \begin{DoxyCode} 436 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell width [H]} 437 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: u0\_l\textcolor{comment}{ !< left edge value in arbitrary units [A]} 438 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: u0\_r\textcolor{comment}{ !< right edge value [A]} 439 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: sigma\_l\textcolor{comment}{ !< left 2nd-order slopes [A H-1]} 440 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: sigma\_r\textcolor{comment}{ !< right 2nd-order slopes [A H-1]} 441 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: slope\textcolor{comment}{ !< limited PLM slope [A H-1]} 442 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(inout)} :: u1\_l\textcolor{comment}{ !< left edge slopes [A H-1]} 443 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(inout)} :: u1\_r\textcolor{comment}{ !< right edge slopes [A H-1]} 444 \textcolor{comment}{! Local variables} 445 \textcolor{keywordtype}{logical} :: found\_ip 446 \textcolor{keywordtype}{logical} :: inflexion\_l \textcolor{comment}{! bool telling if inflex. pt must be on left} 447 \textcolor{keywordtype}{logical} :: inflexion\_r \textcolor{comment}{! bool telling if inflex. pt must be on right} 448 \textcolor{keywordtype}{real} :: a1, a2, a3 \textcolor{comment}{! Temporary slopes times the cell width [A]} 449 \textcolor{keywordtype}{real} :: u1\_l\_tmp \textcolor{comment}{! trial left edge slope [A H-1]} 450 \textcolor{keywordtype}{real} :: u1\_r\_tmp \textcolor{comment}{! trial right edge slope [A H-1]} 451 \textcolor{keywordtype}{real} :: xi\_ip \textcolor{comment}{! location of inflexion point in cell coordinates (0,1) [nondim]} 452 \textcolor{keywordtype}{real} :: slope\_ip \textcolor{comment}{! slope at inflexion point times cell width [A]} 453 454 found\_ip = .false. 455 inflexion\_l = .false. 456 inflexion\_r = .false. 457 458 \textcolor{comment}{! If the edge slopes are inconsistent w.r.t. the limited PLM slope,} 459 \textcolor{comment}{! set them to zero} 460 \textcolor{keywordflow}{if} ( u1\_l*slope <= 0.0 ) \textcolor{keywordflow}{then} 461 u1\_l = 0.0 462 \textcolor{keywordflow}{ endif} 463 464 \textcolor{keywordflow}{if} ( u1\_r*slope <= 0.0 ) \textcolor{keywordflow}{then} 465 u1\_r = 0.0 466 \textcolor{keywordflow}{ endif} 467 468 \textcolor{comment}{! Compute the location of the inflexion point, which is the root} 469 \textcolor{comment}{! of the second derivative} 470 a1 = h * u1\_l 471 a2 = 3.0 * ( u0\_r - u0\_l ) - h*(u1\_r + 2.0*u1\_l) 472 a3 = h*(u1\_r + u1\_l) + 2.0*(u0\_l - u0\_r) 473 474 \textcolor{comment}{! There is a possible root (and inflexion point) only if a3 is nonzero.} 475 \textcolor{comment}{! When a3 is zero, the second derivative of the cubic is constant (the} 476 \textcolor{comment}{! cubic degenerates into a parabola) and no inflexion point exists.} 477 \textcolor{keywordflow}{if} ( a3 /= 0.0 ) \textcolor{keywordflow}{then} 478 \textcolor{comment}{! Location of inflexion point} 479 xi\_ip = - a2 / (3.0 * a3) 480 481 \textcolor{comment}{! If the inflexion point lies in [0,1], change boolean value} 482 \textcolor{keywordflow}{if} ( (xi\_ip >= 0.0) .AND. (xi\_ip <= 1.0) ) \textcolor{keywordflow}{then} 483 found\_ip = .true. 484 \textcolor{keywordflow}{ endif} 485 \textcolor{keywordflow}{ endif} 486 487 \textcolor{comment}{! When there is an inflexion point within [0,1], check the slope} 488 \textcolor{comment}{! to see if it is consistent with the limited PLM slope. If not,} 489 \textcolor{comment}{! decide on which side we want to collapse the inflexion point.} 490 \textcolor{comment}{! If the inflexion point lies on one of the edges, the cubic is} 491 \textcolor{comment}{! guaranteed to be monotonic} 492 \textcolor{keywordflow}{if} ( found\_ip ) \textcolor{keywordflow}{then} 493 slope\_ip = a1 + 2.0*a2*xi\_ip + 3.0*a3*xi\_ip*xi\_ip 494 495 \textcolor{comment}{! Check whether slope is consistent} 496 \textcolor{keywordflow}{if} ( slope\_ip*slope < 0.0 ) \textcolor{keywordflow}{then} 497 \textcolor{keywordflow}{if} ( abs(sigma\_l) < abs(sigma\_r) ) \textcolor{keywordflow}{then} 498 inflexion\_l = .true. 499 \textcolor{keywordflow}{else} 500 inflexion\_r = .true. 501 \textcolor{keywordflow}{ endif} 502 \textcolor{keywordflow}{ endif} 503 \textcolor{keywordflow}{ endif} \textcolor{comment}{! found\_ip} 504 505 \textcolor{comment}{! At this point, if the cubic is not monotonic, we know where the} 506 \textcolor{comment}{! inflexion point should lie. When the cubic is monotonic, both} 507 \textcolor{comment}{! 'inflexion\_l' and 'inflexion\_r' are false and nothing is to be done.} 508 509 \textcolor{comment}{! Move inflexion point on the left} 510 \textcolor{keywordflow}{if} ( inflexion\_l ) \textcolor{keywordflow}{then} 511 512 u1\_l\_tmp = 1.5*(u0\_r-u0\_l)/h - 0.5*u1\_r 513 u1\_r\_tmp = 3.0*(u0\_r-u0\_l)/h - 2.0*u1\_l 514 515 \textcolor{keywordflow}{if} ( (u1\_l\_tmp*slope < 0.0) .AND. (u1\_r\_tmp*slope < 0.0) ) \textcolor{keywordflow}{then} 516 517 u1\_l = 0.0 518 u1\_r = 3.0 * (u0\_r - u0\_l) / h 519 520 \textcolor{keywordflow}{elseif} (u1\_l\_tmp*slope < 0.0) \textcolor{keywordflow}{then} 521 522 u1\_r = u1\_r\_tmp 523 u1\_l = 1.5*(u0\_r - u0\_l)/h - 0.5*u1\_r 524 525 \textcolor{keywordflow}{elseif} (u1\_r\_tmp*slope < 0.0) \textcolor{keywordflow}{then} 526 527 u1\_l = u1\_l\_tmp 528 u1\_r = 3.0*(u0\_r - u0\_l)/h - 2.0*u1\_l 529 530 \textcolor{keywordflow}{else} 531 532 u1\_l = u1\_l\_tmp 533 u1\_r = u1\_r\_tmp 534 535 \textcolor{keywordflow}{ endif} 536 537 \textcolor{keywordflow}{ endif} \textcolor{comment}{! end treating case with inflexion point on the left} 538 539 \textcolor{comment}{! Move inflexion point on the right} 540 \textcolor{keywordflow}{if} ( inflexion\_r ) \textcolor{keywordflow}{then} 541 542 u1\_l\_tmp = 3.0*(u0\_r-u0\_l)/h - 2.0*u1\_r 543 u1\_r\_tmp = 1.5*(u0\_r-u0\_l)/h - 0.5*u1\_l 544 545 \textcolor{keywordflow}{if} ( (u1\_l\_tmp*slope < 0.0) .AND. (u1\_r\_tmp*slope < 0.0) ) \textcolor{keywordflow}{then} 546 547 u1\_l = 3.0 * (u0\_r - u0\_l) / h 548 u1\_r = 0.0 549 550 \textcolor{keywordflow}{elseif} (u1\_l\_tmp*slope < 0.0) \textcolor{keywordflow}{then} 551 552 u1\_r = u1\_r\_tmp 553 u1\_l = 3.0*(u0\_r - u0\_l)/h - 2.0*u1\_r 554 555 \textcolor{keywordflow}{elseif} (u1\_r\_tmp*slope < 0.0) \textcolor{keywordflow}{then} 556 557 u1\_l = u1\_l\_tmp 558 u1\_r = 1.5*(u0\_r - u0\_l)/h - 0.5*u1\_l 559 560 \textcolor{keywordflow}{else} 561 562 u1\_l = u1\_l\_tmp 563 u1\_r = u1\_r\_tmp 564 565 \textcolor{keywordflow}{ endif} 566 567 \textcolor{keywordflow}{ endif} \textcolor{comment}{! end treating case with inflexion point on the right} 568 569 \textcolor{comment}{! Zero out negligibly small slopes.} 570 \textcolor{keywordflow}{if} ( abs(u1\_l*h) < epsilon(u0\_l) * (abs(u0\_l) + abs(u0\_r)) ) u1\_l = 0.0 571 \textcolor{keywordflow}{if} ( abs(u1\_r*h) < epsilon(u0\_l) * (abs(u0\_l) + abs(u0\_r)) ) u1\_r = 0.0 572 \end{DoxyCode} \mbox{\Hypertarget{namespacep3m__functions_ab70eb9e69fc6586e1d6a371d2eeb44d1}\label{namespacep3m__functions_ab70eb9e69fc6586e1d6a371d2eeb44d1}} \index{p3m\+\_\+functions@{p3m\+\_\+functions}!p3m\+\_\+boundary\+\_\+extrapolation@{p3m\+\_\+boundary\+\_\+extrapolation}} \index{p3m\+\_\+boundary\+\_\+extrapolation@{p3m\+\_\+boundary\+\_\+extrapolation}!p3m\+\_\+functions@{p3m\+\_\+functions}} \subsubsection{\texorpdfstring{p3m\+\_\+boundary\+\_\+extrapolation()}{p3m\_boundary\_extrapolation()}} {\footnotesize\ttfamily subroutine, public p3m\+\_\+functions\+::p3m\+\_\+boundary\+\_\+extrapolation (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(\+:), intent(in)}]{h, }\item[{real, dimension(\+:), intent(in)}]{u, }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+values, }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+S, }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+coef, }\item[{real, intent(in), optional}]{h\+\_\+neglect, }\item[{real, intent(in), optional}]{h\+\_\+neglect\+\_\+edge }\end{DoxyParamCaption})} Calculate the edge values and slopes at boundary cells as part of building a piecewise cubic sub-\/grid scale profiles. The following explanations apply to the left boundary cell. The same reasoning holds for the right boundary cell. A cubic needs to be built in the cell and requires four degrees of freedom, which are the edge values and slopes. The right edge values and slopes are taken to be that of the neighboring cell (i.\+e., the left edge value and slope of the neighboring cell). The left edge value and slope are determined by computing the parabola based on the cell average and the right edge value and slope. The resulting cubic is not necessarily monotonic and the slopes are subsequently modified to yield a monotonic cubic. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em n} & Number of cells\\ \hline \mbox{\tt in} & {\em h} & cell widths (size N) \mbox{[}H\mbox{]}\\ \hline \mbox{\tt in} & {\em u} & cell averages (size N) in arbitrary units \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in,out} & {\em edge\+\_\+values} & Edge value of polynomial \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in,out} & {\em ppoly\+\_\+s} & Edge slope of polynomial \mbox{[}A H-\/1\mbox{]}\\ \hline \mbox{\tt in,out} & {\em ppoly\+\_\+coef} & Coefficients of polynomial \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in} & {\em h\+\_\+neglect} & A negligibly small width for the purpose of cell reconstructions \mbox{[}H\mbox{]}\\ \hline \mbox{\tt in} & {\em h\+\_\+neglect\+\_\+edge} & A negligibly small width for the purpose of finding edge values \mbox{[}H\mbox{]} \\ \hline \end{DoxyParams} Definition at line 193 of file P3\+M\+\_\+functions.\+F90. \begin{DoxyCode} 193 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells} 194 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths (size N) [H]} 195 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell averages (size N) in arbitrary units [A]} 196 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_values\textcolor{comment}{ !< Edge value of polynomial [A]} 197 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_S\textcolor{comment}{ !< Edge slope of polynomial [A H-1]} 198 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_coef\textcolor{comment}{ !< Coefficients of polynomial [A]} 199 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\textcolor{comment}{ !< A negligibly small width for the} 200 \textcolor{comment}{ !! purpose of cell reconstructions [H]} 201 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\_edge\textcolor{comment}{ !< A negligibly small width} 202 \textcolor{comment}{ !! for the purpose of finding edge values [H]} 203 \textcolor{comment}{! Local variables} 204 \textcolor{keywordtype}{integer} :: i0, i1 205 \textcolor{keywordtype}{logical} :: monotonic \textcolor{comment}{! boolean indicating whether the cubic is monotonic} 206 \textcolor{keywordtype}{real} :: u0, u1 \textcolor{comment}{! Values of u in two adjacent cells [A]} 207 \textcolor{keywordtype}{real} :: h0, h1 \textcolor{comment}{! Values of h in two adjacent cells, plus a smal increment [H]} 208 \textcolor{keywordtype}{real} :: b, c, d \textcolor{comment}{! Temporary variables [A]} 209 \textcolor{keywordtype}{real} :: u0\_l, u0\_r \textcolor{comment}{! Left and right edge values [A]} 210 \textcolor{keywordtype}{real} :: u1\_l, u1\_r \textcolor{comment}{! Left and right edge slopes [A H-1]} 211 \textcolor{keywordtype}{real} :: slope \textcolor{comment}{! The cell center slope [A H-1]} 212 \textcolor{keywordtype}{real} :: hNeglect, hNeglect\_edge \textcolor{comment}{! Negligibly small thickness [H]} 213 214 hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect 215 hneglect\_edge = hneglect\_edge\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect\_edge)) hneglect\_edge = h\_neglect\_edge 216 217 \textcolor{comment}{! ----- Left boundary -----} 218 i0 = 1 219 i1 = 2 220 h0 = h(i0) + hneglect\_edge 221 h1 = h(i1) + hneglect\_edge 222 u0 = u(i0) 223 u1 = u(i1) 224 225 \textcolor{comment}{! Compute the left edge slope in neighboring cell and express it in} 226 \textcolor{comment}{! the global coordinate system} 227 b = ppoly\_coef(i1,2) 228 u1\_r = b / h1 \textcolor{comment}{! derivative evaluated at xi = 0.0, expressed w.r.t. x} 229 230 \textcolor{comment}{! Limit the right slope by the PLM limited slope} 231 slope = 2.0 * ( u1 - u0 ) / ( h0 + hneglect ) 232 \textcolor{keywordflow}{if} ( abs(u1\_r) > abs(slope) ) \textcolor{keywordflow}{then} 233 u1\_r = slope 234 \textcolor{keywordflow}{ endif} 235 236 \textcolor{comment}{! The right edge value in the boundary cell is taken to be the left} 237 \textcolor{comment}{! edge value in the neighboring cell} 238 u0\_r = edge\_values(i1,1) 239 240 \textcolor{comment}{! Given the right edge value and slope, we determine the left} 241 \textcolor{comment}{! edge value and slope by computing the parabola as determined by} 242 \textcolor{comment}{! the right edge value and slope and the boundary cell average} 243 u0\_l = 3.0 * u0 + 0.5 * h0*u1\_r - 2.0 * u0\_r 244 u1\_l = ( - 6.0 * u0 - 2.0 * h0*u1\_r + 6.0 * u0\_r) / ( h0 + hneglect ) 245 246 \textcolor{comment}{! Check whether the edge values are monotonic. For example, if the left edge} 247 \textcolor{comment}{! value is larger than the right edge value while the slope is positive, the} 248 \textcolor{comment}{! edge values are inconsistent and we need to modify the left edge value} 249 \textcolor{keywordflow}{if} ( (u0\_r-u0\_l) * slope < 0.0 ) \textcolor{keywordflow}{then} 250 u0\_l = u0\_r 251 u1\_l = 0.0 252 u1\_r = 0.0 253 \textcolor{keywordflow}{ endif} 254 255 \textcolor{comment}{! Store edge values and slope, build cubic and check monotonicity} 256 edge\_values(i0,1) = u0\_l 257 edge\_values(i0,2) = u0\_r 258 ppoly\_s(i0,1) = u1\_l 259 ppoly\_s(i0,2) = u1\_r 260 261 \textcolor{comment}{! Store edge values and slope, build cubic and check monotonicity} 262 \textcolor{keyword}{call }build\_cubic\_interpolant( h, i0, edge\_values, ppoly\_s, ppoly\_coef ) 263 monotonic = is\_cubic\_monotonic( ppoly\_coef, i0 ) 264 265 \textcolor{keywordflow}{if} ( .not.monotonic ) \textcolor{keywordflow}{then} 266 \textcolor{keyword}{call }monotonize\_cubic( h0, u0\_l, u0\_r, 0.0, slope, slope, u1\_l, u1\_r ) 267 268 \textcolor{comment}{! Rebuild cubic after monotonization} 269 ppoly\_s(i0,1) = u1\_l 270 ppoly\_s(i0,2) = u1\_r 271 \textcolor{keyword}{call }build\_cubic\_interpolant( h, i0, edge\_values, ppoly\_s, ppoly\_coef ) 272 273 \textcolor{keywordflow}{ endif} 274 275 \textcolor{comment}{! ----- Right boundary -----} 276 i0 = n-1 277 i1 = n 278 h0 = h(i0) + hneglect\_edge 279 h1 = h(i1) + hneglect\_edge 280 u0 = u(i0) 281 u1 = u(i1) 282 283 \textcolor{comment}{! Compute the right edge slope in neighboring cell and express it in} 284 \textcolor{comment}{! the global coordinate system} 285 b = ppoly\_coef(i0,2) 286 c = ppoly\_coef(i0,3) 287 d = ppoly\_coef(i0,4) 288 u1\_l = (b + 2*c + 3*d) / ( h0 + hneglect ) \textcolor{comment}{! derivative evaluated at xi = 1.0} 289 290 \textcolor{comment}{! Limit the left slope by the PLM limited slope} 291 slope = 2.0 * ( u1 - u0 ) / ( h1 + hneglect ) 292 \textcolor{keywordflow}{if} ( abs(u1\_l) > abs(slope) ) \textcolor{keywordflow}{then} 293 u1\_l = slope 294 \textcolor{keywordflow}{ endif} 295 296 \textcolor{comment}{! The left edge value in the boundary cell is taken to be the right} 297 \textcolor{comment}{! edge value in the neighboring cell} 298 u0\_l = edge\_values(i0,2) 299 300 \textcolor{comment}{! Given the left edge value and slope, we determine the right} 301 \textcolor{comment}{! edge value and slope by computing the parabola as determined by} 302 \textcolor{comment}{! the left edge value and slope and the boundary cell average} 303 u0\_r = 3.0 * u1 - 0.5 * h1*u1\_l - 2.0 * u0\_l 304 u1\_r = ( 6.0 * u1 - 2.0 * h1*u1\_l - 6.0 * u0\_l) / ( h1 + hneglect ) 305 306 \textcolor{comment}{! Check whether the edge values are monotonic. For example, if the right edge} 307 \textcolor{comment}{! value is smaller than the left edge value while the slope is positive, the} 308 \textcolor{comment}{! edge values are inconsistent and we need to modify the right edge value} 309 \textcolor{keywordflow}{if} ( (u0\_r-u0\_l) * slope < 0.0 ) \textcolor{keywordflow}{then} 310 u0\_r = u0\_l 311 u1\_l = 0.0 312 u1\_r = 0.0 313 \textcolor{keywordflow}{ endif} 314 315 \textcolor{comment}{! Store edge values and slope, build cubic and check monotonicity} 316 edge\_values(i1,1) = u0\_l 317 edge\_values(i1,2) = u0\_r 318 ppoly\_s(i1,1) = u1\_l 319 ppoly\_s(i1,2) = u1\_r 320 321 \textcolor{keyword}{call }build\_cubic\_interpolant( h, i1, edge\_values, ppoly\_s, ppoly\_coef ) 322 monotonic = is\_cubic\_monotonic( ppoly\_coef, i1 ) 323 324 \textcolor{keywordflow}{if} ( .not.monotonic ) \textcolor{keywordflow}{then} 325 \textcolor{keyword}{call }monotonize\_cubic( h1, u0\_l, u0\_r, slope, 0.0, slope, u1\_l, u1\_r ) 326 327 \textcolor{comment}{! Rebuild cubic after monotonization} 328 ppoly\_s(i1,1) = u1\_l 329 ppoly\_s(i1,2) = u1\_r 330 \textcolor{keyword}{call }build\_cubic\_interpolant( h, i1, edge\_values, ppoly\_s, ppoly\_coef ) 331 332 \textcolor{keywordflow}{ endif} 333 \end{DoxyCode} \mbox{\Hypertarget{namespacep3m__functions_a685d3ef292536dae810b2059ccaa5819}\label{namespacep3m__functions_a685d3ef292536dae810b2059ccaa5819}} \index{p3m\+\_\+functions@{p3m\+\_\+functions}!p3m\+\_\+interpolation@{p3m\+\_\+interpolation}} \index{p3m\+\_\+interpolation@{p3m\+\_\+interpolation}!p3m\+\_\+functions@{p3m\+\_\+functions}} \subsubsection{\texorpdfstring{p3m\+\_\+interpolation()}{p3m\_interpolation()}} {\footnotesize\ttfamily subroutine, public p3m\+\_\+functions\+::p3m\+\_\+interpolation (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(\+:), intent(in)}]{h, }\item[{real, dimension(\+:), intent(in)}]{u, }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+values, }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+S, }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+coef, }\item[{real, intent(in), optional}]{h\+\_\+neglect, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})} Set up a piecewise cubic interpolation from cell averages and estimated edge slopes and values. Cubic interpolation between edges. The edge values and slopes M\+U\+ST have been estimated prior to calling this routine. It is assumed that the size of the array \textquotesingle{}u\textquotesingle{} is equal to the number of cells defining \textquotesingle{}grid\textquotesingle{} and \textquotesingle{}ppoly\textquotesingle{}. No consistency check is performed here. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em n} & Number of cells\\ \hline \mbox{\tt in} & {\em h} & cell widths (size N) \mbox{[}H\mbox{]}\\ \hline \mbox{\tt in} & {\em u} & cell averages (size N) in arbitrary units \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in,out} & {\em edge\+\_\+values} & Edge value of polynomial \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in,out} & {\em ppoly\+\_\+s} & Edge slope of polynomial \mbox{[}A H-\/1\mbox{]}.\\ \hline \mbox{\tt in,out} & {\em ppoly\+\_\+coef} & Coefficients of polynomial \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in} & {\em h\+\_\+neglect} & A negligibly small width for the purpose of cell reconstructions \mbox{[}H\mbox{]}\\ \hline \mbox{\tt in} & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\ \hline \end{DoxyParams} Definition at line 29 of file P3\+M\+\_\+functions.\+F90. \begin{DoxyCode} 29 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells} 30 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths (size N) [H]} 31 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell averages (size N) in arbitrary units [A]} 32 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_values\textcolor{comment}{ !< Edge value of polynomial [A]} 33 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_S\textcolor{comment}{ !< Edge slope of polynomial [A H-1].} 34 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_coef\textcolor{comment}{ !< Coefficients of polynomial [A]} 35 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\textcolor{comment}{ !< A negligibly small width for the} 36 \textcolor{comment}{ !! purpose of cell reconstructions [H]} 37 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.} 38 39 \textcolor{comment}{! Call the limiter for p3m, which takes care of everything from} 40 \textcolor{comment}{! computing the coefficients of the cubic to monotonizing it.} 41 \textcolor{comment}{! This routine could be called directly instead of having to call} 42 \textcolor{comment}{! 'P3M\_interpolation' first but we do that to provide an homogeneous} 43 \textcolor{comment}{! interface.} 44 \textcolor{keyword}{call }p3m\_limiter( n, h, u, edge\_values, ppoly\_s, ppoly\_coef, h\_neglect, answers\_2018 ) 45 \end{DoxyCode} \mbox{\Hypertarget{namespacep3m__functions_a44db59cb5df7f139e05b745746342fcf}\label{namespacep3m__functions_a44db59cb5df7f139e05b745746342fcf}} \index{p3m\+\_\+functions@{p3m\+\_\+functions}!p3m\+\_\+limiter@{p3m\+\_\+limiter}} \index{p3m\+\_\+limiter@{p3m\+\_\+limiter}!p3m\+\_\+functions@{p3m\+\_\+functions}} \subsubsection{\texorpdfstring{p3m\+\_\+limiter()}{p3m\_limiter()}} {\footnotesize\ttfamily subroutine p3m\+\_\+functions\+::p3m\+\_\+limiter (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(\+:), intent(in)}]{h, }\item[{real, dimension(\+:), intent(in)}]{u, }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+values, }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+S, }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+coef, }\item[{real, intent(in), optional}]{h\+\_\+neglect, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}} Adust a piecewise cubic reconstruction with a limiter that adjusts the edge values and slopes. The p3m limiter operates as follows\+: \begin{DoxyEnumerate} \item Edge values are bounded \item Discontinuous edge values are systematically averaged \item Loop on cells and do the following a. Build cubic curve b. Check if cubic curve is monotonic c. If not, monotonize cubic curve and rebuild it \end{DoxyEnumerate} Step 3 of the monotonization process leaves all edge values unchanged. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em n} & Number of cells\\ \hline \mbox{\tt in} & {\em h} & cell widths (size N) \mbox{[}H\mbox{]}\\ \hline \mbox{\tt in} & {\em u} & cell averages (size N) in arbitrary units \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in,out} & {\em edge\+\_\+values} & Edge value of polynomial \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in,out} & {\em ppoly\+\_\+s} & Edge slope of polynomial \mbox{[}A H-\/1\mbox{]}\\ \hline \mbox{\tt in,out} & {\em ppoly\+\_\+coef} & Coefficients of polynomial \mbox{[}A\mbox{]}\\ \hline \mbox{\tt in} & {\em h\+\_\+neglect} & A negligibly small width for the purpose of cell reconstructions \mbox{[}H\mbox{]}\\ \hline \mbox{\tt in} & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\ \hline \end{DoxyParams} Definition at line 62 of file P3\+M\+\_\+functions.\+F90. \begin{DoxyCode} 62 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells} 63 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths (size N) [H]} 64 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell averages (size N) in arbitrary units [A]} 65 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_values\textcolor{comment}{ !< Edge value of polynomial [A]} 66 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_S\textcolor{comment}{ !< Edge slope of polynomial [A H-1]} 67 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_coef\textcolor{comment}{ !< Coefficients of polynomial [A]} 68 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\textcolor{comment}{ !< A negligibly small width for} 69 \textcolor{comment}{ !! the purpose of cell reconstructions [H]} 70 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.} 71 72 \textcolor{comment}{! Local variables} 73 \textcolor{keywordtype}{integer} :: k \textcolor{comment}{! loop index} 74 \textcolor{keywordtype}{logical} :: monotonic \textcolor{comment}{! boolean indicating whether the cubic is monotonic} 75 \textcolor{keywordtype}{real} :: u0\_l, u0\_r \textcolor{comment}{! edge values [A]} 76 \textcolor{keywordtype}{real} :: u1\_l, u1\_r \textcolor{comment}{! edge slopes [A H-1]} 77 \textcolor{keywordtype}{real} :: u\_l, u\_c, u\_r \textcolor{comment}{! left, center and right cell averages [A]} 78 \textcolor{keywordtype}{real} :: h\_l, h\_c, h\_r \textcolor{comment}{! left, center and right cell widths [H]} 79 \textcolor{keywordtype}{real} :: sigma\_l, sigma\_c, sigma\_r \textcolor{comment}{! left, center and right van Leer slopes [A H-1]} 80 \textcolor{keywordtype}{real} :: slope \textcolor{comment}{! retained PLM slope [A H-1]} 81 \textcolor{keywordtype}{real} :: eps 82 \textcolor{keywordtype}{real} :: hNeglect \textcolor{comment}{! A negligibly small thickness [H]} 83 84 hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect 85 86 eps = 1e-10 87 88 \textcolor{comment}{! 1. Bound edge values (boundary cells are assumed to be local extrema)} 89 \textcolor{keyword}{call }bound\_edge\_values( n, h, u, edge\_values, hneglect, answers\_2018 ) 90 91 \textcolor{comment}{! 2. Systematically average discontinuous edge values} 92 \textcolor{keyword}{call }average\_discontinuous\_edge\_values( n, edge\_values ) 93 94 95 \textcolor{comment}{! 3. Loop on cells and do the following} 96 \textcolor{comment}{! (a) Build cubic curve} 97 \textcolor{comment}{! (b) Check if cubic curve is monotonic} 98 \textcolor{comment}{! (c) If not, monotonize cubic curve and rebuild it} 99 \textcolor{keywordflow}{do} k = 1,n 100 101 \textcolor{comment}{! Get edge values, edge slopes and cell width} 102 u0\_l = edge\_values(k,1) 103 u0\_r = edge\_values(k,2) 104 u1\_l = ppoly\_s(k,1) 105 u1\_r = ppoly\_s(k,2) 106 107 \textcolor{comment}{! Get cell widths and cell averages (boundary cells are assumed to} 108 \textcolor{comment}{! be local extrema for the sake of slopes)} 109 u\_c = u(k) 110 h\_c = h(k) 111 112 \textcolor{keywordflow}{if} ( k == 1 ) \textcolor{keywordflow}{then} 113 h\_l = h(k) 114 u\_l = u(k) 115 \textcolor{keywordflow}{else} 116 h\_l = h(k-1) 117 u\_l = u(k-1) 118 \textcolor{keywordflow}{ endif} 119 120 \textcolor{keywordflow}{if} ( k == n ) \textcolor{keywordflow}{then} 121 h\_r = h(k) 122 u\_r = u(k) 123 \textcolor{keywordflow}{else} 124 h\_r = h(k+1) 125 u\_r = u(k+1) 126 \textcolor{keywordflow}{ endif} 127 128 \textcolor{comment}{! Compute limited slope} 129 sigma\_l = 2.0 * ( u\_c - u\_l ) / ( h\_c + hneglect ) 130 sigma\_c = 2.0 * ( u\_r - u\_l ) / ( h\_l + 2.0*h\_c + h\_r + hneglect ) 131 sigma\_r = 2.0 * ( u\_r - u\_c ) / ( h\_c + hneglect ) 132 133 \textcolor{keywordflow}{if} ( (sigma\_l * sigma\_r) > 0.0 ) \textcolor{keywordflow}{then} 134 slope = sign( min(abs(sigma\_l),abs(sigma\_c),abs(sigma\_r)), sigma\_c ) 135 \textcolor{keywordflow}{else} 136 slope = 0.0 137 \textcolor{keywordflow}{ endif} 138 139 \textcolor{comment}{! If the slopes are small, set them to zero to prevent asymmetric representation near extrema.} 140 \textcolor{keywordflow}{if} ( abs(u1\_l*h\_c) < epsilon(u\_c)*abs(u\_c) ) u1\_l = 0.0 141 \textcolor{keywordflow}{if} ( abs(u1\_r*h\_c) < epsilon(u\_c)*abs(u\_c) ) u1\_r = 0.0 142 143 \textcolor{comment}{! The edge slopes are limited from above by the respective} 144 \textcolor{comment}{! one-sided slopes} 145 \textcolor{keywordflow}{if} ( abs(u1\_l) > abs(sigma\_l) ) \textcolor{keywordflow}{then} 146 u1\_l = sigma\_l 147 \textcolor{keywordflow}{ endif} 148 149 \textcolor{keywordflow}{if} ( abs(u1\_r) > abs(sigma\_r) ) \textcolor{keywordflow}{then} 150 u1\_r = sigma\_r 151 \textcolor{keywordflow}{ endif} 152 153 \textcolor{comment}{! Build cubic interpolant (compute the coefficients)} 154 \textcolor{keyword}{call }build\_cubic\_interpolant( h, k, edge\_values, ppoly\_s, ppoly\_coef ) 155 156 \textcolor{comment}{! Check whether cubic is monotonic} 157 monotonic = is\_cubic\_monotonic( ppoly\_coef, k ) 158 159 \textcolor{comment}{! If cubic is not monotonic, monotonize it by modifiying the} 160 \textcolor{comment}{! edge slopes, store the new edge slopes and recompute the} 161 \textcolor{comment}{! cubic coefficients} 162 \textcolor{keywordflow}{if} ( .not.monotonic ) \textcolor{keywordflow}{then} 163 \textcolor{keyword}{call }monotonize\_cubic( h\_c, u0\_l, u0\_r, sigma\_l, sigma\_r, slope, u1\_l, u1\_r ) 164 \textcolor{keywordflow}{ endif} 165 166 \textcolor{comment}{! Store edge slopes} 167 ppoly\_s(k,1) = u1\_l 168 ppoly\_s(k,2) = u1\_r 169 170 \textcolor{comment}{! Recompute coefficients of cubic} 171 \textcolor{keyword}{call }build\_cubic\_interpolant( h, k, edge\_values, ppoly\_s, ppoly\_coef ) 172 173 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! loop on cells} 174 \end{DoxyCode}