Index: /branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7167)
+++ /branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7168)
@@ -23,4 +23,5 @@
    USE sbcapr
    USE sbcdcy          ! surface boundary condition: diurnal cycle
+   USE sbcwave         ! surface boundary condition: waves
    USE phycst          ! physical constants
 #if defined key_lim3
@@ -105,5 +106,14 @@
    INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness 
    INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received
+   INTEGER, PARAMETER ::   jpr_mslp   = 43            ! mean sea level pressure 
+   INTEGER, PARAMETER ::   jpr_hsig   = 44            ! Hsig 
+   INTEGER, PARAMETER ::   jpr_phioc  = 45            ! Wave=>ocean energy flux 
+   INTEGER, PARAMETER ::   jpr_sdrftx = 46            ! Stokes drift on grid 1 
+   INTEGER, PARAMETER ::   jpr_sdrfty = 47            ! Stokes drift on grid 2 
+   INTEGER, PARAMETER ::   jpr_wper   = 48            ! Mean wave period
+   INTEGER, PARAMETER ::   jpr_wnum   = 49            ! Mean wavenumber
+   INTEGER, PARAMETER ::   jpr_wstrf  = 50            ! Stress fraction adsorbed by waves
+   INTEGER, PARAMETER ::   jpr_wdrag  = 51            ! Neutral surface drag coefficient
+   INTEGER, PARAMETER ::   jprcv      = 51            ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
@@ -135,5 +145,9 @@
    INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl)
    INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jpsnd      = 28            ! total number of fields sended
+   INTEGER, PARAMETER ::   jps_ficet  = 29            ! total ice fraction  
+   INTEGER, PARAMETER ::   jps_ocxw   = 30            ! currents on grid 1  
+   INTEGER, PARAMETER ::   jps_ocyw   = 31            ! currents on grid 2
+   INTEGER, PARAMETER ::   jps_wlev   = 32            ! water level 
+   INTEGER, PARAMETER ::   jpsnd      = 32            ! total number of fields sent 
 
    !                                                         !!** namelist namsbc_cpl **
@@ -149,5 +163,9 @@
    ! Received from the atmosphere                     !
    TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
-   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                        
+   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp                           
+   ! Send to waves 
+   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 
+   ! Received from waves 
+   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_wstrf,sn_rcv_wdrag
    ! Other namelist parameters                        !
    INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -161,4 +179,7 @@
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
+
+   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2] 
+   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0) 
 
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
@@ -179,5 +200,5 @@
       !!             ***  FUNCTION sbc_cpl_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(3)
+      INTEGER :: ierr(4)
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -190,4 +211,6 @@
       ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
       !
+      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) ) 
+
       sbc_cpl_alloc = MAXVAL( ierr )
       IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
@@ -216,8 +239,10 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
       !!
-      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      &
-         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      &
-         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   &
-         &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask
+      NAMELIST/namsbc_cpl/  sn_snd_temp , sn_snd_alb  , sn_snd_thick , sn_snd_crt   , sn_snd_co2,      & 
+         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
+         &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   & 
+         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wstrf ,   &
+         &                  sn_rcv_wdrag, sn_rcv_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   &
+         &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp 
       !!---------------------------------------------------------------------
       !
@@ -260,8 +285,17 @@
          WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
+         WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')' 
+         WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')' 
+         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')' 
+         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')' 
+         WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')' 
+         WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')' 
+         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')' 
+         WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
          WRITE(numout,*)'  sent fields (multiple ice categories)'
          WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
          WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
          WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
+         WRITE(numout,*)'      total ice fraction              = ', TRIM(sn_snd_ifrac%cldes ), ' (', TRIM(sn_snd_ifrac%clcat ), ')' 
          WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
          WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
@@ -269,4 +303,10 @@
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
          WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
