Index: /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 7339)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 7340)
@@ -39,5 +39,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
-   USE sbcwave,  ONLY:  usd3dt, vsd3dt,wsd3d
+   USE sbcwave         ! Stokes velocities
 
    IMPLICIT NONE
@@ -162,5 +162,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
-      REAL(wp) ::  dsshnu, dsshnv
+      REAL(wp) ::  zdsshu, zdsshv
       !!----------------------------------------------------------------------
       !
@@ -209,21 +209,17 @@
          END DO
       ENDIF
-!
-!     In case ln_wave and ln_sdw, the surface vertical velocity is modified
-!     accounting for Sokes Drift velocity
-!
-      IF ( ln_wave .AND. ln_sdw )  THEN
-      ! Compute d(ssh)/dx  and d(ssh)/dy  
-      ! Compute the surface vertical velocity accounting for the Stokes Drift
-      !---------------------------------
-       DO jj = 2 , jpjm1
-        DO ji = 2 , jpim1
-          dsshnu = ( ssha(ji+1,jj) - ssha(ji,jj) ) / e1u(ji,jj)
-          dsshnv = ( ssha(ji,jj+1) - ssha(ji,jj) ) / e2v(ji,jj)
-          wn(ji,jj,1) = wn(ji,jj,1) +( usd3dt(ji,jj,1) * dsshnu     &
-                &     + vsd3dt(ji,jj,1) * dsshnv                    &
-                &     - ( wsd3d (ji,jj,1) )) * tmask(ji,jj,1)
-        ENDDO
-       ENDDO
+
+      IF( ln_wave .AND. ln_sdw ) THEN
+         ! Compute d(ssh)/dx  and d(ssh)/dy  
+         ! Compute the surface vertical velocity accounting for the Stokes Drift
+         DO jj = 1 , jpjm1
+            DO ji = 1 , fs_jpim1
+              zdsshu = ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
+              zdsshv = ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
+              wn(ji,jj,1) = wn(ji,jj,1) + ( usd3dt(ji,jj,1) * zdsshu   &
+                 &                        + vsd3dt(ji,jj,1) * zdsshv   &
+                 &        - wsd3d (ji,jj,1) ) * tmask(ji,jj,1)
+            END DO
+         END DO
       ENDIF
       !
Index: /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7339)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7340)
@@ -4,7 +4,7 @@
    !! Wave module 
    !!======================================================================
-   !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient 
-   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift 
-   !!            3.6  !   2014-09  (Clementi E, Oddo P)New Stokes Drift Computation
+   !! History :  3.3  !   2011-09  (M. Adani)  Original code: Drag Coefficient 
+   !!         :  3.4  !   2012-10  (M. Adani)  Stokes Drift 
+   !!            3.6  !   2014-09  (E. Clementi,P. Oddo) New Stokes Drift Computation
    !!----------------------------------------------------------------------
 
@@ -13,12 +13,11 @@
    !!----------------------------------------------------------------------
    USE oce            ! 
-   USE sbc_oce	      ! Surface boundary condition: ocean fields
+   USE sbc_oce	       ! Surface boundary condition: ocean fields
    USE bdy_oce        !
    USE domvvl         !
-   !
    USE iom            ! I/O manager library
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! distribued memory computing library
-   USE fldread	      ! read input fields
+   USE fldread	       ! read input fields
    USE wrk_nemo       !
    USE phycst         ! physical constants 
@@ -29,9 +28,9 @@
    PUBLIC   sbc_wave    ! routine called in sbcmod
    
-   INTEGER , PARAMETER ::   jpfld  = 4           ! number of files to read for stokes drift
-   INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_swh = 3           ! index of significant wave hight      (m)      at T-point
-   INTEGER , PARAMETER ::   jp_wmp = 4           ! index of mean wave period            (s)      at T-point
+   INTEGER, PARAMETER ::   jpfld  = 4           ! number of files to read for stokes drift
+   INTEGER, PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point
+   INTEGER, PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point
+   INTEGER, PARAMETER ::   jp_swh = 3           ! index of significant wave hight      (m)      at T-point
+   INTEGER, PARAMETER ::   jp_wmp = 4           ! index of mean wave period            (s)      at T-point
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient
@@ -43,8 +42,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tauoc_wave
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tsd2d
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: usd2d, vsd2d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     :: usd3dt, vsd3dt
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3dt , vsd3dt 
 
    !! * Substitutions
