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 7166)
+++ 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
