Index: branches/2014/dev_r4822_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/2014/dev_r4822_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 4823)
+++ branches/2014/dev_r4822_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 4824)
@@ -4,15 +4,16 @@
    !! Wave module 
    !!======================================================================
-   !! History :  3.3.1  !   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 
+   !!         :  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
    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 oce
-   USE sbc_oce	       ! Surface boundary condition: ocean fields
+   USE sbc_oce         ! Surface boundary condition: ocean fields
    USE domvvl
-
+   USE phycst
    
    !!----------------------------------------------------------------------
@@ -25,12 +26,17 @@
    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
    
-   INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift
+   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_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
+   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
    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
-   REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum 
+   REAL(wp),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 
 
@@ -54,5 +60,6 @@
       !!              - Read stokes drift 2d in netcdf files 
       !!              - Read wave number      in netcdf files 
-      !!              - Compute 3d stokes drift using monochromatic
+      !!              - Compute 3d stokes drift using Breivik et al.,2014
+      !!                formulation
       !! ** action  :   
       !!               
@@ -61,4 +68,5 @@
       USE divcur
       USE wrk_nemo
+      USE zdf_oce,  ONLY : ln_zdfqiao
 #if defined key_bdy
       USE bdy_oce, ONLY : bdytmask
@@ -69,10 +77,13 @@
       INTEGER                ::   ios     ! Local integer output status for namelist read
       REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy
-      REAL                                          ::  z2dt,z1_2dt
+      REAL(wp)                                      ::  z2dt,z1_2dt
+      REAL(wp)                                      ::  ztransp,zsp0, zk, zus,zvs
+      REAL(wp), DIMENSION(jpi,jpj)                  ::  zfac 
       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, sn_usd, sn_vsd, sn_wn
+      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
+                             &   sn_swh, sn_wmp, sn_wnum         ! informations about the fields to be read
+      !!---------------------------------------------------------------------
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum
       !!---------------------------------------------------------------------
 
@@ -103,7 +114,9 @@
             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
+            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
             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
             !
@@ -112,12 +125,25 @@
                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( usd2d(jpi,jpj),vsd2d(jpi,jpj),usd2dt(jpi,jpj),vsd2dt(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
+            ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) )
+            usd2d(:,:)   = 0.0 ;  vsd2d(:,:)   = 0.0 ;
+            usd2dt(:,:)  = 0.0 ;  vsd2dt(:,:)  = 0.0 ;
+            usd3d(:,:,:) = 0.0 ;  vsd3d(:,:,:) = 0.0 ;
+            wsd3d(:,:,:) = 0.0 ;
+            swh(:,:)     = 0.0 ;  wmp(:,:)     = 0.0 ;
+            IF ( ln_zdfqiao ) THEN
+               ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fillsf_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( 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.0; tsd2d(:,:) = 0.0
+            ENDIF
          ENDIF
       ENDIF
-         !
          !
       IF ( ln_cdgw ) THEN
@@ -126,46 +152,78 @@
       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
+          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing
+          swh(:,:)    = sf_sd(jp_swh)%fnow(:,:,1)
+          wmp(:,:)    = sf_sd(jp_wmp)%fnow(:,:,1)
+          usd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1) 
+          vsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)
          !-------------------------------------------------
-
+         ! Interpolate 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) )
+               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) )
             END DO
          END DO
-
-          !Computation of the 3d Stokes Drift
+          !Computation of the 3d Stokes Drift according to Breivik et al.,2014
+          !(DOI: 10.1175/JPO-D-14-0020.1)
           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( sf_sd(jp_usd)%fnow(ji,jj,1)**2 +  sf_sd(jp_vsd)%fnow(ji,jj,1)**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-1
                 DO ji = 1, jpi-1
-                   usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk))))
-                   vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk))))
+                 ! 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) &
+                 &                                + zfac(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
-             usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) )
-             vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,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(:,:,:)
+      !-------------------------------------------------------------------          
+      ! Store horiz. velocity divergence and Rot in TMP array 
+      ! -------------------------------------------------------------------
+
+          udummy(:,:,:)     = un(:,:,:)
+          vdummy(:,:,:)     = vn(:,:,:)
+          hdivdummy(:,:,:)  = hdivn(:,:,:)
+          rotdummy(:,:,:)   = rotn(:,:,:)
+          un(:,:,:)         = usd3d(:,:,:)
+          vn(:,:,:)         = vsd3d(:,:,:)
+
+      ! Compute divergence using 3d stokes drift to calculate vertical SD
+      ! velocity
+      ! -------------------------------------------------------------------
+
           CALL div_cur(kt)
+
+      ! -------------------------------------------------------------------
       !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
@@ -177,16 +235,33 @@
              ! - 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
+                &                            - ( 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(:,:,:)
+      ! Back to state variables  
+      ! -------------------------------------------------------------------
+          hdivn(:,:,:)  = hdivdummy(:,:,:)
+          rotn(:,:,:)   = rotdummy(:,:,:)
+          vn(:,:,:)     = vdummy(:,:,:)
+          un(:,:,:)     = udummy(:,:,:)
+
           CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
+!
+        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
       ENDIF
+
+
    END SUBROUTINE sbc_wave
       
