Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/CONFIG/GYRE/cpp_GYRE.fcm
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/CONFIG/GYRE/cpp_GYRE.fcm	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/CONFIG/GYRE/cpp_GYRE.fcm	(revision 2926)
@@ -1,1 +1,1 @@
- bld::tool::fppkeys key_gyre key_dynspg_flt key_ldfslp key_zdftke   key_traldf_c3d key_traldf_smag key_dynldf_c3d key_dynldf_smag key_mpp_mpi
+ bld::tool::fppkeys key_gyre key_dynspg_flt key_ldfslp key_zdftke key_vectopt_loop key_iomput
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 2926)
@@ -4,4 +4,11 @@
    !! Ocean dynamics:  lateral viscosity trend
    !!======================================================================
+   !! History :  OPA  ! 1990-09  (G. Madec)  Original code
+   !!            4.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
+   !!            6.0  ! 1996-01  (G. Madec)  statement function for e3
+   !!            8.0  ! 1997-07  (G. Madec)  lbc calls
+   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
+   !!            2.0  ! 2004-08  (C. Talandier) New trends organization
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -9,5 +16,4 @@
    !!                   using an iso-level bilaplacian operator
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -21,6 +27,5 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dyn_ldf_bilap  ! called by step.F90
+   PUBLIC   dyn_ldf_bilap   ! called by step.F90
 
    !! * Substitutions
@@ -31,7 +36,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -69,32 +73,24 @@
       !! ** Action : - Update (ua,va) with the before iso-level biharmonic
       !!               mixing trend.
-      !!
-      !! History :
-      !!        !  90-09  (G. Madec)  Original code
-      !!        !  91-11  (G. Madec)
-      !!        !  93-03  (M. Guyon)  symetrical conditions (M. Guyon)
-      !!        !  96-01  (G. Madec)  statement function for e3
-      !!        !  97-07  (G. Madec)  lbc calls
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!   9.0  !  04-08  (C. Talandier) New trends organization
       !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk                ! dummy loop indices
-      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v ! temporary scalar
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zcu, zcv                             ! temporary workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         zuf, zut, zlu, zlv                   ! temporary workspace
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   zcu => wrk_2d_1 , zcv => wrk_2d_2   ! 3D workspace
+      USE wrk_nemo, ONLY:   zuf => wrk_3d_3 , zut => wrk_3d_4   ! 3D workspace
+      USE wrk_nemo, ONLY:   zlu => wrk_3d_5 , zlv => wrk_3d_6
+      !
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
+      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v   ! temporary scalar
       !!----------------------------------------------------------------------
-      !!  OPA 8.5, LODYC-IPSL (2002)
-      !!----------------------------------------------------------------------
-
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+
+      IF( wrk_in_use(2, 1,2) .OR. wrk_in_use(3, 3,4,5,6) ) THEN
+         CALL ctl_stop('dyn_ldf_bilap : requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
       ENDIF
 
@@ -104,8 +100,8 @@
 !!$      zlu(:,:,jpk) = 0.e0
 !!$      zlv(:,:,jpk) = 0.e0
-      zuf(:,:,:) = 0.e0
-      zut(:,:,:) = 0.e0
-      zlu(:,:,:) = 0.e0
-      zlv(:,:,:) = 0.e0
+      zuf(:,:,:) = 0._wp
+      zut(:,:,:) = 0._wp
+      zlu(:,:,:) = 0._wp
+      zlv(:,:,:) = 0._wp
 
       !                                                ! ===============
@@ -137,10 +133,7 @@
             END DO  
          ENDIF
-      ENDDO
-
-      ! Boundary conditions on the laplacian  (zlu,zlv)
-      CALL lbc_lnk( zlu, 'U', -1. )
-      CALL lbc_lnk( zlv, 'V', -1. )
-         
+      END DO
+      CALL lbc_lnk( zlu, 'U', -1. )   ;   CALL lbc_lnk( zlv, 'V', -1. )   ! Boundary conditions
+
          
       DO jk = 1, jpkm1
@@ -150,8 +143,7 @@
          
          ! Multiply by the eddy viscosity coef. (at u- and v-points)
-#if ! defined key_dynldf_smag
          zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk)
          zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk)
-#endif         
+         
          ! Contravariant "laplacian"
          zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
@@ -195,5 +187,5 @@
          ! Bilaplacian
          ! -----------
-#if !defined key_dynldf_smag
+
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -211,27 +203,11 @@
             END DO
          END DO
-#else
-
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
-               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
-               ! horizontal biharmonic diffusive trends
-               zua = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
-                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj)
-
-               zva = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
-                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)
-               ! add it to the general momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zua*fsahmu(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) + zva*fsahmv(ji,jj,jk)
-            END DO
-         END DO
-
-#endif
+
          !                                             ! ===============
       END DO                                           !   End of slab
       !                                                ! ===============
-
+      IF( wrk_not_released(2, 1,2)     .OR.   &
+          wrk_not_released(3, 3,4,5,6) )   CALL ctl_stop('dyn_ldf_bilap: failed to release workspace arrays')
+      !
    END SUBROUTINE dyn_ldf_bilap
 
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 2926)
@@ -4,4 +4,8 @@
    !! Ocean dynamics:  lateral viscosity trend
    !!======================================================================
