\hypertarget{namespaceregrid__edge__values}{}\doxysection{regrid\+\_\+edge\+\_\+values Module Reference} \label{namespaceregrid__edge__values}\index{regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsection{Detailed Description} Edge value estimation for high-\/order resconstruction. \doxysubsection*{Functions/\+Subroutines} \begin{DoxyCompactItemize} \item subroutine, public \mbox{\hyperlink{namespaceregrid__edge__values_a8480933738145ee2525e1bf449f14097}{bound\+\_\+edge\+\_\+values}} (N, h, u, edge\+\_\+val, h\+\_\+neglect, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Bound edge values by neighboring cell averages. \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespaceregrid__edge__values_ad43eb7fa3a284e2b33068f47232521ca}{average\+\_\+discontinuous\+\_\+edge\+\_\+values}} (N, edge\+\_\+val) \begin{DoxyCompactList}\small\item\em Replace discontinuous collocated edge values with their average. \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespaceregrid__edge__values_a4a496536a77bef8467c441731619761d}{check\+\_\+discontinuous\+\_\+edge\+\_\+values}} (N, u, edge\+\_\+val) \begin{DoxyCompactList}\small\item\em Check discontinuous edge values and replace them with their average if not monotonic. \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespaceregrid__edge__values_a8dc053bb28658c7091f1f5a82184019a}{edge\+\_\+values\+\_\+explicit\+\_\+h2}} (N, h, u, edge\+\_\+val) \begin{DoxyCompactList}\small\item\em Compute h2 edge values (explicit second order accurate) in the same units as u. \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespaceregrid__edge__values_a2fa3a70b208abe68f29bef9589081209}{edge\+\_\+values\+\_\+explicit\+\_\+h4}} (N, h, u, edge\+\_\+val, h\+\_\+neglect, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Compute h4 edge values (explicit fourth order accurate) in the same units as u. \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespaceregrid__edge__values_a0c268712aaf87b3597cef51c85fb32cb}{edge\+\_\+values\+\_\+implicit\+\_\+h4}} (N, h, u, edge\+\_\+val, h\+\_\+neglect, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Compute ih4 edge values (implicit fourth order accurate) in the same units as u. \end{DoxyCompactList}\item subroutine \mbox{\hyperlink{namespaceregrid__edge__values_a7082604865f485ac18af8b44a7aa0224}{end\+\_\+value\+\_\+h4}} (dz, u, Csys) \begin{DoxyCompactList}\small\item\em Determine a one-\/sided 4th order polynomial fit of u to the data points for the purposes of specifying edge values, as described in the appendix of White and Adcroft J\+CP 2008. \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespaceregrid__edge__values_ac33bba0097d34e462abda17e78847862}{edge\+\_\+slopes\+\_\+implicit\+\_\+h3}} (N, h, u, edge\+\_\+slopes, h\+\_\+neglect, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Compute ih3 edge slopes (implicit third order accurate) in the same units as h. \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespaceregrid__edge__values_a7c3fb27d99df9e3594ff32d1d03e1b34}{edge\+\_\+slopes\+\_\+implicit\+\_\+h5}} (N, h, u, edge\+\_\+slopes, h\+\_\+neglect, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Compute ih5 edge slopes (implicit fifth order accurate) \end{DoxyCompactList}\item subroutine, public \mbox{\hyperlink{namespaceregrid__edge__values_a9955c45dcd1bfec32fbf5602315cb5b1}{edge\+\_\+values\+\_\+implicit\+\_\+h6}} (N, h, u, edge\+\_\+val, h\+\_\+neglect, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Compute ih6 edge values (implicit sixth order accurate) in the same units as u. \end{DoxyCompactList}\item subroutine \mbox{\hyperlink{namespaceregrid__edge__values_a9f0ab1f167eb57cbaca668f7a11b31a1}{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag}} (Al, Ac, Au, R, X, N) \begin{DoxyCompactList}\small\item\em Solve the tridiagonal system AX = R. \end{DoxyCompactList}\item subroutine \mbox{\hyperlink{namespaceregrid__edge__values_a09a4c9f369606d81580116620975352e}{linear\+\_\+solver}} (N, A, R, X) \begin{DoxyCompactList}\small\item\em Solve the linear system AX = R by Gaussian elimination. \end{DoxyCompactList}\item subroutine \mbox{\hyperlink{namespaceregrid__edge__values_aef5443593d35636eb95de2b580878a20}{test\+\_\+line}} (msg, N, A, C, R, mag, tol) \begin{DoxyCompactList}\small\item\em Test that A$\ast$C = R to within a tolerance, issuing a fatal error with an explanatory message if they do not. \end{DoxyCompactList}\end{DoxyCompactItemize} \doxysubsection*{Variables} \begin{DoxyCompactItemize} \item \mbox{\Hypertarget{namespaceregrid__edge__values_a241fc99294d87afd8146f4b3e48dc021}\label{namespaceregrid__edge__values_a241fc99294d87afd8146f4b3e48dc021}} real, parameter \mbox{\hyperlink{namespaceregrid__edge__values_a241fc99294d87afd8146f4b3e48dc021}{hneglect\+\_\+edge\+\_\+dflt}} = 1.e-\/10 \begin{DoxyCompactList}\small\item\em The default value for cut-\/off minimum thickness for sum(h) in edge value inversions. \end{DoxyCompactList}\item \mbox{\Hypertarget{namespaceregrid__edge__values_a2800ef300c67ac316d81e5fef12907b5}\label{namespaceregrid__edge__values_a2800ef300c67ac316d81e5fef12907b5}} real, parameter \mbox{\hyperlink{namespaceregrid__edge__values_a2800ef300c67ac316d81e5fef12907b5}{hneglect\+\_\+dflt}} = 1.e-\/30 \begin{DoxyCompactList}\small\item\em The default value for cut-\/off minimum thickness for sum(h) in other calculations. \end{DoxyCompactList}\item \mbox{\Hypertarget{namespaceregrid__edge__values_a6e31d54d55af5650f72574cf82eb0f1f}\label{namespaceregrid__edge__values_a6e31d54d55af5650f72574cf82eb0f1f}} real, parameter \mbox{\hyperlink{namespaceregrid__edge__values_a6e31d54d55af5650f72574cf82eb0f1f}{hminfrac}} = 1.e-\/5 \begin{DoxyCompactList}\small\item\em A minimum fraction for min(h)/sum(h) \end{DoxyCompactList}\end{DoxyCompactItemize} \doxysubsection{Function/\+Subroutine Documentation} \mbox{\Hypertarget{namespaceregrid__edge__values_ad43eb7fa3a284e2b33068f47232521ca}\label{namespaceregrid__edge__values_ad43eb7fa3a284e2b33068f47232521ca}} \index{regrid\_edge\_values@{regrid\_edge\_values}!average\_discontinuous\_edge\_values@{average\_discontinuous\_edge\_values}} \index{average\_discontinuous\_edge\_values@{average\_discontinuous\_edge\_values}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{average\_discontinuous\_edge\_values()}{average\_discontinuous\_edge\_values()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+edge\+\_\+values\+::average\+\_\+discontinuous\+\_\+edge\+\_\+values (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n,2), intent(inout)}]{edge\+\_\+val }\end{DoxyParamCaption})} Replace discontinuous collocated edge values with their average. For each interior edge, check whether the edge values are discontinuous. If so, compute the average and replace the edge values by the average. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & Number of cells \\ \hline \mbox{\texttt{ in,out}} & {\em edge\+\_\+val} & Edge values that may be modified \mbox{[}A\mbox{]}; the second index is for the two edges of each cell. \\ \hline \end{DoxyParams} Definition at line 116 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{116 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells}} \DoxyCodeLine{117 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Edge values that may be modified [A]; the}} \DoxyCodeLine{118 \textcolor{comment}{ !! second index is for the two edges of each cell.}} \DoxyCodeLine{119 \textcolor{comment}{! Local variables}} \DoxyCodeLine{120 \textcolor{keywordtype}{integer} :: k \textcolor{comment}{! loop index}} \DoxyCodeLine{121 \textcolor{keywordtype}{ real} :: u0\_avg \textcolor{comment}{! avg value at given edge}} \DoxyCodeLine{122 } \DoxyCodeLine{123 \textcolor{comment}{! Loop on interior edges}} \DoxyCodeLine{124 \textcolor{keywordflow}{do} k = 1,n-\/1} \DoxyCodeLine{125 \textcolor{comment}{! Compare edge values on the right and left sides of the edge}} \DoxyCodeLine{126 \textcolor{keywordflow}{if} ( edge\_val(k,2) /= edge\_val(k+1,1) ) \textcolor{keywordflow}{then}} \DoxyCodeLine{127 u0\_avg = 0.5 * ( edge\_val(k,2) + edge\_val(k+1,1) )} \DoxyCodeLine{128 edge\_val(k,2) = u0\_avg} \DoxyCodeLine{129 edge\_val(k+1,1) = u0\_avg} \DoxyCodeLine{130 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{131 } \DoxyCodeLine{132 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on interior edges}} \DoxyCodeLine{133 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a8480933738145ee2525e1bf449f14097}\label{namespaceregrid__edge__values_a8480933738145ee2525e1bf449f14097}} \index{regrid\_edge\_values@{regrid\_edge\_values}!bound\_edge\_values@{bound\_edge\_values}} \index{bound\_edge\_values@{bound\_edge\_values}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{bound\_edge\_values()}{bound\_edge\_values()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+edge\+\_\+values\+::bound\+\_\+edge\+\_\+values (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n), intent(in)}]{h, }\item[{real, dimension(n), intent(in)}]{u, }\item[{real, dimension(n,2), intent(inout)}]{edge\+\_\+val, }\item[{real, intent(in), optional}]{h\+\_\+neglect, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})} Bound edge values by neighboring cell averages. In this routine, we loop on all cells to bound their left and right edge values by the cell averages. That is, the left edge value must lie between the left cell average and the central cell average. A similar reasoning applies to the right edge values. Both boundary edge values are set equal to the boundary cell averages. Any extrapolation scheme is applied after this routine has been called. Therefore, boundary cells are treated as if they were local extrama. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & Number of cells \\ \hline \mbox{\texttt{ in}} & {\em h} & cell widths \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em edge\+\_\+val} & Potentially modified edge values \mbox{[}A\mbox{]}; the second index is for the two edges of each cell. \\ \hline \mbox{\texttt{ in}} & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\ \hline \end{DoxyParams} Definition at line 45 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{45 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells}} \DoxyCodeLine{46 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths [H]}} \DoxyCodeLine{47 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}} \DoxyCodeLine{48 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Potentially modified edge values [A]; the}} \DoxyCodeLine{49 \textcolor{comment}{ !! second index is for the two edges of each cell.}} \DoxyCodeLine{50 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}} \DoxyCodeLine{51 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}} \DoxyCodeLine{52 \textcolor{comment}{! Local variables}} \DoxyCodeLine{53 \textcolor{keywordtype}{ real} :: sigma\_l, sigma\_c, sigma\_r \textcolor{comment}{! left, center and right van Leer slopes [A H-\/1] or [A]}} \DoxyCodeLine{54 \textcolor{keywordtype}{ real} :: slope\_x\_h \textcolor{comment}{! retained PLM slope times half grid step [A]}} \DoxyCodeLine{55 \textcolor{keywordtype}{ real} :: hNeglect \textcolor{comment}{! A negligible thickness [H].}} \DoxyCodeLine{56 \textcolor{keywordtype}{logical} :: use\_2018\_answers \textcolor{comment}{! If true use older, less acccurate expressions.}} \DoxyCodeLine{57 \textcolor{keywordtype}{integer} :: k, km1, kp1 \textcolor{comment}{! Loop index and the values to either side.}} \DoxyCodeLine{58 } \DoxyCodeLine{59 use\_2018\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) use\_2018\_answers = answers\_2018} \DoxyCodeLine{60 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{61 hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect} \DoxyCodeLine{62 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{63 } \DoxyCodeLine{64 \textcolor{comment}{! Loop on cells to bound edge value}} \DoxyCodeLine{65 \textcolor{keywordflow}{do} k = 1,n} \DoxyCodeLine{66 } \DoxyCodeLine{67 \textcolor{comment}{! For the sake of bounding boundary edge values, the left neighbor of the left boundary cell}} \DoxyCodeLine{68 \textcolor{comment}{! is assumed to be the same as the left boundary cell and the right neighbor of the right}} \DoxyCodeLine{69 \textcolor{comment}{! boundary cell is assumed to be the same as the right boundary cell. This effectively makes}} \DoxyCodeLine{70 \textcolor{comment}{! boundary cells look like extrema.}} \DoxyCodeLine{71 km1 = max(1,k-\/1) ; kp1 = min(k+1,n)} \DoxyCodeLine{72 } \DoxyCodeLine{73 slope\_x\_h = 0.0} \DoxyCodeLine{74 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{75 sigma\_l = 2.0 * ( u(k) -\/ u(km1) ) / ( h(k) + hneglect )} \DoxyCodeLine{76 sigma\_c = 2.0 * ( u(kp1) -\/ u(km1) ) / ( h(km1) + 2.0*h(k) + h(kp1) + hneglect )} \DoxyCodeLine{77 sigma\_r = 2.0 * ( u(kp1) -\/ u(k) ) / ( h(k) + hneglect )} \DoxyCodeLine{78 } \DoxyCodeLine{79 \textcolor{comment}{! The limiter is used in the local coordinate system to each cell, so for convenience store}} \DoxyCodeLine{80 \textcolor{comment}{! the slope times a half grid spacing. (See White and Adcroft JCP 2008 Eqs 19 and 20)}} \DoxyCodeLine{81 \textcolor{keywordflow}{if} ( (sigma\_l * sigma\_r) > 0.0 ) \&} \DoxyCodeLine{82 slope\_x\_h = 0.5 * h(k) * sign( min(abs(sigma\_l),abs(sigma\_c),abs(sigma\_r)), sigma\_c )} \DoxyCodeLine{83 \textcolor{keywordflow}{elseif} ( ((h(km1) + h(kp1)) + 2.0*h(k)) > 0.0 ) \textcolor{keywordflow}{then}} \DoxyCodeLine{84 sigma\_l = ( u(k) -\/ u(km1) )} \DoxyCodeLine{85 sigma\_c = ( u(kp1) -\/ u(km1) ) * ( h(k) / ((h(km1) + h(kp1)) + 2.0*h(k)) )} \DoxyCodeLine{86 sigma\_r = ( u(kp1) -\/ u(k) )} \DoxyCodeLine{87 } \DoxyCodeLine{88 \textcolor{comment}{! The limiter is used in the local coordinate system to each cell, so for convenience store}} \DoxyCodeLine{89 \textcolor{comment}{! the slope times a half grid spacing. (See White and Adcroft JCP 2008 Eqs 19 and 20)}} \DoxyCodeLine{90 \textcolor{keywordflow}{if} ( (sigma\_l * sigma\_r) > 0.0 ) \&} \DoxyCodeLine{91 slope\_x\_h = sign( min(abs(sigma\_l),abs(sigma\_c),abs(sigma\_r)), sigma\_c )} \DoxyCodeLine{92 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{93 } \DoxyCodeLine{94 \textcolor{comment}{! Limit the edge values}} \DoxyCodeLine{95 \textcolor{keywordflow}{if} ( (u(km1)-\/edge\_val(k,1)) * (edge\_val(k,1)-\/u(k)) < 0.0 ) \textcolor{keywordflow}{then}} \DoxyCodeLine{96 edge\_val(k,1) = u(k) -\/ sign( min( abs(slope\_x\_h), abs(edge\_val(k,1)-\/u(k)) ), slope\_x\_h )} \DoxyCodeLine{97 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{98 } \DoxyCodeLine{99 \textcolor{keywordflow}{if} ( (u(kp1)-\/edge\_val(k,2)) * (edge\_val(k,2)-\/u(k)) < 0.0 ) \textcolor{keywordflow}{then}} \DoxyCodeLine{100 edge\_val(k,2) = u(k) + sign( min( abs(slope\_x\_h), abs(edge\_val(k,2)-\/u(k)) ), slope\_x\_h )} \DoxyCodeLine{101 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{102 } \DoxyCodeLine{103 \textcolor{comment}{! Finally bound by neighboring cell means in case of roundoff}} \DoxyCodeLine{104 edge\_val(k,1) = max( min( edge\_val(k,1), max(u(km1), u(k)) ), min(u(km1), u(k)) )} \DoxyCodeLine{105 edge\_val(k,2) = max( min( edge\_val(k,2), max(u(kp1), u(k)) ), min(u(kp1), u(k)) )} \DoxyCodeLine{106 } \DoxyCodeLine{107 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! loop on interior edges}} \DoxyCodeLine{108 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a4a496536a77bef8467c441731619761d}\label{namespaceregrid__edge__values_a4a496536a77bef8467c441731619761d}} \index{regrid\_edge\_values@{regrid\_edge\_values}!check\_discontinuous\_edge\_values@{check\_discontinuous\_edge\_values}} \index{check\_discontinuous\_edge\_values@{check\_discontinuous\_edge\_values}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{check\_discontinuous\_edge\_values()}{check\_discontinuous\_edge\_values()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+edge\+\_\+values\+::check\+\_\+discontinuous\+\_\+edge\+\_\+values (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n), intent(in)}]{u, }\item[{real, dimension(n,2), intent(inout)}]{edge\+\_\+val }\end{DoxyParamCaption})} Check discontinuous edge values and replace them with their average if not monotonic. For each interior edge, check whether the edge values are discontinuous. If so and if they are not monotonic, replace each edge value by their average. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & Number of cells \\ \hline \mbox{\texttt{ in}} & {\em u} & cell averages in arbitrary units \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em edge\+\_\+val} & Cell edge values \mbox{[}A\mbox{]}; the second index is for the two edges of each cell. \\ \hline \end{DoxyParams} Definition at line 141 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{141 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells}} \DoxyCodeLine{142 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell averages in arbitrary units [A]}} \DoxyCodeLine{143 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Cell edge values [A]; the}} \DoxyCodeLine{144 \textcolor{comment}{ !! second index is for the two edges of each cell.}} \DoxyCodeLine{145 \textcolor{comment}{! Local variables}} \DoxyCodeLine{146 \textcolor{keywordtype}{integer} :: k \textcolor{comment}{! loop index}} \DoxyCodeLine{147 \textcolor{keywordtype}{ real} :: u0\_avg \textcolor{comment}{! avg value at given edge [A]}} \DoxyCodeLine{148 } \DoxyCodeLine{149 \textcolor{keywordflow}{do} k = 1,n-\/1} \DoxyCodeLine{150 \textcolor{keywordflow}{if} ( (edge\_val(k+1,1) -\/ edge\_val(k,2)) * (u(k+1) -\/ u(k)) < 0.0 ) \textcolor{keywordflow}{then}} \DoxyCodeLine{151 u0\_avg = 0.5 * ( edge\_val(k,2) + edge\_val(k+1,1) )} \DoxyCodeLine{152 u0\_avg = max( min( u0\_avg, max(u(k), u(k+1)) ), min(u(k), u(k+1)) )} \DoxyCodeLine{153 edge\_val(k,2) = u0\_avg} \DoxyCodeLine{154 edge\_val(k+1,1) = u0\_avg} \DoxyCodeLine{155 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{156 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on interior edges}} \DoxyCodeLine{157 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_ac33bba0097d34e462abda17e78847862}\label{namespaceregrid__edge__values_ac33bba0097d34e462abda17e78847862}} \index{regrid\_edge\_values@{regrid\_edge\_values}!edge\_slopes\_implicit\_h3@{edge\_slopes\_implicit\_h3}} \index{edge\_slopes\_implicit\_h3@{edge\_slopes\_implicit\_h3}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{edge\_slopes\_implicit\_h3()}{edge\_slopes\_implicit\_h3()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+edge\+\_\+values\+::edge\+\_\+slopes\+\_\+implicit\+\_\+h3 (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n), intent(in)}]{h, }\item[{real, dimension(n), intent(in)}]{u, }\item[{real, dimension(n,2), intent(inout)}]{edge\+\_\+slopes, }\item[{real, intent(in), optional}]{h\+\_\+neglect, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})} Compute ih3 edge slopes (implicit third order accurate) in the same units as h. Compute edge slopes based on third-\/order implicit estimates. Note that the estimates are fourth-\/order accurate on uniform grids Third-\/order implicit estimates of edge slopes are based on a two-\/cell stencil. A tridiagonal system is set up and is based on expressing the edge slopes in terms of neighboring cell averages. The generic relationship is \[ \alpha u\textnormal{\textquotesingle}_{i-1/2} + u\textnormal{\textquotesingle}_{i+1/2} + \beta u\textnormal{\textquotesingle}_{i+3/2} = a \bar{u}_i + b \bar{u}_{i+1} \] and the stencil looks like this \begin{DoxyVerb} i i+1 \end{DoxyVerb} ..--o-\/-\/-\/---o-\/-\/-\/---o--.. i-\/1/2 i+1/2 i+3/2 In this routine, the coefficients $\alpha$, $\beta$, a and b are computed, the tridiagonal system is built, boundary conditions are prescribed and the system is solved to yield edge-\/slope estimates. There are N+1 unknowns and we are able to write N-\/1 equations. The boundary conditions close the system. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & Number of cells \\ \hline \mbox{\texttt{ in}} & {\em h} & cell widths \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em edge\+\_\+slopes} & Returned edge slopes \mbox{[}A H-\/1\mbox{]}; the second index is for the two edges of each cell. \\ \hline \mbox{\texttt{ in}} & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\ \hline \end{DoxyParams} Definition at line 697 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{697 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells}} \DoxyCodeLine{698 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths [H]}} \DoxyCodeLine{699 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}} \DoxyCodeLine{700 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_slopes\textcolor{comment}{ !< Returned edge slopes [A H-\/1]; the}} \DoxyCodeLine{701 \textcolor{comment}{ !! second index is for the two edges of each cell.}} \DoxyCodeLine{702 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}} \DoxyCodeLine{703 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}} \DoxyCodeLine{704 \textcolor{comment}{! Local variables}} \DoxyCodeLine{705 \textcolor{keywordtype}{integer} :: i, j \textcolor{comment}{! loop indexes}} \DoxyCodeLine{706 \textcolor{keywordtype}{ real} :: h0, h1 \textcolor{comment}{! cell widths [H or nondim]}} \DoxyCodeLine{707 \textcolor{keywordtype}{ real} :: h0\_2, h1\_2, h0h1 \textcolor{comment}{! products of cell widths [H2 or nondim]}} \DoxyCodeLine{708 \textcolor{keywordtype}{ real} :: h0\_3, h1\_3 \textcolor{comment}{! products of three cell widths [H3 or nondim]}} \DoxyCodeLine{709 \textcolor{keywordtype}{ real} :: h\_min \textcolor{comment}{! A minimal cell width [H]}} \DoxyCodeLine{710 \textcolor{keywordtype}{ real} :: d \textcolor{comment}{! A temporary variable [H3]}} \DoxyCodeLine{711 \textcolor{keywordtype}{ real} :: I\_d \textcolor{comment}{! A temporary variable [nondim]}} \DoxyCodeLine{712 \textcolor{keywordtype}{ real} :: I\_h \textcolor{comment}{! Inverses of thicknesses [H-\/1]}} \DoxyCodeLine{713 \textcolor{keywordtype}{ real} :: alpha, beta \textcolor{comment}{! stencil coefficients [nondim]}} \DoxyCodeLine{714 \textcolor{keywordtype}{ real} :: a, b \textcolor{comment}{! weights of cells [H-\/1]}} \DoxyCodeLine{715 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{parameter} :: C1\_12 = 1.0 / 12.0} \DoxyCodeLine{716 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)} :: dz \textcolor{comment}{! A temporary array of limited layer thicknesses [H]}} \DoxyCodeLine{717 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)} :: u\_tmp \textcolor{comment}{! A temporary array of cell average properties [A]}} \DoxyCodeLine{718 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(5)} :: x \textcolor{comment}{! Coordinate system with 0 at edges [H]}} \DoxyCodeLine{719 \textcolor{keywordtype}{ real} :: dx, xavg \textcolor{comment}{! Differences and averages of successive values of x [H]}} \DoxyCodeLine{720 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4,4)} :: Asys \textcolor{comment}{! matrix used to find boundary conditions}} \DoxyCodeLine{721 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)} :: Bsys, Csys} \DoxyCodeLine{722 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(3)} :: Dsys} \DoxyCodeLine{723 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N+1)} :: tri\_l, \& \textcolor{comment}{! tridiagonal system (lower diagonal) [nondim]}} \DoxyCodeLine{724 tri\_d, \& \textcolor{comment}{! tridiagonal system (middle diagonal) [nondim]}} \DoxyCodeLine{725 tri\_c, \& \textcolor{comment}{! tridiagonal system central value, with tri\_d = tri\_c+tri\_l+tri\_u}} \DoxyCodeLine{726 tri\_u, \& \textcolor{comment}{! tridiagonal system (upper diagonal) [nondim]}} \DoxyCodeLine{727 tri\_b, \& \textcolor{comment}{! tridiagonal system (right hand side) [A H-\/1]}} \DoxyCodeLine{728 tri\_x \textcolor{comment}{! tridiagonal system (solution vector) [A H-\/1]}} \DoxyCodeLine{729 \textcolor{keywordtype}{ real} :: hNeglect \textcolor{comment}{! A negligible thickness [H].}} \DoxyCodeLine{730 \textcolor{keywordtype}{ real} :: hNeglect3 \textcolor{comment}{! hNeglect\string^3 [H3].}} \DoxyCodeLine{731 \textcolor{keywordtype}{logical} :: use\_2018\_answers \textcolor{comment}{! If true use older, less acccurate expressions.}} \DoxyCodeLine{732 } \DoxyCodeLine{733 hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect} \DoxyCodeLine{734 hneglect3 = hneglect**3} \DoxyCodeLine{735 use\_2018\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) use\_2018\_answers = answers\_2018} \DoxyCodeLine{736 } \DoxyCodeLine{737 \textcolor{comment}{! Loop on cells (except last one)}} \DoxyCodeLine{738 \textcolor{keywordflow}{do} i = 1,n-\/1} \DoxyCodeLine{739 } \DoxyCodeLine{740 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{741 \textcolor{comment}{! Get cell widths}} \DoxyCodeLine{742 h0 = h(i)} \DoxyCodeLine{743 h1 = h(i+1)} \DoxyCodeLine{744 } \DoxyCodeLine{745 \textcolor{comment}{! Auxiliary calculations}} \DoxyCodeLine{746 h0h1 = h0 * h1} \DoxyCodeLine{747 h0\_2 = h0 * h0} \DoxyCodeLine{748 h1\_2 = h1 * h1} \DoxyCodeLine{749 h0\_3 = h0\_2 * h0} \DoxyCodeLine{750 h1\_3 = h1\_2 * h1} \DoxyCodeLine{751 } \DoxyCodeLine{752 d = 4.0 * h0h1 * ( h0 + h1 ) + h1\_3 + h0\_3} \DoxyCodeLine{753 } \DoxyCodeLine{754 \textcolor{comment}{! Coefficients}} \DoxyCodeLine{755 alpha = h1 * (h0\_2 + h0h1 -\/ h1\_2) / ( d + hneglect3 )} \DoxyCodeLine{756 beta = h0 * (h1\_2 + h0h1 -\/ h0\_2) / ( d + hneglect3 )} \DoxyCodeLine{757 a = -\/12.0 * h0h1 / ( d + hneglect3 )} \DoxyCodeLine{758 b = -\/a} \DoxyCodeLine{759 } \DoxyCodeLine{760 tri\_l(i+1) = alpha} \DoxyCodeLine{761 tri\_d(i+1) = 1.0} \DoxyCodeLine{762 tri\_u(i+1) = beta} \DoxyCodeLine{763 } \DoxyCodeLine{764 tri\_b(i+1) = a * u(i) + b * u(i+1)} \DoxyCodeLine{765 \textcolor{keywordflow}{else}} \DoxyCodeLine{766 \textcolor{comment}{! Get cell widths}} \DoxyCodeLine{767 h0 = max(h(i), hneglect)} \DoxyCodeLine{768 h1 = max(h(i+1), hneglect)} \DoxyCodeLine{769 } \DoxyCodeLine{770 i\_h = 1.0 / (h0 + h1)} \DoxyCodeLine{771 h0 = h0 * i\_h ; h1 = h1 * i\_h} \DoxyCodeLine{772 } \DoxyCodeLine{773 h0h1 = h0 * h1 ; h0\_2 = h0 * h0 ; h1\_2 = h1 * h1} \DoxyCodeLine{774 h0\_3 = h0\_2 * h0 ; h1\_3 = h1\_2 * h1} \DoxyCodeLine{775 } \DoxyCodeLine{776 \textcolor{comment}{! Set the tridiagonal coefficients}} \DoxyCodeLine{777 i\_d = 1.0 / (4.0 * h0h1 * ( h0 + h1 ) + h1\_3 + h0\_3) \textcolor{comment}{! = 1 / ((h0 + h1)**3 + h0*h1*(h0 + h1))}} \DoxyCodeLine{778 tri\_l(i+1) = (h1 * ((h0\_2 + h0h1) -\/ h1\_2)) * i\_d} \DoxyCodeLine{779 \textcolor{comment}{! tri\_d(i+1) = 1.0}} \DoxyCodeLine{780 tri\_c(i+1) = 2.0 * ((h0\_2 + h1\_2) * (h0 + h1)) * i\_d} \DoxyCodeLine{781 tri\_u(i+1) = (h0 * ((h1\_2 + h0h1) -\/ h0\_2)) * i\_d} \DoxyCodeLine{782 \textcolor{comment}{! The following expressions have been simplified using the nondimensionalization above:}} \DoxyCodeLine{783 \textcolor{comment}{! I\_d = 1.0 / (1.0 + h0h1)}} \DoxyCodeLine{784 \textcolor{comment}{! tri\_l(i+1) = (h0h1 -\/ h1\_3) * I\_d}} \DoxyCodeLine{785 \textcolor{comment}{! tri\_c(i+1) = 2.0 * (h0\_2 + h1\_2) * I\_d}} \DoxyCodeLine{786 \textcolor{comment}{! tri\_u(i+1) = (h0h1 -\/ h0\_3) * I\_d}} \DoxyCodeLine{787 } \DoxyCodeLine{788 tri\_b(i+1) = 12.0 * (h0h1 * i\_d) * ((u(i+1) -\/ u(i)) * i\_h)} \DoxyCodeLine{789 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{790 } \DoxyCodeLine{791 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on cells}} \DoxyCodeLine{792 } \DoxyCodeLine{793 \textcolor{comment}{! Boundary conditions: set the first edge slope}} \DoxyCodeLine{794 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{795 x(1) = 0.0} \DoxyCodeLine{796 \textcolor{keywordflow}{do} i = 1,4} \DoxyCodeLine{797 dx = h(i)} \DoxyCodeLine{798 x(i+1) = x(i) + dx} \DoxyCodeLine{799 \textcolor{keywordflow}{do} j = 1,4 ; asys(i,j) = ( (x(i+1)**j) -\/ (x(i)**j) ) / j ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{800 bsys(i) = u(i) * dx} \DoxyCodeLine{801 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{802 } \DoxyCodeLine{803 \textcolor{keyword}{call }solve\_linear\_system( asys, bsys, csys, 4 )} \DoxyCodeLine{804 } \DoxyCodeLine{805 dsys(1) = csys(2) ; dsys(2) = 2.0 * csys(3) ; dsys(3) = 3.0 * csys(4)} \DoxyCodeLine{806 tri\_b(1) = evaluation\_polynomial( dsys, 3, x(1) ) \textcolor{comment}{! Set the first edge slope}} \DoxyCodeLine{807 tri\_d(1) = 1.0} \DoxyCodeLine{808 \textcolor{keywordflow}{else} \textcolor{comment}{! Use expressions with less sensitivity to roundoff}} \DoxyCodeLine{809 \textcolor{keywordflow}{do} i=1,4 ; dz(i) = max(hneglect, h(i) ) ; u\_tmp(i) = u(i) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{810 \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, csys)} \DoxyCodeLine{811 } \DoxyCodeLine{812 \textcolor{comment}{! Set the first edge slope}} \DoxyCodeLine{813 tri\_b(1) = csys(2)} \DoxyCodeLine{814 tri\_c(1) = 1.0} \DoxyCodeLine{815 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{816 tri\_u(1) = 0.0 \textcolor{comment}{! tri\_l(1) = 0.0}} \DoxyCodeLine{817 } \DoxyCodeLine{818 \textcolor{comment}{! Boundary conditions: set the last edge slope}} \DoxyCodeLine{819 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{820 x(1) = 0.0} \DoxyCodeLine{821 \textcolor{keywordflow}{do} i = 1,4} \DoxyCodeLine{822 dx = h(n-\/4+i)} \DoxyCodeLine{823 x(i+1) = x(i) + dx} \DoxyCodeLine{824 \textcolor{keywordflow}{do} j = 1,4 ; asys(i,j) = ( (x(i+1)**j) -\/ (x(i)**j) ) / j ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{825 bsys(i) = u(n-\/4+i) * dx} \DoxyCodeLine{826 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{827 } \DoxyCodeLine{828 \textcolor{keyword}{call }solve\_linear\_system( asys, bsys, csys, 4 )} \DoxyCodeLine{829 } \DoxyCodeLine{830 dsys(1) = csys(2) ; dsys(2) = 2.0 * csys(3) ; dsys(3) = 3.0 * csys(4)} \DoxyCodeLine{831 \textcolor{comment}{! Set the last edge slope}} \DoxyCodeLine{832 tri\_b(n+1) = evaluation\_polynomial( dsys, 3, x(5) )} \DoxyCodeLine{833 tri\_d(n+1) = 1.0} \DoxyCodeLine{834 \textcolor{keywordflow}{else}} \DoxyCodeLine{835 \textcolor{comment}{! Use expressions with less sensitivity to roundoff, including using a coordinate}} \DoxyCodeLine{836 \textcolor{comment}{! system that sets the origin at the last interface in the domain.}} \DoxyCodeLine{837 \textcolor{keywordflow}{do} i=1,4 ; dz(i) = max(hneglect, h(n+1-\/i) ) ; u\_tmp(i) = u(n+1-\/i) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{838 } \DoxyCodeLine{839 \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, csys)} \DoxyCodeLine{840 } \DoxyCodeLine{841 \textcolor{comment}{! Set the last edge slope}} \DoxyCodeLine{842 tri\_b(n+1) = -\/csys(2)} \DoxyCodeLine{843 tri\_c(n+1) = 1.0} \DoxyCodeLine{844 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{845 tri\_l(n+1) = 0.0 \textcolor{comment}{! tri\_u(N+1) = 0.0}} \DoxyCodeLine{846 } \DoxyCodeLine{847 \textcolor{comment}{! Solve tridiagonal system and assign edge slopes}} \DoxyCodeLine{848 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{849 \textcolor{keyword}{call }solve\_tridiagonal\_system( tri\_l, tri\_d, tri\_u, tri\_b, tri\_x, n+1 )} \DoxyCodeLine{850 \textcolor{keywordflow}{else}} \DoxyCodeLine{851 \textcolor{keyword}{call }solve\_diag\_dominant\_tridiag( tri\_l, tri\_c, tri\_u, tri\_b, tri\_x, n+1 )} \DoxyCodeLine{852 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{853 } \DoxyCodeLine{854 \textcolor{keywordflow}{do} i = 2,n} \DoxyCodeLine{855 edge\_slopes(i,1) = tri\_x(i)} \DoxyCodeLine{856 edge\_slopes(i-\/1,2) = tri\_x(i)} \DoxyCodeLine{857 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{858 edge\_slopes(1,1) = tri\_x(1)} \DoxyCodeLine{859 edge\_slopes(n,2) = tri\_x(n+1)} \DoxyCodeLine{860 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a7c3fb27d99df9e3594ff32d1d03e1b34}\label{namespaceregrid__edge__values_a7c3fb27d99df9e3594ff32d1d03e1b34}} \index{regrid\_edge\_values@{regrid\_edge\_values}!edge\_slopes\_implicit\_h5@{edge\_slopes\_implicit\_h5}} \index{edge\_slopes\_implicit\_h5@{edge\_slopes\_implicit\_h5}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{edge\_slopes\_implicit\_h5()}{edge\_slopes\_implicit\_h5()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+edge\+\_\+values\+::edge\+\_\+slopes\+\_\+implicit\+\_\+h5 (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n), intent(in)}]{h, }\item[{real, dimension(n), intent(in)}]{u, }\item[{real, dimension(n,2), intent(inout)}]{edge\+\_\+slopes, }\item[{real, intent(in), optional}]{h\+\_\+neglect, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})} Compute ih5 edge slopes (implicit fifth order accurate) \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & Number of cells \\ \hline \mbox{\texttt{ in}} & {\em h} & cell widths \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em edge\+\_\+slopes} & Returned edge slopes \mbox{[}A H-\/1\mbox{]}; the second index is for the two edges of each cell. \\ \hline \mbox{\texttt{ in}} & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\ \hline \end{DoxyParams} Definition at line 867 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{867 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells}} \DoxyCodeLine{868 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths [H]}} \DoxyCodeLine{869 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}} \DoxyCodeLine{870 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_slopes\textcolor{comment}{ !< Returned edge slopes [A H-\/1]; the}} \DoxyCodeLine{871 \textcolor{comment}{ !! second index is for the two edges of each cell.}} \DoxyCodeLine{872 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}} \DoxyCodeLine{873 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}} \DoxyCodeLine{874 \textcolor{comment}{! -\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/}} \DoxyCodeLine{875 \textcolor{comment}{! Fifth-\/order implicit estimates of edge slopes are based on a four-\/cell,}} \DoxyCodeLine{876 \textcolor{comment}{! three-\/edge stencil. A tridiagonal system is set up and is based on}} \DoxyCodeLine{877 \textcolor{comment}{! expressing the edge slopes in terms of neighboring cell averages.}} \DoxyCodeLine{878 \textcolor{comment}{!}} \DoxyCodeLine{879 \textcolor{comment}{! The generic relationship is}} \DoxyCodeLine{880 \textcolor{comment}{!}} \DoxyCodeLine{881 \textcolor{comment}{! \(\backslash\)alpha u'\_\{i-\/1/2\} + u'\_\{i+1/2\} + \(\backslash\)beta u'\_\{i+3/2\} =}} \DoxyCodeLine{882 \textcolor{comment}{! a \(\backslash\)bar\{u\}\_\{i-\/1\} + b \(\backslash\)bar\{u\}\_i + c \(\backslash\)bar\{u\}\_\{i+1\} + d \(\backslash\)bar\{u\}\_\{i+2\}}} \DoxyCodeLine{883 \textcolor{comment}{!}} \DoxyCodeLine{884 \textcolor{comment}{! and the stencil looks like this}} \DoxyCodeLine{885 \textcolor{comment}{!}} \DoxyCodeLine{886 \textcolor{comment}{! i-\/1 i i+1 i+2}} \DoxyCodeLine{887 \textcolor{comment}{! ..-\/-\/o-\/-\/-\/-\/-\/-\/o-\/-\/-\/-\/-\/-\/o-\/-\/-\/-\/-\/-\/o-\/-\/-\/-\/-\/-\/o-\/-\/..}} \DoxyCodeLine{888 \textcolor{comment}{! i-\/1/2 i+1/2 i+3/2}} \DoxyCodeLine{889 \textcolor{comment}{!}} \DoxyCodeLine{890 \textcolor{comment}{! In this routine, the coefficients \(\backslash\)alpha, \(\backslash\)beta, a, b, c and d are}} \DoxyCodeLine{891 \textcolor{comment}{! computed, the tridiagonal system is built, boundary conditions are}} \DoxyCodeLine{892 \textcolor{comment}{! prescribed and the system is solved to yield edge-\/value estimates.}} \DoxyCodeLine{893 \textcolor{comment}{!}} \DoxyCodeLine{894 \textcolor{comment}{! Note that the centered stencil only applies to edges 3 to N-\/1 (edges are}} \DoxyCodeLine{895 \textcolor{comment}{! numbered 1 to n+1), which yields N-\/3 equations for N+1 unknowns. Two other}} \DoxyCodeLine{896 \textcolor{comment}{! equations are written by using a right-\/biased stencil for edge 2 and a}} \DoxyCodeLine{897 \textcolor{comment}{! left-\/biased stencil for edge N. The prescription of boundary conditions}} \DoxyCodeLine{898 \textcolor{comment}{! (using sixth-\/order polynomials) closes the system.}} \DoxyCodeLine{899 \textcolor{comment}{!}} \DoxyCodeLine{900 \textcolor{comment}{! CAUTION: For each edge, in order to determine the coefficients of the}} \DoxyCodeLine{901 \textcolor{comment}{! implicit expression, a 6x6 linear system is solved. This may}} \DoxyCodeLine{902 \textcolor{comment}{! become computationally expensive if regridding is carried out}} \DoxyCodeLine{903 \textcolor{comment}{! often. Figuring out closed-\/form expressions for these coefficients}} \DoxyCodeLine{904 \textcolor{comment}{! on nonuniform meshes turned out to be intractable.}} \DoxyCodeLine{905 \textcolor{comment}{! -\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/}} \DoxyCodeLine{906 } \DoxyCodeLine{907 \textcolor{comment}{! Local variables}} \DoxyCodeLine{908 \textcolor{keywordtype}{ real} :: h0, h1, h2, h3 \textcolor{comment}{! cell widths [H]}} \DoxyCodeLine{909 \textcolor{keywordtype}{ real} :: hMin \textcolor{comment}{! The minimum thickness used in these calculations [H]}} \DoxyCodeLine{910 \textcolor{keywordtype}{ real} :: h01, h01\_2 \textcolor{comment}{! Summed thicknesses to various powers [H\string^n \string~> m\string^n or kg\string^n m-\/2n]}} \DoxyCodeLine{911 \textcolor{keywordtype}{ real} :: h23, h23\_2 \textcolor{comment}{! Summed thicknesses to various powers [H\string^n \string~> m\string^n or kg\string^n m-\/2n]}} \DoxyCodeLine{912 \textcolor{keywordtype}{ real} :: hNeglect \textcolor{comment}{! A negligible thickness [H].}} \DoxyCodeLine{913 \textcolor{keywordtype}{ real} :: h1\_2, h2\_2 \textcolor{comment}{! the coefficients of the}} \DoxyCodeLine{914 \textcolor{keywordtype}{ real} :: h1\_3, h2\_3 \textcolor{comment}{! tridiagonal system}} \DoxyCodeLine{915 \textcolor{keywordtype}{ real} :: h1\_4, h2\_4 \textcolor{comment}{! ...}} \DoxyCodeLine{916 \textcolor{keywordtype}{ real} :: h1\_5, h2\_5 \textcolor{comment}{! ...}} \DoxyCodeLine{917 \textcolor{keywordtype}{ real} :: alpha, beta \textcolor{comment}{! stencil coefficients}} \DoxyCodeLine{918 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(7)} :: x \textcolor{comment}{! Coordinate system with 0 at edges [same units as h]}} \DoxyCodeLine{919 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{parameter} :: C1\_12 = 1.0 / 12.0} \DoxyCodeLine{920 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{parameter} :: C5\_6 = 5.0 / 6.0} \DoxyCodeLine{921 \textcolor{keywordtype}{ real} :: dx, xavg \textcolor{comment}{! Differences and averages of successive values of x [same units as h]}} \DoxyCodeLine{922 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(6,6)} :: Asys \textcolor{comment}{! matrix used to find boundary conditions}} \DoxyCodeLine{923 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(6)} :: Bsys, Csys \textcolor{comment}{! ...}} \DoxyCodeLine{924 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N+1)} :: tri\_l, \& \textcolor{comment}{! trid. system (lower diagonal)}} \DoxyCodeLine{925 tri\_d, \& \textcolor{comment}{! trid. system (middle diagonal)}} \DoxyCodeLine{926 tri\_u, \& \textcolor{comment}{! trid. system (upper diagonal)}} \DoxyCodeLine{927 tri\_b, \& \textcolor{comment}{! trid. system (unknowns vector)}} \DoxyCodeLine{928 tri\_x \textcolor{comment}{! trid. system (rhs)}} \DoxyCodeLine{929 \textcolor{keywordtype}{ real} :: h\_Min\_Frac = 1.0e-\/4} \DoxyCodeLine{930 \textcolor{keywordtype}{integer} :: i, j, k \textcolor{comment}{! loop indexes}} \DoxyCodeLine{931 } \DoxyCodeLine{932 hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect} \DoxyCodeLine{933 } \DoxyCodeLine{934 \textcolor{comment}{! Loop on cells (except the first and last ones)}} \DoxyCodeLine{935 \textcolor{keywordflow}{do} k = 2,n-\/2} \DoxyCodeLine{936 \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}} \DoxyCodeLine{937 hmin = max(hneglect, h\_min\_frac*((h(k-\/1) + h(k)) + (h(k+1) + h(k+2))))} \DoxyCodeLine{938 h0 = max(h(k-\/1), hmin) ; h1 = max(h(k), hmin)} \DoxyCodeLine{939 h2 = max(h(k+1), hmin) ; h3 = max(h(k+2), hmin)} \DoxyCodeLine{940 } \DoxyCodeLine{941 \textcolor{comment}{! Auxiliary calculations}} \DoxyCodeLine{942 h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2} \DoxyCodeLine{943 h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2} \DoxyCodeLine{944 } \DoxyCodeLine{945 \textcolor{comment}{! Compute matrix entries as described in Eq. (52) of White and Adcroft (2009). The last 4 rows are}} \DoxyCodeLine{946 \textcolor{comment}{! Asys(1:6,n) = (/ -\/n*(n-\/1)*(-\/h1)**(n-\/2), -\/n*(n-\/1)*h1**(n-\/2), (-\/1)**(n-\/1) * ((h0+h1)**n -\/ h0**n) / h0, \&}} \DoxyCodeLine{947 \textcolor{comment}{! (-\/h1)**(n-\/1), h2**(n-\/1), ((h2+h3)**n -\/ h2**n) / h3 /)}} \DoxyCodeLine{948 } \DoxyCodeLine{949 asys(1:6,1) = (/ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /)} \DoxyCodeLine{950 asys(1:6,2) = (/ 2.0, 2.0, (2.0*h1 + h0), h1, -\/h2, -\/(2.0*h2 + h3) /)} \DoxyCodeLine{951 asys(1:6,3) = (/ 6.0*h1, -\/6.0* h2, (3.0*h1\_2 + h0*(3.0*h1 + h0)), \&} \DoxyCodeLine{952 h1\_2, h2\_2, (3.0*h2\_2 + h3*(3.0*h2 + h3)) /)} \DoxyCodeLine{953 asys(1:6,4) = (/ -\/12.0*h1\_2, -\/12.0*h2\_2, -\/(4.0*h1\_3 + h0*(6.0*h1\_2 + h0*(4.0*h1 + h0))), \&} \DoxyCodeLine{954 -\/h1\_3, h2\_3, (4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)} \DoxyCodeLine{955 asys(1:6,5) = (/ 20.0*h1\_3, -\/20.0*h2\_3, (5.0*h1\_4 + h0*(10.0*h1\_3 + h0*(10.0*h1\_2 + h0*(5.0*h1 + h0)))), \&} \DoxyCodeLine{956 h1\_4, h2\_4, (5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))) /)} \DoxyCodeLine{957 asys(1:6,6) = (/ -\/30.0*h1\_4, -\/30.0*h2\_4, \&} \DoxyCodeLine{958 -\/(6.0*h1\_5 + h0*(15.0*h1\_4 + h0*(20.0*h1\_3 + h0*(15.0*h1\_2 + h0*(6.0*h1 + h0))))), \&} \DoxyCodeLine{959 -\/h1\_5, h2\_5, \&} \DoxyCodeLine{960 (6.0*h2\_5 + h3*(15.0*h2\_4 + h3*(20.0*h2\_3 + h3*(15.0*h2\_2 + h3*(6.0*h2 + h3))))) /)} \DoxyCodeLine{961 bsys(1:6) = (/ 0.0, -\/2.0, 0.0, 0.0, 0.0, 0.0 /)} \DoxyCodeLine{962 } \DoxyCodeLine{963 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{964 } \DoxyCodeLine{965 alpha = csys(1)} \DoxyCodeLine{966 beta = csys(2)} \DoxyCodeLine{967 } \DoxyCodeLine{968 tri\_l(k+1) = alpha} \DoxyCodeLine{969 tri\_d(k+1) = 1.0} \DoxyCodeLine{970 tri\_u(k+1) = beta} \DoxyCodeLine{971 tri\_b(k+1) = csys(3) * u(k-\/1) + csys(4) * u(k) + csys(5) * u(k+1) + csys(6) * u(k+2)} \DoxyCodeLine{972 } \DoxyCodeLine{973 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on cells}} \DoxyCodeLine{974 } \DoxyCodeLine{975 \textcolor{comment}{! Use a right-\/biased stencil for the second row, as described in Eq. (53) of White and Adcroft (2009).}} \DoxyCodeLine{976 } \DoxyCodeLine{977 \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}} \DoxyCodeLine{978 hmin = max(hneglect, h\_min\_frac*((h(1) + h(2)) + (h(3) + h(4))))} \DoxyCodeLine{979 h0 = max(h(1), hmin) ; h1 = max(h(2), hmin)} \DoxyCodeLine{980 h2 = max(h(3), hmin) ; h3 = max(h(4), hmin)} \DoxyCodeLine{981 } \DoxyCodeLine{982 \textcolor{comment}{! Auxiliary calculations}} \DoxyCodeLine{983 h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2} \DoxyCodeLine{984 h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2} \DoxyCodeLine{985 h01 = h0 + h1 ; h01\_2 = h01 * h01} \DoxyCodeLine{986 } \DoxyCodeLine{987 \textcolor{comment}{! Compute matrix entries}} \DoxyCodeLine{988 asys(1:6,1) = (/ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /)} \DoxyCodeLine{989 asys(1:6,2) = (/ 2.0, 2.0, (2.0*h1 + h0), h1, -\/h2, -\/(2.0*h2 + h3) /)} \DoxyCodeLine{990 asys(1:6,3) = (/ 6.0*h01, 0.0, (3.0*h1\_2 + h0*(3.0*h1 + h0)), \&} \DoxyCodeLine{991 h1\_2, h2\_2, (3.0*h2\_2 + h3*(3.0*h2 + h3)) /)} \DoxyCodeLine{992 asys(1:6,4) = (/ -\/12.0*h01\_2, 0.0, -\/(4.0*h1\_3 + h0*(6.0*h1\_2 + h0*(4.0*h1 + h0))), \&} \DoxyCodeLine{993 -\/h1\_3, h2\_3, (4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)} \DoxyCodeLine{994 asys(1:6,5) = (/ 20.0*(h01*h01\_2), 0.0, (5.0*h1\_4 + h0*(10.0*h1\_3 + h0*(10.0*h1\_2 + h0*(5.0*h1 + h0)))), \&} \DoxyCodeLine{995 h1\_4, h2\_4, (5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))) /)} \DoxyCodeLine{996 asys(1:6,6) = (/ -\/30.0*(h01\_2*h01\_2), 0.0, \&} \DoxyCodeLine{997 -\/(6.0*h1\_5 + h0*(15.0*h1\_4 + h0*(20.0*h1\_3 + h0*(15.0*h1\_2 + h0*(6.0*h1 + h0))))), \&} \DoxyCodeLine{998 -\/h1\_5, h2\_5, \&} \DoxyCodeLine{999 (6.0*h2\_5 + h3*(15.0*h2\_4 + h3*(20.0*h2\_3 + h3*(15.0*h2\_2 + h3*(6.0*h2 + h3))))) /)} \DoxyCodeLine{1000 bsys(1:6) = (/ 0.0, -\/2.0, -\/6.0*h1, 12.0*h1\_2, -\/20.0*h1\_3, 30.0*h1\_4 /)} \DoxyCodeLine{1001 } \DoxyCodeLine{1002 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{1003 } \DoxyCodeLine{1004 alpha = csys(1)} \DoxyCodeLine{1005 beta = csys(2)} \DoxyCodeLine{1006 } \DoxyCodeLine{1007 tri\_l(2) = alpha} \DoxyCodeLine{1008 tri\_d(2) = 1.0} \DoxyCodeLine{1009 tri\_u(2) = beta} \DoxyCodeLine{1010 tri\_b(2) = csys(3) * u(1) + csys(4) * u(2) + csys(5) * u(3) + csys(6) * u(4)} \DoxyCodeLine{1011 } \DoxyCodeLine{1012 \textcolor{comment}{! Boundary conditions: left boundary}} \DoxyCodeLine{1013 x(1) = 0.0} \DoxyCodeLine{1014 \textcolor{keywordflow}{do} i = 1,6} \DoxyCodeLine{1015 dx = h(i)} \DoxyCodeLine{1016 xavg = x(i) + 0.5 * dx} \DoxyCodeLine{1017 asys(1:6,i) = (/ 1.0, xavg, (xavg**2 + c1\_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), \&} \DoxyCodeLine{1018 (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), \&} \DoxyCodeLine{1019 xavg * (xavg**4 + c5\_6*xavg**2*dx**2 + 0.0625*dx**4) /)} \DoxyCodeLine{1020 bsys(i) = u(i)} \DoxyCodeLine{1021 x(i+1) = x(i) + dx} \DoxyCodeLine{1022 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1023 } \DoxyCodeLine{1024 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{1025 } \DoxyCodeLine{1026 tri\_d(1) = 0.0} \DoxyCodeLine{1027 tri\_d(1) = 1.0} \DoxyCodeLine{1028 tri\_u(1) = 0.0} \DoxyCodeLine{1029 tri\_b(1) = csys(2) \textcolor{comment}{! first edge value}} \DoxyCodeLine{1030 } \DoxyCodeLine{1031 \textcolor{comment}{! Use a left-\/biased stencil for the second to last row, as described in Eq. (54) of White and Adcroft (2009).}} \DoxyCodeLine{1032 } \DoxyCodeLine{1033 \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}} \DoxyCodeLine{1034 hmin = max(hneglect, h\_min\_frac*((h(n-\/3) + h(n-\/2)) + (h(n-\/1) + h(n))))} \DoxyCodeLine{1035 h0 = max(h(n-\/3), hmin) ; h1 = max(h(n-\/2), hmin)} \DoxyCodeLine{1036 h2 = max(h(n-\/1), hmin) ; h3 = max(h(n), hmin)} \DoxyCodeLine{1037 } \DoxyCodeLine{1038 \textcolor{comment}{! Auxiliary calculations}} \DoxyCodeLine{1039 h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2} \DoxyCodeLine{1040 h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2} \DoxyCodeLine{1041 } \DoxyCodeLine{1042 h23 = h2 + h3 ; h23\_2 = h23 * h23} \DoxyCodeLine{1043 } \DoxyCodeLine{1044 \textcolor{comment}{! Compute matrix entries}} \DoxyCodeLine{1045 asys(1:6,1) = (/ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /)} \DoxyCodeLine{1046 asys(1:6,2) = (/ 2.0, 2.0, (2.0*h1 + h0), h1, -\/h2, -\/(2.0*h2 + h3) /)} \DoxyCodeLine{1047 asys(1:6,3) = (/ 0.0, -\/6.0*h23, (3.0*h1\_2 + h0*(3.0*h1 + h0)), \&} \DoxyCodeLine{1048 h1\_2, h2\_2, (3.0*h2\_2 + h3*(3.0*h2 + h3)) /)} \DoxyCodeLine{1049 asys(1:6,4) = (/ 0.0, -\/12.0*h23\_2, -\/(4.0*h1\_3 + h0*(6.0*h1\_2 + h0*(4.0*h1 + h0))), \&} \DoxyCodeLine{1050 -\/h1\_3, h2\_3, (4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)} \DoxyCodeLine{1051 asys(1:6,5) = (/ 0.0, -\/20.0*(h23*h23\_2), (5.0*h1\_4 + h0*(10.0*h1\_3 + h0*(10.0*h1\_2 + h0*(5.0*h1 + h0)))), \&} \DoxyCodeLine{1052 h1\_4, h2\_4, (5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))) /)} \DoxyCodeLine{1053 asys(1:6,6) = (/ 0.0, -\/30.0*(h23\_2*h23\_2), \&} \DoxyCodeLine{1054 -\/(6.0*h1\_5 + h0*(15.0*h1\_4 + h0*(20.0*h1\_3 + h0*(15.0*h1\_2 + h0*(6.0*h1 + h0))))), \&} \DoxyCodeLine{1055 -\/h1\_5, h2\_5, \&} \DoxyCodeLine{1056 (6.0*h2\_5 + h3*(15.0*h2\_4 + h3*(20.0*h2\_3 + h3*(15.0*h2\_2 + h3*(6.0*h2 + h3))))) /)} \DoxyCodeLine{1057 bsys(1:6) = (/ 0.0, -\/2.0, 6.0*h2, 12.0*h2\_2, 20.0*h2\_3, 30.0*h2\_4 /)} \DoxyCodeLine{1058 } \DoxyCodeLine{1059 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{1060 } \DoxyCodeLine{1061 alpha = csys(1)} \DoxyCodeLine{1062 beta = csys(2)} \DoxyCodeLine{1063 } \DoxyCodeLine{1064 tri\_l(n) = alpha} \DoxyCodeLine{1065 tri\_d(n) = 1.0} \DoxyCodeLine{1066 tri\_u(n) = beta} \DoxyCodeLine{1067 tri\_b(n) = csys(3) * u(n-\/3) + csys(4) * u(n-\/2) + csys(5) * u(n-\/1) + csys(6) * u(n)} \DoxyCodeLine{1068 } \DoxyCodeLine{1069 \textcolor{comment}{! Boundary conditions: right boundary}} \DoxyCodeLine{1070 x(1) = 0.0} \DoxyCodeLine{1071 \textcolor{keywordflow}{do} i = 1,6} \DoxyCodeLine{1072 dx = h(n+1-\/i)} \DoxyCodeLine{1073 xavg = x(i) + 0.5*dx} \DoxyCodeLine{1074 asys(1:6,i) = (/ 1.0, xavg, (xavg**2 + c1\_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), \&} \DoxyCodeLine{1075 (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), \&} \DoxyCodeLine{1076 xavg * (xavg**4 + c5\_6*xavg**2*dx**2 + 0.0625*dx**4) /)} \DoxyCodeLine{1077 bsys(i) = u(n+1-\/i)} \DoxyCodeLine{1078 x(i+1) = x(i) + dx} \DoxyCodeLine{1079 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1080 } \DoxyCodeLine{1081 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{1082 } \DoxyCodeLine{1083 tri\_l(n+1) = 0.0} \DoxyCodeLine{1084 tri\_d(n+1) = 1.0} \DoxyCodeLine{1085 tri\_u(n+1) = 0.0} \DoxyCodeLine{1086 tri\_b(n+1) = -\/csys(2)} \DoxyCodeLine{1087 } \DoxyCodeLine{1088 \textcolor{comment}{! Solve tridiagonal system and assign edge values}} \DoxyCodeLine{1089 \textcolor{keyword}{call }solve\_tridiagonal\_system( tri\_l, tri\_d, tri\_u, tri\_b, tri\_x, n+1 )} \DoxyCodeLine{1090 } \DoxyCodeLine{1091 \textcolor{keywordflow}{do} i = 2,n} \DoxyCodeLine{1092 edge\_slopes(i,1) = tri\_x(i)} \DoxyCodeLine{1093 edge\_slopes(i-\/1,2) = tri\_x(i)} \DoxyCodeLine{1094 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1095 edge\_slopes(1,1) = tri\_x(1)} \DoxyCodeLine{1096 edge\_slopes(n,2) = tri\_x(n+1)} \DoxyCodeLine{1097 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a8dc053bb28658c7091f1f5a82184019a}\label{namespaceregrid__edge__values_a8dc053bb28658c7091f1f5a82184019a}} \index{regrid\_edge\_values@{regrid\_edge\_values}!edge\_values\_explicit\_h2@{edge\_values\_explicit\_h2}} \index{edge\_values\_explicit\_h2@{edge\_values\_explicit\_h2}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{edge\_values\_explicit\_h2()}{edge\_values\_explicit\_h2()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+edge\+\_\+values\+::edge\+\_\+values\+\_\+explicit\+\_\+h2 (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n), intent(in)}]{h, }\item[{real, dimension(n), intent(in)}]{u, }\item[{real, dimension(n,2), intent(inout)}]{edge\+\_\+val }\end{DoxyParamCaption})} Compute h2 edge values (explicit second order accurate) in the same units as u. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & Number of cells \\ \hline \mbox{\texttt{ in}} & {\em h} & cell widths \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em edge\+\_\+val} & Returned edge values \mbox{[}A\mbox{]}; the second index is for the two edges of each cell. \\ \hline \end{DoxyParams} Definition at line 175 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{175 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells}} \DoxyCodeLine{176 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths [H]}} \DoxyCodeLine{177 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}} \DoxyCodeLine{178 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Returned edge values [A]; the}} \DoxyCodeLine{179 \textcolor{comment}{ !! second index is for the two edges of each cell.}} \DoxyCodeLine{180 } \DoxyCodeLine{181 \textcolor{comment}{! Local variables}} \DoxyCodeLine{182 \textcolor{keywordtype}{integer} :: k \textcolor{comment}{! loop index}} \DoxyCodeLine{183 } \DoxyCodeLine{184 \textcolor{comment}{! Boundary edge values are simply equal to the boundary cell averages}} \DoxyCodeLine{185 edge\_val(1,1) = u(1)} \DoxyCodeLine{186 edge\_val(n,2) = u(n)} \DoxyCodeLine{187 } \DoxyCodeLine{188 \textcolor{keywordflow}{do} k = 2,n} \DoxyCodeLine{189 \textcolor{comment}{! Compute left edge value}} \DoxyCodeLine{190 \textcolor{keywordflow}{if} (h(k-\/1) + h(k) == 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! Avoid singularities when h0+h1=0}} \DoxyCodeLine{191 edge\_val(k,1) = 0.5 * (u(k-\/1) + u(k))} \DoxyCodeLine{192 \textcolor{keywordflow}{else}} \DoxyCodeLine{193 edge\_val(k,1) = ( u(k-\/1)*h(k) + u(k)*h(k-\/1) ) / ( h(k-\/1) + h(k) )} \DoxyCodeLine{194 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{195 } \DoxyCodeLine{196 \textcolor{comment}{! Left edge value of the current cell is equal to right edge value of left cell}} \DoxyCodeLine{197 edge\_val(k-\/1,2) = edge\_val(k,1)} \DoxyCodeLine{198 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{199 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a2fa3a70b208abe68f29bef9589081209}\label{namespaceregrid__edge__values_a2fa3a70b208abe68f29bef9589081209}} \index{regrid\_edge\_values@{regrid\_edge\_values}!edge\_values\_explicit\_h4@{edge\_values\_explicit\_h4}} \index{edge\_values\_explicit\_h4@{edge\_values\_explicit\_h4}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{edge\_values\_explicit\_h4()}{edge\_values\_explicit\_h4()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+edge\+\_\+values\+::edge\+\_\+values\+\_\+explicit\+\_\+h4 (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n), intent(in)}]{h, }\item[{real, dimension(n), intent(in)}]{u, }\item[{real, dimension(n,2), intent(inout)}]{edge\+\_\+val, }\item[{real, intent(in), optional}]{h\+\_\+neglect, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})} Compute h4 edge values (explicit fourth order accurate) in the same units as u. Compute edge values based on fourth-\/order explicit estimates. These estimates are based on a cubic interpolant spanning four cells and evaluated at the location of the middle edge. An interpolant spanning cells i-\/2, i-\/1, i and i+1 is evaluated at edge i-\/1/2. The estimate for each edge is unique. \begin{DoxyVerb} i-2 i-1 i i+1 \end{DoxyVerb} ..--o-\/-\/-\/---o-\/-\/-\/---o-\/-\/-\/---o-\/-\/-\/---o--.. i-\/1/2 The first two edge values are estimated by evaluating the first available cubic interpolant, i.\+e., the interpolant spanning cells 1, 2, 3 and 4. Similarly, the last two edge values are estimated by evaluating the last available interpolant. For this fourth-\/order scheme, at least four cells must exist. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & Number of cells \\ \hline \mbox{\texttt{ in}} & {\em h} & cell widths \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em edge\+\_\+val} & Returned edge values \mbox{[}A\mbox{]}; the second index is for the two edges of each cell. \\ \hline \mbox{\texttt{ in}} & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\ \hline \end{DoxyParams} Definition at line 222 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{222 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells}} \DoxyCodeLine{223 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths [H]}} \DoxyCodeLine{224 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}} \DoxyCodeLine{225 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Returned edge values [A]; the second index}} \DoxyCodeLine{226 \textcolor{comment}{ !! is for the two edges of each cell.}} \DoxyCodeLine{227 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}} \DoxyCodeLine{228 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}} \DoxyCodeLine{229 } \DoxyCodeLine{230 \textcolor{comment}{! Local variables}} \DoxyCodeLine{231 \textcolor{keywordtype}{ real} :: h0, h1, h2, h3 \textcolor{comment}{! temporary thicknesses [H]}} \DoxyCodeLine{232 \textcolor{keywordtype}{ real} :: h\_sum \textcolor{comment}{! A sum of adjacent thicknesses [H]}} \DoxyCodeLine{233 \textcolor{keywordtype}{ real} :: h\_min \textcolor{comment}{! A minimal cell width [H]}} \DoxyCodeLine{234 \textcolor{keywordtype}{ real} :: f1, f2, f3 \textcolor{comment}{! auxiliary variables with various units}} \DoxyCodeLine{235 \textcolor{keywordtype}{ real} :: et1, et2, et3 \textcolor{comment}{! terms the expresson for edge values [A H]}} \DoxyCodeLine{236 \textcolor{keywordtype}{ real} :: I\_h12 \textcolor{comment}{! The inverse of the sum of the two central thicknesses [H-\/1]}} \DoxyCodeLine{237 \textcolor{keywordtype}{ real} :: I\_h012, I\_h123 \textcolor{comment}{! Inverses of sums of three succesive thicknesses [H-\/1]}} \DoxyCodeLine{238 \textcolor{keywordtype}{ real} :: I\_den\_et2, I\_den\_et3 \textcolor{comment}{! Inverses of denominators in edge value terms [H-\/2]}} \DoxyCodeLine{239 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(5)} :: x \textcolor{comment}{! Coordinate system with 0 at edges [H]}} \DoxyCodeLine{240 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)} :: dz \textcolor{comment}{! A temporary array of limited layer thicknesses [H]}} \DoxyCodeLine{241 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)} :: u\_tmp \textcolor{comment}{! A temporary array of cell average properties [A]}} \DoxyCodeLine{242 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{parameter} :: C1\_12 = 1.0 / 12.0} \DoxyCodeLine{243 \textcolor{keywordtype}{ real} :: dx, xavg \textcolor{comment}{! Differences and averages of successive values of x [H]}} \DoxyCodeLine{244 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4,4)} :: A \textcolor{comment}{! values near the boundaries}} \DoxyCodeLine{245 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)} :: B, C} \DoxyCodeLine{246 \textcolor{keywordtype}{ real} :: hNeglect \textcolor{comment}{! A negligible thickness in the same units as h.}} \DoxyCodeLine{247 \textcolor{keywordtype}{integer} :: i, j} \DoxyCodeLine{248 \textcolor{keywordtype}{logical} :: use\_2018\_answers \textcolor{comment}{! If true use older, less acccurate expressions.}} \DoxyCodeLine{249 } \DoxyCodeLine{250 use\_2018\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) use\_2018\_answers = answers\_2018} \DoxyCodeLine{251 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{252 hneglect = hneglect\_edge\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect} \DoxyCodeLine{253 \textcolor{keywordflow}{else}} \DoxyCodeLine{254 hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect} \DoxyCodeLine{255 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{256 } \DoxyCodeLine{257 \textcolor{comment}{! Loop on interior cells}} \DoxyCodeLine{258 \textcolor{keywordflow}{do} i = 3,n-\/1} \DoxyCodeLine{259 } \DoxyCodeLine{260 h0 = h(i-\/2)} \DoxyCodeLine{261 h1 = h(i-\/1)} \DoxyCodeLine{262 h2 = h(i)} \DoxyCodeLine{263 h3 = h(i+1)} \DoxyCodeLine{264 } \DoxyCodeLine{265 \textcolor{comment}{! Avoid singularities when consecutive pairs of h vanish}} \DoxyCodeLine{266 \textcolor{keywordflow}{if} (h0+h1==0.0 .or. h1+h2==0.0 .or. h2+h3==0.0) \textcolor{keywordflow}{then}} \DoxyCodeLine{267 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{268 h\_min = hminfrac*max( hneglect, h0+h1+h2+h3 )} \DoxyCodeLine{269 \textcolor{keywordflow}{else}} \DoxyCodeLine{270 h\_min = hminfrac*max( hneglect, (h0+h1)+(h2+h3) )} \DoxyCodeLine{271 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{272 h0 = max( h\_min, h(i-\/2) )} \DoxyCodeLine{273 h1 = max( h\_min, h(i-\/1) )} \DoxyCodeLine{274 h2 = max( h\_min, h(i) )} \DoxyCodeLine{275 h3 = max( h\_min, h(i+1) )} \DoxyCodeLine{276 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{277 } \DoxyCodeLine{278 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{279 f1 = (h0+h1) * (h2+h3) / (h1+h2)} \DoxyCodeLine{280 f2 = h2 * u(i-\/1) + h1 * u(i)} \DoxyCodeLine{281 f3 = 1.0 / (h0+h1+h2) + 1.0 / (h1+h2+h3)} \DoxyCodeLine{282 et1 = f1 * f2 * f3} \DoxyCodeLine{283 et2 = ( h2 * (h2+h3) / ( (h0+h1+h2)*(h0+h1) ) ) * \&} \DoxyCodeLine{284 ((h0+2.0*h1) * u(i-\/1) -\/ h1 * u(i-\/2))} \DoxyCodeLine{285 et3 = ( h1 * (h0+h1) / ( (h1+h2+h3)*(h2+h3) ) ) * \&} \DoxyCodeLine{286 ((2.0*h2+h3) * u(i) -\/ h2 * u(i+1))} \DoxyCodeLine{287 edge\_val(i,1) = (et1 + et2 + et3) / ( h0 + h1 + h2 + h3)} \DoxyCodeLine{288 \textcolor{keywordflow}{else}} \DoxyCodeLine{289 i\_h12 = 1.0 / (h1+h2)} \DoxyCodeLine{290 i\_den\_et2 = 1.0 / ( ((h0+h1)+h2)*(h0+h1) ) ; i\_h012 = (h0+h1) * i\_den\_et2} \DoxyCodeLine{291 i\_den\_et3 = 1.0 / ( (h1+(h2+h3))*(h2+h3) ) ; i\_h123 = (h2+h3) * i\_den\_et3} \DoxyCodeLine{292 } \DoxyCodeLine{293 et1 = ( 1.0 + (h1 * i\_h012 + (h0+h1) * i\_h123) ) * i\_h12 * (h2*(h2+h3)) * u(i-\/1) + \&} \DoxyCodeLine{294 ( 1.0 + (h2 * i\_h123 + (h2+h3) * i\_h012) ) * i\_h12 * (h1*(h0+h1)) * u(i)} \DoxyCodeLine{295 et2 = ( h1 * (h2*(h2+h3)) * i\_den\_et2 ) * (u(i-\/1)-\/u(i-\/2))} \DoxyCodeLine{296 et3 = ( h2 * (h1*(h0+h1)) * i\_den\_et3 ) * (u(i) -\/ u(i+1))} \DoxyCodeLine{297 edge\_val(i,1) = (et1 + (et2 + et3)) / ((h0 + h1) + (h2 + h3))} \DoxyCodeLine{298 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{299 edge\_val(i-\/1,2) = edge\_val(i,1)} \DoxyCodeLine{300 } \DoxyCodeLine{301 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on interior cells}} \DoxyCodeLine{302 } \DoxyCodeLine{303 \textcolor{comment}{! Determine first two edge values}} \DoxyCodeLine{304 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{305 h\_min = max( hneglect, hminfrac*sum(h(1:4)) )} \DoxyCodeLine{306 x(1) = 0.0} \DoxyCodeLine{307 \textcolor{keywordflow}{do} i = 1,4} \DoxyCodeLine{308 dx = max(h\_min, h(i) )} \DoxyCodeLine{309 x(i+1) = x(i) + dx} \DoxyCodeLine{310 \textcolor{keywordflow}{do} j = 1,4 ; a(i,j) = ( (x(i+1)**j) -\/ (x(i)**j) ) / real(j) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{311 b(i) = u(i) * dx} \DoxyCodeLine{312 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{313 } \DoxyCodeLine{314 \textcolor{keyword}{call }solve\_linear\_system( a, b, c, 4 )} \DoxyCodeLine{315 } \DoxyCodeLine{316 \textcolor{comment}{! Set the edge values of the first cell}} \DoxyCodeLine{317 edge\_val(1,1) = evaluation\_polynomial( c, 4, x(1) )} \DoxyCodeLine{318 edge\_val(1,2) = evaluation\_polynomial( c, 4, x(2) )} \DoxyCodeLine{319 \textcolor{keywordflow}{else} \textcolor{comment}{! Use expressions with less sensitivity to roundoff}} \DoxyCodeLine{320 \textcolor{keywordflow}{do} i=1,4 ; dz(i) = max(hneglect, h(i) ) ; u\_tmp(i) = u(i) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{321 \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, c)} \DoxyCodeLine{322 } \DoxyCodeLine{323 \textcolor{comment}{! Set the edge values of the first cell}} \DoxyCodeLine{324 edge\_val(1,1) = c(1)} \DoxyCodeLine{325 edge\_val(1,2) = c(1) + dz(1)*(c(2) + dz(1)*(c(3) + dz(1)*c(4)))} \DoxyCodeLine{326 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{327 edge\_val(2,1) = edge\_val(1,2)} \DoxyCodeLine{328 } \DoxyCodeLine{329 \textcolor{comment}{! Determine two edge values of the last cell}} \DoxyCodeLine{330 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{331 h\_min = max( hneglect, hminfrac*sum(h(n-\/3:n)) )} \DoxyCodeLine{332 } \DoxyCodeLine{333 x(1) = 0.0} \DoxyCodeLine{334 \textcolor{keywordflow}{do} i = 1,4} \DoxyCodeLine{335 dx = max(h\_min, h(n-\/4+i) )} \DoxyCodeLine{336 x(i+1) = x(i) + dx} \DoxyCodeLine{337 \textcolor{keywordflow}{do} j = 1,4 ; a(i,j) = ( (x(i+1)**j) -\/ (x(i)**j) ) / real(j) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{338 b(i) = u(n-\/4+i) * dx} \DoxyCodeLine{339 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{340 } \DoxyCodeLine{341 \textcolor{keyword}{call }solve\_linear\_system( a, b, c, 4 )} \DoxyCodeLine{342 } \DoxyCodeLine{343 \textcolor{comment}{! Set the last and second to last edge values}} \DoxyCodeLine{344 edge\_val(n,2) = evaluation\_polynomial( c, 4, x(5) )} \DoxyCodeLine{345 edge\_val(n,1) = evaluation\_polynomial( c, 4, x(4) )} \DoxyCodeLine{346 \textcolor{keywordflow}{else}} \DoxyCodeLine{347 \textcolor{comment}{! Use expressions with less sensitivity to roundoff, including using a coordinate}} \DoxyCodeLine{348 \textcolor{comment}{! system that sets the origin at the last interface in the domain.}} \DoxyCodeLine{349 \textcolor{keywordflow}{do} i=1,4 ; dz(i) = max(hneglect, h(n+1-\/i) ) ; u\_tmp(i) = u(n+1-\/i) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{350 \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, c)} \DoxyCodeLine{351 } \DoxyCodeLine{352 \textcolor{comment}{! Set the last and second to last edge values}} \DoxyCodeLine{353 edge\_val(n,2) = c(1)} \DoxyCodeLine{354 edge\_val(n,1) = c(1) + dz(1)*(c(2) + dz(1)*(c(3) + dz(1)*c(4)))} \DoxyCodeLine{355 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{356 edge\_val(n-\/1,2) = edge\_val(n,1)} \DoxyCodeLine{357 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a0c268712aaf87b3597cef51c85fb32cb}\label{namespaceregrid__edge__values_a0c268712aaf87b3597cef51c85fb32cb}} \index{regrid\_edge\_values@{regrid\_edge\_values}!edge\_values\_implicit\_h4@{edge\_values\_implicit\_h4}} \index{edge\_values\_implicit\_h4@{edge\_values\_implicit\_h4}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{edge\_values\_implicit\_h4()}{edge\_values\_implicit\_h4()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+edge\+\_\+values\+::edge\+\_\+values\+\_\+implicit\+\_\+h4 (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n), intent(in)}]{h, }\item[{real, dimension(n), intent(in)}]{u, }\item[{real, dimension(n,2), intent(inout)}]{edge\+\_\+val, }\item[{real, intent(in), optional}]{h\+\_\+neglect, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})} Compute ih4 edge values (implicit fourth order accurate) in the same units as u. Compute edge values based on fourth-\/order implicit estimates. Fourth-\/order implicit estimates of edge values are based on a two-\/cell stencil. A tridiagonal system is set up and is based on expressing the edge values in terms of neighboring cell averages. The generic relationship is \[ \alpha u_{i-1/2} + u_{i+1/2} + \beta u_{i+3/2} = a \bar{u}_i + b \bar{u}_{i+1} \] and the stencil looks like this \begin{DoxyVerb} i i+1 \end{DoxyVerb} ..--o-\/-\/-\/---o-\/-\/-\/---o--.. i-\/1/2 i+1/2 i+3/2 In this routine, the coefficients $\alpha$, $\beta$, $a$ and $b$ are computed, the tridiagonal system is built, boundary conditions are prescribed and the system is solved to yield edge-\/value estimates. There are N+1 unknowns and we are able to write N-\/1 equations. The boundary conditions close the system. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & Number of cells \\ \hline \mbox{\texttt{ in}} & {\em h} & cell widths \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em edge\+\_\+val} & Returned edge values \mbox{[}A\mbox{]}; the second index is for the two edges of each cell. \\ \hline \mbox{\texttt{ in}} & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\ \hline \end{DoxyParams} Definition at line 387 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{387 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells}} \DoxyCodeLine{388 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths [H]}} \DoxyCodeLine{389 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}} \DoxyCodeLine{390 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Returned edge values [A]; the second index}} \DoxyCodeLine{391 \textcolor{comment}{ !! is for the two edges of each cell.}} \DoxyCodeLine{392 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}} \DoxyCodeLine{393 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}} \DoxyCodeLine{394 } \DoxyCodeLine{395 \textcolor{comment}{! Local variables}} \DoxyCodeLine{396 \textcolor{keywordtype}{integer} :: i, j \textcolor{comment}{! loop indexes}} \DoxyCodeLine{397 \textcolor{keywordtype}{ real} :: h0, h1, h2 \textcolor{comment}{! cell widths [H]}} \DoxyCodeLine{398 \textcolor{keywordtype}{ real} :: h\_min \textcolor{comment}{! A minimal cell width [H]}} \DoxyCodeLine{399 \textcolor{keywordtype}{ real} :: h\_sum \textcolor{comment}{! A sum of adjacent thicknesses [H]}} \DoxyCodeLine{400 \textcolor{keywordtype}{ real} :: h0\_2, h1\_2, h0h1} \DoxyCodeLine{401 \textcolor{keywordtype}{ real} :: h0ph1\_2, h0ph1\_4} \DoxyCodeLine{402 \textcolor{keywordtype}{ real} :: alpha, beta \textcolor{comment}{! stencil coefficients [nondim]}} \DoxyCodeLine{403 \textcolor{keywordtype}{ real} :: I\_h2, abmix \textcolor{comment}{! stencil coefficients [nondim]}} \DoxyCodeLine{404 \textcolor{keywordtype}{ real} :: a, b} \DoxyCodeLine{405 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(5)} :: x \textcolor{comment}{! Coordinate system with 0 at edges [H]}} \DoxyCodeLine{406 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{parameter} :: C1\_12 = 1.0 / 12.0} \DoxyCodeLine{407 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{parameter} :: C1\_3 = 1.0 / 3.0} \DoxyCodeLine{408 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)} :: dz \textcolor{comment}{! A temporary array of limited layer thicknesses [H]}} \DoxyCodeLine{409 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)} :: u\_tmp \textcolor{comment}{! A temporary array of cell average properties [A]}} \DoxyCodeLine{410 \textcolor{keywordtype}{ real} :: dx, xavg \textcolor{comment}{! Differences and averages of successive values of x [H]}} \DoxyCodeLine{411 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4,4)} :: Asys \textcolor{comment}{! boundary conditions}} \DoxyCodeLine{412 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)} :: Bsys, Csys} \DoxyCodeLine{413 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N+1)} :: tri\_l, \& \textcolor{comment}{! tridiagonal system (lower diagonal) [nondim]}} \DoxyCodeLine{414 tri\_d, \& \textcolor{comment}{! tridiagonal system (middle diagonal) [nondim]}} \DoxyCodeLine{415 tri\_c, \& \textcolor{comment}{! tridiagonal system central value, with tri\_d = tri\_c+tri\_l+tri\_u}} \DoxyCodeLine{416 tri\_u, \& \textcolor{comment}{! tridiagonal system (upper diagonal) [nondim]}} \DoxyCodeLine{417 tri\_b, \& \textcolor{comment}{! tridiagonal system (right hand side) [A]}} \DoxyCodeLine{418 tri\_x \textcolor{comment}{! tridiagonal system (solution vector) [A]}} \DoxyCodeLine{419 \textcolor{keywordtype}{ real} :: hNeglect \textcolor{comment}{! A negligible thickness [H]}} \DoxyCodeLine{420 \textcolor{keywordtype}{logical} :: use\_2018\_answers \textcolor{comment}{! If true use older, less acccurate expressions.}} \DoxyCodeLine{421 } \DoxyCodeLine{422 use\_2018\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) use\_2018\_answers = answers\_2018} \DoxyCodeLine{423 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{424 hneglect = hneglect\_edge\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect} \DoxyCodeLine{425 \textcolor{keywordflow}{else}} \DoxyCodeLine{426 hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect} \DoxyCodeLine{427 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{428 } \DoxyCodeLine{429 \textcolor{comment}{! Loop on cells (except last one)}} \DoxyCodeLine{430 \textcolor{keywordflow}{do} i = 1,n-\/1} \DoxyCodeLine{431 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{432 \textcolor{comment}{! Get cell widths}} \DoxyCodeLine{433 h0 = h(i)} \DoxyCodeLine{434 h1 = h(i+1)} \DoxyCodeLine{435 \textcolor{comment}{! Avoid singularities when h0+h1=0}} \DoxyCodeLine{436 \textcolor{keywordflow}{if} (h0+h1==0.) \textcolor{keywordflow}{then}} \DoxyCodeLine{437 h0 = hneglect} \DoxyCodeLine{438 h1 = hneglect} \DoxyCodeLine{439 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{440 } \DoxyCodeLine{441 \textcolor{comment}{! Auxiliary calculations}} \DoxyCodeLine{442 h0ph1\_2 = (h0 + h1)**2} \DoxyCodeLine{443 h0ph1\_4 = h0ph1\_2**2} \DoxyCodeLine{444 h0h1 = h0 * h1} \DoxyCodeLine{445 h0\_2 = h0 * h0} \DoxyCodeLine{446 h1\_2 = h1 * h1} \DoxyCodeLine{447 } \DoxyCodeLine{448 \textcolor{comment}{! Coefficients}} \DoxyCodeLine{449 alpha = h1\_2 / h0ph1\_2} \DoxyCodeLine{450 beta = h0\_2 / h0ph1\_2} \DoxyCodeLine{451 a = 2.0 * h1\_2 * ( h1\_2 + 2.0 * h0\_2 + 3.0 * h0h1 ) / h0ph1\_4} \DoxyCodeLine{452 b = 2.0 * h0\_2 * ( h0\_2 + 2.0 * h1\_2 + 3.0 * h0h1 ) / h0ph1\_4} \DoxyCodeLine{453 } \DoxyCodeLine{454 tri\_d(i+1) = 1.0} \DoxyCodeLine{455 \textcolor{keywordflow}{else} \textcolor{comment}{! Use expressions with less sensitivity to roundoff}} \DoxyCodeLine{456 \textcolor{comment}{! Get cell widths}} \DoxyCodeLine{457 h0 = max(h(i), hneglect)} \DoxyCodeLine{458 h1 = max(h(i+1), hneglect)} \DoxyCodeLine{459 \textcolor{comment}{! The 1e-\/12 here attempts to balance truncation errors from the differences of}} \DoxyCodeLine{460 \textcolor{comment}{! large numbers against errors from approximating thin layers as non-\/vanishing.}} \DoxyCodeLine{461 \textcolor{keywordflow}{if} (abs(h0) < 1.0e-\/12*abs(h1)) h0 = 1.0e-\/12*h1} \DoxyCodeLine{462 \textcolor{keywordflow}{if} (abs(h1) < 1.0e-\/12*abs(h0)) h1 = 1.0e-\/12*h0} \DoxyCodeLine{463 i\_h2 = 1.0 / ((h0 + h1)**2)} \DoxyCodeLine{464 alpha = (h1 * h1) * i\_h2} \DoxyCodeLine{465 beta = (h0 * h0) * i\_h2} \DoxyCodeLine{466 abmix = (h0 * h1) * i\_h2} \DoxyCodeLine{467 a = 2.0 * alpha * ( alpha + 2.0 * beta + 3.0 * abmix )} \DoxyCodeLine{468 b = 2.0 * beta * ( beta + 2.0 * alpha + 3.0 * abmix )} \DoxyCodeLine{469 } \DoxyCodeLine{470 tri\_c(i+1) = 2.0*abmix \textcolor{comment}{! = 1.0 -\/ alpha -\/ beta}} \DoxyCodeLine{471 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{472 } \DoxyCodeLine{473 tri\_l(i+1) = alpha} \DoxyCodeLine{474 tri\_u(i+1) = beta} \DoxyCodeLine{475 } \DoxyCodeLine{476 tri\_b(i+1) = a * u(i) + b * u(i+1)} \DoxyCodeLine{477 } \DoxyCodeLine{478 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on cells}} \DoxyCodeLine{479 } \DoxyCodeLine{480 \textcolor{comment}{! Boundary conditions: set the first boundary value}} \DoxyCodeLine{481 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{482 h\_min = max( hneglect, hminfrac*sum(h(1:4)) )} \DoxyCodeLine{483 x(1) = 0.0} \DoxyCodeLine{484 \textcolor{keywordflow}{do} i = 1,4} \DoxyCodeLine{485 dx = max(h\_min, h(i) )} \DoxyCodeLine{486 x(i+1) = x(i) + dx} \DoxyCodeLine{487 \textcolor{keywordflow}{do} j = 1,4 ; asys(i,j) = ( (x(i+1)**j) -\/ (x(i)**j) ) / j ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{488 bsys(i) = u(i) * dx} \DoxyCodeLine{489 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{490 } \DoxyCodeLine{491 \textcolor{keyword}{call }solve\_linear\_system( asys, bsys, csys, 4 )} \DoxyCodeLine{492 } \DoxyCodeLine{493 tri\_b(1) = evaluation\_polynomial( csys, 4, x(1) ) \textcolor{comment}{! Set the first edge value}} \DoxyCodeLine{494 tri\_d(1) = 1.0} \DoxyCodeLine{495 \textcolor{keywordflow}{else} \textcolor{comment}{! Use expressions with less sensitivity to roundoff}} \DoxyCodeLine{496 \textcolor{keywordflow}{do} i=1,4 ; dz(i) = max(hneglect, h(i) ) ; u\_tmp(i) = u(i) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{497 \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, csys)} \DoxyCodeLine{498 } \DoxyCodeLine{499 tri\_b(1) = csys(1) \textcolor{comment}{! Set the first edge value.}} \DoxyCodeLine{500 tri\_c(1) = 1.0} \DoxyCodeLine{501 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{502 tri\_u(1) = 0.0 \textcolor{comment}{! tri\_l(1) = 0.0}} \DoxyCodeLine{503 } \DoxyCodeLine{504 \textcolor{comment}{! Boundary conditions: set the last boundary value}} \DoxyCodeLine{505 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{506 h\_min = max( hneglect, hminfrac*sum(h(n-\/3:n)) )} \DoxyCodeLine{507 x(1) = 0.0} \DoxyCodeLine{508 \textcolor{keywordflow}{do} i=1,4} \DoxyCodeLine{509 dx = max(h\_min, h(n-\/4+i) )} \DoxyCodeLine{510 x(i+1) = x(i) + dx} \DoxyCodeLine{511 \textcolor{keywordflow}{do} j = 1,4 ; asys(i,j) = ( (x(i+1)**j) -\/ (x(i)**j) ) / j ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{512 bsys(i) = u(n-\/4+i) * dx} \DoxyCodeLine{513 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{514 } \DoxyCodeLine{515 \textcolor{keyword}{call }solve\_linear\_system( asys, bsys, csys, 4 )} \DoxyCodeLine{516 } \DoxyCodeLine{517 \textcolor{comment}{! Set the last edge value}} \DoxyCodeLine{518 tri\_b(n+1) = evaluation\_polynomial( csys, 4, x(5) )} \DoxyCodeLine{519 tri\_d(n+1) = 1.0} \DoxyCodeLine{520 } \DoxyCodeLine{521 \textcolor{keywordflow}{else}} \DoxyCodeLine{522 \textcolor{comment}{! Use expressions with less sensitivity to roundoff, including using a coordinate}} \DoxyCodeLine{523 \textcolor{comment}{! system that sets the origin at the last interface in the domain.}} \DoxyCodeLine{524 \textcolor{keywordflow}{do} i=1,4 ; dz(i) = max(hneglect, h(n+1-\/i) ) ; u\_tmp(i) = u(n+1-\/i) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{525 \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, csys)} \DoxyCodeLine{526 } \DoxyCodeLine{527 tri\_b(n+1) = csys(1) \textcolor{comment}{! Set the last edge value}} \DoxyCodeLine{528 tri\_c(n+1) = 1.0} \DoxyCodeLine{529 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{530 tri\_l(n+1) = 0.0 \textcolor{comment}{! tri\_u(N+1) = 0.0}} \DoxyCodeLine{531 } \DoxyCodeLine{532 \textcolor{comment}{! Solve tridiagonal system and assign edge values}} \DoxyCodeLine{533 \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}} \DoxyCodeLine{534 \textcolor{keyword}{call }solve\_tridiagonal\_system( tri\_l, tri\_d, tri\_u, tri\_b, tri\_x, n+1 )} \DoxyCodeLine{535 \textcolor{keywordflow}{else}} \DoxyCodeLine{536 \textcolor{keyword}{call }solve\_diag\_dominant\_tridiag( tri\_l, tri\_c, tri\_u, tri\_b, tri\_x, n+1 )} \DoxyCodeLine{537 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{538 } \DoxyCodeLine{539 edge\_val(1,1) = tri\_x(1)} \DoxyCodeLine{540 \textcolor{keywordflow}{do} i=2,n} \DoxyCodeLine{541 edge\_val(i,1) = tri\_x(i)} \DoxyCodeLine{542 edge\_val(i-\/1,2) = tri\_x(i)} \DoxyCodeLine{543 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{544 edge\_val(n,2) = tri\_x(n+1)} \DoxyCodeLine{545 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a9955c45dcd1bfec32fbf5602315cb5b1}\label{namespaceregrid__edge__values_a9955c45dcd1bfec32fbf5602315cb5b1}} \index{regrid\_edge\_values@{regrid\_edge\_values}!edge\_values\_implicit\_h6@{edge\_values\_implicit\_h6}} \index{edge\_values\_implicit\_h6@{edge\_values\_implicit\_h6}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{edge\_values\_implicit\_h6()}{edge\_values\_implicit\_h6()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+edge\+\_\+values\+::edge\+\_\+values\+\_\+implicit\+\_\+h6 (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n), intent(in)}]{h, }\item[{real, dimension(n), intent(in)}]{u, }\item[{real, dimension(n,2), intent(inout)}]{edge\+\_\+val, }\item[{real, intent(in), optional}]{h\+\_\+neglect, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})} Compute ih6 edge values (implicit sixth order accurate) in the same units as u. Sixth-\/order implicit estimates of edge values are based on a four-\/cell, three-\/edge stencil. A tridiagonal system is set up and is based on expressing the edge values in terms of neighboring cell averages. The generic relationship is \[ \alpha u_{i-1/2} + u_{i+1/2} + \beta u_{i+3/2} = a \bar{u}_{i-1} + b \bar{u}_i + c \bar{u}_{i+1} + d \bar{u}_{i+2} \] and the stencil looks like this \begin{DoxyVerb} i-1 i i+1 i+2 \end{DoxyVerb} ..--o-\/-\/-\/---o-\/-\/-\/---o-\/-\/-\/---o-\/-\/-\/---o--.. i-\/1/2 i+1/2 i+3/2 In this routine, the coefficients $\alpha$, $\beta$, a, b, c and d are computed, the tridiagonal system is built, boundary conditions are prescribed and the system is solved to yield edge-\/value estimates. This scheme is described in detail by White and Adcroft, 2009, J. Comp. Phys, \href{https://doi.org/10.1016/j.jcp.2008.04.026}{\texttt{ https\+://doi.\+org/10.\+1016/j.\+jcp.\+2008.\+04.\+026}} Note that the centered stencil only applies to edges 3 to N-\/1 (edges are numbered 1 to n+1), which yields N-\/3 equations for N+1 unknowns. Two other equations are written by using a right-\/biased stencil for edge 2 and a left-\/biased stencil for edge N. The prescription of boundary conditions (using sixth-\/order polynomials) closes the system. C\+A\+U\+T\+I\+ON\+: For each edge, in order to determine the coefficients of the implicit expression, a 6x6 linear system is solved. This may become computationally expensive if regridding is carried out often. Figuring out closed-\/form expressions for these coefficients on nonuniform meshes turned out to be intractable. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & Number of cells \\ \hline \mbox{\texttt{ in}} & {\em h} & cell widths \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em u} & cell average properties (size N) in arbitrary units \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em edge\+\_\+val} & Returned edge values \mbox{[}A\mbox{]}; the second index is for the two edges of each cell. \\ \hline \mbox{\texttt{ in}} & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]} \\ \hline \mbox{\texttt{ in}} & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\ \hline \end{DoxyParams} Definition at line 1137 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{1137 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< Number of cells}} \DoxyCodeLine{1138 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{ !< cell widths [H]}} \DoxyCodeLine{1139 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< cell average properties (size N) in arbitrary units [A]}} \DoxyCodeLine{1140 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Returned edge values [A]; the second index}} \DoxyCodeLine{1141 \textcolor{comment}{ !! is for the two edges of each cell.}} \DoxyCodeLine{1142 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}} \DoxyCodeLine{1143 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}} \DoxyCodeLine{1144 } \DoxyCodeLine{1145 \textcolor{comment}{! Local variables}} \DoxyCodeLine{1146 \textcolor{keywordtype}{ real} :: h0, h1, h2, h3 \textcolor{comment}{! cell widths [H]}} \DoxyCodeLine{1147 \textcolor{keywordtype}{ real} :: hMin \textcolor{comment}{! The minimum thickness used in these calculations [H]}} \DoxyCodeLine{1148 \textcolor{keywordtype}{ real} :: h01, h01\_2, h01\_3 \textcolor{comment}{! Summed thicknesses to various powers [H\string^n \string~> m\string^n or kg\string^n m-\/2n]}} \DoxyCodeLine{1149 \textcolor{keywordtype}{ real} :: h23, h23\_2, h23\_3 \textcolor{comment}{! Summed thicknesses to various powers [H\string^n \string~> m\string^n or kg\string^n m-\/2n]}} \DoxyCodeLine{1150 \textcolor{keywordtype}{ real} :: hNeglect \textcolor{comment}{! A negligible thickness [H].}} \DoxyCodeLine{1151 \textcolor{keywordtype}{ real} :: h1\_2, h2\_2, h1\_3, h2\_3 \textcolor{comment}{! Cell widths raised to the 2nd and 3rd powers [H2] or [H3]}} \DoxyCodeLine{1152 \textcolor{keywordtype}{ real} :: h1\_4, h2\_4, h1\_5, h2\_5 \textcolor{comment}{! Cell widths raised to the 4th and 5th powers [H4] or [H5]}} \DoxyCodeLine{1153 \textcolor{keywordtype}{ real} :: alpha, beta \textcolor{comment}{! stencil coefficients}} \DoxyCodeLine{1154 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(7)} :: x \textcolor{comment}{! Coordinate system with 0 at edges [same units as h]}} \DoxyCodeLine{1155 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{parameter} :: C1\_12 = 1.0 / 12.0} \DoxyCodeLine{1156 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{parameter} :: C5\_6 = 5.0 / 6.0} \DoxyCodeLine{1157 \textcolor{keywordtype}{ real} :: dx, xavg \textcolor{comment}{! Differences and averages of successive values of x [same units as h]}} \DoxyCodeLine{1158 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(6,6)} :: Asys \textcolor{comment}{! boundary conditions}} \DoxyCodeLine{1159 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(6)} :: Bsys, Csys \textcolor{comment}{! ...}} \DoxyCodeLine{1160 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N+1)} :: tri\_l, \& \textcolor{comment}{! trid. system (lower diagonal)}} \DoxyCodeLine{1161 tri\_d, \& \textcolor{comment}{! trid. system (middle diagonal)}} \DoxyCodeLine{1162 tri\_u, \& \textcolor{comment}{! trid. system (upper diagonal)}} \DoxyCodeLine{1163 tri\_b, \& \textcolor{comment}{! trid. system (unknowns vector)}} \DoxyCodeLine{1164 tri\_x \textcolor{comment}{! trid. system (rhs)}} \DoxyCodeLine{1165 \textcolor{keywordtype}{integer} :: i, j, k \textcolor{comment}{! loop indexes}} \DoxyCodeLine{1166 } \DoxyCodeLine{1167 hneglect = hneglect\_edge\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect} \DoxyCodeLine{1168 } \DoxyCodeLine{1169 \textcolor{comment}{! Loop on interior cells}} \DoxyCodeLine{1170 \textcolor{keywordflow}{do} k = 2,n-\/2} \DoxyCodeLine{1171 \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}} \DoxyCodeLine{1172 hmin = max(hneglect, hminfrac*((h(k-\/1) + h(k)) + (h(k+1) + h(k+2))))} \DoxyCodeLine{1173 h0 = max(h(k-\/1), hmin) ; h1 = max(h(k), hmin)} \DoxyCodeLine{1174 h2 = max(h(k+1), hmin) ; h3 = max(h(k+2), hmin)} \DoxyCodeLine{1175 } \DoxyCodeLine{1176 \textcolor{comment}{! Auxiliary calculations}} \DoxyCodeLine{1177 h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2} \DoxyCodeLine{1178 h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2} \DoxyCodeLine{1179 } \DoxyCodeLine{1180 \textcolor{comment}{! Compute matrix entries as described in Eq. (48) of White and Adcroft (2009)}} \DoxyCodeLine{1181 asys(1:6,1) = (/ 1.0, 1.0, -\/1.0, -\/1.0, -\/1.0, -\/1.0 /)} \DoxyCodeLine{1182 asys(1:6,2) = (/ -\/2.0*h1, 2.0*h2, (2.0*h1 + h0), h1, -\/h2, -\/(2.0*h2 + h3) /)} \DoxyCodeLine{1183 asys(1:6,3) = (/ 3.0*h1\_2, 3.0*h2\_2, -\/(3.0*h1\_2 + h0*(3.0*h1 + h0)), \& \textcolor{comment}{! = -\/((h0+h1)**3 -\/ h1**3) / h0}} \DoxyCodeLine{1184 -\/h1\_2, -\/h2\_2, -\/(3.0*h2\_2 + h3*(3.0*h2 + h3)) /) \textcolor{comment}{! = -\/((h2+h3)**3 -\/ h2**3) / h3}} \DoxyCodeLine{1185 asys(1:6,4) = (/ -\/4.0*h1\_3, 4.0*h2\_3, (4.0*h1\_3 + h0*(6.0*h1\_2 + h0*(4.0*h1 + h0))), \&} \DoxyCodeLine{1186 h1\_3, -\/h2\_3, -\/(4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)} \DoxyCodeLine{1187 asys(1:6,5) = (/ 5.0*h1\_4, 5.0*h2\_4, -\/(5.0*h1\_4 + h0*(10.0*h1\_3 + h0*(10.0*h1\_2 + h0*(5.0*h1 + h0)))), \&} \DoxyCodeLine{1188 -\/h1\_4, -\/h2\_4, -\/(5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))) /)} \DoxyCodeLine{1189 asys(1:6,6) = (/ -\/6.0*h1\_5, 6.0*h2\_5, \&} \DoxyCodeLine{1190 (6.0*h1\_5 + h0*(15.0*h1\_4 + h0*(20.0*h1\_3 + h0*(15.0*h1\_2 + h0*(6.0*h1 + h0))))), \&} \DoxyCodeLine{1191 h1\_5, -\/h2\_5, \&} \DoxyCodeLine{1192 -\/(6.0*h2\_5 + h3*(15.0*h2\_4 + h3*(20.0*h2\_3 + h3*(15.0*h2\_2 + h3*(6.0*h2 + h3))))) /)} \DoxyCodeLine{1193 bsys(1:6) = (/ -\/1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)} \DoxyCodeLine{1194 } \DoxyCodeLine{1195 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{1196 } \DoxyCodeLine{1197 alpha = csys(1)} \DoxyCodeLine{1198 beta = csys(2)} \DoxyCodeLine{1199 } \DoxyCodeLine{1200 tri\_l(k+1) = alpha} \DoxyCodeLine{1201 tri\_d(k+1) = 1.0} \DoxyCodeLine{1202 tri\_u(k+1) = beta} \DoxyCodeLine{1203 tri\_b(k+1) = csys(3) * u(k-\/1) + csys(4) * u(k) + csys(5) * u(k+1) + csys(6) * u(k+2)} \DoxyCodeLine{1204 } \DoxyCodeLine{1205 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on cells}} \DoxyCodeLine{1206 } \DoxyCodeLine{1207 \textcolor{comment}{! Use a right-\/biased stencil for the second row, as described in Eq. (49) of White and Adcroft (2009).}} \DoxyCodeLine{1208 } \DoxyCodeLine{1209 \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}} \DoxyCodeLine{1210 hmin = max(hneglect, hminfrac*((h(1) + h(2)) + (h(3) + h(4))))} \DoxyCodeLine{1211 h0 = max(h(1), hmin) ; h1 = max(h(2), hmin)} \DoxyCodeLine{1212 h2 = max(h(3), hmin) ; h3 = max(h(4), hmin)} \DoxyCodeLine{1213 } \DoxyCodeLine{1214 \textcolor{comment}{! Auxiliary calculations}} \DoxyCodeLine{1215 h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2} \DoxyCodeLine{1216 h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2} \DoxyCodeLine{1217 h01 = h0 + h1 ; h01\_2 = h01 * h01 ; h01\_3 = h01 * h01\_2} \DoxyCodeLine{1218 } \DoxyCodeLine{1219 \textcolor{comment}{! Compute matrix entries}} \DoxyCodeLine{1220 asys(1:6,1) = (/ 1.0, 1.0, -\/1.0, -\/1.0, -\/1.0, -\/1.0 /)} \DoxyCodeLine{1221 asys(1:6,2) = (/ -\/2.0*h01, 0.0, (2.0*h1 + h0), h1, -\/h2, -\/(2.0*h2 + h3) /)} \DoxyCodeLine{1222 asys(1:6,3) = (/ 3.0*h01\_2, 0.0, -\/(3.0*h1\_2 + h0*(3.0*h1 + h0)), \&} \DoxyCodeLine{1223 -\/h1\_2, -\/h2\_2, -\/(3.0*h2\_2 + h3*(3.0*h2 + h3)) /)} \DoxyCodeLine{1224 asys(1:6,4) = (/ -\/4.0*h01\_3, 0.0, (4.0*h1\_3 + h0*(6.0*h1\_2 + h0*(4.0*h1 + h0))), \&} \DoxyCodeLine{1225 h1\_3, -\/h2\_3, -\/(4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)} \DoxyCodeLine{1226 asys(1:6,5) = (/ 5.0*(h01\_2*h01\_2), 0.0, -\/(5.0*h1\_4 + h0*(10.0*h1\_3 + h0*(10.0*h1\_2 + h0*(5.0*h1 + h0)))), \&} \DoxyCodeLine{1227 -\/h1\_4, -\/h2\_4, -\/(5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))) /)} \DoxyCodeLine{1228 asys(1:6,6) = (/ -\/6.0*(h01\_3*h01\_2), 0.0, \&} \DoxyCodeLine{1229 (6.0*h1\_5 + h0*(15.0*h1\_4 + h0*(20.0*h1\_3 + h0*(15.0*h1\_2 + h0*(6.0*h1 + h0))))), \&} \DoxyCodeLine{1230 h1\_5, -\/ h2\_5, \&} \DoxyCodeLine{1231 -\/(6.0*h2\_5 + h3*(15.0*h2\_4 + h3*(20.0*h2\_3 + h3*(15.0*h2\_2 + h3*(6.0*h2 + h3))))) /)} \DoxyCodeLine{1232 bsys(1:6) = (/ -\/1.0, 2.0*h1, -\/3.0*h1\_2, 4.0*h1\_3, -\/5.0*h1\_4, 6.0*h1\_5 /)} \DoxyCodeLine{1233 } \DoxyCodeLine{1234 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{1235 } \DoxyCodeLine{1236 alpha = csys(1)} \DoxyCodeLine{1237 beta = csys(2)} \DoxyCodeLine{1238 } \DoxyCodeLine{1239 tri\_l(2) = alpha} \DoxyCodeLine{1240 tri\_d(2) = 1.0} \DoxyCodeLine{1241 tri\_u(2) = beta} \DoxyCodeLine{1242 tri\_b(2) = csys(3) * u(1) + csys(4) * u(2) + csys(5) * u(3) + csys(6) * u(4)} \DoxyCodeLine{1243 } \DoxyCodeLine{1244 \textcolor{comment}{! Boundary conditions: left boundary}} \DoxyCodeLine{1245 hmin = max( hneglect, hminfrac*((h(1)+h(2)) + (h(5)+h(6)) + (h(3)+h(4))) )} \DoxyCodeLine{1246 x(1) = 0.0} \DoxyCodeLine{1247 \textcolor{keywordflow}{do} i = 1,6} \DoxyCodeLine{1248 dx = max( hmin, h(i) )} \DoxyCodeLine{1249 xavg = x(i) + 0.5*dx} \DoxyCodeLine{1250 asys(1:6,i) = (/ 1.0, xavg, (xavg**2 + c1\_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), \&} \DoxyCodeLine{1251 (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), \&} \DoxyCodeLine{1252 xavg * (xavg**4 + c5\_6*xavg**2*dx**2 + 0.0625*dx**4) /)} \DoxyCodeLine{1253 bsys(i) = u(i)} \DoxyCodeLine{1254 x(i+1) = x(i) + dx} \DoxyCodeLine{1255 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1256 } \DoxyCodeLine{1257 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{1258 } \DoxyCodeLine{1259 tri\_l(1) = 0.0} \DoxyCodeLine{1260 tri\_d(1) = 1.0} \DoxyCodeLine{1261 tri\_u(1) = 0.0} \DoxyCodeLine{1262 tri\_b(1) = evaluation\_polynomial( csys, 6, x(1) ) \textcolor{comment}{! first edge value}} \DoxyCodeLine{1263 } \DoxyCodeLine{1264 \textcolor{comment}{! Use a left-\/biased stencil for the second to last row, as described in Eq. (50) of White and Adcroft (2009).}} \DoxyCodeLine{1265 } \DoxyCodeLine{1266 \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}} \DoxyCodeLine{1267 hmin = max(hneglect, hminfrac*((h(n-\/3) + h(n-\/2)) + (h(n-\/1) + h(n))))} \DoxyCodeLine{1268 h0 = max(h(n-\/3), hmin) ; h1 = max(h(n-\/2), hmin)} \DoxyCodeLine{1269 h2 = max(h(n-\/1), hmin) ; h3 = max(h(n), hmin)} \DoxyCodeLine{1270 } \DoxyCodeLine{1271 \textcolor{comment}{! Auxiliary calculations}} \DoxyCodeLine{1272 h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2} \DoxyCodeLine{1273 h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2} \DoxyCodeLine{1274 h23 = h2 + h3 ; h23\_2 = h23 * h23 ; h23\_3 = h23 * h23\_2} \DoxyCodeLine{1275 } \DoxyCodeLine{1276 \textcolor{comment}{! Compute matrix entries}} \DoxyCodeLine{1277 asys(1:6,1) = (/ 1.0, 1.0, -\/1.0, -\/1.0, -\/1.0, -\/1.0 /)} \DoxyCodeLine{1278 asys(1:6,2) = (/ 0.0, 2.0*h23, (2.0*h1 + h0), h1, -\/h2, -\/(2.0*h2 + h3) /)} \DoxyCodeLine{1279 asys(1:6,3) = (/ 0.0, 3.0*h23\_2, -\/(3.0*h1\_2 + h0*(3.0*h1 + h0)), \&} \DoxyCodeLine{1280 -\/h1\_2, -\/h2\_2, -\/(3.0*h2\_2 + h3*(3.0*h2 + h3)) /)} \DoxyCodeLine{1281 asys(1:6,4) = (/ 0.0, 4.0*h23\_3, (4.0*h1\_3 + h0*(6.0*h1\_2 + h0*(4.0*h1 + h0))), \&} \DoxyCodeLine{1282 h1\_3, -\/h2\_3, -\/(4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)} \DoxyCodeLine{1283 asys(1:6,5) = (/ 0.0, 5.0*(h23\_2*h23\_2), -\/(5.0*h1\_4 + h0*(10.0*h1\_3 + h0*(10.0*h1\_2 + h0*(5.0*h1 + h0)))), \&} \DoxyCodeLine{1284 -\/h1\_4, -\/h2\_4, -\/(5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))) /)} \DoxyCodeLine{1285 asys(1:6,6) = (/ 0.0, 6.0*(h23\_3*h23\_2), \&} \DoxyCodeLine{1286 (6.0*h1\_5 + h0*(15.0*h1\_4 + h0*(20.0*h1\_3 + h0*(15.0*h1\_2 + h0*(6.0*h1 + h0))))), \&} \DoxyCodeLine{1287 h1\_5, -\/h2\_5, \&} \DoxyCodeLine{1288 -\/(6.0*h2\_5 + h3*(15.0*h2\_4 + h3*(20.0*h2\_3 + h3*(15.0*h2\_2 + h3*(6.0*h2 + h3))))) /)} \DoxyCodeLine{1289 bsys(1:6) = (/ -\/1.0, -\/2.0*h2, -\/3.0*h2\_2, -\/4.0*h2\_3, -\/5.0*h2\_4, -\/6.0*h2\_5 /)} \DoxyCodeLine{1290 } \DoxyCodeLine{1291 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{1292 } \DoxyCodeLine{1293 alpha = csys(1)} \DoxyCodeLine{1294 beta = csys(2)} \DoxyCodeLine{1295 } \DoxyCodeLine{1296 tri\_l(n) = alpha} \DoxyCodeLine{1297 tri\_d(n) = 1.0} \DoxyCodeLine{1298 tri\_u(n) = beta} \DoxyCodeLine{1299 tri\_b(n) = csys(3) * u(n-\/3) + csys(4) * u(n-\/2) + csys(5) * u(n-\/1) + csys(6) * u(n)} \DoxyCodeLine{1300 } \DoxyCodeLine{1301 \textcolor{comment}{! Boundary conditions: right boundary}} \DoxyCodeLine{1302 hmin = max( hneglect, hminfrac*(h(n-\/3) + h(n-\/2)) + ((h(n-\/1) + h(n)) + (h(n-\/5) + h(n-\/4))) )} \DoxyCodeLine{1303 x(1) = 0.0} \DoxyCodeLine{1304 \textcolor{keywordflow}{do} i = 1,6} \DoxyCodeLine{1305 dx = max( hmin, h(n+1-\/i) )} \DoxyCodeLine{1306 xavg = x(i) + 0.5 * dx} \DoxyCodeLine{1307 asys(1:6,i) = (/ 1.0, xavg, (xavg**2 + c1\_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), \&} \DoxyCodeLine{1308 (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), \&} \DoxyCodeLine{1309 xavg * (xavg**4 + c5\_6*xavg**2*dx**2 + 0.0625*dx**4) /)} \DoxyCodeLine{1310 bsys(i) = u(n+1-\/i)} \DoxyCodeLine{1311 x(i+1) = x(i) + dx} \DoxyCodeLine{1312 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1313 } \DoxyCodeLine{1314 \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )} \DoxyCodeLine{1315 } \DoxyCodeLine{1316 tri\_l(n+1) = 0.0} \DoxyCodeLine{1317 tri\_d(n+1) = 1.0} \DoxyCodeLine{1318 tri\_u(n+1) = 0.0} \DoxyCodeLine{1319 tri\_b(n+1) = csys(1)} \DoxyCodeLine{1320 } \DoxyCodeLine{1321 \textcolor{comment}{! Solve tridiagonal system and assign edge values}} \DoxyCodeLine{1322 \textcolor{keyword}{call }solve\_tridiagonal\_system( tri\_l, tri\_d, tri\_u, tri\_b, tri\_x, n+1 )} \DoxyCodeLine{1323 } \DoxyCodeLine{1324 \textcolor{keywordflow}{do} i = 2,n} \DoxyCodeLine{1325 edge\_val(i,1) = tri\_x(i)} \DoxyCodeLine{1326 edge\_val(i-\/1,2) = tri\_x(i)} \DoxyCodeLine{1327 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1328 edge\_val(1,1) = tri\_x(1)} \DoxyCodeLine{1329 edge\_val(n,2) = tri\_x(n+1)} \DoxyCodeLine{1330 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a7082604865f485ac18af8b44a7aa0224}\label{namespaceregrid__edge__values_a7082604865f485ac18af8b44a7aa0224}} \index{regrid\_edge\_values@{regrid\_edge\_values}!end\_value\_h4@{end\_value\_h4}} \index{end\_value\_h4@{end\_value\_h4}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{end\_value\_h4()}{end\_value\_h4()}} {\footnotesize\ttfamily subroutine regrid\+\_\+edge\+\_\+values\+::end\+\_\+value\+\_\+h4 (\begin{DoxyParamCaption}\item[{real, dimension(4), intent(in)}]{dz, }\item[{real, dimension(4), intent(in)}]{u, }\item[{real, dimension(4), intent(out)}]{Csys }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}} Determine a one-\/sided 4th order polynomial fit of u to the data points for the purposes of specifying edge values, as described in the appendix of White and Adcroft J\+CP 2008. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em dz} & The thicknesses of 4 layers, starting at the edge \mbox{[}H\mbox{]}. The values of dz must be positive. \\ \hline \mbox{\texttt{ in}} & {\em u} & The average properties of 4 layers, starting at the edge \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ out}} & {\em csys} & The four coefficients of a 4th order polynomial fit of u as a function of z \mbox{[}A H-\/(n-\/1)\mbox{]} \\ \hline \end{DoxyParams} Definition at line 551 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{551 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)}, \textcolor{keywordtype}{intent(in)} :: dz\textcolor{comment}{ !< The thicknesses of 4 layers, starting at the edge [H].}} \DoxyCodeLine{552 \textcolor{comment}{ !! The values of dz must be positive.}} \DoxyCodeLine{553 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)}, \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{ !< The average properties of 4 layers, starting at the edge [A]}} \DoxyCodeLine{554 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)}, \textcolor{keywordtype}{intent(out)} :: Csys\textcolor{comment}{ !< The four coefficients of a 4th order polynomial fit}} \DoxyCodeLine{555 \textcolor{comment}{ !! of u as a function of z [A H-\/(n-\/1)]}} \DoxyCodeLine{556 } \DoxyCodeLine{557 \textcolor{comment}{! Local variables}} \DoxyCodeLine{558 \textcolor{keywordtype}{ real} :: Wt(3,4) \textcolor{comment}{! The weights of successive u differences in the 4 closed form expressions.}} \DoxyCodeLine{559 \textcolor{comment}{! The units of Wt vary with the second index as [H-\/(n-\/1)].}} \DoxyCodeLine{560 \textcolor{keywordtype}{ real} :: h1, h2, h3, h4 \textcolor{comment}{! Copies of the layer thicknesses [H]}} \DoxyCodeLine{561 \textcolor{keywordtype}{ real} :: h12, h23, h34 \textcolor{comment}{! Sums of two successive thicknesses [H]}} \DoxyCodeLine{562 \textcolor{keywordtype}{ real} :: h123, h234 \textcolor{comment}{! Sums of three successive thicknesses [H]}} \DoxyCodeLine{563 \textcolor{keywordtype}{ real} :: h1234 \textcolor{comment}{! Sums of all four thicknesses [H]}} \DoxyCodeLine{564 \textcolor{comment}{! real :: I\_h1 ! The inverse of the a thickness [H-\/1]}} \DoxyCodeLine{565 \textcolor{keywordtype}{ real} :: I\_h12, I\_h23, I\_h34 \textcolor{comment}{! The inverses of sums of two thicknesses [H-\/1]}} \DoxyCodeLine{566 \textcolor{keywordtype}{ real} :: I\_h123, I\_h234 \textcolor{comment}{! The inverse of the sum of three thicknesses [H-\/1]}} \DoxyCodeLine{567 \textcolor{keywordtype}{ real} :: I\_h1234 \textcolor{comment}{! The inverse of the sum of all four thicknesses [H-\/1]}} \DoxyCodeLine{568 \textcolor{keywordtype}{ real} :: I\_denom \textcolor{comment}{! The inverse of the denominator some expressions [H-\/3]}} \DoxyCodeLine{569 \textcolor{keywordtype}{ real} :: I\_denB3 \textcolor{comment}{! The inverse of the product of three sums of thicknesses [H-\/3]}} \DoxyCodeLine{570 \textcolor{keywordtype}{ real} :: min\_frac = 1.0e-\/6 \textcolor{comment}{! The square of min\_frac should be much larger than roundoff [nondim]}} \DoxyCodeLine{571 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{parameter} :: C1\_3 = 1.0 / 3.0} \DoxyCodeLine{572 \textcolor{keywordtype}{integer} :: i, j, k} \DoxyCodeLine{573 } \DoxyCodeLine{574 \textcolor{comment}{! These are only used for code verification}} \DoxyCodeLine{575 \textcolor{comment}{! real, dimension(4) :: Atest ! The coefficients of an expression that is being tested.}} \DoxyCodeLine{576 \textcolor{comment}{! real :: zavg, u\_mag, c\_mag}} \DoxyCodeLine{577 \textcolor{comment}{! character(len=128) :: mesg}} \DoxyCodeLine{578 \textcolor{comment}{! real, parameter :: C1\_12 = 1.0 / 12.0}} \DoxyCodeLine{579 } \DoxyCodeLine{580 \textcolor{comment}{! if ((dz(1) == dz(2)) .and. (dz(1) == dz(3)) .and. (dz(1) == dz(4))) then}} \DoxyCodeLine{581 \textcolor{comment}{! ! There are simple closed-\/form expressions in this case}} \DoxyCodeLine{582 \textcolor{comment}{! I\_h1 = 0.0 ; if (dz(1) > 0.0) I\_h1 = 1.0 / dz(1)}} \DoxyCodeLine{583 \textcolor{comment}{! Csys(1) = u(1) + (-\/13.0 * (u(2)-\/u(1)) + 10.0 * (u(3)-\/u(2)) -\/ 3.0 * (u(4)-\/u(3))) * (0.25*C1\_3)}} \DoxyCodeLine{584 \textcolor{comment}{! Csys(2) = (35.0 * (u(2)-\/u(1)) -\/ 34.0 * (u(3)-\/u(2)) + 11.0 * (u(4)-\/u(3))) * (0.25*C1\_3 * I\_h1)}} \DoxyCodeLine{585 \textcolor{comment}{! Csys(3) = (-\/5.0 * (u(2)-\/u(1)) + 8.0 * (u(3)-\/u(2)) -\/ 3.0 * (u(4)-\/u(3))) * (0.25 * I\_h1**2)}} \DoxyCodeLine{586 \textcolor{comment}{! Csys(4) = ((u(2)-\/u(1)) -\/ 2.0 * (u(3)-\/u(2)) + (u(4)-\/u(3))) * (0.5*C1\_3)}} \DoxyCodeLine{587 \textcolor{comment}{! else}} \DoxyCodeLine{588 } \DoxyCodeLine{589 \textcolor{comment}{! Express the coefficients as sums of the differences between properties of succesive layers.}} \DoxyCodeLine{590 } \DoxyCodeLine{591 h1 = dz(1) ; h2 = dz(2) ; h3 = dz(3) ; h4 = dz(4)} \DoxyCodeLine{592 \textcolor{comment}{! Some of the weights used below are proportional to (h1/(h2+h3))**2 or (h1/(h2+h3))*(h2/(h3+h4))}} \DoxyCodeLine{593 \textcolor{comment}{! so h2 and h3 should be adjusted to ensure that these ratios are not so large that property}} \DoxyCodeLine{594 \textcolor{comment}{! differences at the level of roundoff are amplified to be of order 1.}} \DoxyCodeLine{595 \textcolor{keywordflow}{if} ((h2+h3) < min\_frac*h1) h3 = min\_frac*h1 -\/ h2} \DoxyCodeLine{596 \textcolor{keywordflow}{if} ((h3+h4) < min\_frac*h1) h4 = min\_frac*h1 -\/ h3} \DoxyCodeLine{597 } \DoxyCodeLine{598 h12 = h1+h2 ; h23 = h2+h3 ; h34 = h3+h4} \DoxyCodeLine{599 h123 = h12 + h3 ; h234 = h2 + h34 ; h1234 = h12 + h34} \DoxyCodeLine{600 \textcolor{comment}{! Find 3 reciprocals with a single division for efficiency.}} \DoxyCodeLine{601 i\_denb3 = 1.0 / (h123 * h12 * h23)} \DoxyCodeLine{602 i\_h12 = (h123 * h23) * i\_denb3} \DoxyCodeLine{603 i\_h23 = (h12 * h123) * i\_denb3} \DoxyCodeLine{604 i\_h123 = (h12 * h23) * i\_denb3} \DoxyCodeLine{605 i\_denom = 1.0 / ( h1234 * (h234 * h34) )} \DoxyCodeLine{606 i\_h34 = (h1234 * h234) * i\_denom} \DoxyCodeLine{607 i\_h234 = (h1234 * h34) * i\_denom} \DoxyCodeLine{608 i\_h1234 = (h234 * h34) * i\_denom} \DoxyCodeLine{609 } \DoxyCodeLine{610 \textcolor{comment}{! Calculation coefficients in the four equations}} \DoxyCodeLine{611 } \DoxyCodeLine{612 \textcolor{comment}{! The expressions for Csys(3) and Csys(4) come from reducing the 4x4 matrix problem into the following 2x2}} \DoxyCodeLine{613 \textcolor{comment}{! matrix problem, then manipulating the analytic solution to avoid any subtraction and simplifying.}} \DoxyCodeLine{614 \textcolor{comment}{! (C1\_3 * h123 * h23) * Csys(3) + (0.25 * h123 * h23 * (h3 + 2.0*h2 + 3.0*h1)) * Csys(4) =}} \DoxyCodeLine{615 \textcolor{comment}{! (u(3)-\/u(1)) -\/ (u(2)-\/u(1)) * (h12 + h23) * I\_h12}} \DoxyCodeLine{616 \textcolor{comment}{! (C1\_3 * ((h23 + h34) * h1234 + h23 * h3)) * Csys(3) +}} \DoxyCodeLine{617 \textcolor{comment}{! (0.25 * ((h1234 + h123 + h12 + h1) * h23 * h3 + (h1234 + h12 + h1) * (h23 + h34) * h1234)) * Csys(4) =}} \DoxyCodeLine{618 \textcolor{comment}{! (u(4)-\/u(1)) -\/ (u(2)-\/u(1)) * (h123 + h234) * I\_h12}} \DoxyCodeLine{619 \textcolor{comment}{! The final expressions for Csys(1) and Csys(2) were derived by algebraically manipulating the following expressions:}} \DoxyCodeLine{620 \textcolor{comment}{! Csys(1) = (C1\_3 * h1 * h12 * Csys(3) + 0.25 * h1 * h12 * (2.0*h1+h2) * Csys(4)) -\/ \&}} \DoxyCodeLine{621 \textcolor{comment}{! (h1*I\_h12)*(u(2)-\/u(1)) + u(1)}} \DoxyCodeLine{622 \textcolor{comment}{! Csys(2) = (-\/2.0*C1\_3 * (2.0*h1+h2) * Csys(3) -\/ 0.5 * (h1**2 + h12 * (2.0*h1+h2)) * Csys(4)) + \&}} \DoxyCodeLine{623 \textcolor{comment}{! 2.0*I\_h12 * (u(2)-\/u(1))}} \DoxyCodeLine{624 \textcolor{comment}{! These expressions are typically evaluated at x=0 and x=h1, so it is important that these are well behaved}} \DoxyCodeLine{625 \textcolor{comment}{! for these values, suggesting that h1/h23 and h1/h34 should not be allowed to be too large.}} \DoxyCodeLine{626 } \DoxyCodeLine{627 wt(1,1) = -\/h1 * (i\_h1234 + i\_h123 + i\_h12) \textcolor{comment}{! > -\/3}} \DoxyCodeLine{628 wt(2,1) = h1 * h12 * ( i\_h234 * i\_h1234 + i\_h23 * (i\_h234 + i\_h123) ) \textcolor{comment}{! < (h1/h234) + (h1/h23)*(2+(h1/h234))}} \DoxyCodeLine{629 wt(3,1) = -\/h1 * h12 * h123 * i\_denom \textcolor{comment}{! > -\/(h1/h34)*(1+(h1/h234))}} \DoxyCodeLine{630 } \DoxyCodeLine{631 wt(1,2) = 2.0 * (i\_h12*(1.0 + (h1+h12) * (i\_h1234 + i\_h123)) + h1 * i\_h1234*i\_h123) \textcolor{comment}{! < 10/h12}} \DoxyCodeLine{632 wt(2,2) = -\/2.0 * ((h1 * h12 * i\_h1234) * (i\_h23 * (i\_h234 + i\_h123)) + \& \textcolor{comment}{! > -\/(10+6*(h1/h234))/h23}} \DoxyCodeLine{633 (h1+h12) * ( i\_h1234*i\_h234 + i\_h23 * (i\_h234 + i\_h123) ) )} \DoxyCodeLine{634 wt(3,2) = 2.0 * ((h1+h12) * h123 + h1*h12 ) * i\_denom \textcolor{comment}{! < (2+(6*h1/h234)) / h34}} \DoxyCodeLine{635 } \DoxyCodeLine{636 wt(1,3) = -\/3.0 * i\_h12 * i\_h123* ( 1.0 + i\_h1234 * ((h1+h12)+h123) ) \textcolor{comment}{! > -\/12 / (h12*h123)}} \DoxyCodeLine{637 wt(2,3) = 3.0 * i\_h23 * ( i\_h123 + i\_h1234 * ((h1+h12)+h123) * (i\_h123 + i\_h234) ) \textcolor{comment}{! < 12 / (h23\string^2)}} \DoxyCodeLine{638 wt(3,3) = -\/3.0 * ((h1+h12)+h123) * i\_denom \textcolor{comment}{! > -\/9 / (h234*h23)}} \DoxyCodeLine{639 } \DoxyCodeLine{640 wt(1,4) = 4.0 * i\_h1234 * i\_h123 * i\_h12 \textcolor{comment}{! Wt*h1\string^3 < 4}} \DoxyCodeLine{641 wt(2,4) = -\/4.0 * i\_h1234 * (i\_h23 * (i\_h123 + i\_h234)) \textcolor{comment}{! Wt*h1\string^3 > -\/4* (h1/h23)*(1+h1/h234)}} \DoxyCodeLine{642 wt(3,4) = 4.0 * i\_denom \textcolor{comment}{! = 4.0*I\_h1234 * I\_h234 * I\_h34 ! Wt*h1\string^3 < 4 * (h1/h234)*(h1/h34)}} \DoxyCodeLine{643 } \DoxyCodeLine{644 csys(1) = ((u(1) + wt(1,1) * (u(2)-\/u(1))) + wt(2,1) * (u(3)-\/u(2))) + wt(3,1) * (u(4)-\/u(3))} \DoxyCodeLine{645 csys(2) = (wt(1,2) * (u(2)-\/u(1)) + wt(2,2) * (u(3)-\/u(2))) + wt(3,2) * (u(4)-\/u(3))} \DoxyCodeLine{646 csys(3) = (wt(1,3) * (u(2)-\/u(1)) + wt(2,3) * (u(3)-\/u(2))) + wt(3,3) * (u(4)-\/u(3))} \DoxyCodeLine{647 csys(4) = (wt(1,4) * (u(2)-\/u(1)) + wt(2,4) * (u(3)-\/u(2))) + wt(3,4) * (u(4)-\/u(3))} \DoxyCodeLine{648 } \DoxyCodeLine{649 \textcolor{comment}{! endif ! End of non-\/uniform layer thickness branch.}} \DoxyCodeLine{650 } \DoxyCodeLine{651 \textcolor{comment}{! To verify that these answers are correct, uncomment the following:}} \DoxyCodeLine{652 \textcolor{comment}{! u\_mag = 0.0 ; do i=1,4 ; u\_mag = max(u\_mag, abs(u(i))) ; enddo}} \DoxyCodeLine{653 \textcolor{comment}{! do i = 1,4}} \DoxyCodeLine{654 \textcolor{comment}{! if (i==1) then ; zavg = 0.5*dz(i) ; else ; zavg = zavg + 0.5*(dz(i-\/1)+dz(i)) ; endif}} \DoxyCodeLine{655 \textcolor{comment}{! Atest(1) = 1.0}} \DoxyCodeLine{656 \textcolor{comment}{! Atest(2) = zavg ! = ( (z(i+1)**2) -\/ (z(i)**2) ) / (2*dz(i))}} \DoxyCodeLine{657 \textcolor{comment}{! Atest(3) = (zavg**2 + 0.25*C1\_3*dz(i)**2) ! = ( (z(i+1)**3) -\/ (z(i)**3) ) / (3*dz(i))}} \DoxyCodeLine{658 \textcolor{comment}{! Atest(4) = zavg * (zavg**2 + 0.25*dz(i)**2) ! = ( (z(i+1)**4) -\/ (z(i)**4) ) / (4*dz(i))}} \DoxyCodeLine{659 \textcolor{comment}{! c\_mag = 1.0 ; do k=0,3 ; do j=1,3 ; c\_mag = c\_mag + abs(Wt(j,k+1) * zavg**k) ; enddo ; enddo}} \DoxyCodeLine{660 \textcolor{comment}{! write(mesg, '("end\_value\_h4 line ", i2, " c\_mag = ", es10.2, " u\_mag = ", es10.2)') i, c\_mag, u\_mag}} \DoxyCodeLine{661 \textcolor{comment}{! call test\_line(mesg, 4, Atest, Csys, u(i), u\_mag*c\_mag, tol=1.0e-\/15)}} \DoxyCodeLine{662 \textcolor{comment}{! enddo}} \DoxyCodeLine{663 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a09a4c9f369606d81580116620975352e}\label{namespaceregrid__edge__values_a09a4c9f369606d81580116620975352e}} \index{regrid\_edge\_values@{regrid\_edge\_values}!linear\_solver@{linear\_solver}} \index{linear\_solver@{linear\_solver}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{linear\_solver()}{linear\_solver()}} {\footnotesize\ttfamily subroutine regrid\+\_\+edge\+\_\+values\+::linear\+\_\+solver (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N, }\item[{real, dimension(n,n), intent(inout)}]{A, }\item[{real, dimension(n), intent(inout)}]{R, }\item[{real, dimension(n), intent(inout)}]{X }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}} Solve the linear system AX = R by Gaussian elimination. This routine uses Gauss\textquotesingle{}s algorithm to transform the system\textquotesingle{}s original matrix into an upper triangular matrix. Back substitution then yields the answer. The matrix A must be square, with the first index varing along the row. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & The size of the system \\ \hline \mbox{\texttt{ in,out}} & {\em a} & The matrix being inverted \mbox{[}nondim\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em r} & system right-\/hand side \mbox{[}A\mbox{]} \\ \hline \mbox{\texttt{ in,out}} & {\em x} & solution vector \mbox{[}A\mbox{]} \\ \hline \end{DoxyParams} Definition at line 1384 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{1384 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< The size of the system}} \DoxyCodeLine{1385 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N,N)}, \textcolor{keywordtype}{intent(inout)} :: A\textcolor{comment}{ !< The matrix being inverted [nondim]}} \DoxyCodeLine{1386 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(inout)} :: R\textcolor{comment}{ !< system right-\/hand side [A]}} \DoxyCodeLine{1387 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(inout)} :: X\textcolor{comment}{ !< solution vector [A]}} \DoxyCodeLine{1388 } \DoxyCodeLine{1389 \textcolor{comment}{! Local variables}} \DoxyCodeLine{1390 \textcolor{keywordtype}{ real} :: factor \textcolor{comment}{! The factor that eliminates the leading nonzero element in a row.}} \DoxyCodeLine{1391 \textcolor{keywordtype}{ real} :: I\_pivot \textcolor{comment}{! The reciprocal of the pivot value [inverse of the input units of a row of A]}} \DoxyCodeLine{1392 \textcolor{keywordtype}{ real} :: swap} \DoxyCodeLine{1393 \textcolor{keywordtype}{integer} :: i, j, k} \DoxyCodeLine{1394 } \DoxyCodeLine{1395 \textcolor{comment}{! Loop on rows to transform the problem into multiplication by an upper-\/right matrix.}} \DoxyCodeLine{1396 \textcolor{keywordflow}{do} i=1,n-\/1} \DoxyCodeLine{1397 \textcolor{comment}{! Seek a pivot for column i starting in row i, and continuing into the remaining rows. If the}} \DoxyCodeLine{1398 \textcolor{comment}{! pivot is in a row other than i, swap them. If no valid pivot is found, i = N+1 after this loop.}} \DoxyCodeLine{1399 \textcolor{keywordflow}{do} k=i,n ; \textcolor{keywordflow}{if} ( abs(a(i,k)) > 0.0 ) \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop to find pivot}} \DoxyCodeLine{1400 \textcolor{keywordflow}{if} ( k > n ) \textcolor{keywordflow}{then} \textcolor{comment}{! No pivot could be found and the system is singular.}} \DoxyCodeLine{1401 \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{' A='},a} \DoxyCodeLine{1402 \textcolor{keyword}{call }mom\_error( fatal, \textcolor{stringliteral}{'The linear system sent to linear\_solver is singular.'} )} \DoxyCodeLine{1403 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{1404 } \DoxyCodeLine{1405 \textcolor{comment}{! If the pivot is in a row that is different than row i, swap those two rows, noting that both}} \DoxyCodeLine{1406 \textcolor{comment}{! rows start with i-\/1 zero values.}} \DoxyCodeLine{1407 \textcolor{keywordflow}{if} ( k /= i ) \textcolor{keywordflow}{then}} \DoxyCodeLine{1408 \textcolor{keywordflow}{do} j=i,n ; swap = a(j,i) ; a(j,i) = a(j,k) ; a(j,k) = swap ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1409 swap = r(i) ; r(i) = r(k) ; r(k) = swap} \DoxyCodeLine{1410 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{1411 } \DoxyCodeLine{1412 \textcolor{comment}{! Transform the pivot to 1 by dividing the entire row (right-\/hand side included) by the pivot}} \DoxyCodeLine{1413 i\_pivot = 1.0 / a(i,i)} \DoxyCodeLine{1414 a(i,i) = 1.0} \DoxyCodeLine{1415 \textcolor{keywordflow}{do} j=i+1,n ; a(j,i) = a(j,i) * i\_pivot ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1416 r(i) = r(i) * i\_pivot} \DoxyCodeLine{1417 } \DoxyCodeLine{1418 \textcolor{comment}{! Put zeros in column for all rows below that contain the pivot (which is row i)}} \DoxyCodeLine{1419 \textcolor{keywordflow}{do} k=i+1,n \textcolor{comment}{! k is the row index}} \DoxyCodeLine{1420 factor = a(i,k)} \DoxyCodeLine{1421 \textcolor{comment}{! A(i,k) = 0.0 ! These elements are not used again, so this line can be skipped for speed.}} \DoxyCodeLine{1422 \textcolor{keywordflow}{do} j=i+1,n ; a(j,k) = a(j,k) -\/ factor * a(j,i) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1423 r(k) = r(k) -\/ factor * r(i)} \DoxyCodeLine{1424 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1425 } \DoxyCodeLine{1426 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on i}} \DoxyCodeLine{1427 } \DoxyCodeLine{1428 \textcolor{comment}{! Solve the system by back substituting into what is now an upper-\/right matrix.}} \DoxyCodeLine{1429 \textcolor{keywordflow}{if} (a(n,n) == 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! No pivot could be found and the system is singular.}} \DoxyCodeLine{1430 \textcolor{comment}{! write(0,*) ' A=',A}} \DoxyCodeLine{1431 \textcolor{keyword}{call }mom\_error( fatal, \textcolor{stringliteral}{'The final pivot in linear\_solver is zero.'} )} \DoxyCodeLine{1432 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{1433 x(n) = r(n) / a(n,n) \textcolor{comment}{! The last row can now be solved trivially.}} \DoxyCodeLine{1434 \textcolor{keywordflow}{do} i=n-\/1,1,-\/1 \textcolor{comment}{! loop on rows, starting from second to last row}} \DoxyCodeLine{1435 x(i) = r(i)} \DoxyCodeLine{1436 \textcolor{keywordflow}{do} j=i+1,n ; x(i) = x(i) -\/ a(j,i) * x(j) ;\textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1437 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1438 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_a9f0ab1f167eb57cbaca668f7a11b31a1}\label{namespaceregrid__edge__values_a9f0ab1f167eb57cbaca668f7a11b31a1}} \index{regrid\_edge\_values@{regrid\_edge\_values}!solve\_diag\_dominant\_tridiag@{solve\_diag\_dominant\_tridiag}} \index{solve\_diag\_dominant\_tridiag@{solve\_diag\_dominant\_tridiag}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{solve\_diag\_dominant\_tridiag()}{solve\_diag\_dominant\_tridiag()}} {\footnotesize\ttfamily subroutine regrid\+\_\+edge\+\_\+values\+::solve\+\_\+diag\+\_\+dominant\+\_\+tridiag (\begin{DoxyParamCaption}\item[{real, dimension(n), intent(in)}]{Al, }\item[{real, dimension(n), intent(in)}]{Ac, }\item[{real, dimension(n), intent(in)}]{Au, }\item[{real, dimension(n), intent(in)}]{R, }\item[{real, dimension(n), intent(out)}]{X, }\item[{integer, intent(in)}]{N }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}} Solve the tridiagonal system AX = R. This routine uses a variant of Thomas\textquotesingle{}s algorithm to solve the tridiagonal system AX = R, in a form that is guaranteed to avoid dividing by a zero pivot. The matrix A is made up of lower (Al) and upper diagonals (Au) and a central diagonal Ad = Ac+\+Al+\+Au, where Al, Au, and Ac are all positive (or negative) definite. However when Ac is smaller than roundoff compared with (Al+\+Au), the answers are prone to inaccuracy. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em n} & The size of the system \\ \hline \mbox{\texttt{ in}} & {\em ac} & Matrix center diagonal offset from Al + Au \\ \hline \mbox{\texttt{ in}} & {\em al} & Matrix lower diagonal \\ \hline \mbox{\texttt{ in}} & {\em au} & Matrix upper diagonal \\ \hline \mbox{\texttt{ in}} & {\em r} & system right-\/hand side \\ \hline \mbox{\texttt{ out}} & {\em x} & solution vector \\ \hline \end{DoxyParams} Definition at line 1342 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{1342 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< The size of the system}} \DoxyCodeLine{1343 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: Ac\textcolor{comment}{ !< Matrix center diagonal offset from Al + Au}} \DoxyCodeLine{1344 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: Al\textcolor{comment}{ !< Matrix lower diagonal}} \DoxyCodeLine{1345 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: Au\textcolor{comment}{ !< Matrix upper diagonal}} \DoxyCodeLine{1346 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: R\textcolor{comment}{ !< system right-\/hand side}} \DoxyCodeLine{1347 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(out)} :: X\textcolor{comment}{ !< solution vector}} \DoxyCodeLine{1348 \textcolor{comment}{! Local variables}} \DoxyCodeLine{1349 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(N)} :: c1 \textcolor{comment}{! Au / pivot for the backward sweep}} \DoxyCodeLine{1350 \textcolor{keywordtype}{ real} :: d1 \textcolor{comment}{! The next value of 1.0 -\/ c1}} \DoxyCodeLine{1351 \textcolor{keywordtype}{ real} :: I\_pivot \textcolor{comment}{! The inverse of the most recent pivot}} \DoxyCodeLine{1352 \textcolor{keywordtype}{ real} :: denom\_t1 \textcolor{comment}{! The first term in the denominator of the inverse of the pivot.}} \DoxyCodeLine{1353 \textcolor{keywordtype}{integer} :: k \textcolor{comment}{! Loop index}} \DoxyCodeLine{1354 } \DoxyCodeLine{1355 \textcolor{comment}{! Factorization and forward sweep, in a form that will never give a division by a}} \DoxyCodeLine{1356 \textcolor{comment}{! zero pivot for positive definite Ac, Al, and Au.}} \DoxyCodeLine{1357 i\_pivot = 1.0 / (ac(1) + au(1))} \DoxyCodeLine{1358 d1 = ac(1) * i\_pivot} \DoxyCodeLine{1359 c1(1) = au(1) * i\_pivot} \DoxyCodeLine{1360 x(1) = r(1) * i\_pivot} \DoxyCodeLine{1361 \textcolor{keywordflow}{do} k=2,n-\/1} \DoxyCodeLine{1362 denom\_t1 = ac(k) + d1 * al(k)} \DoxyCodeLine{1363 i\_pivot = 1.0 / (denom\_t1 + au(k))} \DoxyCodeLine{1364 d1 = denom\_t1 * i\_pivot} \DoxyCodeLine{1365 c1(k) = au(k) * i\_pivot} \DoxyCodeLine{1366 x(k) = (r(k) -\/ al(k) * x(k-\/1)) * i\_pivot} \DoxyCodeLine{1367 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1368 i\_pivot = 1.0 / (ac(n) + d1 * al(n))} \DoxyCodeLine{1369 x(n) = (r(n) -\/ al(n) * x(n-\/1)) * i\_pivot} \DoxyCodeLine{1370 \textcolor{comment}{! Backward sweep}} \DoxyCodeLine{1371 \textcolor{keywordflow}{do} k=n-\/1,1,-\/1} \DoxyCodeLine{1372 x(k) = x(k) -\/ c1(k) * x(k+1)} \DoxyCodeLine{1373 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1374 } \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__edge__values_aef5443593d35636eb95de2b580878a20}\label{namespaceregrid__edge__values_aef5443593d35636eb95de2b580878a20}} \index{regrid\_edge\_values@{regrid\_edge\_values}!test\_line@{test\_line}} \index{test\_line@{test\_line}!regrid\_edge\_values@{regrid\_edge\_values}} \doxysubsubsection{\texorpdfstring{test\_line()}{test\_line()}} {\footnotesize\ttfamily subroutine regrid\+\_\+edge\+\_\+values\+::test\+\_\+line (\begin{DoxyParamCaption}\item[{character(len=$\ast$), intent(in)}]{msg, }\item[{integer, intent(in)}]{N, }\item[{real, dimension(4), intent(in)}]{A, }\item[{real, dimension(4), intent(in)}]{C, }\item[{real, intent(in)}]{R, }\item[{real, intent(in)}]{mag, }\item[{real, intent(in), optional}]{tol }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}} Test that A$\ast$C = R to within a tolerance, issuing a fatal error with an explanatory message if they do not. \begin{DoxyParams}[1]{Parameters} \mbox{\texttt{ in}} & {\em mag} & The magnitude of leading order terms in this line \\ \hline \mbox{\texttt{ in}} & {\em n} & The number of points in the system \\ \hline \mbox{\texttt{ in}} & {\em a} & One of the two vectors being multiplied \\ \hline \mbox{\texttt{ in}} & {\em c} & One of the two vectors being multiplied \\ \hline \mbox{\texttt{ in}} & {\em r} & The expected solution of the equation \\ \hline \mbox{\texttt{ in}} & {\em msg} & An identifying message for this test \\ \hline \mbox{\texttt{ in}} & {\em tol} & The fractional tolerance for the two solutions \\ \hline \end{DoxyParams} Definition at line 1445 of file regrid\+\_\+edge\+\_\+values.\+F90. \begin{DoxyCode}{0} \DoxyCodeLine{1445 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{intent(in)} :: mag\textcolor{comment}{ !< The magnitude of leading order terms in this line}} \DoxyCodeLine{1446 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: N\textcolor{comment}{ !< The number of points in the system}} \DoxyCodeLine{1447 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)}, \textcolor{keywordtype}{intent(in)} :: A\textcolor{comment}{ !< One of the two vectors being multiplied}} \DoxyCodeLine{1448 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{dimension(4)}, \textcolor{keywordtype}{intent(in)} :: C\textcolor{comment}{ !< One of the two vectors being multiplied}} \DoxyCodeLine{1449 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{intent(in)} :: R\textcolor{comment}{ !< The expected solution of the equation}} \DoxyCodeLine{1450 \textcolor{keywordtype}{character(len=*)}, \textcolor{keywordtype}{intent(in)} :: msg\textcolor{comment}{ !< An identifying message for this test}} \DoxyCodeLine{1451 \textcolor{keywordtype}{ real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: tol\textcolor{comment}{ !< The fractional tolerance for the two solutions}} \DoxyCodeLine{1452 } \DoxyCodeLine{1453 \textcolor{keywordtype}{ real} :: sum, sum\_mag} \DoxyCodeLine{1454 \textcolor{keywordtype}{ real} :: tolerance} \DoxyCodeLine{1455 \textcolor{keywordtype}{character(len=128)} :: mesg2} \DoxyCodeLine{1456 \textcolor{keywordtype}{integer} :: i} \DoxyCodeLine{1457 } \DoxyCodeLine{1458 tolerance = 1.0e-\/12 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(tol)) tolerance = tol} \DoxyCodeLine{1459 } \DoxyCodeLine{1460 sum = 0.0 ; sum\_mag = max(0.0,mag)} \DoxyCodeLine{1461 \textcolor{keywordflow}{do} i=1,n} \DoxyCodeLine{1462 sum = sum + a(i) * c(i)} \DoxyCodeLine{1463 sum\_mag = sum\_mag + abs(a(i) * c(i))} \DoxyCodeLine{1464 \textcolor{keywordflow}{ enddo}} \DoxyCodeLine{1465 } \DoxyCodeLine{1466 \textcolor{keywordflow}{if} (abs(sum -\/ r) > tolerance * (sum\_mag + abs(r))) \textcolor{keywordflow}{then}} \DoxyCodeLine{1467 \textcolor{keyword}{write}(mesg2, \textcolor{stringliteral}{'(", Fractional error = ", es12.4,", sum = ", es12.4)'}) (sum -\/ r) / (sum\_mag + abs(r)), sum} \DoxyCodeLine{1468 \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"Failed line test: "}//trim(msg)//trim(mesg2))} \DoxyCodeLine{1469 \textcolor{keywordflow}{ endif}} \DoxyCodeLine{1470 } \end{DoxyCode}