Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcsub.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 2904)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 2910)
@@ -26,5 +26,4 @@
    USE divcur          ! hor. divergence and curl      (div & cur routines)
    USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff 
-   USE trc, ONLY:  nn_dttrc, nittrc000 !TOP_TRC substepping
    USE sbc_oce         ! surface boundary condition: ocean
    USE bdy_oce
@@ -36,11 +35,15 @@
    IMPLICIT NONE
 
-   PUBLIC   trc_sub_stp    ! called by trc_stp
-   PUBLIC   trc_sub_ini    ! called by trc_ini to initialize substepping arrays.
-   PUBLIC   trc_sub_reset  ! called by trc_stp to reset physics variables
-   PUBLIC   trc_sub_ssh    ! called by trc_stp to reset physics variables
-
+   PUBLIC   trc_sub_stp      ! called by trc_stp
+   PUBLIC   trc_sub_ini      ! called by trc_ini to initialize substepping arrays.
+   PUBLIC   trc_sub_reset    ! called by trc_stp to reset physics variables
+   PUBLIC   trc_sub_ssh      ! called by trc_stp to reset physics variables
+
+   !!* Module variables
+   REAL(wp)  :: r1_ndttrc     !    1 /  nn_dttrc 
+   REAL(wp)  :: r1_ndttrcp1   !    1 / (nn_dttrc+1) 
+
+   !!* Substitution
 #  include "top_substitute.h90"
-   
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -59,6 +62,7 @@
       !!              on TOP steps, calculate averages.
       !!-------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::  kt        !ocean time-step index
-      INTEGER ::	        ji,jj,jk  !dummy loop indices
+      INTEGER, INTENT( in ) ::  kt        ! ocean time-step index
+      INTEGER               ::  ji,jj,jk  ! dummy loop indices
+      REAL(wp)              ::  z1_ne3t, z1_ne3u, z1_ne3v, z1_ne3w
       !!-------------------------------------------------------------------
 
@@ -67,506 +71,415 @@
            IF(lwp) WRITE(numout,*) 'trc_sub_stp : substepping of the passive tracers'
            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-
-           sshb_hold(:,:)=sshn(:,:)
-           sshu_b_hold(:,:)=sshu_n(:,:)
-           sshv_b_hold(:,:)=sshv_n(:,:)
-           emp_b_hold(:,:)=emp_b(:,:)
+           !
+           sshb_hold  (:,:) = sshn  (:,:)
+           sshu_b_hold(:,:) = sshu_n(:,:)
+           sshv_b_hold(:,:) = sshv_n(:,:)
+           emp_b_hold (:,:) = emp_b (:,:)
+           !
+           r1_ndttrc        = 1._wp / REAL( nn_dttrc    , wp ) 
+           r1_ndttrcp1      = 1._wp / REAL( nn_dttrc + 1, wp )
+           !
        ENDIF  
+
        IF( MOD( kt , nn_dttrc ) /= 0 ) THEN