+         WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')' 
+         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')' 
+         WRITE(numout,*)'      surface current to waves        = ', TRIM(sn_snd_crtw%cldes   ), ' (', TRIM(sn_snd_crtw%clcat   ), ')' 
+         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref 
+         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor 
+         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd 
          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
@@ -470,4 +510,10 @@
       !                                                      ! ------------------------- !
       srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
+
+      !                                                      ! ------------------------- ! 
+      !                                                      ! Mean Sea Level Pressure   ! 
+      !                                                      ! ------------------------- ! 
+      srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE. 
+
       !                                                      ! ------------------------- !
       !                                                      !   topmelt and botmelt     !   
@@ -483,4 +529,48 @@
          srcv(jpr_topm:jpr_botm)%laction = .TRUE.
       ENDIF
+      !                                                      ! ------------------------- !
+      !                                                      !      Wave breaking        !    
+      !                                                      ! ------------------------- ! 
+      srcv(jpr_hsig)%clname  = 'O_Hsigwa'    ! significant wave height
+      IF( TRIM(sn_rcv_hsig%cldes  ) == 'coupled' )  THEN
+         srcv(jpr_hsig)%laction = .TRUE.
+         cpl_hsig = .TRUE.
+      ENDIF
+      srcv(jpr_phioc)%clname = 'O_PhiOce'    ! wave to ocean energy
+      IF( TRIM(sn_rcv_phioc%cldes ) == 'coupled' )  THEN
+         srcv(jpr_phioc)%laction = .TRUE.
+         cpl_phioc = .TRUE.
+      ENDIF
+      srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction
+      IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' )  THEN
+         srcv(jpr_sdrftx)%laction = .TRUE.
+         cpl_sdrftx = .TRUE.
+      ENDIF
+      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction
+      IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' )  THEN
+         srcv(jpr_sdrfty)%laction = .TRUE.
+         cpl_sdrfty = .TRUE.
+      ENDIF
+      srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period
+      IF( TRIM(sn_rcv_wper%cldes  ) == 'coupled' )  THEN
+         srcv(jpr_wper)%laction = .TRUE.
+         cpl_wper = .TRUE.
+      ENDIF
+      srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number
+      IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN
+         srcv(jpr_wnum)%laction = .TRUE.
+         cpl_wnum = .TRUE.
+      ENDIF
+      srcv(jpr_wstrf)%clname = 'O_WStrf'     ! stress fraction adsorbed by the wave
+      IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' )  THEN
+         srcv(jpr_wstrf)%laction = .TRUE.
+         cpl_wstrf = .TRUE.
+      ENDIF
+      srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient
+      IF( TRIM(sn_rcv_wdrag%cldes ) == 'coupled' )  THEN
+         srcv(jpr_wdrag)%laction = .TRUE.
+         cpl_wdrag = .TRUE.
+      ENDIF
+      ! 
       !                                                      ! ------------------------------- !
       !                                                      !   OPA-SAS coupling - rcv by opa !   
@@ -637,4 +727,5 @@
       !                                                      ! ------------------------- !
       ssnd(jps_fice)%clname = 'OIceFrc'
+      ssnd(jps_ficet)%clname = 'OIceFrcT' 
       ssnd(jps_hice)%clname = 'OIceTck'
       ssnd(jps_hsnw)%clname = 'OSnwTck'
@@ -645,4 +736,6 @@
       ENDIF
       
+      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 
+
       SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
       CASE( 'none'         )       ! nothing to do
@@ -665,4 +758,6 @@
       ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
       ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
+      ssnd(jps_ocxw)%clname = 'O_OCurxw' 
+      ssnd(jps_ocyw)%clname = 'O_OCuryw' 
       !
       ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
@@ -685,8 +780,29 @@
       END SELECT
 
