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