Index: /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6252)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6253)
@@ -1240,5 +1240,5 @@
 !              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
 !              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_cdg      =  'sdw_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
@@ -1248,5 +1248,5 @@
    sn_tauoc    -  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
 !
-   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+   cn_dir  = './'  !  root directory for the location of drag coefficient files
 /
 !-----------------------------------------------------------------------
Index: /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 6252)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 6253)
@@ -423,20 +423,4 @@
       ENDIF  
       !
-      ! Add Stokes Coriolis if defined
-      IF ( ln_stcor ) THEN
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-
-               zy1 = ff(ji  ,jj-1) * ( vsd2d(ji  ,jj-1) + vsd2d(ji+1,jj-1) )
-               zy2 = ff(ji  ,jj  ) * ( vsd2d(ji  ,jj  ) + vsd2d(ji+1,jj  ) )
-               zx1 = ff(ji-1,jj  ) * ( usd2d(ji-1,jj  ) + usd2d(ji-1,jj+1) )
-               zx2 = ff(ji  ,jj  ) * ( usd2d(ji  ,jj  ) + usd2d(ji  ,jj+1) )
-
-               zu_frc(ji,jj) = zu_frc(ji,jj) + 0.25 * ( zy1 + zy2 ) * hur(ji,jj)
-               zv_frc(ji,jj) = zv_frc(ji,jj) - 0.25 * ( zx1 + zx2 ) * hvr(ji,jj)
-           END DO
-         END DO
-      ENDIF
-      !
       IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
          IF (ln_bt_fw) THEN
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 6252)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 6253)
@@ -39,5 +39,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
-   USE sbcwave,  ONLY: usd2dt, vsd2dt,wsd3d
+   USE sbcwave,  ONLY:  zusd2dt, zvsd2dt,wsd3d
 
    IMPLICIT NONE
@@ -240,6 +240,6 @@
        ! Compute the surface vertical velocity accounting for the Stokes Drift 
        !---------------------------------------------------------------------
-       wn(:,:,1) = wn(:,:,1) + usd2dt(:,:) * dsshnu(:,:)     &
-                 &           + vsd2dt(:,:) * dsshnv(:,:)      &
+       wn(:,:,1) = wn(:,:,1) + zusd2dt(:,:) * dsshnu(:,:)     &
+                 &           + zvsd2dt(:,:) * dsshnv(:,:)      &
                  &           - ( wsd3d (:,:,1) ) * tmask(:,:,1)
       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 6252)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 6253)
@@ -6,9 +6,9 @@
    !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient 
    !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift 
-   !! History :  3.6  !2014-09  (Clementi E, Oddo P)New Stokes Drift Computation
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   sbc_wave      : read drag coefficient from wave model in netcdf files 
+   !!            3.6  !   2014-09  (Clementi E, Oddo P)New Stokes Drift Computation
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_wave      : wave data from wave model in netcdf files 
    !!----------------------------------------------------------------------
    USE oce            ! 
@@ -27,5 +27,5 @@
    PRIVATE
 
-   PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
+   PUBLIC   sbc_wave    ! routine called in sbcmod
    
    INTEGER , PARAMETER ::   jpfld  = 4           ! number of files to read for stokes drift
@@ -34,15 +34,16 @@
    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
    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift
    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao
    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: usd2d,vsd2d
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: swh,wmp,wnum
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: usd2dt,vsd2dt,tsd2d
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d 
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: tauoc_wave
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: cdn_wave 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: swh,wmp, wnum
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tauoc_wave
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tsd2d
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: usd2d, vsd2d
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt
 
    !! * Substitutions
@@ -58,7 +59,7 @@
    SUBROUTINE sbc_wave( kt )
       !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE sbc_apr  ***
+      !!                     ***  ROUTINE sbc_wave  ***
       !!
-      !! ** Purpose :   read drag coefficient from wave model  in netcdf files.
+      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
       !!
       !! ** Method  : - Read namelist namsbc_wave
@@ -68,5 +69,5 @@
       !!              - Compute 3d stokes drift using Breivik et al.,2014
       !!                formulation
-      !! ** action  :   
+      !! ** action  
       !!---------------------------------------------------------------------
       USE zdf_oce,  ONLY : ln_zdfqiao
