Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 3625)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 3634)
@@ -67,6 +67,9 @@
       NAMELIST/namdyn_ldf/ ln_dynldf_lap  , ln_dynldf_bilap,                  &
          &                 ln_dynldf_level, ln_dynldf_hor  , ln_dynldf_iso,   &
-         &                 rn_ahm_0_lap   , rn_ahmb_0      , rn_ahm_0_blp
-      !!----------------------------------------------------------------------
+         &                 rn_ahm_0_lap   , rn_ahmb_0      , rn_ahm_0_blp ,   &
+         &                 rn_cmsmag_1    , rn_cmsmag_2    , rn_cmsh,         &
+         &                 rn_ahm_m_lap   , rn_ahm_m_blp
+
+   !!----------------------------------------------------------------------
 
       REWIND( numnam )                  ! Read Namelist namdyn_ldf : Lateral physics
@@ -86,4 +89,7 @@
          WRITE(numout,*) '      background viscosity                    rn_ahmb_0       = ', rn_ahmb_0
          WRITE(numout,*) '      horizontal bilaplacian eddy viscosity   rn_ahm_0_blp    = ', rn_ahm_0_blp
+         WRITE(numout,*) '      upper limit for laplacian eddy visc     rn_ahm_m_lap    = ', rn_ahm_m_lap
+         WRITE(numout,*) '      upper limit for bilap eddy viscosity    rn_ahm_m_blp    = ', rn_ahm_m_blp
+
       ENDIF
 
@@ -142,4 +148,9 @@
       IF(lwp) WRITE(numout,*) '        ahm1 = ahm2 = ahm0 =  ',ahm0
 #endif
+     nkahm_smag = 0
+#if defined key_dynldf_smag
+     nkahm_smag = 1
+#endif
+
       !
    END SUBROUTINE ldf_dyn_init
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90	(revision 3625)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90	(revision 3634)
@@ -23,4 +23,10 @@
    REAL(wp), PUBLIC ::   rn_ahm_0_blp    =     0._wp   !: lateral bilaplacian eddy viscosity (m4/s)
    REAL(wp), PUBLIC ::   ahm0, ahmb0, ahm0_blp         !: OLD namelist names
