Index: trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90
===================================================================
--- trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90	(revision 28)
+++ trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90	(revision 32)
@@ -17,10 +17,4 @@
       !!
       !! ** Method :
-      !!                   **** W A R N I N G ****
-      !!                     ORCA OCEAN VERSION
-      !!    This method is relevant ONLY for the grid build by the method
-      !!               given in the 'Reference' section.
-      !!                   **** W A R N I N G ****
-      !!
       !!      2D eddy viscosity coefficients ( longitude, latitude )
       !!
@@ -55,7 +49,6 @@
 
          zdx_max = MAXVAL( e1t(:,:) )
-#if defined key_mpp
-         CALL mpp_max( zdx_max )
-#endif
+         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
+
          IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
          IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
@@ -96,7 +89,6 @@
 
          zdx_max = MAXVAL( e1u(:,:) )
-#if defined key_mpp
-         CALL mpp_max( zdx_max )
-#endif
+         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
+
          IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
          IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
@@ -141,15 +133,13 @@
 
       !! * Arguments
-      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
+      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
 
       !! * Local variables
-      INTEGER :: ji, jj      ! dummy loop indices
-      INTEGER :: inumcf, iost, iim, ijm
-      INTEGER :: jn
-      INTEGER :: ifreq, il1, il2, ij, ii, inorth, isouth
-      INTEGER :: ipi, ipj, iumout, iwork, icompt, ibtest, ikmax
-      INTEGER :: ijpt0, ijpt1, iipt0, iipt1
-      INTEGER ,DIMENSION(jpidta,jpidta) :: idata
-      INTEGER ,DIMENSION(jpi   ,jpj   ) :: icof
+      INTEGER ::   ji, jj, jn      ! dummy loop indices
+      INTEGER ::   inum = 11       ! temporary logical unit
+      INTEGER ::   iost, iim, ijm
+      INTEGER ::   ifreq, il1, il2, ij, ii
+      INTEGER, DIMENSION(jpidta,jpidta) ::   idata
+      INTEGER, DIMENSION(jpi   ,jpj   ) ::   icof
 
       REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk
@@ -173,11 +163,10 @@
       ! ===================== equatorial strip (20N-20S) defined at t-points
 
-      inumcf = 15
-      OPEN( UNIT=inumcf,FILE='ahmcoef',STATUS='OLD',   &
-            FORM='FORMATTED', ACCESS='SEQUENTIAL', ERR=111 ,   &
-            IOSTAT= iost) 
+      OPEN( UNIT=inum, FILE='ahmcoef', STATUS='OLD',   &
+         &  FORM='FORMATTED', ACCESS='SEQUENTIAL', ERR=111 ,   &
+         &  IOSTAT= iost ) 
       IF( iost == 0 ) THEN
          IF(lwp) WRITE(numout,*) '     file   : ahmcoef open ok'
-         IF(lwp) WRITE(numout,*) '     unit   = ', inumcf
+         IF(lwp) WRITE(numout,*) '     unit   = ', inum
          IF(lwp) WRITE(numout,*) '     status = OLD'
          IF(lwp) WRITE(numout,*) '     form   = FORMATTED'
@@ -195,16 +184,16 @@
       ENDIF
 
-      REWIND inumcf
-      READ(inumcf,9101) clexp, iim, ijm
-      READ(inumcf,'(/)')
+      REWIND inum
+      READ(inum,9101) clexp, iim, ijm
+      READ(inum,'(/)')
       ifreq = 40
       il1 = 1
       DO jn = 1, jpidta/ifreq+1
-         READ(inumcf,'(/)')
+         READ(inum,'(/)')
          il2 = MIN( jpidta, il1+ifreq-1 )
-         READ(inumcf,9201) ( ii, ji = il1, il2, 5 )
-         READ(inumcf,'(/)')
+         READ(inum,9201) ( ii, ji = il1, il2, 5 )
+         READ(inum,'(/)')
          DO jj = jpjdta, 1, -1
