\hypertarget{namespaceregrid__solvers}{}\section{regrid\+\_\+solvers Module Reference} \label{namespaceregrid__solvers}\index{regrid\+\_\+solvers@{regrid\+\_\+solvers}} \subsection{Detailed Description} Solvers of linear systems. Date of creation\+: 2008.\+06.\+12 L. White. This module contains solvers of linear systems. These routines have now been updated for greater efficiency, especially in special cases. \subsection*{Functions/\+Subroutines} \begin{DoxyCompactItemize} \item subroutine, public \hyperlink{namespaceregrid__solvers_a36827238948de49ff0dd0be54cfaf719}{solve\+\_\+linear\+\_\+system} (A, R, X, N, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Solve the linear system AX = R by Gaussian elimination. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespaceregrid__solvers_a2ed09d1e0857610eb3638601e20d2a2c}{linear\+\_\+solver} (N, A, R, X) \begin{DoxyCompactList}\small\item\em Solve the linear system AX = R by Gaussian elimination. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespaceregrid__solvers_aac4382af38975d9cfcfd6b00adafaeab}{solve\+\_\+tridiagonal\+\_\+system} (Al, Ad, Au, R, X, N, answers\+\_\+2018) \begin{DoxyCompactList}\small\item\em Solve the tridiagonal system AX = R. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespaceregrid__solvers_a3f1d27aaab0a19c9abd07c96bb704bb6}{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag} (Al, Ac, Au, R, X, N) \begin{DoxyCompactList}\small\item\em Solve the tridiagonal system AX = R. \end{DoxyCompactList}\end{DoxyCompactItemize} \subsection{Function/\+Subroutine Documentation} \mbox{\Hypertarget{namespaceregrid__solvers_a2ed09d1e0857610eb3638601e20d2a2c}\label{namespaceregrid__solvers_a2ed09d1e0857610eb3638601e20d2a2c}} \index{regrid\+\_\+solvers@{regrid\+\_\+solvers}!linear\+\_\+solver@{linear\+\_\+solver}} \index{linear\+\_\+solver@{linear\+\_\+solver}!regrid\+\_\+solvers@{regrid\+\_\+solvers}} \subsubsection{\texorpdfstring{linear\+\_\+solver()}{linear\_solver()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+solvers\+::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})} 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 112 of file regrid\+\_\+solvers.\+F90. \begin{DoxyCode} 112 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: n\textcolor{comment}{ !< The size of the system} 113 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,N)}, \textcolor{keywordtype}{intent(inout)} :: a\textcolor{comment}{ !< The matrix being inverted [nondim]} 114 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(inout)} :: r\textcolor{comment}{ !< system right-hand side [A]} 115 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(inout)} :: x\textcolor{comment}{ !< solution vector [A]} 116 117 \textcolor{comment}{! Local variables} 118 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: eps = 0.0 \textcolor{comment}{! Minimum pivot magnitude allowed} 119 \textcolor{keywordtype}{real} :: factor \textcolor{comment}{! The factor that eliminates the leading nonzero element in a row.} 120 \textcolor{keywordtype}{real} :: i\_pivot \textcolor{comment}{! The reciprocal of the pivot value [inverse of the input units of a row of A]} 121 \textcolor{keywordtype}{real} :: swap 122 \textcolor{keywordtype}{logical} :: found\_pivot \textcolor{comment}{! If true, a pivot has been found} 123 \textcolor{keywordtype}{integer} :: i, j, k 124 125 \textcolor{comment}{! Loop on rows to transform the problem into multiplication by an upper-right matrix.} 126 \textcolor{keywordflow}{do} i=1,n-1 127 \textcolor{comment}{! Seek a pivot for column i starting in row i, and continuing into the remaining rows. If the} 128 \textcolor{comment}{! pivot is in a row other than i, swap them. If no valid pivot is found, i = N+1 after this loop.} 129 \textcolor{keywordflow}{do} k=i,n ; \textcolor{keywordflow}{if} ( abs( a(i,k) ) > eps ) \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop to find pivot} 130 \textcolor{keywordflow}{if} ( k > n ) \textcolor{keywordflow}{then} \textcolor{comment}{! No pivot could be found and the system is singular.} 131 \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{' A='},a 132 \textcolor{keyword}{call }mom\_error( fatal, \textcolor{stringliteral}{'The linear system is singular !'} ) 133 \textcolor{keywordflow}{ endif} 134 135 \textcolor{comment}{! If the pivot is in a row that is different than row i, swap those two rows, noting that both} 136 \textcolor{comment}{! rows start with i-1 zero values.} 137 \textcolor{keywordflow}{if} ( k /= i ) \textcolor{keywordflow}{then} 138 \textcolor{keywordflow}{do} j=i,n ; swap = a(j,i) ; a(j,i) = a(j,k) ; a(j,k) = swap ;\textcolor{keywordflow}{ enddo} 139 swap = r(i) ; r(i) = r(k) ; r(k) = swap 140 \textcolor{keywordflow}{ endif} 141 142 \textcolor{comment}{! Transform the pivot to 1 by dividing the entire row (right-hand side included) by the pivot} 143 i\_pivot = 1.0 / a(i,i) 144 a(i,i) = 1.0 145 \textcolor{keywordflow}{do} j=i+1,n ; a(j,i) = a(j,i) * i\_pivot ;\textcolor{keywordflow}{ enddo} 146 r(i) = r(i) * i\_pivot 147 148 \textcolor{comment}{! Put zeros in column for all rows below that contain the pivot (which is row i)} 149 \textcolor{keywordflow}{do} k=i+1,n \textcolor{comment}{! k is the row index} 150 factor = a(i,k) 151 \textcolor{comment}{! A(i,k) = 0.0 ! These elements are not used again, so this line can be skipped for speed.} 152 \textcolor{keywordflow}{do} j=i+1,n ; a(j,k) = a(j,k) - factor * a(j,i) ;\textcolor{keywordflow}{ enddo} 153 r(k) = r(k) - factor * r(i) 154 \textcolor{keywordflow}{ enddo} 155 156 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on i} 157 158 \textcolor{comment}{! Solve the system by back substituting into what is now an upper-right matrix.} 159 x(n) = r(n) / a(n,n) \textcolor{comment}{! The last row is now trivially solved.} 160 \textcolor{keywordflow}{do} i=n-1,1,-1 \textcolor{comment}{! loop on rows, starting from second to last row} 161 x(i) = r(i) 162 \textcolor{keywordflow}{do} j=i+1,n ; x(i) = x(i) - a(j,i) * x(j) ;\textcolor{keywordflow}{ enddo} 163 \textcolor{keywordflow}{ enddo} 164 \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__solvers_a3f1d27aaab0a19c9abd07c96bb704bb6}\label{namespaceregrid__solvers_a3f1d27aaab0a19c9abd07c96bb704bb6}} \index{regrid\+\_\+solvers@{regrid\+\_\+solvers}!solve\+\_\+diag\+\_\+dominant\+\_\+tridiag@{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag}} \index{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag@{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag}!regrid\+\_\+solvers@{regrid\+\_\+solvers}} \subsubsection{\texorpdfstring{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag()}{solve\_diag\_dominant\_tridiag()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+solvers\+::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})} 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 235 of file regrid\+\_\+solvers.\+F90. \begin{DoxyCode} 235 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: n\textcolor{comment}{ !< The size of the system} 236 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: ac\textcolor{comment}{ !< Matrix center diagonal offset from Al + Au} 237 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: al\textcolor{comment}{ !< Matrix lower diagonal} 238 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: au\textcolor{comment}{ !< Matrix upper diagonal} 239 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: r\textcolor{comment}{ !< system right-hand side} 240 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(out)} :: x\textcolor{comment}{ !< solution vector} 241 \textcolor{comment}{! Local variables} 242 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)} :: c1 \textcolor{comment}{! Au / pivot for the backward sweep} 243 \textcolor{keywordtype}{real} :: d1 \textcolor{comment}{! The next value of 1.0 - c1} 244 \textcolor{keywordtype}{real} :: i\_pivot \textcolor{comment}{! The inverse of the most recent pivot} 245 \textcolor{keywordtype}{real} :: denom\_t1 \textcolor{comment}{! The first term in the denominator of the inverse of the pivot.} 246 \textcolor{keywordtype}{integer} :: k \textcolor{comment}{! Loop index} 247 248 \textcolor{comment}{! Factorization and forward sweep, in a form that will never give a division by a} 249 \textcolor{comment}{! zero pivot for positive definite Ac, Al, and Au.} 250 i\_pivot = 1.0 / (ac(1) + au(1)) 251 d1 = ac(1) * i\_pivot 252 c1(1) = au(1) * i\_pivot 253 x(1) = r(1) * i\_pivot 254 \textcolor{keywordflow}{do} k=2,n-1 255 denom\_t1 = ac(k) + d1 * al(k) 256 i\_pivot = 1.0 / (denom\_t1 + au(k)) 257 d1 = denom\_t1 * i\_pivot 258 c1(k) = au(k) * i\_pivot 259 x(k) = (r(k) - al(k) * x(k-1)) * i\_pivot 260 \textcolor{keywordflow}{ enddo} 261 i\_pivot = 1.0 / (ac(n) + d1 * al(n)) 262 x(n) = (r(n) - al(n) * x(n-1)) * i\_pivot 263 \textcolor{comment}{! Backward sweep} 264 \textcolor{keywordflow}{do} k=n-1,1,-1 265 x(k) = x(k) - c1(k) * x(k+1) 266 \textcolor{keywordflow}{ enddo} 267 \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__solvers_a36827238948de49ff0dd0be54cfaf719}\label{namespaceregrid__solvers_a36827238948de49ff0dd0be54cfaf719}} \index{regrid\+\_\+solvers@{regrid\+\_\+solvers}!solve\+\_\+linear\+\_\+system@{solve\+\_\+linear\+\_\+system}} \index{solve\+\_\+linear\+\_\+system@{solve\+\_\+linear\+\_\+system}!regrid\+\_\+solvers@{regrid\+\_\+solvers}} \subsubsection{\texorpdfstring{solve\+\_\+linear\+\_\+system()}{solve\_linear\_system()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+solvers\+::solve\+\_\+linear\+\_\+system (\begin{DoxyParamCaption}\item[{real, dimension(n,n), intent(inout)}]{A, }\item[{real, dimension(n), intent(inout)}]{R, }\item[{real, dimension(n), intent(inout)}]{X, }\item[{integer, intent(in)}]{N, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})} 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 yields the answer. The matrix A must be square, with the first index varing down the column. \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 \mbox{\tt in} & {\em answers\+\_\+2018} & If true or absent use older, less efficient expressions. \\ \hline \end{DoxyParams} Definition at line 20 of file regrid\+\_\+solvers.\+F90. \begin{DoxyCode} 20 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: n\textcolor{comment}{ !< The size of the system} 21 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,N)}, \textcolor{keywordtype}{intent(inout)} :: a\textcolor{comment}{ !< The matrix being inverted [nondim]} 22 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(inout)} :: r\textcolor{comment}{ !< system right-hand side [A]} 23 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(inout)} :: x\textcolor{comment}{ !< solution vector [A]} 24 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true or absent use older, less efficient expressions.} 25 \textcolor{comment}{! Local variables} 26 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: eps = 0.0 \textcolor{comment}{! Minimum pivot magnitude allowed} 27 \textcolor{keywordtype}{real} :: factor \textcolor{comment}{! The factor that eliminates the leading nonzero element in a row.} 28 \textcolor{keywordtype}{real} :: pivot, i\_pivot \textcolor{comment}{! The pivot value and its reciprocal [nondim]} 29 \textcolor{keywordtype}{real} :: swap\_a, swap\_b 30 \textcolor{keywordtype}{logical} :: found\_pivot \textcolor{comment}{! If true, a pivot has been found} 31 \textcolor{keywordtype}{logical} :: old\_answers \textcolor{comment}{! If true, use expressions that give the original (2008 through 2018) MOM6 answers} 32 \textcolor{keywordtype}{integer} :: i, j, k 33 34 old\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) old\_answers = answers\_2018 35 36 \textcolor{comment}{! Loop on rows to transform the problem into multiplication by an upper-right matrix.} 37 \textcolor{keywordflow}{do} i = 1,n-1 38 39 40 \textcolor{comment}{! Start to look for a pivot in the current row, i. If the pivot in row i is not valid,} 41 \textcolor{comment}{! keep looking for a valid pivot by searching the entries of column i in rows below row i.} 42 \textcolor{comment}{! Once a valid pivot is found (say in row k), rows i and k are swaped.} 43 found\_pivot = .false. 44 k = i 45 \textcolor{keywordflow}{do} \textcolor{keywordflow}{while} ( ( .NOT. found\_pivot ) .AND. ( k <= n ) ) 46 \textcolor{keywordflow}{if} ( abs( a(k,i) ) > eps ) \textcolor{keywordflow}{then} \textcolor{comment}{! A valid pivot has been found} 47 found\_pivot = .true. 48 \textcolor{keywordflow}{else} \textcolor{comment}{! Seek a valid pivot in the next row} 49 k = k + 1 50 \textcolor{keywordflow}{ endif} 51 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop to find pivot} 52 53 \textcolor{comment}{! If no pivot could be found, the system is singular.} 54 \textcolor{keywordflow}{if} ( .NOT. found\_pivot ) \textcolor{keywordflow}{then} 55 \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{' A='},a 56 \textcolor{keyword}{call }mom\_error( fatal, \textcolor{stringliteral}{'The linear system is singular !'} ) 57 \textcolor{keywordflow}{ endif} 58 59 \textcolor{comment}{! If the pivot is in a row that is different than row i, that is if} 60 \textcolor{comment}{! k is different than i, we need to swap those two rows} 61 \textcolor{keywordflow}{if} ( k /= i ) \textcolor{keywordflow}{then} 62 \textcolor{keywordflow}{do} j = 1,n 63 swap\_a = a(i,j) ; a(i,j) = a(k,j) ; a(k,j) = swap\_a 64 \textcolor{keywordflow}{ enddo} 65 swap\_b = r(i) ; r(i) = r(k) ; r(k) = swap\_b 66 \textcolor{keywordflow}{ endif} 67 68 \textcolor{comment}{! Transform pivot to 1 by dividing the entire row (right-hand side included) by the pivot} 69 \textcolor{keywordflow}{if} (old\_answers) \textcolor{keywordflow}{then} 70 pivot = a(i,i) 71 \textcolor{keywordflow}{do} j = i,n ; a(i,j) = a(i,j) / pivot ;\textcolor{keywordflow}{ enddo} 72 r(i) = r(i) / pivot 73 \textcolor{keywordflow}{else} 74 i\_pivot = 1.0 / a(i,i) 75 a(i,i) = 1.0 76 \textcolor{keywordflow}{do} j = i+1,n ; a(i,j) = a(i,j) * i\_pivot ;\textcolor{keywordflow}{ enddo} 77 r(i) = r(i) * i\_pivot 78 \textcolor{keywordflow}{ endif} 79 80 \textcolor{comment}{! #INV: At this point, A(i,i) is a suitable pivot and it is equal to 1} 81 82 \textcolor{comment}{! Put zeros in column for all rows below that contain the pivot (which is row i)} 83 \textcolor{keywordflow}{do} k = i+1,n \textcolor{comment}{! k is the row index} 84 factor = a(k,i) 85 \textcolor{comment}{! A(k,i) = 0.0 ! These elements are not used again, so this line can be skipped for speed.} 86 \textcolor{keywordflow}{do} j = i+1,n \textcolor{comment}{! j is the column index} 87 a(k,j) = a(k,j) - factor * a(i,j) 88 \textcolor{keywordflow}{ enddo} 89 r(k) = r(k) - factor * r(i) 90 \textcolor{keywordflow}{ enddo} 91 92 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop on i} 93 94 \textcolor{comment}{! Solve system by back substituting in what is now an upper-right matrix.} 95 x(n) = r(n) / a(n,n) \textcolor{comment}{! The last row is now trivially solved.} 96 \textcolor{keywordflow}{do} i = n-1,1,-1 \textcolor{comment}{! loop on rows, starting from second to last row} 97 x(i) = r(i) 98 \textcolor{keywordflow}{do} j = i+1,n 99 x(i) = x(i) - a(i,j) * x(j) 100 \textcolor{keywordflow}{ enddo} 101 \textcolor{keywordflow}{if} (old\_answers) x(i) = x(i) / a(i,i) 102 \textcolor{keywordflow}{ enddo} 103 \end{DoxyCode} \mbox{\Hypertarget{namespaceregrid__solvers_aac4382af38975d9cfcfd6b00adafaeab}\label{namespaceregrid__solvers_aac4382af38975d9cfcfd6b00adafaeab}} \index{regrid\+\_\+solvers@{regrid\+\_\+solvers}!solve\+\_\+tridiagonal\+\_\+system@{solve\+\_\+tridiagonal\+\_\+system}} \index{solve\+\_\+tridiagonal\+\_\+system@{solve\+\_\+tridiagonal\+\_\+system}!regrid\+\_\+solvers@{regrid\+\_\+solvers}} \subsubsection{\texorpdfstring{solve\+\_\+tridiagonal\+\_\+system()}{solve\_tridiagonal\_system()}} {\footnotesize\ttfamily subroutine, public regrid\+\_\+solvers\+::solve\+\_\+tridiagonal\+\_\+system (\begin{DoxyParamCaption}\item[{real, dimension(n), intent(in)}]{Al, }\item[{real, dimension(n), intent(in)}]{Ad, }\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, }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})} Solve the tridiagonal system AX = R. This routine uses Thomas\textquotesingle{}s algorithm to solve the tridiagonal system AX = R. (A is made up of lower, middle and upper diagonals) \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em n} & The size of the system\\ \hline \mbox{\tt in} & {\em ad} & Matrix center diagonal\\ \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 \mbox{\tt in} & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\ \hline \end{DoxyParams} Definition at line 173 of file regrid\+\_\+solvers.\+F90. \begin{DoxyCode} 173 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: n\textcolor{comment}{ !< The size of the system} 174 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: ad\textcolor{comment}{ !< Matrix center diagonal} 175 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: al\textcolor{comment}{ !< Matrix lower diagonal} 176 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: au\textcolor{comment}{ !< Matrix upper diagonal} 177 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)} :: r\textcolor{comment}{ !< system right-hand side} 178 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(out)} :: x\textcolor{comment}{ !< solution vector} 179 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.} 180 \textcolor{comment}{! Local variables} 181 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)} :: pivot, al\_piv 182 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)} :: c1 \textcolor{comment}{! Au / pivot for the backward sweep} 183 \textcolor{keywordtype}{real} :: i\_pivot \textcolor{comment}{! The inverse of the most recent pivot} 184 \textcolor{keywordtype}{integer} :: k \textcolor{comment}{! Loop index} 185 \textcolor{keywordtype}{logical} :: old\_answers \textcolor{comment}{! If true, use expressions that give the original (2008 through 2018) MOM6 answers} 186 187 old\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) old\_answers = answers\_2018 188 189 \textcolor{keywordflow}{if} (old\_answers) \textcolor{keywordflow}{then} 190 \textcolor{comment}{! This version gives the same answers as the original (2008 through 2018) MOM6 code} 191 \textcolor{comment}{! Factorization and forward sweep} 192 pivot(1) = ad(1) 193 x(1) = r(1) 194 \textcolor{keywordflow}{do} k = 2,n 195 al\_piv(k) = al(k) / pivot(k-1) 196 pivot(k) = ad(k) - al\_piv(k) * au(k-1) 197 x(k) = r(k) - al\_piv(k) * x(k-1) 198 \textcolor{keywordflow}{ enddo} 199 200 \textcolor{comment}{! Backward sweep} 201 x(n) = r(n) / pivot(n) \textcolor{comment}{! This should be X(N) / pivot(N), but is OK if Al(N) = 0.} 202 \textcolor{keywordflow}{do} k = n-1,1,-1 203 x(k) = ( x(k) - au(k)*x(k+1) ) / pivot(k) 204 \textcolor{keywordflow}{ enddo} 205 \textcolor{keywordflow}{else} 206 \textcolor{comment}{! This is a more typical implementation of a tridiagonal solver than the one above.} 207 \textcolor{comment}{! It is mathematically equivalent but differs at roundoff, which can cascade up to larger values.} 208 209 \textcolor{comment}{! Factorization and forward sweep} 210 i\_pivot = 1.0 / ad(1) 211 x(1) = r(1) * i\_pivot 212 \textcolor{keywordflow}{do} k = 2,n 213 c1(k-1) = au(k-1) * i\_pivot 214 i\_pivot = 1.0 / (ad(k) - al(k) * c1(k-1)) 215 x(k) = (r(k) - al(k) * x(k-1)) * i\_pivot 216 \textcolor{keywordflow}{ enddo} 217 \textcolor{comment}{! Backward sweep} 218 \textcolor{keywordflow}{do} k = n-1,1,-1 219 x(k) = x(k) - c1(k) * x(k+1) 220 \textcolor{keywordflow}{ enddo} 221 222 \textcolor{keywordflow}{ endif} 223 \end{DoxyCode}