+   REAL(wp), PUBLIC ::   rn_cmsmag_1     =     3._wp   !: constant in laplacian Smagorinsky viscosity
+   REAL(wp), PUBLIC ::   rn_cmsmag_2     =     3._wp   !: constant in bilaplacian Smagorinsky viscosity
+   REAL(wp), PUBLIC ::   rn_cmsh         =     1._wp   !: 1 or 0 , if 0 -use only shear for Smagorinsky viscosity
+   REAL(wp), PUBLIC ::   rn_ahm_m_blp    = -1.e12_wp   !: upper limit for bilap  abs(ahm) < min( dx^4/128rdt, rn_ahm_m_blp)
+   REAL(wp), PUBLIC ::   rn_ahm_m_lap    =  40000_wp   !: upper limit for lap  ahm < min(dx^2/16rdt, rn_ahm_m_lap)
+   INTEGER , PUBLIC ::   nkahm_smag      =  0          !: 
 
    !                                                                                  !!! eddy coeff. at U-,V-,W-pts [m2/s]
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_smag.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_smag.F90	(revision 3634)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_smag.F90	(revision 3634)
@@ -0,0 +1,301 @@
+MODULE ldfdyn_smag
+   !!======================================================================
+   !!                     ***  MODULE  ldftrasmag  ***
+   !! Ocean physics:  variable eddy induced velocity coefficients
+   !!======================================================================
+#if   defined key_dynldf_smag   &&   defined key_dynldf_c3d
+   !!----------------------------------------------------------------------
+   !!   'key_dynldf_smag'      and           smagorinsky  diffusivity
+   !!   'key_dynldf_c3d'                    3D tracer lateral  mixing coef.
+   !!----------------------------------------------------------------------
+   !!   ldf_eiv      : compute the eddy induced velocity coefficients
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE sbcrnf          ! river runoffs
+   USE ldfdyn_oce      ! ocean tracer   lateral physics
+   USE phycst          ! physical constants
+   USE ldfslp          ! iso-neutral slopes
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
+   USE iom
+   USE wrk_nemo
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC ldf_dyn_smag               ! routine called by step.F90
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Id: ldf_tra_smag.F90 1482 2010-06-13 15:28:06Z  $
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !!----------------------------------------------------------------------
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+
+
+
+
+   !!----------------------------------------------------------------------
+   !!                        ***  ldfdyn_smag.F90  ***
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Id: ldfdyn_c3d.h90 1581 2009-08-05 14:53:12Z smasson $ 
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   'key_dynldf_smag'             3D lateral eddy viscosity coefficients
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE ldf_dyn_smag( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ldf_dyn_smag  ***
+      !!                   
+      !! ** Purpose :   initializations of the horizontal ocean physics
+      !!
+      !! ** Method  :   3D eddy viscosity coef. 
+      !!    M.Griffies, R.Hallberg AMS, 2000
+      !! for laplacian:
+      !!   Asmag=(C/pi)^2*dx*dy sqrt(D^2), C=3-4
+      !! for bilaplacian:
+      !!   Bsmag=Asmag*dx*dy/8
+      !!   D^2=(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates
+      !!  in general case du/dx ==> e2 d(u/e2)/dx;  du/dy ==> e1 d(u/e1)/dy; 
+      !!                  dv/dx ==> e2 d(v/e2)/dx;  dv/dy ==> e1 d(v/e1)/dy
+      !!
+      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
+      !!                              ahm3, ahm4 never used
+      !!       bilaplacian operator : ahm1, ahm2 never used
+      !!                           :  ahm3, ahm4 defined at U- and V-points
+      !!       explanation of the default is missingi
+      !!  last modified : Maria Luneva, September 2011
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      !! ahm0 here is a background viscosity
+
+      !! * Arguments
+
+      !! * local variables
+
+      INTEGER              :: kt                   ! timestep
+
+      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
+      REAL (wp):: zdeltat,zdeltaf,zdeltau,zdeltav  ! temporary scalars
+      REAL (wp), POINTER, DIMENSION (:,:) ::   zux, zuy , zvx ,zvy, zue1, zue2, zve1, zve2 
+      REAL (wp)::  zcmsmag_1, zcmsmag_2 , zcmsh
+
+
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi,jpj,zux,zuy,zvx,zvy )
+      CALL wrk_alloc( jpi,jpj,zue1,zue2,zve1,zve2 )
+
+
+      IF(  kt == nit000 ) THEN
+
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'ldf_dyn_smag : 3D lateral eddy viscosity coefficient'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+     
+      ENDIF
+      
+      zcmsmag_1 = rn_cmsmag_1
+      zcmsmag_2 = rn_cmsmag_2
+      zcmsh     = rn_cmsh
+
+   
+
+
+      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operators
+      ! =================                       whatever its orientation is)
+      IF( ln_dynldf_lap ) THEN
+         ! define ahm1 and ahm2 at the right grid point position
+         ! (USER: modify ahm1 and ahm2 following your desiderata)
+         
+         DO jk=1,jpk
+           zue2(:,:)=un(:,:,jk)/e2u(:,:)
+           zve1(:,:)=vn(:,:,jk)/e1v(:,:)
+           zue1(:,:)=un(:,:,jk)/e1u(:,:)
+           zve2(:,:)=vn(:,:,jk)/e2v(:,:)
+
+ 
+           DO jj=2,jpj
+            DO ji=2,jpi
+            zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk) * zcmsh
+            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk) * zcmsh
+            ENDDO
+           ENDDO
+
+           DO jj=1,jpjm1
+            DO ji=1,jpim1
+            zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
+            zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
+            ENDDO
+           ENDDO
+             
+          DO jj=2,jpjm1
+           DO ji=2,jpim1
+
+            zdeltat=2._wp /(e1t(ji,jj)**(-2)+e2t(ji,jj)**(-2))
+            zdeltaf=2._wp /(e1f(ji,jj)**(-2)+e2f(ji,jj)**(-2))   
+            ahm1(ji,jj,jk)=(zcmsmag_1/rpi)**2*zdeltat*                                         &
+                            sqrt( (zux(ji,jj)-zvy(ji,jj))**2+                                    &
+                     0.0625_wp*(zuy(ji,jj)+zuy(ji,jj-1)+zuy(ji-1,jj)+zuy(ji-1,jj-1)+             &
+                            zvx(ji,jj)+zvx(ji,jj-1)+zvx(ji-1,jj)+zvx(ji-1,jj-1))**2)
+
+            ahm2(ji,jj,jk)=(zcmsmag_1/rpi)**2*zdeltaf*                                         &
+                            sqrt( (zuy(ji,jj)+zvx(ji,jj))**2+                                    &
+                     0.0625_wp*(zux(ji,jj)+zux(ji,jj+1)+zux(ji+1,jj)+zux(ji+1,jj+1)-             &
+                             zvy(ji,jj)-zvy(ji,jj+1)-zvy(ji+1,jj)-zvy(ji+1,jj+1))**2)
+
+            ahm1(ji,jj,jk)=MAX(ahm1(ji,jj,jk),rn_ahm_0_lap)
+            ahm2(ji,jj,jk)=MAX(ahm2(ji,jj,jk),rn_ahm_0_lap)
+
+! stability criteria or upper limit set from namelist
+            ahm1(ji,jj,jk)=MIN(ahm1(ji,jj,jk),zdeltat / (16_wp*rdt),rn_ahm_m_lap)
+            ahm2(ji,jj,jk)=MIN(ahm2(ji,jj,jk),zdeltaf / (16_wp*rdt),rn_ahm_m_lap)
+
+            ENDDO
+           ENDDO
+
+         ENDDO ! jpk
+            ahm1(:,:,jpk) = ahm1(:,:,jpkm1)
+            ahm2(:,:,jpk) = ahm2(:,:,jpkm1)
+            IF(lwp.and.kt==nit000) WRITE(numout,'(36x," ahm ", 7x)')
+            DO jk = 1, jpk
+
+               IF(lwp.and.kt==nit000) WRITE(numout,'(30x,E10.2,8x,i3)') ahm1(jpi/2,jpj/2,jk), jk
+            END DO
+      CALL lbc_lnk( ahm1, 'T', 1. )   ! Lateral boundary conditions on ( ahtt )
+      CALL lbc_lnk( ahm2, 'F', 1. )   ! Lateral boundary conditions on ( ahtt )
+
+      ENDIF    ! ln_dynldf_lap
+      
+
+
+      ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
+      ! ================================  whatever its orientation is)
+      ! Here: ahm is proportional to the cube of the maximum of the grid spacing
+      !       in the to horizontal direction
+
+      IF( ln_dynldf_bilap ) THEN
+         DO jk=1,jpk
+           zue2(:,:) = un(:,:,jk)/e2u(:,:)
+           zve1(:,:) = vn(:,:,jk)/e1v(:,:)
+           zue1(:,:) = un(:,:,jk)/e1u(:,:)
+           zve2(:,:) = vn(:,:,jk)/e2v(:,:)
+
+
+           DO jj=2,jpj
+            DO ji=2,jpi
+            zux(ji,jj) = (zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
+            zvy(ji,jj) = (zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
+            ENDDO
+           ENDDO
+
+           DO jj=1,jpjm1
+            DO ji=1,jpim1
+            zuy(ji,jj) = (zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
+            zvx(ji,jj) = (zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
+            ENDDO
+           ENDDO
+
+
+          DO jj=2,jpjm1
+           DO ji=2,jpim1
+            zdeltau = 2._wp/(e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2))
+            zdeltav = 2._wp/(e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2))
+
+            ahm3(ji,jj,jk) = -(zcmsmag_2/rpi)**2/8.0_wp*zdeltau**2*                          &
+
+                         sqrt(0.25_wp*(zux(ji,jj)+zux(ji+1,jj)-zvy(ji,jj)-zvy(ji+1,jj))**2+    &
+                              0.25_wp*(zuy(ji,jj)+zuy(ji,jj-1)+zvx(ji,jj)+zvx(ji,jj-1))**2)
+
+            ahm4(ji,jj,jk) = -(zcmsmag_2/rpi)**2/8.0_wp*zdeltav**2*                           &
+
+                         sqrt(0.25_wp*(zux(ji,jj)+zux(ji,jj+1)-zvy(ji,jj)-zvy(ji,jj+1))**2+    &
+                              0.25_wp*(zuy(ji,jj)+zuy(ji-1,jj)+zvx(ji-1,jj)+zvx(ji,jj))**2)
+
+            ahm3(ji,jj,jk) = MIN (rn_ahm_0_blp , ahm3(ji,jj,jk) )
+            ahm4(ji,jj,jk) = MIN (rn_ahm_0_blp , ahm4(ji,jj,jk) ) 
+
+! stability criteria or upper limit set in namelist
+
+            ahm3(ji,jj,jk) = MAX( ahm3(ji,jj,jk),-zdeltau**2/( 128._wp*rdt ),rn_ahm_m_blp )
+            ahm4(ji,jj,jk) = MAX( ahm4(ji,jj,jk),-zdeltav**2/( 128._wp*rdt ),rn_ahm_m_blp )
+            
+
+            ENDDO
+           ENDDO
+
+         ENDDO
+            ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
+            ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
+
+      DO jk = 1, jpk
+      IF(  kt == nit000 ) THEN
+
+               IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(jpi/2,jpj/2,jk), jk
+      ENDIF   
+      END DO
+      CALL lbc_lnk( ahm3, 'U', 1. )   ! Lateral boundary conditions
+      CALL lbc_lnk( ahm4, 'V', 1. )
+   ENDIF
+
+      CALL wrk_dealloc( jpi,jpj,zux,zuy,zvx,zvy )
+      CALL wrk_dealloc( jpi,jpj,zue1,zue2,zve1,zve2 )
+   !    zumax = MAXVAL( ABS( ahm3(:,:,:) ) )                ! slower than the following loop on NEC SX5
+      zdeltat = 0._wp
+   If(ln_dynldf_lap)THEN
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zdeltat = MAX(zdeltat,ABS(ahm1(ji,jj,jk)),ABS(ahm2(ji,jj,jk)) )
+          END DO
+        END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_max( zdeltat )                 ! max over the global domain
+      !
+      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,'dynlap:  abs(ahm) max: ', zdeltat
+    ENDIF
+    If(ln_dynldf_bilap)THEN
+     zdeltat = 0._wp
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zdeltat = MAX(zdeltat,ABS(ahm3(ji,jj,jk)),ABS(ahm3(ji,jj,jk)) )
+          END DO
+        END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_max( zdeltat )                 ! max over the global domain
+      !
+      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,'dyn_bilap abs(ahm) max: ', zdeltat
+      !
+   ENDIF
+      !
+
+END SUBROUTINE ldf_dyn_smag
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                                         Dummy module
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE ldf_dyn_smag( kt )       ! Empty routine
+      WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
+   END SUBROUTINE ldf_dyn_smag
+#endif
+
+   END MODULE ldfdyn_smag
+
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 3625)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 3634)
@@ -66,8 +66,9 @@
       NAMELIST/namtra_ldf/ ln_traldf_lap  , ln_traldf_bilap,                  &
          &                 ln_traldf_level, ln_traldf_hor  , ln_traldf_iso,   &
