\hypertarget{namespaceregrid__edge__values}{}\section{regrid\+\_\+edge\+\_\+values Module Reference}
\label{namespaceregrid__edge__values}\index{regrid\+\_\+edge\+\_\+values@{regrid\+\_\+edge\+\_\+values}}


\subsection{Detailed Description}
Edge value estimation for high-\/order resconstruction. \subsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \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 \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 \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 \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 \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 \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 \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 \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 \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 \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 \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 \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 \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}
\subsection*{Variables}
\begin{DoxyCompactItemize}
\item 
\mbox{\Hypertarget{namespaceregrid__edge__values_a241fc99294d87afd8146f4b3e48dc021}\label{namespaceregrid__edge__values_a241fc99294d87afd8146f4b3e48dc021}} 
real, parameter \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 \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 \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}


\subsection{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}}
\subsubsection{\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{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt 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}
116   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of cells}
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}
118 \textcolor{comment}{                                           !! second index is for the two edges of each cell.}
119   \textcolor{comment}{! Local variables}
120   \textcolor{keywordtype}{integer}       :: k            \textcolor{comment}{! loop index}
121   \textcolor{keywordtype}{real}          :: u0\_avg       \textcolor{comment}{! avg value at given edge}
122 
123   \textcolor{comment}{! Loop on interior edges}
124   \textcolor{keywordflow}{do} k = 1,n-1
125     \textcolor{comment}{! Compare edge values on the right and left sides of the edge}
126     \textcolor{keywordflow}{if} ( edge\_val(k,2) /= edge\_val(k+1,1) ) \textcolor{keywordflow}{then}
127       u0\_avg = 0.5 * ( edge\_val(k,2) + edge\_val(k+1,1) )
128       edge\_val(k,2) = u0\_avg
129       edge\_val(k+1,1) = u0\_avg
130 \textcolor{keywordflow}{    endif}
131 
132 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on interior edges}
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}}
\subsubsection{\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{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt 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{\tt in}  & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt 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}
45   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of cells}
46   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths [H]}
47   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}
48   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Potentially modified edge values [A]; the}
49 \textcolor{comment}{                                           !! second index is for the two edges of each cell.}
50   \textcolor{keywordtype}{real},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}
51   \textcolor{keywordtype}{logical},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}
52   \textcolor{comment}{! Local variables}
53   \textcolor{keywordtype}{real}    :: sigma\_l, sigma\_c, sigma\_r    \textcolor{comment}{! left, center and right van Leer slopes [A H-1] or [A]}
54   \textcolor{keywordtype}{real}    :: slope\_x\_h     \textcolor{comment}{! retained PLM slope times  half grid step [A]}
55   \textcolor{keywordtype}{real}    :: hneglect      \textcolor{comment}{! A negligible thickness [H].}
56   \textcolor{keywordtype}{logical} :: use\_2018\_answers  \textcolor{comment}{! If true use older, less acccurate expressions.}
57   \textcolor{keywordtype}{integer} :: k, km1, kp1   \textcolor{comment}{! Loop index and the values to either side.}
58 
59   use\_2018\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) use\_2018\_answers = answers\_2018
60   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
61     hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect
62 \textcolor{keywordflow}{  endif}
63 
64   \textcolor{comment}{! Loop on cells to bound edge value}
65   \textcolor{keywordflow}{do} k = 1,n
66 
67     \textcolor{comment}{! For the sake of bounding boundary edge values, the left neighbor of the left boundary cell}
68     \textcolor{comment}{! is assumed to be the same as the left boundary cell and the right neighbor of the right}
69     \textcolor{comment}{! boundary cell is assumed to be the same as the right boundary cell. This effectively makes}
70     \textcolor{comment}{! boundary cells look like extrema.}
71     km1 = max(1,k-1) ; kp1 = min(k+1,n)
72 
73     slope\_x\_h = 0.0
74     \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
75       sigma\_l = 2.0 * ( u(k) - u(km1) ) / ( h(k) + hneglect )
76       sigma\_c = 2.0 * ( u(kp1) - u(km1) ) / ( h(km1) + 2.0*h(k) + h(kp1) + hneglect )
77       sigma\_r = 2.0 * ( u(kp1) - u(k) ) / ( h(k) + hneglect )
78 
79       \textcolor{comment}{! The limiter is used in the local coordinate system to each cell, so for convenience store}
80       \textcolor{comment}{! the slope times a half grid spacing.  (See White and Adcroft JCP 2008 Eqs 19 and 20)}
81       \textcolor{keywordflow}{if} ( (sigma\_l * sigma\_r) > 0.0 ) &
82         slope\_x\_h = 0.5 * h(k) * sign( min(abs(sigma\_l),abs(sigma\_c),abs(sigma\_r)), sigma\_c )
83     \textcolor{keywordflow}{elseif} ( ((h(km1) + h(kp1)) + 2.0*h(k)) > 0.0 ) \textcolor{keywordflow}{then}
84       sigma\_l = ( u(k) - u(km1) )
85       sigma\_c = ( u(kp1) - u(km1) ) * ( h(k) / ((h(km1) + h(kp1)) + 2.0*h(k)) )
86       sigma\_r = ( u(kp1) - u(k) )
87 
88       \textcolor{comment}{! The limiter is used in the local coordinate system to each cell, so for convenience store}
89       \textcolor{comment}{! the slope times a half grid spacing.  (See White and Adcroft JCP 2008 Eqs 19 and 20)}
90       \textcolor{keywordflow}{if} ( (sigma\_l * sigma\_r) > 0.0 ) &
91         slope\_x\_h = sign( min(abs(sigma\_l),abs(sigma\_c),abs(sigma\_r)), sigma\_c )
92 \textcolor{keywordflow}{    endif}
93 
94     \textcolor{comment}{! Limit the edge values}
95     \textcolor{keywordflow}{if} ( (u(km1)-edge\_val(k,1)) * (edge\_val(k,1)-u(k)) < 0.0 ) \textcolor{keywordflow}{then}
96       edge\_val(k,1) = u(k) - sign( min( abs(slope\_x\_h), abs(edge\_val(k,1)-u(k)) ), slope\_x\_h )
97 \textcolor{keywordflow}{    endif}
98 
99     \textcolor{keywordflow}{if} ( (u(kp1)-edge\_val(k,2)) * (edge\_val(k,2)-u(k)) < 0.0 ) \textcolor{keywordflow}{then}
100       edge\_val(k,2) = u(k) + sign( min( abs(slope\_x\_h), abs(edge\_val(k,2)-u(k)) ), slope\_x\_h )
101 \textcolor{keywordflow}{    endif}
102 
103     \textcolor{comment}{! Finally bound by neighboring cell means in case of roundoff}
104     edge\_val(k,1) = max( min( edge\_val(k,1), max(u(km1), u(k)) ), min(u(km1), u(k)) )
105     edge\_val(k,2) = max( min( edge\_val(k,2), max(u(kp1), u(k)) ), min(u(kp1), u(k)) )
106 
107 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! loop on interior edges}
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}}
\subsubsection{\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{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em u} & cell averages in arbitrary units \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt 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}
141   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of cells}
142   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell averages in arbitrary units [A]}
143   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Cell edge values [A]; the}
144 \textcolor{comment}{                                           !! second index is for the two edges of each cell.}
145   \textcolor{comment}{! Local variables}
146   \textcolor{keywordtype}{integer}       :: k            \textcolor{comment}{! loop index}
147   \textcolor{keywordtype}{real}          :: u0\_avg       \textcolor{comment}{! avg value at given edge [A]}
148 
149   \textcolor{keywordflow}{do} k = 1,n-1
150     \textcolor{keywordflow}{if} ( (edge\_val(k+1,1) - edge\_val(k,2)) * (u(k+1) - u(k)) < 0.0 ) \textcolor{keywordflow}{then}
151       u0\_avg = 0.5 * ( edge\_val(k,2) + edge\_val(k+1,1) )
152       u0\_avg = max( min( u0\_avg, max(u(k), u(k+1)) ), min(u(k), u(k+1)) )
153       edge\_val(k,2) = u0\_avg
154       edge\_val(k+1,1) = u0\_avg
155 \textcolor{keywordflow}{    endif}
156 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on interior edges}
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}}
\subsubsection{\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'_{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-\/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{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt 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{\tt in}  & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt 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}
697   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of cells}
698   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths [H]}
699   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}
700   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_slopes\textcolor{comment}{ !< Returned edge slopes [A H-1]; the}
701 \textcolor{comment}{                                           !! second index is for the two edges of each cell.}
702   \textcolor{keywordtype}{real},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}
703   \textcolor{keywordtype}{logical},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}
704   \textcolor{comment}{! Local variables}
705   \textcolor{keywordtype}{integer}               :: i, j                 \textcolor{comment}{! loop indexes}
706   \textcolor{keywordtype}{real}                  :: h0, h1               \textcolor{comment}{! cell widths [H or nondim]}
707   \textcolor{keywordtype}{real}                  :: h0\_2, h1\_2, h0h1     \textcolor{comment}{! products of cell widths [H2 or nondim]}
708   \textcolor{keywordtype}{real}                  :: h0\_3, h1\_3           \textcolor{comment}{! products of three cell widths [H3 or nondim]}
709   \textcolor{keywordtype}{real}                  :: h\_min                \textcolor{comment}{! A minimal cell width [H]}
710   \textcolor{keywordtype}{real}                  :: d                    \textcolor{comment}{! A temporary variable [H3]}
711   \textcolor{keywordtype}{real}                  :: i\_d                  \textcolor{comment}{! A temporary variable [nondim]}
712   \textcolor{keywordtype}{real}                  :: i\_h                  \textcolor{comment}{! Inverses of thicknesses [H-1]}
713   \textcolor{keywordtype}{real}                  :: alpha, beta          \textcolor{comment}{! stencil coefficients [nondim]}
714   \textcolor{keywordtype}{real}                  :: a, b                 \textcolor{comment}{! weights of cells [H-1]}
715   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}       :: c1\_12 = 1.0 / 12.0
716   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}    :: dz                   \textcolor{comment}{! A temporary array of limited layer thicknesses [H]}
717   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}    :: u\_tmp                \textcolor{comment}{! A temporary array of cell average properties [A]}
718   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(5)}    :: x          \textcolor{comment}{! Coordinate system with 0 at edges [H]}
719   \textcolor{keywordtype}{real}                  :: dx, xavg   \textcolor{comment}{! Differences and averages of successive values of x [H]}
720   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4,4)}  :: asys       \textcolor{comment}{! matrix used to find boundary conditions}
721   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}    :: bsys, csys
722   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(3)}    :: dsys
723   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N+1)}  :: tri\_l, &     \textcolor{comment}{! tridiagonal system (lower diagonal) [nondim]}
724                            tri\_d, &     \textcolor{comment}{! tridiagonal system (middle diagonal) [nondim]}
725                            tri\_c, &     \textcolor{comment}{! tridiagonal system central value, with tri\_d = tri\_c+tri\_l+tri\_u}
726                            tri\_u, &     \textcolor{comment}{! tridiagonal system (upper diagonal) [nondim]}
727                            tri\_b, &     \textcolor{comment}{! tridiagonal system (right hand side) [A H-1]}
728                            tri\_x        \textcolor{comment}{! tridiagonal system (solution vector) [A H-1]}
729   \textcolor{keywordtype}{real}      :: hneglect  \textcolor{comment}{! A negligible thickness [H].}
730   \textcolor{keywordtype}{real}      :: hneglect3 \textcolor{comment}{! hNeglect^3 [H3].}
731   \textcolor{keywordtype}{logical}   :: use\_2018\_answers  \textcolor{comment}{! If true use older, less acccurate expressions.}
732 
733   hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect))  hneglect = h\_neglect
734   hneglect3 = hneglect**3
735   use\_2018\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) use\_2018\_answers = answers\_2018
736 
737   \textcolor{comment}{! Loop on cells (except last one)}
738   \textcolor{keywordflow}{do} i = 1,n-1
739 
740     \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
741       \textcolor{comment}{! Get cell widths}
742       h0 = h(i)
743       h1 = h(i+1)
744 
745       \textcolor{comment}{! Auxiliary calculations}
746       h0h1 = h0 * h1
747       h0\_2 = h0 * h0
748       h1\_2 = h1 * h1
749       h0\_3 = h0\_2 * h0
750       h1\_3 = h1\_2 * h1
751 
752       d = 4.0 * h0h1 * ( h0 + h1 ) + h1\_3 + h0\_3
753 
754       \textcolor{comment}{! Coefficients}
755       alpha = h1 * (h0\_2 + h0h1 - h1\_2) / ( d + hneglect3 )
756       beta  = h0 * (h1\_2 + h0h1 - h0\_2) / ( d + hneglect3 )
757       a = -12.0 * h0h1 / ( d + hneglect3 )
758       b = -a
759 
760       tri\_l(i+1) = alpha
761       tri\_d(i+1) = 1.0
762       tri\_u(i+1) = beta
763 
764       tri\_b(i+1) = a * u(i) + b * u(i+1)
765     \textcolor{keywordflow}{else}
766       \textcolor{comment}{! Get cell widths}
767       h0 = max(h(i), hneglect)
768       h1 = max(h(i+1), hneglect)
769 
770       i\_h = 1.0 / (h0 + h1)
771       h0 = h0 * i\_h ; h1 = h1 * i\_h
772 
773       h0h1 = h0 * h1 ; h0\_2 = h0 * h0 ; h1\_2 = h1 * h1
774       h0\_3 = h0\_2 * h0 ; h1\_3 = h1\_2 * h1
775 
776       \textcolor{comment}{! Set the tridiagonal coefficients}
777       i\_d = 1.0 / (4.0 * h0h1 * ( h0 + h1 ) + h1\_3 + h0\_3) \textcolor{comment}{! = 1 / ((h0 + h1)**3 + h0*h1*(h0 + h1))}
778       tri\_l(i+1) = (h1 * ((h0\_2 + h0h1) - h1\_2)) * i\_d
779       \textcolor{comment}{! tri\_d(i+1) = 1.0}
780       tri\_c(i+1) = 2.0 * ((h0\_2 + h1\_2) * (h0 + h1)) * i\_d
781       tri\_u(i+1) = (h0 * ((h1\_2 + h0h1) - h0\_2)) * i\_d
782       \textcolor{comment}{! The following expressions have been simplified using the nondimensionalization above:}
783       \textcolor{comment}{! I\_d = 1.0 / (1.0 + h0h1)}
784       \textcolor{comment}{! tri\_l(i+1) = (h0h1 - h1\_3) * I\_d}
785       \textcolor{comment}{! tri\_c(i+1) = 2.0 * (h0\_2 + h1\_2) * I\_d}
786       \textcolor{comment}{! tri\_u(i+1) = (h0h1 - h0\_3) * I\_d}
787 
788       tri\_b(i+1) = 12.0 * (h0h1 * i\_d) * ((u(i+1) - u(i)) * i\_h)
789 \textcolor{keywordflow}{    endif}
790 
791 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on cells}
792 
793   \textcolor{comment}{! Boundary conditions: set the first edge slope}
794   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
795     x(1) = 0.0
796     \textcolor{keywordflow}{do} i = 1,4
797       dx = h(i)
798       x(i+1) = x(i) + dx
799       \textcolor{keywordflow}{do} j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ;\textcolor{keywordflow}{ enddo}
800       bsys(i) = u(i) * dx
801 \textcolor{keywordflow}{    enddo}
802 
803     \textcolor{keyword}{call }solve\_linear\_system( asys, bsys, csys, 4 )
804 
805     dsys(1) = csys(2) ; dsys(2) = 2.0 * csys(3) ; dsys(3) = 3.0 * csys(4)
806     tri\_b(1) = evaluation\_polynomial( dsys, 3, x(1) )  \textcolor{comment}{! Set the first edge slope}
807     tri\_d(1) = 1.0
808   \textcolor{keywordflow}{else} \textcolor{comment}{! Use expressions with less sensitivity to roundoff}
809     \textcolor{keywordflow}{do} i=1,4 ; dz(i) = max(hneglect, h(i) ) ; u\_tmp(i) = u(i) ;\textcolor{keywordflow}{ enddo}
810     \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, csys)
811 
812     \textcolor{comment}{! Set the first edge slope}
813     tri\_b(1) = csys(2)
814     tri\_c(1) = 1.0
815 \textcolor{keywordflow}{  endif}
816   tri\_u(1) = 0.0 \textcolor{comment}{! tri\_l(1) = 0.0}
817 
818   \textcolor{comment}{! Boundary conditions: set the last edge slope}
819   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
820     x(1) = 0.0
821     \textcolor{keywordflow}{do} i = 1,4
822       dx = h(n-4+i)
823       x(i+1) = x(i) + dx
824       \textcolor{keywordflow}{do} j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ;\textcolor{keywordflow}{ enddo}
825       bsys(i) = u(n-4+i) * dx
826 \textcolor{keywordflow}{    enddo}
827 
828     \textcolor{keyword}{call }solve\_linear\_system( asys, bsys, csys, 4 )
829 
830     dsys(1) = csys(2) ; dsys(2) = 2.0 * csys(3) ; dsys(3) = 3.0 * csys(4)
831     \textcolor{comment}{! Set the last edge slope}
832     tri\_b(n+1) = evaluation\_polynomial( dsys, 3, x(5) )
833     tri\_d(n+1) = 1.0
834   \textcolor{keywordflow}{else}
835     \textcolor{comment}{! Use expressions with less sensitivity to roundoff, including using a coordinate}
836     \textcolor{comment}{! system that sets the origin at the last interface in the domain.}
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}
838 
839     \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, csys)
840 
841     \textcolor{comment}{! Set the last edge slope}
842     tri\_b(n+1) = -csys(2)
843     tri\_c(n+1) = 1.0
844 \textcolor{keywordflow}{  endif}
845   tri\_l(n+1) = 0.0 \textcolor{comment}{! tri\_u(N+1) = 0.0}
846 
847   \textcolor{comment}{! Solve tridiagonal system and assign edge slopes}
848   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
849     \textcolor{keyword}{call }solve\_tridiagonal\_system( tri\_l, tri\_d, tri\_u, tri\_b, tri\_x, n+1 )
850   \textcolor{keywordflow}{else}
851     \textcolor{keyword}{call }solve\_diag\_dominant\_tridiag( tri\_l, tri\_c, tri\_u, tri\_b, tri\_x, n+1 )
852 \textcolor{keywordflow}{  endif}
853 
854   \textcolor{keywordflow}{do} i = 2,n
855     edge\_slopes(i,1)   = tri\_x(i)
856     edge\_slopes(i-1,2) = tri\_x(i)
857 \textcolor{keywordflow}{  enddo}
858   edge\_slopes(1,1) = tri\_x(1)
859   edge\_slopes(n,2) = tri\_x(n+1)
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}}
\subsubsection{\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{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt 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{\tt in}  & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt 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}
867   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of cells}
868   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths [H]}
869   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}
870   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_slopes\textcolor{comment}{ !< Returned edge slopes [A H-1]; the}
871 \textcolor{comment}{                                           !! second index is for the two edges of each cell.}
872   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional},       \textcolor{keywordtype}{intent(in)}    :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}
873   \textcolor{keywordtype}{logical},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}
874 \textcolor{comment}{! -----------------------------------------------------------------------------}
875 \textcolor{comment}{! Fifth-order implicit estimates of edge slopes are based on a four-cell,}
876 \textcolor{comment}{! three-edge stencil. A tridiagonal system is set up and is based on}
877 \textcolor{comment}{! expressing the edge slopes in terms of neighboring cell averages.}
878 \textcolor{comment}{!}
879 \textcolor{comment}{! The generic relationship is}
880 \textcolor{comment}{!}
881 \textcolor{comment}{! \(\backslash\)alpha u'\_\{i-1/2\} + u'\_\{i+1/2\} + \(\backslash\)beta u'\_\{i+3/2\} =}
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\}}
883 \textcolor{comment}{!}
884 \textcolor{comment}{! and the stencil looks like this}
885 \textcolor{comment}{!}
886 \textcolor{comment}{!         i-1     i     i+1    i+2}
887 \textcolor{comment}{!   ..--o------o------o------o------o--..}
888 \textcolor{comment}{!            i-1/2  i+1/2  i+3/2}
889 \textcolor{comment}{!}
890 \textcolor{comment}{! In this routine, the coefficients \(\backslash\)alpha, \(\backslash\)beta, a, b, c and d are}
891 \textcolor{comment}{! computed, the tridiagonal system is built, boundary conditions are}
892 \textcolor{comment}{! prescribed and the system is solved to yield edge-value estimates.}
893 \textcolor{comment}{!}
894 \textcolor{comment}{! Note that the centered stencil only applies to edges 3 to N-1 (edges are}
895 \textcolor{comment}{! numbered 1 to n+1), which yields N-3 equations for N+1 unknowns. Two other}
896 \textcolor{comment}{! equations are written by using a right-biased stencil for edge 2 and a}
897 \textcolor{comment}{! left-biased stencil for edge N. The prescription of boundary conditions}
898 \textcolor{comment}{! (using sixth-order polynomials) closes the system.}
899 \textcolor{comment}{!}
900 \textcolor{comment}{! CAUTION: For each edge, in order to determine the coefficients of the}
901 \textcolor{comment}{!          implicit expression, a 6x6 linear system is solved. This may}
902 \textcolor{comment}{!          become computationally expensive if regridding is carried out}
903 \textcolor{comment}{!          often. Figuring out closed-form expressions for these coefficients}
904 \textcolor{comment}{!          on nonuniform meshes turned out to be intractable.}
905 \textcolor{comment}{! -----------------------------------------------------------------------------}
906 
907   \textcolor{comment}{! Local variables}
908   \textcolor{keywordtype}{real} :: h0, h1, h2, h3       \textcolor{comment}{! cell widths [H]}
909   \textcolor{keywordtype}{real} :: hmin                 \textcolor{comment}{! The minimum thickness used in these calculations [H]}
910   \textcolor{keywordtype}{real} :: h01, h01\_2           \textcolor{comment}{! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n]}
911   \textcolor{keywordtype}{real} :: h23, h23\_2           \textcolor{comment}{! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n]}
912   \textcolor{keywordtype}{real} :: hneglect             \textcolor{comment}{! A negligible thickness [H].}
913   \textcolor{keywordtype}{real}                  :: h1\_2, h2\_2           \textcolor{comment}{! the coefficients of the}
914   \textcolor{keywordtype}{real}                  :: h1\_3, h2\_3           \textcolor{comment}{! tridiagonal system}
915   \textcolor{keywordtype}{real}                  :: h1\_4, h2\_4           \textcolor{comment}{! ...}
916   \textcolor{keywordtype}{real}                  :: h1\_5, h2\_5           \textcolor{comment}{! ...}
917   \textcolor{keywordtype}{real}                  :: alpha, beta          \textcolor{comment}{! stencil coefficients}
918   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(7)}    :: x                    \textcolor{comment}{! Coordinate system with 0 at edges [same units as h]}
919   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}       :: c1\_12 = 1.0 / 12.0
920   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}       :: c5\_6 = 5.0 / 6.0
921   \textcolor{keywordtype}{real}                  :: dx, xavg             \textcolor{comment}{! Differences and averages of successive values of x [same
       units as h]}