@@ -77,13 +78,12 @@
       INTEGER                ::   ifpr, jj,ji,jk 
       INTEGER                ::   ios      ! Local integer output status for namelist read
-
-      REAL(wp)                       ::  ztransp,zsp0, zk, zus,zvs
+      !
+      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
-      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
                              &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   zusd_t, zvsd_t, ze3hdiv   ! 3D workspace
       !!
       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc
@@ -112,5 +112,5 @@
             cdn_wave(:,:) = 0.0
          ENDIF
-!
+
          IF ( ln_tauoc ) THEN
             ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc
@@ -123,5 +123,5 @@
             tauoc_wave(:,:) = 0.0
         ENDIF
-!
+
          IF ( ln_sdw ) THEN
             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd;
@@ -134,11 +134,11 @@
                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
             END DO
-!
+
             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),usd2dt(jpi,jpj),vsd2dt(jpi,jpj))
+            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj) )
             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
             ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) )
-            usd3d(:,:,:) = 0._wp   ;   usd2d(:,:) = 0._wp   ;    usd2dt(:,:) = 0._wp ;
-            vsd3d(:,:,:) = 0._wp   ;   vsd2d(:,:) = 0._wp   ;    vsd2dt(:,:) = 0._wp ;
+            usd3d(:,:,:) = 0._wp   ;   usd2d(:,:) = 0._wp   ;   
+            vsd3d(:,:,:) = 0._wp   ;   vsd2d(:,:) = 0._wp   ; 
             wsd3d(:,:,:) = 0._wp   ;
             swh  (:,:)   = 0._wp   ;    wmp (:,:) = 0._wp ;
@@ -159,17 +159,17 @@
          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
          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
-         swh(:,:)    = sf_sd(jp_swh)%fnow(:,:,1)  ! significant wave height
-         wmp(:,:)    = sf_sd(jp_wmp)%fnow(:,:,1)  ! wave mean period
-         usd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)  ! 2D zonal Stokes Drift
-         vsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)  ! 2D meridional Stokes Drift
+         swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height
+         wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
+         zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
+         zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
          !
          !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014
@@ -177,7 +177,7 @@
          !
          CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv )
-          DO jk = 1, jpk
-             DO jj = 1, jpj
-                DO ji = 1, jpi
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
                ! On T grid
                ! Stokes transport speed estimated from Hs and Tmean
@@ -189,32 +189,32 @@
                ! 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, jpjm1
-                DO ji = 1, jpim1
-                 ! Into the U and V Grid 
-                 zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) &
+               END DO
+            END DO
+
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  ! Into the U and V Grid 
+                  zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) &
                  &                                + zfac(ji+1,jj) * tmask(ji+1,jj,1) )
-!
-                 zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) &
+
+                  zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) &
                  &                                + zfac(ji,jj+1) * tmask(ji,jj+1,1) )
-!
-                 usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( usd2dt(ji,jj) * tmask(ji,jj,1) &
-                 &                                             +  usd2dt(ji+1,jj) * tmask(ji+1,jj,1) )
-!
-                 vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vsd2dt(ji,jj) * tmask(ji,jj,1) &
-                 &                                              + vsd2dt(ji,jj+1) * tmask(ji,jj+1,1) )
-!
-                 usd3d(ji,jj,jk) = usd2d(ji,jj)*zus
-                 vsd3d(ji,jj,jk) = vsd2d(ji,jj)*zvs
-                END DO
-             END DO
-          END DO
+
+                  usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zusd2dt(ji,jj) * tmask(ji,jj,1) &
+                 &                                             +  zusd2dt(ji+1,jj) * tmask(ji+1,jj,1) )
+
+                  vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zvsd2dt(ji,jj) * tmask(ji,jj,1) &
+                 &                                              + zvsd2dt(ji,jj+1) * tmask(ji,jj+1,1) )
+
+                  usd3d(ji,jj,jk) = usd2d(ji,jj) * zus
+                  vsd3d(ji,jj,jk) = vsd2d(ji,jj) * zvs
+               END DO
+            END DO
+         END DO
          !
          CALL lbc_lnk( usd3d(:,:,:), 'U', -1. )
          CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. )
          !
-         DO jk = 1, jpkm1                       !* e3t * Horizontal divergence  ==!
+         DO jk = 1, jpkm1                       ! e3t * Horizontal divergence
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -244,19 +244,19 @@
          ENDIF
 #endif
