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
|
Rev | Line | |
---|
[3493] | 1 | !================================================ |
---|
| 2 | ! Tridiagonal solver |
---|
| 3 | !================================================ |
---|
[3470] | 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 |
---|
[3493] | 13 | ! if(n.gt.NMAX) then |
---|
| 14 | ! print *, 'tridag: too many points, set NMAX>',n |
---|
| 15 | ! stop |
---|
| 16 | ! endif |
---|
[3470] | 17 | bet=b(1) |
---|
| 18 | u(1)=r(1)/bet |
---|
[3493] | 19 | do j=2,n |
---|
[3470] | 20 | gam(j)=c(j-1)/bet |
---|
| 21 | bet=b(j)-a(j)*gam(j) |
---|
[3493] | 22 | ! if(bet.eq.0.)pause 'tridag failed' |
---|
[3470] | 23 | u(j)=(r(j)-a(j)*u(j-1))/bet |
---|
[3493] | 24 | enddo |
---|
| 25 | do j=n-1,1,-1 |
---|
[3470] | 26 | u(j)=u(j)-gam(j+1)*u(j+1) |
---|
[3493] | 27 | enddo |
---|
| 28 | ! return |
---|
[3470] | 29 | END |
---|
[3493] | 30 | ! (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31(+. |
---|
Note: See
TracBrowser
for help on using the repository browser.