@@ -74,18 +71,16 @@
       USE zdf_oce,  ONLY : ln_zdfqiao
 
-      INTEGER, INTENT( in  ) ::   kt       ! ocean time step
-      !
-      INTEGER                ::   ierror   ! return error code
-      INTEGER                ::   ifpr, jj,ji,jk 
-      INTEGER                ::   ios      ! Local integer output status for namelist read
-      !
-      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
-      REAL(wp)                       ::  ztransp, zsp0, zk, zus, zvs
-      REAL(wp), DIMENSION(jpi,jpj)   ::  zfac 
-      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace
-      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
-      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
-                             &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read
-      !!
+      INTEGER, INTENT( in  ) :: kt       ! ocean time step
+      !
+      INTEGER            ::  ierror           ! return error code
+      INTEGER            ::  ifpr, jj,ji,jk   ! dummy loop indice 
+      INTEGER            ::  ios              ! Local integer output status for namelist read
+      CHARACTER(len=100) ::  cn_dir           ! Root directory for location of drag coefficient files
+      REAL(wp)           ::  ztransp, zfac, zsp0, zk, zus, zvs
+      REAL(wp), DIMENSION(jpi,jpj) :: zusd2dt, zvsd2dt
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3hdiv   ! 3D workspace
+      TYPE(FLD_N) ::  sn_cdg, sn_usd, sn_vsd            ! informations about the fields to be read
+      TYPE(FLD_N) ::  sn_swh, sn_wmp, sn_wnum, sn_tauoc !     "          "    "     "    "  "   "
+      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i          ! array of namelist informations on the fields to read
       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc
       !!---------------------------------------------------------------------
@@ -103,33 +98,32 @@
          IF(lwm) WRITE ( numond, namsbc_wave )
          !
-         IF ( ln_cdgw ) THEN
-            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
+         IF( ln_cdgw ) THEN
+            ALLOCATE( sf_cd(1), STAT=ierror )           ! allocate and fill sf_wave with sn_cdg
             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
             !
-                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
-            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
+            ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1) )
+            IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
             CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
             ALLOCATE( cdn_wave(jpi,jpj) )
-            cdn_wave(:,:) = 0.0
-         ENDIF
-
-         IF ( ln_tauoc ) THEN
-            ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc
+         ENDIF
+
+         IF( ln_tauoc ) THEN
+            ALLOCATE( sf_tauoc(1), STAT=ierror )           ! allocate and fill sf_wave with sn_tauoc
             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
             !
-                                   ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   )
-            IF( sn_tauoc%ln_tint )   ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) )
+            ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1) )
+            IF( sn_tauoc%ln_tint ) ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) )
             CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
             ALLOCATE( tauoc_wave(jpi,jpj) )
             tauoc_wave(:,:) = 0.0
-        ENDIF
-
-         IF ( ln_sdw ) THEN
+         ENDIF
+
+         IF( ln_sdw ) THEN
             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd;
             slf_i(jp_swh) = sn_swh ; slf_i(jp_wmp) = sn_wmp;
-            ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift
+            ALLOCATE( sf_sd(jpfld), STAT=ierror )           ! allocate and fill sf_sd with stokes drift
             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
             !
-            DO ifpr= 1, jpfld
+            DO ifpr = 1, jpfld
                ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
@@ -137,39 +131,33 @@
 
             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj) )
-            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
-            ALLOCATE( usd3dt(jpi,jpj,jpk),vsd3dt(jpi,jpj,jpk) )
+            ALLOCATE( usd3dt(jpi,jpj,jpk), vsd3dt(jpi,jpj,jpk), wsd3d(jpi,jpj,jpk) )
+            ALLOCATE( usd3d (jpi,jpj,jpk), vsd3d (jpi,jpj,jpk) )
             ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) )