+   !! History :  OPA  !  1997-07  (G. Madec)  Original code
+   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!            2.0  !  2004-08  (C. Talandier) New trends organization
+   !!----------------------------------------------------------------------
 #if defined key_ldfslp   ||   defined key_esopa
    !!----------------------------------------------------------------------
@@ -12,13 +16,13 @@
    !!   ldfguv         : 
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
    USE ldfdyn_oce      ! ocean dynamics lateral physics
    USE zdf_oce         ! ocean vertical physics
-   USE in_out_manager  ! I/O manager
    USE trdmod          ! ocean dynamics trends 
    USE trdmod_oce      ! ocean variables trends
    USE ldfslp          ! iso-neutral slopes available
+   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
@@ -27,6 +31,8 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dyn_ldf_bilapg ! called by step.F90
+   PUBLIC   dyn_ldf_bilapg       ! called by step.F90
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zfuw, zfvw , zdiu, zdiv   ! 2D workspace (ldfguv)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zdju, zdj1u, zdjv, zdj1v  ! 2D workspace (ldfguv)
 
    !! * Substitutions
@@ -36,8 +42,18 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
+
+   INTEGER FUNCTION dyn_ldf_bilapg_alloc()
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE dyn_ldf_bilapg_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) ,     &
+         &      zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc )
+         !
+      IF( dyn_ldf_bilapg_alloc /= 0 )   CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays')
+   END FUNCTION dyn_ldf_bilapg_alloc
+
 
    SUBROUTINE dyn_ldf_bilapg( kt )
@@ -67,23 +83,17 @@
       !!                biharmonic mixing trend.
       !!              - save the trend in (zwk3,zwk4) ('key_trddyn')
-      !!
-      !! History :
-      !!   8.0  !  97-07  (G. Madec)  Original code
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!   9.0  !  04-08  (C. Talandier) New trends organization
-      !!----------------------------------------------------------------------
-      !! * Modules used     
-      USE oce, ONLY :    zwk3 => ta,   & ! use ta as 3D workspace   
-                         zwk4 => sa      ! use sa as 3D workspace   
-
-      !! * Arguments
+      !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   zwk1 => wrk_3d_3 , zwk2 => wrk_3d_4   ! 3D workspace
+      USE oce     , ONLY:   zwk3 => ta       , zwk4 => sa         ! ta, sa used as 3D workspace   
+      !
       INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
-
-      !! * Local declarations
+      !
       INTEGER ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         zwk1, zwk2                ! work array used for rotated biharmonic
-         !                         ! operator on tracers and/or momentum
-      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(3, 3,4) ) THEN
+         CALL ctl_stop('dyn_ldf_bilapg: requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
 
       IF( kt == nit000 ) THEN
@@ -93,42 +103,32 @@
          zwk1(:,:,:) = 0.e0   ;   zwk3(:,:,:) = 0.e0
          zwk2(:,:,:) = 0.e0   ;   zwk4(:,:,:) = 0.e0
+         !                                      ! allocate dyn_ldf_bilapg arrays
+         IF( dyn_ldf_bilapg_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays')
       ENDIF
 
       ! Laplacian of (ub,vb) multiplied by ahm
       ! --------------------------------------  
-      ! rotated harmonic operator applied to (ub,vb)
-      !     and multiply by ahmu, ahmv (output in (zwk1,zwk2) )
-
-      CALL ldfguv ( ub, vb, zwk1, zwk2, 1 )
-
-
-      ! Lateral boundary conditions on (zwk1,zwk2)
-      CALL lbc_lnk( zwk1, 'U', -1. )
-      CALL lbc_lnk( zwk2, 'V', -1. )
-
+      CALL ldfguv( ub, vb, zwk1, zwk2, 1 )      ! rotated harmonic operator applied to (ub,vb)
+      !                                         ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) )
+      CALL lbc_lnk( zwk1, 'U', -1. )   ;   CALL lbc_lnk( zwk2, 'V', -1. )     ! Lateral boundary conditions
 
       ! Bilaplacian of (ub,vb)
       ! ---------------------- 
-      ! rotated harmonic operator applied to (zwk1,zwk2) (output in (zwk3,zwk4) )
-
-      CALL ldfguv ( zwk1, zwk2, zwk3, zwk4, 2 )
-
-
-      ! Update the momentum trends			   (j-slab :   2, jpj-1)
+      CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 )  ! rotated harmonic operator applied to (zwk1,zwk2) 
+      !                                         ! (output in (zwk3,zwk4) )
+
+      ! Update the momentum trends
       ! --------------------------
-      !                                                ! ===============
-      DO jj = 2, jpjm1                                 !  Vertical slab
-         !                                             ! ===============
+      DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
-               ! add the diffusive trend to the general momentum trends
                ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk)
                va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk)
             END DO
          END DO
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-
+      END DO
+      !
+      IF( wrk_not_released(3, 3,4) )   CALL ctl_stop('dyn_ldf_bilapg: failed to release workspace arrays')
+      !
    END SUBROUTINE dyn_ldf_bilapg
 
@@ -174,37 +174,27 @@
       !!                          second order vertical derivative term)
       !!      'key_trddyn' defined: the trend is saved for diagnostics.
