New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdf.matrixsolver.h90 in tags/nemo_v3_2/nemo_v3_2/NEMO/OFF_SRC/ZDF – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OFF_SRC/ZDF/zdf.matrixsolver.h90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm

File size: 2.2 KB
Line 
1   !!
2   !!         ZDF.MATRIXSOLVER
3   !!       ********************
4   !!
5   !! Matrix inversion
6   !!----------------------------------------------------------------------
7   !!   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
8   !!
9   !!        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
10   !!        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
11   !!        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
12   !!        (        ...               )( ...  ) ( ...  )
13   !!        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
14   !!
15   !!   m is decomposed in the product of an upper and lower triangular
16   !!   matrix
17   !!   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
18   !!   The second member is in 2d array zwy
19   !!   The solution is in 2d array zwx
20   !!   The 2d arry zwt and zwz are work space arrays
21   !!
22   !!   N.B. the starting vertical index (ikst) is equal to 1 except for
23   !!   the resolution of tke matrix where surface tke value is prescribed
24   !!   so that ikstrt=2.
25   !!----------------------------------------------------------------------
26   !!----------------------------------------------------------------------
27   !!   OPA 9.0 , LOCEAN-IPSL  (2005)
28   !!   $Id: zdf.matrixsolver.h90 1152 2008-06-26 14:11:13Z rblod $
29   !!   This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
30   !!----------------------------------------------------------------------
31
32      ikstp1 = ikst + 1
33      ikenm2 = jpk - 2
34      DO ji = 2, jpim1
35         zwt(ji,ikst) = zwd(ji,ikst)
36      END DO
37      DO jk = ikstp1, jpkm1
38         DO ji = 2, jpim1
39            zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
40         END DO
41      END DO
42      DO ji = 2, jpim1
43         zwz(ji,ikst) = zwy(ji,ikst)
44      END DO
45      DO jk = ikstp1, jpkm1
46         DO ji = 2, jpim1
47            zwz(ji,jk) = zwy(ji,jk) - zwi(ji,jk) / zwt(ji,jk-1) * zwz(ji,jk-1)
48         END DO
49      END DO
50      DO ji = 2, jpim1
51         zwx(ji,jpkm1) = zwz(ji,jpkm1) / zwt(ji,jpkm1)
52      END DO
53      DO jk = ikenm2, ikst, -1
54         DO ji = 2, jpim1
55            zwx(ji,jk) =( zwz(ji,jk) - zws(ji,jk) * zwx(ji,jk+1) ) / zwt(ji,jk)
56         END DO
57      END DO
58     
Note: See TracBrowser for help on using the repository browser.