-            ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) )
-            usd3d(:,:,:) = 0._wp   ;   usd2d(:,:) = 0._wp   ;   
-            vsd3d(:,:,:) = 0._wp   ;   vsd2d(:,:) = 0._wp   ; 
-            wsd3d(:,:,:) = 0._wp   ;
-            usd3dt(:,:,:) = 0._wp  ;   vsd3dt(:,:,:) = 0._wp   ;
-            swh  (:,:)   = 0._wp   ;    wmp (:,:) = 0._wp ;
-            IF ( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==!
+            usd3d(:,:,:) = 0._wp        
+            vsd3d(:,:,:) = 0._wp     
+            wsd3d(:,:,:) = 0._wp   
+            IF( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==!
                ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
-               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' )
-                                      ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
+               ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1) )
                IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
                CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
                ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) )
-               wnum(:,:) = 0._wp ; tsd2d(:,:) = 0._wp
             ENDIF
          ENDIF
       ENDIF
       !
-      IF ( ln_cdgw ) THEN              !==  Neutral drag coefficient  ==!
-         CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing
+      IF( ln_cdgw ) THEN              !==  Neutral drag coefficient  ==!
+         CALL fld_read( kt, nn_fsbc, sf_cd )        ! read from external forcing
          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
       ENDIF
 
-      IF ( ln_tauoc ) THEN             !==  Wave induced stress  ==!
-         CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing
+      IF( ln_tauoc ) THEN             !==  Wave induced stress  ==!
+         CALL fld_read( kt, nn_fsbc, sf_tauoc )     ! read wave norm stress from external forcing
          tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1)
       ENDIF
 
-      IF ( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==! 
-         !
-         CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing
+      IF( ln_sdw )  THEN              !==  Computation of the 3d Stokes Drift  ==! 
+         CALL fld_read( kt, nn_fsbc, sf_sd )        ! read wave parameters from external forcing
          swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height
          wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
@@ -180,24 +168,19 @@
          !(DOI: 10.1175/JPO-D-14-0020.1)==! 
          !
-         CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv )
          DO jk = 1, jpk
             DO jj = 1, jpj
                DO ji = 1, jpi
-               ! On T grid
-               ! Stokes transport speed estimated from Hs and Tmean
-               ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp))
-               ! Stokes surface speed
-               zsp0 = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2)
-               ! Wavenumber scale
-               zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp)
-               ! Depth attenuation
-               zfac(ji,jj) = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk))
-               END DO
-            END DO
-            !
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                 usd3dt(ji,jj,jk) = zfac(ji,jj) * zusd2dt (ji,jj) * tmask(ji,jj,jk)
-                 vsd3dt(ji,jj,jk) = zfac(ji,jj) * zvsd2dt (ji,jj) * tmask(ji,jj,jk)
+                  ! On T grid
+                  ! Stokes transport speed estimated from Hs and Tmean
+                  ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp))
+                  ! Stokes surface speed
+                  zsp0 = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2 )
+                  ! Wavenumber scale
+                  zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp)
+                  ! Depth attenuation
+                  zfac = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk))
+                  !
+                  usd3dt(ji,jj,jk) = zfac * zusd2dt(ji,jj) * tmask(ji,jj,jk)
+                  vsd3dt(ji,jj,jk) = zfac * zvsd2dt(ji,jj) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -207,8 +190,8 @@
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-               usd3d(ji,jj,jk) =   0.5 *  umask(ji,jj,jk)  *        &
-                               &  (usd3dt(ji,jj,jk) + usd3dt(ji+1,jj,jk))
-               vsd3d(ji,jj,jk) =  0.5 *  vmask(ji,jj,jk)  *           &
-                               &  (vsd3dt(ji,jj,jk) + vsd3dt(ji,jj+1,jk))
+                  usd3d(ji,jj,jk) = 0.5 *  umask(ji,jj,jk) *   &
+                                  &  ( usd3dt(ji,jj,jk) + usd3dt(ji+1,jj,jk) )
+                  vsd3d(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk) *   &
+                                  &  ( vsd3dt(ji,jj,jk) + vsd3dt(ji,jj+1,jk) )
                END DO
             END DO
@@ -218,24 +201,26 @@
          CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. )
          !
