source: tags/nemo_dev_x3/NEMO/OPA_SRC/ZDF/zdf.matrixsolver.h90 @ 105

!! 
!!         ZDF.MATRIXSOLVER
!!       ********************
!! 
!! Matrix inversion
!!----------------------------------------------------------------------
!   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 2d arry zwt and zwz are work space arrays
!
!   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.
!!----------------------------------------------------------------------
!!  OPA 8.5, LODYC-IPSL (2002)
!!----------------------------------------------------------------------

      ikstp1 = ikst + 1
      ikenm2 = jpk - 2
      DO ji = 2, jpim1
         zwt(ji,ikst) = zwd(ji,ikst)
      END DO
      DO jk = ikstp1, jpkm1
         DO ji = 2, jpim1
            zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
         END DO
      END DO
      DO ji = 2, jpim1
         zwz(ji,ikst) = zwy(ji,ikst)
      END DO
      DO jk = ikstp1, jpkm1
         DO ji = 2, jpim1
            zwz(ji,jk) = zwy(ji,jk) - zwi(ji,jk) / zwt(ji,jk-1) * zwz(ji,jk-1)
         END DO
      END DO
      DO ji = 2, jpim1
         zwx(ji,jpkm1) = zwz(ji,jpkm1) / zwt(ji,jpkm1)
      END DO
      DO jk = ikenm2, ikst, -1
         DO ji = 2, jpim1
            zwx(ji,jk) =( zwz(ji,jk) - zws(ji,jk) * zwx(ji,jk+1) ) / zwt(ji,jk)
         END DO
      END DO