-      !!
-      !! History :
-      !!   8.0  !  97-07  (G. Madec)  Original code
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   &
-         pu, pv     ! momentum fields (before u and v for the 1st call, and
-      !             ! laplacian of these fields multiplied by ahm for the 2nd
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   &
-         plu, plv   ! partial harmonic operator applied to
-      !             ! pu and pv (all the components except
-      !             ! second order vertical derivative term)
-      INTEGER, INTENT( in ) ::   &
-         kahm       ! =1 the laplacian is multiplied by the eddy diffusivity coef.
-      !             ! =2 no multiplication
-
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk       ! dummy loop indices
-      REAL(wp) ::   &
-         zabe1, zabe2, zcof1, zcof2,    &  ! temporary scalars
-         zcoef0, zcoef3, zcoef4
-      REAL(wp) ::   &
-         zbur, zbvr, zmkt, zmkf, zuav, zvav,    &
-         zuwslpi, zuwslpj, zvwslpi, zvwslpj
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         ziut, zjuf , zjvt, zivf,       &  ! workspace
-         zdku, zdk1u, zdkv, zdk1v
-      REAL(wp), DIMENSION(jpi,jpk) ::   &
-         zfuw, zfvw, zdiu, zdiv,        &  ! workspace
-         zdju, zdj1u, zdjv, zdj1v 
-      !!----------------------------------------------------------------------
-
+      !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   ziut => wrk_2d_1 , zjuf  => wrk_2d_2 , zjvt  => wrk_2d_3
+      USE wrk_nemo, ONLY:   zivf => wrk_2d_4 , zdku  => wrk_2d_5 , zdk1u => wrk_2d_6
+      USE wrk_nemo, ONLY:   zdkv => wrk_2d_7 , zdk1v => wrk_2d_8
+      !!
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu , pv    ! 1st call: before horizontal velocity 
+      !                                                               ! 2nd call: ahm x these fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   plu, plv   ! partial harmonic operator applied to
+      !                                                               ! pu and pv (all the components except
+      !                                                               ! second order vertical derivative term)
+      INTEGER                         , INTENT(in   ) ::   kahm       ! =1 1st call ; =2 2nd call
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zabe1 , zabe2 , zcof1 , zcof2        ! local scalar
+      REAL(wp) ::   zcoef0, zcoef3, zcoef4               !   -      -
+      REAL(wp) ::   zbur, zbvr, zmkt, zmkf, zuav, zvav   !   -      -
+      REAL(wp) ::   zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN
+         CALL ctl_stop('dyn:ldfguv: requested workspace arrays unavailable')   ;   RETURN
+      END IF
       !                               ! ********** !   ! ===============
       DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
@@ -422,5 +412,5 @@
          ! II.3 Divergence of vertical fluxes added to the horizontal divergence
          ! ---------------------------------------------------------------------
-#if defined key_dynldf_smag 
+
          IF( kahm == 1 ) THEN
             ! multiply the laplacian by the eddy viscosity coefficient
@@ -458,45 +448,10 @@
             STOP 'ldfguv'
          ENDIF
-#else
-         IF( kahm == 2 ) THEN
-            ! multiply the laplacian by the eddy viscosity coefficient
-            DO jk = 1, jpkm1
-               DO ji = 2, jpim1
-                  ! eddy coef. divided by the volume element
-                  zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
-                  zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
-                  ! vertical divergence
-                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
-                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
-                  ! harmonic operator applied to (pu,pv) and multiply by ahm
-                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
-                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
-               END DO
-            END DO
-         ELSEIF( kahm == 1 ) THEN
-            ! second call, no multiplication
-            DO jk = 1, jpkm1
-               DO ji = 2, jpim1
-                  ! inverse of the volume element
-                  zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
-                  zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
-                  ! vertical divergence
-                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
-                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
-                  ! harmonic operator applied to (pu,pv) 
-                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
-                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
-               END DO
-            END DO
-         ELSE
-            IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm
-            IF(lwp)WRITE(numout,*) '         We stop'
-            STOP 'ldfguv'
-         ENDIF
-
-#endif
          !                                             ! ===============
       END DO                                           !   End of slab
       !                                                ! ===============
+
+      IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn:ldfguv: failed to release workspace arrays')
+      !
    END SUBROUTINE ldfguv
 
@@ -507,4 +462,5 @@
 CONTAINS
    SUBROUTINE dyn_ldf_bilapg( kt )               ! Dummy routine
+      INTEGER, INTENT(in) :: kt
       WRITE(*,*) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt
    END SUBROUTINE dyn_ldf_bilapg
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 2926)
@@ -38,7 +38,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -63,9 +62,9 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ioptio         ! ???
-      LOGICAL :: ll_print = .FALSE.    ! Logical flag for printing viscosity coef.
+      LOGICAL ::   ll_print = .FALSE.    ! Logical flag for printing viscosity coef.
       !! 
       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 ,cmsmag1,cmsmag2
+         &                 rn_ahm_0_lap   , rn_ahmb_0      , rn_ahm_0_blp
       !!----------------------------------------------------------------------
 
@@ -86,6 +85,4 @@
          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,*) '      smagorinsky coefficient for laplacian       cmsmag1    = ' , cmsmag1
-         WRITE(numout,*) '      smagorinsky coefficient for bilaplacian     cmsmag2    = ' , cmsmag2
       ENDIF
 
@@ -131,8 +128,4 @@
       ! Lateral eddy viscosity
       ! ======================