+      ssnd(jps_ocxw:jps_ocyw)%nsgn = -1.   ! vectors: change of the sign at the north fold 
+        
+      IF( sn_snd_crtw%clvgrd == 'U,V' ) THEN 
+         ssnd(jps_ocxw)%clgrid = 'U' ; ssnd(jps_ocyw)%clgrid = 'V' 
+      ELSE IF( sn_snd_crtw%clvgrd /= 'T' ) THEN 
+         CALL ctl_stop( 'sn_snd_crtw%clvgrd must be equal to T' ) 
+      ENDIF 
+      IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) ssnd(jps_ocxw:jps_ocyw)%nsgn = 1. 
+      SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
+         CASE( 'none'                 )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .FALSE. 
+         CASE( 'oce only'             )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .TRUE. 
+         CASE( 'weighted oce and ice' )   !   nothing to do 
+         CASE( 'mixed oce-ice'        )   ; ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE. 
+         CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crtw%cldes' ) 
+      END SELECT 
+
       !                                                      ! ------------------------- !
       !                                                      !          CO2 flux         !
       !                                                      ! ------------------------- !
       ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
+
+      !                                                      ! ------------------------- ! 
+      !                                                      !     Sea surface height    ! 
+      !                                                      ! ------------------------- ! 
+      ssnd(jps_wlev)%clname = 'O_Wlevel' ;  IF( TRIM(sn_snd_wlev%cldes) == 'coupled' )   ssnd(jps_wlev)%laction = .TRUE. 
 
       !                                                      ! ------------------------------- !
@@ -837,4 +953,8 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
+      USE zdf_oce,  ONLY : ln_zdfqiao
+
+      IMPLICIT NONE
+
       INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
       INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
@@ -992,4 +1112,62 @@
       IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
 #endif
+      ! 
+      !                                                      ! ========================= ! 
+      !                                                      ! Mean Sea Level Pressure   !   (taum) 
+      !                                                      ! ========================= ! 
+      ! 
+      IF( srcv(jpr_mslp)%laction ) THEN                    ! UKMO SHELF effect of atmospheric pressure on SSH 
+          IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields 
+
+          r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization 
+          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer) 
+          apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1)                         !atmospheric pressure 
+    
+          IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible) 
+      END IF 
+      !
+      IF( ln_sdw ) THEN  ! Stokes Drift correction activated
+      !                                                      ! ========================= ! 
+      !                                                      !       Stokes drift u      !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_sdrftx)%laction ) zusd2dt(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)
+      !
+      !                                                      ! ========================= ! 
+      !                                                      !       Stokes drift v      !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_sdrfty)%laction ) zvsd2dt(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)
+      !
+      !                                                      ! ========================= ! 
+      !                                                      !      Wave mean period     !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)
+      !
+      !                                                      ! ========================= ! 
+      !                                                      !  Significant wave height  !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_hsig)%laction ) swh(:,:) = frcv(jpr_hsig)%z3(:,:,1)
+      !
+      !                                                      ! ========================= ! 
+      !                                                      !    Vertical mixing Qiao   !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_wnum)%laction .AND. ln_zdfqiao ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1)
+
+         ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode
+         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &
+                                                                    .OR. srcv(jpr_hsig)%laction ) THEN
+            CALL sbc_stokes()
+            IF( ln_zdfqiao .AND. .NOT. srcv(jpr_wnum)%laction ) CALL sbc_qiao()
+         ENDIF
+         IF( ln_zdfqiao .AND. srcv(jpr_wnum)%laction ) CALL sbc_qiao()
+      ENDIF
+      !                                                      ! ========================= ! 
+      !                                                      ! Stress adsorbed by waves  !
+      !                                                      ! ========================= ! 
+      IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1)
+
+      !                                                      ! ========================= ! 
+      !                                                      !   Wave drag coefficient   !
+      !                                                      ! ========================= ! 
+      IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw ) cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1)
 
       !  Fields received by SAS when OASIS coupling
@@ -1984,4 +2162,148 @@
       ENDIF
       !
