Changeset 490 for trunk/NEMO/TOP_SRC/TRP/trczdf_imp.F90

Timestamp:

2006-09-01T15:50:43+02:00 (18 years ago)

Author:

opalod

Message:

nemo_v1_bugfix_053:CE+RB: enable use of ndttrc

File:

• trunk/NEMO/TOP_SRC/TRP/trczdf_imp.F90 (modified) (2 diffs)

trunk/NEMO/TOP_SRC/TRP/trczdf_imp.F90

@@ -97 +97 @@
          IF( neuler == 0 .AND. kt == nittrc000 ) THEN
             rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)             ! restarting with Euler time stepping
-         ELSEIF( kt <= nittrc000 + 1 ) THEN
+         ELSEIF( kt <= nittrc000 + ndttrc ) THEN
             rdttrc(:) = 2. * rdttra(:) * FLOAT(ndttrc)         ! leapfrog
          ENDIF
@@ -145 +145 @@
          END DO
          
- 
-   ! Matrix inversion from the first level
-   ikst = 1
-
-#   include "zdf.matrixsolver.vopt.h90"        
- 
+         
+         ! Matrix inversion from the first level
+         !----------------------------------------------------------------------
+         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
+         !
+         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
+         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
+         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
+         !        (        ...               )( ...  ) ( ...  )
+         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
+         !
+         !   m is decomposed in the product of an upper and lower triangular matrix
+         !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
+         !   The second member is in 2d array zwy
+         !   The solution is in 2d array zwx
+         !   The 3d arry zwt is a work space array
+         !   zwy is used and then used as a work space array : its value is modified!
+         !
+         !   N.B. the starting vertical index (ikst) is equal to 1 except for
+         !   the resolution of tke matrix where surface tke value is prescribed
+         !   so that ikstrt=2.
+        ikst = 1
+        ikstp1 = ikst + 1
+        ikenm2 = jpk - 2
+        
+        DO jj = 2, jpjm1
+           DO ji = fs_2, fs_jpim1
+              zwt(ji,jj,ikst) = zwd(ji,jj,ikst)
+           END DO
+        END DO
+        
+        DO jk = ikstp1, jpkm1
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1
+                 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
+              END DO
+           END DO
+        END DO
+        
+        DO jk = ikstp1, jpkm1
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1
+                 zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * zwy(ji,jj,jk-1)
+              END DO
+           END DO
+        END DO
+        
+        DO jj = 2, jpjm1
+           DO ji = fs_2, fs_jpim1
+              zwx(ji,jj,jpkm1) = zwy(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
+           END DO
+        END DO
+        
+        DO jk = ikenm2, ikst, -1
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1
+                 zwx(ji,jj,jk) = ( zwy(ji,jj,jk) - zws(ji,jj,jk) * zwx(ji,jj,jk+1) ) / zwt(ji,jj,jk)
+              END DO
+           END DO
+        END DO
+
          
 #if defined key_trc_diatrd