-#if defined key_dynldf_smag && ! defined key_dynldf_c3d
-      IF(lwp) WRITE(numout,*) '   key_dynldf_smag can not be used without key_dynldf_c3d'
-#endif
-
 #if defined key_dynldf_c3d
       CALL ldf_dyn_c3d( ll_print )   ! ahm = 3D coef. = F( longitude, latitude, depth )
@@ -213,6 +206,6 @@
       REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
       REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
-      REAL(wp), INTENT(in   ), DIMENSION        (jpk) ::   pdep       ! depth of the gridpoint (T, U, V, F)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pah        ! adimensional vertical profile
+      REAL(wp), INTENT(in   ), DIMENSION          (:) ::   pdep       ! depth of the gridpoint (T, U, V, F)
+      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
       !!
       INTEGER  ::   jk           ! dummy loop indices
@@ -255,6 +248,6 @@
       REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
       REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pdep       ! dep of the gridpoint (T, U, V, F)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pah        ! adimensional vertical profile
+      REAL(wp), INTENT(in   ), DIMENSION      (:,:,:) ::   pdep       ! dep of the gridpoint (T, U, V, F)
+      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
       !!
       INTEGER  ::   jk           ! dummy loop indices
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90	(revision 2926)
@@ -23,6 +23,5 @@
    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 ::   cmsmag1          = 3._wp       !: constant in laplacian Smagorinsky viscosity
-   REAL(wp), PUBLIC ::   cmsmag2          = 3._wp       !: constant in bilaplacian Smagorinsky viscosity
+
    !                                                                                  !!! eddy coeff. at U-,V-,W-pts [m2/s]
 #if defined key_dynldf_c3d
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_smag.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_smag.F90	(revision 2925)
+++ 	(revision )
@@ -1,262 +1,0 @@
-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
-
-   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
-      USE ldftra_oce, ONLY : aht0
-
-      REAL (wp), DIMENSION (:,:), ALLOCATABLE::    ux,uy,vx,vy     ! local variables
-      REAL (wp), DIMENSION (:,:), ALLOCATABLE::    ue1,ue2,ve1,ve2 ! local variables
-      !! * Arguments
-      INTEGER              :: kt         ! timestep
-
-      !! * local variables
-      INTEGER  ::   ji, jj, jk      ! dummy loop indices
-      REAL (wp):: deltat,deltaf,deltau,deltav 
-      
-
-      !!----------------------------------------------------------------------
-      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
-   
-      ALLOCATE ( ux(jpi,jpj) ); ux(:,:)=0_wp
-      ALLOCATE ( uy(jpi,jpj) ); uy(:,:)=0_wp
-      ALLOCATE ( vx(jpi,jpj) ); vx(:,:)=0_wp
-      ALLOCATE ( vy(jpi,jpj) ); vy(:,:)=0_wp 
-      ALLOCATE(ue1(jpi,jpj)); ALLOCATE(ue2(jpi,jpj))
-      ALLOCATE(ve1(jpi,jpj)); ALLOCATE(ve2(jpi,jpj))
-
-      
-      ! 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
-           ue2(:,:)=un(:,:,jk)/e2u(:,:)
-           ve1(:,:)=vn(:,:,jk)/e1v(:,:)
-           ue1(:,:)=un(:,:,jk)/e1u(:,:)
-           ve2(:,:)=vn(:,:,jk)/e2v(:,:)
-
- 
-           DO jj=2,jpj
-            DO ji=2,jpi
-            ux(ji,jj)=(ue2(ji,jj)-ue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
-            vy(ji,jj)=(ve1(ji,jj)-ve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
-            ENDDO
-           ENDDO
-
-           DO jj=1,jpjm1
-            DO ji=1,jpi
-            uy(ji,jj)=(ue1(ji,jj+1)-ue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
-            vx(ji,jj)=(ve2(ji+1,jj)-ve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
-            ENDDO
-           ENDDO
-             
-          DO jj=2,jpjm1
-           DO ji=2,jpim1
-            deltat=2./(e1t(ji,jj)**(-2)+e2t(ji,jj)**(-2))
-            deltaf=2./(e1f(ji,jj)**(-2)+e2f(ji,jj)**(-2))   
-            ahm1(ji,jj,jk)=(cmsmag1/3.14)**2*deltat*                                       &
-                            sqrt( (ux(ji,jj)-vy(ji,jj))**2+                               &
-                     0.0625*(uy(ji,jj)+uy(ji,jj-1)+uy(ji-1,jj)+uy(ji-1,jj-1)+             &
-                            vx(ji,jj)+vx(ji,jj-1)+vx(ji-1,jj)+vx(ji-1,jj-1))**2)
-
-            ahm2(ji,jj,jk)=(cmsmag1/3.14)**2*deltaf*                                       &
-                            sqrt( (uy(ji,jj)+vx(ji,jj))**2+                               &
-                     0.0625*(ux(ji,jj)+ux(ji,jj+1)+ux(ji+1,jj)+ux(ji+1,jj+1)-             &
-                             vy(ji,jj)-vy(ji,jj+1)-vy(ji+1,jj)-vy(ji+1,jj+1))**2)
-
-            ahm1(ji,jj,jk)=MAX(ahm1(ji,jj,jk),ahm0)
-            ahm2(ji,jj,jk)=MAX(ahm2(ji,jj,jk),ahm0)
-
-! stability criteria 
-            ahm1(ji,jj,jk)=MIN(ahm1(ji,jj,jk),deltat**2/(16*rdt))
-            ahm2(ji,jj,jk)=MIN(ahm2(ji,jj,jk),deltaf**2/(16*rdt))
-
-            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
-      
-
-
-      ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
-      ! ================================  whatever its orientation is)
-      ! (USER: modify ahm3 and ahm4 following your desiderata)
-      ! Here: ahm is proportional to the cube of the maximum of the gridspacing
-      !       in the to horizontal direction
-
-      IF( ln_dynldf_bilap ) THEN
-         DO jk=1,jpk
-           ue2(:,:)=un(:,:,jk)/e2u(:,:)
-           ve1(:,:)=vn(:,:,jk)/e1v(:,:)
-           ue1(:,:)=un(:,:,jk)/e1u(:,:)
-           ve2(:,:)=vn(:,:,jk)/e2v(:,:)
-
-
-           DO jj=2,jpj
-            DO ji=2,jpi
-            ux(ji,jj)=(ue2(ji,jj)-ue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
-            vy(ji,jj)=(ve1(ji,jj)-ve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
-            ENDDO
-           ENDDO
-
-           DO jj=1,jpjm1
-            DO ji=1,jpim1
-            uy(ji,jj)=(ue1(ji,jj+1)-ue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
-            vx(ji,jj)=(ve2(ji+1,jj)-ve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
-            ENDDO
-           ENDDO
-
-
-          DO jj=2,jpjm1
-           DO ji=2,jpim1
-            deltau=2./(e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2))
-            deltav=2./(e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2))
-
-             ahm3(ji,jj,jk)=MIN(ahm0_blp,(cmsmag2/3.14)**2/8*deltau**2*            &
-
-                         sqrt(0.25*(ux(ji,jj)+ux(ji+1,jj)-vy(ji,jj)-vy(ji+1,jj))**2+    &
-                              0.25*(uy(ji,jj)+uy(ji,jj-1)+vx(ji,jj)+vx(ji,jj-1))**2))
-
-            ahm4(ji,jj,jk)=MIN(ahm0_blp ,(cmsmag2/3.14)**2/8*deltav**2*            &
-
-                         sqrt(0.25*(ux(ji,jj)+ux(ji,jj+1)-vy(ji,jj)-vy(ji,jj+1))**2+    &
-                              0.25*(uy(ji,jj)+uy(ji-1,jj)+vx(ji-1,jj)+vx(ji,jj))**2))
-! stability criteria 
-            ahm3(ji,jj,jk)=MAX(ahm3(ji,jj,jk),-deltau**2/(16*rdt))
-            ahm4(ji,jj,jk)=MAX(ahm4(ji,jj,jk),-deltav**2/(16*rdt))
-            
-
-            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
-     DEALLOCATE ( ux ); DEALLOCATE ( uy ); DEALLOCATE ( vx ); DEALLOCATE ( vy )
-     DEALLOCATE ( ue1 ); DEALLOCATE ( ue2 ); DEALLOCATE ( ve1 ); DEALLOCATE ( ve2 )
-
-   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/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2926)
@@ -31,5 +31,4 @@
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
-   USE lib_mpp
 
    IMPLICIT NONE
@@ -692,9 +691,9 @@
             !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
             !                                kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-            zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,ik)* ABS( zau )  )
-            zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ik)* ABS( zav )  )
+            zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau )  )
+            zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav )  )
             !                          !- Slope at u- & v-points (uslpml, vslpml)
