Index: /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/AMM12/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/AMM12/EXP00/namelist	(revision 3528)
+++ /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/AMM12/EXP00/namelist	(revision 3529)
@@ -144,5 +144,7 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
+   ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave)
    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
+   ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -971,4 +973,7 @@
 !              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
    sn_cdg      =  'cdg_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'  ,''         , ''
+   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   , .true.   , .false. , 'daily'  ,''         , ''
 !
    cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
Index: /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/AMM12_PISCES/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/AMM12_PISCES/EXP00/namelist	(revision 3528)
+++ /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/AMM12_PISCES/EXP00/namelist	(revision 3529)
@@ -144,5 +144,7 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
+   ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave)
    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
+   ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -972,4 +974,7 @@
 !              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
    sn_cdg      =  'cdg_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'  ,''         , ''
+   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   , .true.   , .false. , 'daily'  ,''         , ''
 !
    cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
Index: /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/GYRE/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 3528)
+++ /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 3529)
@@ -144,5 +144,7 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
+   ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave)
    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
+   ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -963,4 +965,7 @@
 !              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
    sn_cdg      =  'cdg_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'  ,''         , ''
+   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   , .true.   , .false. , 'daily'  ,''         , ''
 !
    cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
Index: /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 3528)
+++ /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 3529)
@@ -144,5 +144,7 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
+   ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave)
    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
+   ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -967,4 +969,7 @@
 !              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
    sn_cdg      =  'cdg_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'  ,''         , ''
+   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   , .true.   , .false. , 'daily'  ,''         , ''
 !
    cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
Index: /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 3528)
+++ /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 3529)
@@ -144,5 +144,7 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
-   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave )
+   ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave)
+   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
+   ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -964,4 +966,7 @@
 !              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
    sn_cdg      =  'cdg_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'  ,''         , ''
+   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   , .true.   , .false. , 'daily'  ,''         , ''
 !
    cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
Index: /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 3528)
+++ /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 3529)
@@ -44,5 +44,7 @@
    !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step
    !                                             !:  = 2 annual global mean of e-p-r set to zero
-   LOGICAL , PUBLIC ::   ln_cdgw     = .FALSE.   !: true if neutral drag coefficient read from wave model
+   LOGICAL , PUBLIC ::   ln_wave     = .FALSE.   !: true if some coupling with wave model
+   LOGICAL , PUBLIC ::   ln_cdgw     = .FALSE.   !: true if neutral drag coefficient from wave model
+   LOGICAL , PUBLIC ::   ln_sdw      = .FALSE.   !: true if 3d stokes drift from wave model
 
    !!----------------------------------------------------------------------
Index: /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3528)
+++ /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3529)
@@ -83,5 +83,6 @@
       !!
       NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   &
-         &             ln_blk_mfs, ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr     , nn_fwb, ln_cdgw
+         &             ln_blk_mfs, ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr, nn_fwb, &
+         &             ln_wave, ln_cdgw, ln_sdw
       !!----------------------------------------------------------------------
 
@@ -92,6 +93,8 @@
       ENDIF
 
+      call flush(numout)
       REWIND( numnam )           ! Read Namelist namsbc
       READ  ( numnam, namsbc )
+      call flush(numout)
 
       !                          ! overwrite namelist parameter using CPP key information
@@ -166,7 +169,16 @@
          &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
 
-       !drag coefficient read from wave model definable only with mfs bulk formulae and core 
-       IF(ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )              &
-          &   CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
+      IF ( ln_wave ) THEN
+      !Activated wave module but neither drag nor stokes drift activated
+         IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN
+            CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
+      !drag coefficient read from wave model definable only with mfs bulk formulae and core 
+         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
+             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
+         ENDIF
+      ELSE
+      IF ( ln_cdgw .OR. ln_sdw  )                                         & 
+         &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ')
+      ENDIF 
       
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
@@ -249,5 +261,5 @@
       !                                                  ! averaged over nf_sbc time-step
 
-      IF (ln_cdgw) CALL sbc_wave( kt )
+      IF (ln_wave) CALL sbc_wave( kt )
                                                    !==  sbc formulation  ==!
                                                             
Index: /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 3528)
+++ /branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 3529)
@@ -4,5 +4,6 @@
    !! Wave module 
    !!======================================================================
-   !! History :  3.3.1  !   2011-09  (Adani M)  Original code
+   !! History :  3.3.1  !   2011-09  (Adani M)  Original code: Drag Coefficient 
+   !!         :  3.4    !   2012-10  (Adani M)                 Stokes Drift 
    !!----------------------------------------------------------------------
    USE iom             ! I/O manager library
@@ -22,7 +23,16 @@
    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
    
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wave	  ! structure of input fields (file informations, fields read)
+   INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes 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_wn  = 3           ! index of wave number                 (1/m)    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
    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
+   REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum 
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d 
 
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
@@ -40,13 +50,26 @@
       !! ** Method  : - Read namelist namsbc_wave
       !!              - Read Cd_n10 fields in netcdf files 