922   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(6,6)}  :: asys                 \textcolor{comment}{! matrix used to find  boundary conditions}
923   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(6)}    :: bsys, csys           \textcolor{comment}{! ...}
924   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N+1)}  :: tri\_l, &             \textcolor{comment}{! trid. system (lower diagonal)}
925                            tri\_d, &             \textcolor{comment}{! trid. system (middle diagonal)}
926                            tri\_u, &             \textcolor{comment}{! trid. system (upper diagonal)}
927                            tri\_b, &             \textcolor{comment}{! trid. system (unknowns vector)}
928                            tri\_x                \textcolor{comment}{! trid. system (rhs)}
929   \textcolor{keywordtype}{real} :: h\_min\_frac = 1.0e-4
930   \textcolor{keywordtype}{integer} :: i, j, k              \textcolor{comment}{! loop indexes}
931 
932   hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect
933 
934   \textcolor{comment}{! Loop on cells (except the first and last ones)}
935   \textcolor{keywordflow}{do} k = 2,n-2
936     \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}
937     hmin = max(hneglect, h\_min\_frac*((h(k-1) + h(k)) + (h(k+1) + h(k+2))))
938     h0 = max(h(k-1), hmin) ; h1 = max(h(k), hmin)
939     h2 = max(h(k+1), hmin) ; h3 = max(h(k+2), hmin)
940 
941     \textcolor{comment}{! Auxiliary calculations}
942     h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2
943     h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2
944 
945     \textcolor{comment}{! Compute matrix entries as described in Eq. (52) of White and Adcroft (2009).  The last 4 rows are}
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,
       &}