-            READ(inumcf,9202) ij, ( idata(ji,jj), ji = il1, il2 )
+            READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
          END DO
          il1 = il1 + ifreq
Index: trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90
===================================================================
--- trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90	(revision 28)
+++ trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90	(revision 32)
@@ -49,7 +49,6 @@
 
          zdx_max = MAXVAL( e1t(:,:) )
-#if defined key_mpp
-         CALL mpp_max( zdx_max )
-#endif
+         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
+
          IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
          IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
@@ -105,7 +104,6 @@
 
          zdx_max = MAXVAL( e1u(:,:) )
-#if defined key_mpp
-         CALL mpp_max( zdx_max )
-#endif
+         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
+
          IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
          IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
@@ -133,7 +131,6 @@
          ELSE                            ! partial steps or s-ccordinate
             zc = MAXVAL( fsdept(:,:,jpkm1) )
-#if defined key_mpp
-            CALL mpp_max( zc )
-#endif
+            IF( lk_mpp )   CALL mpp_max( zc )   ! max over the global domain
+
             zc = 1. / (  1. - EXP( ( zc - zh ) / zh )  )
             DO jk = 2, jpkm1
@@ -188,12 +185,13 @@
 
       !! * local variables
-      INTEGER ::inumcf, iost, iim, ijm
-      INTEGER ::ji,jj,jk, jn
-      INTEGER ::ifreq, il1, il2, ij, ii
-      INTEGER ,DIMENSION(jpidta, jpjdta) :: idata
-      INTEGER ,DIMENSION(jpi   , jpj   ) :: icof
+      INTEGER ::   ji, jj, jk, jn      ! dummy loop indices
+      INTEGER ::   inum = 11           ! temporary logical unit
+      INTEGER ::   iost, iim, ijm
+      INTEGER ::   ifreq, il1, il2, ij, ii
+      INTEGER, DIMENSION(jpidta, jpjdta) ::   idata
+      INTEGER, DIMENSION(jpi   , jpj   ) ::   icof
 
       REAL(wp) ::   zahmeq, zcoff, zcoft, zmsk
-      REAL(wp) ::   zcoef(jpk)
+      REAL(wp), DIMENSION(jpk) ::   zcoef
 
       CHARACTER (len=15) ::   clexp
@@ -211,12 +209,11 @@
       ! ===================== equatorial strip (20N-20S) defined at t-points
 
-      inumcf = 15
-      OPEN( UNIT=inumcf,FILE='ahmcoef',STATUS='OLD',   &
-            FORM='FORMATTED', ACCESS='SEQUENTIAL', ERR=111 ,   &
-            IOSTAT= iost)
+      OPEN( UNIT=inum, FILE='ahmcoef', STATUS='OLD',   &
+         &  FORM='FORMATTED', ACCESS='SEQUENTIAL', ERR=111 ,   &
+         &  IOSTAT= iost)
       IF( iost == 0 ) THEN
          IF(lwp) THEN
             WRITE(numout,*) '     file   : ahmcoef open ok'
-            WRITE(numout,*) '     unit   = ', inumcf
+            WRITE(numout,*) '     unit   = ', inum
             WRITE(numout,*) '     status = OLD'
             WRITE(numout,*) '     form   = FORMATTED'
@@ -235,16 +232,16 @@
       ENDIF
 
-      REWIND inumcf
-      READ(inumcf,9101) clexp, iim, ijm
-      READ(inumcf,'(/)')
+      REWIND inum
+      READ(inum,9101) clexp, iim, ijm
+      READ(inum,'(/)')
       ifreq = 40
       il1 = 1
       DO jn = 1, jpidta/ifreq+1
-         READ(inumcf,'(/)')
+         READ(inum,'(/)')
          il2 = MIN( jpidta, il1+ifreq-1 )
-         READ(inumcf,9201) ( ii, ji = il1, il2, 5 )
-         READ(inumcf,'(/)')
+         READ(inum,9201) ( ii, ji = il1, il2, 5 )
+         READ(inum,'(/)')
          DO jj = jpjdta, 1, -1
