Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 9098)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 9115)
@@ -45,5 +45,5 @@
    LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_tauwoc  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tauwoc = .FALSE.
    LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE.
@@ -59,5 +59,5 @@
    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_tauwoc   ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauwoc  ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauw    ! structure of input fields (file informations, fields read) ocean stress components from wave model
 
@@ -100,16 +100,15 @@
       REAL(wp) ::  zstokes_psi_u_bot, zstokes_psi_v_bot
       REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v
-      REAL(wp), DIMENSION(:,:)  , POINTER ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd   ! 2D workspace
-      REAL(wp), DIMENSION(:,:)  , POINTER ::   zstokes_psi_u_top, zstokes_psi_v_top   ! 2D workspace
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   ze3divh                            ! 3D workspace
-      !!---------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   ze3divh )
-      CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
-      CALL wrk_alloc( jpi,jpj,       zstokes_psi_u_top, zstokes_psi_v_top)
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd     ! 2D workspace
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zstokes_psi_u_top, zstokes_psi_v_top ! 2D workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ze3divh                              ! 3D workspace
+      !!---------------------------------------------------------------------
+      !
+      ALLOCATE( ze3divh(jpi,jpj,jpk) )
+      ALLOCATE( zk_t(jpi,jpj), zk_u(jpi,jpj), zk_v(jpi,jpj), zu0_sd(jpi,jpj), zv0_sd(jpi,jpj) )
       !
       ! select parameterization for the calculation of vertical Stokes drift
       ! exp. wave number at t-point
-      IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN   ! (Eq. (19) in Breivick et al. (2014) )
+      IF( ll_st_bv_li ) THEN   ! (Eq. (19) in Breivik et al. (2014) )
          zfac = 2.0_wp * rpi / 16.0_wp
          DO jj = 1, jpj
@@ -132,5 +131,5 @@
             END DO
          END DO
-      ELSE IF( nn_sdrift==jp_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
+      ELSE IF( ll_st_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
@@ -145,5 +144,5 @@
       !
       !                       !==  horizontal Stokes Drift 3D velocity  ==!
-      IF( nn_sdrift==jp_breivik ) THEN
+      IF( ll_st_bv2014 ) THEN
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -163,16 +162,42 @@
             END DO
          END DO
-      ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakfr ) THEN
+      ELSE IF( ll_st_li2017 .OR. ll_st_peakfr ) THEN
+         ALLOCATE( zstokes_psi_u_top(jpi,jpj), zstokes_psi_v_top(jpi,jpj) )
+         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
+            DO ji = 1, jpim1
+               zstokes_psi_u_top(ji,jj) = 0._wp
+               zstokes_psi_v_top(ji,jj) = 0._wp
+            END DO
+         END DO
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
-                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
-                  !                          
-                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
-                  zkh_v = zk_v(ji,jj) * zdep_v
-                  !                                ! Depth attenuation
-                  zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u))
-                  zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v))
+                  zbot_u = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji+1,jj,jk+1) )
+                  zbot_v = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji,jj+1,jk+1) )
+                  zkb_u  = 2.0_wp * zk_u(ji,jj) * zbot_u                             ! 2k * bottom depth
+                  zkb_v  = 2.0_wp * zk_v(ji,jj) * zbot_v                             ! 2k * bottom depth
+                  !
+                  zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u_n(ji,jj,jk))     ! 2k * thickness
+                  zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v_n(ji,jj,jk))     ! 2k * thickness
+
+                  ! Depth attenuation .... do u component first..
+                  zdepth      = zkb_u
+                  zsqrt_depth = SQRT(zdepth)
+                  zexp_depth  = EXP(-zdepth)
+                  zstokes_psi_u_bot = 1.0_wp - zexp_depth  &
+                       &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
+                       &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
+                  zda_u                    = ( zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj) ) / zke3_u
+                  zstokes_psi_u_top(ji,jj) =   zstokes_psi_u_bot
+
+                  !         ... and then v component
+                  zdepth      =zkb_v
+                  zsqrt_depth = SQRT(zdepth)
+                  zexp_depth  = EXP(-zdepth)
+                  zstokes_psi_v_bot = 1.0_wp - zexp_depth  &
+                       &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
+                       &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
+                  zda_v                    = ( zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj) ) / zke3_v
+                  zstokes_psi_v_top(ji,jj) =   zstokes_psi_v_bot
                   !
                   usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