+          !
+          un_tm   (:,:,:)        = un_tm   (:,:,:)        + un   (:,:,:)        * fse3u(:,:,:) 
+          vn_tm   (:,:,:)        = vn_tm   (:,:,:)        + vn   (:,:,:)        * fse3v(:,:,:) 
+          wn_tm   (:,:,:)        = wn_tm   (:,:,:)        + wn   (:,:,:)        * fse3w(:,:,:) 
+          tsn_tm  (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_tem) * fse3t(:,:,:)  
+          tsn_tm  (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_sal) * fse3t(:,:,:)  
+          avt_tm  (:,:,:)        = avt_tm  (:,:,:)        + avt  (:,:,:)        * fse3w(:,:,:)  
+# if defined key_zdfddm
+          avs_tm  (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * fse3w(:,:,:)  
+# endif
+#if defined key_ldfslp
+          wslpi_tm(:,:,:)        = wslpi_tm(:,:,:)        + wslpi(:,:,:)        * fse3w(:,:,:) 
+          wslpj_tm(:,:,:)        = wslpj_tm(:,:,:)        + wslpj(:,:,:)        * fse3w(:,:,:) 
+          uslp_tm (:,:,:)        = uslp_tm (:,:,:)        + uslp (:,:,:)        * fse3u(:,:,:) 
+          vslp_tm (:,:,:)        = vslp_tm (:,:,:)        + vslp (:,:,:)        * fse3v(:,:,:) 
+#endif
+          !
+          sshn_tm  (:,:)         = sshn_tm  (:,:)         + sshn  (:,:) 
+          sshu_n_tm(:,:)         = sshu_n_tm(:,:)         + sshu_n(:,:) 
+          sshv_n_tm(:,:)         = sshv_n_tm(:,:)         + sshv_n(:,:) 
+          rnf_tm   (:,:)         = rnf_tm   (:,:)         + rnf   (:,:) 
+          h_rnf_tm (:,:)         = h_rnf_tm (:,:)         + h_rnf (:,:) 
+          hmld_tm  (:,:)         = hmld_tm  (:,:)         + hmld  (:,:)
+          fr_i_tm  (:,:)         = fr_i_tm  (:,:)         + fr_i  (:,:)
+          emp_tm   (:,:)         = emp_tm   (:,:)         + emp   (:,:) 
+          emps_tm  (:,:)         = emps_tm  (:,:)         + emps  (:,:)
+          qsr_tm   (:,:)         = qsr_tm   (:,:)         + qsr   (:,:)
+          wndm_tm  (:,:)         = wndm_tm  (:,:)         + wndm  (:,:)
+          !
+#if defined key_traldf_c3d
+          ahtt_tm  (:,:,:)       = ahtt_tm  (:,:,:)       + ahtt(:,:,:)         * fse3t(:,:,:)
+          ahtu_tm  (:,:,:)       = ahtu_tm  (:,:,:)       + ahtu(:,:,:)         * fse3u(:,:,:)
+          ahtv_tm  (:,:,:)       = ahtv_tm  (:,:,:)       + ahtv(:,:,:)         * fse3v(:,:,:)
+          ahtw_tm  (:,:,:)       = ahtw_tm  (:,:,:)       + ahtw(:,:,:)         * fse3w(:,:,:)
+#elif defined key_traldf_c2d
+          ahtt_tm  (:,:)         = ahtt_tm  (:,:)         + ahtt(:,:)
+          ahtu_tm  (:,:)         = ahtu_tm  (:,:)         + ahtu(:,:)
+          ahtv_tm  (:,:)         = ahtv_tm  (:,:)         + ahtv(:,:)
+          ahtw_tm  (:,:)         = ahtw_tm  (:,:)         + ahtw(:,:)
+#elif defined key_traldf_c1d
+          ahtt_tm  (:)           = ahtt_tm  (:,:)         + ahtt(:)
+          ahtu_tm  (:)           = ahtu_tm  (:,:)         + ahtu(:)
+          ahtv_tm  (:)           = ahtv_tm  (:,:)         + ahtv(:)
+          ahtw_tm  (:)           = ahtw_tm  (:,:)         + ahtw(:)
+#else
+          ahtt_tm                = ahtt_tm                + ahtt
+          ahtu_tm                = ahtu_tm                + ahtu
+          ahtv_tm                = ahtv_tm                + ahtv
+          ahtw_tm                = ahtw_tm                + ahtw
+#endif
+#if defined key_traldf_eiv
+#  if defined key_traldf_c3d
+          aieu_tm  (:,:,:)       = aieu_tm  (:,:,:)       + aieu(:,:,:)         * fse3u(:,:,:)
+          aiev_tm  (:,:,:)       = aiev_tm  (:,:,:)       + aiev(:,:,:)         * fse3v(:,:,:)
+          aiew_tm  (:,:,:)       = aiew_tm  (:,:,:)       + aiew(:,:,:)         * fse3w(:,:,:)
+#  elif defined key_traldf_c2d
+          aieu_tm  (:,:)         = aieu_tm  (:,:)         + aieu(:,:)
+          aiev_tm  (:,:)         = aiev_tm  (:,:)         + aiev(:,:)
+          aiew_tm  (:,:)         = aiew_tm  (:,:)         + aiew(:,:)
+#  elif defined key_traldf_c1d
+          aieu_tm  (:)           = aieu_tm  (:,:)         + aieu(:)
+          aiev_tm  (:)           = aiev_tm  (:,:)         + aiev(:)
+          aiew_tm  (:)           = aiew_tm  (:,:)         + aiew(:)
+#  else
+          aieu_tm                = aieu_tm                + aieu
+          aiev_tm                = aiev_tm                + aiev
+          aiew_tm                = aiew_tm                + aiew
+#  endif
+#endif
+
+      ELSE                           !  It is time to substep 
+         !   1. set temporary arrays to hold physics variables
+         sshn_temp  (:,:)        = sshn  (:,:)
+         sshu_n_temp(:,:)        = sshu_n(:,:)
+         sshv_n_temp(:,:)        = sshv_n(:,:)
+         sshf_n_temp(:,:)        = sshf_n(:,:)
+         sshb_temp  (:,:)        = sshb  (:,:)
+         sshu_b_temp(:,:)        = sshu_b(:,:)
+         sshv_b_temp(:,:)        = sshv_b(:,:)
+         ssha_temp  (:,:)        = ssha  (:,:)
+         sshu_a_temp(:,:)        = sshu_a(:,:)
+         sshv_a_temp(:,:)        = sshv_a(:,:)
+         rnf_temp   (:,:)        = rnf   (:,:)
+         h_rnf_temp (:,:)        = h_rnf (:,:)
+         rotn_temp  (:,:,:)      = rotn  (:,:,:)
+         hdivn_temp (:,:,:)      = hdivn (:,:,:)
+         rotb_temp  (:,:,:)      = rotb  (:,:,:)
+         hdivb_temp (:,:,:)      = hdivb (:,:,:)
+         un_temp    (:,:,:)      = un    (:,:,:)
+         vn_temp    (:,:,:)      = vn    (:,:,:)
+         wn_temp    (:,:,:)      = wn    (:,:,:)
+         tsn_temp   (:,:,:,:)    = tsn   (:,:,:,:)
+         avt_temp   (:,:,:)      = avt   (:,:,:)
+# if defined key_zdfddm
+         avs_temp(:,:,:)         = avs   (:,:,:)
+# endif
+         hmld_temp  (:,:)        = hmld  (:,:)
+         qsr_temp   (:,:)        = qsr   (:,:)
+         wndm_temp  (:,:)        = wndm  (:,:)
+         emp_temp   (:,:)        = emp   (:,:)
+         emp_b_temp (:,:)        = emp_b (:,:)
+         emps_temp  (:,:)        = emps  (:,:)
+#if defined key_ldfslp
+         wslpi_temp (:,:,:)      = wslpi (:,:,:)
+         wslpj_temp (:,:,:)      = wslpj (:,:,:)
+         uslp_temp  (:,:,:)      = uslp  (:,:,:)
+         vslp_temp  (:,:,:)      = vslp  (:,:,:)
+#endif
+         hu_temp    (:,:)        = hu    (:,:)
+         hv_temp    (:,:)        = hv    (:,:)
+         hur_temp   (:,:)        = hur   (:,:)
+         hvr_temp   (:,:)        = hvr   (:,:)
+#if defined key_traldf_c3d
+         ahtu_temp  (:,:,:)      = ahtu  (:,:,:)
+         ahtv_temp  (:,:,:)      = ahtv  (:,:,:)
+         ahtw_temp  (:,:,:)      = ahtw  (:,:,:)
+         ahtt_temp  (:,:,:)      = ahtt  (:,:,:)
+#elif defined key_traldf_c2d
+         ahtu_temp  (:,:)        = ahtu  (:,:)
+         ahtv_temp  (:,:)        = ahtv  (:,:)
+         ahtw_temp  (:,:)        = ahtw  (:,:)
+         ahtt_temp  (:,:)        = ahtt  (:,:)
+#elif defined key_traldf_c1d
+         ahtu_temp  (:)          = ahtu  (:)
+         ahtv_temp  (:)          = ahtv  (:)
+         ahtw_temp  (:)          = ahtw  (:)
+         ahtt_temp  (:)          = ahtt  (:)
+#else
+         ahtu_temp               = ahtu
+         ahtv_temp               = ahtv
+         ahtw_temp               = ahtw
+         ahtt_temp               = ahtt
+#endif
+
+#if defined key_traldf_eiv
+# if defined key_traldf_c3d
+         aeiu_temp  (:,:,:)      = aeiu(:,:,:)
+         aeiv_temp  (:,:,:)      = aeiv(:,:,:)
+         aeiw_temp  (:,:,:)      = aeiw(:,:,:)
+# elif defined key_traldf_c2d
+         aeiu_temp  (:,:)        = aeiu(:,:)
+         aeiv_temp  (:,:)        = aeiv(:,:)
+         aeiw_temp  (:,:)        = aeiw(:,:)
+# elif defined key_traldf_c1d
+         aeiu_temp  (:)          = aeiu(:)
+         aeiv_temp  (:)          = aeiv(:)
+         aeiw_temp  (:)          = aeiw(:)
+# else
+         aeiu_temp               = aeiu
+         aeiv_temp               = aeiv
+         aeiw_temp               = aeiw
+# endif
+#endif
+         !
+         DO jk = 1, jpk
+            e3t_temp(:,:,jk)     = fse3t(:,:,jk)
+            e3u_temp(:,:,jk)     = fse3u(:,:,jk)
+            e3v_temp(:,:,jk)     = fse3v(:,:,jk)
+            e3w_temp(:,:,jk)     = fse3w(:,:,jk)
+         ENDDO
+         IF( lk_vvl ) THEN                      !  Update Now Vertical coord.  !   (only in vvl case)
+           !                                    !------------------------------!
+           DO jk = 1, jpk
+              fse3t (:,:,jk)     = fse3t_n (:,:,jk)   ! vertical scale factors stored in fse3. arrays
+              fse3u (:,:,jk)     = fse3u_n (:,:,jk)
+              fse3v (:,:,jk)     = fse3v_n (:,:,jk)
+              fse3w (:,:,jk)     = fse3w_n (:,:,jk)
+           END DO
+         ENDIF
+
+         ! 2. Create averages and reassign variables
+         sshn_tm  (:,:)          = sshn_tm    (:,:)       + sshn  (:,:) 
+         sshu_n_tm(:,:)          = sshu_n_tm  (:,:)       + sshu_n(:,:) 
+         sshv_n_tm(:,:)          = sshv_n_tm  (:,:)       + sshv_n(:,:) 
+         rnf_tm   (:,:)          = rnf_tm     (:,:)       + rnf   (:,:) 
+         h_rnf_tm (:,:)          = h_rnf_tm   (:,:)       + h_rnf (:,:) 
+         hmld_tm  (:,:)          = hmld_tm    (:,:)       + hmld  (:,:)
+         fr_i_tm  (:,:)          = fr_i_tm    (:,:)       + fr_i  (:,:)
+         emp_tm   (:,:)          = emp_tm     (:,:)       + emp   (:,:) 
+         emps_tm  (:,:)          = emps_tm    (:,:)       + emps  (:,:)
+         qsr_tm   (:,:)          = qsr_tm     (:,:)       + qsr   (:,:)
+         wndm_tm  (:,:)          = wndm_tm    (:,:)       + wndm  (:,:)
+         !
+         sshn     (:,:)          = sshn_tm    (:,:) * r1_ndttrcp1 
+         sshu_n   (:,:)          = sshu_n_tm  (:,:) * r1_ndttrcp1  
+         sshv_n   (:,:)          = sshv_n_tm  (:,:) * r1_ndttrcp1  
+         sshb     (:,:)          = sshb_hold  (:,:)
+         sshu_b   (:,:)          = sshu_b_hold(:,:)
+         sshv_b   (:,:)          = sshv_b_hold(:,:)
+         rnf      (:,:)          = rnf_tm     (:,:) * r1_ndttrcp1 
+         h_rnf    (:,:)          = h_rnf_tm   (:,:) * r1_ndttrcp1 
+         hmld     (:,:)          = hmld_tm    (:,:) * r1_ndttrcp1 
+         fr_i     (:,:)          = fr_i_tm    (:,:) * r1_ndttrcp1 
+         !  variables that are initialized after averages initialized
+         emp_b    (:,:) = emp_b_hold (:,:)
+         IF( kt == nittrc000 ) THEN
+            wndm  (:,:)          = wndm_tm    (:,:) * r1_ndttrc 
+            qsr   (:,:)          = qsr_tm     (:,:) * r1_ndttrc 
+            emp   (:,:)          = emp_tm     (:,:) * r1_ndttrc 
+            emps  (:,:)          = emps_tm    (:,:) * r1_ndttrc 
+         ELSE
+            wndm  (:,:)          = wndm_tm    (:,:) * r1_ndttrcp1 
+            qsr   (:,:)          = qsr_tm     (:,:) * r1_ndttrcp1 
+            emp   (:,:)          = emp_tm     (:,:) * r1_ndttrcp1 
+            emps  (:,:)          = emps_tm    (:,:) * r1_ndttrcp1 
+         ENDIF
+
+         un_tm    (:,:,:)        = un_tm   (:,:,:)        + un   (:,:,:)        * e3u_temp(:,:,:) 
+         vn_tm    (:,:,:)        = vn_tm   (:,:,:)        + vn   (:,:,:)        * e3v_temp(:,:,:) 
+         wn_tm    (:,:,:)        = wn_tm   (:,:,:)        + wn   (:,:,:)        * e3w_temp(:,:,:) 
+         tsn_tm   (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_tem) * e3t_temp(:,:,:)  
+         tsn_tm   (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_sal) * e3t_temp(:,:,:)  
+         avt_tm   (:,:,:)        = avt_tm  (:,:,:)        + avt  (:,:,:)        * e3w_temp(:,:,:)  
+# if defined key_zdfddm
+         avs_tm   (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * e3w_temp(:,:,:)  
+# endif
+#if defined key_ldfslp
+         wslpi_tm (:,:,:)        = wslpi_tm(:,:,:)        + wslpi(:,:,:)        * e3w_temp(:,:,:) 
+         wslpj_tm (:,:,:)        = wslpj_tm(:,:,:)        + wslpj(:,:,:)        * e3w_temp(:,:,:) 
+         uslp_tm  (:,:,:)        = uslp_tm (:,:,:)        + uslp (:,:,:)        * e3u_temp(:,:,:) 
+         vslp_tm  (:,:,:)        = vslp_tm (:,:,:)        + vslp (:,:,:)        * e3v_temp(:,:,:) 
+#endif
+         !
          DO jk = 1, jpk
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  un_tm(ji,jj,jk)= un_tm(ji,jj,jk)+un(ji,jj,jk)*fse3u(ji,jj,jk)
-                  vn_tm(ji,jj,jk)= vn_tm(ji,jj,jk)+vn(ji,jj,jk)*fse3v(ji,jj,jk)
-		  wn_tm(ji,jj,jk)= wn_tm(ji,jj,jk)+wn(ji,jj,jk)*fse3w(ji,jj,jk)
-                  avt_tm(ji,jj,jk)= avt_tm(ji,jj,jk)+avt(ji,jj,jk)*fse3t(ji,jj,jk)
-                  tsn_tm(ji,jj,jk,jp_sal)= tsn_tm(ji,jj,jk,jp_sal)+tsn(ji,jj,jk,jp_sal)*fse3t(ji,jj,jk)
-                  tsn_tm(ji,jj,jk,jp_tem)= tsn_tm(ji,jj,jk,jp_tem)+tsn(ji,jj,jk,jp_tem)*fse3t(ji,jj,jk)
+                  z1_ne3t = r1_ndttrcp1  / fse3t(ji,jj,jk)
+                  z1_ne3u = r1_ndttrcp1  / fse3u(ji,jj,jk)
+                  z1_ne3v = r1_ndttrcp1  / fse3v(ji,jj,jk)
+                  z1_ne3w = r1_ndttrcp1  / fse3w(ji,jj,jk)
+                  !
+                  un   (ji,jj,jk)        = un_tm   (ji,jj,jk)        * z1_ne3u
+                  vn   (ji,jj,jk)        = vn_tm   (ji,jj,jk)        * z1_ne3v
+                  wn   (ji,jj,jk)        = wn_tm   (ji,jj,jk)        * z1_ne3w
+                  tsn  (ji,jj,jk,jp_tem) = tsn_tm  (ji,jj,jk,jp_tem) * z1_ne3t
+                  tsn  (ji,jj,jk,jp_sal) = tsn_tm  (ji,jj,jk,jp_sal) * z1_ne3t
+                  avt  (ji,jj,jk)        = avt_tm  (ji,jj,jk)        * z1_ne3w
 # if defined key_zdfddm
-                  avs_tm(ji,jj,jk)= avs_tm(ji,jj,jk)+avs(ji,jj,jk)*fse3w(ji,jj,jk)
+                  avs  (ji,jj,jk)        = avs_tm  (ji,jj,jk)        * z1_ne3w
 # endif
 #if defined key_ldfslp
-		  wslpi_tm(ji,jj,jk)= wslpi_tm(ji,jj,jk)+wslpi(ji,jj,jk)*fse3w(ji,jj,jk)
-		  wslpj_tm(ji,jj,jk)= wslpj_tm(ji,jj,jk)+wslpj(ji,jj,jk)*fse3w(ji,jj,jk)
-		  uslp_tm(ji,jj,jk)= uslp_tm(ji,jj,jk)+uslp(ji,jj,jk)*fse3u(ji,jj,jk)
-		  vslp_tm(ji,jj,jk)= vslp_tm(ji,jj,jk)+vslp(ji,jj,jk)*fse3v(ji,jj,jk)
-#endif
-               END DO
-            END DO
-         END DO
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-         	sshn_tm(ji,jj)=sshn_tm(ji,jj)+sshn(ji,jj)
-                sshu_n_tm(ji,jj)=sshu_n_tm(ji,jj)+sshu_n(ji,jj)
-                sshv_n_tm(ji,jj)=sshv_n_tm(ji,jj)+sshv_n(ji,jj)
-	        rnf_tm(ji,jj)= rnf_tm(ji,jj)+rnf(ji,jj)
-	        h_rnf_tm(ji,jj)= h_rnf_tm(ji,jj)+h_rnf(ji,jj)
-         	hmld_tm(ji,jj)=hmld_tm(ji,jj)+hmld(ji,jj)
-         	qsr_tm(ji,jj)=qsr_tm(ji,jj)+qsr(ji,jj)
-         	fr_i_tm(ji,jj)=fr_i_tm(ji,jj)+fr_i(ji,jj)
-         	wndm_tm(ji,jj)=wndm_tm(ji,jj)+wndm(ji,jj)
-         	emp_tm(ji,jj)=emp_tm(ji,jj)+emp(ji,jj)
-         	emps_tm(ji,jj)=emps_tm(ji,jj)+emps(ji,jj)
-             END DO
-         END DO
+                  wslpi(ji,jj,jk)        = wslpi_tm(ji,jj,jk)        * z1_ne3w 
+                  wslpj(ji,jj,jk)        = wslpj_tm(ji,jj,jk)        * z1_ne3w 
+                  uslp (ji,jj,jk)        = uslp_tm (ji,jj,jk)        * z1_ne3u 
+                  vslp (ji,jj,jk)        = vslp_tm (ji,jj,jk)        * z1_ne3v 
+#endif
+               ENDDO
+            ENDDO
+         ENDDO
+
 #if defined key_traldf_c3d
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ahtt_tm(ji,jj,jk)= ahtt_tm(ji,jj,jk)+ahtt(ji,jj,jk)*fse3t(ji,jj,jk)
-                  ahtu_tm(ji,jj,jk)= ahtu_tm(ji,jj,jk)+ahtu(ji,jj,jk)*fse3u(ji,jj,jk)
-                  ahtv_tm(ji,jj,jk)= ahtv_tm(ji,jj,jk)+ahtv(ji,jj,jk)*fse3v(ji,jj,jk)
-                  ahtw_tm(ji,jj,jk)= ahtw_tm(ji,jj,jk)+ahtw(ji,jj,jk)*fse3w(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+          ahtt_tm  (:,:,:)       = ahtt_tm  (:,:,:)       + ahtt(:,:,:)         * e3t_temp(:,:,:)
+          ahtu_tm  (:,:,:)       = ahtu_tm  (:,:,:)       + ahtu(:,:,:)         * e3u_temp(:,:,:)
+          ahtv_tm  (:,:,:)       = ahtv_tm  (:,:,:)       + ahtv(:,:,:)         * e3v_temp(:,:,:)
+          ahtw_tm  (:,:,:)       = ahtw_tm  (:,:,:)       + ahtw(:,:,:)         * e3w_temp(:,:,:)
+          !
+          ahtt     (:,:,:)       = ahtt_tm  (:,:,:) * r1_ndttrcp1  / fse3t(:,:,:)
+          ahtu     (:,:,:)       = ahtu_tm  (:,:,:) * r1_ndttrcp1  / fse3u(:,:,:)
+          ahtv     (:,:,:)       = ahtv_tm  (:,:,:) * r1_ndttrcp1  / fse3v(:,:,:)
+          ahtw     (:,:,:)       = ahtw_tm  (:,:,:) * r1_ndttrcp1  / fse3w(:,:,:)
 #elif defined key_traldf_c2d
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ahtt_tm(ji,jj)= ahtt_tm(ji,jj)+ahtt(ji,jj)
-               ahtu_tm(ji,jj)= ahtu_tm(ji,jj)+ahtu(ji,jj)
-               ahtv_tm(ji,jj)= ahtv_tm(ji,jj)+ahtv(ji,jj)
-               ahtw_tm(ji,jj)= ahtw_tm(ji,jj)+ahtw(ji,jj)
-            END DO
-         END DO
+          ahtt_tm  (:,:)         = ahtt_tm  (:,:)         + ahtt(:,:)
+          ahtu_tm  (:,:)         = ahtu_tm  (:,:)         + ahtu(:,:)
+          ahtv_tm  (:,:)         = ahtv_tm  (:,:)         + ahtv(:,:)
+          ahtw_tm  (:,:)         = ahtw_tm  (:,:)         + ahtw(:,:)
+          !
+          ahtt     (:,:)         = ahtt_tm  (:,:)   * r1_ndttrcp1
+          ahtu     (:,:)         = ahtu_tm  (:,:)   * r1_ndttrcp1
+          ahtv     (:,:)         = ahtv_tm  (:,:)   * r1_ndttrcp1
+          ahtw     (:,:)         = ahtw_tm  (:,:)   * r1_ndttrcp1
 #elif defined key_traldf_c1d
-         DO jk = 1, jpk
-            ahtt_tm(jk)= ahtt_tm(jk)+ahtt(jk)
-            ahtu_tm(jk)= ahtu_tm(jk)+ahtu(jk)
-            ahtv_tm(jk)= ahtv_tm(jk)+ahtv(jk)
-            ahtw_tm(jk)= ahtw_tm(jk)+ahtw(jk)
-         END DO
+          ahtt_tm  (:)           = ahtt_tm  (:,:)         + ahtt(:)
+          ahtu_tm  (:)           = ahtu_tm  (:,:)         + ahtu(:)
+          ahtv_tm  (:)           = ahtv_tm  (:,:)         + ahtv(:)
+          ahtw_tm  (:)           = ahtw_tm  (:,:)         + ahtw(:)
+          !
+          ahtt     (:)           = ahtt_tm  (:)     * r1_ndttrcp1
+          ahtu     (:)           = ahtu_tm  (:)     * r1_ndttrcp1
+          ahtv     (:)           = ahtv_tm  (:)     * r1_ndttrcp1
+          ahtw     (:)           = ahtw_tm  (:)     * r1_ndttrcp1
 #else
-         ahtt_tm= ahtt_tm+ahtt
-         ahtu_tm= ahtu_tm+ahtu
-         ahtv_tm= ahtv_tm+ahtv
-         ahtw_tm= ahtw_tm+ahtw
-#endif
+          ahtt_tm                = ahtt_tm                + ahtt
+          ahtu_tm                = ahtu_tm                + ahtu
+          ahtv_tm                = ahtv_tm                + ahtv
+          ahtw_tm                = ahtw_tm                + ahtw
+          !
+          ahtt                   = ahtt_tm          * r1_ndttrcp1
+          ahtu                   = ahtu_tm          * r1_ndttrcp1
+          ahtv                   = ahtv_tm          * r1_ndttrcp1
+          ahtw                   = ahtw_tm          * r1_ndttrcp1
+#endif
+
 #if defined key_traldf_eiv
 # if defined key_traldf_c3d
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  aeiu_tm(ji,jj,jk)= aeiu_tm(ji,jj,jk)+aeiu(ji,jj,jk)*fse3u(ji,jj,jk)
-                  aeiv_tm(ji,jj,jk)= aeiv_tm(ji,jj,jk)+aeiv(ji,jj,jk)*fse3v(ji,jj,jk)
-                  aeiw_tm(ji,jj,jk)= aeiw_tm(ji,jj,jk)+aeiw(ji,jj,jk)*fse3w(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+          aeiu_tm  (:,:,:)       = aeiu_tm  (:,:,:)       + aeiu(:,:,:)         * e3u_temp(:,:,:)
+          aeiv_tm  (:,:,:)       = aeiv_tm  (:,:,:)       + aeiv(:,:,:)         * e3v_temp(:,:,:)
+          aeiw_tm  (:,:,:)       = aeiw_tm  (:,:,:)       + aeiw(:,:,:)         * e3w_temp(:,:,:)
+          !
+          aeiu     (:,:,:)       = aeiu_tm  (:,:,:) * r1_ndttrcp1  / fse3u(:,:,:)
+          aeiv     (:,:,:)       = aeiv_tm  (:,:,:) * r1_ndttrcp1  / fse3v(:,:,:)
+          aeiw     (:,:,:)       = aeiw_tm  (:,:,:) * r1_ndttrcp1  / fse3w(:,:,:)
 # elif defined key_traldf_c2d
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               aeiu_tm(ji,jj)= aeiu_tm(ji,jj)+aeiu(ji,jj)
-               aeiv_tm(ji,jj)= aeiv_tm(ji,jj)+aeiv(ji,jj)
-               aeiw_tm(ji,jj)= aeiw_tm(ji,jj)+aeiw(ji,jj)
-            END DO
-         END DO
+          aeiu_tm  (:,:)         = aeiu_tm  (:,:)         + aeiu(:,:)
+          aeiv_tm  (:,:)         = aeiv_tm  (:,:)         + aeiv(:,:)
+          aeiw_tm  (:,:)         = aeiw_tm  (:,:)         + aeiw(:,:)
+          !
+          aeiu     (:,:)         = aeiu_tm  (:,:)   * r1_ndttrcp1
+          aeiv     (:,:)         = aeiv_tm  (:,:)   * r1_ndttrcp1
+          aeiw     (:,:)         = aeiw_tm  (:,:)   * r1_ndttrcp1
 # elif defined key_traldf_c1d
-         DO jk = 1, jpk
-            aeiu_tm(jk)= aeiu_tm(jk)+aeiu(jk)
-            aeiv_tm(jk)= aeiv_tm(jk)+aeiv(jk)
-            aeiw_tm(jk)= aeiw_tm(jk)+aeiw(jk)
-         END DO
+          aeiu_tm  (:)           = aeiu_tm  (:,:)         + aeiu(:)
+          aeiv_tm  (:)           = aeiv_tm  (:,:)         + aeiv(:)
+          aeiw_tm  (:)           = aeiw_tm  (:,:)         + aeiw(:)
+          !
+          aeiu     (:)           = aeiu_tm  (:)     * r1_ndttrcp1
+          aeiv     (:)           = aeiv_tm  (:)     * r1_ndttrcp1
+          aeiw     (:)           = aeiw_tm  (:)     * r1_ndttrcp1
 # else
-         aeiu_tm= aeiu_tm+aeiu
-         aeiv_tm= aeiv_tm+aeiv
-         aeiw_tm= aeiw_tm+aeiw
-# endif
-# endif
-
-      ELSE
-!   it is time to substep
-!   1. set temporary arrays to hold physics variables
-         sshn_temp(:,:)=sshn(:,:)
-         sshu_n_temp(:,:)=sshu_n(:,:)
-         sshv_n_temp(:,:)=sshv_n(:,:)
-         sshf_n_temp(:,:)=sshf_n(:,:)
-         sshb_temp(:,:)=sshb(:,:)
-         sshu_b_temp(:,:)=sshu_b(:,:)
-         sshv_b_temp(:,:)=sshv_b(:,:)
-         ssha_temp(:,:)=ssha(:,:)
-         sshu_a_temp(:,:)=sshu_a(:,:)
-         sshv_a_temp(:,:)=sshv_a(:,:)
-         rnf_temp(:,:)=rnf(:,:)
-         h_rnf_temp(:,:)=h_rnf(:,:)
-         rotn_temp(:,:,:)=rotn(:,:,:)
-         hdivn_temp(:,:,:)=hdivn(:,:,:)
-         rotb_temp(:,:,:)=rotb(:,:,:)
-         hdivb_temp(:,:,:)=hdivb(:,:,:)
-         un_temp(:,:,:)=un(:,:,:)
-         vn_temp(:,:,:)=vn(:,:,:)
-         wn_temp(:,:,:)=wn(:,:,:)
-         tsn_temp(:,:,:,:)=tsn(:,:,:,:)
-         avt_temp(:,:,:)=avt(:,:,:)
-         hmld_temp(:,:)=hmld(:,:)
-         qsr_temp(:,:)=qsr(:,:)
-         wndm_temp(:,:)=wndm(:,:)
-         emp_temp(:,:)=emp(:,:)
-         emp_b_temp(:,:)=emp_b(:,:)
-         emps_temp(:,:)=emps(:,:)
+          aeiu_tm                = aeiu_tm                + aeiu
+          aeiv_tm                = aeiv_tm                + aeiv
+          aeiw_tm                = aeiw_tm                + aeiw
+          !
+          aeiu                   = aeiu_tm          * r1_ndttrcp1
+          aeiv                   = aeiv_tm          * r1_ndttrcp1
+          aeiw                   = aeiw_tm          * r1_ndttrcp1
+# endif
+#endif
+
+         CALL lbc_lnk( sshn  (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( sshu_n(:,:)         , 'U', 1. ) 
+         CALL lbc_lnk( sshv_n(:,:)         , 'V', 1. ) 
+         CALL lbc_lnk( sshf_n(:,:)         , 'F', 1. ) 
+         CALL lbc_lnk( sshb  (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( sshu_b(:,:)         , 'U', 1. ) 
+         CALL lbc_lnk( sshv_b(:,:)         , 'V', 1. ) 
+         CALL lbc_lnk( ssha  (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( sshu_a(:,:)         , 'U', 1. ) 
+         CALL lbc_lnk( sshv_a(:,:)         , 'V', 1. ) 
+         CALL lbc_lnk( rnf   (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( h_rnf (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( rotn  (:,:,:)       , 'F', 1. ) 
+         CALL lbc_lnk( hdivn (:,:,:)       , 'T', 1. ) 
+         CALL lbc_lnk( rotb  (:,:,:)       , 'F', 1. ) 
+         CALL lbc_lnk( hdivb (:,:,:)       , 'T', 1. ) 
+         CALL lbc_lnk( un    (:,:,:)       , 'U', 1. ) 
+         CALL lbc_lnk( vn    (:,:,:)       , 'V', 1. ) 
+         CALL lbc_lnk( wn    (:,:,:)       , 'V', 1. ) 
+         CALL lbc_lnk( tsn   (:,:,:,jp_tem), 'T', 1. ) 
+         CALL lbc_lnk( tsn   (:,:,:,jp_sal), 'T', 1. ) 
+         CALL lbc_lnk( avt   (:,:,:)       , 'W', 1. ) 
+         CALL lbc_lnk( hmld  (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( qsr   (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( wndm  (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( emp   (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( emp_b (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( emps  (:,:)         , 'T', 1. ) 
 #if defined key_traldf_c3d
-         ahtu_temp(:,:,:)=ahtu(:,:,:)
-         ahtv_temp(:,:,:)=ahtv(:,:,:)
-         ahtw_temp(:,:,:)=ahtw(:,:,:)
-         ahtt_temp(:,:,:)=ahtt(:,:,:)
+         CALL lbc_lnk( ahtt  (:,:,:)       , 'T', 1. ) 
+         CALL lbc_lnk( ahtu  (:,:,:)       , 'U', 1. ) 
+         CALL lbc_lnk( ahtv  (:,:,:)       , 'V', 1. ) 
+         CALL lbc_lnk( ahtw  (:,:,:)       , 'W', 1. ) 
 #elif defined key_traldf_c2d
-         ahtu_temp(:,:)=ahtu(:,:)
-         ahtv_temp(:,:)=ahtv(:,:)
-         ahtw_temp(:,:)=ahtw(:,:)
-         ahtt_temp(:,:)=ahtt(:,:)
-#elif defined key_traldf_c1d
-         ahtu_temp(:)=ahtu(:)
-         ahtv_temp(:)=ahtv(:)
-         ahtw_temp(:)=ahtw(:)
-         ahtt_temp(:)=ahtt(:)
-#else
-         ahtu_temp=ahtu
-         ahtv_temp=ahtv
-         ahtw_temp=ahtw
-         ahtt_temp=ahtt
+         CALL lbc_lnk( ahtt  (:,:)         , 'T', 1. ) 
+         CALL lbc_lnk( ahtu  (:,:)         , 'U', 1. ) 
+         CALL lbc_lnk( ahtv  (:,:)         , 'V', 1. ) 
+         CALL lbc_lnk( ahtw  (:,:)         , 'W', 1. ) 
 #endif
 # if defined key_zdfddm
-         avs_temp(:,:,:)=avs(:,:,:)
+          CALL lbc_lnk( avs  (:,:,:)       , 'W', 1. ) 
 # endif
 #if defined key_traldf_eiv
 #if defined key_traldf_c3d
-         aeiu_temp(:,:,:)=aeiu(:,:,:)
-         aeiv_temp(:,:,:)=aeiv(:,:,:)
-         aeiw_temp(:,:,:)=aeiw(:,:,:)
+         CALL lbc_lnk( aeiu  (:,:,:)       , 'U', 1. ) 
+         CALL lbc_lnk( aeiv  (:,:,:)       , 'V', 1. ) 
+         CALL lbc_lnk( aeiw  (:,:,:)       , 'W', 1. ) 
 #elif defined key_traldf_c2d
-         aeiu_temp(:,:)=aeiu(:,:)
-         aeiv_temp(:,:)=aeiv(:,:)
-         aeiw_temp(:,:)=aeiw(:,:)
-#elif defined key_traldf_c1d
-         aeiu_temp(:)=aeiu(:)
-         aeiv_temp(:)=aeiv(:)
-         aeiw_temp(:)=aeiw(:)
-#else
-         aeiu_temp=aeiu
-         aeiv_temp=aeiv
-         aeiw_temp=aeiw
+         CALL lbc_lnk( aeiu  (:,:)         , 'U', 1. ) 
+         CALL lbc_lnk( aeiv  (:,:)         , 'V', 1. ) 
+         CALL lbc_lnk( aeiw  (:,:)         , 'W', 1. ) 
 #endif
 #endif
 #if defined key_ldfslp
-         wslpi_temp(:,:,:)=wslpi(:,:,:)
-         wslpj_temp(:,:,:)=wslpj(:,:,:)
-         uslp_temp(:,:,:)=uslp(:,:,:)
-         vslp_temp(:,:,:)=vslp(:,:,:)
-#endif
-         hu_temp(:,:)=hu(:,:)
-         hv_temp(:,:)=hv(:,:)
-         hur_temp(:,:)=hur(:,:)
-         hvr_temp(:,:)=hvr(:,:)
-
-!   2. create averages and reassign variables
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               sshn_tm(ji,jj)= sshn_tm(ji,jj)+sshn(ji,jj)
-               sshn(ji,jj)=sshn_tm(ji,jj)/(nn_dttrc+1)
-		  
-               sshu_n_tm(ji,jj)= sshu_n_tm(ji,jj)+sshu_n(ji,jj)
-               sshu_n(ji,jj)=sshu_n_tm(ji,jj)/(nn_dttrc+1)
-
-               sshv_n_tm(ji,jj)= sshv_n_tm(ji,jj)+sshv_n(ji,jj)
-               sshv_n(ji,jj)=sshv_n_tm(ji,jj)/(nn_dttrc+1)
-
-               sshb(ji,jj)=sshb_hold(ji,jj)
-
-               sshu_b(ji,jj)=sshu_b_hold(ji,jj)
-
-	       sshv_b(ji,jj)=sshv_b_hold(ji,jj)
-
-               rnf_tm(ji,jj)= rnf_tm(ji,jj)+rnf(ji,jj)
-               rnf(ji,jj)=rnf_tm(ji,jj)/(nn_dttrc+1)
- 
-               h_rnf_tm(ji,jj)= h_rnf_tm(ji,jj)+h_rnf(ji,jj)
-               h_rnf(ji,jj)=h_rnf_tm(ji,jj)/(nn_dttrc+1)
-
-               hmld_tm(ji,jj)= hmld_tm(ji,jj)+hmld(ji,jj)
-               hmld(ji,jj)=hmld_tm(ji,jj)/(nn_dttrc+1)
- 
-               fr_i_tm(ji,jj)= fr_i_tm(ji,jj)+fr_i(ji,jj)
-               fr_i(ji,jj)=fr_i_tm(ji,jj)/(nn_dttrc+1)
-
-!  variables that are initialized after averages initialized
-               wndm_tm(ji,jj)= wndm_tm(ji,jj)+wndm(ji,jj)
-               IF(kt == nittrc000)THEN
-                   wndm(ji,jj)=wndm_tm(ji,jj)/(nn_dttrc)
-               ELSE
-                   wndm(ji,jj)=wndm_tm(ji,jj)/(nn_dttrc+1)
-               ENDIF 
-
-               qsr_tm(ji,jj)= qsr_tm(ji,jj)+qsr(ji,jj)
-               IF(kt == nittrc000)THEN
-                   qsr(ji,jj)=qsr_tm(ji,jj)/(nn_dttrc)
-               ELSE
-                   qsr(ji,jj)=qsr_tm(ji,jj)/(nn_dttrc+1)
-               ENDIF 
-
-	       emp_b(ji,jj)=emp_b_hold(ji,jj)
-
-               emp_tm(ji,jj)= emp_tm(ji,jj)+emp(ji,jj)
-               IF(kt == nittrc000)THEN
-                  emp(ji,jj)=emp_tm(ji,jj)/(nn_dttrc)
-               ELSE
-                  emp(ji,jj)=emp_tm(ji,jj)/(nn_dttrc+1)
-               ENDIF 
-
-               emps_tm(ji,jj)= emps_tm(ji,jj)+emps(ji,jj)
-               IF(kt == nittrc000)THEN
-                  emps(ji,jj)=emps_tm(ji,jj)/(nn_dttrc)
-               ELSE
-                  emps(ji,jj)=emps_tm(ji,jj)/(nn_dttrc+1)
-               ENDIF
-
-            END DO
-         END DO
-         DO jk = 1, jpk
-            e3t_temp(:,:,jk)= fse3t(:,:,jk) 
-            e3u_temp(:,:,jk)= fse3u(:,:,jk) 
-            e3v_temp(:,:,jk)= fse3v(:,:,jk)
-            e3w_temp(:,:,jk)= fse3w(:,:,jk)  
-         END DO
-         IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
-         !                                           !------------------------------!
-            DO jk = 1, jpk
-                 fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
-                 fse3u (:,:,jk) = fse3u_n (:,:,jk)
-                 fse3v (:,:,jk) = fse3v_n (:,:,jk)
-                 fse3w (:,:,jk) = fse3w_n (:,:,jk)
-            END DO
-         ENDIF
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-
-                  un_tm(ji,jj,jk)= un_tm(ji,jj,jk)+un(ji,jj,jk)*e3u_temp(ji,jj,jk)
-                  un(ji,jj,jk)=(un_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)
-		  
-                  vn_tm(ji,jj,jk)= vn_tm(ji,jj,jk)+vn(ji,jj,jk)*e3v_temp(ji,jj,jk)
-                  vn(ji,jj,jk)=(vn_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)
-		  
-                  wn_tm(ji,jj,jk)= wn_tm(ji,jj,jk)+wn(ji,jj,jk)*e3w_temp(ji,jj,jk)
-                  wn(ji,jj,jk)=(wn_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
-		  
-                  avt_tm(ji,jj,jk)= avt_tm(ji,jj,jk)+avt(ji,jj,jk)*e3t_temp(ji,jj,jk)
-                  avt(ji,jj,jk)=(avt_tm(ji,jj,jk)/(nn_dttrc+1))/fse3t(ji,jj,jk)	
-	  	  
-                  tsn_tm(ji,jj,jk,jp_tem)= tsn_tm(ji,jj,jk,jp_tem)+tsn(ji,jj,jk,jp_tem)*e3t_temp(ji,jj,jk)
-                  tsn(ji,jj,jk,jp_tem)=(tsn_tm(ji,jj,jk,jp_tem)/(nn_dttrc+1))/fse3t(ji,jj,jk)
-	  
-                  tsn_tm(ji,jj,jk,jp_sal)= tsn_tm(ji,jj,jk,jp_sal)+tsn(ji,jj,jk,jp_sal)*e3t_temp(ji,jj,jk)
-                  tsn(ji,jj,jk,jp_sal)=(tsn_tm(ji,jj,jk,jp_sal)/(nn_dttrc+1))/fse3t(ji,jj,jk)
-	
-# if defined key_zdfddm
-                  avs_tm(ji,jj,jk)= avs_tm(ji,jj,jk)+avs(ji,jj,jk)*e3w_temp(ji,jj,jk)
-                  avs(ji,jj,jk)=(avs_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
-
-# endif
-	  
-#if defined key_ldfslp
-                  wslpi_tm(ji,jj,jk)= wslpi_tm(ji,jj,jk)+wslpi(ji,jj,jk)*e3w_temp(ji,jj,jk)
-                  wslpi(ji,jj,jk)=(wslpi_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
-
-                  wslpj_tm(ji,jj,jk)= wslpj_tm(ji,jj,jk)+wslpj(ji,jj,jk)*e3w_temp(ji,jj,jk)
-                  wslpj(ji,jj,jk)=(wslpj_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
-
-                  uslp_tm(ji,jj,jk)= uslp_tm(ji,jj,jk)+uslp(ji,jj,jk)*e3u_temp(ji,jj,jk)
-                  uslp(ji,jj,jk)=(uslp_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)
-
-                  vslp_tm(ji,jj,jk)= vslp_tm(ji,jj,jk)+vslp(ji,jj,jk)*e3v_temp(ji,jj,jk)
-                  vslp(ji,jj,jk)=(vslp_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)
-#endif
-	       END DO
-            END DO
-         END DO
-#if defined key_traldf_c3d
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ahtt_tm(ji,jj,jk)= ahtt_tm(ji,jj,jk)+ahtt(ji,jj,jk)*e3t_temp(ji,jj,jk)
-                  ahtt(ji,jj,jk)=(ahtt_tm(ji,jj,jk)/(nn_dttrc+1))/fse3t(ji,jj,jk)	
-
-                  ahtu_tm(ji,jj,jk)= ahtu_tm(ji,jj,jk)+ahtu(ji,jj,jk)*e3u_temp(ji,jj,jk)
-                  ahtu(ji,jj,jk)=(ahtu_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)	
-
-                  ahtv_tm(ji,jj,jk)= ahtv_tm(ji,jj,jk)+ahtv(ji,jj,jk)*e3v_temp(ji,jj,jk)
-                  ahtv(ji,jj,jk)=(ahtv_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)	
-
-                  ahtw_tm(ji,jj,jk)= ahtw_tm(ji,jj,jk)+ahtw(ji,jj,jk)*e3w_temp(ji,jj,jk)
-                  ahtw(ji,jj,jk)=(ahtw_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)	
-
-               END DO
-            END DO
-         END DO
-#elif defined key_traldf_c2d
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ahtt_tm(ji,jj)= ahtt_tm(ji,jj)+ahtt(ji,jj)
-               ahtt(ji,jj)=ahtt_tm(ji,jj)/(nn_dttrc+1)
-
-               ahtu_tm(ji,jj)= ahtu_tm(ji,jj)+ahtu(ji,jj)
-               ahtu(ji,jj)=ahtu_tm(ji,jj)/(nn_dttrc+1)
-
-               ahtv_tm(ji,jj)= ahtv_tm(ji,jj)+ahtv(ji,jj)
-               ahtv(ji,jj)=ahtv_tm(ji,jj)/(nn_dttrc+1)
-
-               ahtw_tm(ji,jj)= ahtw_tm(ji,jj)+ahtw(ji,jj)
-               ahtw(ji,jj)=ahtw_tm(ji,jj)/(nn_dttrc+1)
-            END DO
-         END DO
-#elif defined key_traldf_c1d
-         DO jk = 1, jpk
-            ahtt_tm(jk)= ahtt_tm(jk)+ahtt(jk)
-            ahtt(jk)=ahtt_tm(jk)/(nn_dttrc+1)
-
-            ahtu_tm(jk)= ahtu_tm(jk)+ahtu(jk)
-            ahtu(jk)=ahtu_tm(jk)/(nn_dttrc+1)
-
-            ahtv_tm(jk)= ahtv_tm(jk)+ahtv(jk)
-            ahtv(jk)=ahtv_tm(jk)/(nn_dttrc+1)
-
-            ahtw_tm(jk)= ahtw_tm(jk)+ahtw(jk)
-            ahtw(jk)=ahtw_tm(jk)/(nn_dttrc+1)
-         END DO
-#else
-         ahtt_tm= ahtt_tm+ahtt
-         ahtt=ahtt_tm/(nn_dttrc+1)
-
-         ahtu_tm= ahtu_tm+ahtu
-         ahtu=ahtu_tm/(nn_dttrc+1)
-
-         ahtv_tm= ahtv_tm+ahtv
-         ahtv=ahtv_tm/(nn_dttrc+1)
-
-         ahtw_tm= ahtw_tm+ahtw
-         ahtw=ahtw_tm/(nn_dttrc+1)
-
-#endif
-#if defined key_traldf_eiv
-# if defined key_traldf_c3d
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  aeiu_tm(ji,jj,jk)= aeiu_tm(ji,jj,jk)+aeiu(ji,jj,jk)*e3u_temp(ji,jj,jk)
-                  aeiu(ji,jj,jk)=(aeiu_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)	
-
-                  aeiv_tm(ji,jj,jk)= aeiv_tm(ji,jj,jk)+aeiv(ji,jj,jk)*e3v_temp(ji,jj,jk)
-                  aeiv(ji,jj,jk)=(aeiv_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)	
-
-                  aeiw_tm(ji,jj,jk)= aeiw_tm(ji,jj,jk)+aeiw(ji,jj,jk)*e3w_temp(ji,jj,jk)
-                  aeiw(ji,jj,jk)=(aeiw_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)	
-
-               END DO
-            END DO
-         END DO
-# elif defined key_traldf_c2d
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               aeiu_tm(ji,jj)= aeiu_tm(ji,jj)+aeiu(ji,jj)
-               aeiu(ji,jj)=aeiu_tm(ji,jj)/(nn_dttrc+1)
-
-               aeiv_tm(ji,jj)= aeiv_tm(ji,jj)+aeiv(ji,jj)
-               aeiv(ji,jj)=aeiv_tm(ji,jj)/(nn_dttrc+1)
-
-               aeiw_tm(ji,jj)= aeiw_tm(ji,jj)+aeiw(ji,jj)
-               aeiw(ji,jj)=aeiw_tm(ji,jj)/(nn_dttrc+1)
-            END DO
-         END DO
-# elif defined key_traldf_c1d
-         DO jk = 1, jpk
-            aeiu_tm(jk)= aeiu_tm(jk)+aeiu(jk)
-            aeiu(jk)=aeiu_tm(jk)/(nn_dttrc+1)
-
-            aeiv_tm(jk)= aeiv_tm(jk)+aeiv(jk)
-            aeiv(jk)=aeiv_tm(jk)/(nn_dttrc+1)
-
-            aeiw_tm(jk)= aeiw_tm(jk)+aeiw(jk)
-            aeiw(jk)=aeiw_tm(jk)/(nn_dttrc+1)
-         END DO
-# else
-         aeiu_tm= aeiu_tm+aeiu
-         aeiu=aeiu_tm/(nn_dttrc+1)
-
-         aeiv_tm= aeiv_tm+aeiv
-         aeiv=aeiv_tm/(nn_dttrc+1)
-
-         aeiw_tm= aeiw_tm+aeiw
-         aeiw=aeiw_tm/(nn_dttrc+1)
-# endif
-# endif
-
-         CALL lbc_lnk( sshn(:,:), 'T', 1. ) 
-         CALL lbc_lnk( sshu_n(:,:), 'U', 1. ) 
-         CALL lbc_lnk( sshv_n(:,:), 'V', 1. ) 
-         CALL lbc_lnk( sshf_n(:,:), 'F', 1. ) 
-         CALL lbc_lnk( sshb(:,:), 'T', 1. ) 
-         CALL lbc_lnk( sshu_b(:,:), 'U', 1. ) 
-         CALL lbc_lnk( sshv_b(:,:), 'V', 1. ) 
-         CALL lbc_lnk( ssha(:,:), 'T', 1. ) 
-         CALL lbc_lnk( sshu_a(:,:), 'U', 1. ) 
-         CALL lbc_lnk( sshv_a(:,:), 'V', 1. ) 
-         CALL lbc_lnk( rnf(:,:), 'T', 1. ) 
-         CALL lbc_lnk( h_rnf(:,:), 'T', 1. ) 
-         CALL lbc_lnk( rotn(:,:,:), 'F', 1. ) 
-         CALL lbc_lnk( hdivn(:,:,:), 'T', 1. ) 
-         CALL lbc_lnk( rotb(:,:,:), 'F', 1. ) 
-         CALL lbc_lnk( hdivb(:,:,:), 'T', 1. ) 
-         CALL lbc_lnk( un(:,:,:), 'U', 1. ) 
-         CALL lbc_lnk( vn(:,:,:), 'V', 1. ) 
-         CALL lbc_lnk( wn(:,:,:), 'V', 1. ) 
-         CALL lbc_lnk( tsn(:,:,:,jp_tem), 'T', 1. ) 
-         CALL lbc_lnk( tsn(:,:,:,jp_sal), 'T', 1. ) 
-         CALL lbc_lnk( avt(:,:,:), 'T', 1. ) 
-         CALL lbc_lnk( hmld(:,:), 'T', 1. ) 
-         CALL lbc_lnk( qsr(:,:), 'T', 1. ) 
-         CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
-         CALL lbc_lnk( emp(:,:), 'T', 1. ) 
-         CALL lbc_lnk( emp_b(:,:), 'T', 1. ) 
-         CALL lbc_lnk( emps(:,:), 'T', 1. ) 
-#if defined key_traldf_c3d
-         CALL lbc_lnk( ahtt(:,:,:), 'T', 1. ) 
-         CALL lbc_lnk( ahtu(:,:,:), 'U', 1. ) 
-         CALL lbc_lnk( ahtv(:,:,:), 'V', 1. ) 
-         CALL lbc_lnk( ahtw(:,:,:), 'W', 1. ) 
-#elif defined key_traldf_c2d
-         CALL lbc_lnk( ahtt(:,:), 'T', 1. ) 
-         CALL lbc_lnk( ahtu(:,:), 'U', 1. ) 
-         CALL lbc_lnk( ahtv(:,:), 'V', 1. ) 
-         CALL lbc_lnk( ahtw(:,:), 'W', 1. ) 
-#endif
-# if defined key_zdfddm
-         CALL lbc_lnk( avs, 'W', 1. ) 
-# endif
-#if defined key_traldf_eiv
-#if defined key_traldf_c3d
-         CALL lbc_lnk( aeiu(:,:,:), 'U', 1. ) 
-         CALL lbc_lnk( aeiv(:,:,:), 'V', 1. ) 
-         CALL lbc_lnk( aeiw(:,:,:), 'W', 1. ) 
-#elif defined key_traldf_c2d
-         CALL lbc_lnk( aeiu(:,:), 'U', 1. ) 
-         CALL lbc_lnk( aeiv(:,:), 'V', 1. ) 
-         CALL lbc_lnk( aeiw(:,:), 'W', 1. ) 
-#endif
-#endif
-#if defined key_ldfslp
-         CALL lbc_lnk( uslp(:,:,:), 'U', 1. ) 
-         CALL lbc_lnk( vslp(:,:,:), 'V', 1. ) 
-         CALL lbc_lnk( wslpi(:,:,:), 'W', 1. ) 
-         CALL lbc_lnk( wslpj(:,:,:), 'W', 1. ) 
-#endif
-         CALL lbc_lnk( hu(:,:), 'U', 1. ) 
-         CALL lbc_lnk( hv(:,:), 'V', 1. ) 
-         CALL lbc_lnk( hur(:,:), 'U', 1. ) 
-         CALL lbc_lnk( hvr(:,:), 'V', 1. ) 
-
+         CALL lbc_lnk( uslp  (:,:,:)       , 'U', 1. ) 
+         CALL lbc_lnk( vslp  (:,:,:)       , 'V', 1. ) 
+         CALL lbc_lnk( wslpi (:,:,:)       , 'W', 1. ) 
+         CALL lbc_lnk( wslpj (:,:,:)       , 'W', 1. ) 
+#endif
+         CALL lbc_lnk( hu    (:,:)         , 'U', 1. ) 
+         CALL lbc_lnk( hv    (:,:)         , 'V', 1. ) 
+         CALL lbc_lnk( hur   (:,:)         , 'U', 1. ) 
+         CALL lbc_lnk( hvr   (:,:)         , 'V', 1. ) 
+         !
          CALL trc_sub_ssh( kt )         ! after ssh & vertical velocity
+         !
       ENDIF
 
 
    END SUBROUTINE trc_sub_stp
-  !!-------------------------------------------------------------------
-
-  !!-------------------------------------------------------------------
+
    SUBROUTINE trc_sub_ini
       !!-------------------------------------------------------------------
@@ -578,5 +491,4 @@
       !!              Compute the averages for sub-stepping
       !!-------------------------------------------------------------------
-      INTEGER ::   jk, jn ,ji,jj   ! dummy loop indices
       INTEGER ::   ierr
       !!-------------------------------------------------------------------
@@ -590,108 +502,79 @@
       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'top_sub_alloc : unable to allocate standard ocean arrays' )
 
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               un_tm(ji,jj,jk)= un(ji,jj,jk)*fse3u(ji,jj,jk)
-               vn_tm(ji,jj,jk)= vn(ji,jj,jk)*fse3v(ji,jj,jk)
-	       wn_tm(ji,jj,jk)= wn(ji,jj,jk)*fse3w(ji,jj,jk)
-               avt_tm(ji,jj,jk)=avt(ji,jj,jk)*fse3t(ji,jj,jk)
-               tsn_tm(ji,jj,jk,jp_tem)= tsn_tm(ji,jj,jk,jp_tem)+tsn(ji,jj,jk,jp_tem)*fse3t(ji,jj,jk)
-               tsn_tm(ji,jj,jk,jp_sal)= tsn_tm(ji,jj,jk,jp_sal)+tsn(ji,jj,jk,jp_sal)*fse3t(ji,jj,jk)
-#if defined key_zdfddm
-               avs_tm(ji,jj,jk)= avs_tm(ji,jj,jk)+avs(ji,jj,jk)*fse3w(ji,jj,jk)
-#endif
+      un_tm   (:,:,:)        = un   (:,:,:)        * fse3u(:,:,:) 
+      vn_tm   (:,:,:)        = vn   (:,:,:)        * fse3v(:,:,:) 
+      wn_tm   (:,:,:)        = wn   (:,:,:)        * fse3w(:,:,:) 
+      tsn_tm  (:,:,:,jp_tem) = tsn  (:,:,:,jp_tem) * fse3t(:,:,:)  
+      tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * fse3t(:,:,:)  
+      avt_tm  (:,:,:)        = avt  (:,:,:)        * fse3w(:,:,:)  
+# if defined key_zdfddm
+      avs_tm  (:,:,:)        = avs  (:,:,:)        * fse3w(:,:,:)  
+# endif
 #if defined key_ldfslp
-               wslpi_tm(ji,jj,jk)= wslpi(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
-               wslpj_tm(ji,jj,jk)= wslpj(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
-               uslp_tm(ji,jj,jk)= uslp(ji,jj,jk)*fse3u(ji,jj,jk) !0. 
-               vslp_tm(ji,jj,jk)= vslp(ji,jj,jk)*fse3v(ji,jj,jk) !0. 
-#endif
-           END DO
-         END DO
-      END DO
-
-!   physics variables that are set after initialization:
-      wndm_tm(:,:)=0.  !wndm(:,:)
-      qsr_tm(:,:)= 0.  !qsr(:,:)
-      emp_tm(:,:)= 0.  !emp(:,:)
-      emps_tm(:,:)=0.  !emps(:,:)
- 
-      sshn_tm(:,:) = sshn(:,:) !0.
-      sshu_n_tm(:,:) = sshu_n(:,:) !0.
-      sshv_n_tm(:,:) = sshv_n(:,:) !0.
-      rnf_tm(:,:)= rnf(:,:) !0. 
-      h_rnf_tm(:,:)= h_rnf(:,:) !0.
-      hmld_tm(:,:)=hmld(:,:)
-      fr_i_tm(:,:)=fr_i(:,:)
+      wslpi_tm(:,:,:)        = wslpi(:,:,:)        * fse3w(:,:,:) 
+      wslpj_tm(:,:,:)        = wslpj(:,:,:)        * fse3w(:,:,:) 
+      uslp_tm (:,:,:)        = uslp (:,:,:)        * fse3u(:,:,:) 
+      vslp_tm (:,:,:)        = vslp (:,:,:)        * fse3v(:,:,:) 
+#endif
+
+      ! Physics variables that are set after initialization:
+      emp_tm (:,:) = 0.
+      emps_tm(:,:) = 0.
+      qsr_tm (:,:) = 0.
+      wndm_tm(:,:) = 0.
+      !
+      sshn_tm  (:,:) = sshn  (:,:) 
+      sshu_n_tm(:,:) = sshu_n(:,:) 
+      sshv_n_tm(:,:) = sshv_n(:,:) 
+      rnf_tm   (:,:) = rnf   (:,:) 
+      h_rnf_tm (:,:) = h_rnf (:,:) 
+      hmld_tm  (:,:) = hmld  (:,:)
+      fr_i_tm  (:,:) = fr_i  (:,:)
+      !
 #if defined key_traldf_c3d
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ahtt_tm(ji,jj,jk)= ahtt(ji,jj,jk)*fse3t(ji,jj,jk)
-                  ahtu_tm(ji,jj,jk)= ahtu(ji,jj,jk)*fse3u(ji,jj,jk)
-                  ahtv_tm(ji,jj,jk)= ahtv(ji,jj,jk)*fse3v(ji,jj,jk)
-                  ahtw_tm(ji,jj,jk)= ahtw(ji,jj,jk)*fse3w(ji,jj,jk)
-             END DO
-            END DO
-         END DO
+      ahtt_tm(:,:,:) = ahtt(:,:,:) * fse3t(:,:,:)
+      ahtu_tm(:,:,:) = ahtu(:,:,:) * fse3u(:,:,:)
+      ahtv_tm(:,:,:) = ahtv(:,:,:) * fse3v(:,:,:)
+      ahtw_tm(:,:,:) = ahtw(:,:,:) * fse3w(:,:,:)
 #elif defined key_traldf_c2d
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ahtt_tm(ji,jj)= ahtt(ji,jj)
-                  ahtu_tm(ji,jj)= ahtu(ji,jj)
-                  ahtv_tm(ji,jj)= ahtv(ji,jj)
-                  ahtw_tm(ji,jj)= ahtw(ji,jj)
-             END DO
-            END DO
+      ahtt_tm(:,:)   = ahtt(:,:)
+      ahtu_tm(:,:)   = ahtu(:,:)
+      ahtv_tm(:,:)   = ahtv(:,:)
+      ahtw_tm(:,:)   = ahtw(:,:)
 #elif defined key_traldf_c1d
-          DO jk = 1, jpk
-             ahtt_tm(jk)= ahtt(jk)
-             ahtu_tm(jk)= ahtu(jk)
-             ahtv_tm(jk)= ahtv(jk)
-             ahtw_tm(jk)= ahtw(jk)
-          END DO
+      ahtt_tm(:)     = ahtt(:)
+      ahtu_tm(:)     = ahtu(:)
+      ahtv_tm(:)     = ahtv(:)
+      ahtw_tm(:)     = ahtw(:)
 #else
-             ahtt_tm= ahtt
-             ahtu_tm= ahtu
-             ahtv_tm= ahtv
-             ahtw_tm= ahtw
-#endif
+      ahtt_tm        = ahtt
+      ahtu_tm        = ahtu
+      ahtv_tm        = ahtv
+      ahtw_tm        = ahtw
+#endif
+      !
 #if defined key_traldf_eiv
-# if defined key_traldf_c3d
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  aeiu_tm(ji,jj,jk)= aeiu(ji,jj,jk)*fse3u(ji,jj,jk)
-                  aeiv_tm(ji,jj,jk)= aeiv(ji,jj,jk)*fse3v(ji,jj,jk)
-                  aeiw_tm(ji,jj,jk)= aeiw(ji,jj,jk)*fse3w(ji,jj,jk)
-             END DO
-            END DO
-         END DO
-# elif defined key_traldf_c2d
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  aeiu_tm(ji,jj)= aeiu(ji,jj)
-                  aeiv_tm(ji,jj)= aeiv(ji,jj)
-                  aeiw_tm(ji,jj)= aeiw(ji,jj)
-             END DO
-            END DO
-# elif defined key_traldf_c1d
-          DO jk = 1, jpk
-             aeiu_tm(jk)= aeiu(jk)
-             aeiv_tm(jk)= aeiv(jk)
-             aeiw_tm(jk)= aeiw(jk)
-          END DO
-# else
-             aeiu_tm= aeiu
-             aeiv_tm= aeiv
-             aeiw_tm= aeiw
-# endif
-# endif
-
+#  if defined key_traldf_c3d
+      aieu_tm(:,:,:) = aieu(:,:,:) * fse3u(:,:,:)
+      aiev_tm(:,:,:) = aiev(:,:,:) * fse3v(:,:,:)
+      aiew_tm(:,:,:) = aiew(:,:,:) * fse3w(:,:,:)
+#  elif defined key_traldf_c2d
+      aieu_tm(:,:)   = aeiu(:,:)
+      aiev_tm(:,:)   = aeiv(:,:)
+      aiew_tm(:,:)   = aeiw(:,:)
+#  elif defined key_traldf_c1d
+      aeiu_tm(:)     = aeiu(:)
+      aeiv_tm(:)     = aeiv(:)
+      aeiw_tm(:)     = aeiw(:)
+#  else
+      aeiu_tm        = aeiu
+      aeiv_tm        = aeiv
+      aeiw_tm        = aeiw
+#  endif
+#endif
+      !
    END SUBROUTINE trc_sub_ini
 
-  !!-------------------------------------------------------------------
-   SUBROUTINE trc_sub_reset(kt)
+   SUBROUTINE trc_sub_reset( kt )
       !!-------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sub_reset  ***
@@ -702,102 +585,103 @@
       !!              Compute the averages for sub-stepping
       !!-------------------------------------------------------------------
-      INTEGER ::   jk, jn ,ji,jj    ! dummy loop indices
-      INTEGER ::   ierr
       INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
-      !!-------------------------------------------------------------------
-
-!   restore physics variables
-        avt(:,:,:)=avt_temp(:,:,:)
-        un(:,:,:)=un_temp(:,:,:)
-        vn(:,:,:)=vn_temp(:,:,:)
-        wn(:,:,:)=wn_temp(:,:,:)
-        sshn(:,:)=sshn_temp(:,:)
-        sshb(:,:)=sshb_temp(:,:)
-        ssha(:,:)=ssha_temp(:,:)
-        sshu_n(:,:)=sshu_n_temp(:,:)
-        sshu_b(:,:)=sshu_b_temp(:,:)
-        sshu_a(:,:)=sshu_a_temp(:,:)
-        sshv_n(:,:)=sshv_n_temp(:,:)
-        sshv_b(:,:)=sshv_b_temp(:,:)
-        sshv_a(:,:)=sshv_a_temp(:,:)
-        sshf_n(:,:)=sshf_n_temp(:,:)
+      INTEGER :: jk                 ! dummy loop indices
+      !!-------------------------------------------------------------------
+
+      !   restore physics variables
+      avt   (:,:,:)   =  avt_temp   (:,:,:)
+      un    (:,:,:)   =  un_temp    (:,:,:)
+      vn    (:,:,:)   =  vn_temp    (:,:,:)
+      wn    (:,:,:)   =  wn_temp    (:,:,:)
+      sshn  (:,:)     =  sshn_temp  (:,:)
+      sshb  (:,:)     =  sshb_temp  (:,:)
+      ssha  (:,:)     =  ssha_temp  (:,:)
+      sshu_n(:,:)     =  sshu_n_temp(:,:)
+      sshu_b(:,:)     =  sshu_b_temp(:,:)
+      sshu_a(:,:)     =  sshu_a_temp(:,:)
+      sshv_n(:,:)     =  sshv_n_temp(:,:)
+      sshv_b(:,:)     =  sshv_b_temp(:,:)
+      sshv_a(:,:)     =  sshv_a_temp(:,:)
+      sshf_n(:,:)     =  sshf_n_temp(:,:)
 #if defined key_ldfslp
-        wslpi(:,:,:)=wslpi_temp(:,:,:)
-        wslpj(:,:,:)=wslpj_temp(:,:,:)
-        uslp(:,:,:)=uslp_temp(:,:,:)
-        vslp(:,:,:)=vslp_temp(:,:,:)
-#endif
-        rnf(:,:)=rnf_temp(:,:)
-        h_rnf(:,:)=h_rnf_temp(:,:)
-        hdivn(:,:,:)=hdivn_temp(:,:,:)
-        rotn(:,:,:)=rotn_temp(:,:,:)
-        hdivb(:,:,:)=hdivb_temp(:,:,:)
-        rotb(:,:,:)=rotb_temp(:,:,:)
-        hu(:,:)=hu_temp(:,:)
-        hv(:,:)=hv_temp(:,:)
-        hur(:,:)=hur_temp(:,:)
-        hvr(:,:)=hvr_temp(:,:)
-        tsn(:,:,:,:)=tsn_temp(:,:,:,:)
+      wslpi (:,:,:)   =  wslpi_temp (:,:,:)
+      wslpj (:,:,:)   =  wslpj_temp (:,:,:)
+      uslp  (:,:,:)   =  uslp_temp  (:,:,:)
+      vslp  (:,:,:)   =  vslp_temp  (:,:,:)
+#endif
+      rnf   (:,:)     =  rnf_temp   (:,:)
+      h_rnf (:,:)     =  h_rnf_temp (:,:)
+      hdivn (:,:,:)   =  hdivn_temp (:,:,:)
+      rotn  (:,:,:)   =  rotn_temp  (:,:,:)
+      hdivb (:,:,:)   =  hdivb_temp (:,:,:)
+      rotb  (:,:,:)   =  rotb_temp  (:,:,:)
+      hu    (:,:)     =  hu_temp    (:,:)
+      hv    (:,:)     =  hv_temp    (:,:)
+      hur   (:,:)     =  hur_temp   (:,:)
+      hvr   (:,:)     =  hvr_temp   (:,:)
+      !
+      tsn   (:,:,:,:) =  tsn_temp   (:,:,:,:)
 # if defined key_zdfddm
-        avs(:,:,:)=avs_temp(:,:,:)
-# endif
-       hmld(:,:)=hmld_temp(:,:)
-       qsr(:,:)=qsr_temp(:,:)
-       fr_i(:,:)=fr_i_temp(:,:)
-       wndm(:,:)=wndm_temp(:,:)
-       emp(:,:)=emp_temp(:,:)
-       emps(:,:)=emps_temp(:,:)
-       emp_b(:,:)=emp_b_temp(:,:)
+      avs   (:,:,:)   =  avs_temp   (:,:,:)
+# endif
+      hmld  (:,:)     =  hmld_temp  (:,:)
+      qsr   (:,:)     =  qsr_temp   (:,:)
+      fr_i  (:,:)     =  fr_i_temp  (:,:)
+      wndm  (:,:)     =  wndm_temp  (:,:)
+      emp   (:,:)     =  emp_temp   (:,:)
+      emps  (:,:)     =  emps_temp  (:,:)
+      emp_b (:,:)     =  emp_b_temp (:,:)
+      !
 #if defined key_traldf_c3d
-       ahtu_temp(:,:,:)=ahtu(:,:,:)
-       ahtv_temp(:,:,:)=ahtv(:,:,:)
-       ahtw_temp(:,:,:)=ahtw(:,:,:)
-       ahtt_temp(:,:,:)=ahtt(:,:,:)
+      ahtu  (:,:,:)   =  ahtu_temp  (:,:,:)
+      ahtv  (:,:,:)   =  ahtv_temp  (:,:,:)
+      ahtw  (:,:,:)   =  ahtw_temp  (:,:,:)
+      ahtt  (:,:,:)   =  ahtt_temp  (:,:,:)
 #elif defined key_traldf_c2d
-       ahtu_temp(:,:)=ahtu(:,:)
-       ahtv_temp(:,:)=ahtv(:,:)
-       ahtw_temp(:,:)=ahtw(:,:)
-       ahtt_temp(:,:)=ahtt(:,:)
+      ahtu  (:,:)     =  ahtu_temp  (:,:)
+      ahtv  (:,:)     =  ahtv_temp  (:,:)
+      ahtw  (:,:)     =  ahtw_temp  (:,:)
+      ahtt  (:,:)     =  ahtt_temp  (:,:)
 #elif defined key_traldf_c1d
-       ahtu_temp(:)=ahtu(:)
-       ahtv_temp(:)=ahtv(:)
-       ahtw_temp(:)=ahtw(:)
-       ahtt_temp(:)=ahtt(:)
+      ahtu  (:)       =  ahtu_temp  (:)
+      ahtv  (:)       =  ahtv_temp  (:)
+      ahtw  (:)       =  ahtw_temp  (:)
+      ahtt  (:)       =  ahtt_temp  (:)
 #else
-       ahtu_temp=ahtu
-       ahtv_temp=ahtv
-       ahtw_temp=ahtw
-       ahtt_temp=ahtt
-#endif
+      ahtu            =  ahtu_temp
+      ahtv            =  ahtv_temp
+      ahtw            =  ahtw_temp
+      ahtt            =  ahtt_temp
+#endif
+      !
 #if defined key_traldf_eiv
 #if defined key_traldf_c3d
-       aeiu_temp(:,:,:)=aeiu(:,:,:)
-       aeiv_temp(:,:,:)=aeiv(:,:,:)
-       aeiw_temp(:,:,:)=aeiw(:,:,:)
+      aeiu  (:,:,:)  =  aeiu_temp(:,:,:)
+      aeiv  (:,:,:)  =  aeiv_temp(:,:,:)
+      aeiw  (:,:,:)  =  aeiw_temp(:,:,:)
 #elif defined key_traldf_c2d
-       aeiu_temp(:,:)=aeiu(:,:)
-       aeiv_temp(:,:)=aeiv(:,:)
-       aeiw_temp(:,:)=aeiw(:,:)
+      aeiu  (:,:)    =  aeiu_temp(:,:)
+      aeiv  (:,:)    =  aeiv_temp(:,:)
+      aeiw  (:,:)    =  aeiw_temp(:,:)
 #elif defined key_traldf_c1d
-       aeiu_temp(:)=aeiu(:)
-       aeiv_temp(:)=aeiv(:)
-       aeiw_temp(:)=aeiw(:)
+      aeiu  (:)      =  aeiu_temp(:)
+      aeiv  (:)      =  aeiv_temp(:)
+      aeiw  (:)      =  aeiw_temp(:)
 #else
-       aeiu_temp=aeiu
-       aeiv_temp=aeiv
-       aeiw_temp=aeiw
-#endif
-#endif        
-
-
-       !                                           !------------------------------!
-       DO jk = 1, jpk
-          fse3t(:,:,jk)= e3t_temp(:,:,jk) !0. 
-          fse3u(:,:,jk)= e3u_temp(:,:,jk) !0. 
-          fse3v(:,:,jk)= e3v_temp(:,:,jk)!0. 
-          fse3w(:,:,jk)= e3w_temp(:,:,jk) !0. 
-       END DO
-       IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
-       !                                           !------------------------------!
+      aeiu           =  aeiu_temp
+      aeiv           =  aeiv_temp
+      aeiw           =  aeiw_temp
+#endif
+#endif 
+      !                                      
+      DO jk = 1, jpk
+         fse3t(:,:,jk)= e3t_temp(:,:,jk) 
+         fse3u(:,:,jk)= e3u_temp(:,:,jk) 
+         fse3v(:,:,jk)= e3v_temp(:,:,jk) 
+         fse3w(:,:,jk)= e3w_temp(:,:,jk) 
+      END DO
+      !                                           !------------------------------!
+      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
+        !                                           !------------------------------!
          DO jk = 1, jpkm1
             fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
@@ -813,111 +697,91 @@
             fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
          END DO
+         !
       ENDIF
 
-!   start new averages
+      ! Start new averages
       DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               un_tm(ji,jj,jk)=un(ji,jj,jk)*fse3u(ji,jj,jk)  !0.0
-               vn_tm(ji,jj,jk)=vn(ji,jj,jk)*fse3v(ji,jj,jk)  !0.0
-               wn_tm(ji,jj,jk)=wn(ji,jj,jk)*fse3w(ji,jj,jk)  !0.0
-               avt_tm(ji,jj,jk)=avt(ji,jj,jk)*fse3t(ji,jj,jk)  !0.0
-               tsn_tm(ji,jj,jk,jp_tem)=tsn(ji,jj,jk,jp_tem)*fse3t(ji,jj,jk)  !0.0
-               tsn_tm(ji,jj,jk,jp_sal)=tsn(ji,jj,jk,jp_sal)*fse3t(ji,jj,jk)  !0.0
+         un_tm   (:,:,jk)        = un   (:,:,jk)        * fse3u(:,:,jk) 
+         vn_tm   (:,:,jk)        = vn   (:,:,jk)        * fse3v(:,:,jk) 
+         wn_tm   (:,:,jk)        = wn   (:,:,jk)        * fse3w(:,:,jk) 
+         tsn_tm  (:,:,jk,jp_tem) = tsn  (:,:,jk,jp_tem) * fse3t(:,:,jk)  
+         tsn_tm  (:,:,jk,jp_sal) = tsn  (:,:,jk,jp_sal) * fse3t(:,:,jk)  
+         avt_tm  (:,:,jk)        = avt  (:,:,jk)        * fse3w(:,:,jk)  
 # if defined key_zdfddm
-               avs_tm==avs(ji,jj,jk)*fse3w(ji,jj,jk)  !0.0
+         avs_tm  (:,:,jk)        = avs  (:,:,jk)        * fse3w(:,:,jk)  
 # endif
 #if defined key_ldfslp
-               wslpi_tm(ji,jj,jk)= wslpi(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
-               wslpj_tm(ji,jj,jk)=wslpj(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
-               uslp_tm(ji,jj,jk)=uslp(ji,jj,jk)*fse3u(ji,jj,jk) !0. 
-               vslp_tm(ji,jj,jk)=vslp(ji,jj,jk)*fse3v(ji,jj,jk) !0. 
-#endif
-            END DO
-         END DO
+         wslpi_tm(:,:,jk)        = wslpi(:,:,jk)        * fse3w(:,:,jk) 
+         wslpj_tm(:,:,jk)        = wslpj(:,:,jk)        * fse3w(:,:,jk) 
+         uslp_tm (:,:,jk)        = uslp (:,:,jk)        * fse3u(:,:,jk) 
+         vslp_tm (:,:,jk)        = vslp (:,:,jk)        * fse3v(:,:,jk) 
+#endif
       END DO
-      sshn_tm(:,:) = sshn(:,:) !0.
-      sshu_n_tm(:,:) = sshu_n(:,:) !0.
-      sshv_n_tm(:,:) = sshv_n(:,:) !0.
-      rnf_tm(:,:)= rnf(:,:) !0. 
-      h_rnf_tm(:,:)= h_rnf(:,:) !0.
-      hmld_tm(:,:)=hmld(:,:)
-      fr_i_tm(:,:)=fr_i(:,:)
-      sshb_hold(:,:)=sshn(:,:)
-      sshu_b_hold(:,:)=sshu_n(:,:)
-      sshv_b_hold(:,:)=sshv_n(:,:)
-      emp_b_hold(:,:)=emp(:,:)
-      emp_tm(:,:)=emp(:,:)
-      emps_tm(:,:)=emps(:,:)
-      qsr_tm(:,:)=qsr(:,:)
-      wndm_tm(:,:)=wndm(:,:)
+      !
+      sshn_tm    (:,:) = sshn  (:,:) 
+      sshu_n_tm  (:,:) = sshu_n(:,:) 
+      sshv_n_tm  (:,:) = sshv_n(:,:) 
+      rnf_tm     (:,:) = rnf   (:,:) 
+      h_rnf_tm   (:,:) = h_rnf (:,:) 
+      hmld_tm    (:,:) = hmld  (:,:)
+      fr_i_tm    (:,:) = fr_i  (:,:)
+      sshb_hold  (:,:) = sshn  (:,:)
+      sshu_b_hold(:,:) = sshu_n(:,:)
+      sshv_b_hold(:,:) = sshv_n(:,:)
+      emp_b_hold (:,:) = emp   (:,:)
+      emp_tm     (:,:) = emp   (:,:)
+      emps_tm    (:,:) = emps  (:,:)
+      qsr_tm     (:,:) = qsr   (:,:)
+      wndm_tm    (:,:) = wndm  (:,:)
+      !
 #if defined key_traldf_c3d
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ahtt_tm(ji,jj,jk)= ahtt(ji,jj,jk)*fse3t(ji,jj,jk)
-               ahtu_tm(ji,jj,jk)= ahtu(ji,jj,jk)*fse3u(ji,jj,jk)
-               ahtv_tm(ji,jj,jk)= ahtv(ji,jj,jk)*fse3v(ji,jj,jk)
-               ahtw_tm(ji,jj,jk)= ahtw(ji,jj,jk)*fse3w(ji,jj,jk)
-            END DO
-         END DO
+      DO jk = 1, jpkm1
+         ahtt_tm(:,:,jk) = ahtt(:,:,jk) * fse3t(:,:,jk)
+         ahtu_tm(:,:,jk) = ahtu(:,:,jk) * fse3u(:,:,jk)
+         ahtv_tm(:,:,jk) = ahtv(:,:,jk) * fse3v(:,:,jk)
+         ahtw_tm(:,:,jk) = ahtw(:,:,jk) * fse3w(:,:,jk)
       END DO
 #elif defined key_traldf_c2d
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ahtt_tm(ji,jj)= ahtt(ji,jj)
-            ahtu_tm(ji,jj)= ahtu(ji,jj)
-            ahtv_tm(ji,jj)= ahtv(ji,jj)
-            ahtw_tm(ji,jj)= ahtw(ji,jj)
-         END DO
+      ahtt_tm(:,:)   = ahtt(:,:)
+      ahtu_tm(:,:)   = ahtu(:,:)
+      ahtv_tm(:,:)   = ahtv(:,:)
+      ahtw_tm(:,:)   = ahtw(:,:)
+#elif defined key_traldf_c1d
+      ahtt_tm(:)     = ahtt(:)
+      ahtu_tm(:)     = ahtu(:)
+      ahtv_tm(:)     = ahtv(:)
+      ahtw_tm(:)     = ahtw(:)
+#else
+      ahtt_tm        = ahtt
+      ahtu_tm        = ahtu
+      ahtv_tm        = ahtv
+      ahtw_tm        = ahtw
+#endif
+      !
+#if defined key_traldf_eiv
+#  if defined key_traldf_c3d
+      DO jk = 1, jpk
+         aieu_tm(:,:,jk) = aieu(:,:,jk) * fse3u(:,:,jk)
+         aiev_tm(:,:,jk) = aiev(:,:,jk) * fse3v(:,:,jk)
+         aiew_tm(:,:,jk) = aiew(:,:,jk) * fse3w(:,:,jk)
       END DO
-#elif defined key_traldf_c1d
-      DO jk = 1, jpk
-         ahtt_tm(jk)= ahtt(jk)
-         ahtu_tm(jk)= ahtu(jk)
-         ahtv_tm(jk)= ahtv(jk)
-         ahtw_tm(jk)= ahtw(jk)
-      END DO
-#else
-      ahtt_tm= ahtt
-      ahtu_tm= ahtu
-      ahtv_tm= ahtv
-      ahtw_tm= ahtw
-#endif
-#if defined key_traldf_eiv
-# if defined key_traldf_c3d
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               aeiu_tm(ji,jj,jk)= aeiu(ji,jj,jk)*fse3u(ji,jj,jk)
-               aeiv_tm(ji,jj,jk)= aeiv(ji,jj,jk)*fse3v(ji,jj,jk)
-               aeiw_tm(ji,jj,jk)= aeiw(ji,jj,jk)*fse3w(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-# elif defined key_traldf_c2d
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            aeiu_tm(ji,jj)= aeiu(ji,jj)
-            aeiv_tm(ji,jj)= aeiv(ji,jj)
-            aeiw_tm(ji,jj)= aeiw(ji,jj)
-         END DO
-      END DO
-# elif defined key_traldf_c1d
-      DO jk = 1, jpk
-         aeiu_tm(jk)= aeiu(jk)
-         aeiv_tm(jk)= aeiv(jk)
-         aeiw_tm(jk)= aeiw(jk)
-      END DO
-# else
-      aeiu_tm= aeiu
-      aeiv_tm= aeiv
-      aeiw_tm= aeiw
-# endif
-# endif
-
-
+#  elif defined key_traldf_c2d
+      aieu_tm(:,:)   = aeiu(:,:)
+      aiev_tm(:,:)   = aeiv(:,:)
+      aiew_tm(:,:)   = aeiw(:,:)
+#  elif defined key_traldf_c1d
+      aeiu_tm(:)     = aeiu(:)
+      aeiv_tm(:)     = aeiv(:)
+      aeiw_tm(:)     = aeiw(:)
+#  else
+      aeiu_tm        = aeiu
+      aeiv_tm        = aeiv
+      aeiw_tm        = aeiw
+#  endif
+#endif
+      !
    END SUBROUTINE trc_sub_reset
-  !!-------------------------------------------------------------------
+
+
    SUBROUTINE trc_sub_ssh( kt ) 
       !!----------------------------------------------------------------------
@@ -940,6 +804,6 @@
       !!----------------------------------------------------------------------
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace
       USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1 , z2d => wrk_2d_2   ! 2D workspace
+      USE wrk_nemo, ONLY:   z3d => wrk_3d_10                      ! 3D workspace
       !
       INTEGER, INTENT(in) ::   kt   ! time step
@@ -949,7 +813,8 @@
       !!----------------------------------------------------------------------
 
-      IF( wrk_in_use(2, 1,2) ) THEN
-         CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN
+      IF( wrk_in_use(2, 1,2) .OR. wrk_in_use(3,10) ) THEN
+         CALL ctl_stop('trc_sub_ssh: requested workspace arrays unavailable')   ;   RETURN
       ENDIF
+
 
       IF( kt == nittrc000 ) THEN
@@ -960,5 +825,4 @@
          !
          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
          !
       ENDIF
@@ -1052,10 +916,10 @@
 
       !
-      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
+      IF( wrk_not_released(2, 1,2) .OR. wrk_not_released(3, 10) ) THEN
+          CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
+      ENDIF
       !
    END SUBROUTINE trc_sub_ssh
-  !!-------------------------------------------------------------------
-
-  !!-------------------------------------------------------------------
+
    INTEGER FUNCTION trc_sub_alloc()
       !!-------------------------------------------------------------------
@@ -1063,5 +927,5 @@
       !!-------------------------------------------------------------------
       USE lib_mpp, ONLY: ctl_warn
-      INTEGER ::         ierr
+      INTEGER ::  ierr
       !!-------------------------------------------------------------------
       !
