Index: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 4990)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 5836)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                     ***  MODULE  dynldf_iso  ***
-   !! Ocean dynamics:  lateral viscosity trend
+   !! Ocean dynamics:   lateral viscosity trend (rotated laplacian operator)
    !!======================================================================
    !! History :  OPA  !  97-07  (G. Madec)  Original code
@@ -8,8 +8,8 @@
    !!             -   !  2004-08  (C. Talandier) New trends organization
    !!            2.0  !  2005-11  (G. Madec)  s-coordinate: horizontal diffusion
+   !!            3.7  !  2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
+   !!                 !                                   add velocity dependent coefficient and optional read in file
    !!----------------------------------------------------------------------
-#if defined key_ldfslp   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_ldfslp'                      slopes of the direction of mixing
+
    !!----------------------------------------------------------------------
    !!   dyn_ldf_iso  : update the momentum trend with the horizontal part
@@ -19,12 +19,12 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE ldfdyn_oce      ! ocean dynamics lateral physics
-   USE ldftra_oce      ! ocean tracer   lateral physics
+   USE ldfdyn          ! lateral diffusion: eddy viscosity coef.
+   USE ldftra          ! lateral physics: eddy diffusivity
    USE zdf_oce         ! ocean vertical physics
    USE ldfslp          ! iso-neutral slopes 
    !
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
    USE wrk_nemo        ! Memory Allocation
@@ -42,5 +42,4 @@
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-#  include "ldfdyn_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -83,13 +82,13 @@
       !!      horizontal fluxes associated with the rotated lateral mixing:
       !!      u-component:
-      !!         ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t  di[ ub ]
-      !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
-      !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
-      !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ]
+      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ ub ]
+      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
+      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
+      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(ub)) ]
       !!      v-component:
-      !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ]
-      !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ]
-      !!         zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
-      !!               -      ahmt       e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
+      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vb ]
+      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vb)) ]
+      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
+      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
       !!      take the horizontal divergence of the fluxes:
       !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
@@ -106,10 +105,9 @@
       !!      of the rotated operator in dynzdf module
       !!----------------------------------------------------------------------
-      !
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zabe1, zabe2, zcof1, zcof2                       ! local scalars
-      REAL(wp) ::   zmskt, zmskf, zbu, zbv, zuah, zvah               !   -      -
+      REAL(wp) ::   zmskt, zmskf                                     !   -      -
       REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      -
       REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
@@ -130,5 +128,6 @@
       ENDIF
 
-      ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum
+!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED
+      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
       IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
          !
@@ -136,8 +135,8 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
-                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
-                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
-                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
+                  uslp (ji,jj,jk) = - ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
+                  vslp (ji,jj,jk) = - ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
+                  wslpi(ji,jj,jk) = - ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
+                  wslpj(ji,jj,jk) = - ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
                END DO
             END DO
@@ -184,11 +183,26 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, jpi   ! vector opt.
-                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
-
-                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
-                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
-
-                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
+                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) * r1_e1t(ji,jj)
+
+                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
+                     &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
+
+                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
    
+                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    &
+                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
+                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         ELSE                   ! other coordinate system (zco or sco) : e3t
+            DO jj = 2, jpjm1
+               DO ji = fs_2, jpi   ! vector opt.
+                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * fse3t(ji,jj,jk) * r1_e1t(ji,jj)
+
+                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
+                     &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
+
+                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
+
                   ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
                      &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
@@ -196,19 +210,4 @@
                END DO
             END DO
-         ELSE                   ! other coordinate system (zco or sco) : e3t
-            DO jj = 2, jpjm1
-               DO ji = fs_2, jpi   ! vector opt.
-                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
-
-                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
-                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
-
-                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
-
-                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
-                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
-                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
-               END DO
-            END DO
          ENDIF
 
@@ -216,10 +215,10 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
-
-               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
-                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
-
-               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
+               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * fse3f(ji,jj,jk) * r1_e2f(ji,jj)
+
+               zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
+                  &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp )
+
+               zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
 
                zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
@@ -237,14 +236,14 @@
          DO jj = 2, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
-
-               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
-                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
-
-               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
-
-               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
-                  &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     &
-                  &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
+               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * fse3f(ji,jj,jk) * r1_e1f(ji,jj)
+
+               zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
+                  &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
+
+               zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
+
+               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    &
+                  &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
+                  &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
             END DO
          END DO
@@ -254,10 +253,25 @@
             DO jj = 2, jpj
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
+                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) * r1_e2t(ji,jj)
+
+                  zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
+                     &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
+
+                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
+
+                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    &
+                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
+                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         ELSE                   ! other coordinate system (zco or sco) : e3t
+            DO jj = 2, jpj
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * fse3t(ji,jj,jk) * r1_e2t(ji,jj)
 
                   zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
                      &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
 
-                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
+                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
 
                   zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
@@ -266,19 +280,4 @@
                END DO
             END DO
-         ELSE                   ! other coordinate system (zco or sco) : e3t
-            DO jj = 2, jpj
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
-
-                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
-                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
-
-                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
-
-                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
-                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
-                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
-               END DO
-            END DO
          ENDIF
 
@@ -286,18 +285,10 @@
          ! Second derivative (divergence) and add to the general trend
          ! -----------------------------------------------------------
-
          DO jj = 2, jpjm1
-            DO ji = 2, jpim1          !! Question vectop possible??? !!bug
-               ! volume elements
-               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-               ! horizontal component of isopycnal momentum diffusive trends
-               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   &
-                  &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu
-               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   &
-                  &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv
-               ! add the trends to the general trends
-               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
-               va (ji,jj,jk) = va (ji,jj,jk) + zvah
+            DO ji = 2, jpim1          !!gm Question vectop possible??? !!bug
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj  )    &
+                  &                          + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj  ) - zivf(ji-1,jj)    &
+                  &                          + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
             END DO
          END DO
@@ -358,13 +349,13 @@
          DO jk = 2, jpkm1
             DO ji = 2, jpim1
-               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
-
+               zcoef0= 0.5 * rn_aht_0 * umask(ji,jj,jk)
+               !
                zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
                zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
-
+               !
                zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
                              + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
-               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
-                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
+               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)   &
+                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ), 1. )
 
                zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
@@ -376,6 +367,6 @@
                                        +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
                ! update avmu (add isopycnal vertical coefficient to avmu)
-               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
-               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
+               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
+               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0
             END DO
          END DO
@@ -384,5 +375,5 @@
          DO jk = 2, jpkm1
             DO ji = 2, jpim1
-               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
+               zcoef0 = 0.5 * rn_aht_0 * vmask(ji,jj,jk)
 
                zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
@@ -398,10 +389,10 @@
                ! vertical flux on v field
                zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
-                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
-                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
-                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
+                  &                    +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
+                  &        + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
+                  &                    +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
                ! update avmv (add isopycnal vertical coefficient to avmv)
-               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
-               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
+               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
+               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0
             END DO
          END DO
@@ -412,13 +403,6 @@
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
-               ! volume elements
-               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-               ! part of the k-component of isopycnal momentum diffusive trends
-               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
-               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
-               ! add the trends to the general trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
-               va(ji,jj,jk) = va(ji,jj,jk) + zvav
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
             END DO
          END DO
@@ -432,15 +416,4 @@
    END SUBROUTINE dyn_ldf_iso
 
-# else
-   !!----------------------------------------------------------------------
-   !!   Dummy module                           NO rotation of mixing tensor
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt
-   END SUBROUTINE dyn_ldf_iso
-#endif
-
    !!======================================================================
 END MODULE dynldf_iso