-            READ(inumcf,9202) ij, ( idata(ji,jj), ji = il1, il2 )
+            READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
          END DO
          il1 = il1 + ifreq
Index: trunk/NEMO/OPA_SRC/LDF/ldfeiv_substitute.h90
===================================================================
--- trunk/NEMO/OPA_SRC/LDF/ldfeiv_substitute.h90	(revision 28)
+++ trunk/NEMO/OPA_SRC/LDF/ldfeiv_substitute.h90	(revision 32)
@@ -1,3 +1,3 @@
-#if defined key_traldf_eiv
+#if defined key_traldf_eiv   ||   defined key_esopa
    !!----------------------------------------------------------------------
    !!                   ***  ldfeiv_substitute.h90  ***
Index: trunk/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 28)
+++ trunk/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 32)
@@ -29,8 +29,8 @@
 
    !! * Share module variables
-   LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp = .TRUE.     ! slopes flag
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &
-      uslp, wslpi,         &  ! i_slope at U- and W-points
-      vslp, wslpj             ! j-slope at V- and W-points
+   LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp = .TRUE.     !: slopes flag
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
+      uslp, wslpi,         &  !: i_slope at U- and W-points
+      vslp, wslpj             !: j-slope at V- and W-points
    
    !! * Module variables
@@ -95,5 +95,6 @@
 
       !! * Local declarations
-      INTEGER  ::   ji, jj, jk      ! dummy loop indices
+      INTEGER  ::   ji, jj, jk              ! dummy loop indices
+      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integer
 #if defined key_partial_steps
       INTEGER  ::   iku, ikv  ! temporary integers
@@ -120,6 +121,6 @@
       
       zeps  =  1.e-20
-      zmg   = -1.0 / g
-      zm05g = -0.5 / g
+      zmg   = -1.0 / grav
+      zm05g = -0.5 / grav
 
       zww(:,:,:) = 0.e0
@@ -177,5 +178,5 @@
          
          ! Local vertical density gradient evaluated from N^2
-         ! zwy = d/dz(prd)= - ( prd ) / g * mk(pn2) -- at t point
+         ! zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
 
          DO jj = 1, jpj
@@ -323,5 +324,5 @@
          
          ! Local vertical density gradient evaluated from N^2
-         !     zwy = d/dz(prd)= - mk ( prd ) / g * pn2 -- at w point
+         !     zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -450,14 +451,18 @@
             !
             !                                             ! Gibraltar Strait
-            uslp ( mi0(69):mi1(71) , mj0(50):mj1(53) , jk ) = 0.0e0
-            vslp ( mi0(68):mi1(71) , mj0(51):mj1(53) , jk ) = 0.0e0
-            wslpi( mi0(69):mi1(71) , mj0(51):mj1(53) , jk ) = 0.0e0
-            wslpj( mi0(69):mi1(71) , mj0(51):mj1(53) , jk ) = 0.0e0
+            ij0 =  50   ;   ij1 =  53
+            ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
+            ij0 =  51   ;   ij1 =  53
+            ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
+            ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
+            ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
 
             !                                             ! Mediterrannean Sea
-            uslp ( mi0(71):mi1(90) , mj0(49):mj1(56) , jk ) = 0.0e0
-            vslp ( mi0(70):mi1(90) , mj0(50):mj1(56) , jk ) = 0.0e0
-            wslpi( mi0(71):mi1(90) , mj0(50):mj1(56) , jk ) = 0.0e0
-            wslpj( mi0(71):mi1(90) , mj0(50):mj1(56) , jk ) = 0.0e0
+            ij0 =  49   ;   ij1 =  56
+            ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
+            ij0 =  50   ;   ij1 =  56
+            ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
+            ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
+            ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk ) = 0.e0
          ENDIF
          !                                             ! ===============
@@ -527,6 +532,6 @@
 
       zeps  =  1.e-20
