Index: NEMO/trunk/src/OCE/LDF/ldftra.F90
===================================================================
--- NEMO/trunk/src/OCE/LDF/ldftra.F90	(revision 13982)
+++ NEMO/trunk/src/OCE/LDF/ldftra.F90	(revision 14072)
@@ -2,9 +2,9 @@
    !!======================================================================
    !!                       ***  MODULE  ldftra  ***
-   !! Ocean physics:  lateral diffusivity coefficients 
+   !! Ocean physics:  lateral diffusivity coefficients
    !!=====================================================================
    !! History :       ! 1997-07  (G. Madec)  from inimix.F split in 2 routines
    !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
-   !!            2.0  ! 2005-11  (G. Madec)  
+   !!            2.0  ! 2005-11  (G. Madec)
    !!            3.7  ! 2013-12  (F. Lemarie, G. Madec)  restructuration/simplification of aht/aeiv specification,
    !!                 !                                  add velocity dependent coefficient and optional read in file
@@ -13,6 +13,6 @@
    !!----------------------------------------------------------------------
    !!   ldf_tra_init : initialization, namelist read, and parameters control
-   !!   ldf_tra      : update lateral eddy diffusivity coefficients at each time step 
-   !!   ldf_eiv_init : initialization of the eiv coeff. from namelist choices 
+   !!   ldf_tra      : update lateral eddy diffusivity coefficients at each time step
+   !!   ldf_eiv_init : initialization of the eiv coeff. from namelist choices
    !!   ldf_eiv      : time evolution of the eiv coefficients (function of the growth rate of baroclinic instability)
    !!   ldf_eiv_trp  : add to the input ocean transport the contribution of the EIV parametrization
@@ -23,5 +23,5 @@
    USE phycst          ! physical constants
    USE ldfslp          ! lateral diffusion: slope of iso-neutral surfaces
-   USE ldfc1d_c2d      ! lateral diffusion: 1D & 2D cases 
+   USE ldfc1d_c2d      ! lateral diffusion: 1D & 2D cases
    USE diaptr
    !
@@ -40,6 +40,6 @@
    PUBLIC   ldf_eiv_trp    ! called by traadv.F90
    PUBLIC   ldf_eiv_dia    ! called by traldf_iso and traldf_iso_triad.F90
-   
-   !                                   !!* Namelist namtra_ldf : lateral mixing on tracers * 
+
+   !                                   !!* Namelist namtra_ldf : lateral mixing on tracers *
    !                                    != Operator type =!
    LOGICAL , PUBLIC ::   ln_traldf_OFF       !: no operator: No explicit diffusion
@@ -52,5 +52,5 @@
    !		       	                      != iso-neutral options =!
 !  LOGICAL , PUBLIC ::   ln_traldf_triad     !: griffies triad scheme                    (see ldfslp)
-   LOGICAL , PUBLIC ::   ln_traldf_msc       !: Method of Stabilizing Correction 
+   LOGICAL , PUBLIC ::   ln_traldf_msc       !: Method of Stabilizing Correction
 !  LOGICAL , PUBLIC ::   ln_triad_iso        !: pure horizontal mixing in ML             (see ldfslp)
 !  LOGICAL , PUBLIC ::   ln_botmix_triad     !: mixing on bottom                         (see ldfslp)
@@ -59,5 +59,5 @@
    !                                    !=  Coefficients =!
    INTEGER , PUBLIC ::   nn_aht_ijk_t        !: choice of time & space variations of the lateral eddy diffusivity coef.
-   !                                            !  time invariant coefficients:  aht_0 = 1/2  Ud*Ld   (lap case) 
+   !                                            !  time invariant coefficients:  aht_0 = 1/2  Ud*Ld   (lap case)
    !                                            !                                bht_0 = 1/12 Ud*Ld^3 (blp case)
    REAL(wp), PUBLIC ::      rn_Ud               !: lateral diffusive velocity  [m/s]
@@ -72,5 +72,5 @@
    REAL(wp), PUBLIC ::      rn_Ue               !: lateral diffusive velocity  [m/s]
    REAL(wp), PUBLIC ::      rn_Le               !: lateral diffusive length    [m]
