Index: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 2023)
+++ branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 2024)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                 ***  MODULE  trazdf_imp  ***
-   !! Ocean active tracers:  vertical component of the tracer mixing trend
+   !! Ocean  tracers:  vertical component of the tracer mixing trend
    !!======================================================================
    !! History :  OPA  !  1990-10  (B. Blanke)  Original code
@@ -14,4 +14,5 @@
    !!            2.0  !  2006-11  (G. Madec) New step reorganisation
    !!            3.2  !  2009-03  (G. Madec)  heat and salt content trends
+   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
    !!----------------------------------------------------------------------
   
@@ -25,10 +26,7 @@
    USE ldftra_oce      ! ocean active tracers: lateral physics
    USE ldfslp          ! lateral physics: slope of diffusion
-   USE trdmod          ! ocean active tracers trends 
-   USE trdmod_oce      ! ocean variables trends
    USE zdfddm          ! ocean vertical physics: double diffusion
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl          ! Print control
    USE domvvl          ! variable volume
    USE ldftra          ! lateral mixing type
@@ -50,6 +48,7 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   
-   SUBROUTINE tra_zdf_imp( kt, p2dt )
+ 
+   SUBROUTINE tra_zdf_imp( kt    , cdtype, p2dt,    &
+      &                    ptrab , ptraa , kjpt     )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_zdf_imp  ***
@@ -71,5 +70,5 @@
       !!                  associated with the lateral mixing, through the
       !!                  update of avt)
-      !!      The vertical diffusion of tracers (t & s) is given by:
+      !!      The vertical diffusion of the tracer t  is given by:
       !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
       !!      It is computed using a backward time scheme (t=ta).
@@ -78,25 +77,31 @@
       !!      Add this trend to the general trend ta,sa :
       !!         ta = ta + dz( avt dz(t) )
-      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
+      !!         if lk_zdfddm=T, use avs for salinity or for passive tracers
+      !!         (sa = sa + dz( avs dz(t) ) 
       !!
       !!      Third part: recover avt resulting from the vertical physics
       !!      ==========  alone, for further diagnostics (for example to
       !!                  compute the turbocline depth in zdfmxl.F90).
-      !!         avt = zavt
+      !!         if lk_zdfddm=T, use avt = zavt
       !!         (avs = zavs if lk_zdfddm=T )
       !!
-      !! ** Action  : - Update (ta,sa) with before vertical diffusion trend
+      !! ** Action  : - Update (ta) with before vertical diffusion trend
       !!
       !!---------------------------------------------------------------------
+      !! * Modules used
       USE oce    , ONLY :   zwd   => ua   ! ua used as workspace
       USE oce    , ONLY :   zws   => va   ! va  -          -
-      !!
-      INTEGER                 , INTENT(in) ::   kt     ! ocean time-step index
-      REAL(wp), DIMENSION(jpk), INTENT(in) ::   p2dt   ! vertical profile of tracer time-step
-      !!
-      INTEGER  ::   ji, jj, jk            ! dummy loop indices
-      REAL(wp) ::   zavi, zrhs, znvvl     ! temporary scalars
-      REAL(wp) ::   ze3tb, ze3tn, ze3ta   ! variable vertical scale factors
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi, zwt, zavsi   ! workspace arrays
+      !! * Arguments
+      INTEGER         , INTENT(in   )                                ::   kt             ! ocean time-step index
+      CHARACTER(len=3), INTENT(in   )                                ::   cdtype         ! =TRA or TRC (tracer indicator)
+      INTEGER         , INTENT(in   )                                ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)                ::   p2dt        ! vertical profile of tracer time-step
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)   ::   ptrab          ! before and now tracer fields
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)   ::   ptraa          ! tracer trend 
+      !!
+      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices
+      REAL(wp) ::  zavi, zrhs, znvvl     ! temporary scalars
+      REAL(wp) ::  ze3tb, ze3tn, ze3ta   ! variable vertical scale factors
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi, zwt   ! workspace arrays
       !!---------------------------------------------------------------------
 
@@ -107,14 +112,12 @@
          zavi = 0.e0      ! avoid warning at compilation phase when lk_ldfslp=F
       ENDIF
-
+      !
       ! I. Local initialization
       ! -----------------------