-         DO jk = 1, jpkm1                       ! e3t * Horizontal divergence
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * fse3u_n(ji  ,jj,jk) * usd3d(ji  ,jj,jk)     &
-                     &                 - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk)     &
-                     &                 + e1v(ji,jj  ) * fse3v_n(ji,jj  ,jk) * vsd3d(ji,jj  ,jk)     &
-                     &                 - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj)
+         DO jk = 1, jpkm1               ! Horizontal divergence
+            DO jj = 2, jpj
+               DO ji = fs_2, jpi   
+                  ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * usd3d(ji  ,jj,jk)     &
+                     &                 - e2u(ji-1,jj) * usd3d(ji-1,jj,jk)     &
+                     &                 + e1v(ji,jj  ) * vsd3d(ji,jj  ,jk)     &
+                     &                 - e1v(ji,jj-1) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj)
                END DO  
             END DO  
-            IF( .NOT. AGRIF_Root() ) THEN
-               IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,jk) = 0._wp      ! east
-               IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,jk) = 0._wp      ! west
-               IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,jk) = 0._wp      ! north
-               IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,jk) = 0._wp      ! south
-            ENDIF
-         END DO
+         END DO
+         !
+         IF( .NOT. AGRIF_Root() ) THEN
+            IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,:) = 0._wp      ! east
+            IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,:) = 0._wp      ! west
+            IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,:) = 0._wp      ! north
+            IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,:) = 0._wp      ! south
+         ENDIF
+         !
          CALL lbc_lnk( ze3hdiv, 'T', 1. ) 
          !
-         DO jk = jpkm1, 1, -1                   !* integrate from the bottom the e3t * hor. divergence
-            wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk)
+         DO jk = jpkm1, 1, -1                   ! integrate from the bottom the e3t * hor. divergence
+            wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - fse3t_n(:,:,jk) * ze3hdiv(:,:,jk)
          END DO
 #if defined key_bdy
@@ -246,8 +231,7 @@
          ENDIF
 #endif
-         CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv )
 
          IF ( ln_zdfqiao ) THEN
-            CALL fld_read( kt, nn_fsbc, sf_wn )      !* read wave parameters from external forcing
+            CALL fld_read( kt, nn_fsbc, sf_wn )      ! read wave parameters from external forcing
             wnum(:,:) = sf_wn(1)%fnow(:,:,1)
            
@@ -256,5 +240,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                   tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2)
+                  tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2 ) 
                END DO
             END DO
Index: /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7339)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7340)
@@ -37,8 +37,5 @@
    USE sbcwave        ! wave module
    USE sbc_oce        ! surface boundary condition: ocean
-
    USE diaptr         ! Poleward heat transport 
-   USE sbcwave        ! wave module
-   USE sbc_oce        ! surface boundary condition: ocean
 
    IMPLICIT NONE
@@ -111,8 +108,8 @@
       !
       !                                         !==  effective transport  ==!
-      IF (ln_wave .AND. ln_sdw) THEN
+      IF( ln_wave .AND. ln_sdw )  THEN
          DO jk = 1, jpkm1
             zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) *      &
-                        &  ( un(:,:,jk) + usd3d(:,:,jk) )                     !eulerian transport + Stokes Drift
+                        &  ( un(:,:,jk) + usd3d(:,:,jk) )                       ! eulerian transport + Stokes Drift
             zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) *      &
                         &  ( vn(:,:,jk) + vsd3d(:,:,jk) )
@@ -121,9 +118,9 @@
          END DO
       ELSE
-      DO jk = 1, jpkm1
-         zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
-         zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk)
-         zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
-      END DO
+         DO jk = 1, jpkm1
+            zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)               ! eulerian transport only
+            zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk)
+            zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
+         END DO
       ENDIF
       !
Index: /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90
===================================================================
--- /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7339)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7340)
@@ -6,10 +6,7 @@
    !! History :  3.6  !  2014-10  (E. Clementi)  Original code
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   qiao_init      
    !!   zdf_qiao        : compute Qiao parameters
    !!----------------------------------------------------------------------
 
-   USE iom             ! I/O manager library
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distribued memory computing library
@@ -18,19 +15,16 @@
    USE sbcwave         ! wave module
    USE dom_oce
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)  
    
-   !!----------------------------------------------------------------------
-   !!   qiao_init       : compute QBv: Qiao terms to be added to vertical eddy
-   !!                     diffusivity and viscosity coefficients 
-   !!----------------------------------------------------------------------
-
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   zdf_qiao    ! routine called in zdf_ric
+   PUBLIC zdf_qiao    ! routine called in step
 
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: QBv, QBvu, QBvv
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: qbv, qbvu, qbvv
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
@@ -45,8 +39,8 @@
       !!                     ***  ROUTINE zdf_qiao ***
       !!
