!! 
   !!         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 9.0 , LOCEAN-IPSL  (2005)
   !!   $Header$
   !!   This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
   !!----------------------------------------------------------------------

      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