-            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,ik)
-            vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ik)
+            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku)
+            vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv)
             !
             !                    !==   i- & j-slopes at w-points just below the Mixed Layer   ==!
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 2926)
@@ -68,5 +68,5 @@
          &                 ln_traldf_grif , ln_traldf_gdia,                   &
          &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0,       &
-         &                 rn_slpmax , chsmag
+         &                 rn_slpmax
       !!----------------------------------------------------------------------
 
@@ -157,8 +157,4 @@
       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
-
       !
    END SUBROUTINE ldf_tra_init
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c3d.h90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c3d.h90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c3d.h90	(revision 2926)
@@ -89,5 +89,5 @@
       ENDIF
 
-# if defined key_traldf_eiv 
+# if defined key_traldf_eiv
       aeiu(:,:,:) = aeiv0            ! set aeiu = aeiv at u- and v-points,
       aeiv(:,:,:) = aeiv0            ! and aeiw at w-point
@@ -108,4 +108,5 @@
       CALL lbc_lnk( aeiv, 'V', 1. )
       CALL lbc_lnk( aeiw, 'W', 1. )
+# endif
 
       IF(lwp .AND. ld_print ) THEN
@@ -120,5 +121,4 @@
          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/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90	(revision 2926)
@@ -30,5 +30,5 @@
    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 ::   chsmag          = 1._wp     !: constant in Smagorinsky diffusivity