-      !! ** Purpose :Compute the Qiao term (QBv) to be added to
+      !! ** Purpose :Compute the Qiao term (qbv) to be added to
       !!             vertical viscosity and diffusivity coeffs.  
       !!
-      !! ** Method  :QBv = alpha * A * Us(0) * exp (3 * k * z)
+      !! ** Method  :qbv = alpha * A * Us(0) * exp (3 * k * z)
       !!             
       !! ** action  :Compute the Qiao wave dependent term 
@@ -56,68 +50,79 @@
       INTEGER, INTENT( in  ) ::  kt   ! ocean time step
       !
-      INTEGER                ::  jj, ji, jk
+      INTEGER :: jj, ji, jk   ! dummy loop indices
       !!---------------------------------------------------------------------
       !
-      !
-      !                                         ! -------------------- !
       IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
-         ALLOCATE(QBv(jpi,jpj,jpk))             ! -------------------- !
-         ALLOCATE(QBvu(jpi,jpj,jpk))
-         ALLOCATE(QBvv(jpi,jpj,jpk))
+         IF( .NOT. ( ln_wave .AND. ln_sdw ) )   &
+            &   CALL ctl_stop ( 'Ask for wave Qiao enhanced turbulence but ln_wave   &
+            &                    and ln_sdw have to be activated')
+         IF( zdf_qiao_alloc() /= 0 )   &
+            &   CALL ctl_stop( 'STOP', 'zdf_qiao : unable to allocate arrays' )
       ENDIF
 
-      QBv (:,:,:) = 0.0
-      QBvu(:,:,:) = 0.0
-      QBvv(:,:,:) = 0.0
-
       !
-      ! Compute the Qiao term Bv (QBv) to be added to
+      ! Compute the Qiao term Bv (qbv) to be added to
       ! vertical viscosity and diffusivity
-      ! QBv = alpha * A * Us(0) * exp (3 * k * z)
+      ! qbv = alpha * A * Us(0) * exp (3 * k * z)
       ! alpha here is set to 1
       !---------------------------------------------------------------------------------
       !
-      IF ( ln_wave ) THEN
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  QBv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *       &
-               &              exp(3.0 * wnum(ji,jj) *                                &                     
-               &              (-MIN( fsdept(ji  ,jj  ,jk) , fsdept(ji+1,jj  ,jk),    &
-               &                     fsdept(ji  ,jj+1,jk) , fsdept(ji+1,jj+1,jk))))
-               END DO
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               qbv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *           &
+            &                  EXP(3.0 * wnum(ji,jj) *                                &                     
+            &                  (-MIN( fsdepw(ji  ,jj  ,jk), fsdepw(ji+1,jj  ,jk),     &
+            &                         fsdepw(ji  ,jj+1,jk), fsdepw(ji+1,jj+1,jk))))   &
+            &                          * wmask(ji,jj,jk)
             END DO
          END DO
-
-         QBv(jpi,:,:)=QBv(jpim1,:,:)
-         QBv(:,jpj,:)=QBv(:,jpjm1,:)
-
-         !
-         ! Interpolate Qiao parameter QBv into the grid_U and grid_V
-         !-------------------------------------------------
-         !
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  QBvu(ji,jj,jk) = 0.5 *  umask(ji,jj,jk)  *               &
-               &           ( QBv(ji  ,jj,jk) * tmask(ji  ,jj,jk)           &
-               &           + QBv(ji+1,jj,jk) * tmask(ji+1,jj,jk) )
-                  QBvv(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk)  *               &
-               &           ( QBv(ji,jj  ,jk) * tmask(ji,jj  ,jk)           &
-               &           + QBv(ji,jj+1,jk) * tmask(ji,jj+1,jk) )
-               END DO
+      END DO
+      !
+      CALL lbc_lnk( qbv, 'W', 1. )   ! Lateral boundary conditions
+         
+      !
+      ! Interpolate Qiao parameter qbv into the grid_U and grid_V
+      !----------------------------------------------------------
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               qbvu(ji,jj,jk) = 0.5 * wumask(ji,jj,jk)  *              &  
+            &                  ( qbv(ji,jj,jk) + qbv(ji+1,jj  ,jk) )
+               qbvv(ji,jj,jk) = 0.5 * wvmask(ji,jj,jk)  *              &
+            &                  ( qbv(ji,jj,jk) + qbv(ji  ,jj+1,jk) )
             END DO
          END DO
-         ! 
-         QBvu(jpi,:,:)=QBvu(jpim1,:,:)
-         QBvu(:,jpj,:)=QBvu(:,jpjm1,:)
-         QBvv(jpi,:,:)=QBvv(jpim1,:,:)
-         QBvv(:,jpj,:)=QBvv(:,jpjm1,:)
+      END DO
+      ! 
+      CALL lbc_lnk( qbvu, 'U', 1. ) ; CALL lbc_lnk( qbvv, 'V', 1. )   ! Lateral boundary conditions
 
-        ELSE
-           CALL ctl_stop( 'STOP', 'To use Qiao formulation you have to set: ln_wave=.true.')
-        ENDIF
-        !
+      ! Enhance vertical mixing coeff.         
+      !-------------------------------
+      !
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               avmu(ji,jj,jk) = ( avmu(ji,jj,jk) + qbvu(ji,jj,jk) ) * umask(ji,jj,jk)
+               avmv(ji,jj,jk) = ( avmv(ji,jj,jk) + qbvv(ji,jj,jk) ) * vmask(ji,jj,jk)
+               avt (ji,jj,jk) = ( avt (ji,jj,jk) + qbv (ji,jj,jk) ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
    END SUBROUTINE zdf_qiao
+
+   INTEGER FUNCTION zdf_qiao_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_qiao_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( qbv(jpi,jpj,jpk), qbvu(jpi,jpj,jpk), qbvv(jpi,jpj,jpk),   &
+         &      STAT = zdf_qiao_alloc )
+      !
+      IF( lk_mpp             )  CALL mpp_sum ( zdf_qiao_alloc )
+      IF( zdf_qiao_alloc > 0 )  CALL ctl_warn('zdf_qiao_alloc: allocation of arrays failed.')
+      !
+   END FUNCTION zdf_qiao_alloc
       
    !!======================================================================
Index: /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7339)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7340)
@@ -130,26 +130,10 @@
                          CALL zdf_bfr( kstp )         ! bottom friction (if quadratic)
       !                                               ! Vertical eddy viscosity and diffusivity coefficients