+      !                                                      ! ------------------------- ! 
+      !                                                      !  Surface current to waves ! 
+      !                                                      ! ------------------------- ! 
+      IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN 
+          !     
+          !                                                  j+1  j     -----V---F 
+          ! surface velocity always sent from T point                    !       | 
+          !                                                       j      |   T   U 
+          !                                                              |       | 
+          !                                                   j   j-1   -I-------| 
+          !                                               (for I)        |       | 
+          !                                                             i-1  i   i 
+          !                                                              i      i+1 (for I) 
+          SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
+          CASE( 'oce only'             )      ! C-grid ==> T 
+             DO jj = 2, jpjm1 
+                DO ji = fs_2, fs_jpim1   ! vector opt. 
+                   zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) ) 
+                   zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) )  
+                END DO 
+             END DO 
+          CASE( 'weighted oce and ice' )    
+             SELECT CASE ( cp_ice_msh ) 
+             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = fs_2, fs_jpim1   ! vector opt. 
+                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   
+                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj) 
+                      zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
+                      zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = 2, jpim1   ! NO vector opt. 
+                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   
+                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   
+                      zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     & 
+                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                      zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     & 
+                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = 2, jpim1   ! NO vector opt. 
+                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   
+                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   
+                      zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     & 
+                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                      zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     & 
+                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             END SELECT 
+             CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. ) 
+          CASE( 'mixed oce-ice'        ) 
+             SELECT CASE ( cp_ice_msh ) 
+             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = fs_2, fs_jpim1   ! vector opt. 
+                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   & 
+                         &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
+                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
+                         &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = 2, jpim1   ! NO vector opt. 
+                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &    
+                         &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     & 
+                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &  
+                         &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     & 
+                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = 2, jpim1   ! NO vector opt. 
+                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &    
+                         &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     & 
+                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &  
+                         &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     & 
+                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             END SELECT 
+          END SELECT 
+         CALL lbc_lnk( zotx1, ssnd(jps_ocxw)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocyw)%clgrid, -1. ) 
+         ! 
+         ! 
+         IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components 
+         !                                                                        ! Ocean component 
+            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component  
+            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component  
+            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components  
+            zoty1(:,:) = ztmp2(:,:)  
+            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component 
+               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component  
+               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component  
+               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components  
+               zity1(:,:) = ztmp2(:,:) 
+            ENDIF 
+         ENDIF 
+         ! 
+!         ! spherical coordinates to cartesian -> 2 components to 3 components 
+!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN 
+!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents 
+!            ztmp2(:,:) = zoty1(:,:) 
+!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 ) 
+!            ! 
+!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities 
+!               ztmp1(:,:) = zitx1(:,:) 
+!               ztmp1(:,:) = zity1(:,:) 
+!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 ) 
+!            ENDIF 
+!         ENDIF 
+         ! 
+         IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid 
+         IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid 
+         !  
+      ENDIF 
+      ! 
+      IF( ssnd(jps_ficet)%laction ) THEN 
+         CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info ) 
+      END IF 
+      !                                                      ! ------------------------- ! 
+      !                                                      !   Water levels to waves   ! 
+      !                                                      ! ------------------------- ! 
+      IF( ssnd(jps_wlev)%laction ) THEN 
+         IF( ln_apr_dyn ) THEN  
+            IF( kt /= nit000 ) THEN  
+               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
+            ELSE  
+               ztmp1(:,:) = sshb(:,:)  
+            ENDIF  
+         ELSE  
+            ztmp1(:,:) = sshn(:,:)  
+         ENDIF  
+         CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
+      END IF 
       !
       !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
Index: /branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7167)
+++ /branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7168)
@@ -27,11 +27,22 @@
    PRIVATE
 
+   PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl
    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