-   
+
    !                                  ! Flag to control the type of lateral diffusive operator
    INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   ! error in specification of lateral diffusion
@@ -106,5 +106,5 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_tra_init  ***
-      !! 
+      !!
       !! ** Purpose :   initializations of the tracer lateral mixing coeff.
       !!
@@ -116,5 +116,5 @@
       !!    nn_aht_ijk_t  =  0 => = constant
       !!                  !
-      !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth 
+      !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth
       !!                  !
       !!                  =-20 => = F(i,j)   = shape read in 'eddy_diffusivity.nc' file
@@ -126,6 +126,6 @@
       !!                  = 31    = F(i,j,k,t) = F(local velocity) (  1/2  |u|e     laplacian operator
       !!                                                           or 1/12 |u|e^3 bilaplacian operator )
-      !!              * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init 
-      !!            
+      !!              * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init
+      !!
       !! ** action  : ahtu, ahtv initialized one for all or l_ldftra_time set to true
       !!              aeiu, aeiv initialized one for all or l_ldfeiv_time set to true
@@ -148,5 +148,5 @@
          WRITE(numout,*) '~~~~~~~~~~~~ '
       ENDIF
-      
+
       !
       !  Choice of lateral tracer physics
@@ -182,5 +182,5 @@
       !
       !
-      ! Operator and its acting direction   (set nldf_tra)  
+      ! Operator and its acting direction   (set nldf_tra)
       ! =================================
       !
@@ -210,5 +210,5 @@
             ENDIF
             IF ( ln_zps ) THEN                  ! z-coordinate with partial step
-               IF ( ln_traldf_lev   )   ierr     = 1          ! iso-level not allowed 
+               IF ( ln_traldf_lev   )   ierr     = 1          ! iso-level not allowed
                IF ( ln_traldf_hor   )   nldf_tra = np_lap     ! horizontal             (no rotation)
                IF ( ln_traldf_iso   )   nldf_tra = np_lap_i   ! iso-neutral: standard     (rotation)
@@ -231,5 +231,5 @@
             ENDIF
             IF ( ln_zps ) THEN                  ! z-coordinate with partial step
-               IF ( ln_traldf_lev   )   ierr     = 1          ! iso-level not allowed 
+               IF ( ln_traldf_lev   )   ierr     = 1          ! iso-level not allowed
                IF ( ln_traldf_hor   )   nldf_tra = np_blp     ! horizontal             (no rotation)
                IF ( ln_traldf_iso   )   nldf_tra = np_blp_i   ! iso-neutral: standard  (   rotation)