-      zwd  (1,:, : ) = 0.e0     ;     zwd  (jpi,:,:) = 0.e0
-      zws  (1,:, : ) = 0.e0     ;     zws  (jpi,:,:) = 0.e0
-      zwi  (1,:, : ) = 0.e0     ;     zwi  (jpi,:,:) = 0.e0
-      zwt  (1,:, : ) = 0.e0     ;     zwt  (jpi,:,:) = 0.e0
-      zavsi(1,:, : ) = 0.e0     ;     zavsi(jpi,:,:) = 0.e0
-      zwt  (:,:,jpk) = 0.e0     ;     zwt  ( : ,:,1) = 0.e0
-      zavsi(:,:,jpk) = 0.e0     ;     zavsi( : ,:,1) = 0.e0
+      zwd(1,:, : ) = 0.e0     ;     zwd(jpi,:,:) = 0.e0
+      zws(1,:, : ) = 0.e0     ;     zws(jpi,:,:) = 0.e0
+      zwi(1,:, : ) = 0.e0     ;     zwi(jpi,:,:) = 0.e0
+      zwt(1,:, : ) = 0.e0     ;     zwt(jpi,:,:) = 0.e0
+      zwt(:,:,jpk) = 0.e0     ;     zwt( : ,:,1) = 0.e0
 
       ! I.1 Variable volume : to take into account vertical variable vertical scale factors
@@ -130,233 +133,167 @@
       !     dk[ avt dk[ (t,s) ] ] diffusive trends
 
-
+      !
       ! II.0 Matrix construction
       ! ------------------------
-
+      DO jn = 1, kjpt
+         !
+         !  Matrix construction
+         ! ------------------------
+         IF( cdtype == 'TRA' .AND. jn == jp_tem )  THEN 
 #if defined key_ldfslp
-      ! update and save of avt (and avs if double diffusive mixing)
-      IF( l_traldf_rot ) THEN
-         DO jk = 2, jpkm1
+            ! update and save of avt (and avs if double diffusive mixing)
+            IF( l_traldf_rot ) THEN
+               DO jk = 2, jpkm1
+                  DO jj = 2, jpjm1
+                     DO ji = fs_2, fs_jpim1   ! vector opt.
+                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
+                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
+                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
+                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
+                     END DO
+                  END DO
+               END DO
+            ELSE                         ! no rotation but key_ldfslp defined
+               zwt  (:,:,:) = avt(:,:,:)
+            ENDIF
+#else
+            ! No isopycnal diffusion
+            zwt(:,:,:) = avt(:,:,:)           
+#endif
+            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
+                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
+                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
+                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
+                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
+                 END DO
+               END DO
+            END DO
+            ! Surface boudary conditions
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
-                     & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
-                     &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
-                  zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
-# if defined key_zdfddm
-                  zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi          ! dd mixing: zavsi = total vertical mixing coef. on salinity
-# endif
-               END DO
-            END DO
-         END DO
-      ELSE                         ! no rotation but key_ldfslp defined
-         zwt  (:,:,:) = avt(:,:,:)
-# if defined key_zdfddm
-         zavsi(:,:,:) = avs(:,:,:)       ! avs /= avt (double diffusion mixing)
-# endif
-      ENDIF   
+                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
+                 zwi(ji,jj,1) = 0.e0
+                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
+               END DO
+            END DO
+            !
+         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
+#if defined key_ldfslp
+            ! update and save of avt (and avs if double diffusive mixing)
+            IF( l_traldf_rot ) THEN
+               DO jk = 2, jpkm1
+                  DO jj = 2, jpjm1
+                     DO ji = fs_2, fs_jpim1   ! vector opt.
+                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
+                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
+                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
+                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
+                     END DO
+                  END DO
+               END DO
+            ELSE                         ! no rotation but key_ldfslp defined
+               zwt  (:,:,:) = fsavs(:,:,:)
+            ENDIF
 #else
-      ! No isopycnal diffusion
-      zwt(:,:,:) = avt(:,:,:)
-# if defined key_zdfddm
-      zavsi(:,:,:) = avs(:,:,:)
-# endif
-
+            ! No isopycnal diffusion
+            zwt(:,:,:) = fsavs(:,:,:)           
 #endif