947     \textcolor{comment}{!                  (-h1)**(n-1), h2**(n-1), ((h2+h3)**n - h2**n) / h3  /)}
948 
949     asys(1:6,1) = (/  0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /)
950     asys(1:6,2) = (/  2.0, 2.0, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /)
951     asys(1:6,3) = (/  6.0*h1, -6.0* h2, (3.0*h1\_2 + h0*(3.0*h1 + h0)), &
952                       h1\_2, h2\_2, (3.0*h2\_2 + h3*(3.0*h2 + h3)) /)
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))), &
954                      -h1\_3,       h2\_3,       (4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)
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)))
      ), &
956                       h1\_4,      h2\_4,       (5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))
      ) /)
957     asys(1:6,6) = (/ -30.0*h1\_4, -30.0*h2\_4, &
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))))), &
959                      -h1\_5, h2\_5, &
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))))) /)
961     bsys(1:6) = (/ 0.0, -2.0, 0.0, 0.0, 0.0, 0.0 /)
962 
963     \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )
964 
965     alpha = csys(1)
966     beta  = csys(2)
967 
968     tri\_l(k+1) = alpha
969     tri\_d(k+1) = 1.0
970     tri\_u(k+1) = beta
971     tri\_b(k+1) = csys(3) * u(k-1) + csys(4) * u(k) + csys(5) * u(k+1) + csys(6) * u(k+2)
972 
973 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on cells}
974 
975   \textcolor{comment}{! Use a right-biased stencil for the second row, as described in Eq. (53) of White and Adcroft (2009).}
976 
977   \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}
978   hmin = max(hneglect, h\_min\_frac*((h(1) + h(2)) + (h(3) + h(4))))
979   h0 = max(h(1), hmin) ; h1 = max(h(2), hmin)
980   h2 = max(h(3), hmin) ; h3 = max(h(4), hmin)
981 
982   \textcolor{comment}{! Auxiliary calculations}
983   h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2
984   h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2
985   h01 = h0 + h1 ; h01\_2 = h01 * h01
986 
987   \textcolor{comment}{! Compute matrix entries}
988   asys(1:6,1) = (/  0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /)
989   asys(1:6,2) = (/  2.0,     2.0, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /)
990   asys(1:6,3) = (/  6.0*h01, 0.0, (3.0*h1\_2 + h0*(3.0*h1 + h0)), &
991                     h1\_2,   h2\_2, (3.0*h2\_2 + h3*(3.0*h2 + h3)) /)
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))), &
993                    -h1\_3,       h2\_3,  (4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)
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)))),
       &
995                     h1\_4,      h2\_4,       (5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))) 
      /)
996   asys(1:6,6) = (/ -30.0*(h01\_2*h01\_2), 0.0, &
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))))), &
998                    -h1\_5, h2\_5, &
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))))) /)
1000   bsys(1:6) = (/ 0.0, -2.0, -6.0*h1, 12.0*h1\_2, -20.0*h1\_3, 30.0*h1\_4 /)
1001 
1002   \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )
1003 
1004   alpha = csys(1)
1005   beta  = csys(2)
1006 
1007   tri\_l(2) = alpha
1008   tri\_d(2) = 1.0
1009   tri\_u(2) = beta
1010   tri\_b(2) = csys(3) * u(1) + csys(4) * u(2) + csys(5) * u(3) + csys(6) * u(4)
1011 
1012   \textcolor{comment}{! Boundary conditions: left boundary}
1013   x(1) = 0.0
1014   \textcolor{keywordflow}{do} i = 1,6
1015     dx = h(i)
1016     xavg = x(i) + 0.5 * dx
1017     asys(1:6,i) =  (/ 1.0, xavg, (xavg**2 + c1\_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), &
1018                       (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), &
1019                        xavg * (xavg**4 + c5\_6*xavg**2*dx**2 + 0.0625*dx**4) /)
1020     bsys(i) = u(i)
1021     x(i+1) = x(i) + dx
1022 \textcolor{keywordflow}{  enddo}
1023 
1024   \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )
1025 
1026   tri\_d(1) = 0.0
1027   tri\_d(1) = 1.0
1028   tri\_u(1) = 0.0
1029   tri\_b(1) = csys(2) \textcolor{comment}{! first edge value}
1030 
1031   \textcolor{comment}{! Use a left-biased stencil for the second to last row, as described in Eq. (54) of White and Adcroft
       (2009).}
1032 
1033   \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}
1034   hmin = max(hneglect, h\_min\_frac*((h(n-3) + h(n-2)) + (h(n-1) + h(n))))
1035   h0 = max(h(n-3), hmin) ; h1 = max(h(n-2), hmin)
1036   h2 = max(h(n-1), hmin) ; h3 = max(h(n), hmin)
1037 
1038   \textcolor{comment}{! Auxiliary calculations}
1039   h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2
1040   h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2
1041 
1042   h23 = h2 + h3 ; h23\_2 = h23 * h23
1043 
1044   \textcolor{comment}{! Compute matrix entries}
1045   asys(1:6,1) = (/  0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /)
1046   asys(1:6,2) = (/  2.0,     2.0, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /)
1047   asys(1:6,3) = (/  0.0, -6.0*h23, (3.0*h1\_2 + h0*(3.0*h1 + h0)), &
1048                     h1\_2,    h2\_2, (3.0*h2\_2 + h3*(3.0*h2 + h3)) /)
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))), &
1050                    -h1\_3,       h2\_3,  (4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)
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))))
      , &
1052                     h1\_4,       h2\_4,       (5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3))))
       /)