+   ! Variables checking if the wave parameters are coupled (if not, they are read from file)
+   LOGICAL, PUBLIC     ::   cpl_hsig=.FALSE.
+   LOGICAL, PUBLIC     ::   cpl_phioc=.FALSE.
+   LOGICAL, PUBLIC     ::   cpl_sdrftx=.FALSE.
+   LOGICAL, PUBLIC     ::   cpl_sdrfty=.FALSE.
+   LOGICAL, PUBLIC     ::   cpl_wper=.FALSE.
+   LOGICAL, PUBLIC     ::   cpl_wnum=.FALSE.
+   LOGICAL, PUBLIC     ::   cpl_wstrf=.FALSE.
+   LOGICAL, PUBLIC     ::   cpl_wdrag=.FALSE.
+
+   INTEGER             ::   jpfld                ! number of files to read for stokes drift
+   INTEGER             ::   jp_usd               ! index of stokes drift  (i-component) (m/s)    at T-point
+   INTEGER             ::   jp_vsd               ! index of stokes drift  (j-component) (m/s)    at T-point
+   INTEGER             ::   jp_swh               ! index of significant wave hight      (m)      at T-point
+   INTEGER             ::   jp_wmp               ! 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
@@ -57,4 +68,113 @@
 CONTAINS
 
+   SUBROUTINE sbc_stokes( )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_stokes  ***
+      !!
+      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al.,
+      !!                2014 (DOI: 10.1175/JPO-D-14-0020.1)
+      !!
+      !! ** Method  : - Calculate Stokes transport speed 
+      !!              - Calculate horizontal divergence 
+      !!              - Integrate the horizontal divergenze from the bottom 
+      !! ** action  
+      !!---------------------------------------------------------------------
+      INTEGER                ::   jj,ji,jk 
+      REAL(wp)                       ::  ztransp, zsp0, zk, zus, zvs
+      REAL(wp), DIMENSION(jpi,jpj)   ::  zfac 
+      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace
+
+
+      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( 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, 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) &
+                  &                                + zfac(ji,jj+1) * tmask(ji,jj+1,1) )
+               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 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)
+            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
+      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)
+      END DO
+#if defined key_bdy
+      IF( lk_bdy ) THEN
+         DO jk = 1, jpkm1
+            wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
+         END DO
+      ENDIF
+#endif
+      CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv )
+      !
+   END SUBROUTINE sbc_stokes
+
+   SUBROUTINE sbc_qiao( )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_qiao  ***
+      !!
+      !! ** Purpose :   Qiao formulation for wave enhanced turbulence
+      !!                2010 (DOI: 10.1007/s10236-010-0326) 
+      !!
+      !! ** Method  : - 
+      !! ** action  
+      !!---------------------------------------------------------------------
+      INTEGER                ::   jj,ji
+
+      ! Calculate the module of the stokes drift on T grid
+      !-------------------------------------------------
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            tsd2d(ji,jj) = ((zusd2dt(ji,jj) * tmask(ji,jj,1))**2.0  + &
+            &               (zvsd2dt(ji,jj) * tmask(ji,jj,1))**2.0)**0.5
+         END DO
+      END DO
+      !
+   END SUBROUTINE sbc_qiao
+
    SUBROUTINE sbc_wave( kt )
       !!---------------------------------------------------------------------
@@ -71,17 +191,16 @@
       !! ** action  
       !!---------------------------------------------------------------------
-      USE zdf_oce,  ONLY : ln_zdfqiao
+      USE zdf_oce    ,  ONLY : ln_zdfqiao
+
+      IMPLICIT NONE
 
       INTEGER, INTENT( in  ) ::   kt       ! ocean time step
       !
       INTEGER                ::   ierror   ! return error code
-      INTEGER                ::   ifpr, jj,ji,jk 
+      INTEGER                ::   ifpr
       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), ALLOCATABLE, DIMENSION(:) ::   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
@@ -103,10 +222,12 @@
          !
          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' )
+            IF ( .NOT. cpl_wdrag ) 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' )
+            ENDIF
             ALLOCATE( cdn_wave(jpi,jpj) )
             cdn_wave(:,:) = 0.0
@@ -114,26 +235,54 @@
 
          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) )