@@ -249,5 +249,5 @@
            !
       IF(  nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. &
-         & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required 
+         & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required
       !
       IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN     ! iso-neutral bilaplacian need MSC
@@ -270,5 +270,5 @@
 
       !
-      !  Space/time variation of eddy coefficients 
+      !  Space/time variation of eddy coefficients
       ! ===========================================
       !
@@ -286,5 +286,5 @@
          IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays')
          !
-         ahtu(:,:,jpk) = 0._wp                     ! last level always 0  
+         ahtu(:,:,jpk) = 0._wp                     ! last level always 0
          ahtv(:,:,jpk) = 0._wp
          !.
@@ -363,5 +363,5 @@
          END SELECT
          !
-         IF( .NOT.l_ldftra_time ) THEN             !* No time variation 
+         IF( .NOT.l_ldftra_time ) THEN             !* No time variation
             IF(     ln_traldf_lap ) THEN                 !   laplacian operator (mask only)
                ahtu(:,:,1:jpkm1) =       ahtu(:,:,1:jpkm1)   * umask(:,:,1:jpkm1)
@@ -381,5 +381,5 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_tra  ***
-      !! 
+      !!
       !! ** Purpose :   update at kt the tracer lateral mixing coeff. (aht and aeiv)
       !!
@@ -395,6 +395,6 @@
       !!              * time varying EIV coefficients: call to ldf_eiv routine
       !!
-      !! ** action  :   ahtu, ahtv   update at each time step   
-      !!                aeiu, aeiv      -       -     -    -   (if ln_ldfeiv=T) 
+      !! ** action  :   ahtu, ahtv   update at each time step
+      !!                aeiu, aeiv      -       -     -    -   (if ln_ldfeiv=T)
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step
@@ -420,11 +420,11 @@
             ahtu(:,:,1) = aeiu(:,:,1)
             ahtv(:,:,1) = aeiv(:,:,1)
-         ELSE                                            ! compute aht. 
+         ELSE                                            ! compute aht.
             CALL ldf_eiv( kt, aht0, ahtu, ahtv, Kmm )
          ENDIF
          !
-         z1_f20   = 1._wp / (  2._wp * omega * SIN( rad * 20._wp )  )   ! 1 / ff(20 degrees)   
+         z1_f20   = 1._wp / (  2._wp * omega * SIN( rad * 20._wp )  )   ! 1 / ff(20 degrees)
          zaht_min = 0.2_wp * aht0                                       ! minimum value for aht
-         zDaht    = aht0 - zaht_min                                      
+         zDaht    = aht0 - zaht_min
          ! NOTE: [tiling-comms-merge] Change needed to preserve results with respect to the trunk
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
@@ -480,5 +480,5 @@
       !!    nn_aei_ijk_t  =  0 => = constant
       !!                  !
-      !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth 
+      !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth
       !!                  !
       !!                  =-20 => = F(i,j)   = shape read in 'eddy_diffusivity.nc' file
@@ -547,5 +547,5 @@
          !                                != Specification of space-time variations of eaiu, aeiv
          !
-         aeiu(:,:,jpk) = 0._wp               ! last level always 0  
+         aeiu(:,:,jpk) = 0._wp               ! last level always 0
          aeiv(:,:,jpk) = 0._wp
          !                                   ! value of EIV coef. (laplacian operator)
@@ -609,5 +609,5 @@
          END SELECT
          !
-         IF( .NOT.l_ldfeiv_time ) THEN             !* mask if No time variation 
+         IF( .NOT.l_ldfeiv_time ) THEN             !* mask if No time variation
             DO jk = 1, jpkm1
                aeiu(:,:,jk) = aeiu(:,:,jk) * umask(:,:,jk)
@@ -617,5 +617,5 @@
          !
       ENDIF
-      !                    
+      !
    END SUBROUTINE ldf_eiv_init
 
@@ -649,11 +649,11 @@
       IF( ln_traldf_triad ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            ! Take the max of N^2 and zero then take the vertical sum 
-            ! of the square root of the resulting N^2 ( required to compute 
-            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
+            ! Take the max of N^2 and zero then take the vertical sum
+            ! of the square root of the resulting N^2 ( required to compute
+            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
             zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
             zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm)
             ! Compute elements required for the inverse time scale of baroclinic
-            ! eddies using the isopycnal slopes calculated in ldfslp.F : 
+            ! eddies using the isopycnal slopes calculated in ldfslp.F :
             ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
             ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
@@ -663,11 +663,11 @@
       ELSE
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            ! Take the max of N^2 and zero then take the vertical sum 
-            ! of the square root of the resulting N^2 ( required to compute 
-            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
+            ! Take the max of N^2 and zero then take the vertical sum
+            ! of the square root of the resulting N^2 ( required to compute
+            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
             zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
             zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm)
             ! Compute elements required for the inverse time scale of baroclinic
-            ! eddies using the isopycnal slopes calculated in ldfslp.F : 
+            ! eddies using the isopycnal slopes calculated in ldfslp.F :
             ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
             ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
@@ -693,5 +693,5 @@
       END_2D
       CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1.0_wp )       ! lateral boundary condition
-      !               
+      !
       DO_2D( 0, 0, 0, 0 )                       !== aei at u- and v-points  ==!
          paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj  ) ) * umask(ji,jj,1)
@@ -704,5 +704,5 @@
          paeiv(:,:,jk) = paeiv(:,:,1) * vmask(:,:,jk)
       END DO
-      !  
+      !
    END SUBROUTINE ldf_eiv
 
@@ -711,11 +711,11 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_eiv_trp  ***
-      !! 
-      !! ** Purpose :   add to the input ocean transport the contribution of 
+      !!
+      !! ** Purpose :   add to the input ocean transport the contribution of
       !!              the eddy induced velocity parametrization.
       !!
       !! ** Method  :   The eddy induced transport is computed from a flux stream-
       !!              function which depends on the slope of iso-neutral surfaces
-      !!              (see ldf_slp). For example, in the i-k plan : 
+      !!              (see ldf_slp). For example, in the i-k plan :
       !!                   psi_uw = mk(aeiu) e2u mi(wslpi)   [in m3/s]
       !!                   Utr_eiv = - dk[psi_uw]
@@ -748,5 +748,5 @@
       ENDIF
 
-      
+
       zpsi_uw(:,:, 1 ) = 0._wp   ;   zpsi_vw(:,:, 1 ) = 0._wp
       zpsi_uw(:,:,jpk) = 0._wp   ;   zpsi_vw(:,:,jpk) = 0._wp
@@ -794,5 +794,5 @@
       !
 !!gm I don't like this routine....   Crazy  way of doing things, not optimal at all...
-!!gm     to be redesigned....   
+!!gm     to be redesigned....
       !                                                  !==  eiv stream function: output  ==!
 !!gm      CALL iom_put( "psi_eiv_uw", psi_uw )                 ! output
@@ -826,5 +826,5 @@
             zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:)
          END DO
