Index: branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
===================================================================
--- branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90	(revision 2783)
+++ branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90	(revision 2840)
@@ -4,5 +4,5 @@
    !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
    !!======================================================================
-   !! History : 3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  
+   !! History : 3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)
    !!                !          Griffies operator version adapted from traldf_iso.F90
    !!----------------------------------------------------------------------
@@ -11,6 +11,6 @@
    !!   'key_ldfslp'               slope of the lateral diffusive direction
    !!----------------------------------------------------------------------
-   !!   tra_ldf_iso_grif  : update the tracer trend with the horizontal component  
-   !!                       of the Griffies iso-neutral laplacian operator 
+   !!   tra_ldf_iso_grif  : update the tracer trend with the horizontal component
+   !!                       of the Griffies iso-neutral laplacian operator
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and active tracers
@@ -53,9 +53,9 @@
       !!                  ***  ROUTINE tra_ldf_iso_grif  ***
       !!
-      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
-      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 
+      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
+      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
       !!      add it to the general trend of tracer equation.
       !!
-      !! ** Method  :   The horizontal component of the lateral diffusive trends 
+      !! ** Method  :   The horizontal component of the lateral diffusive trends
       !!      is provided by a 2nd order operator rotated along neural or geopo-
       !!      tential surfaces to which an eddy induced advection can be added
@@ -67,5 +67,5 @@
       !!
       !!      2nd part :  horizontal fluxes of the lateral mixing operator
-      !!      ========    
+      !!      ========
       !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ]
       !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ]
@@ -99,5 +99,5 @@
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
       REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef
       !
@@ -143,5 +143,5 @@
 
       !!----------------------------------------------------------------------
-      !!   0 - calculate  ah_wslp2, psix_eiv, psiy_eiv 
+      !!   0 - calculate  ah_wslp2, psix_eiv, psiy_eiv
       !!----------------------------------------------------------------------
 
@@ -159,9 +159,12 @@
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1
+                     ze1ur = 1._wp / e1u(ji,jj)
                      ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
                      zbu   = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
                      zah   = fsahtu(ji,jj,jk)                                  !  fsaht(ji+ip,jj,jk)
                      zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
-                     zslope2 = zslope_skew - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+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) ) * ze1ur * umask(ji,jj,jk+kp)
                      zslope2 = zslope2 *zslope2
                      ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp)    &
@@ -182,9 +185,12 @@
                DO jj = 1, jpjm1
                   DO ji=1,fs_jpim1
+                     ze2vr = 1._wp / e2v(ji,jj)
                      ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp)
                      zbv   = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-                     zah   = fsahtu(ji,jj,jk)                                       !fsaht(ji,jj+jp,jk)
+                     zah   = fsahtv(ji,jj,jk)                                       !fsaht(ji,jj+jp,jk)
                      zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
-                     zslope2 = zslope_skew - ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+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) ) * ze2vr * vmask(ji,jj,jk+kp)
                      zslope2 = zslope2 * zslope2
                      ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp)   &
@@ -207,5 +213,5 @@
          zftu(:,:,:) = 0._wp
          zftv(:,:,:) = 0._wp
-         !                                               
+         !
          DO jk = 1, jpkm1                          !==  before lateral T & S gradients at T-level jk  ==!
             DO jj = 1, jpjm1
@@ -216,5 +222,5 @@
             END DO
          END DO
-         IF( ln_zps ) THEN                               ! partial steps: correction at the last level
+         IF( ln_zps.and.ln_grad_zps ) THEN              ! partial steps: correction at the last level
 # if defined key_vectopt_loop
             DO jj = 1, 1
@@ -224,6 +230,6 @@
                DO ji = 1, jpim1
 # endif
-                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)          
-                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)      
+                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
+                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
                END DO
             END DO
@@ -244,52 +250,96 @@
             ENDIF
 