-         &                 ln_traldf_grif , ln_traldf_gdia,                   &
-         &                 ln_triad_iso   , ln_botmix_grif,                   &
+         &                 ln_traldf_grif , ln_traldf_gdia ,                  &
+         &                 ln_triad_iso   , ln_botmix_grif ,                  &
          &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0,       &
-         &                 rn_slpmax
+         &                 rn_slpmax      , rn_chsmag      ,    rn_smsh,      &
+         &                 rn_aht_m
       !!----------------------------------------------------------------------
 
@@ -153,8 +154,18 @@
       IF(lwp)WRITE(numout,*) '      constant eddy diffusivity coef.   ahtu = ahtv = ahtw = aht0 = ', aht0
       IF( lk_traldf_eiv ) THEN
-         IF(lwp)WRITE(numout,*)
          IF(lwp)WRITE(numout,*) '      constant eddy induced velocity coef.   aeiu = aeiv = aeiw = aeiv0 = ', aeiv0
+      
       ENDIF
 #endif
+
+#if defined key_traldf_smag && ! defined key_traldf_c3d
+        CALL ctl_stop( 'key_traldf_smag can only be used with key_traldf_c3d' )
+#endif
+#if defined key_traldf_smag
+        IF(lwp) WRITE(numout,*)' SMAGORINSKY DIFFUSION'
+        IF(lwp .AND. rn_smsh < 1)  WRITE(numout,*)' only  shear is used '
+        IF(lwp.and.ln_traldf_bilap) CALL ctl_stop(' SMAGORINSKY + BILAPLACIAN - UNSTABLE OR NON_CONSERVATIVE' )
+#endif
+
       !
    END SUBROUTINE ldf_tra_init
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c3d.h90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c3d.h90	(revision 3625)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c3d.h90	(revision 3634)
@@ -108,6 +108,4 @@
       CALL lbc_lnk( aeiv, 'V', 1. )
       CALL lbc_lnk( aeiw, 'W', 1. )