1053   asys(1:6,6) = (/  0.0, -30.0*(h23\_2*h23\_2), &
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))))), &
1055                    -h1\_5, h2\_5, &
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))))) /)
1057   bsys(1:6) = (/ 0.0, -2.0, 6.0*h2, 12.0*h2\_2, 20.0*h2\_3, 30.0*h2\_4 /)
1058 
1059   \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )
1060 
1061   alpha = csys(1)
1062   beta  = csys(2)
1063 
1064   tri\_l(n) = alpha
1065   tri\_d(n) = 1.0
1066   tri\_u(n) = beta
1067   tri\_b(n) = csys(3) * u(n-3) + csys(4) * u(n-2) + csys(5) * u(n-1) + csys(6) * u(n)
1068 
1069   \textcolor{comment}{! Boundary conditions: right boundary}
1070   x(1) = 0.0
1071   \textcolor{keywordflow}{do} i = 1,6
1072     dx = h(n+1-i)
1073     xavg = x(i) + 0.5*dx
1074     asys(1:6,i) =  (/ 1.0, xavg, (xavg**2 + c1\_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), &
1075                       (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), &
1076                        xavg * (xavg**4 + c5\_6*xavg**2*dx**2 + 0.0625*dx**4) /)
1077     bsys(i) = u(n+1-i)
1078     x(i+1) = x(i) + dx
1079 \textcolor{keywordflow}{  enddo}
1080 
1081   \textcolor{keyword}{call }linear\_solver( 6,  asys, bsys, csys )
1082 
1083   tri\_l(n+1) = 0.0
1084   tri\_d(n+1) = 1.0
1085   tri\_u(n+1) = 0.0
1086   tri\_b(n+1) = -csys(2)
1087 
1088   \textcolor{comment}{! Solve tridiagonal system and assign edge values}
1089   \textcolor{keyword}{call }solve\_tridiagonal\_system( tri\_l, tri\_d, tri\_u, tri\_b, tri\_x, n+1 )
1090 
1091   \textcolor{keywordflow}{do} i = 2,n
1092     edge\_slopes(i,1)   = tri\_x(i)
1093     edge\_slopes(i-1,2) = tri\_x(i)
1094 \textcolor{keywordflow}{  enddo}
1095   edge\_slopes(1,1) = tri\_x(1)
1096   edge\_slopes(n,2) = tri\_x(n+1)
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}}
\subsubsection{\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{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt 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}
175   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of cells}
176   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths [H]}
177   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}
178   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Returned edge values [A]; the}
179 \textcolor{comment}{                                           !! second index is for the two edges of each cell.}
180 
181   \textcolor{comment}{! Local variables}
182   \textcolor{keywordtype}{integer}   :: k        \textcolor{comment}{! loop index}
183 
184   \textcolor{comment}{! Boundary edge values are simply equal to the boundary cell averages}
185   edge\_val(1,1) = u(1)
186   edge\_val(n,2) = u(n)
187 
188   \textcolor{keywordflow}{do} k = 2,n
189     \textcolor{comment}{! Compute left edge value}
190     \textcolor{keywordflow}{if} (h(k-1) + h(k) == 0.0) \textcolor{keywordflow}{then}    \textcolor{comment}{! Avoid singularities when h0+h1=0}
191       edge\_val(k,1) = 0.5 * (u(k-1) + u(k))
192     \textcolor{keywordflow}{else}
193       edge\_val(k,1) = ( u(k-1)*h(k) + u(k)*h(k-1) ) / ( h(k-1) + h(k) )
194 \textcolor{keywordflow}{    endif}
195 
196     \textcolor{comment}{! Left edge value of the current cell is equal to right edge value of left cell}
197     edge\_val(k-1,2) = edge\_val(k,1)
198 \textcolor{keywordflow}{  enddo}
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}}
\subsubsection{\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{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+val} & Returned edge values \mbox{[}A\mbox{]}; the second index is for the two edges of each cell.\\
\hline
\mbox{\tt in}  & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt 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}
222   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of cells}
223   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths [H]}
224   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}
225   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Returned edge values [A]; the second index}
226 \textcolor{comment}{                                           !! is for the two edges of each cell.}
227   \textcolor{keywordtype}{real},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}
228   \textcolor{keywordtype}{logical},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}
229 
230   \textcolor{comment}{! Local variables}
231   \textcolor{keywordtype}{real} :: h0, h1, h2, h3        \textcolor{comment}{! temporary thicknesses [H]}
232   \textcolor{keywordtype}{real} :: h\_sum                 \textcolor{comment}{! A sum of adjacent thicknesses [H]}
233   \textcolor{keywordtype}{real} :: h\_min                 \textcolor{comment}{! A minimal cell width [H]}
234   \textcolor{keywordtype}{real} :: f1, f2, f3            \textcolor{comment}{! auxiliary variables with various units}
235   \textcolor{keywordtype}{real} :: et1, et2, et3         \textcolor{comment}{! terms the expresson for edge values [A H]}
236   \textcolor{keywordtype}{real} :: i\_h12                 \textcolor{comment}{! The inverse of the sum of the two central thicknesses [H-1]}
237   \textcolor{keywordtype}{real} :: i\_h012, i\_h123        \textcolor{comment}{! Inverses of sums of three succesive thicknesses [H-1]}
238   \textcolor{keywordtype}{real} :: i\_den\_et2, i\_den\_et3  \textcolor{comment}{! Inverses of denominators in edge value terms [H-2]}
239   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(5)}    :: x          \textcolor{comment}{! Coordinate system with 0 at edges [H]}
240   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}    :: dz               \textcolor{comment}{! A temporary array of limited layer thicknesses [H]}
241   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}    :: u\_tmp            \textcolor{comment}{! A temporary array of cell average properties [A]}
242   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}       :: c1\_12 = 1.0 / 12.0
243   \textcolor{keywordtype}{real}                  :: dx, xavg         \textcolor{comment}{! Differences and averages of successive values of x [H]}
244   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4,4)}  :: a                \textcolor{comment}{! values near the boundaries}
245   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}    :: b, c
246   \textcolor{keywordtype}{real}      :: hneglect \textcolor{comment}{! A negligible thickness in the same units as h.}
247   \textcolor{keywordtype}{integer}               :: i, j
248   \textcolor{keywordtype}{logical}   :: use\_2018\_answers  \textcolor{comment}{! If true use older, less acccurate expressions.}
249 
250   use\_2018\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) use\_2018\_answers = answers\_2018
251   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
252     hneglect = hneglect\_edge\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect
253   \textcolor{keywordflow}{else}
254     hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect
255 \textcolor{keywordflow}{  endif}
256 
257   \textcolor{comment}{! Loop on interior cells}
258   \textcolor{keywordflow}{do} i = 3,n-1
259 
260     h0 = h(i-2)
261     h1 = h(i-1)
262     h2 = h(i)
263     h3 = h(i+1)
264 
265     \textcolor{comment}{! Avoid singularities when consecutive pairs of h vanish}
266     \textcolor{keywordflow}{if} (h0+h1==0.0 .or. h1+h2==0.0 .or. h2+h3==0.0) \textcolor{keywordflow}{then}
267       \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
268         h\_min = hminfrac*max( hneglect, h0+h1+h2+h3 )
269       \textcolor{keywordflow}{else}
270         h\_min = hminfrac*max( hneglect, (h0+h1)+(h2+h3) )
271 \textcolor{keywordflow}{      endif}
272       h0 = max( h\_min, h(i-2) )
273       h1 = max( h\_min, h(i-1) )
274       h2 = max( h\_min, h(i) )
275       h3 = max( h\_min, h(i+1) )
276 \textcolor{keywordflow}{    endif}
277 
278     \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
279       f1 = (h0+h1) * (h2+h3) / (h1+h2)
280       f2 = h2 * u(i-1) + h1 * u(i)
281       f3 = 1.0 / (h0+h1+h2) + 1.0 / (h1+h2+h3)
282       et1 = f1 * f2 * f3
283       et2 = ( h2 * (h2+h3) / ( (h0+h1+h2)*(h0+h1) ) ) * &
284             ((h0+2.0*h1) * u(i-1) - h1 * u(i-2))
285       et3 = ( h1 * (h0+h1) / ( (h1+h2+h3)*(h2+h3) ) ) * &
286             ((2.0*h2+h3) * u(i) - h2 * u(i+1))
287       edge\_val(i,1) = (et1 + et2 + et3) / ( h0 + h1 + h2 + h3)
288     \textcolor{keywordflow}{else}
289       i\_h12 = 1.0 / (h1+h2)
290       i\_den\_et2 = 1.0 / ( ((h0+h1)+h2)*(h0+h1) ) ; i\_h012 = (h0+h1) * i\_den\_et2
291       i\_den\_et3 = 1.0 / ( (h1+(h2+h3))*(h2+h3) ) ; i\_h123 = (h2+h3) * i\_den\_et3
292 
293       et1 = ( 1.0 + (h1 * i\_h012 + (h0+h1) * i\_h123) ) * i\_h12 * (h2*(h2+h3)) * u(i-1) + &
294             ( 1.0 + (h2 * i\_h123 + (h2+h3) * i\_h012) ) * i\_h12 * (h1*(h0+h1)) * u(i)
295       et2 = ( h1 * (h2*(h2+h3)) * i\_den\_et2 ) * (u(i-1)-u(i-2))
296       et3 = ( h2 * (h1*(h0+h1)) * i\_den\_et3 ) * (u(i) - u(i+1))
297       edge\_val(i,1) = (et1 + (et2 + et3)) / ((h0 + h1) + (h2 + h3))
298 \textcolor{keywordflow}{    endif}
299     edge\_val(i-1,2) = edge\_val(i,1)
300 
301 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on interior cells}
302 
303   \textcolor{comment}{! Determine first two edge values}
304   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
305     h\_min = max( hneglect, hminfrac*sum(h(1:4)) )
306     x(1) = 0.0
307     \textcolor{keywordflow}{do} i = 1,4
308       dx = max(h\_min, h(i) )
309       x(i+1) = x(i) + dx
310       \textcolor{keywordflow}{do} j = 1,4 ; a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / \textcolor{keywordtype}{real(j)} ; enddo
311       b(i) = u(i) * dx
312 \textcolor{keywordflow}{    enddo}
313 
314     \textcolor{keyword}{call }solve\_linear\_system( a, b, c, 4 )
315 
316     \textcolor{comment}{! Set the edge values of the first cell}
317     edge\_val(1,1) = evaluation\_polynomial( c, 4, x(1) )
318     edge\_val(1,2) = evaluation\_polynomial( c, 4, x(2) )
319   \textcolor{keywordflow}{else}  \textcolor{comment}{! Use expressions with less sensitivity to roundoff}
320     \textcolor{keywordflow}{do} i=1,4 ; dz(i) = max(hneglect, h(i) ) ; u\_tmp(i) = u(i) ;\textcolor{keywordflow}{ enddo}
321     \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, c)
322 
323     \textcolor{comment}{! Set the edge values of the first cell}
324     edge\_val(1,1) = c(1)
325     edge\_val(1,2) = c(1) + dz(1)*(c(2) + dz(1)*(c(3) + dz(1)*c(4)))
326 \textcolor{keywordflow}{  endif}
327   edge\_val(2,1) = edge\_val(1,2)
328 
329   \textcolor{comment}{! Determine two edge values of the last cell}
330   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
331     h\_min = max( hneglect, hminfrac*sum(h(n-3:n)) )
332 
333     x(1) = 0.0
334     \textcolor{keywordflow}{do} i = 1,4
335       dx = max(h\_min, h(n-4+i) )
336       x(i+1) = x(i) + dx
337       \textcolor{keywordflow}{do} j = 1,4 ; a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / \textcolor{keywordtype}{real(j)} ; enddo
338       b(i) = u(n-4+i) * dx
339 \textcolor{keywordflow}{    enddo}
340 
341     \textcolor{keyword}{call }solve\_linear\_system( a, b, c, 4 )
342 
343     \textcolor{comment}{! Set the last and second to last edge values}
344     edge\_val(n,2) = evaluation\_polynomial( c, 4, x(5) )
345     edge\_val(n,1) = evaluation\_polynomial( c, 4, x(4) )
346   \textcolor{keywordflow}{else}
347     \textcolor{comment}{! Use expressions with less sensitivity to roundoff, including using a coordinate}
348     \textcolor{comment}{! system that sets the origin at the last interface in the domain.}
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}
350     \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, c)
351 
352     \textcolor{comment}{! Set the last and second to last edge values}
353     edge\_val(n,2) = c(1)
354     edge\_val(n,1) = c(1) + dz(1)*(c(2) + dz(1)*(c(3) + dz(1)*c(4)))
355 \textcolor{keywordflow}{  endif}
356   edge\_val(n-1,2) = edge\_val(n,1)
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}}
\subsubsection{\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{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell average properties in arbitrary units \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+val} & Returned edge values \mbox{[}A\mbox{]}; the second index is for the two edges of each cell.\\
\hline
\mbox{\tt in}  & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt 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}
387   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of cells}
388   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths [H]}
389   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell average properties in arbitrary units [A]}
390   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{ !< Returned edge values [A]; the second index}
391 \textcolor{comment}{                                           !! is for the two edges of each cell.}
392   \textcolor{keywordtype}{real},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}
393   \textcolor{keywordtype}{logical},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}
394 
395   \textcolor{comment}{! Local variables}
396   \textcolor{keywordtype}{integer}               :: i, j                 \textcolor{comment}{! loop indexes}
397   \textcolor{keywordtype}{real}                  :: h0, h1, h2           \textcolor{comment}{! cell widths [H]}
398   \textcolor{keywordtype}{real}                  :: h\_min                \textcolor{comment}{! A minimal cell width [H]}
399   \textcolor{keywordtype}{real}                  :: h\_sum                \textcolor{comment}{! A sum of adjacent thicknesses [H]}
400   \textcolor{keywordtype}{real}                  :: h0\_2, h1\_2, h0h1
401   \textcolor{keywordtype}{real}                  :: h0ph1\_2, h0ph1\_4
402   \textcolor{keywordtype}{real}                  :: alpha, beta          \textcolor{comment}{! stencil coefficients [nondim]}
403   \textcolor{keywordtype}{real}                  :: i\_h2, abmix          \textcolor{comment}{! stencil coefficients [nondim]}
404   \textcolor{keywordtype}{real}                  :: a, b
405   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(5)}    :: x                    \textcolor{comment}{! Coordinate system with 0 at edges [H]}
406   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}       :: c1\_12 = 1.0 / 12.0
407   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}       :: c1\_3 = 1.0 / 3.0
408   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}    :: dz                   \textcolor{comment}{! A temporary array of limited layer thicknesses [H]}
409   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}    :: u\_tmp                \textcolor{comment}{! A temporary array of cell average properties [A]}
410   \textcolor{keywordtype}{real}                  :: dx, xavg             \textcolor{comment}{! Differences and averages of successive values of x [H]}
411   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4,4)}  :: asys                 \textcolor{comment}{! boundary conditions}
412   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}    :: bsys, csys
413   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N+1)}  :: tri\_l, &     \textcolor{comment}{! tridiagonal system (lower diagonal) [nondim]}
414                            tri\_d, &     \textcolor{comment}{! tridiagonal system (middle diagonal) [nondim]}
415                            tri\_c, &     \textcolor{comment}{! tridiagonal system central value, with tri\_d = tri\_c+tri\_l+tri\_u}
416                            tri\_u, &     \textcolor{comment}{! tridiagonal system (upper diagonal) [nondim]}
417                            tri\_b, &     \textcolor{comment}{! tridiagonal system (right hand side) [A]}
418                            tri\_x        \textcolor{comment}{! tridiagonal system (solution vector) [A]}
419   \textcolor{keywordtype}{real}      :: hneglect          \textcolor{comment}{! A negligible thickness [H]}
420   \textcolor{keywordtype}{logical}   :: use\_2018\_answers  \textcolor{comment}{! If true use older, less acccurate expressions.}
421 
422   use\_2018\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) use\_2018\_answers = answers\_2018
423   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
424     hneglect = hneglect\_edge\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect
425   \textcolor{keywordflow}{else}
426     hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect
427 \textcolor{keywordflow}{  endif}
428 
429   \textcolor{comment}{! Loop on cells (except last one)}
430   \textcolor{keywordflow}{do} i = 1,n-1
431     \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
432       \textcolor{comment}{! Get cell widths}
433       h0 = h(i)
434       h1 = h(i+1)
435       \textcolor{comment}{! Avoid singularities when h0+h1=0}
436       \textcolor{keywordflow}{if} (h0+h1==0.) \textcolor{keywordflow}{then}
437         h0 = hneglect
438         h1 = hneglect
439 \textcolor{keywordflow}{      endif}
440 
441       \textcolor{comment}{! Auxiliary calculations}
442       h0ph1\_2 = (h0 + h1)**2
443       h0ph1\_4 = h0ph1\_2**2
444       h0h1 = h0 * h1
445       h0\_2 = h0 * h0
446       h1\_2 = h1 * h1
447 
448       \textcolor{comment}{! Coefficients}
449       alpha = h1\_2 / h0ph1\_2
450       beta = h0\_2 / h0ph1\_2
451       a = 2.0 * h1\_2 * ( h1\_2 + 2.0 * h0\_2 + 3.0 * h0h1 ) / h0ph1\_4
452       b = 2.0 * h0\_2 * ( h0\_2 + 2.0 * h1\_2 + 3.0 * h0h1 ) / h0ph1\_4
453 
454       tri\_d(i+1) = 1.0
455     \textcolor{keywordflow}{else}  \textcolor{comment}{! Use expressions with less sensitivity to roundoff}
456       \textcolor{comment}{! Get cell widths}
457       h0 = max(h(i), hneglect)
458       h1 = max(h(i+1), hneglect)
459       \textcolor{comment}{! The 1e-12 here attempts to balance truncation errors from the differences of}
460       \textcolor{comment}{! large numbers against errors from approximating thin layers as non-vanishing.}
461       \textcolor{keywordflow}{if} (abs(h0) < 1.0e-12*abs(h1)) h0 = 1.0e-12*h1
462       \textcolor{keywordflow}{if} (abs(h1) < 1.0e-12*abs(h0)) h1 = 1.0e-12*h0
463       i\_h2 = 1.0 / ((h0 + h1)**2)
464       alpha = (h1 * h1) * i\_h2
465       beta = (h0 * h0) * i\_h2
466       abmix = (h0 * h1) * i\_h2
467       a = 2.0 * alpha * ( alpha + 2.0 * beta + 3.0 * abmix )
468       b = 2.0 * beta * ( beta + 2.0 * alpha + 3.0 * abmix )
469 
470       tri\_c(i+1) = 2.0*abmix  \textcolor{comment}{! = 1.0 - alpha - beta}
471 \textcolor{keywordflow}{    endif}
472 
473     tri\_l(i+1) = alpha
474     tri\_u(i+1) = beta
475 
476     tri\_b(i+1) = a * u(i) + b * u(i+1)
477 
478 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on cells}
479 
480   \textcolor{comment}{! Boundary conditions: set the first boundary value}
481   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
482     h\_min = max( hneglect, hminfrac*sum(h(1:4)) )
483     x(1) = 0.0
484     \textcolor{keywordflow}{do} i = 1,4
485       dx = max(h\_min, h(i) )
486       x(i+1) = x(i) + dx
487       \textcolor{keywordflow}{do} j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ;\textcolor{keywordflow}{ enddo}
488       bsys(i) = u(i) * dx
489 \textcolor{keywordflow}{    enddo}
490 
491     \textcolor{keyword}{call }solve\_linear\_system( asys, bsys, csys, 4 )
492 
493     tri\_b(1) = evaluation\_polynomial( csys, 4, x(1) )  \textcolor{comment}{! Set the first edge value}
494     tri\_d(1) = 1.0
495   \textcolor{keywordflow}{else} \textcolor{comment}{! Use expressions with less sensitivity to roundoff}
496     \textcolor{keywordflow}{do} i=1,4 ; dz(i) = max(hneglect, h(i) ) ; u\_tmp(i) = u(i) ;\textcolor{keywordflow}{ enddo}
497     \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, csys)
498 
499     tri\_b(1) = csys(1)  \textcolor{comment}{! Set the first edge value.}
500     tri\_c(1) = 1.0
501 \textcolor{keywordflow}{  endif}
502   tri\_u(1) = 0.0 \textcolor{comment}{! tri\_l(1) = 0.0}
503 
504   \textcolor{comment}{! Boundary conditions: set the last boundary value}
505   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
506     h\_min = max( hneglect, hminfrac*sum(h(n-3:n)) )
507     x(1) = 0.0
508     \textcolor{keywordflow}{do} i=1,4
509       dx = max(h\_min, h(n-4+i) )
510       x(i+1) = x(i) + dx
511       \textcolor{keywordflow}{do} j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ;\textcolor{keywordflow}{ enddo}
512       bsys(i) = u(n-4+i) * dx
513 \textcolor{keywordflow}{    enddo}
514 
515     \textcolor{keyword}{call }solve\_linear\_system( asys, bsys, csys, 4 )
516 
517     \textcolor{comment}{! Set the last edge value}
518     tri\_b(n+1) = evaluation\_polynomial( csys, 4, x(5) )
519     tri\_d(n+1) = 1.0
520 
521   \textcolor{keywordflow}{else}
522     \textcolor{comment}{! Use expressions with less sensitivity to roundoff, including using a coordinate}
523     \textcolor{comment}{! system that sets the origin at the last interface in the domain.}
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}
525     \textcolor{keyword}{call }end\_value\_h4(dz, u\_tmp, csys)
526 
527     tri\_b(n+1) = csys(1) \textcolor{comment}{! Set the last edge value}
528     tri\_c(n+1) = 1.0
529 \textcolor{keywordflow}{  endif}
530   tri\_l(n+1) = 0.0 \textcolor{comment}{! tri\_u(N+1) = 0.0}
531 
532   \textcolor{comment}{! Solve tridiagonal system and assign edge values}
533   \textcolor{keywordflow}{if} (use\_2018\_answers) \textcolor{keywordflow}{then}
534     \textcolor{keyword}{call }solve\_tridiagonal\_system( tri\_l, tri\_d, tri\_u, tri\_b, tri\_x, n+1 )
535   \textcolor{keywordflow}{else}
536     \textcolor{keyword}{call }solve\_diag\_dominant\_tridiag( tri\_l, tri\_c, tri\_u, tri\_b, tri\_x, n+1 )
537 \textcolor{keywordflow}{  endif}
538 
539   edge\_val(1,1) = tri\_x(1)
540   \textcolor{keywordflow}{do} i=2,n
541     edge\_val(i,1)   = tri\_x(i)
542     edge\_val(i-1,2) = tri\_x(i)
543 \textcolor{keywordflow}{  enddo}
544   edge\_val(n,2) = tri\_x(n+1)
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}}
\subsubsection{\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}{\tt 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{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell average properties (size N) in arbitrary units \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+val} & Returned edge values \mbox{[}A\mbox{]}; the second index is for the two edges of each cell.\\
\hline
\mbox{\tt in}  & {\em h\+\_\+neglect} & A negligibly small width \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt 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}
1137   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of cells}
1138   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths [H]}
1139   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell average properties (size N) in arbitrary units [A]}
1140   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,2)}, \textcolor{keywordtype}{intent(inout)} :: edge\_val\textcolor{comment}{  !< Returned edge values [A]; the second index}
1141 \textcolor{comment}{                                           !! is for the two edges of each cell.}
1142   \textcolor{keywordtype}{real},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: h\_neglect\textcolor{comment}{ !< A negligibly small width [H]}
1143   \textcolor{keywordtype}{logical},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}
1144 
1145   \textcolor{comment}{! Local variables}
1146   \textcolor{keywordtype}{real} :: h0, h1, h2, h3       \textcolor{comment}{! cell widths [H]}
1147   \textcolor{keywordtype}{real} :: hmin                 \textcolor{comment}{! The minimum thickness used in these calculations [H]}
1148   \textcolor{keywordtype}{real} :: h01, h01\_2, h01\_3  \textcolor{comment}{! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n]}
1149   \textcolor{keywordtype}{real} :: h23, h23\_2, h23\_3  \textcolor{comment}{! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n]}
1150   \textcolor{keywordtype}{real} :: hneglect             \textcolor{comment}{! A negligible thickness [H].}
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]}
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]}
1153   \textcolor{keywordtype}{real}                  :: alpha, beta          \textcolor{comment}{! stencil coefficients}
1154   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(7)}    :: x          \textcolor{comment}{! Coordinate system with 0 at edges [same units as h]}
1155   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}       :: c1\_12 = 1.0 / 12.0
1156   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}       :: c5\_6 = 5.0 / 6.0
1157   \textcolor{keywordtype}{real}                  :: dx, xavg   \textcolor{comment}{! Differences and averages of successive values of x [same units as
       h]}