-!         CALL wrk_dealloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv )
          CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv )
-!
-        IF ( ln_zdfqiao )  THEN
-          wnum(:,:) = sf_wn(1)%fnow(:,:,1)
-         ! Calculate the module of the stokes drift on T grid
-         !-------------------------------------------------
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-                tsd2d(ji,jj) = ((sf_sd(jp_usd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0  +     &
+
+         IF ( ln_zdfqiao ) THEN
+            wnum(:,:) = sf_wn(1)%fnow(:,:,1)
+           
+            ! Calculate the module of the stokes drift on T grid
+            !-------------------------------------------------
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                   tsd2d(ji,jj) = ((sf_sd(jp_usd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0  +     &
                 &               (sf_sd(jp_vsd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0)**0.5
-            END DO
-         END DO
-        ENDIF
-      !
+               END DO
+            END DO
+         ENDIF
+         !
       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 6252)
+++ /branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 6253)
@@ -53,22 +53,20 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT( in  ) ::  kt   ! ocean time step
-      INTEGER                ::  ierror   ! return error code
-      INTEGER                ::  jj,ji,jk
+      !
+      INTEGER                ::  jj, ji, jk
       !!---------------------------------------------------------------------
-      !!----------------------------------------------------------------------
       !
       !
       !                                         ! -------------------- !
       IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
-      !                                         ! -------------------- !
-         ALLOCATE(QBv(jpi,jpj,jpk))
+         ALLOCATE(QBv(jpi,jpj,jpk))             ! -------------------- !
          ALLOCATE(QBvu(jpi,jpj,jpk))
          ALLOCATE(QBvv(jpi,jpj,jpk))
-
       ENDIF
 
-         QBv (:,:,:) = 0.0
-         QBvu(:,:,:) = 0.0
-         QBvv(:,:,:) = 0.0
+      QBv (:,:,:) = 0.0
+      QBvu(:,:,:) = 0.0
+      QBvv(:,:,:) = 0.0
+
       !
       ! Compute the Qiao term Bv (QBv) to be added to
@@ -80,15 +78,13 @@
       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( gdept_0(ji  ,jj  ,jk) , gdept_0(ji+1,jj  ,jk),     & 
+            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( gdept_0(ji  ,jj  ,jk) , gdept_0(ji+1,jj  ,jk),  & 
                &               gdept_0(ji  ,jj+1,jk),gdept_0(ji+1,jj+1,jk))))
-
-            ENDDO
-          ENDDO
-         ENDDO
+               END DO
+            END DO
+         END DO
 
          QBv(jpi,:,:)=QBv(jpim1,:,:)
@@ -98,31 +94,28 @@
          ! 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 * ( 2. - umask(ji,jj,jk) ) *                 &
-               &           ( QBv(ji,jj,jk) * tmask(ji,jj,jk)                   & 
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  QBvu(ji,jj,jk)=0.5 * ( 2. - 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 * ( 2. - vmask(ji,jj,jk) ) *                 &
-               &           ( QBv(ji,jj,jk) * tmask(ji,jj,jk)                   &
+                  QBvv(ji,jj,jk)=0.5 * ( 2. - vmask(ji,jj,jk) ) *                 &
+               &           ( QBv(ji,jj,jk) * tmask(ji,jj,jk)                      &
                &           + QBv(ji,jj+1,jk) * tmask(ji,jj+1,jk) )
 
-             ENDDO
-            ENDDO
-          ENDDO
-!
+               END DO
+            END DO
+         END DO
+         ! 
          QBvu(jpi,:,:)=QBvu(jpim1,:,:)
          QBvu(:,jpj,:)=QBvu(:,jpjm1,:)
          QBvv(jpi,:,:)=QBvv(jpim1,:,:)
          QBvv(:,jpj,:)=QBvv(:,jpjm1,:)
-!
+
         ELSE
+           CALL ctl_stop( 'STOP', 'To use Qiao formulation you have to set: ln_wave=.true.')
+        ENDIF
         !
-        CALL ctl_stop( 'STOP', 'To use Qiao formulation you have to set: ln_wave=.true.')
-        !
-        ENDIF
-!
    END SUBROUTINE zdf_qiao
       