@@ -181,4 +206,5 @@
             END DO
          END DO
+         DEALLOCATE( zstokes_psi_u_top, zstokes_psi_v_top )
       ENDIF
 
@@ -232,7 +258,6 @@
       CALL iom_put( "wstokes",  wsd  )
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   ze3divh )
-      CALL wrk_dealloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
-      CALL wrk_dealloc( jpi,jpj,       zstokes_psi_u_top, zstokes_psi_v_top)
+      DEALLOCATE( ze3divh )
+      DEALLOCATE( zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
       !
    END SUBROUTINE sbc_stokes
@@ -298,6 +323,6 @@
       ENDIF
 
-      IF( ln_tauwoc .AND. .NOT. cpl_tauwoc ) THEN    !==  Wave induced stress  ==!
-         CALL fld_read( kt, nn_fsbc, sf_tauwoc )          ! read wave norm stress from external forcing
+      IF( ln_tauwoc .AND. .NOT. cpl_tauwoc ) THEN  !==  Wave induced stress  ==!
+         CALL fld_read( kt, nn_fsbc, sf_tauwoc )         ! read wave norm stress from external forcing
          tauoc_wave(:,:) = sf_tauwoc(1)%fnow(:,:,1)
       ENDIF
@@ -326,11 +351,9 @@
          ENDIF
            
-         !                                         !==  Computation of the 3d Stokes Drift  ==! 
+         ! Calculate only if required fields have been read
+         ! In coupled wave model-NEMO case the call is done after coupling
          !
-         IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. &
-                          jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. &
-             (nn_sdrift==jp_peakfr .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) &
-            CALL sbc_stokes()            ! Calculate only if required fields are read
-         !                               ! In coupled wave model-NEMO case the call is done after coupling
+         IF( ( ll_st_bv_li   .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0 ) .OR. &
+           & ( ll_st_peakfr  .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0                ) ) CALL sbc_stokes()
          !
       ENDIF
@@ -356,9 +379,9 @@
       INTEGER ::   ifpr
       !!
-      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
+      CHARACTER(len=100)     ::  cn_dir                            ! Root directory for location of drag coefficient files
       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i, slf_j     ! array of namelist informations on the fields to read
       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
                              &   sn_hsw, sn_wmp, sn_wfr, sn_wnum, &
-                             &   sn_tauwoc, sn_tauwx, sn_tauwy      ! informations about the fields to be read
+                             &   sn_tauwoc, sn_tauwx, sn_tauwy     ! informations about the fields to be read
       !
       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, &
@@ -377,5 +400,5 @@
       IF( ln_cdgw ) THEN
          IF( .NOT. cpl_wdrag ) THEN
-            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
+            ALLOCATE( sf_cd(1), STAT=ierror )               !* allocate and fill sf_wave with sn_cdg
             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
             !
@@ -392,5 +415,5 @@
             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
             !
-                                    ALLOCATE( sf_tauwoc(1)%fnow(jpi,jpj,1)   )
+                                     ALLOCATE( sf_tauwoc(1)%fnow(jpi,jpj,1)   )
             IF( sn_tauwoc%ln_tint )  ALLOCATE( sf_tauwoc(1)%fdta(jpi,jpj,1,2) )
             CALL fld_fill( sf_tauwoc, (/ sn_tauwoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
@@ -428,13 +451,13 @@
             jp_vsd = jpfld
          ENDIF
-         IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
+         IF( .NOT. cpl_hsig  .AND. ll_st_bv_li  ) THEN
             jpfld  = jpfld + 1
             jp_hsw = jpfld
          ENDIF
-         IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
+         IF( .NOT. cpl_wper  .AND. ll_st_bv_li  ) THEN
             jpfld  = jpfld + 1
             jp_wmp = jpfld
          ENDIF
-         IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakfr ) THEN
+         IF( .NOT. cpl_wfreq .AND. ll_st_peakfr ) THEN
             jpfld  = jpfld + 1
             jp_wfr = jpfld
