|
Last change
on this file since 3493 was
3493,
checked in by jbclement, 3 weeks ago
|
|
PEM:
- Renaming of Norbert Schorghofer's subroutines with the prefix 'NS_';
- Making the extension of all NS's subroutines as '.F90';
- Deletion of the wrapper subroutine;
- Making the initialization, variables management and arguments of the main subroutine for the dynamic computation of ice table to be more suitable.
JBC
|
|
File size:
831 bytes
|
| Line | |
|---|
| 1 | !================================================ |
|---|
| 2 | ! Tridiagonal solver |
|---|
| 3 | !================================================ |
|---|
| 4 | SUBROUTINE tridag(a,b,c,r,u,n) |
|---|
| 5 | INTEGER n,NMAX |
|---|
| 6 | REAL*8 a(n),b(n),c(n),r(n),u(n) |
|---|
| 7 | PARAMETER (NMAX=1000) |
|---|
| 8 | INTEGER j |
|---|
| 9 | REAL*8 bet,gam(NMAX) |
|---|
| 10 | if(b(1).eq.0.) then |
|---|
| 11 | stop 'tridag: rewrite equations' |
|---|
| 12 | endif |
|---|
| 13 | ! if(n.gt.NMAX) then |
|---|
| 14 | ! print *, 'tridag: too many points, set NMAX>',n |
|---|
| 15 | ! stop |
|---|
| 16 | ! endif |
|---|
| 17 | bet=b(1) |
|---|
| 18 | u(1)=r(1)/bet |
|---|
| 19 | do j=2,n |
|---|
| 20 | gam(j)=c(j-1)/bet |
|---|
| 21 | bet=b(j)-a(j)*gam(j) |
|---|
| 22 | ! if(bet.eq.0.)pause 'tridag failed' |
|---|
| 23 | u(j)=(r(j)-a(j)*u(j-1))/bet |
|---|
| 24 | enddo |
|---|
| 25 | do j=n-1,1,-1 |
|---|
| 26 | u(j)=u(j)-gam(j+1)*u(j+1) |
|---|
| 27 | enddo |
|---|
| 28 | ! return |
|---|
| 29 | END |
|---|
| 30 | ! (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31(+. |
|---|
Note: See
TracBrowser
for help on using the repository browser.
