Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90	(revision 6060)
@@ -11,21 +11,21 @@
    !!   tra_ldf_triad : update the tracer trend with the iso-neutral laplacian triad-operator
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE trc_oce         ! share passive tracers/Ocean variables
-   USE zdf_oce         ! ocean vertical physics
-   USE ldftra          ! lateral physics: eddy diffusivity
-   USE ldfslp          ! lateral physics: iso-neutral slopes
-   USE traldf_iso      ! lateral diffusion (Madec operator)         (tra_ldf_iso routine)
-   USE diaptr          ! poleward transport diagnostics
-   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
+   USE oce            ! ocean dynamics and active tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE trc_oce        ! share passive tracers/Ocean variables
+   USE zdf_oce        ! ocean vertical physics
+   USE ldftra         ! lateral physics: eddy diffusivity
+   USE ldfslp         ! lateral physics: iso-neutral slopes
+   USE traldf_iso     ! lateral diffusion (Madec operator)         (tra_ldf_iso routine)
+   USE diaptr         ! poleward transport diagnostics
+   USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
    !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O library
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -37,5 +37,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -113,5 +112,4 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
       ENDIF
-      !
       !                                               ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdttra(1)      ! at nit000   (Euler)
@@ -123,5 +121,5 @@
       ELSE                    ;   zsign = -1._wp
       ENDIF
-                  
+      !    
       !!----------------------------------------------------------------------
       !!   0 - calculate  ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw
@@ -142,15 +140,15 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
-                        zbu   = e1e2u(ji,jj) * fse3u(ji,jj,jk)
+                        ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
+                        zbu   = e1e2u(ji,jj) * e3u_n(ji,jj,jk)
                         zah   = 0.25_wp * pahu(ji,jj,jk)
                         zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                         ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
-                        zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+                        zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
                         zslope2 = zslope2 *zslope2
                         ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2
                         akz     (ji+ip,jj,jk+kp) = akz     (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj)       &
                            &                                                      * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
-                        !
+                           !
                        IF( ln_ldfeiv_dia )   zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)   &
                            &                                       + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew
@@ -166,11 +164,11 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp)
-                        zbv   = e1e2v(ji,jj) * fse3v(ji,jj,jk)
+                        ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp)
+                        zbv   = e1e2v(ji,jj) * e3v_n(ji,jj,jk)
                         zah   = 0.25_wp * pahv(ji,jj,jk)
                         zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                         ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
                         !    (do this by *adding* gradient of depth)
-                        zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+                        zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
                         zslope2 = zslope2 * zslope2
                         ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2
@@ -193,5 +191,5 @@
                      DO ji = 1, fs_jpim1
                         akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )  )
+                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  )
                      END DO
                   END DO
@@ -201,5 +199,5 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
+                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
@@ -250,5 +248,5 @@
             ENDIF
          ENDIF
-
+         !
          !!----------------------------------------------------------------------
          !!   II - horizontal trend  (full)
@@ -256,5 +254,4 @@
          !
          DO jk = 1, jpkm1
-            !
             !                    !==  Vertical tracer gradient at level jk and jk+1
             zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
@@ -274,10 +271,10 @@
                            ze1ur = r1_e1u(ji,jj)
                            zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
+                           ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
                            zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
                            zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                            zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp)
-
-                           zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk)
+                           !
+                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk)
                            ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
                            zah = pahu(ji,jj,jk)
@@ -290,5 +287,5 @@
                   END DO
                END DO
-
+               !
                DO jp = 0, 1
                   DO kp = 0, 1
@@ -297,9 +294,9 @@
                            ze2vr = r1_e2v(ji,jj)
                            zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
+                           ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp)
                            zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
                            zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                            zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk)
                            ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked...
                            zah = pahv(ji,jj,jk)
@@ -312,7 +309,7 @@
                   END DO
                END DO
-
+               !
             ELSE
-
+               !
                DO ip = 0, 1               !==  Horizontal & vertical fluxes
                   DO kp = 0, 1
@@ -321,14 +318,14 @@
                            ze1ur = r1_e1u(ji,jj)
                            zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
+                           ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
                            zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
                            zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                            zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
-
-                           zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk)
+                           !
+                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk)
                            ! ln_botmix_triad is .F. mask zah for bottom half cells
                            zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
                            zah_slp  = zah * zslope_iso
-                           IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! fsaeit(ji+ip,jj,jk)*zslope_skew
+                           IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! aeit(ji+ip,jj,jk)*zslope_skew
                            zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
                            ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
@@ -337,5 +334,5 @@
                   END DO
                END DO
-
+               !
                DO jp = 0, 1
                   DO kp = 0, 1
@@ -344,13 +341,13 @@
                            ze2vr = r1_e2v(ji,jj)
                            zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
+                           ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp)
                            zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
                            zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                            zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk)
                            ! ln_botmix_triad is .F. mask zah for bottom half cells
                            zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ????
                            zah_slp = zah * zslope_iso
-                           IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! fsaeit(ji,jj+jp,jk)*zslope_skew
+                           IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! aeit(ji,jj+jp,jk)*zslope_skew
                            zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
                            ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
@@ -365,5 +362,5 @@
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
                      &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
-                     &                                        / (  e1e2t(ji,jj) * fse3t(ji,jj,jk)  )
+                     &                                        / (  e1e2t(ji,jj) * e3t_n(ji,jj,jk)  )
                END DO
             END DO
@@ -376,5 +373,5 @@
                DO jj = 1, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)   &
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)   &
                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
@@ -388,5 +385,5 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)             &
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)             &
                            &                            * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
                      END DO
@@ -397,5 +394,5 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)                      &
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)                      &
                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
@@ -410,5 +407,5 @@
                DO ji = fs_2, fs_jpim1  ! vector opt.
                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
-                     &                                        / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
+                     &                                        / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
                END DO
             END DO