1158   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(6,6)}  :: asys                 \textcolor{comment}{! boundary conditions}
1159   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(6)}    :: bsys, csys           \textcolor{comment}{! ...}
1160   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N+1)}  :: tri\_l, &             \textcolor{comment}{! trid. system (lower diagonal)}
1161                            tri\_d, &             \textcolor{comment}{! trid. system (middle diagonal)}
1162                            tri\_u, &             \textcolor{comment}{! trid. system (upper diagonal)}
1163                            tri\_b, &             \textcolor{comment}{! trid. system (unknowns vector)}
1164                            tri\_x                \textcolor{comment}{! trid. system (rhs)}
1165   \textcolor{keywordtype}{integer} :: i, j, k              \textcolor{comment}{! loop indexes}
1166 
1167   hneglect = hneglect\_edge\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect
1168 
1169   \textcolor{comment}{! Loop on interior cells}
1170   \textcolor{keywordflow}{do} k = 2,n-2
1171     \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}
1172     hmin = max(hneglect, hminfrac*((h(k-1) + h(k)) + (h(k+1) + h(k+2))))
1173     h0 = max(h(k-1), hmin) ; h1 = max(h(k), hmin)
1174     h2 = max(h(k+1), hmin) ; h3 = max(h(k+2), hmin)
1175 
1176     \textcolor{comment}{! Auxiliary calculations}
1177     h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2
1178     h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2
1179 
1180     \textcolor{comment}{! Compute matrix entries as described in Eq. (48) of White and Adcroft (2009)}
1181     asys(1:6,1) = (/ 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 /)
1182     asys(1:6,2) = (/ -2.0*h1, 2.0*h2, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /)
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}
1184                      -h1\_2,    -h2\_2,    -(3.0*h2\_2 + h3*(3.0*h2 + h3)) /) \textcolor{comment}{! = -((h2+h3)**3 - h2**3) / h3}
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))), &
1186                       h1\_3,    -h2\_3,    -(4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)
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)))),
       &