-      IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
-      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
-      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
-      IF( ln_zdfqiao )   THEN
-      !Activated Qiao enhanced turbulence but neither ln_wave or ln_sdw are activated
-       IF ( .NOT.( ln_wave .AND. ln_sdw ) )   THEN
-        CALL ctl_stop ( 'Ask for wave Qiao enhanced turbulence but ln_wave and ln_sdw have to be activated')
-       ELSE
-        CALL zdf_qiao(kstp )                             ! Qiao vertical mixing 
-         DO jk = 1, jpkm1
-          DO jj = 1, jpj
-            DO ji = 1, jpi
-               avmu(ji,jj,jk) = (avmu(ji,jj,jk) + QBvu(ji,jj,jk)) * umask(ji,jj,jk)
-               avmv(ji,jj,jk) = (avmv(ji,jj,jk) + QBvv(ji,jj,jk)) * vmask(ji,jj,jk)
-               avt( ji,jj,jk) = (avt( ji,jj,jk) + QBv(ji,jj,jk))  * tmask(ji,jj,jk)
-            END DO
-          END DO
-         END DO
-       ENDIF 
-      ENDIF
-      !
-      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
+      IF( lk_zdfric  )   CALL zdf_ric ( kstp )             ! Richardson number dependent Kz
+      IF( lk_zdftke  )   CALL zdf_tke ( kstp )             ! TKE closure scheme for Kz
+      IF( lk_zdfgls  )   CALL zdf_gls ( kstp )             ! GLS closure scheme for Kz
+      IF( ln_zdfqiao )   CALL zdf_qiao( kstp )             ! Qiao vertical mixing 
+      !
+      IF( lk_zdfcst  ) THEN                                ! Constant Kz (reset avt, avm[uv] to the background value)
          avt (:,:,:) = rn_avt0 * wmask (:,:,:)
          avmu(:,:,:) = rn_avm0 * wumask(:,:,:)