+      !!              - Read stokes drift 2d in netcdf files 
+      !!              - Read wave number      in netcdf files 
+      !!              - Compute 3d stokes drift using monochromatic
       !! ** action  :   
       !!               
       !!---------------------------------------------------------------------
-      INTEGER, INTENT( in  ) ::  kt   ! ocean time step
+      USE oce,  ONLY : un,vn,hdivn,rotn
+      USE divcur
+      USE wrk_nemo
+#if defined key_bdy
+      USE bdy_oce, ONLY : bdytmask
+#endif
+      INTEGER, INTENT( in  ) ::  kt       ! ocean time step
       INTEGER                ::  ierror   ! return error code
-      CHARACTER(len=100)     ::  cn_dir_cdg                       ! Root directory for location of drag coefficient files
-      TYPE(FLD_N)            ::  sn_cdg                          ! informations about the fields to be read
+      INTEGER                ::  ifpr, jj,ji,jk 
+      REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy
+      REAL                                          ::  z2dt,z1_2dt
+      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_wn   ! informations about the fields to be read
       !!---------------------------------------------------------------------
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir_cdg
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn
       !!---------------------------------------------------------------------
 
@@ -62,5 +85,8 @@
          !              !   name   !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
          sn_cdg = FLD_N('cdg_wave'  ,    1     ,'drag_coeff',  .true.    , .false. ,   'daily'   , ''       , ''       )
-         cn_dir_cdg = './'          ! directory in which the Patm data are 
+         sn_usd = FLD_N('sdw_wave'  ,    1     ,'u_sd2d',      .true.    , .false. ,   'daily'   , ''       , ''       )
+         sn_vsd = FLD_N('sdw_wave'  ,    1     ,'v_sd2d',      .true.    , .false. ,   'daily'   , ''       , ''       )
+         sn_wn = FLD_N( 'sdw_wave'  ,    1     ,'wave_num',    .true.    , .false. ,   'daily'   , ''       , ''       )
+         cn_dir = './'          ! directory in which the wave data are 
          
 
@@ -69,18 +95,100 @@
          !
 
-         ALLOCATE( sf_wave(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' )
-         !
-         CALL fld_fill( sf_wave, (/ sn_cdg /), cn_dir_cdg, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
-                                ALLOCATE( sf_wave(1)%fnow(jpi,jpj,1)   )
-         IF( sn_cdg%ln_tint )   ALLOCATE( sf_wave(1)%fdta(jpi,jpj,1,2) )
-         ALLOCATE( cdn_wave(jpi,jpj) )
-         cdn_wave(:,:) = 0.0
+         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) )
+            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_sdw ) THEN
+            slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn
+            ALLOCATE( sf_sd(3), 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' )
+            !
+            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) )
+            END DO
+            CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
+            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) )
+            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
+            usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0
+            usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0
+         ENDIF
       ENDIF
          !
          !
-      CALL fld_read( kt, nn_fsbc, sf_wave )      !* read drag coefficient from external forcing
-      cdn_wave(:,:) = sf_wave(1)%fnow(:,:,1)
+      IF ( ln_cdgw ) THEN
+         CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing
+         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
+      ENDIF
+      IF ( ln_sdw )  THEN
+          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing
 
+         ! Interpolate wavenumber, stokes drift into the grid_V and grid_V
+         !-------------------------------------------------
+
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
+               &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
+
+               vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
+               &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
+
+               usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) &
+               &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
+
+               vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) &
+               &                                + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
+            END DO
+         END DO
+
+          !Computation of the 3d Stokes Drift
+          DO jk = 1, jpk
+             DO jj = 1, jpj-1
+                DO ji = 1, jpi-1
+                   usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk))))
+                   vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk))))
+                END DO
+             END DO
+             usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept(jpi,:,jk)) )
+             vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept(:,jpj,jk)) )
+          END DO
+
+          CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
+          
+          udummy(:,:,:)=un(:,:,:)
+          vdummy(:,:,:)=vn(:,:,:)
+          hdivdummy(:,:,:)=hdivn(:,:,:)
+          rotdummy(:,:,:)=rotn(:,:,:)
+          un(:,:,:)=usd3d(:,:,:)
+          vn(:,:,:)=vsd3d(:,:,:)
+          CALL div_cur(kt)
+      !                                           !------------------------------!
+      !                                           !     Now Vertical Velocity    !
+      !                                           !------------------------------!
+          z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
+
+          z1_2dt = 1.e0 / z2dt
+          DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
+             ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
+             wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
+                &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
+                &                         * tmask(:,:,jk) * z1_2dt
+#if defined key_bdy
+             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
+#endif
+          END DO
+          hdivn(:,:,:)=hdivdummy(:,:,:)
+          rotn(:,:,:)=rotdummy(:,:,:)
+          vn(:,:,:)=vdummy(:,:,:)
+          un(:,:,:)=udummy(:,:,:)
+          CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
+      ENDIF
    END SUBROUTINE sbc_wave
       