1188                      -h1\_4,    -h2\_4,     -(5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))) 
      /)
1189     asys(1:6,6) = (/ -6.0*h1\_5, 6.0*h2\_5, &
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))))), &
1191                       h1\_5, -h2\_5, &
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))))) /)
1193     bsys(1:6) = (/ -1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
1194 
1195     \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )
1196 
1197     alpha = csys(1)
1198     beta  = csys(2)
1199 
1200     tri\_l(k+1) = alpha
1201     tri\_d(k+1) = 1.0
1202     tri\_u(k+1) = beta
1203     tri\_b(k+1) = csys(3) * u(k-1) + csys(4) * u(k) + csys(5) * u(k+1) + csys(6) * u(k+2)
1204 
1205 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on cells}
1206 
1207   \textcolor{comment}{! Use a right-biased stencil for the second row, as described in Eq. (49) of White and Adcroft (2009).}
1208 
1209   \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}
1210   hmin = max(hneglect, hminfrac*((h(1) + h(2)) + (h(3) + h(4))))
1211   h0 = max(h(1), hmin) ; h1 = max(h(2), hmin)
1212   h2 = max(h(3), hmin) ; h3 = max(h(4), hmin)
1213 
1214   \textcolor{comment}{! Auxiliary calculations}
1215   h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2
1216   h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2
1217   h01 = h0 + h1 ; h01\_2 = h01 * h01 ; h01\_3 = h01 * h01\_2
1218 
1219   \textcolor{comment}{! Compute matrix entries}
1220   asys(1:6,1) = (/ 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 /)
1221   asys(1:6,2) = (/  -2.0*h01, 0.0, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /)
1222   asys(1:6,3) = (/ 3.0*h01\_2, 0.0, -(3.0*h1\_2 + h0*(3.0*h1 + h0)), &
1223                    -h1\_2, -h2\_2, -(3.0*h2\_2 + h3*(3.0*h2 + h3)) /)
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))), &
1225                    h1\_3, -h2\_3, -(4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)
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))))
      , &
1227                    -h1\_4, -h2\_4, -(5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3)))) /)
1228   asys(1:6,6) = (/ -6.0*(h01\_3*h01\_2), 0.0, &
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))))), &
1230                    h1\_5, - h2\_5, &
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))))) /)
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 /)
1233 
1234   \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )
1235 
1236   alpha = csys(1)
1237   beta  = csys(2)
1238 
1239   tri\_l(2) = alpha
1240   tri\_d(2) = 1.0
1241   tri\_u(2) = beta
1242   tri\_b(2) = csys(3) * u(1) + csys(4) * u(2) + csys(5) * u(3) + csys(6) * u(4)
1243 
1244   \textcolor{comment}{! Boundary conditions: left boundary}
1245   hmin = max( hneglect, hminfrac*((h(1)+h(2)) + (h(5)+h(6)) + (h(3)+h(4))) )
1246   x(1) = 0.0
1247   \textcolor{keywordflow}{do} i = 1,6
1248     dx = max( hmin, h(i) )
1249     xavg = x(i) + 0.5*dx
1250     asys(1:6,i) =  (/ 1.0, xavg, (xavg**2 + c1\_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), &
1251                       (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), &
1252                        xavg * (xavg**4 + c5\_6*xavg**2*dx**2 + 0.0625*dx**4) /)
1253     bsys(i) = u(i)
1254     x(i+1) = x(i) + dx
1255 \textcolor{keywordflow}{  enddo}
1256 
1257   \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )
1258 
1259   tri\_l(1) = 0.0
1260   tri\_d(1) = 1.0
1261   tri\_u(1) = 0.0
1262   tri\_b(1) = evaluation\_polynomial( csys, 6, x(1) )        \textcolor{comment}{! first edge value}
1263 
1264   \textcolor{comment}{! Use a left-biased stencil for the second to last row, as described in Eq. (50) of White and Adcroft
       (2009).}
1265 
1266   \textcolor{comment}{! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.}
1267   hmin = max(hneglect, hminfrac*((h(n-3) + h(n-2)) + (h(n-1) + h(n))))
1268   h0 = max(h(n-3), hmin) ; h1 = max(h(n-2), hmin)
1269   h2 = max(h(n-1), hmin) ; h3 = max(h(n), hmin)
1270 
1271   \textcolor{comment}{! Auxiliary calculations}
1272   h1\_2 = h1 * h1 ; h1\_3 = h1\_2 * h1 ; h1\_4 = h1\_2 * h1\_2 ; h1\_5 = h1\_3 * h1\_2
1273   h2\_2 = h2 * h2 ; h2\_3 = h2\_2 * h2 ; h2\_4 = h2\_2 * h2\_2 ; h2\_5 = h2\_3 * h2\_2
1274   h23 = h2 + h3 ; h23\_2 = h23 * h23 ; h23\_3 = h23 * h23\_2
1275 
1276   \textcolor{comment}{! Compute matrix entries}
1277   asys(1:6,1) = (/ 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 /)
1278   asys(1:6,2) = (/ 0.0, 2.0*h23, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /)
1279   asys(1:6,3) = (/ 0.0, 3.0*h23\_2, -(3.0*h1\_2 + h0*(3.0*h1 + h0)), &
1280                  -h1\_2,     -h2\_2, -(3.0*h2\_2 + h3*(3.0*h2 + h3)) /)
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))), &
1282                   h1\_3,     -h2\_3, -(4.0*h2\_3 + h3*(6.0*h2\_2 + h3*(4.0*h2 + h3))) /)
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))))
      , &
1284                  -h1\_4,     -h2\_4,         -(5.0*h2\_4 + h3*(10.0*h2\_3 + h3*(10.0*h2\_2 + h3*(5.0*h2 + h3))))
       /)