-            IF( .NOT. l_triad_iso ) THEN
-               triadi = triadi_g
-               triadj = triadj_g
-            ENDIF
-
-            DO ip = 0, 1              !==  Horizontal & vertical fluxes
-               DO kp = 0, 1
-                  DO jj = 1, jpjm1
-                     DO ji = 1, fs_jpim1
-                        ze1ur = 1._wp / e1u(ji,jj)
-                        zdxt = zdit(ji,jj,jk) * ze1ur
-                        ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
-                        zdzt  = zdkt(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 * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-                        zah = fsahtu(ji,jj,jk)   !*umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)           ===>>  ????
-                        zah_slp  = zah * zslope_iso
-                        zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew    !fsaeit(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
+
+            IF (ln_botmix) THEN
+               DO ip = 0, 1              !==  Horizontal & vertical fluxes
+                  DO kp = 0, 1
+                     DO jj = 1, jpjm1
+                        DO ji = 1, fs_jpim1
+                           ze1ur = 1._wp / e1u(ji,jj)
+                           zdxt = zdit(ji,jj,jk) * ze1ur
+                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
+                           zdzt  = zdkt(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 * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
+                           ! ln_botmix is .T. don't mask zah for bottom half cells
+                           zah = fsahtu(ji,jj,jk)   !*umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)           ===>>  ????
+                           zah_slp  = zah * zslope_iso
+                           zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew    !fsaeit(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
+                        END DO
                      END DO
                   END DO
                END DO
-            END DO
-
-            DO jp = 0, 1
-               DO kp = 0, 1
-                  DO jj = 1, jpjm1
-                     DO ji = 1, fs_jpim1
-                        ze2vr = 1._wp / e2v(ji,jj)
-                        zdyt = zdjt(ji,jj,jk) * ze2vr
-                        ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
-                        zdzt = zdkt(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 * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-                        zah = fsahtv(ji,jj,jk)        !*vmask(ji,jj,jk+kp)         !fsaht(ji,jj+jp,jk)
-                        zah_slp = zah * zslope_iso
-                        zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew    !fsaeit(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
+
+               DO jp = 0, 1
+                  DO kp = 0, 1
+                     DO jj = 1, jpjm1
+                        DO ji = 1, fs_jpim1
+                           ze2vr = 1._wp / e2v(ji,jj)
+                           zdyt = zdjt(ji,jj,jk) * ze2vr
+                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
+                           zdzt = zdkt(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 * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
+                           ! ln_botmix is .T. don't mask zah for bottom half cells
+                           zah = fsahtv(ji,jj,jk)        !*vmask(ji,jj,jk+kp)         !fsaht(ji,jj+jp,jk)
+                           zah_slp = zah * zslope_iso
+                           zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew    !fsaeit(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
+                        END DO
                      END DO
                   END DO
                END DO
-            END DO
-
+            ELSE
+               DO ip = 0, 1              !==  Horizontal & vertical fluxes
+                  DO kp = 0, 1
+                     DO jj = 1, jpjm1
+                        DO ji = 1, fs_jpim1
+                           ze1ur = 1._wp / e1u(ji,jj)
+                           zdxt = zdit(ji,jj,jk) * ze1ur
+                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
+                           zdzt  = zdkt(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 * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
+                           ! ln_botmix is .F. mask zah for bottom half cells
+                           zah = fsahtu(ji,jj,jk) * umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)           ===>>  ????
+                           zah_slp  = zah * zslope_iso
+                           zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew    !fsaeit(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
+                        END DO
+                     END DO
+                  END DO
+               END DO
+
+               DO jp = 0, 1
+                  DO kp = 0, 1
+                     DO jj = 1, jpjm1
+                        DO ji = 1, fs_jpim1
+                           ze2vr = 1._wp / e2v(ji,jj)
+                           zdyt = zdjt(ji,jj,jk) * ze2vr
+                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
+                           zdzt = zdkt(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 * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
+                           ! ln_botmix is .F. mask zah for bottom half cells
+                           zah = fsahtv(ji,jj,jk) * vmask(ji,jj,jk+kp)         !fsaht(ji,jj+jp,jk)
+                           zah_slp = zah * zslope_iso
+                           zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew    !fsaeit(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
+                        END DO
+                     END DO
+                  END DO
+               END DO
+            END IF
             !                        !==  divergence and add to the general trend  ==!
             DO jj = 2 , jpjm1
@@ -320,10 +370,10 @@
 #if defined key_diaar5
          IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
-            z2d(:,:) = 0._wp 
-            zztmp = rau0 * rcp 
+            z2d(:,:) = 0._wp
+            zztmp = rau0 * rcp
             DO jk = 1, jpkm1
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
+                     z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)
                   END DO
                END DO
@@ -332,9 +382,9 @@
             CALL lbc_lnk( z2d, 'U', -1. )
             CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
-            z2d(:,:) = 0._wp 
+            z2d(:,:) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
+                     z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)
                   END DO
                END DO