-
-      ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ze3ta =  ( 1. - znvvl )   &                            ! after scale factor at T-point
-                  &   +        znvvl    * fse3t_a(ji,jj,jk) 
-               ze3tn =    znvvl          &                            ! now   scale factor at T-point
-                  &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)
-               zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
-               zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
-               zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! Surface boudary conditions
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3ta = ( 1. - znvvl )        &                           ! after scale factor at T-point
-               &   +       znvvl      * fse3t_a(ji,jj,1) 
-            zwi(ji,jj,1) = 0.e0
-            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
-         END DO
-      END DO
-
-
-      ! II.1. Vertical diffusion on t
-      ! ---------------------------
-
-      !! 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!
-
-      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            zwt(ji,jj,1) = zwd(ji,jj,1)
-         END DO
-      END DO
-      DO jk = 2, jpkm1
+            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
+                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
+                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
+                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
+                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
+                 END DO
+               END DO
+            END DO
+            ! Surface boudary conditions
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
+                 zwi(ji,jj,1) = 0.e0
+                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
+               END DO
+            END DO
+            !
+         END IF
+
+         ! II.1. Vertical diffusion on tracer
+         ! ---------------------------
+
+         !! 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!
+         
+         ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
          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
-
-      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
-            ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
-            ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1)
-         END DO
-      END DO
-      DO jk = 2, jpkm1
+               zwt(ji,jj,1) = zwd(ji,jj,1)
+            END DO
+         END DO
+         DO jk = 2, 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
+         
+         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
-               ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
-               zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk)   ! zrhs=right hand side 
-               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-
-      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
-      ! Save the masked temperature after in ta
-      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 1, -1
+               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
+               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
+               ptraa(ji,jj,1,jn) = ze3tb * ptrab(ji,jj,1,jn) + p2dt(1) * ze3tn * ptraa(ji,jj,1,jn)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
+                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
+                  zrhs = ze3tb * ptrab(ji,jj,jk,jn) + p2dt(jk) * ze3tn * ptraa(ji,jj,jk,jn)   ! zrhs=right hand side 
+                  ptraa(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ptraa(ji,jj,jk-1,jn)
+               END DO
+            END DO
+         END DO
+         
+         ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
+         ! Save the masked temperature after in ta
+         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               ta(ji,jj,jk) = ( ta(ji,jj,jk) - zws(ji,jj,jk) * ta(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! II.2 Vertical diffusion on salinity
-      ! -----------------------------------
-
-#if defined key_zdfddm
-      ! Rebuild the Matrix as avt /= avs
-
-      ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ze3ta =  ( 1. - znvvl )                        &         ! after scale factor at T-point
-                  &   +        znvvl   * fse3t_a(ji,jj,jk)           
-               ze3tn =    znvvl                               &         ! now   scale factor at T-point
-                  &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)
-               zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
-               zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
-               zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! Surface boudary conditions
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1)
-            zwi(ji,jj,1) = 0.e0
-            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
-         END DO
-      END DO
-#endif
-
-
-      !! 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!
-
-      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            zwt(ji,jj,1) = zwd(ji,jj,1)
-         END DO
-      END DO
-      DO jk = 2, 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
-
-      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            ze3tb = ( 1. - znvvl )   &                                           ! before scale factor at T-point
-               &   +  znvvl       * fse3t_b(ji,jj,1)
-            ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,1)                    ! now    scale factor at T-point
-            sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1)
-         END DO
-      END DO
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               ze3tb = ( 1. - znvvl )   &                                        ! before scale factor at T-point
-                  &   +  znvvl       * fse3t_b(ji,jj,jk)
-               ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)                ! now    scale factor at T-point
-               zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk)     ! zrhs=right hand side
-               sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-
-      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
-      ! Save the masked temperature after in ta
-      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 1, -1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               sa(ji,jj,jk) = ( sa(ji,jj,jk) - zws(ji,jj,jk) * sa(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
-            END DO
-         END DO
+               ptraa(ji,jj,jpkm1,jn) = ptraa(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
+            END DO
+         END DO
+         DO jk = jpk-2, 1, -1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  ptraa(ji,jj,jk,jn) = ( ptraa(ji,jj,jk,jn) - zws(ji,jj,jk) * ptraa(ji,jj,jk+1,jn) ) &
+                  &                    / zwt(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
       END DO
       !