-      zmg   = -1.0 / g
-      zm05g = -0.5 / g
+      zmg   = -1.0 / grav
+      zm05g = -0.5 / grav
 
 
@@ -567,5 +572,5 @@
 
       ! Local vertical density gradient evaluated from N^2
-      ! zwy = d/dz(prd)= - ( prd ) / g * mk(pn2) -- at t point
+      ! zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
 
       !-----------------------------------------------------------------------
@@ -619,5 +624,5 @@
 
       ! Local vertical density gradient evaluated from N^2
-      !     zwy = d/dz(prd)= - ( prd ) / g * mk(pn2) -- at t point
+      !     zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
       zwy ( :, jpj) = 0.0e0
       zwy ( jpi, :) = 0.0e0
@@ -673,5 +678,5 @@
 
       ! Local vertical density gradient evaluated from N^2
-      ! zwy = d/dz(prd)= - mk ( prd ) / g * pn2 -- at w point
+      ! zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
 # if defined key_vectopt_loop   &&   ! defined key_autotasking
       jj = 1
@@ -809,10 +814,10 @@
    !!   Dummy module :                 NO Rotation of lateral mixing tensor
    !!------------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    ! slopes flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag
 CONTAINS
    SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine
       INTEGER, INTENT(in) :: kt 
       REAL,DIMENSION(:,:,:), INTENT(in) :: prd, pn2
-      WRITE(*,*) kt, prd, pn2
+      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
    END SUBROUTINE ldf_slp
 #endif
Index: trunk/NEMO/OPA_SRC/LDF/ldftra_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/LDF/ldftra_oce.F90	(revision 28)
+++ trunk/NEMO/OPA_SRC/LDF/ldftra_oce.F90	(revision 32)
@@ -22,35 +22,35 @@
    !!----------------------------------------------------------------------
 
-   LOGICAL , PUBLIC ::              & !!! ** lateral mixing namelist (nam_traldf) **
-      ln_traldf_lap   = .TRUE.  ,   &  ! laplacian operator
-      ln_traldf_bilap = .FALSE. ,   &  ! bilaplacian operator
-      ln_traldf_level = .FALSE. ,   &  ! iso-level direction
-      ln_traldf_hor   = .FALSE. ,   &  ! horizontal (geopotential) direction
-      ln_traldf_iso   = .TRUE.         ! iso-neutral direction
+   LOGICAL , PUBLIC ::              & !!: ** lateral mixing namelist (nam_traldf) **
+      ln_traldf_lap   = .TRUE.  ,   &  !: laplacian operator
+      ln_traldf_bilap = .FALSE. ,   &  !: bilaplacian operator
+      ln_traldf_level = .FALSE. ,   &  !: iso-level direction
+      ln_traldf_hor   = .FALSE. ,   &  !: horizontal (geopotential) direction
+      ln_traldf_iso   = .TRUE.         !: iso-neutral direction
 
-   REAL(wp), PUBLIC ::              & !!! ** lateral mixing namelist (namldf) **
-      aht0  = 2000._wp     ,        &  ! lateral eddy diffusivity (m2/s)
-      ahtb0 =    0._wp     ,        &  ! lateral background eddy diffusivity (m2/s)
-      aeiv0 = 2000._wp                 ! eddy induced velocity coefficient (m2/s)
+   REAL(wp), PUBLIC ::              & !!: ** lateral mixing namelist (namldf) **
+      aht0  = 2000._wp     ,        &  !: lateral eddy diffusivity (m2/s)
+      ahtb0 =    0._wp     ,        &  !: lateral background eddy diffusivity (m2/s)
+      aeiv0 = 2000._wp                 !: eddy induced velocity coefficient (m2/s)
 