+
    REAL(wp), PUBLIC ::   aht0, ahtb0, aeiv0         !!: OLD namelist names
    LOGICAL , PUBLIC ::   l_triad_iso     = .FALSE.   !: calculate triads twice
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_smag.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_smag.F90	(revision 2925)
+++ 	(revision )
@@ -1,263 +1,0 @@
-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
-
-   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
-      !! 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 missing
-      !!  last modified : Maria Luneva, September 2011
-      !!----------------------------------------------------------------------
-      !!
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE ioipsl
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx
-      !! * Arguments
-      INTEGER               :: ji,jj,jk
-      REAL (wp), DIMENSION (:,:), ALLOCATABLE::    ux,uy,vx,vy     ! local velocity derivatives
-      REAL (wp), DIMENSION (:,:), ALLOCATABLE::    ue1,ue2,ve1,ve2 ! local variables
-      REAL (wp)                     ::    deltaxy, deltau, deltav  ! local step 
-
-      !!----------------------------------------------------------------------
-      IF(  kt == nit000 ) THEN
-      IF( lk_traldf_eiv ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity and eddy induced velocity coefficients'
-         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  '
-         IF(lwp) WRITE(numout,*) '               Coefficients are computed'
-         IF(lwp) WRITE(numout,*)
-      ELSE
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity coefficient'
-         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   -- '
-         IF(lwp) WRITE(numout,*) '               Coefficients are computed'
-         IF(lwp) WRITE(numout,*)
-      ENDIF
-      ENDIF
-
-      ALLOCATE(ux(jpi,jpj)) ; ux(:,:)=0._wp
-      ALLOCATE(uy(jpi,jpj)) ; uy(:,:)=0._wp
-      ALLOCATE(vx(jpi,jpj)) ; vx(:,:)=0._wp
-      ALLOCATE(vy(jpi,jpj)) ; vy(:,:)=0._wp
-      ALLOCATE(ue1(jpi,jpj)); ALLOCATE(ue2(jpi,jpj)) 
-      ALLOCATE(ve1(jpi,jpj)); ALLOCATE(ve2(jpi,jpj)) 
-
-      ! biharmonic operator   (T-point)
-      ! -------------------
-     
-      ahtt(:,:,:) = aht0              ! set ahtt at T-point (here no space variation)
-        IF( ln_traldf_bilap ) THEN
-         ! define ahm1 and ahm2 at the right grid point position
-         ! (USER: modify ahm1 and ahm2 following your desiderata)
-         DO jk=1,jpk
-           ue2(:,:)=un(:,:,jk)/e2u(:,:)
-           ve1(:,:)=vn(:,:,jk)/e1v(:,:)
-           ue1(:,:)=un(:,:,jk)/e1u(:,:)
-           ve2(:,:)=vn(:,:,jk)/e2v(:,:)
-
-
-           DO jj=2,jpj
-            DO ji=2,jpi
-            ux(ji,jj)=(ue2(ji,jj)-ue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
-            vy(ji,jj)=(ve1(ji,jj)-ve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
-            ENDDO
-           ENDDO
-
-           DO jj=1,jpjm1
-            DO ji=1,jpim1
-            uy(ji,jj)=(ue1(ji,jj+1)-ue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
-            vx(ji,jj)=(ve2(ji+1,jj)-ve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
-            ENDDO
-           ENDDO
-
-
-          DO jj=2,jpjm1
-           DO ji=2,jpim1
-            deltaxy=2./(e1t(ji,jj)**(-2)+e2t(ji,jj)**(-2))
-             ahtt(ji,jj,jk)=(chsmag/3.14)**2/8*deltaxy**2*                      &
-                            sqrt( (ux(ji,jj)-vy(ji,jj))**2+                              &
-                     0.0625*(uy(ji,jj)+uy(ji,jj-1)+uy(ji-1,jj)+uy(ji-1,jj-1)+            &
-                        vx(ji,jj)+vx(ji,jj-1)+vx(ji-1,jj)+vx(ji-1,jj-1) )**2)
-   
-              !!! stabulity criteria: abs(aht)<delta**4/(4*8*dt)   dt=2*rdt 
-             ahtt(ji,jj,jk)=MIN( ahtt(ji,jj,jk) , aht0 )      
-             ahtt(ji,jj,jk)=MAX( ahtt(ji,jj,jk) ,- deltaxy**4/(8*16*rdt) ) 
-            ENDDO
-           ENDDO
-         ENDDO
-            ahtt(:,:,jpk) = ahtt(:,:,jpkm1)
-
-      CALL lbc_lnk( ahtt, 'T', 1. )   ! Lateral boundary conditions on ( ahtt )
-IF(  kt == nit000 ) THEN
-
-      IF(lwp ) THEN                    ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'inildf: ahtt at k = 1'
-         CALL prihre( ahtt(:,:,1), jpi, jpj, 1, jpi, 1,   &
-            &                                1, jpj, 1, 1.e-1, numout )
-
-      ENDIF
-
-      ENDIF
-ENDIF
-
-
-
-      ! harmonic operator   (U-, V-, W-points)
-      ! ----------------- 
-
-      ahtu(:,:,:) = aht0              ! set ahtu = ahtv at u- and v-points,
-      ahtv(:,:,:) = aht0              ! and ahtw at w-point
-      ahtw(:,:,:) = aht0              ! (here example: no space variation)
-
-      IF( ln_traldf_lap ) THEN
-         DO jk=1,jpk
-           ue2(:,:)=un(:,:,jk)/e2u(:,:)
-           ve1(:,:)=vn(:,:,jk)/e1v(:,:)
-           ue1(:,:)=un(:,:,jk)/e1u(:,:)
-           ve2(:,:)=vn(:,:,jk)/e2v(:,:)
-
-
-           DO jj=2,jpj
-            DO ji=2,jpi
-            ux(ji,jj)=(ue2(ji,jj)-ue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
-            vy(ji,jj)=(ve1(ji,jj)-ve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
-            ENDDO
-           ENDDO
-
-           DO jj=1,jpjm1
-            DO ji=1,jpim1
-            uy(ji,jj)=(ue1(ji,jj+1)-ue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
-            vx(ji,jj)=(ve2(ji+1,jj)-ve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
-            ENDDO
-           ENDDO
-
-
-          DO jj=2,jpjm1
-           DO ji=2,jpim1
-           deltau=2./(e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2))
-           deltav=2./(e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2))
-
-           ahtu(ji,jj,jk)=MAX(aht0 , (chsmag/3.14)**2*deltau*            &
-                         sqrt(0.25*(ux(ji,jj)+ux(ji+1,jj)-vy(ji,jj)-vy(ji+1,jj))**2+    &
-                              0.25*(uy(ji,jj)+uy(ji,jj-1)+vx(ji,jj)+vx(ji,jj-1))**2))
-
-           ahtv(ji,jj,jk)=MAX(aht0 , (chsmag/3.14)**2*deltav*            &
-                         sqrt(0.25*(ux(ji,jj)+ux(ji,jj+1)-vy(ji,jj)-vy(ji,jj+1))**2+    &
-                              0.25*(uy(ji,jj)+uy(ji-1,jj)+vx(ji-1,jj)+vx(ji,jj))**2))
-
-
-         !!! stabulity criteria: aht<delta**2/(4*dt)   dt=2*rdt
-
-             ahtu(ji,jj,jk)=MIN(ahtu(ji,jj,jk),deltau**2/(16*rdt) )
-             ahtv(ji,jj,jk)=MIN(ahtv(ji,jj,jk),deltav**2/(16*rdt) )
-
-            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. )
-
-
-          DEALLOCATE( ux ) ;DEALLOCATE( uy );DEALLOCATE( vx );DEALLOCATE( vy )     
-          DEALLOCATE( ue1 ) ;DEALLOCATE( ue2 );DEALLOCATE( ve1 );DEALLOCATE( ve2 )
-
-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
-
-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
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 2926)
@@ -28,4 +28,5 @@
    USE diaptr          ! poleward transport diagnostics
    USE trc_oce         ! share passive tracers/Ocean variables
+   USE lib_mpp         ! MPP library
 
    IMPLICIT NONE
@@ -73,6 +74,7 @@
       !!               biharmonic mixing trend.
       !!----------------------------------------------------------------------
-      USE oce         , ztu => ua   ! use ua as workspace
-      USE oce         , ztv => va   ! use va as workspace
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE oce     , ONLY:   ztu  => ua       , ztv  => va                           ! (ua,va) used as workspace
+      USE wrk_nemo, ONLY:   zeeu => wrk_2d_1 , zeev => wrk_2d_2 , zlt => wrk_2d_3   ! 2D workspace
       !!
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
@@ -85,6 +87,9 @@
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
       REAL(wp) ::  zbtr, ztra       ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zeeu, zeev, zlt   ! 2D workspace
       !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 1,2,3) ) THEN
+         CALL ctl_stop('tra_ldf_bilap: requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
 
       IF( kt == nit000 )  THEN
@@ -123,5 +128,4 @@
                END DO
             ENDIF
-#if ! defined key_traldf_smag
             DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -131,13 +135,4 @@
                END DO
             END DO
-#else
-            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  zlt(ji,jj) = zbtr * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
-                     &                                     + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   )
-               END DO
-            END DO
-#endif
             CALL lbc_lnk( zlt, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn)
 
@@ -150,5 +145,4 @@
                END DO
             END DO
-#if ! defined key_traldf_smag
             DO jj = 2, jpjm1                 ! fourth derivative (divergence) and add to the general tracer trend
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -160,16 +154,4 @@
                END DO
             END DO
-#else
-            DO jj = 2, jpjm1                 ! fourth derivative (divergence) and add to the general tracer trend
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! horizontal diffusive trends
-                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  ztra = zbtr * fsahtt(ji,jj,jk)* (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
-                  ! add it to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
-               END DO
-            END DO
-
-#endif
             !                                             
          END DO                                           ! Horizontal slab
@@ -183,4 +165,6 @@
       END DO                                              ! tracer loop
       !                                                   ! ===========
+      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('tra_ldf_bilap: failed to release workspace arrays')
+      !
    END SUBROUTINE tra_ldf_bilap
 
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 2926)
@@ -24,4 +24,5 @@
    USE diaptr          ! poleward transport diagnostics 
    USE trc_oce         ! share passive tracers/Ocean variables
+   USE lib_mpp         ! MPP library
 
    IMPLICIT NONE
@@ -65,4 +66,7 @@
       !!               biharmonic mixing trend.
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   wk1 => wrk_4d_1 , wk2 => wrk_4d_2     ! 4D workspace
+      !
       INTEGER         , INTENT(in   )                      ::   kt       ! ocean time-step index
       CHARACTER(len=3), INTENT(in   )                      ::   cdtype   ! =TRA or TRC (tracer indicator)
@@ -70,8 +74,11 @@
       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 
-      !!
-      INTEGER ::   ji, jj, jk, jn                 ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt) ::   wk1, wk2   ! 4D workspace
-      !!----------------------------------------------------------------------
+      !
+      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(4, 1,2) ) THEN
+         CALL ctl_stop('tra_ldf_bilapg: requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
 
       IF( kt == nit000 )  THEN
@@ -107,4 +114,6 @@
          END DO
       END DO
+      !
+      IF( wrk_not_released(4, 1,2) )   CALL ctl_stop('tra_ldf_bilapg : failed to release workspace arrays.')
       !
    END SUBROUTINE tra_ldf_bilapg
@@ -149,6 +158,10 @@
       !!
       !!----------------------------------------------------------------------
-      USE oce         , zftv => ua     ! use ua as workspace
-      !!
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released, wrk_in_use_xz, wrk_not_released_xz
+      USE oce     , ONLY:   zftv => ua       ! ua used as workspace
+      USE wrk_nemo, ONLY:   zftu => wrk_2d_1 , zdkt  => wrk_2d_2 , zdk1t => wrk_2d_3
+      USE wrk_nemo, ONLY:   zftw => wrk_xz_1 , zdit  => wrk_xz_2 
+      USE wrk_nemo, ONLY:   zdjt => wrk_xz_3 , zdj1t => wrk_xz_4
+      !
       INTEGER         , INTENT(in )                              ::  kt      ! ocean time-step index
       CHARACTER(len=3), INTENT(in )                              ::  cdtype  ! =TRA or TRC (tracer indicator) 
@@ -166,8 +179,9 @@
       REAL(wp) ::  zbtr, ztah, ztav
       REAL(wp) ::  zcof0, zcof1, zcof2, zcof3, zcof4
-      REAL(wp), DIMENSION(jpi,jpj) ::  zftu,  zdkt, zdk1t       ! workspace
-      REAL(wp), DIMENSION(jpi,jpk) ::  zftw, zdit, zdjt, zdj1t  ! 
-      !!----------------------------------------------------------------------
-
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 1,2,3) .OR. wrk_in_use_xz(1,2,3,4) )THEN
+         CALL ctl_stop('ldfght : requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
       !
       DO jn = 1, kjpt
@@ -285,6 +299,6 @@
 
             ! II.3 Divergence of vertical fluxes added to the horizontal divergence
-            ! --------------------------------------------------------------------
-#if ! defined key_traldf_smag            
+            ! ---------------------------------------------------------------------
+            
             IF( kaht == 1 ) THEN
                ! multiply the laplacian by the eddy diffusivity coefficient
@@ -320,41 +334,8 @@
          !                                                ! ===============
       END DO
-#else
-            IF( kaht == 2 ) THEN
-               ! multiply the laplacian by the eddy diffusivity coefficient
-               DO jk = 1, jpkm1
-                  DO ji = 2, jpim1  
-                     ! eddy coef. divided by the volume element
-                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                     ! vertical divergence
-                     ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) )
-                     ! harmonic operator applied to (pt,ps) and multiply by aht
-                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
-                  END DO
-               END DO
-            ELSEIF( kaht == 1 ) THEN
-               ! second call, no multiplication
-               DO jk = 1, jpkm1
-                  DO ji = 2, jpim1
-                     ! inverse of the volume element
-                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                     ! vertical divergence
-                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
-                     ! harmonic operator applied to (pt,ps) 
-                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
-                  END DO
-               END DO
-            ELSE
-               IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
-               IF(lwp) WRITE(numout,*) '         We stop'
-               STOP 'ldfght'
-            ENDIF
-            !                                             ! ===============
-         END DO                                           !   End of slab
-         !                                                ! ===============
-      END DO
-
-
-#endif      
+      !
+      IF( wrk_not_released(2, 1,2,3)   .OR.   &
+          wrk_not_released_xz(1,2,3,4) )   CALL ctl_stop('ldfght : failed to release workspace arrays.')
+      !
    END SUBROUTINE ldfght
 
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2926)
@@ -156,11 +156,4 @@
       IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
 #endif
