|
Last change
on this file since 3470 was
3470,
checked in by evos, 4 weeks ago
|
|
we added the option to use NS dynamical subsurface ice in the model to more realisticly calculate the amount of ice in the subsurface and therfore the subsurface thermal inertia
|
|
File size:
842 bytes
|
| Line | |
|---|
| 1 | C================================================ |
|---|
| 2 | C Tridiagonal solver |
|---|
| 3 | C================================================ |
|---|
| 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 | c if(n.gt.NMAX) then |
|---|
| 14 | c print *, 'tridag: too many points, set NMAX>',n |
|---|
| 15 | c stop |
|---|
| 16 | c endif |
|---|
| 17 | bet=b(1) |
|---|
| 18 | u(1)=r(1)/bet |
|---|
| 19 | do 11 j=2,n |
|---|
| 20 | gam(j)=c(j-1)/bet |
|---|
| 21 | bet=b(j)-a(j)*gam(j) |
|---|
| 22 | c if(bet.eq.0.)pause 'tridag failed' |
|---|
| 23 | u(j)=(r(j)-a(j)*u(j-1))/bet |
|---|
| 24 | 11 continue |
|---|
| 25 | do 12 j=n-1,1,-1 |
|---|
| 26 | u(j)=u(j)-gam(j+1)*u(j+1) |
|---|
| 27 | 12 continue |
|---|
| 28 | return |
|---|
| 29 | END |
|---|
| 30 | C (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31(+. |
|---|
Note: See
TracBrowser
for help on using the repository browser.