-   LOGICAL , PUBLIC ::              &  ! flag of the lateral diff. scheme used 
-      l_traldf_lap         ,        &  ! iso-level laplacian operator
-      l_traldf_bilap       ,        &  ! iso-level bilaplacian operator
-      l_traldf_bilapg      ,        &  ! geopotential bilap. (s-coord)
-      l_traldf_iso         ,        &  ! iso-neutral laplacian or horizontal lapacian (s-coord)
-      l_trazdf_iso         ,        &  ! idem for the vertical component
-      l_trazdf_iso_vo      ,        &  ! idem with vectopt_memory
-      l_traldf_iso_zps                 ! iso-neutral laplacian (partial steps)
+   LOGICAL , PUBLIC ::              &  !: flag of the lateral diff. scheme used 
+      l_traldf_lap         ,        &  !: iso-level laplacian operator
+      l_traldf_bilap       ,        &  !: iso-level bilaplacian operator
+      l_traldf_bilapg      ,        &  !: geopotential bilap. (s-coord)
+      l_traldf_iso         ,        &  !: iso-neutral laplacian or horizontal lapacian (s-coord)
+      l_trazdf_iso         ,        &  !: idem for the vertical component
+      l_trazdf_iso_vo      ,        &  !: idem with vectopt_memory
+      l_traldf_iso_zps                 !: iso-neutral laplacian (partial steps)
 
 #if defined key_traldf_c3d
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  ! ** 3D coefficients **
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !: ** 3D coefficients **
 #elif defined key_traldf_c2d
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   &  ! ** 2D coefficients **
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   &  !: ** 2D coefficients **
 #elif defined key_traldf_c1d
-   REAL(wp), PUBLIC, DIMENSION(jpk)         ::   &  ! ** 1D coefficients **
+   REAL(wp), PUBLIC, DIMENSION(jpk)         ::   &  !: ** 1D coefficients **
 #else
-   REAL(wp), PUBLIC                         ::   &  ! ** 0D coefficients **
+   REAL(wp), PUBLIC                         ::   &  !: ** 0D coefficients **
 #endif
-      ahtt, ahtu, ahtv, ahtw                ! T-, U-, V-, W-points coefficients
+      ahtt, ahtu, ahtv, ahtw                !: T-, U-, V-, W-points coefficients
 
 
@@ -59,22 +59,20 @@
    !!   'key_traldf_eiv'                              eddy induced velocity
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_traldf_eiv   = .TRUE.   ! eddy induced velocity flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_traldf_eiv   = .TRUE.   !: eddy induced velocity flag
       
 # if defined key_traldf_c3d
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  ! ** 3D coefficients **
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !: ** 3D coefficients **
 # elif defined key_traldf_c2d
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   &  ! ** 2D coefficients **
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   &  !: ** 2D coefficients **
 # elif defined key_traldf_c1d
-   REAL(wp), PUBLIC, DIMENSION(jpk)         ::   &  ! ** 1D coefficients **
+   REAL(wp), PUBLIC, DIMENSION(jpk)         ::   &  !: ** 1D coefficients **
 # else
-   REAL(wp), PUBLIC                         ::   &  ! ** 0D coefficients **
+   REAL(wp), PUBLIC                         ::   &  !: ** 0D coefficients **
 # endif
-      aeiu, aeiv, aeiw                              ! U-, V-, W-points  induced velocity coef. (m2/s)
+      aeiu, aeiv, aeiw                              !: U-, V-, W-points  induced velocity coef. (m2/s)
 
 # if defined key_diaeiv
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::    &
-      u_eiv,      &  ! The three component of the eddy induced velocity (m/s)
-      v_eiv,      &  ! saved for diagnostics and/or outputs
-      w_eiv          ! 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::    &  !:
+      u_eiv, v_eiv, w_eiv     !: The three component of the eddy induced velocity (m/s)
 # endif
 
@@ -83,5 +81,5 @@
    !!   Default option :                           NO eddy induced velocity
    !!----------------------------------------------------------------------
-   LOGICAL , PUBLIC, PARAMETER ::   lk_traldf_eiv   = .FALSE.   ! eddy induced velocity flag
+   LOGICAL , PUBLIC, PARAMETER ::   lk_traldf_eiv   = .FALSE.   !: eddy induced velocity flag
    REAL(wp), PUBLIC ::   aeiu, aeiv, aeiw
 #endif