1285   asys(1:6,6) = (/ 0.0, 6.0*(h23\_3*h23\_2), &
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))))), &
1287                    h1\_5, -h2\_5, &
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))))) /)
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 /)
1290 
1291   \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )
1292 
1293   alpha = csys(1)
1294   beta  = csys(2)
1295 
1296   tri\_l(n) = alpha
1297   tri\_d(n) = 1.0
1298   tri\_u(n) = beta
1299   tri\_b(n) = csys(3) * u(n-3) + csys(4) * u(n-2) + csys(5) * u(n-1) + csys(6) * u(n)
1300 
1301   \textcolor{comment}{! Boundary conditions: right boundary}
1302   hmin = max( hneglect, hminfrac*(h(n-3) + h(n-2)) + ((h(n-1) + h(n)) + (h(n-5) + h(n-4))) )
1303   x(1) = 0.0
1304   \textcolor{keywordflow}{do} i = 1,6
1305     dx = max( hmin, h(n+1-i) )
1306     xavg = x(i) + 0.5 * dx
1307     asys(1:6,i) =  (/ 1.0, xavg, (xavg**2 + c1\_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), &
1308                       (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), &
1309                        xavg * (xavg**4 + c5\_6*xavg**2*dx**2 + 0.0625*dx**4) /)
1310     bsys(i) = u(n+1-i)
1311     x(i+1) = x(i) + dx
1312 \textcolor{keywordflow}{  enddo}
1313 
1314   \textcolor{keyword}{call }linear\_solver( 6, asys, bsys, csys )
1315 
1316   tri\_l(n+1) = 0.0
1317   tri\_d(n+1) = 1.0
1318   tri\_u(n+1) = 0.0
1319   tri\_b(n+1) = csys(1)
1320 
1321   \textcolor{comment}{! Solve tridiagonal system and assign edge values}
1322   \textcolor{keyword}{call }solve\_tridiagonal\_system( tri\_l, tri\_d, tri\_u, tri\_b, tri\_x, n+1 )
1323 
1324   \textcolor{keywordflow}{do} i = 2,n
1325     edge\_val(i,1)   = tri\_x(i)
1326     edge\_val(i-1,2) = tri\_x(i)
1327 \textcolor{keywordflow}{  enddo}
1328   edge\_val(1,1) = tri\_x(1)
1329   edge\_val(n,2) = tri\_x(n+1)
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}}
\subsubsection{\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{\tt in}  & {\em dz} & The thicknesses of 4 layers, starting at the edge \mbox{[}H\mbox{]}. The values of dz must be positive.\\
\hline
\mbox{\tt in}  & {\em u} & The average properties of 4 layers, starting at the edge \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt 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}
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].}
552 \textcolor{comment}{                                           !! The values of dz must be positive.}
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]}
554   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}, \textcolor{keywordtype}{intent(out)} :: csys\textcolor{comment}{  !< The four coefficients of a 4th order polynomial fit}
555 \textcolor{comment}{                                           !! of u as a function of z [A H-(n-1)]}
556 
557   \textcolor{comment}{! Local variables}
558   \textcolor{keywordtype}{real} :: wt(3,4)         \textcolor{comment}{! The weights of successive u differences in the 4 closed form expressions.}
559                           \textcolor{comment}{! The units of Wt vary with the second index as [H-(n-1)].}
560   \textcolor{keywordtype}{real} :: h1, h2, h3, h4  \textcolor{comment}{! Copies of the layer thicknesses [H]}
561   \textcolor{keywordtype}{real} :: h12, h23, h34   \textcolor{comment}{! Sums of two successive thicknesses [H]}
562   \textcolor{keywordtype}{real} :: h123, h234      \textcolor{comment}{! Sums of three successive thicknesses [H]}
563   \textcolor{keywordtype}{real} :: h1234           \textcolor{comment}{! Sums of all four thicknesses [H]}
564   \textcolor{comment}{! real :: I\_h1          ! The inverse of the a thickness [H-1]}
565   \textcolor{keywordtype}{real} :: i\_h12, i\_h23, i\_h34 \textcolor{comment}{! The inverses of sums of two thicknesses [H-1]}
566   \textcolor{keywordtype}{real} :: i\_h123, i\_h234  \textcolor{comment}{! The inverse of the sum of three thicknesses [H-1]}
567   \textcolor{keywordtype}{real} :: i\_h1234         \textcolor{comment}{! The inverse of the sum of all four thicknesses [H-1]}
568   \textcolor{keywordtype}{real} :: i\_denom         \textcolor{comment}{! The inverse of the denominator some expressions [H-3]}
569   \textcolor{keywordtype}{real} :: i\_denb3         \textcolor{comment}{! The inverse of the product of three sums of thicknesses [H-3]}
570   \textcolor{keywordtype}{real} :: min\_frac = 1.0e-6  \textcolor{comment}{! The square of min\_frac should be much larger than roundoff [nondim]}
571   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: c1\_3 = 1.0 / 3.0
572   \textcolor{keywordtype}{integer} :: i, j, k
573 
574   \textcolor{comment}{! These are only used for code verification}
575   \textcolor{comment}{! real, dimension(4) :: Atest  ! The  coefficients of an expression that is being tested.}
576   \textcolor{comment}{! real :: zavg, u\_mag, c\_mag}
577   \textcolor{comment}{! character(len=128) :: mesg}
578   \textcolor{comment}{! real, parameter :: C1\_12 = 1.0 / 12.0}
579 
580  \textcolor{comment}{! if ((dz(1) == dz(2)) .and. (dz(1) == dz(3)) .and. (dz(1) == dz(4))) then}
581  \textcolor{comment}{!   ! There are simple closed-form expressions in this case}
582  \textcolor{comment}{!   I\_h1 = 0.0 ; if (dz(1) > 0.0) I\_h1 = 1.0 / dz(1)}
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)}
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)}
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)}
586  \textcolor{comment}{!   Csys(4) = ((u(2)-u(1)) - 2.0 * (u(3)-u(2)) + (u(4)-u(3))) * (0.5*C1\_3)}
587  \textcolor{comment}{! else}
588 
589   \textcolor{comment}{! Express the coefficients as sums of the differences between properties of succesive layers.}
590 
591   h1 = dz(1) ; h2 = dz(2) ; h3 = dz(3) ; h4 = dz(4)
592   \textcolor{comment}{! Some of the weights used below are proportional to (h1/(h2+h3))**2 or (h1/(h2+h3))*(h2/(h3+h4))}
593   \textcolor{comment}{! so h2 and h3 should be adjusted to ensure that these ratios are not so large that property}
594   \textcolor{comment}{! differences at the level of roundoff are amplified to be of order 1.}
595   \textcolor{keywordflow}{if} ((h2+h3) < min\_frac*h1) h3 = min\_frac*h1 - h2
596   \textcolor{keywordflow}{if} ((h3+h4) < min\_frac*h1) h4 = min\_frac*h1 - h3
597 
598   h12 = h1+h2 ; h23 = h2+h3 ; h34 = h3+h4
599   h123 = h12 + h3 ; h234 = h2 + h34 ; h1234 = h12 + h34
600   \textcolor{comment}{! Find 3 reciprocals with a single division for efficiency.}
601   i\_denb3 = 1.0 / (h123 * h12 * h23)
602   i\_h12 = (h123 * h23) * i\_denb3
603   i\_h23 = (h12 * h123) * i\_denb3
604   i\_h123 = (h12 * h23) * i\_denb3
605   i\_denom = 1.0 / ( h1234 * (h234 * h34) )
606   i\_h34 = (h1234 * h234) * i\_denom
607   i\_h234 = (h1234 * h34) * i\_denom
608   i\_h1234 = (h234 * h34) * i\_denom
609 
610   \textcolor{comment}{! Calculation coefficients in the four equations}
611 
612   \textcolor{comment}{! The expressions for Csys(3) and Csys(4) come from reducing the 4x4 matrix problem into the following
       2x2}
613   \textcolor{comment}{! matrix problem, then manipulating the analytic solution to avoid any subtraction and simplifying.}
614   \textcolor{comment}{!  (C1\_3 * h123 * h23) * Csys(3) + (0.25 * h123 * h23 * (h3 + 2.0*h2 + 3.0*h1)) * Csys(4) =}
615   \textcolor{comment}{!            (u(3)-u(1)) - (u(2)-u(1)) * (h12 + h23) * I\_h12}
616   \textcolor{comment}{!  (C1\_3 * ((h23 + h34) * h1234 + h23 * h3)) * Csys(3) +}
617   \textcolor{comment}{!  (0.25 * ((h1234 + h123 + h12 + h1) * h23 * h3 + (h1234 + h12 + h1) * (h23 + h34) * h1234)) * Csys(4) =}
618   \textcolor{comment}{!            (u(4)-u(1)) - (u(2)-u(1)) * (h123 + h234) * I\_h12}
619   \textcolor{comment}{! The final expressions for Csys(1) and Csys(2) were derived by algebraically manipulating the following
       expressions:}
620   \textcolor{comment}{!  Csys(1) = (C1\_3 * h1 * h12 * Csys(3) + 0.25 * h1 * h12 * (2.0*h1+h2) * Csys(4)) - &}
621   \textcolor{comment}{!            (h1*I\_h12)*(u(2)-u(1)) + u(1)}
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)) + &}
623   \textcolor{comment}{!            2.0*I\_h12 * (u(2)-u(1))}
624   \textcolor{comment}{! These expressions are typically evaluated at x=0 and x=h1, so it is important that these are well
       behaved}
625   \textcolor{comment}{! for these values, suggesting that h1/h23 and h1/h34 should not be allowed to be too large.}
626 
627   wt(1,1) = -h1 * (i\_h1234 + i\_h123 + i\_h12)                              \textcolor{comment}{! > -3}
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))}
629   wt(3,1) = -h1 * h12 * h123 * i\_denom                                    \textcolor{comment}{! > -(h1/h34)*(1+(h1/h234))}
630 
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}
632   wt(2,2) = -2.0 * ((h1 * h12 * i\_h1234) *       (i\_h23 * (i\_h234 + i\_h123)) + &       \textcolor{comment}{! >
       -(10+6*(h1/h234))/h23}
633                     (h1+h12) * ( i\_h1234*i\_h234 + i\_h23 * (i\_h234 + i\_h123) ) )
634   wt(3,2) =  2.0 * ((h1+h12) * h123 + h1*h12 ) * i\_denom                               \textcolor{comment}{! < (2+(6*h1/h234))
       / h34}