-            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
+            IF ( .NOT. cpl_wstrf ) 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) )
+               CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
+            ENDIF
             ALLOCATE( tauoc_wave(jpi,jpj) )
             tauoc_wave(:,:) = 0.0
-        ENDIF
+         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
-            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' )
+            ! Find out how many fields have to be read from file if not coupled
+            jpfld=0
+            jp_usd=0; jp_vsd=0; jp_swh=0; jp_wmp=0
+            IF( .NOT. cpl_sdrftx ) THEN
+               jpfld=jpfld+1
+               jp_usd=jpfld
+            ENDIF
+            IF( .NOT. cpl_sdrfty ) THEN
+               jpfld=jpfld+1
+               jp_vsd=jpfld
+            ENDIF
+            IF( .NOT. cpl_hsig ) THEN
+               jpfld=jpfld+1
+               jp_swh=jpfld
+            ENDIF
+            IF( .NOT. cpl_wper ) THEN
+               jpfld=jpfld+1
+               jp_wmp=jpfld
+            ENDIF
+
+            ! Read from file only the non-coupled fields 
+            IF( jpfld > 0 ) THEN
+               ALLOCATE( slf_i(jpfld) )
+               IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd
+               IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd
+               IF( jp_swh > 0 ) slf_i(jp_swh) = sn_swh
+               IF( jp_wmp > 0 ) 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' )
+               !
+               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' )
+            ENDIF
             ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj) )
             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
@@ -145,9 +294,11 @@
             swh  (:,:)   = 0._wp   ;    wmp (:,:) = 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( 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' )
+               IF ( .NOT. cpl_wnum ) THEN
+                  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( 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' )
+               ENDIF
                ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) )
                wnum(:,:) = 0._wp ; tsd2d(:,:) = 0._wp
@@ -156,10 +307,10 @@
       ENDIF
       !
-      IF ( ln_cdgw ) THEN              !==  Neutral drag coefficient  ==!
+      IF ( ln_cdgw .AND. .NOT. cpl_wdrag ) 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  ==!
+      IF ( ln_tauoc .AND. .NOT. cpl_wstrf ) 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)
@@ -168,96 +319,29 @@
       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
-         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
+         ! Read from file only if the field is not coupled
+         IF( jpfld > 0 ) THEN
+            CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing
+            IF( jp_swh > 0 ) swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height
+            IF( jp_wmp > 0 ) wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
+            IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
+            IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
+         ENDIF
          !
+         ! Read also wave number if needed, so that it is available in coupling routines
+         IF ( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN
+            CALL fld_read( kt, nn_fsbc, sf_wn )      !* read wave parameters from external forcing
+            wnum(:,:) = sf_wn(1)%fnow(:,:,1)
+         ENDIF
+           
          !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014
          !(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( 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, 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) &
-                 &                                + zfac(ji,jj+1) * tmask(ji,jj+1,1) )
-
-                  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 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)
-               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
-         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)
-         END DO
-#if defined key_bdy
-         IF( lk_bdy ) THEN
-            DO jk = 1, jpkm1
-               wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
-            END DO
-         ENDIF
-#endif
-         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  +     &
-                &               (sf_sd(jp_vsd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0)**0.5
-               END DO
-            END DO
-         ENDIF
-         !
+         ! Calculate only if no necessary fields are coupled, if not calculate later after coupling
+         IF( jpfld == 4 ) THEN
+            CALL sbc_stokes()
+            IF ( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN
+               CALL sbc_qiao()
+            ENDIF
+         ENDIF
       ENDIF
       !
Index: /branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7167)
+++ /branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7168)
@@ -231,6 +231,6 @@
                          CALL dyn_adv       ( kstp )  ! advection (vector or flux form)
                          CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis
+      IF( ln_stcor    )  CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing
                          CALL dyn_ldf       ( kstp )  ! lateral mixing
-      IF( ln_stcor    )  CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing
                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
                          CALL dyn_spg       ( kstp )  ! surface pressure gradient