-# endif
-
       IF(lwp .AND. ld_print ) THEN
          WRITE(numout,*)
@@ -121,4 +119,5 @@
          CALL prihre(aeiw(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
       ENDIF
-      !
+
+# endif
    END SUBROUTINE ldf_tra_c3d
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90	(revision 3625)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90	(revision 3634)
@@ -30,6 +30,10 @@
    REAL(wp), PUBLIC ::   rn_aeiv_0       = 2000._wp  !: eddy induced velocity coefficient (m2/s)
    REAL(wp), PUBLIC ::   rn_slpmax       = 0.01_wp   !: slope limit
+   REAL(wp), PUBLIC ::   rn_chsmag       = 1._wp     !:  multiplicative factor in Smagorinsky diffusivity
+   REAL(wp), PUBLIC ::   rn_smsh         = 1._wp     !:  Smagorinsky diffusivity: = 0 - use only sheer
+   REAL(wp), PUBLIC ::   rn_aht_m        = 2000._wp  !:  upper limit or stability criteria for lateral eddy diffusivity (m2/s)
 
    REAL(wp), PUBLIC ::   aht0, ahtb0, aeiv0         !!: OLD namelist names
+
    LOGICAL , PUBLIC ::   ln_triad_iso    = .FALSE.   !: calculate triads twice
    LOGICAL , PUBLIC ::   ln_botmix_grif  = .FALSE.   !: mixing on bottom
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_smag.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_smag.F90	(revision 3634)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_smag.F90	(revision 3634)
@@ -0,0 +1,210 @@
+MODULE ldftra_smag
+   !!======================================================================
+   !!                     ***  MODULE  ldftrasmag  ***
+   !! Ocean physics:  variable eddy induced velocity coefficients
+   !!======================================================================
+#if   defined key_traldf_smag   &&   defined key_traldf_c3d
+   !!----------------------------------------------------------------------
+   !!   'key_traldf_smag'      and           smagorinsky  diffusivity
+   !!   'key_traldf_c3d'                    3D tracer lateral  mixing coef.
+   !!----------------------------------------------------------------------
+   !!   ldf_eiv      : compute the eddy induced velocity coefficients
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE sbcrnf          ! river runoffs
+   USE ldftra_oce      ! ocean tracer   lateral physics
+   USE phycst          ! physical constants
+   USE ldfslp          ! iso-neutral slopes
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
+   USE iom
+   USE wrk_nemo
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC ldf_tra_smag               ! routine called by step.F90
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Id: ldf_tra_smag.F90 1482 2010-06-13 15:28:06Z  $
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !!----------------------------------------------------------------------
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+ 
+
+
+
+  !!----------------------------------------------------------------------
+  !!                        ***  ldf_tra_smag.F90  ***
+  !!----------------------------------------------------------------------
+
+
+   SUBROUTINE ldf_tra_smag( kt )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ldf_tra_smag  ***
+      !!                   
+      !! ** Purpose :   initializations of the horizontal ocean physics
+      !!
+      !! ** Method  :   3D eddy viscosity coef. 
+      !!    M.Griffies, R.Hallberg AMS, 2000
+      !!  for laplacian:
+      !!   Asmag=(C/pi)^2*dx*dy sqrt(D^2), C=1 for tracers, C=3-4 for viscosity
+      !!   D^2= rm_smsh*(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates
+      !!  IF rm_smsh = 0 , only shear is used, recommended for tidal flows
+      !!  in general case du/dx ==> e2 d(u/e2)/dx;  du/dy ==> e1 d(u/e1)/dy; 
+      !!                  dv/dx ==> e2 d(v/e2)/dx;  dv/dy ==> e1 d(v/e1)/dy
+      !!  for bilaplacian: now this option is deleted as unstable or non-conservative
+      !!   - delta{Bsmag (delta(T)} = -Bsmag* delta{delta(T)} - delta(Bsmag)*delta( T )
+      !!  second term is of arbitrary sign on the edge of fronts and can induce instability 
+      !!   Bsmag=Asmag*dx*dy/8
+      !!  
+      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
+      !!                              ahm3, ahm4 never used
+      !!       bilaplacian operator : ahm1, ahm2 never used
+      !!                           :  ahm3, ahm4 defined at U- and V-points
+      !!       ??? explanation of the default is missing
+      !!  last modified : Maria Luneva, October 2012
+      !!----------------------------------------------------------------------
+      !!
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      USE ioipsl
+      REAL ( wp), POINTER , DIMENSION (:,:) :: zux, zvx , zuy , zvy 
+      REAL ( wp), POINTER , DIMENSION (:,:) :: zue1, zue2 , zve1 , zve2 
+      INTEGER, INTENT( in )                 ::   kt                             ! ocean time-step inedx
+      !! * Arguments
+      INTEGER                               :: ji,jj,jk
+
+      REAL (wp)                             ::     zdeltau, zdeltav, zhsmag ,zsmsh    ! temporary scalars
+      
+      CALL wrk_alloc (jpi,jpj,zux, zvx , zuy , zvy )
+      CALL wrk_alloc (jpi,jpj,zue1, zue2 , zve1 , zve2 )
+      !!----------------------------------------------------------------------
+      IF(  kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity '
+         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  '
+         IF(lwp) WRITE(numout,*) '               Coefficients are computed'
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      zhsmag = rn_chsmag
+      zsmsh  = rn_smsh 
+      zux(:,:)=0._wp ; zuy(:,:)=0._wp ; zvx(:,:)=0._wp ; zvy(:,:)=0._wp
+
+      ! -------------------
+      ahtt(:,:,:) = rn_aht_0
+       IF( ln_traldf_bilap ) THEN
+        IF( lwp .AND. kt == nit000) WRITE(numout,* )'ldf_tra_smag :no bilaplacian Smagorinsky diffusivity'
+        IF( lwp .AND. kt == nit000) WRITE(numout,* )'ldf_tra_smag :bilaplacian diffusivity set to constant'  
+       ENDIF
+
+
+
+      ! harmonic operator   (U-, V-, W-points)
+      ! ----------------- 
+
+      ahtu(:,:,:) = rn_aht_0                  ! set ahtu , ahtv at u- and v-points,
+      ahtv(:,:,:) = rn_aht_0                  ! and ahtw at w-point
+      ahtw(:,:,:) = rn_aht_0                  ! (here example: no space variation)
+      
+      IF( ln_traldf_lap ) THEN
+
+         DO jk=1,jpk
+
+           zue2(:,:)=un(:,:,jk)/e2u(:,:)
+           zve1(:,:)=vn(:,:,jk)/e1v(:,:)
+           zue1(:,:)=un(:,:,jk)/e1u(:,:)
+           zve2(:,:)=vn(:,:,jk)/e2v(:,:)
+
+
+           DO jj=2,jpj
+            DO ji=2,jpi
+            zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk) * zsmsh 
+            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk) * zsmsh 
+            ENDDO
+           ENDDO
+
+           DO jj=1,jpjm1
+            DO ji=1,jpim1
+            zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
+            zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
+            ENDDO
+           ENDDO
+
+
+          DO jj=2,jpjm1
+           DO ji=2,jpim1
+           zdeltau=2._wp/( e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2) )
+           zdeltav=2._wp/( e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2) )
+
+           ahtu(ji,jj,jk)=MAX( rn_aht_0 , (zhsmag/rpi)**2*zdeltau*                                &
+                          SQRT(0.25_wp*( zux(ji,jj)+zux(ji+1,jj)-zvy(ji,jj)-zvy(ji+1,jj) )**2+    &
+                               0.25_wp*( zuy(ji,jj)+zuy(ji,jj-1)+zvx(ji,jj)+zvx(ji,jj-1)  )**2) )
+
+           ahtv(ji,jj,jk)=MAX( rn_aht_0 ,  (zhsmag/rpi)**2*zdeltav*                               &
+                          SQRT(0.25_wp*( zux(ji,jj)+zux(ji,jj+1)-zvy(ji,jj)-zvy(ji,jj+1) )**2+    &
+                               0.25_wp*( zuy(ji,jj)+zuy(ji-1,jj)+zvx(ji-1,jj)+zvx(ji,jj)  )**2) )
+
+
+         !!! stability criteria: aht<delta**2/(4*dt)   dt=2*rdt , positiveness require aht<delta**2/(8*dt)
+             ahtu(ji,jj,jk)=MIN(ahtu(ji,jj,jk),zdeltau/(16*rdt) ,rn_aht_m)
+             ahtv(ji,jj,jk)=MIN(ahtv(ji,jj,jk),zdeltav/(16*rdt) ,rn_aht_m)
+         ! so...
+
+
+            ENDDO
+           ENDDO
+         ENDDO
+        ENDIF
+            ahtu(:,:,jpk) = ahtu(:,:,jpkm1)
+            ahtv(:,:,jpk) = ahtv(:,:,jpkm1)
+        CALL lbc_lnk( ahtu, 'U', 1. )   ! Lateral boundary conditions
+        CALL lbc_lnk( ahtv, 'V', 1. )
+
+IF(  kt == nit000 ) THEN
+
+      IF(lwp ) THEN                    ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'inildf: ahtu at k = 1'
+         CALL prihre( ahtu(:,:,1), jpi, jpj, 1, jpi, 1,   &
+            &                                1, jpj, 1, 1.e-1, numout )
+         WRITE(numout,*)
+         WRITE(numout,*) 'inildf: ahtv at k = 1'
+         CALL prihre( ahtv(:,:,1), jpi, jpj, 1, jpi, 1,   &
+            &                                1, jpj, 1, 1.e-1, numout )
+         WRITE(numout,*)
+         WRITE(numout,*) 'inildf: ahtw at k = 1'
+         CALL prihre( ahtw(:,:,1), jpi, jpj, 1, jpi, 1,   &
+            &                                1, jpj, 1, 1.e-1, numout )
+      ENDIF
+ENDIF
+
+      CALL wrk_dealloc ( jpi,jpj,zux, zvx , zuy , zvy     )
+      CALL wrk_dealloc ( jpi,jpj,zue1, zue2 , zve1 , zve2 )
+
+
+END SUBROUTINE ldf_tra_smag
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                                         Dummy module
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE ldf_tra_smag( kt )       ! Empty routine
+      WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
+   END SUBROUTINE ldf_tra_smag
+#endif
+
+END MODULE ldftra_smag