635 
636   wt(1,3) = -3.0 * i\_h12 * i\_h123* ( 1.0 + i\_h1234 * ((h1+h12)+h123) )                 \textcolor{comment}{! > -12 / (h12*h123)}
637   wt(2,3) =  3.0 * i\_h23 * ( i\_h123 + i\_h1234 * ((h1+h12)+h123) * (i\_h123 + i\_h234) )  \textcolor{comment}{! < 12 / (h23^2)}
638   wt(3,3) = -3.0 * ((h1+h12)+h123) * i\_denom                                           \textcolor{comment}{! > -9 / (h234*h23)}
639 
640   wt(1,4) =  4.0 * i\_h1234 * i\_h123 * i\_h12                          \textcolor{comment}{! Wt*h1^3 < 4}
641   wt(2,4) = -4.0 * i\_h1234 * (i\_h23 * (i\_h123 + i\_h234))             \textcolor{comment}{! Wt*h1^3 > -4* (h1/h23)*(1+h1/h234)}
642   wt(3,4) =  4.0 * i\_denom  \textcolor{comment}{! = 4.0*I\_h1234 * I\_h234 * I\_h34         ! Wt*h1^3 < 4 * (h1/h234)*(h1/h34)}
643 
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))
645   csys(2) = (wt(1,2) * (u(2)-u(1)) + wt(2,2) * (u(3)-u(2))) + wt(3,2) * (u(4)-u(3))
646   csys(3) = (wt(1,3) * (u(2)-u(1)) + wt(2,3) * (u(3)-u(2))) + wt(3,3) * (u(4)-u(3))
647   csys(4) = (wt(1,4) * (u(2)-u(1)) + wt(2,4) * (u(3)-u(2))) + wt(3,4) * (u(4)-u(3))
648 
649   \textcolor{comment}{! endif ! End of non-uniform layer thickness branch.}
650 
651   \textcolor{comment}{! To verify that these answers are correct, uncomment the following:}
652 \textcolor{comment}{!  u\_mag = 0.0 ; do i=1,4 ; u\_mag = max(u\_mag, abs(u(i))) ; enddo}
653 \textcolor{comment}{!  do i = 1,4}
654 \textcolor{comment}{!    if (i==1) then ; zavg = 0.5*dz(i) ; else ; zavg = zavg + 0.5*(dz(i-1)+dz(i)) ; endif}
655 \textcolor{comment}{!    Atest(1) = 1.0}
656 \textcolor{comment}{!    Atest(2) = zavg                             ! = ( (z(i+1)**2) - (z(i)**2) ) / (2*dz(i))}
657 \textcolor{comment}{!    Atest(3) = (zavg**2 + 0.25*C1\_3*dz(i)**2)   ! = ( (z(i+1)**3) - (z(i)**3) ) / (3*dz(i))}
658 \textcolor{comment}{!    Atest(4) = zavg * (zavg**2 + 0.25*dz(i)**2) ! = ( (z(i+1)**4) - (z(i)**4) ) / (4*dz(i))}
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}
660 \textcolor{comment}{!    write(mesg, '("end\_value\_h4 line ", i2, " c\_mag = ", es10.2, " u\_mag = ", es10.2)') i, c\_mag, u\_mag}
661 \textcolor{comment}{!    call test\_line(mesg, 4, Atest, Csys, u(i), u\_mag*c\_mag, tol=1.0e-15)}
662 \textcolor{comment}{!  enddo}
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}}
\subsubsection{\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{\tt in}  & {\em n} & The size of the system\\
\hline
\mbox{\tt in,out}  & {\em a} & The matrix being inverted \mbox{[}nondim\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em r} & system right-\/hand side \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em x} & solution vector \mbox{[}A\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 1384 of file regrid\+\_\+edge\+\_\+values.\+F90.


\begin{DoxyCode}
1384   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{  !< The size of the system}
1385   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,N)}, \textcolor{keywordtype}{intent(inout)} :: a\textcolor{comment}{  !< The matrix being inverted [nondim]}
1386   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(inout)} :: r\textcolor{comment}{  !< system right-hand side [A]}
1387   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(inout)} :: x\textcolor{comment}{  !< solution vector [A]}
1388 
1389   \textcolor{comment}{! Local variables}
1390   \textcolor{keywordtype}{real}    :: factor       \textcolor{comment}{! The factor that eliminates the leading nonzero element in a row.}
1391   \textcolor{keywordtype}{real}    :: i\_pivot      \textcolor{comment}{! The reciprocal of the pivot value [inverse of the input units of a row of A]}
1392   \textcolor{keywordtype}{real}    :: swap
1393   \textcolor{keywordtype}{integer} :: i, j, k
1394 
1395   \textcolor{comment}{! Loop on rows to transform the problem into multiplication by an upper-right matrix.}
1396   \textcolor{keywordflow}{do} i=1,n-1
1397     \textcolor{comment}{! Seek a pivot for column i starting in row i, and continuing into the remaining rows.  If the}
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.}
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}
1400     \textcolor{keywordflow}{if} ( k > n ) \textcolor{keywordflow}{then}  \textcolor{comment}{! No pivot could be found and the system is singular.}
1401       \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{' A='},a
1402       \textcolor{keyword}{call }mom\_error( fatal, \textcolor{stringliteral}{'The linear system sent to linear\_solver is singular.'} )
1403 \textcolor{keywordflow}{    endif}
1404 
1405     \textcolor{comment}{! If the pivot is in a row that is different than row i, swap those two rows, noting that both}
1406     \textcolor{comment}{! rows start with i-1 zero values.}
1407     \textcolor{keywordflow}{if} ( k /= i ) \textcolor{keywordflow}{then}
1408       \textcolor{keywordflow}{do} j=i,n ; swap = a(j,i) ; a(j,i) = a(j,k) ; a(j,k) = swap ;\textcolor{keywordflow}{ enddo}
1409       swap = r(i) ; r(i) = r(k) ; r(k) = swap
1410 \textcolor{keywordflow}{    endif}
1411 
1412     \textcolor{comment}{! Transform the pivot to 1 by dividing the entire row (right-hand side included) by the pivot}
1413     i\_pivot = 1.0 / a(i,i)
1414     a(i,i) = 1.0
1415     \textcolor{keywordflow}{do} j=i+1,n ; a(j,i) = a(j,i) * i\_pivot ;\textcolor{keywordflow}{ enddo}
1416     r(i) = r(i) * i\_pivot
1417 
1418     \textcolor{comment}{! Put zeros in column for all rows below that contain the pivot (which is row i)}
1419     \textcolor{keywordflow}{do} k=i+1,n    \textcolor{comment}{! k is the row index}
1420       factor = a(i,k)
1421       \textcolor{comment}{! A(i,k) = 0.0  ! These elements are not used again, so this line can be skipped for speed.}
1422       \textcolor{keywordflow}{do} j=i+1,n ; a(j,k) = a(j,k) - factor * a(j,i) ;\textcolor{keywordflow}{ enddo}
1423       r(k) = r(k) - factor * r(i)
1424 \textcolor{keywordflow}{    enddo}
1425 
1426 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on i}
1427 
1428   \textcolor{comment}{! Solve the system by back substituting into what is now an upper-right matrix.}
1429   \textcolor{keywordflow}{if} (a(n,n) == 0.0) \textcolor{keywordflow}{then}  \textcolor{comment}{! No pivot could be found and the system is singular.}
1430     \textcolor{comment}{! write(0,*) ' A=',A}
1431     \textcolor{keyword}{call }mom\_error( fatal, \textcolor{stringliteral}{'The final pivot in linear\_solver is zero.'} )
1432 \textcolor{keywordflow}{  endif}
1433   x(n) = r(n) / a(n,n)  \textcolor{comment}{! The last row can now be solved trivially.}
1434   \textcolor{keywordflow}{do} i=n-1,1,-1 \textcolor{comment}{! loop on rows, starting from second to last row}
1435     x(i) = r(i)
1436     \textcolor{keywordflow}{do} j=i+1,n ; x(i) = x(i) - a(j,i) * x(j) ;\textcolor{keywordflow}{ enddo}
1437 \textcolor{keywordflow}{  enddo}
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}}
\subsubsection{\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{\tt in}  & {\em n} & The size of the system\\
\hline
\mbox{\tt in}  & {\em ac} & Matrix center diagonal offset from Al + Au\\
\hline
\mbox{\tt in}  & {\em al} & Matrix lower diagonal\\
\hline
\mbox{\tt in}  & {\em au} & Matrix upper diagonal\\
\hline
\mbox{\tt in}  & {\em r} & system right-\/hand side\\
\hline
\mbox{\tt out}  & {\em x} & solution vector \\
\hline
\end{DoxyParams}


Definition at line 1342 of file regrid\+\_\+edge\+\_\+values.\+F90.


\begin{DoxyCode}
1342   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)}  :: n\textcolor{comment}{   !< The size of the system}
1343   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: ac\textcolor{comment}{  !< Matrix center diagonal offset from Al + Au}
1344   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: al\textcolor{comment}{  !< Matrix lower diagonal}
1345   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: au\textcolor{comment}{  !< Matrix upper diagonal}
1346   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: r\textcolor{comment}{   !< system right-hand side}
1347   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(out)} :: x\textcolor{comment}{   !< solution vector}
1348   \textcolor{comment}{! Local variables}
1349   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)} :: c1       \textcolor{comment}{! Au / pivot for the backward sweep}
1350   \textcolor{keywordtype}{real}               :: d1       \textcolor{comment}{! The next value of 1.0 - c1}
1351   \textcolor{keywordtype}{real}               :: i\_pivot  \textcolor{comment}{! The inverse of the most recent pivot}
1352   \textcolor{keywordtype}{real}               :: denom\_t1 \textcolor{comment}{! The first term in the denominator of the inverse of the pivot.}
1353   \textcolor{keywordtype}{integer}            :: k        \textcolor{comment}{! Loop index}
1354 
1355   \textcolor{comment}{! Factorization and forward sweep, in a form that will never give a division by a}
1356   \textcolor{comment}{! zero pivot for positive definite Ac, Al, and Au.}
1357   i\_pivot = 1.0 / (ac(1) + au(1))
1358   d1 = ac(1) * i\_pivot
1359   c1(1) = au(1) * i\_pivot
1360   x(1) = r(1) * i\_pivot
1361   \textcolor{keywordflow}{do} k=2,n-1
1362     denom\_t1 = ac(k) + d1 * al(k)
1363     i\_pivot = 1.0 / (denom\_t1 + au(k))
1364     d1 = denom\_t1 * i\_pivot
1365     c1(k) = au(k) * i\_pivot
1366     x(k) = (r(k) - al(k) * x(k-1)) * i\_pivot
1367 \textcolor{keywordflow}{  enddo}
1368   i\_pivot = 1.0 / (ac(n) + d1 * al(n))
1369   x(n) = (r(n) - al(n) * x(n-1)) * i\_pivot
1370   \textcolor{comment}{! Backward sweep}
1371   \textcolor{keywordflow}{do} k=n-1,1,-1
1372     x(k) = x(k) - c1(k) * x(k+1)
1373 \textcolor{keywordflow}{  enddo}
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}}
\subsubsection{\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{\tt in}  & {\em mag} & The magnitude of leading order terms in this line\\
\hline
\mbox{\tt in}  & {\em n} & The number of points in the system\\
\hline
\mbox{\tt in}  & {\em a} & One of the two vectors being multiplied\\
\hline
\mbox{\tt in}  & {\em c} & One of the two vectors being multiplied\\
\hline
\mbox{\tt in}  & {\em r} & The expected solution of the equation\\
\hline
\mbox{\tt in}  & {\em msg} & An identifying message for this test\\
\hline
\mbox{\tt 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}
1445   \textcolor{keywordtype}{real},               \textcolor{keywordtype}{intent(in)} :: mag\textcolor{comment}{  !< The magnitude of leading order terms in this line}
1446   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)} :: n\textcolor{comment}{    !< The number of points in the system}
1447   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}, \textcolor{keywordtype}{intent(in)} :: a\textcolor{comment}{    !< One of the two vectors being multiplied}
1448   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(4)}, \textcolor{keywordtype}{intent(in)} :: c\textcolor{comment}{    !< One of the two vectors being multiplied}
1449   \textcolor{keywordtype}{real},               \textcolor{keywordtype}{intent(in)} :: r\textcolor{comment}{    !< The expected solution of the equation}
1450   \textcolor{keywordtype}{character(len=*)},   \textcolor{keywordtype}{intent(in)} :: msg\textcolor{comment}{  !< An identifying message for this test}
1451   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional},     \textcolor{keywordtype}{intent(in)} :: tol\textcolor{comment}{  !< The fractional tolerance for the two solutions}
1452 
1453   \textcolor{keywordtype}{real} :: sum, sum\_mag
1454   \textcolor{keywordtype}{real} :: tolerance
1455   \textcolor{keywordtype}{character(len=128)} :: mesg2
1456   \textcolor{keywordtype}{integer} :: i
1457 
1458   tolerance = 1.0e-12 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(tol)) tolerance = tol
1459 
1460   sum = 0.0 ; sum\_mag = max(0.0,mag)
1461   \textcolor{keywordflow}{do} i=1,n
1462     sum = sum + a(i) * c(i)
1463     sum\_mag = sum\_mag + abs(a(i) * c(i))
1464 \textcolor{keywordflow}{  enddo}
1465 
1466   \textcolor{keywordflow}{if} (abs(sum - r) > tolerance * (sum\_mag + abs(r))) \textcolor{keywordflow}{then}
1467     \textcolor{keyword}{write}(mesg2, \textcolor{stringliteral}{'(", Fractional error = ", es12.4,", sum = ", es12.4)'}) (sum - r) / (sum\_mag + abs(r)), 
      sum
1468     \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"Failed line test: "}//trim(msg)//trim(mesg2))
1469 \textcolor{keywordflow}{  endif}
1470 
\end{DoxyCode}