-         CALL iom_put( "weiv_masstr" , zw3d )  
+         CALL iom_put( "weiv_masstr" , zw3d )
       ENDIF
       !
@@ -832,16 +832,16 @@
          zw3d(:,:,:) = 0.e0
          DO jk = 1, jpkm1
-            zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 
+            zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) )
          END DO
          CALL iom_put( "ueiv_masstr", zw3d )                  ! mass transport in i-direction
       ENDIF
       !
-      zztmp = 0.5_wp * rho0 * rcp 
+      zztmp = 0.5_wp * rho0 * rcp
       IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
-        zw2d(:,:)   = 0._wp 
-        zw3d(:,:,:) = 0._wp 
+        zw2d(:,:)   = 0._wp
+        zw3d(:,:,:) = 0._wp
         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
            zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
-              &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji+1,jj,jk,jp_tem,Kmm) ) 
+              &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji+1,jj,jk,jp_tem,Kmm) )
            zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
         END_3D
@@ -853,14 +853,14 @@
          zw3d(:,:,:) = 0.e0
          DO jk = 1, jpkm1
-            zw3d(:,:,jk) = rho0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) 
+            zw3d(:,:,jk) = rho0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) )
          END DO
          CALL iom_put( "veiv_masstr", zw3d )                  ! mass transport in i-direction
       ENDIF
       !
-      zw2d(:,:)   = 0._wp 
-      zw3d(:,:,:) = 0._wp 
+      zw2d(:,:)   = 0._wp
+      zw3d(:,:,:) = 0._wp
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
-            &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji,jj+1,jk,jp_tem,Kmm) ) 
+            &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji,jj+1,jk,jp_tem,Kmm) )
          zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
       END_3D
@@ -872,9 +872,9 @@
       zztmp = 0.5_wp * 0.5
       IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
-        zw2d(:,:) = 0._wp 
-        zw3d(:,:,:) = 0._wp 
+        zw2d(:,:) = 0._wp
+        zw3d(:,:,:) = 0._wp
         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
            zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
-              &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji+1,jj,jk,jp_sal,Kmm) ) 
+              &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji+1,jj,jk,jp_sal,Kmm) )
            zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
         END_3D
@@ -882,9 +882,9 @@
         CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                ! salt transport in i-direction
       ENDIF
-      zw2d(:,:) = 0._wp 
-      zw3d(:,:,:) = 0._wp 
+      zw2d(:,:) = 0._wp
+      zw3d(:,:,:) = 0._wp
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
-            &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji,jj+1,jk,jp_sal,Kmm) ) 
+            &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji,jj+1,jk,jp_sal,Kmm) )
          zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
       END_3D