-#if defined key_traldf_c3d && key_traldf_smag
-       CALL ldf_tra_smag( kstp )      ! Smagorinsky diffusivity
-#  endif
-#if defined key_dynldf_c3d && key_dynldf_smag
-       CALL ldf_dyn_smag( kstp )      ! Smagorinsky viscosity
-#  endif
-
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 2926)
@@ -60,6 +60,4 @@
    USE ldfslp           ! iso-neutral slopes               (ldf_slp routine)
    USE ldfeiv           ! eddy induced velocity coef.      (ldf_eiv routine)
-   USE ldftra_smag      ! smagorinsky diffusivity   .      (ldf_tra_smag routine)
-   USE ldfdyn_smag      ! smagorinsky viscosity            (ldf_dyn_smag routine)
 
    USE zdftmx           ! tide-induced vertical mixing     (zdf_tmx routine)
Index: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/TOOLS/COMPILE/cfg.txt
===================================================================
--- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/TOOLS/COMPILE/cfg.txt	(revision 2925)
+++ branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/TOOLS/COMPILE/cfg.txt	(revision 2926)
@@ -1,2 +1,3 @@
+GYRE OPA_SRC
 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
 GYRE_LOBSTER OPA_SRC TOP_SRC
@@ -5,4 +6,2 @@
 ORCA2_LIM3 OPA_SRC LIM_SRC_3
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
-GYRE OPA_SRC
-GYRE_5 OPA_SRC
