Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 12448)
@@ -374,4 +374,6 @@
          <field id="uocet"        long_name="ocean transport along i-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_U_3D" />
          <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="1e-3*m/s"   grid_ref="grid_U_3D" />
+	 <field id="ustokes"      long_name="Stokes Drift Velocity i-axis"                           standard_name="StokesDrift_x_velocity"      unit="m/s"        grid_ref="grid_U_3D" />  
+	 <field id="ustokes_e3u"  long_name="Stokes Drift Velocity i-axis  (thickness weighted)"                                                 unit="m/s"        grid_ref="grid_U_3D"  > ustokes * e3u </field>
 
          <!-- variables available with MLE -->
@@ -409,4 +411,6 @@
          <field id="vocet"        long_name="ocean transport along j-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_V_3D" />
          <field id="voces"        long_name="ocean transport along j-axis times salinity (CRS)"                                                  unit="1e-3*m/s"   grid_ref="grid_V_3D" />
+	 <field id="vstokes"      long_name="Stokes Drift Velocity j-axis"                           standard_name="StokesDrift_y_velocity"      unit="m/s"        grid_ref="grid_V_3D" />  
+	 <field id="vstokes_e3v"  long_name="Stokes Drift Velocity j-axis  (thickness weighted)"                                                 unit="m/s"        grid_ref="grid_V_3D"  > vstokes * e3v </field>
 
          <!-- variables available with MLE -->
@@ -437,4 +441,5 @@
         <field id="woce"         long_name="ocean vertical velocity"              standard_name="upward_sea_water_velocity"   unit="m/s"  />
         <field id="wocetr_eff"   long_name="effective ocean vertical transport"                                               unit="m3/s" />
+        <field id="wstokes"      long_name="Stokes Drift vertical velocity"       standard_name="upward_StokesDrift_velocity" unit="m/s" />
 
         <!-- woce_eiv: available with key_traldf_eiv and key_diaeiv -->
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 12448)
@@ -240,6 +240,7 @@
    ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core)
    ln_blk_mfs  = .false.   !  MFS bulk formulation                      (T => fill namsbc_mfs )
-   ln_cpl      = .false.   !  atmosphere coupled   formulation          ( requires key_oasis3 )
-   ln_mixcpl   = .false.   !  forced-coupled mixed formulation          ( requires key_oasis3 )
+   ln_cpl      = .false.   !  coupled   formulation                       ( requires key_oasis3 ) 
+   ln_mixcpl   = .false.   !  forced-coupled mixed atmosphere formulation ( requires key_oasis3 ) 
+   ln_wavcpl   = .false.   !  forced-coupled mixed wave       formulation ( requires key_oasis3 )
    nn_components = 0       !  configuration of the opa-sas OASIS coupling
                            !  =0 no opa-sas OASIS coupling: default single executable configuration
@@ -264,7 +265,5 @@
                            !     =1 global mean of e-p-r set to zero at each time step
                            !     =2 annual global mean of e-p-r set to zero
-   ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave)
-   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
-   ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
+   ln_wave = .false.       !  Activate coupling with wave (T => fill namsbc_wave)   
    nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) ,
                            !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field)
@@ -274,4 +273,9 @@
                            !  = 1  Average and redistribute per-category fluxes, forced mode only, not yet implemented coupled
                            !  = 2  Redistribute a single flux over categories (coupled mode only)
+   nn_drag   = 0      !  formula to calculate momentum from the wind components 
+                           ! = 0 UKMO SHELF formulation 
+                           ! = 1 standard formulation with forced of coupled drag coefficient 
+                           ! = 2 standard formulation with constant drag coefficient 
+                           ! = 3 momentum calculated from core forcing fields
 /
 !-----------------------------------------------------------------------
@@ -296,5 +300,8 @@
    sn_emp      = 'emp'       ,        24         , 'emp'     , .false.      , .false., 'yearly'  , ''       , ''       , ''
 
-   cn_dir      = './'      !  root directory for the location of the flux files
+   cn_dir       = './'      !  root directory for the location of the flux files 
+   ln_shelf_flx = .false.   !  UKMO SHELF specific flux flag - read from file wind components instead of momentum  
+   ln_rel_wind  = .false.   !  UKMO SHELF - calculate momentum from wind speed relative to currents  
+   rn_wfac      = 1.0       !  UKMO SHELF - multiplicative factor for ocean/wind velocity
 /
 !-----------------------------------------------------------------------
@@ -363,4 +370,7 @@
    sn_snd_crt    =       'none'                 ,    'no'    , 'spherical' , 'eastward-northward' ,  'T'
    sn_snd_co2    =       'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_crtw   =       'none'                 ,    'no'    ,     ''      ,         ''           , 'U,V'  
+   sn_snd_ifrac  =       'none'                 ,    'no'    ,     ''      ,         ''           ,   ''  
+   sn_snd_wlev   =       'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
 ! receive
    sn_rcv_w10m   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
@@ -374,8 +384,21 @@
    sn_rcv_cal    =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
    sn_rcv_co2    =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_hsig   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_iceflx =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_mslp   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_phioc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_sdrft  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_wfreq  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_wnum   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_tauoc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_tauw   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_wdrag  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
 !
    nn_cplmodel   =     1     !  Maximum number of models to/from which NEMO is potentialy sending/receiving data
    ln_usecplmask = .false.   !  use a coupling mask file to merge data received from several models
                              !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   ln_usernfmask = .false.   !  use a runoff mask file to merge data received from several models 
+                             !   -> file rnfmask.nc with the float variable called rnfmask (jpi,jpj,nn_cplmodel)
 /
 !-----------------------------------------------------------------------
@@ -984,9 +1007,10 @@
    rn_clim_galp  = 0.267   !  galperin limit
    ln_sigpsi     = .true.  !  Activate or not Burchard 2001 mods on psi schmidt number in the wb case
-   rn_crban      = 100.    !  Craig and Banner 1994 constant for wb tke flux
+   rn_crban_default = 100.    !  Craig and Banner 1994 constant for wb tke flux
    rn_charn      = 70000.  !  Charnock constant for wb induced roughness length
    rn_hsro       =  0.02   !  Minimum surface roughness
    rn_frac_hs    =   1.3   !  Fraction of wave height as roughness (if nn_z0_met=2)
-   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2)
+   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3) 
+      !                             ! =3 requires ln_wave=T 
    nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum)
    nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum)
@@ -1300,14 +1324,35 @@
 /
 !-----------------------------------------------------------------------
-&namsbc_wave   ! External fields from wave model
+&namsbc_wave   ! External fields from wave model                        (ln_wave=T)
 !-----------------------------------------------------------------------
 !              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
 !              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model 
+   ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift 
+   ln_stcor    = .false.   !  Activate Stokes Coriolis term 
+   ln_tauoc    = .false.   !  Activate ocean stress modified by external wave induced stress 
+   ln_tauw     = .false.   !  Activate ocean stress components from wave model 
+   ln_phioc    = .false.   !  Activate wave to ocean energy 
+   ln_rough    = .false.   !  Wave roughness equals the significant wave height 
+   ln_zdfqiao  = .false.   !  Enhanced wave vertical mixing Qiao (2010) 
+   nn_sdrift   =  0   !  Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift 
+                           ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 
+                           ! = 1 Phillips:                      v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))] 
+   nn_wmix     =  0   !  type of wave breaking mixing 
+                           ! = 0 Craig and Banner formulation (original NEMO formulation)  
+                           ! = 1 Janssen formulation (no assumption of direct energy conversion)
+   sn_cdg      =  'cdw_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
-   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
-!
-   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+   sn_hsw      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_wfr      =  'sdw_wave' ,        1          , 'wave_freq'  ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_tauwx    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_tauwy    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_phioc    =  'sdw_wave' ,        1          , 'wave_energy',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   cn_dir      = './'  !  root directory for the location of wave forcing files
+!  
 /
 !-----------------------------------------------------------------------
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/CONFIG/cfg.txt	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/CONFIG/cfg.txt	(revision 12448)
@@ -11,2 +11,3 @@
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
+GYRE_LONG OPA_SRC
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90	(revision 12448)
@@ -124,7 +124,6 @@
 
     CASE DEFAULT
-       IF(lwp) WRITE(numout,*) ' E R R O R : bad cd_type in dia_wri_dimg '
-       STOP 'dia_wri_dimg'
-
+       WRITE(numout,*) ' E R R O R : bad cd_type in dia_wri_dimg' 
+       CALL ctl_stop( 'STOP', 'dia_wri_dimg :bad cd_type in dia_wri_dimg ' ) 
     END SELECT
 
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 12448)
@@ -41,4 +41,5 @@
    USE zdfddm          ! vertical  physics: double diffusion
    USE diahth          ! thermocline diagnostics
+   USE sbcwave         ! wave parameters
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
@@ -709,4 +710,8 @@
          CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
             &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
+         IF( ln_wave .AND. ln_sdw) THEN  
+            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd  
+               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )  
+         ENDIF 
          IF( ln_traldf_gdia ) THEN
             CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
@@ -727,4 +732,8 @@
          CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
             &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
+         IF( ln_wave .AND. ln_sdw) THEN  
+            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd  
+               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )  
+         ENDIF  
          IF( ln_traldf_gdia ) THEN
             CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
@@ -763,4 +772,9 @@
                &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
          ENDIF
+           
+         IF( ln_wave .AND. ln_sdw) THEN  
+            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd  
+               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )  
+         ENDIF  
          !                                                                                      !!! nid_W : 2D
 #if defined key_traldf_c2d
@@ -943,4 +957,10 @@
 #endif
 
+      IF( ln_wave .AND. ln_sdw ) THEN  
+         CALL histwrite( nid_U, "sdzocrtx", it, usd           , ndim_U , ndex_U )    ! i-StokesDrift-current  
+         CALL histwrite( nid_V, "sdmecrty", it, vsd           , ndim_V , ndex_V )    ! j-StokesDrift-current  
+         CALL histwrite( nid_W, "sdvecrtz", it, wsd           , ndim_T , ndex_T )    ! StokesDrift vert. current  
+      ENDIF  
+       
       ! 3. Close all files
       ! ---------------------------------------
@@ -1048,4 +1068,13 @@
             &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
       END IF
+      !  
+      IF( ln_wave .AND. ln_sdw ) THEN  
+         CALL histdef( id_i, "sdzocrtx", "Stokes Drift Zonal", "m/s"    , &   ! StokesDrift zonal current  
+            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )  
+         CALL histdef( id_i, "sdmecrty", "Stokes Drift Merid", "m/s"    , &   ! StokesDrift meridonal current  
+            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )  
+         CALL histdef( id_i, "sdvecrtz", "Stokes Drift Vert", "m/s"    , &   ! StokesDrift vertical current  
+            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )  
+      ENDIF  
 
 #if defined key_lim2
@@ -1065,4 +1094,9 @@
 
       ! Write all fields on T grid
+      IF( ln_wave .AND. ln_sdw ) THEN  
+         CALL histwrite( id_i, "sdzocrtx", kt, usd           , jpi*jpj*jpk, idex)     ! now StokesDrift i-velocity  
+         CALL histwrite( id_i, "sdmecrty", kt, vsd           , jpi*jpj*jpk, idex)     ! now StokesDrift j-velocity  
+         CALL histwrite( id_i, "sdvecrtz", kt, wsd           , jpi*jpj*jpk, idex)     ! now StokesDrift k-velocity  
+      ENDIF
       CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
       CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 12448)
@@ -11,4 +11,5 @@
    !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
    !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
+   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence
    !!---------------------------------------------------------------------
 #if defined key_dynspg_ts   ||   defined key_esopa
@@ -31,4 +32,6 @@
    USE sbctide         ! tides
    USE updtide         ! tide potential
+   USE sbcwave         ! surface wave 
+   ! 
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -45,7 +48,4 @@
 #if defined key_agrif
    USE agrif_opa_interp ! agrif
-#endif
-#if defined key_asminc   
-   USE asminc          ! Assimilation increment
 #endif
 
@@ -333,4 +333,9 @@
          END DO
       END DO
+              
+      !!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....  
+      !!gm  Is it correct to do so ?   I think so...  
+              
+              
       !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
       !                                   ! --------------------------------------------------------
@@ -462,10 +467,4 @@
                 &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
       ENDIF
-#if defined key_asminc
-      !                                         ! Include the IAU weighted SSH increment
-      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
-         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
-      ENDIF
-#endif
       !                                   !* Fill boundary data arrays with AGRIF
       !                                   ! -------------------------------------
@@ -475,10 +474,10 @@
       !
       ! -----------------------------------------------------------------------
-      !  Phase 2 : Integration of the barotropic equations 
+      !  Phase 2 : Integration of the barotropic equations
       ! -----------------------------------------------------------------------
       !
       !                                             ! ==================== !
       !                                             !    Initialisations   !
-      !                                             ! ==================== !  
+      !                                             ! ==================== !
       ! Initialize barotropic variables:      
       IF( ll_init )THEN
@@ -489,4 +488,9 @@
          ub_e   (:,:) = 0._wp
          vb_e   (:,:) = 0._wp
+      ENDIF
+      ! 
+      ! Moved here to allow merging with other branch  
+      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary  
+         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)  
       ENDIF
       !
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 12448)
@@ -16,4 +16,5 @@
    !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
    !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 
+   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
    !!----------------------------------------------------------------------
 
@@ -32,4 +33,7 @@
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
+   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force) 
+   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force 
+   ! 
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE prtctl         ! Print control
@@ -91,8 +95,8 @@
       !
       CASE ( -1 )                                      ! esopa: test all possibility with control print
-         CALL vor_ene( kt, ntot, ua, va )
+         CALL vor_ene( kt, ntot, un, vn, ua, va )
          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL vor_ens( kt, ntot, ua, va )
+         CALL vor_ens( kt, ntot, un, vn, ua, va )
          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
@@ -100,5 +104,5 @@
          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL vor_een( kt, ntot, ua, va )
+         CALL vor_een( kt, ntot, un, vn, ua, va )
          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
@@ -108,5 +112,5 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend
+            CALL vor_ene( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -114,10 +118,11 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend
+            CALL vor_ene( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity
+                             CALL vor_ene( kt, ntot, un, vn, ua, va )    ! total vorticity trend 
+            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
          ENDIF
          !
@@ -126,5 +131,5 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend
+            CALL vor_ens( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -132,10 +137,11 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend
+            CALL vor_ens( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity
+                             CALL vor_ens( kt, ntot, un, vn, ua, va )    ! total vorticity 
+            IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
          ENDIF
          !
@@ -144,5 +150,5 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
+            CALL vor_ens( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend (ens)
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -150,10 +156,12 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
+            CALL vor_ene( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend (ene)
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
+                             CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens) 
+                             CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene) 
+            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
          ENDIF
          !
@@ -162,5 +170,5 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend
+            CALL vor_een( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -168,10 +176,11 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend
+            CALL vor_een( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_een( kt, ntot, ua, va )                ! total vorticity
+                             CALL vor_een( kt, ntot, un, vn, ua, va )    ! total vorticity 
+            IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
          ENDIF
          !
@@ -189,5 +198,5 @@
 
 
-   SUBROUTINE vor_ene( kt, kvor, pua, pva )
+   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE vor_ene  ***
@@ -214,9 +223,10 @@
       !!----------------------------------------------------------------------
       !
-      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-      !                                                           ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
+      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index 
+      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ; 
+      !                                                                ! =nrvm (relative vorticity or metric) 
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend 
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend 
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -250,6 +260,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
+                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
                END DO
@@ -260,6 +270,6 @@
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = ( ff (ji,jj)                                                                       &
-                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                       &       + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
+                       &           - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
                        &       )
@@ -270,9 +280,9 @@
          IF( ln_sco ) THEN
             zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
-            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * pun(:,:,jk) 
+            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * pvn(:,:,jk)
          ELSE
-            zwx(:,:) = e2u(:,:) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
          ENDIF
 
@@ -418,5 +428,5 @@
 
 
-   SUBROUTINE vor_ens( kt, kvor, pua, pva )
+   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_ens  ***
@@ -443,9 +453,10 @@
       !!----------------------------------------------------------------------
       !
-      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-         !                                                        ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
+      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index 
+      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ; 
+         !                                                             ! =nrvm (relative vorticity or metric) 
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend 
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend 
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
       !
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
@@ -479,6 +490,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
+                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
                END DO
@@ -489,6 +500,6 @@
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = ( ff (ji,jj)                                                                       &
-                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                       &       + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
+                       &           - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
                        &       )
@@ -501,6 +512,6 @@
                DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
                   zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
-                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
-                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
+                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
+                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
                END DO
             END DO
@@ -508,6 +519,6 @@
             DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
                DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
-                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
-                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
+                  zwx(ji,jj) = e2u(ji,jj) * pun(ji,jj,jk)
+                  zwy(ji,jj) = e1v(ji,jj) * pvn(ji,jj,jk)
                END DO
             END DO
@@ -536,5 +547,5 @@
 
 
-   SUBROUTINE vor_een( kt, kvor, pua, pva )
+   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_een  ***
@@ -554,9 +565,10 @@
       !!----------------------------------------------------------------------
       !
-      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-      !                                                           ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
+      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index 
+      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ; 
+      !                                                                ! =nrvm (relative vorticity or metric) 
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend 
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend 
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
       !!
       INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
@@ -604,5 +616,7 @@
                      ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
                         &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
-                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3
+                     IF( ze3 /= 0._wp ) THEN  ;   ze3f(ji,jj,jk) = 4.0_wp / ze3 
+                     ELSE                     ;   ze3f(ji,jj,jk) = 0._wp 
+                     ENDIF
                   END DO
                END DO
@@ -616,5 +630,7 @@
                      zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
                         &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
-                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3
+                     IF( ze3 /= 0._wp ) THEN  ;  ze3f(ji,jj,jk) = zmsk / ze3 
+                     ELSE                     ;  ze3f(ji,jj,jk) = 0._wp 
+                     ENDIF
                   END DO
                END DO
@@ -643,6 +659,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
+                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
                END DO
@@ -655,6 +671,6 @@
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = ( ff (ji,jj)                                                                       &
-                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                       &       + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
+                       &           - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
                        &       ) * ze3f(ji,jj,jk)
@@ -664,6 +680,6 @@
          END SELECT
 
-         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
-         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
+         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * pun(:,:,jk)
+         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * pvn(:,:,jk)
 
          ! Compute and add the vorticity term trend
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 12448)
@@ -62,4 +62,5 @@
    USE lbcnfd         ! north fold treatment
    USE in_out_manager ! I/O manager
+   USE mod_oasis      ! coupling routines
 
    IMPLICIT NONE
@@ -2023,6 +2024,13 @@
       !!----------------------------------------------------------------------
       !
+#if defined key_oasis3 || defined key_oasis3mct   
+      ! If we're trying to shut down cleanly then we need to consider the fact   
+      ! that this could be part of an MPMD configuration - we don't want to   
+      ! leave other components deadlocked.   
+      CALL oasis_abort(nproc,"mppstop","NEMO initiated abort")   
+#else   
       CALL mppsync
       CALL mpi_finalize( info )
+#endif 
       !
    END SUBROUTINE mppstop
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 12448)
@@ -24,5 +24,5 @@
    !!   cpl_finalize : finalize the coupled mode communication
    !!----------------------------------------------------------------------
-#if defined key_oasis3
+#if defined key_oasis3 || defined key_oasis3mct
    USE mod_oasis                    ! OASIS3-MCT module
 #endif
@@ -32,4 +32,9 @@
    USE lbclnk                       ! ocean lateral boundary conditions (or mpp link)
 
+#if defined key_cpl_rootexchg   
+   USE lib_mpp, only : mppsync   
+   USE lib_mpp, only : mppscatter,mppgather   
+#endif    
+
    IMPLICIT NONE
    PRIVATE
@@ -41,4 +46,15 @@
    PUBLIC   cpl_freq
    PUBLIC   cpl_finalize
+#if defined key_mpp_mpi   
+   INCLUDE 'mpif.h'   
+#endif   
+      
+   INTEGER, PARAMETER         :: localRoot  = 0   
+   LOGICAL                    :: commRank            ! true for ranks doing OASIS communication   
+#if defined key_cpl_rootexchg   
+   LOGICAL                    :: rootexchg =.true.   ! logical switch    
+#else   
+   LOGICAL                    :: rootexchg =.false.  ! logical switch    
+#endif    
 
    INTEGER, PUBLIC            ::   OASIS_Rcv  = 1    !: return code if received field
@@ -46,5 +62,5 @@
    INTEGER                    ::   ncomp_id          ! id returned by oasis_init_comp
    INTEGER                    ::   nerror            ! return error code
-#if ! defined key_oasis3
+#if ! defined key_oasis3 && ! defined key_oasis3mct
    ! OASIS Variables not used. defined only for compilation purpose
    INTEGER                    ::   OASIS_Out         = -1
@@ -65,5 +81,5 @@
    INTEGER                    ::   nsnd         ! total number of fields sent 
    INTEGER                    ::   ncplmodel    ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
-   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=50   ! Maximum number of coupling fields
+   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=55   ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields
@@ -82,4 +98,6 @@
 
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   exfld   ! Temporary buffer for receiving
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tbuf  ! Temporary buffer for sending / receiving    
+   INTEGER, PUBLIC :: localComm    
 
    !!----------------------------------------------------------------------
@@ -120,4 +138,5 @@
       IF ( nerror /= OASIS_Ok ) &
          CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' )
+      localComm = kl_comm
       !
    END SUBROUTINE cpl_init
@@ -149,4 +168,11 @@
       IF(lwp) WRITE(numout,*)
 
+      commRank = .false.   
+      IF ( rootexchg ) THEN   
+         IF ( nproc == localRoot ) commRank = .true.   
+      ELSE   
+         commRank = .true.   
+      ENDIF   
+
       ncplmodel = kcplmodel
       IF( kcplmodel > nmaxcpl ) THEN
@@ -172,10 +198,4 @@
       ishape(:,2) = (/ 1, nlej-nldj+1 /)
       !
-      ! ... Allocate memory for data exchange
-      !
-      ALLOCATE(exfld(nlei-nldi+1, nlej-nldj+1), stat = nerror)
-      IF( nerror > 0 ) THEN
-         CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld')   ;   RETURN
-      ENDIF
       !
       ! -----------------------------------------------------------------
@@ -183,4 +203,13 @@
       ! -----------------------------------------------------------------
       
+      IF ( rootexchg ) THEN   
+       
+         paral(1) = 2              ! box partitioning   
+         paral(2) = 0              ! NEMO lower left corner global offset        
+         paral(3) = jpiglo         ! local extent in i    
+         paral(4) = jpjglo         ! local extent in j   
+         paral(5) = jpiglo         ! global extent in x   
+       
+      ELSE    
       paral(1) = 2                                              ! box partitioning
       paral(2) = jpiglo * (nldj-1+njmpp-1) + (nldi-1+nimpp-1)   ! NEMO lower left corner global offset    
@@ -196,5 +225,22 @@
       ENDIF
       
-      CALL oasis_def_partition ( id_part, paral, nerror )
+      ENDIF 
+      IF ( commRank )  CALL oasis_def_partition ( id_part, paral, nerror )   
+       
+      ! ... Allocate memory for data exchange   
+      !   
+      ALLOCATE(exfld(paral(3), paral(4)), stat = nerror)   
+      IF( nerror > 0 ) THEN   
+         CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld')   ;   RETURN   
+      ENDIF   
+      IF ( rootexchg ) THEN   
+         ! Should possibly use one of the work arrays for tbuf really   
+         ALLOCATE(tbuf(jpi, jpj, jpnij), stat = nerror)   
+         IF( nerror > 0 ) THEN   
+            CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating tbuf') ; RETURN   
+         ENDIF   
+      ENDIF                
+      !   
+      IF (commRank ) THEN 
       !
       ! ... Announce send variables. 
@@ -288,4 +334,6 @@
          ENDIF
       END DO
+      !  
+      ENDIF  ! commRank=true  
       
       !------------------------------------------------------------------
@@ -293,6 +341,8 @@
       !------------------------------------------------------------------
       
-      CALL oasis_enddef(nerror)
-      IF( nerror /= OASIS_Ok )   CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef')
+      IF ( commRank ) THEN        
+         CALL oasis_enddef(nerror)   
+         IF( nerror /= OASIS_Ok )   CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef')   
+      ENDIF   
       !
    END SUBROUTINE cpl_define
@@ -311,5 +361,5 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdata
       !!
-      INTEGER                                   ::   jc,jm     ! local loop index
+      INTEGER                                   ::   jn,jc,jm  ! local loop index
       !!--------------------------------------------------------------------
       !
@@ -320,5 +370,23 @@
         
             IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN
-               CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(nldi:nlei, nldj:nlej,jc), kinfo )
+               IF ( rootexchg ) THEN   
+                  !   
+                  ! collect data on the local root process   
+                  !   
+                  CALL mppgather (pdata(:,:,jc),localRoot,tbuf)    
+                  CALL mppsync    
+                         
+                  IF ( nproc == localRoot ) THEN   
+                     DO jn = 1, jpnij   
+                        exfld(nimppt(jn)-1+nldit(jn):nimppt(jn)+nleit(jn)-1,njmppt(jn)-1+nldjt(jn):njmppt(jn)+nlejt(jn)-1)= &   
+                          tbuf(nldit(jn):nleit(jn),nldjt(jn):nlejt(jn),jn)   
+                     ENDDO   
+                     ! snd data to OASIS3   
+                     CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, exfld, kinfo )   
+                  ENDIF   
+               ELSE   
+                  ! snd data to OASIS3   
+                  CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(nldi:nlei, nldj:nlej,jc), kinfo )   
+               ENDIF   
                
                IF ( ln_ctl ) THEN        
@@ -358,5 +426,5 @@
       INTEGER                   , INTENT(  out) ::   kinfo     ! OASIS3 info argument
       !!
-      INTEGER                                   ::   jc,jm     ! local loop index
+      INTEGER                                   ::   jn,jc,jm  ! local loop index
       LOGICAL                                   ::   llaction, llfisrt
       !!--------------------------------------------------------------------
@@ -372,6 +440,9 @@
 
             IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN
-
-               CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld, kinfo )         
+               !   
+               ! receive data from OASIS3   
+               !   
+               IF ( commRank )   CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld, kinfo )   
+               IF ( rootexchg )  CALL MPI_BCAST ( kinfo, 1, MPI_INTEGER, localRoot, localComm, nerror ) 
                
                llaction =  kinfo == OASIS_Recvd   .OR. kinfo == OASIS_FromRest .OR.   &
@@ -384,5 +455,19 @@
                   kinfo = OASIS_Rcv
                   IF( llfisrt ) THEN 
-                     pdata(nldi:nlei,nldj:nlej,jc) =                                 exfld(:,:) * pmask(nldi:nlei,nldj:nlej,jm)
+                     IF ( rootexchg ) THEN   
+                        ! distribute data to processes   
+                        !   
+                        IF ( nproc == localRoot ) THEN   
+                           DO jn = 1, jpnij   
+                              tbuf(nldit(jn):nleit(jn),nldjt(jn):nlejt(jn),jn)=          &   
+                              exfld(nimppt(jn)-1+nldit(jn):nimppt(jn)+nleit(jn)-1,njmppt(jn)-1+nldjt(jn):njmppt(jn)+nlejt(jn)-1)   
+                              ! NOTE: we are missing combining this with pmask (see else below)   
+                           ENDDO   
+                        ENDIF   
+                        CALL mppscatter(tbuf,localRoot,pdata(:,:,jc))    
+                        CALL mppsync   
+                     ELSE   
+                        pdata(nldi:nlei, nldj:nlej, jc) = exfld(:,:) * pmask(nldi:nlei,nldj:nlej,jm)   
+                     ENDIF   
                      llfisrt = .FALSE.
                   ELSE
@@ -462,8 +547,12 @@
 #if defined key_oa3mct_v3
          CALL oasis_get_freqs(id, mop, 1, itmp, info)
-#else
+#endif
+#if defined key_oasis3
          CALL oasis_get_freqs(id,      1, itmp, info)
 #endif
          cpl_freq = itmp(1)
+#if defined key_oasis3mct  
+         cpl_freq = namflddti( id )  
+#endif 
       ENDIF
       !
@@ -481,4 +570,5 @@
       !
       DEALLOCATE( exfld )
+      IF ( rootexchg ) DEALLOCATE ( tbuf )
       IF (nstop == 0) THEN
          CALL oasis_terminate( nerror )         
@@ -489,5 +579,5 @@
    END SUBROUTINE cpl_finalize
 
-#if ! defined key_oasis3
+#if ! defined key_oasis3 && ! defined key_oasis3mct
 
    !!----------------------------------------------------------------------
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 12448)
@@ -51,5 +51,5 @@
 
    SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   &
-                       px2 , py2 )
+                       px2 , py2, kchoix )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE repcmo  ***
@@ -67,10 +67,23 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2          ! i-componante (defined at u-point)
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   py2          ! j-componante (defined at v-point)
-      !!----------------------------------------------------------------------
-      
-      ! Change from geographic to stretched coordinate
-      ! ----------------------------------------------
-      CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )
-      CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )
+      INTEGER, INTENT( IN )                       ::   kchoix       ! type of transformation   
+                                                                    ! = 1 change from geographic to model grid.   
+                                                                    ! =-1 change from model to geographic grid   
+      !!----------------------------------------------------------------------  
+        
+      SELECT CASE (kchoix)   
+      CASE ( 1)   
+        ! Change from geographic to stretched coordinate   
+        ! ----------------------------------------------   
+        
+        CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )   
+        CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )   
+      CASE (-1)   
+        ! Change from stretched to geographic coordinate   
+        ! ----------------------------------------------   
+        
+        CALL rot_rep( pxu1, pyu1, 'U', 'ij->e',px2 )   
+        CALL rot_rep( pxv1, pyv1, 'V', 'ij->n',py2 )   
+      END SELECT   
       
    END SUBROUTINE repcmo
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 12448)
@@ -35,5 +35,5 @@
    LOGICAL , PUBLIC ::   ln_blk_core    !: CORE bulk formulation
    LOGICAL , PUBLIC ::   ln_blk_mfs     !: MFS  bulk formulation
-#if defined key_oasis3
+#if defined key_oasis3 || defined key_oasis3mct
    LOGICAL , PUBLIC ::   lk_oasis = .TRUE.  !: OASIS used
 #else
@@ -42,8 +42,12 @@
    LOGICAL , PUBLIC ::   ln_cpl         !: ocean-atmosphere coupled formulation
    LOGICAL , PUBLIC ::   ln_mixcpl      !: ocean-atmosphere forced-coupled mixed formulation
+   LOGICAL , PUBLIC ::   ln_wavcpl      !: ocean-wave forced-coupled mixed formulation 
+   LOGICAL , PUBLIC ::   ll_purecpl     !: ocean-atmosphere or ocean-wave pure coupled formulation
    LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr)
    LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths
    LOGICAL , PUBLIC ::   ln_ssr         !: Sea Surface restoring on SST and/or SSS      
    LOGICAL , PUBLIC ::   ln_apr_dyn     !: Atmospheric pressure forcing used on dynamics (ocean & ice)
+   LOGICAL, PUBLIC  ::   ln_usernfmask = .false.   !  use a runoff mask file to merge data received from several models 
+                                                   !   -> file rnfmask.nc with the float variable called rnfmask (jpi,jpj,nn_cplmodel)
    INTEGER , PUBLIC ::   nn_ice         !: flag for ice in the surface boundary condition (=0/1/2/3)
    INTEGER , PUBLIC ::   nn_isf         !: flag for isf in the surface boundary condition (=0/1/2/3/4) 
@@ -64,5 +68,13 @@
    LOGICAL , PUBLIC ::   ln_wave        !: true if some coupling with wave model
    LOGICAL , PUBLIC ::   ln_cdgw        !: true if neutral drag coefficient from wave model
-   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model
+   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d Stokes drift from wave model 
+   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used
+   LOGICAL , PUBLIC ::   ln_tauoc       !: true if normalized stress from wave is used  
+   LOGICAL , PUBLIC ::   ln_tauw        !: true if ocean stress components from wave is used  
+   LOGICAL , PUBLIC ::   ln_phioc       !: true if wave energy to ocean is used  
+   LOGICAL , PUBLIC ::   ln_rough       !: true if wave roughness equals significant wave height 
+   LOGICAL , PUBLIC ::   ln_zdfqiao     !: Enhanced wave vertical mixing Qiao(2010) formulation flag 
+   INTEGER , PUBLIC ::   nn_drag        ! type of formula to calculate wind stress from wind components 
+   INTEGER , PUBLIC ::   nn_wmix        ! type of wave breaking mixing
    !
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
@@ -91,4 +103,18 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_iam_sas  = 2      !: Multi executable configuration - SAS component
                                                          !  (internal OASIS coupling)
+   !!---------------------------------------------------------------------- 
+   !!           wind stress definition 
+   !!---------------------------------------------------------------------- 
+   INTEGER, PUBLIC, PARAMETER ::   jp_ukmo  = 0        ! UKMO SHELF formulation 
+   INTEGER, PUBLIC, PARAMETER ::   jp_std   = 1        ! standard formulation with forced or coupled drag coefficient  
+   INTEGER, PUBLIC, PARAMETER ::   jp_const = 2        ! standard formulation with constant drag coefficient  
+   INTEGER, PUBLIC, PARAMETER ::   jp_mcore = 3        ! momentum calculated from core forcing fields  
+
+   !!---------------------------------------------------------------------- 
+   !!           Wave mixing vertical parameterization 
+   !!---------------------------------------------------------------------- 
+   INTEGER, PUBLIC, PARAMETER ::   jp_craigbanner = 0  ! Craig and Banner formulation (original NEMO formulation - 
+                                                       ! direct conversion of mechanical to turbulent energy) 
+   INTEGER, PUBLIC, PARAMETER ::   jp_janssen     = 1  ! Janssen formulation - no assumption on direct energy conversion
    !!----------------------------------------------------------------------
    !!              Ocean Surface Boundary Condition fields
@@ -128,4 +154,5 @@
 #endif
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xrnfmask          !: coupling mask for ln_usernfmask (warning: allocated in sbcrnf)
 
    !!----------------------------------------------------------------------
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 12448)
@@ -26,7 +26,8 @@
    
    !                                !!* namsbc_apr namelist (Atmospheric PRessure) *
+   LOGICAL, PUBLIC ::   cpl_mslp = .FALSE. ! Is the pressure read from coupling?
    LOGICAL, PUBLIC ::   ln_apr_obc   !: inverse barometer added to OBC ssh data 
    LOGICAL, PUBLIC ::   ln_ref_apr   !: ref. pressure: global mean Patm (F) or a constant (F)
-   REAL(wp)        ::   rn_pref      !  reference atmospheric pressure   [N/m2]
+   REAL(wp),PUBLIC ::   rn_pref      !  reference atmospheric pressure   [N/m2]
 
    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   ssh_ib    ! Inverse barometer now    sea surface height   [m]
@@ -34,6 +35,6 @@
    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   apr       ! atmospheric pressure at kt                 [N/m2]
    
-   REAL(wp) ::   tarea                ! whole domain mean masked ocean surface
-   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0)
+   REAL(wp), PUBLIC ::   tarea                ! whole domain mean masked ocean surface
+   REAL(wp), PUBLIC ::   r1_grau              ! = 1.e0 / (grav * rau0)
    
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_apr   ! structure of input fields (file informations, fields read)
@@ -85,10 +86,12 @@
          IF(lwm) WRITE ( numond, namsbc_apr )
          !
-         ALLOCATE( sf_apr(1), STAT=ierror )           !* allocate and fill sf_sst (forcing structure) with sn_sst
-         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_apr: unable to allocate sf_apr structure' )
-         !
-         CALL fld_fill( sf_apr, (/ sn_apr /), cn_dir, 'sbc_apr', 'Atmospheric pressure ', 'namsbc_apr' )
-                                ALLOCATE( sf_apr(1)%fnow(jpi,jpj,1)   )
-         IF( sn_apr%ln_tint )   ALLOCATE( sf_apr(1)%fdta(jpi,jpj,1,2) )
+         IF( .NOT. cpl_mslp ) THEN 
+            ALLOCATE( sf_apr(1), STAT=ierror )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_apr: unable to allocate sf_apr structure' ) 
+            ! 
+            CALL fld_fill( sf_apr, (/ sn_apr /), cn_dir, 'sbc_apr', 'Atmospheric pressure ', 'namsbc_apr' ) 
+                                   ALLOCATE( sf_apr(1)%fnow(jpi,jpj,1)   ) 
+            IF( sn_apr%ln_tint )   ALLOCATE( sf_apr(1)%fdta(jpi,jpj,1,2) ) 
+         ENDIF
                                 ALLOCATE( ssh_ib(jpi,jpj) , ssh_ibb(jpi,jpj) )
                                 ALLOCATE( apr (jpi,jpj) )
@@ -96,5 +99,9 @@
          IF(lwp) THEN                                 !* control print
             WRITE(numout,*)
-            WRITE(numout,*) '   Namelist namsbc_apr : Atmospheric PRessure as extrenal forcing'
+            IF( cpl_mslp ) THEN 
+               WRITE(numout,*) '   Namelist namsbc_apr : Atmospheric Pressure as extrenal coupling' 
+            ELSE 
+               WRITE(numout,*) '   Namelist namsbc_apr : Atmospheric Pressure as extrenal forcing' 
+            ENDIF
             WRITE(numout,*) '      ref. pressure: global mean Patm (T) or a constant (F)  ln_ref_apr = ', ln_ref_apr
          ENDIF
@@ -119,42 +126,43 @@
       ENDIF
 
-      !                                         ! ========================== !
-      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    At each sbc time-step   !
-         !                                      ! ===========+++============ !
-         !
-         IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields
-         !
-         CALL fld_read( kt, nn_fsbc, sf_apr )               !* input Patm provided at kt + nn_fsbc/2
-         !
-         !                                                  !* update the reference atmospheric pressure (if necessary)
-         IF( ln_ref_apr )   rn_pref = glob_sum( sf_apr(1)%fnow(:,:,1) * e1e2t(:,:) ) / tarea
-         !
-         !                                                  !* Patm related forcing at kt
-         ssh_ib(:,:) = - ( sf_apr(1)%fnow(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer)
-         apr   (:,:) =     sf_apr(1)%fnow(:,:,1)                        ! atmospheric pressure
-         !
-         CALL iom_put( "ssh_ib", ssh_ib )                   !* output the inverse barometer ssh
-      ENDIF
-
-      !                                         ! ---------------------------------------- !
-      IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    !
-         !                                      ! ---------------------------------------- !
-         !                                            !* Restart: read in restart file
-         IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN 
-            IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb read in the restart file'
-            CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh
+      IF( .NOT. cpl_mslp ) THEN                    ! ========================== ! 
+         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    At each sbc time-step   ! 
+            !                                      ! ===========+++============ !
             !
-         ELSE                                         !* no restart: set from nit000 values
-            IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb set to nit000 values'
-            ssh_ibb(:,:) = ssh_ib(:,:)
+            IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields 
+            ! 
+            CALL fld_read( kt, nn_fsbc, sf_apr )               !* input Patm provided at kt + nn_fsbc/2 
+            ! 
+            !                                                  !* update the reference atmospheric pressure (if necessary) 
+            IF( ln_ref_apr )   rn_pref = glob_sum( sf_apr(1)%fnow(:,:,1) * e1e2t(:,:) ) / tarea 
+            ! 
+            !                                                  !* Patm related forcing at kt 
+            ssh_ib(:,:) = - ( sf_apr(1)%fnow(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer) 
+            apr   (:,:) =     sf_apr(1)%fnow(:,:,1)                        ! atmospheric pressure 
+            ! 
+            CALL iom_put( "ssh_ib", ssh_ib )                   !* output the inverse barometer ssh
          ENDIF
-      ENDIF
-      !                                         ! ---------------------------------------- !
-      IF( lrst_oce ) THEN                       !      Write in the ocean restart file     !
-         !                                      ! ---------------------------------------- !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'sbc_apr : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp
-         IF(lwp) WRITE(numout,*) '~~~~'
-         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib )
+ 
+         !                                         ! ---------------------------------------- ! 
+         IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    ! 
+            !                                      ! ---------------------------------------- ! 
+            !                                            !* Restart: read in restart file 
+            IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN  
+               IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb read in the restart file' 
+               CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh 
+               ! 
+            ELSE                                         !* no restart: set from nit000 values 
+               IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb set to nit000 values' 
+               ssh_ibb(:,:) = ssh_ib(:,:) 
+            ENDIF 
+         ENDIF 
+         !                                         ! ---------------------------------------- ! 
+         IF( lrst_oce ) THEN                       !      Write in the ocean restart file     ! 
+            !                                      ! ---------------------------------------- ! 
+            IF(lwp) WRITE(numout,*) 
+            IF(lwp) WRITE(numout,*) 'sbc_apr : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp 
+            IF(lwp) WRITE(numout,*) '~~~~' 
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib ) 
+         ENDIF
       ENDIF
       !
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 12448)
@@ -319,4 +319,5 @@
          &               Cd, Ch, Ce, zt_zu, zq_zu )
     
+      IF ( ln_cdgw .AND. nn_drag==jp_std ) Cd(:,:) = cdn_wave(:,:)
       ! ... tau module, i and j component
       DO jj = 1, jpj
@@ -737,5 +738,5 @@
 
       !! Neutral coefficients at 10m:
-      IF( ln_cdgw ) THEN      ! wave drag case
+      IF( ln_cdgw .AND. nn_drag==jp_mcore ) THEN      ! wave drag case
          cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
          ztmp0   (:,:) = cdn_wave(:,:)
@@ -783,5 +784,5 @@
          END IF
        
-         IF( ln_cdgw ) THEN      ! surface wave case
+         IF( ln_cdgw .AND. nn_drag==jp_mcore ) THEN      ! surface wave case
             sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
             Cd      = sqrt_Cd * sqrt_Cd
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 12448)
@@ -24,5 +24,4 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
-   USE sbcwave,ONLY : cdn_wave !wave module
 
    IMPLICIT NONE
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 12448)
@@ -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
@@ -41,4 +42,7 @@
    USE timing          ! Timing
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+#if defined key_oasis3 || defined key_oasis3mct   
+   USE mod_oasis                    ! OASIS3-MCT module   
+#endif  
    USE eosbn2
    USE sbcrnf   , ONLY : l_rnfcpl
@@ -105,5 +109,17 @@
    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_tauoc  = 50            ! Stress fraction adsorbed by waves
+   INTEGER, PARAMETER ::   jpr_wdrag  = 51            ! Neutral surface drag coefficient  
+   INTEGER, PARAMETER ::   jpr_wfreq  = 52            ! Wave peak frequency  
+   INTEGER, PARAMETER ::   jpr_tauwx  = 53            ! x component of the ocean stress from waves 
+   INTEGER, PARAMETER ::   jpr_tauwy  = 54            ! y component of the ocean stress from waves 
+   INTEGER, PARAMETER ::   jprcv      = 54            ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
@@ -135,5 +151,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 +169,11 @@
    ! 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_sdrft,sn_rcv_wper, & 
+                    sn_rcv_wfreq,sn_rcv_wnum,sn_rcv_tauoc,sn_rcv_tauw, & 
+                    sn_rcv_wdrag
    ! Other namelist parameters                        !
    INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -161,5 +187,5 @@
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
-
+    
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
 
@@ -179,5 +205,5 @@
       !!             ***  FUNCTION sbc_cpl_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(3)
+      INTEGER :: ierr(4)
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -188,6 +214,10 @@
       ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init)
 #endif
-      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
-      !
+!      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
+      ! Hardwire three models as nn_cplmodel has not been read in from the namelist yet.   
+      ALLOCATE( xcplmask(jpi,jpj,0:3) , 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 +246,11 @@
       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_sdrft, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wfreq, sn_rcv_tauoc,    & 
+         &                  sn_rcv_wdrag, sn_rcv_qns   , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal ,     & 
+         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp , sn_rcv_tauw ,                  & 
+         &                  nn_cplmodel, ln_usecplmask, ln_usernfmask
       !!---------------------------------------------------------------------
       !
@@ -260,8 +293,19 @@
          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,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%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 u,v        = ', TRIM(sn_rcv_sdrft%cldes ), ' (', TRIM(sn_rcv_sdrft%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,*)'      Wave peak frequency             = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')'   
+         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_tauoc%cldes ), ' (', TRIM(sn_rcv_tauoc%clcat ), ')'   
+         WRITE(numout,*)'      Stress components by waves      = ', TRIM(sn_rcv_tauw%cldes  ), ' (', TRIM(sn_rcv_tauw%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,10 +313,16 @@
          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,*)'      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
+         WRITE(numout,*)'  ln_usernfmask                       = ', ln_usernfmask
       ENDIF
 
       !                                   ! allocate sbccpl arrays
-      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
+      !IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
      
       ! ================================ !
@@ -307,4 +357,5 @@
       ! 
       ! Vectors: change of sign at north fold ONLY if on the local grid
+      IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled
       IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
       
@@ -337,5 +388,9 @@
          srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
          srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
-         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
+         !srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2  
+         ! Currently needed for HadGEM3 - but shouldn't affect anyone else for the moment  
+         srcv(jpr_otx1)%laction = .TRUE.   
+         srcv(jpr_oty1)%laction = .TRUE.  
+         !  
          srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
       CASE( 'T,I' ) 
@@ -374,4 +429,5 @@
          srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
       ENDIF
+      ENDIF
        
       !                                                      ! ------------------------- !
@@ -403,5 +459,5 @@
       IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
          srcv(jpr_rnf)%laction = .TRUE.
-         l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf
+         l_rnfcpl              = .NOT. ln_usernfmask         ! -> no need to read runoffs in sbcrnf
          ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
          IF(lwp) WRITE(numout,*)
@@ -470,4 +526,14 @@
       !                                                      ! ------------------------- !
       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' ) THEN 
+         srcv(jpr_mslp)%laction = .TRUE. 
+         cpl_mslp = .TRUE. 
+      ENDIF
+
       !                                                      ! ------------------------- !
       !                                                      !   topmelt and botmelt     !   
@@ -483,4 +549,62 @@
          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  
+      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction  
+      IF( TRIM(sn_rcv_sdrft%cldes ) == 'coupled' )  THEN
+         srcv(jpr_sdrftx)%laction = .TRUE.  
+         srcv(jpr_sdrfty)%laction = .TRUE.  
+         cpl_sdrft = .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_wfreq)%clname = 'O_WFreq'     ! wave peak frequency  
+      IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' )  THEN  
+         srcv(jpr_wfreq)%laction = .TRUE.  
+         cpl_wfreq = .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_tauoc)%clname = 'O_TauOce'     ! stress fraction adsorbed by the wave  
+      IF( TRIM(sn_rcv_tauoc%cldes ) == 'coupled' )  THEN  
+         srcv(jpr_tauoc)%laction = .TRUE.  
+         cpl_tauoc = .TRUE.  
+      ENDIF  
+      srcv(jpr_tauwx)%clname = 'O_Tauwx'      ! ocean stress from wave in the x direction 
+      srcv(jpr_tauwy)%clname = 'O_Tauwy'      ! ocean stress from wave in the y direction 
+      IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' )  THEN  
+         srcv(jpr_tauwx)%laction = .TRUE.  
+         srcv(jpr_tauwy)%laction = .TRUE.  
+         cpl_tauw = .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  
+      ! 
+      IF( srcv(jpr_tauoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) & 
+            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 
+                                     '(sn_rcv_tauoc=coupled and sn_rcv_tauw=coupled)' ) 
+      !
+      !   
       !                                                      ! ------------------------------- !
       !                                                      !   OPA-SAS coupling - rcv by opa !   
@@ -637,4 +761,5 @@
       !                                                      ! ------------------------- !
       ssnd(jps_fice)%clname = 'OIceFrc'
+      ssnd(jps_ficet)%clname = 'OIceFrcT'
       ssnd(jps_hice)%clname = 'OIceTck'
       ssnd(jps_hsnw)%clname = 'OSnwTck'
@@ -645,4 +770,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 +792,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 +814,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.
 
       !                                                      ! ------------------------------- !
@@ -783,5 +933,5 @@
       IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
          &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
-      ncpl_qsr_freq = 86400 / ncpl_qsr_freq
+      IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq
 
       CALL wrk_dealloc( jpi,jpj, zacs, zaos )
@@ -837,4 +987,8 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
+      USE sbcflx ,  ONLY : ln_shelf_flx 
+      USE sbcssm ,  ONLY : sbc_ssm_cpl 
+      USE lib_fortran     ! distributed memory computing library
+
       INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
       INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
@@ -845,4 +999,5 @@
       INTEGER  ::   ji, jj, jn             ! dummy loop indices
       INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
+      INTEGER  ::   ikchoix
       REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
       REAL(wp) ::   zcoef                  ! temporary scalar
@@ -850,10 +1005,10 @@
       REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
       REAL(wp) ::   zzx, zzy               ! temporary variables
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, ztx2, zty2, zmsk, zemp, zqns, zqsr
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
       !
-      CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
+      CALL wrk_alloc( jpi,jpj, ztx, zty, ztx2, zty2, zmsk, zemp, zqns, zqsr )
       !
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -893,12 +1048,34 @@
             IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
                !                                                       ! (geographical to local grid -> rotate the components)
-               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
-               IF( srcv(jpr_otx2)%laction ) THEN
-                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
-               ELSE	
-                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )  
+               IF( srcv(jpr_otx1)%clgrid == 'U' .AND. (.NOT. srcv(jpr_otx2)%laction) ) THEN   
+                  ! Temporary code for HadGEM3 - will be removed eventually.   
+                  ! Only applies when we have only taux on U grid and tauy on V grid   
+                  DO jj=2,jpjm1   
+                     DO ji=2,jpim1   
+                        ztx(ji,jj)=0.25*vmask(ji,jj,1) &   
+                           *(frcv(jpr_otx1)%z3(ji,jj,1)+frcv(jpr_otx1)%z3(ji-1,jj,1)    &   
+                           +frcv(jpr_otx1)%z3(ji,jj+1,1)+frcv(jpr_otx1)%z3(ji-1,jj+1,1))   
+                        zty(ji,jj)=0.25*umask(ji,jj,1) &   
+                           *(frcv(jpr_oty1)%z3(ji,jj,1)+frcv(jpr_oty1)%z3(ji+1,jj,1)    &   
+                           +frcv(jpr_oty1)%z3(ji,jj-1,1)+frcv(jpr_oty1)%z3(ji+1,jj-1,1))   
+                     ENDDO   
+                  ENDDO   
+                              
+                  ikchoix = 1   
+                  CALL repcmo(frcv(jpr_otx1)%z3(:,:,1),zty,ztx,frcv(jpr_oty1)%z3(:,:,1),ztx2,zty2,ikchoix)   
+                  CALL lbc_lnk (ztx2,'U', -1. )   
+                  CALL lbc_lnk (zty2,'V', -1. )   
+                  frcv(jpr_otx1)%z3(:,:,1)=ztx2(:,:)   
+                  frcv(jpr_oty1)%z3(:,:,1)=zty2(:,:)   
+               ELSE   
+                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )      
+                  frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid   
+                  IF( srcv(jpr_otx2)%laction ) THEN   
+                     CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )      
+                  ELSE   
+                     CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )    
+                  ENDIF   
+                  frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
                ENDIF
-               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
-               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
             ENDIF
             !                              
@@ -919,7 +1096,26 @@
       ELSE                                                   !   No dynamical coupling   !
          !                                                   ! ========================= !
+         ! it is possible that the momentum is calculated from the winds (ln_shelf_flx) and a coupled drag coefficient 
+         IF( srcv(jpr_wdrag)%laction .AND. ln_shelf_flx .AND. ln_cdgw .AND. nn_drag == jp_std ) THEN 
+            DO jj = 1, jpjm1 
+               DO ji = 1, jpim1 
+                  ! here utau and vtau should contain the wind components as read from the forcing files, 
+                  ! and the wind module should already be properly calculated 
+                  frcv(jpr_otx1)%z3(ji,jj,1) = zrhoa * 0.5 * ( frcv(jpr_wdrag)%z3(ji,jj,1) + frcv(jpr_wdrag)%z3(ji+1,jj,1) ) * & 
+                                                                         utau(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) ) 
+                  frcv(jpr_oty1)%z3(ji,jj,1) = zrhoa * 0.5 * ( frcv(jpr_wdrag)%z3(ji,jj,1) + frcv(jpr_wdrag)%z3(ji,jj+1,1) ) * & 
+                                                                         vtau(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) ) 
+                  utau(ji,jj) = frcv(jpr_otx1)%z3(ji,jj,1) 
+                  vtau(ji,jj) = frcv(jpr_oty1)%z3(ji,jj,1) 
+               END DO 
+            END DO 
+            CALL lbc_lnk_multi( frcv(jpr_otx1)%z3(:,:,1), 'U', -1. , frcv(jpr_oty1)%z3(:,:,1), 'V', -1. , & 
+                                                             utau(:,:), 'U', -1. , vtau(:,:), 'V',  -1. ) 
+            llnewtx = .TRUE. 
+         ELSE
          frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero 
          frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
          llnewtx = .TRUE.
+         ENDIF
          !
       ENDIF
@@ -941,4 +1137,7 @@
             END DO
             CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
+            IF( .NOT. srcv(jpr_otx1)%laction .AND. srcv(jpr_wdrag)%laction .AND. & 
+                                ln_shelf_flx .AND. ln_cdgw .AND. nn_drag == jp_std ) & 
+               taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
             llnewtau = .TRUE.
          ELSE
@@ -955,4 +1154,5 @@
       !                                                      ! ========================= !
       !                                                      !      10 m wind speed      !   (wndm)
+      !                                                      !   include wave drag coef  !   (wndm)
       !                                                      ! ========================= !
       !
@@ -965,5 +1165,9 @@
 !CDIR NOVERRCHK
                DO ji = 1, jpi 
+                  IF( ln_shelf_flx ) THEN   ! the 10 wind module is properly calculated before if ln_shelf_flx 
+                     frcv(jpr_w10m)%z3(ji,jj,1) = wndm(ji,jj) 
+                  ELSE
                   frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
+                  ENDIF
                END DO
             END DO
@@ -975,10 +1179,15 @@
       IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
          !
+         ! if ln_wavcpl, the fields already contain the right information from forcing even if not ln_mixcpl
          IF( ln_mixcpl ) THEN
-            utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)
-            vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)
-            taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)
-            wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)
-         ELSE
+            IF( srcv(jpr_otx1)%laction ) THEN 
+               utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:) 
+               vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:) 
+            ENDIF 
+            IF( srcv(jpr_taum)%laction )   & 
+               taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:) 
+            IF( srcv(jpr_w10m)%laction )   & 
+               wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:) 
+         ELSE IF( ll_purecpl ) THEN
             utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
             vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
@@ -996,4 +1205,111 @@
       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( ln_apr_dyn ) THEN
+             IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields   
+           
+             !                                                  !* update the reference atmospheric pressure (if necessary) 
+             IF( ln_ref_apr )  rn_pref = glob_sum( frcv(jpr_mslp)%z3(:,:,1) * e1e2t(:,:) ) / tarea 
+           
+             ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer)   
+             apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1) !atmospheric pressure   
+             ! 
+             CALL iom_put( "ssh_ib", ssh_ib )                                    !* output the inverse barometer ssh
+           
+             !                                         ! ---------------------------------------- ! 
+             IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    ! 
+                !                                      ! ---------------------------------------- ! 
+                !* Restart: read in restart file 
+                IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN 
+                   IF(lwp) WRITE(numout,*) 'sbc_cpl:   ssh_ibb read in the restart file' 
+                   CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh 
+                ELSE                                         !* no restart: set from nit000 values 
+                   IF(lwp) WRITE(numout,*) 'sbc_cpl:   ssh_ibb set to nit000 values' 
+                   ssh_ibb(:,:) = ssh_ib(:,:) 
+                ENDIF 
+             ENDIF 
+             !                                         ! ---------------------------------------- ! 
+             IF( lrst_oce ) THEN                       !      Write in the ocean restart file     ! 
+                !                                      ! ---------------------------------------- ! 
+                IF(lwp) WRITE(numout,*) 
+                IF(lwp) WRITE(numout,*) 'sbc_cpl : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp 
+                IF(lwp) WRITE(numout,*) '~~~~' 
+                CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib ) 
+             ENDIF 
+          ENDIF 
+        
+          ! Update mean ssh 
+          CALL sbc_ssm_cpl( kt )
+      END IF   
+      !  
+      IF( ln_sdw ) THEN  ! Stokes Drift correction activated  
+      !                                                      ! ========================= !   
+      !                                                      !       Stokes drift u,v    !  
+      !                                                      ! ========================= !   
+      IF( srcv(jpr_sdrftx)%laction .AND. srcv(jpr_sdrfty)%laction ) THEN 
+                                     ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)  
+                                     vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)  
+      ENDIF
+      !  
+      !                                                      ! ========================= !   
+      !                                                      !      Wave mean period     !  
+      !                                                      ! ========================= !   
+         IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)  
+      !  
+      !                                                      ! ========================= !   
+      !                                                      !  Significant wave height  !  
+      !                                                      ! ========================= !   
+         IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)  
+      !  
+      !                                                      ! ========================= !   
+      !                                                      !    Wave peak frequency    !  
+      !                                                      ! ========================= !   
+         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%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 .AND. srcv(jpr_sdrfty)%laction) .OR. srcv(jpr_wper)%laction &  
+                                        .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) &
+            CALL sbc_stokes()  
+      ENDIF  
+      !                                                      ! ========================= !   
+      !                                                      ! Stress adsorbed by waves  !  
+      !                                                      ! ========================= !   
+      IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) THEN 
+         tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1) 
+         ! cap the value of tauoc 
+         WHERE(tauoc_wave <   0.0 ) tauoc_wave = 1.0 
+         WHERE(tauoc_wave > 100.0 ) tauoc_wave = 1.0 
+      ENDIF 
+      !                                                      ! ========================= !   
+      !                                                      ! Stress component by waves !  
+      !                                                      ! ========================= !   
+      IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN 
+         tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1) 
+         tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1) 
+         ! cap the value of tauoc 
+         WHERE(tauw_x < -100.0 ) tauw_x = 0.0 
+         WHERE(tauw_x >  100.0 ) tauw_x = 0.0 
+         WHERE(tauw_y < -100.0 ) tauw_y = 0.0 
+         WHERE(tauw_y >  100.0 ) tauw_y = 0.0 
+      ENDIF
+       
+      !                                                      ! ========================= !   
+      !                                                      !   Wave to ocean energy    !
+      !                                                      ! ========================= !   
+      IF( srcv(jpr_phioc)%laction .AND. ln_phioc ) THEN 
+         rn_crban(:,:) = 29.0 * frcv(jpr_phioc)%z3(:,:,1) 
+         WHERE( rn_crban <    0.0 ) rn_crban = 0.0 
+         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0 
+      ENDIF 
 
       !  Fields received by SAS when OASIS coupling
@@ -1067,5 +1383,5 @@
                CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
             END SELECT
-         ELSE
+         ELSE IF( ll_purecpl ) THEN
             zemp(:,:) = 0._wp
          ENDIF
@@ -1075,6 +1391,7 @@
          IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
          
-         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
-         ELSE                   ;   emp(:,:) =                              zemp(:,:)
+         IF( ln_mixcpl .AND. ( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction )) THEN 
+                                         emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:) 
+         ELSE IF( ll_purecpl ) THEN  ;   emp(:,:) = zemp(:,:)
          ENDIF
          !
@@ -1091,6 +1408,7 @@
             ENDIF
          ENDIF
-         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
-         ELSE                   ;   qns(:,:) =                              zqns(:,:)
+         IF( ln_mixcpl .AND. ( srcv(jpr_qnsoce)%laction .OR. srcv(jpr_qnsmix)%laction )) THEN 
+                                          qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:) 
+         ELSE IF( ll_purecpl ) THEN   ;   qns(:,:) = zqns(:,:)
          ENDIF
 
@@ -1101,6 +1419,7 @@
          ENDIF
          IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle
-         IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)
-         ELSE                   ;   qsr(:,:) =                              zqsr(:,:)
+         IF( ln_mixcpl .AND. ( srcv(jpr_qsroce)%laction .OR. srcv(jpr_qsrmix)%laction )) THEN 
+                                          qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:) 
+         ELSE IF( ll_purecpl ) THEN   ;   qsr(:,:) = zqsr(:,:)
          ENDIF
          !
@@ -1113,5 +1432,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
+      CALL wrk_dealloc( jpi,jpj, ztx, zty, ztx2, zty2, zmsk, zemp, zqns, zqsr )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
@@ -1708,4 +2027,5 @@
       !
       INTEGER ::   ji, jj, jl   ! dummy loop indices
+      INTEGER ::   ikchoix
       INTEGER ::   isec, info   ! local integer
       REAL(wp) ::   zumax, zvmax
@@ -1736,5 +2056,5 @@
             !
             SELECT CASE( sn_snd_temp%cldes)
-            CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
+            CASE( 'oce only'             )   ;   ztmp1(:,:) = (ztmp1(:,:) + rt0) * tmask(:,:,1)
             CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
                SELECT CASE( sn_snd_temp%clcat )
@@ -1765,4 +2085,5 @@
                   ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
                ENDDO
+            CASE( 'none'         )       ! nothing to do
             CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
             END SELECT
@@ -1889,5 +2210,5 @@
          !                                                  j+1   j     -----V---F
          ! surface velocity always sent from T point                     !       |
-         !                                                        j      |   T   U
+         ! [except for HadGEM3]                                   j      |   T   U
          !                                                               |       |
          !                                                   j    j-1   -I-------|
@@ -1901,10 +2222,20 @@
             SELECT CASE( TRIM( sn_snd_crt%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
+               IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN   
+                  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   
+               ELSE   
+                  ! Temporarily Changed for UKV   
+                  DO jj = 2, jpjm1   
+                     DO ji = 2, jpim1   
+                        zotx1(ji,jj) = un(ji,jj,1)   
+                        zoty1(ji,jj) = vn(ji,jj,1)   
+                     END DO   
+                  END DO   
+               ENDIF 
             CASE( 'weighted oce and ice' )   
                SELECT CASE ( cp_ice_msh )
@@ -1930,14 +2261,28 @@
                   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
+                  IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN   
+                     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   
+#if defined key_cice   
+                  ELSE   
+                     ! Temporarily Changed for HadGEM3   
+                     DO jj = 2, jpjm1   
+                        DO ji = 2, jpim1   ! NO vector opt.   
+                           zotx1(ji,jj) = (1.0-fr_iu(ji,jj)) * un(ji,jj,1)             &   
+                                &              + fr_iu(ji,jj) * 0.5 * ( u_ice(ji,jj-1) + u_ice(ji,jj) )    
+                           zoty1(ji,jj) = (1.0-fr_iv(ji,jj)) * vn(ji,jj,1)             &   
+                                &              + fr_iv(ji,jj) * 0.5 * ( v_ice(ji-1,jj) + v_ice(ji,jj) )    
+                        END DO   
+                     END DO   
+#endif   
+                  ENDIF 
                END SELECT
                CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
@@ -1984,14 +2329,38 @@
          IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
             !                                                                     ! Ocean component
-            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component 
-            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%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
+             IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN   
+                CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component   
+                CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%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   
+             ELSE   
+                ! Temporary code for HadGEM3 - will be removed eventually.   
+                ! Only applies when we want uvel on U grid and vvel on V grid   
+                ! Rotate U and V onto geographic grid before sending.   
+             
+                DO jj=2,jpjm1   
+                   DO ji=2,jpim1   
+                      ztmp1(ji,jj)=0.25*vmask(ji,jj,1)      &   
+                           *(zotx1(ji,jj)+zotx1(ji-1,jj)    &   
+                           +zotx1(ji,jj+1)+zotx1(ji-1,jj+1))   
+                      ztmp2(ji,jj)=0.25*umask(ji,jj,1)      &   
+                           *(zoty1(ji,jj)+zoty1(ji+1,jj)    &   
+                           +zoty1(ji,jj-1)+zoty1(ji+1,jj-1))   
+                   ENDDO   
+                ENDDO   
+                   
+                ! Ensure any N fold and wrap columns are updated   
+                CALL lbc_lnk(ztmp1, 'V', -1.0)   
+                CALL lbc_lnk(ztmp2, 'U', -1.0)   
+                   
+                ikchoix = -1   
+                CALL repcmo(zotx1,ztmp2,ztmp1,zoty1,zotx1,zoty1,ikchoix)   
+            ENDIF   
          ENDIF
          !
@@ -2019,4 +2388,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/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 12448)
@@ -20,4 +20,5 @@
    USE iom             ! IOM library
    USE in_out_manager  ! I/O manager
+   USE sbcwave         ! wave physics
    USE lib_mpp         ! distribued memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -87,4 +88,5 @@
       REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
+      REAL(wp) ::   totwind               ! UKMO SHELF: Module of wind speed
       REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
       REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface
@@ -115,4 +117,13 @@
          IF(lwm) WRITE ( numond, namsbc_flx ) 
          !
+         IF(lwp) THEN                        ! Namelist print 
+            WRITE(numout,*) 
+            WRITE(numout,*) 'sbc_flx : Flux forcing' 
+            WRITE(numout,*) '~~~~~~~~~~~' 
+            WRITE(numout,*) '       Namelist namsbc_flx : shelf seas configuration (force with winds instead of momentum)' 
+            WRITE(numout,*) '          shelf seas configuration    ln_shelf_flx    = ', ln_shelf_flx 
+            WRITE(numout,*) '          relative wind speed         ln_rel_wind     = ', ln_rel_wind 
+            WRITE(numout,*) '          wind multiplication factor  rn_wfac         = ', rn_wfac 
+         ENDIF
          !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
          IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   &
@@ -210,4 +221,48 @@
             END DO
          END DO
+         !                                                        ! add modification due to drag coefficient read from wave forcing 
+         !                                                        ! this code is inefficient but put here to allow merging with another UKMO branch 
+         IF( ln_shelf_flx ) THEN 
+            ! calculate first the wind module, as it will be used later 
+            DO jj = 2, jpjm1 
+               DO ji = fs_2, fs_jpim1    ! vect. opt. 
+                  ztx = zwnd_i(ji-1,jj  ) + zwnd_i(ji,jj) 
+                  zty = zwnd_j(ji  ,jj-1) + zwnd_j(ji,jj) 
+                  wndm(ji,jj) = 0.5 * SQRT( ztx * ztx + zty * zty ) 
+               END DO 
+            END DO 
+            CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
+        
+            IF( ln_cdgw .AND. nn_drag == jp_std ) THEN 
+               IF( cpl_wdrag ) THEN  
+                  ! reset utau and vtau to the wind components: the momentum will 
+                  ! be calculated from the coupled value of the drag coefficient 
+                  DO jj = 1, jpj 
+                     DO ji = 1, jpi 
+                        utau(ji,jj) = zwnd_i(ji,jj) 
+                        vtau(ji,jj) = zwnd_j(ji,jj) 
+                     END DO 
+                  END DO 
+               ELSE 
+                  DO jj = 1, jpjm1 
+                     DO ji = 1, jpim1 
+                        utau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji+1,jj) ) * zwnd_i(ji,jj) * & 
+                                                                        0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) ) 
+                        vtau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji,jj+1) ) * zwnd_j(ji,jj) * & 
+                                                                        0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) ) 
+                     END DO 
+                  END DO 
+                  CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. ) 
+               ENDIF 
+            ELSE IF( nn_drag == jp_const ) THEN 
+               DO jj = 1, jpjm1 
+                  DO ji = 1, jpim1 
+                     utau(ji,jj) = zrhoa * zcdrag * zwnd_i(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) ) 
+                     vtau(ji,jj) = zrhoa * zcdrag * zwnd_j(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) ) 
+                  END DO 
+               END DO 
+               CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. ) 
+            ENDIF 
+         ENDIF
          !                                                        ! add to qns the heat due to e-p
          qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
@@ -230,5 +285,7 @@
                zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
                taum(ji,jj) = zmod
+               IF ( .NOT. (ln_shelf_flx .AND. ln_cpl)) THEN
                wndm(ji,jj) = SQRT( zmod * zcoef )
+               ENDIF
             END DO
          END DO
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 12448)
@@ -284,4 +284,7 @@
       CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
       !
+      ! In coupled mode get extra fields from CICE for passing back to atmosphere   
+      IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(nit000)   
+      !    
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_init')
       !
@@ -708,9 +711,4 @@
       IF( nn_timing == 1 )  CALL timing_start('cice_sbc_hadgam')
       !
-      IF( kt == nit000 )  THEN
-         IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
-         IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
-      ENDIF
-
       !                                         ! =========================== !
       !                                         !   Prepare Coupling fields   !
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 12448)
@@ -88,9 +88,8 @@
       NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   &
          &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
-         &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
-         &             nn_lsm    , nn_limflx , nn_components, ln_cpl
+         &             ln_ssr    , nn_isf    , nn_fwb, ln_wave, nn_lsm     , nn_limflx,   & 
+                       nn_components, ln_cpl , ln_wavcpl, nn_drag
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
-      LOGICAL  ::   ll_purecpl
       !!----------------------------------------------------------------------
 
@@ -131,8 +130,11 @@
          WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
          WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl
-         WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl
+         WRITE(numout,*) '              forced-coupled atm mixed formulation       ln_mixcpl   = ', ln_mixcpl 
+         WRITE(numout,*) '              forced-coupled wav mixed formulation       ln_wavcpl   = ', ln_wavcpl
+         WRITE(numout,*) '              wave physics                               ln_wave     = ', ln_wave  
          WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis
          WRITE(numout,*) '              components of your executable              nn_components = ', nn_components
          WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx
+         WRITE(numout,*) '              momentum formulation                       nn_drag     = ', nn_drag
          WRITE(numout,*) '           Misc. options of sbc : '
          WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
@@ -164,14 +166,14 @@
       IF ( nn_components == jp_iam_opa .AND. ln_cpl )   &
          &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' )
-      IF ( nn_components == jp_iam_opa .AND. ln_mixcpl )   &
-         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
+      IF ( nn_components == jp_iam_opa .AND. ( ln_mixcpl .OR. ln_wavcpl) )   & 
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl or ln_wavcpl = T in OPA' )
       IF ( ln_cpl .AND. .NOT. lk_oasis )    &
          &      CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' )
-      IF( ln_mixcpl .AND. .NOT. lk_oasis )    &
-         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' )
-      IF( ln_mixcpl .AND. .NOT. ln_cpl )    &
-         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires ln_cpl = T' )
-      IF( ln_mixcpl .AND. nn_components /= jp_iam_nemo )    &
-         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) is not yet working with sas-opa coupling via oasis' )
+      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. .NOT. lk_oasis )    & 
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) requires the cpp key key_oasis3' ) 
+      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. .NOT. ln_cpl )    & 
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) requires ln_cpl = T' ) 
+      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. nn_components /= jp_iam_nemo )    & 
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) is not yet working with sas-opa coupling via oasis' )
 
       !                              ! allocate sbc arrays
@@ -214,20 +216,6 @@
       IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa )   &
          &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
-      
-      IF ( ln_wave ) THEN
-      !Activated wave module but neither drag nor stokes drift activated
-         IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN
-            CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
-      !drag coefficient read from wave model definable only with mfs bulk formulae and core 
-         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
-             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
-         ENDIF
-      ELSE
-      IF ( ln_cdgw .OR. ln_sdw  )                                         & 
-         &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     &
-         & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ')
-      ENDIF 
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
-      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl
+      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl .AND. .NOT. ln_wavcpl
       !
       icpt = 0
@@ -261,12 +249,15 @@
          IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation'
          IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis'
-         IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed formulation'
+         IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed atm formulation' 
+         IF( ln_wavcpl          )   WRITE(numout,*) '              + forced-coupled mixed wav formulation'
          IF( nn_components/= jp_iam_nemo )  &
             &                       WRITE(numout,*) '              + OASIS coupled SAS'
       ENDIF
       !
-      IF( lk_oasis )   CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step
-      !                                                     !                                            (2) the use of nn_fsbc
-
+      IF( lk_oasis ) THEN   
+         IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )            
+         CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step   
+                                      ! (2) the use of nn_fsbc   
+      ENDIF   
 !     nn_fsbc initialization if OPA-SAS coupling via OASIS
 !     sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly
@@ -305,4 +296,7 @@
 
       IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation
+      !  
+      IF( ln_wave          )   CALL sbc_wave_init              ! surface wave initialisation  
+      !
       
    END SUBROUTINE sbc_init
@@ -325,4 +319,6 @@
       !!              - updte the ice fraction : fr_i
       !!----------------------------------------------------------------------
+      USE bdydta, ONLY: bdy_dta  
+      !
       INTEGER, INTENT(in) ::   kt       ! ocean time step
       !!---------------------------------------------------------------------
@@ -345,7 +341,6 @@
       !                                            ! ---------------------------------------- !
       !
-      IF ( .NOT. lk_bdy ) then
-         IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
-      ENDIF
+
+      IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
                                                          ! (caution called before sbc_ssm)
       !
@@ -382,6 +377,9 @@
       END SELECT
 
-      IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
-
+      IF( ln_mixcpl .OR. ln_wavcpl )  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing 
+      
+      IF( ln_wave .AND. (ln_tauoc .OR. ln_tauw) ) CALL sbc_stress( )    ! Wave stress update  
+      IF( lk_bdy )           CALL bdy_dta ( kt, time_offset=+1 )        ! update dynamic & tracer data at open boundaries 
+                                                                        ! (caution called after sbc_ssm[_cpl] and before ice)
 
       !                                            !==  Misc. Options  ==!
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 12448)
@@ -82,5 +82,6 @@
       ALLOCATE( rnfmsk(jpi,jpj)         , rnfmsk_z(jpk)          ,     &
          &      h_rnf (jpi,jpj)         , nk_rnf  (jpi,jpj)      ,     &
-         &      rnf_tsc_b(jpi,jpj,jpts) , rnf_tsc (jpi,jpj,jpts) , STAT=sbc_rnf_alloc )
+         &      rnf_tsc_b(jpi,jpj,jpts) , rnf_tsc (jpi,jpj,jpts) ,     & 
+         &      xrnfmask(jpi,jpj,1)     , STAT=sbc_rnf_alloc )
          !
       IF( lk_mpp            )   CALL mpp_sum ( sbc_rnf_alloc )
@@ -128,5 +129,6 @@
       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
          !
-         IF( .NOT. l_rnfcpl )   rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )       ! updated runoff value at time step kt
+         IF( .NOT. l_rnfcpl )   & 
+             rnf(:,:) = rnf(:,:) * (1. - xrnfmask(:,:,1)) + rn_rfact * sf_rnf(1)%fnow(:,:,1) * xrnfmask(:,:,1)  ! updated runoff value at time step kt
          !
          !                                                     ! set temperature & salinity content of runoffs
@@ -442,4 +444,13 @@
       ENDIF
       !
+      xrnfmask(:,:,:) = 1.    ! default value for points using river forcing 
+      IF (ln_usernfmask) THEN 
+         IF(lwp) WRITE(numout,*) 
+         IF(lwp) WRITE(numout,*) '          runoff mask read in a file' 
+         CALL iom_open( 'rnfmask', inum ) 
+         CALL iom_get( inum, jpdom_data, 'rnfmask', xrnfmask(:,:,1), 1) 
+         CALL iom_close( inum ) 
+      ENDIF 
+      !
       rnf(:,:) =  0._wp                         ! runoff initialisation
       rnf_tsc(:,:,:) = 0._wp                    ! runoffs temperature & salinty contents initilisation
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 12448)
@@ -26,4 +26,5 @@
 
    PUBLIC   sbc_ssm         ! routine called by step.F90
+   PUBLIC   sbc_ssm_cpl     ! routine called by sbccpl.F90
    PUBLIC   sbc_ssm_init    ! routine called by sbcmod.F90
 
@@ -77,6 +78,8 @@
          sss_m(:,:) = zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = sshn(:,:)
+         IF( .NOT. cpl_mslp ) THEN 
+            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
+            ELSE                    ;   ssh_m(:,:) = sshn(:,:) 
+            ENDIF
          ENDIF
          !
@@ -99,6 +102,8 @@
             sss_m(:,:) = zcoef * zts(:,:,jp_sal)
             !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
-            ELSE                    ;   ssh_m(:,:) = zcoef * sshn(:,:)
+            IF( .NOT. cpl_mslp ) THEN 
+               IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) ) 
+               ELSE                    ;   ssh_m(:,:) = zcoef * sshn(:,:) 
+               ENDIF
             ENDIF
             !
@@ -113,5 +118,5 @@
             sst_m(:,:) = 0.e0
             sss_m(:,:) = 0.e0
-            ssh_m(:,:) = 0.e0
+            IF( .NOT. cpl_mslp ) ssh_m(:,:) = 0.e0
             IF( lk_vvl )   e3t_m(:,:) = 0.e0
             frq_m(:,:) = 0.e0
@@ -127,6 +132,8 @@
          sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)
+         IF( .NOT. cpl_mslp ) THEN 
+            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
+            ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) 
+            ENDIF
          ENDIF
          !
@@ -143,5 +150,5 @@
             ssu_m(:,:) = ssu_m(:,:) * zcoef           ! mean suface current  [m/s]
             ssv_m(:,:) = ssv_m(:,:) * zcoef           !
-            ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m]
+            IF( .NOT. cpl_mslp ) ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m]
             IF( lk_vvl )   e3t_m(:,:) = fse3t_m(:,:) * zcoef   ! mean vertical scale factor [m]
             frq_m(:,:) = frq_m(:,:) * zcoef   ! mean fraction of solar net radiation absorbed in the 1st T level [-]
@@ -161,5 +168,5 @@
             CALL iom_rstput( kt, nitrst, numrow, 'sst_m'  , sst_m  )
             CALL iom_rstput( kt, nitrst, numrow, 'sss_m'  , sss_m  )
-            CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  )
+            IF( .NOT. cpl_mslp ) CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  )
             IF( lk_vvl )   CALL iom_rstput( kt, nitrst, numrow, 'e3t_m'  , e3t_m  )
             CALL iom_rstput( kt, nitrst, numrow, 'frq_m'  , frq_m  )
@@ -174,5 +181,5 @@
          CALL iom_put( 'sst_m', sst_m )
          CALL iom_put( 'sss_m', sss_m )
-         CALL iom_put( 'ssh_m', ssh_m )
+         IF( .NOT. cpl_mslp ) CALL iom_put( 'ssh_m', ssh_m )
          IF( lk_vvl )   CALL iom_put( 'e3t_m', e3t_m )
          CALL iom_put( 'frq_m', frq_m )
@@ -180,4 +187,67 @@
       !
    END SUBROUTINE sbc_ssm
+
+   SUBROUTINE sbc_ssm_cpl( kt ) 
+      !!--------------------------------------------------------------------- 
+      !!                   ***  ROUTINE sbc_ssm_cpl  *** 
+      !!                      
+      !! ** Purpose :   provide ocean surface variable to sea-surface boundary 
+      !!                condition computation when pressure is read from coupling 
+      !!                 
+      !! ** Method  :   The inverse barometer ssh (i.e. ssh associated with Patm) 
+      !!                is added to ssh_m when ln_apr_dyn = T. Required for sea-ice dynamics. 
+      !!--------------------------------------------------------------------- 
+      INTEGER, INTENT(in) ::   kt   ! ocean time step 
+      ! 
+      REAL(wp) ::   zcoef       ! local scalar 
+      !!--------------------------------------------------------------------- 
+      ! 
+      IF( nn_fsbc == 1 ) THEN                             !      Instantaneous surface fields        ! 
+         !                                                ! ---------------------------------------- ! 
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
+         ELSE                    ;   ssh_m(:,:) = sshn(:,:) 
+         ENDIF 
+      ELSE 
+         !                                                ! ----------------------------------------------- ! 
+         IF( kt == nit000 .AND. .NOT. l_ssm_mean ) THEN   !   Initialisation: 1st time-step, no input means ! 
+            !                                             ! ----------------------------------------------- ! 
+            IF(lwp) WRITE(numout,*) 
+            IF(lwp) WRITE(numout,*) '~~~~~~~   mean ssh field initialised to instantaneous values' 
+            zcoef = REAL( nn_fsbc - 1, wp ) 
+            zcoef = REAL( nn_fsbc - 1, wp ) 
+            IF( ln_apr_dyn ) THEN    ;  ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) ) 
+            ELSE                     ;  ssh_m(:,:) = zcoef * sshn(:,:) 
+            ENDIF 
+            !                                             ! ---------------------------------------- ! 
+         ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   ! 
+            !                                             ! ---------------------------------------- ! 
+            ssh_m(:,:) = 0.e0 
+         ENDIF 
+   
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
+         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) 
+         ENDIF 
+         !                                                ! ---------------------------------------- ! 
+         IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN          !   Mean value at each nn_fsbc time-step   ! 
+            !                                             ! ---------------------------------------- ! 
+            zcoef = 1. / REAL( nn_fsbc, wp ) 
+            ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH [m] 
+         ENDIF 
+         !                                                ! ---------------------------------------- ! 
+         IF( lrst_oce ) THEN                              !      Write in the ocean restart file     ! 
+            !                                             ! ---------------------------------------- ! 
+            IF(lwp) WRITE(numout,*) 
+            IF(lwp) WRITE(numout,*) 'sbc_ssm_cpl : ssh mean field written in ocean restart file ',   & 
+               &                    'at it= ', kt,' date= ', ndastp 
+            IF(lwp) WRITE(numout,*) '~~~~~~~' 
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  ) 
+         ENDIF 
+      ENDIF 
+      ! 
+      IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN          !   Mean value at each nn_fsbc time-step   ! 
+         CALL iom_put( 'ssh_m', ssh_m ) 
+      ENDIF 
+      ! 
+   END SUBROUTINE sbc_ssm_cpl
 
    SUBROUTINE sbc_ssm_init
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 12448)
@@ -4,37 +4,78 @@
    !! Wave module 
    !!======================================================================
-   !! History :  3.3.1  !   2011-09  (Adani M)  Original code: Drag Coefficient 
-   !!         :  3.4    !   2012-10  (Adani M)                 Stokes Drift 
+   !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient 
+   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift 
+   !!            3.6  !   2014-09  (Clementi E, Oddo P)New Stokes Drift Computation
+   !!             -   !   2016-12  (G. Madec, E. Clementi) update Stoke drift computation 
+   !!                                                     + add sbc_wave_ini routine
    !!----------------------------------------------------------------------
-   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 oce
-   USE sbc_oce	       ! Surface boundary condition: ocean fields
-   USE domvvl
-
-   
+
    !!----------------------------------------------------------------------
-   !!   sbc_wave       : read drag coefficient from wave model in netcdf files 
+   !!   sbc_stokes    : calculate 3D Stokes-drift velocities
+   !!   sbc_wave      : wave data from wave model in netcdf files 
+   !!   sbc_wave_init : initialisation fo surface waves
    !!----------------------------------------------------------------------
+   USE oce            ! ocean variables
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE bdy_oce        ! open boundary condition variables
+   USE domvvl         ! domain: variable volume layers
+   !
+   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 wrk_nemo       !
+   USE phycst         ! physical constants 
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
+   PUBLIC   sbc_stokes      ! routine called in sbccpl 
+   PUBLIC   sbc_stress      ! routine called in sbcmod
+   PUBLIC   sbc_wave        ! routine called in sbcmod 
+   PUBLIC   sbc_wave_init   ! routine called in sbcmod
    
-   INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift
-   INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_wn  = 3           ! index of wave number                 (1/m)    at T-point
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd	  ! structure of input fields (file informations, fields read) Drag Coefficient
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd	  ! structure of input fields (file informations, fields read) Stokes Drift
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
-   REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum 
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d 
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
+   ! 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_sdrft  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tauoc  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE.
+
+   INTEGER ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift
+   INTEGER, PARAMETER ::   jp_breivik  = 0     ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]
+   INTEGER, PARAMETER ::   jp_phillips = 1     ! Phillips:     v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]
+   INTEGER, PARAMETER ::   jp_peakph   = 2     ! Phillips using the peak wave number read from wave model instead of the inverse depth scale
+
+   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_hsw   ! index of significant wave hight      (m)      at T-point
+   INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point
+   INTEGER ::   jp_wfr   ! index of wave peak frequency         (s^-1)   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
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauoc ! 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
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_phioc ! structure of input fields (file informations, fields read) wave to ocean energy 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wfreq               !: 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   rn_crban            !: Craig and Banner constant for surface breaking waves mixing
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_x              !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_y              !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd        !: surface Stokes drift velocities at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp.
+
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
@@ -44,151 +85,524 @@
 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   ! dummy loop argument
+      INTEGER  ::   ik           ! local integer 
+      REAL(wp) ::  ztransp, zfac, ztemp
+      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v
+      REAL(wp), DIMENSION(:,:)  , POINTER ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd   ! 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 ) 
+      !
+      ! 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) )
+         zfac = 2.0_wp * rpi / 16.0_wp
+         DO jj = 1, jpj               
+            DO ji = 1, jpi
+               ! Stokes drift velocity estimated from Hs and Tmean
+               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
+               ! Stokes surface speed
+               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
+               ! Wavenumber scale
+               zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
+            END DO
+         END DO
+         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
+            DO ji = 1, jpim1
+               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
+               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
+               !
+               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
+               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
+            END DO
+         END DO
+      ELSE IF( nn_sdrift==jp_peakph ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
+         DO jj = 1, jpjm1              
+            DO ji = 1, jpim1
+               zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav
+               zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav
+               !
+               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
+               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
+            END DO
+         END DO
+      ENDIF
+      !
+      !                       !==  horizontal Stokes Drift 3D velocity  ==!
+      IF( nn_sdrift==jp_breivik ) THEN
+         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 ) / ( 1.0_wp + 8.0_wp*zkh_u )
+                  zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
+                  !
+                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
+                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph ) THEN
+         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))
+                  !
+                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
+                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. )
+      !
+      !                       !==  vertical Stokes Drift 3D velocity  ==!
+      !
+      DO jk = 1, jpkm1               ! Horizontal e3*divergence
+         DO jj = 2, jpj
+            DO ji = fs_2, jpi
+               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd(ji,  jj,jk)    &
+                  &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd(ji-1,jj,jk)    &
+                  &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd(ji,jj  ,jk)    &
+                  &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd(ji,jj-1,jk)  ) * r1_e12t(ji,jj)
+            END DO
+         END DO
+      END DO
+      !
+      IF( .NOT. AGRIF_Root() ) THEN
+         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3divh(nlci-1,   :  ,:) = 0._wp      ! east
+         IF( nbondi == -1 .OR. nbondi == 2 )   ze3divh(  2   ,   :  ,:) = 0._wp      ! west
+         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3divh(  :   ,nlcj-1,:) = 0._wp      ! north
+         IF( nbondj == -1 .OR. nbondj == 2 )   ze3divh(  :   ,  2   ,:) = 0._wp      ! south
+      ENDIF
+      !
+      CALL lbc_lnk( ze3divh, 'T', 1. )
+      !
+      IF( .NOT. lk_vvl ) THEN   ;   ik = 1   ! none zero velocity through the sea surface
+      ELSE                      ;   ik = 2   ! w=0 at the surface (set one for all in sbc_wave_init)
+      ENDIF
+      DO jk = jpkm1, ik, -1          ! integrate from the bottom the hor. divergence (NB: at k=jpk w is always zero)
+         wsd(:,:,jk) = wsd(:,:,jk+1) - ze3divh(:,:,jk)
+      END DO
+#if defined key_bdy
+      IF( lk_bdy ) THEN
+         DO jk = 1, jpkm1
+            wsd(:,:,jk) = wsd(:,:,jk) * bdytmask(:,:)
+         END DO
+      ENDIF
+#endif
+      !                       !==  Horizontal divergence of barotropic Stokes transport  ==!
+      div_sd(:,:) = 0._wp
+      DO jk = 1, jpkm1                                 ! 
+        div_sd(:,:) = div_sd(:,:) + ze3divh(:,:,jk)
+      END DO
+      !
+      CALL iom_put( "ustokes",  usd  )
+      CALL iom_put( "vstokes",  vsd  )
+      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 )
+      !
+   END SUBROUTINE sbc_stokes
+
+
+   SUBROUTINE sbc_stress( )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_stress  ***
+      !!
+      !! ** Purpose :   Updates the ocean momentum modified by waves
+      !!
+      !! ** Method  : - Calculate u,v components of stress depending on stress
+      !!                model 
+      !!              - Calculate the stress module
+      !!              - The wind module is not modified by waves 
+      !! ** action  
+      !!---------------------------------------------------------------------
+      INTEGER  ::   jj, ji   ! dummy loop argument
+      !
+      IF( ln_tauoc ) THEN
+         utau(:,:) = utau(:,:)*tauoc_wave(:,:)
+         vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
+         taum(:,:) = taum(:,:)*tauoc_wave(:,:)
+      ENDIF
+      !
+      IF( ln_tauw ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               ! Stress components at u- & v-points
+               utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) )
+               vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) )
+               !
+               ! Stress module at t points
+               taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,: ), 'T', -1. )
+      ENDIF
+      !
+   END SUBROUTINE sbc_stress
+
+
    SUBROUTINE sbc_wave( kt )
       !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE sbc_apr  ***
-      !!
-      !! ** Purpose :   read drag coefficient from wave model  in netcdf files.
+      !!                     ***  ROUTINE sbc_wave  ***
+      !!
+      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
       !!
       !! ** Method  : - Read namelist namsbc_wave
       !!              - Read Cd_n10 fields in netcdf files 
       !!              - Read stokes drift 2d in netcdf files 
-      !!              - Read wave number      in netcdf files 
-      !!              - Compute 3d stokes drift using monochromatic
-      !! ** action  :   
-      !!               
-      !!---------------------------------------------------------------------
-      USE oce,  ONLY : un,vn,hdivn,rotn
-      USE divcur
-      USE wrk_nemo
-#if defined key_bdy
-      USE bdy_oce, ONLY : bdytmask
-#endif
-      INTEGER, INTENT( in  ) ::  kt       ! ocean time step
-      INTEGER                ::  ierror   ! return error code
-      INTEGER                ::  ifpr, jj,ji,jk 
-      INTEGER                ::   ios     ! Local integer output status for namelist read
-      REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy
-      REAL                                          ::  z2dt,z1_2dt
-      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
+      !!              - Read wave number in netcdf files 
+      !!              - Compute 3d stokes drift using Breivik et al.,2014
+      !!                formulation
+      !! ** action  
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
+      !!---------------------------------------------------------------------
+      !
+      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)
+         ! check that the drag coefficient contains proper information even if
+         ! the masks do not match - the momentum stress is not masked!
+         WHERE( cdn_wave < 0.0 ) cdn_wave = 1.5e-3
+         WHERE( cdn_wave > 1.0 ) cdn_wave = 1.5e-3
+      ENDIF
+
+      IF( ln_tauoc .AND. .NOT. cpl_tauoc ) 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)
+         WHERE( tauoc_wave <   0.0 ) tauoc_wave = 1.0
+         WHERE( tauoc_wave > 100.0 ) tauoc_wave = 1.0
+      ENDIF
+
+      IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN      !==  Wave induced stress  ==!
+         CALL fld_read( kt, nn_fsbc, sf_tauw )           ! read ocean stress components from external forcing (T grid)
+         tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1)
+         WHERE( tauw_x < -100.0 ) tauw_x = 0.0
+         WHERE( tauw_x >  100.0 ) tauw_x = 0.0
+
+         tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1)
+         WHERE( tauw_y < -100.0 ) tauw_y = 0.0
+         WHERE( tauw_y >  100.0 ) tauw_y = 0.0
+      ENDIF
+
+      IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN    !==  Wave to ocean energy  ==!
+         CALL fld_read( kt, nn_fsbc, sf_phioc )          ! read wave to ocean energy from external forcing
+         rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1)     ! ! Alfa is phioc*sqrt(rau0/zrhoa)  : rau0=water density, zhroa= air density
+         WHERE( rn_crban > 1.e8   ) rn_crban = 0.0    !remove first mask mistmatch points, then cap values in case of low friction velocity
+         WHERE( rn_crban < 0.0    ) rn_crban = 0.0
+         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0
+      ENDIF
+
+      IF( ln_sdw .OR. ln_rough )  THEN             !==  Computation of the 3d Stokes Drift  ==! 
+         !
+         IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled
+            CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing
+            IF( jp_hsw > 0 ) THEN
+               hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height
+               WHERE( hsw > 100.0 ) hsw = 0.0
+               WHERE( hsw <   0.0 ) hsw = 0.0
+            ENDIF
+            IF( jp_wmp > 0 ) THEN
+               wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
+               WHERE( wmp > 100.0 ) wmp = 0.0
+               WHERE( wmp <   0.0 ) wmp = 0.0
+            ENDIF
+            IF( jp_wfr > 0 ) THEN
+               wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1)   ! Peak wave frequency 
+               WHERE( wfreq <    0.0 ) wfreq = 0.0 
+               WHERE( wfreq >  100.0 ) wfreq = 0.0
+            ENDIF
+            IF( jp_usd > 0 ) THEN
+               ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
+               WHERE( ut0sd < -100.0 ) ut0sd = 1.0
+               WHERE( ut0sd >  100.0 ) ut0sd = 1.0
+            ENDIF
+            IF( jp_vsd > 0 ) THEN
+               vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
+               WHERE( vt0sd < -100.0 ) vt0sd = 1.0
+               WHERE( vt0sd >  100.0 ) vt0sd = 1.0
+            ENDIF
+         ENDIF
+      ENDIF
+      !
+      IF( ln_sdw ) THEN
+         ! 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  ==! 
+         !
+         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_peakph .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
+         !
+      ENDIF
+      !
+   END SUBROUTINE sbc_wave
+
+
+   SUBROUTINE sbc_wave_init
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_wave_init  ***
+      !!
+      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
+      !!
+      !! ** Method  : - Read namelist namsbc_wave
+      !!              - Read Cd_n10 fields in netcdf files 
+      !!              - Read stokes drift 2d in netcdf files 
+      !!              - Read wave number in netcdf files 
+      !!              - Compute 3d stokes drift using Breivik et al.,2014
+      !!                formulation
+      !! ** action  
+      !!---------------------------------------------------------------------
+      INTEGER ::   ierror, ios   ! local integer
+      INTEGER ::   ifpr
+      !!
       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
-      !!---------------------------------------------------------------------
-
-      !!----------------------------------------------------------------------
-      !
-      !
-      !                                         ! -------------------- !
-      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
-         !                                      ! -------------------- !
-         REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model
-         READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
-901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp )
-
-         REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model
-         READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
-902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp )
-         IF(lwm) WRITE ( numond, namsbc_wave )
-         !
-
-         IF ( ln_cdgw ) THEN
+      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_phioc, &
+                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum , &
+                             &   sn_tauoc, sn_tauwx, sn_tauwy      ! information about the fields to be read
+      !
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_phioc,          &
+                             sn_tauoc, sn_tauwx, sn_tauwy,                                                       &
+                             ln_cdgw, ln_sdw, ln_stcor, ln_phioc, ln_tauoc, ln_tauw, ln_zdfqiao, ln_rough,       &
+                             nn_sdrift, nn_wmix
+      !!---------------------------------------------------------------------
+      !
+      REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model
+      READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp )
+         
+      REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model
+      READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namsbc_wave )
+      !
+      IF(lwp) THEN               !* Parameter print
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_wave_init: wave physics'
+         WRITE(numout,*) '~~~~~~~~'
+         WRITE(numout,*) '   Namelist namsbc_wave : set wave physics parameters'
+         WRITE(numout,*) '      Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw 
+         WRITE(numout,*) '        vertical parametrization              nn_sdrift   = ', nn_sdrift 
+         WRITE(numout,*) '      Stokes coriolis term                    ln_stcor    = ', ln_stcor 
+         WRITE(numout,*) '      wave modified ocean stress              ln_tauoc    = ', ln_tauoc 
+         WRITE(numout,*) '      wave modified ocean stress components   ln_tauw     = ', ln_tauw 
+         WRITE(numout,*) '      wave to ocean energy                    ln_phioc    = ', ln_phioc
+         WRITE(numout,*) '        vertical mixing parametrization       nn_wmix     = ', nn_wmix 
+         WRITE(numout,*) '      neutral drag coefficient                ln_cdgw     = ', ln_cdgw
+         WRITE(numout,*) '      wave roughness length modification      ln_rough    = ', ln_rough 
+         WRITE(numout,*) '      Qiao vertical mixing formulation        ln_zdfqiao  = ', ln_zdfqiao
+      ENDIF
+
+      IF ( ln_wave ) THEN
+         ! Activated wave physics but no wave physics components activated 
+         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_tauw .OR. ln_stcor .OR. ln_phioc &
+                                                                    .OR. ln_rough .OR. ln_zdfqiao) )   THEN  
+             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_tauw=F, ln_stcor=F ', &
+                                                      'ln_phioc=F, ln_rough=F, ln_zdfqiao=F' ) 
+         ELSE
+         IF (ln_stcor .AND. .NOT. ln_sdw) & 
+             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
+         IF ( ln_cdgw .AND. .NOT.(nn_drag==jp_ukmo .OR. nn_drag==jp_std .OR. nn_drag==jp_const .OR. nn_drag==jp_mcore) ) & 
+             CALL ctl_stop( 'The chosen nn_drag for momentum calculation must be 0, 1, 2, or 3')
+         IF ( ln_cdgw .AND. ln_blk_core .AND. nn_drag==0 ) &
+             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in core forcing must be 1, 2, or 3')
+         IF ( ln_cdgw .AND. ln_flx .AND. nn_drag==3 ) &
+             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in direct forcing must be 0, 1, or 2')
+         IF( ln_phioc .AND. .NOT.(nn_wmix==jp_craigbanner .OR. nn_wmix==jp_janssen) ) & 
+            CALL ctl_stop( 'The chosen nn_wmix for wave vertical mixing must be 0, or 1' )
+         IF( ln_sdw .AND. .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph) ) & 
+            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
+         IF( ln_zdfqiao .AND. .NOT.ln_sdw ) & 
+            CALL ctl_stop( 'Qiao vertical mixing can not be used without Stokes drift (ln_sdw)' )
+         IF( ln_tauoc .AND. ln_tauw ) & 
+            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
+                                     '(ln_tauoc=.true. and ln_tauw=.true.)' )
+         IF( ln_tauoc ) &
+             CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauoc=.true.)' )
+         IF( ln_tauw ) &
+             CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', &
+                                  'This will override any other specification of the ocean stress' ) 
+         ENDIF
+      ELSE
+         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_tauw .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao ) &  
+            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
+            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
+            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
+            &                  'or Stokes-Coriolis term (ln_stcor=T)',                         &
+            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      & 
+            &                  'or ocean stress components from waves (ln_tauw=T) ',          & 
+            &                  'or wave to ocean energy modification (ln_phioc=T) ',           & 
+            &                  'or wave surface roughness (ln_rough=T) ',                      & 
+            &                  'or Qiao vertical mixing formulation (ln_zdfqiao=T) ' )
+      ENDIF 
+      !
+      IF( ln_cdgw ) THEN
+         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' )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_cd structure' )
             !
                                    ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
             IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( cdn_wave(jpi,jpj) )
-            cdn_wave(:,:) = 0.0
-        ENDIF
-         IF ( ln_sdw ) THEN
-            slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn
-            ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
+            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
+         ENDIF
+         ALLOCATE( cdn_wave(jpi,jpj) )
+      ENDIF
+
+      IF( ln_tauoc ) THEN
+         IF( .NOT. cpl_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_init: unable to allocate sf_tauoc structure' )
             !
-            DO ifpr= 1, jpfld
-               ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
-               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
-            END DO
-            CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) )
-            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
-            usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0
-            usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0
-         ENDIF
-      ENDIF
+                                    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_init', 'read wave input', 'namsbc_wave' )
+         ENDIF
+         ALLOCATE( tauoc_wave(jpi,jpj) )
+      ENDIF
+
+      IF( ln_tauw ) THEN
+         IF( .NOT. cpl_tauw ) THEN
+            ALLOCATE( sf_tauw(2), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwx/y
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' )
+            !
+            ALLOCATE( slf_j(2) )
+            slf_j(1) = sn_tauwx
+            slf_j(2) = sn_tauwy
+                                    ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1)   )
+                                    ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1)   )
+            IF( slf_j(1)%ln_tint )  ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) )
+            IF( slf_j(2)%ln_tint )  ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) )
+            CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
+         ENDIF
+         ALLOCATE( tauw_x(jpi,jpj) )
+         ALLOCATE( tauw_y(jpi,jpj) )
+      ENDIF
+
+      IF( ln_phioc ) THEN
+         IF( .NOT. cpl_phioc ) THEN
+            ALLOCATE( sf_phioc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_phioc
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_phioc structure' )
+            !
+                                    ALLOCATE( sf_phioc(1)%fnow(jpi,jpj,1)   )
+            IF( sn_phioc%ln_tint )  ALLOCATE( sf_phioc(1)%fdta(jpi,jpj,1,2) )
+            CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
+         ENDIF
+         ALLOCATE( rn_crban(jpi,jpj) )
+      ENDIF
+
+      ! Find out how many fields have to be read from file if not coupled
+      jpfld=0
+      jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0
+      IF( ln_sdw ) THEN
+         IF( .NOT. cpl_sdrft ) THEN
+            jpfld  = jpfld + 1
+            jp_usd = jpfld
+            jpfld  = jpfld + 1
+            jp_vsd = jpfld
+         ENDIF
+         IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
+            jpfld  = jpfld + 1
+            jp_hsw = jpfld
+         ENDIF
+         IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
+            jpfld  = jpfld + 1
+            jp_wmp = jpfld
+         ENDIF
+         IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakph ) THEN
+            jpfld  = jpfld + 1
+            jp_wfr = jpfld
+         ENDIF
+      ENDIF
+
+      IF( ln_rough .AND. .NOT. cpl_hsig .AND. jp_hsw==0 ) THEN
+         jpfld  = jpfld + 1
+         jp_hsw = 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_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
+         IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
+         IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr
+         ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd 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
          !
-      IF ( ln_cdgw ) THEN
-         CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing
-         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
-      ENDIF
-      IF ( ln_sdw )  THEN
-          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing
-
-         ! Interpolate wavenumber, stokes drift into the grid_V and grid_V
-         !-------------------------------------------------
-
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
-
-               vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
-
-               usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
-
-               vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
-            END DO
-         END DO
-
-          !Computation of the 3d Stokes Drift
-          DO jk = 1, jpk
-             DO jj = 1, jpj-1
-                DO ji = 1, jpi-1
-                   usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_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))))
-                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(:,:,:)
-          CALL div_cur(kt)
-      !                                           !------------------------------!
-      !                                           !     Now Vertical Velocity    !
-      !                                           !------------------------------!
-          z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
-
-          z1_2dt = 1.e0 / z2dt
-          DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
-             ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
-             wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
-                &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
-                &                         * tmask(:,:,jk) * z1_2dt
-#if defined key_bdy
-             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
-#endif
-          END DO
-          hdivn(:,:,:)=hdivdummy(:,:,:)
-          rotn(:,:,:)=rotdummy(:,:,:)
-          vn(:,:,:)=vdummy(:,:,:)
-          un(:,:,:)=udummy(:,:,:)
-          CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
-      ENDIF
-   END SUBROUTINE sbc_wave
-      
+         CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
+      ENDIF
+
+      IF( ln_sdw ) THEN
+         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) )
+         ALLOCATE( wmp  (jpi,jpj)  )
+         ALLOCATE( wfreq (jpi,jpj) )
+         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     )
+         ALLOCATE( div_sd(jpi,jpj) )
+         ALLOCATE( tsd2d (jpi,jpj) )
+         usd(:,:,:) = 0._wp
+         vsd(:,:,:) = 0._wp
+         wsd(:,:,:) = 0._wp
+         ! Wave number needed only if ln_zdfqiao=T
+         IF( ln_zdfqiao .AND. .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_init: unable toallocate sf_wn 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', 'read wave input', 'namsbc_wave' )
+         ENDIF
+         ALLOCATE( wnum(jpi,jpj) )
+      ENDIF
+
+      IF( ln_sdw .OR. ln_rough ) THEN
+         ALLOCATE( hsw  (jpi,jpj) )
+      ENDIF
+      !
+   END SUBROUTINE sbc_wave_init
+
    !!======================================================================
 END MODULE sbcwave
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 12448)
@@ -6,4 +6,5 @@
    !! History :  2.0  !  2005-11  (G. Madec)  Original code
    !!            3.3  !  2010-09  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
+   !!            3.6  !  2015-06  (E. Clementi) Addition of Stokes drift in case of wave coupling
    !!            4.0  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation
    !!----------------------------------------------------------------------
@@ -34,6 +35,7 @@
    USE timing          ! Timing
    USE sbc_oce
-   USE diaptr          ! Poleward heat transport 
-
+   USE sbcwave        ! wave module  
+   USE sbc_oce        ! surface boundary condition: ocean  
+   USE diaptr         ! Poleward heat transport
 
    IMPLICIT NONE
@@ -85,4 +87,8 @@
       CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn )
       !                                          ! set time step
+      zun(:,:,:) = 0.0  
+      zvn(:,:,:) = 0.0  
+      zwn(:,:,:) = 0.0  
+      !
       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
          r2dtra(:) =  rdttra(:)                          ! = rdtra (restarting with Euler time stepping)
@@ -94,9 +100,17 @@
       !
       !                                               !==  effective transport  ==!
-      DO jk = 1, jpkm1
-         zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
-         zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
-         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk)
-      END DO
+      IF(ln_wave .AND. ln_sdw) THEN  
+         DO jk = 1, jpkm1                                                                      ! eulerian transport + Stokes Drift  
+            zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) ) 
+            zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) )  
+            zwn(:,:,jk) = e1e2t(:,:)               * ( wn(:,:,jk) + wsd(:,:,jk) )  
+         END DO  
+      ELSE  
+         DO jk = 1, jpkm1  
+            zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)             ! eulerian transport only  
+            zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk)  
+            zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)  
+         END DO  
+      ENDIF
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 12448)
@@ -10,4 +10,5 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
+   USE sbc_oce, ONLY : ln_zdfqiao
 
    IMPLICIT NONE
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 12448)
@@ -24,4 +24,6 @@
    USE phycst         ! physical constants
    USE zdfmxl         ! mixed layer
+   USE sbcwave
+   !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! MPP manager
@@ -46,4 +48,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_tke1 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_tke3 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_psi1 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rtrans
 
    !                              !! ** Namelist  namzdf_gls  **
@@ -59,7 +66,7 @@
    REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2)
    REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
-   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
+   REAL(wp) ::   rn_crban_default  ! Craig and Banner constant for surface breaking waves mixing
    REAL(wp) ::   rn_hsro           ! Minimum surface roughness
-   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 
+   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met = 1, 3) 
 
    REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
@@ -90,8 +97,8 @@
    REAL(wp) ::   rm7           =  0.0_wp
    REAL(wp) ::   rm8           =  0.318_wp
-   REAL(wp) ::   rtrans        =  0.1_wp
+   REAL(wp) ::   rtrans_default =  0.1_wp
    REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
-   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        -
-   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        -
+   REAL(wp) ::   rsbc_tke1_default, rsbc_tke2, rfact_tke          !     -           -           -        - 
+   REAL(wp) ::   rsbc_psi2, rsbc_psi3, rfact_psi                  !     -           -           -        -
    REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        -
    REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        -
@@ -116,5 +123,7 @@
       !!----------------------------------------------------------------------
       ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
-         &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc )
+         &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     & 
+         &      rsbc_psi1(jpi,jpj), rtrans(jpi,jpj),         &
+         &      ustars2(jpi,jpj), ustarb2(jpi,jpj)                      , STAT= zdf_gls_alloc )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc )
@@ -161,4 +170,17 @@
       ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
 
+      ! variable initialization 
+      IF( ln_phioc ) THEN 
+         IF( nn_wmix==jp_janssen ) THEN 
+            rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.Dirichlet + Wave breaking  
+            rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking 
+            rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking 
+         ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+            rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf 
+            rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf 
+            rtrans(:,:) = 0.2_wp / MAX( rsmall, rsbc_tke1(:,:)**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
+         ENDIF 
+      ENDIF 
+
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
          avt (:,:,:) = avt_k (:,:,:)
@@ -197,5 +219,6 @@
          zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10)
          zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
-      !
+      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)  
+         zhsro = MAX(rn_frac_hs * hsw, rn_hsro)
       END SELECT
 
@@ -311,9 +334,74 @@
       ! ---------------------------------
       !
+      IF( ln_phioc ) THEN 
+         SELECT CASE ( nn_bc_surf ) 
+         ! 
+         CASE ( 0 )             ! Dirichlet case 
+            IF( nn_wmix==jp_janssen ) THEN 
+               ! First level 
+               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 
+               z_elem_a(:,:,1) = en(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+ 
+               ! One level below 
+               en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 
+               z_elem_a(:,:,2) = 0._wp 
+               z_elem_c(:,:,2) = 0._wp 
+               z_elem_b(:,:,2) = 1._wp 
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+               en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 
+               en(:,:,1) = MAX(en(:,:,1), rn_emin)  
+               z_elem_a(:,:,1) = en(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               !  
+               ! One level below 
+               en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:) * ((zhsro(:,:)+fsdepw(:,:,2)) & 
+                   &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
+               en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
+               z_elem_a(:,:,2) = 0._wp  
+               z_elem_c(:,:,2) = 0._wp 
+               z_elem_b(:,:,2) = 1._wp 
+               ! 
+            ENDIF 
+         CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
+            IF( nn_wmix==jp_janssen ) THEN 
+               ! Dirichlet conditions at k=1 
+               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 
+               z_elem_a(:,:,1) = en(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! at k=2, set de/dz=Fw 
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+               z_elem_a(:,:,2) = 0._wp 
+               zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 
+               en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+               ! Dirichlet conditions at k=1 
+               en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 
+               en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
+               z_elem_a(:,:,1) = en(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! at k=2, set de/dz=Fw 
+               !cbr 
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+               z_elem_a(:,:,2) = 0._wp 
+               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:)) )) 
+               zflxs(:,:)      = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * zkar(:,:) & 
+                    &                           * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
+
+               en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
+               ! 
+            ENDIF 
+         END SELECT 
+      ELSE
       SELECT CASE ( nn_bc_surf )
       !
       CASE ( 0 )             ! Dirichlet case
       ! First level
-      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
+      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp)
       en(:,:,1) = MAX(en(:,:,1), rn_emin) 
       z_elem_a(:,:,1) = en(:,:,1)
@@ -322,6 +410,5 @@
       ! 
       ! One level below
-      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) &
-          &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
+      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
       en(:,:,2) = MAX(en(:,:,2), rn_emin )
       z_elem_a(:,:,2) = 0._wp 
@@ -331,7 +418,6 @@
       !
       CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
-      !
       ! Dirichlet conditions at k=1
-      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
+      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp)
       en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
       z_elem_a(:,:,1) = en(:,:,1)
@@ -343,7 +429,6 @@
       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
       z_elem_a(:,:,2) = 0._wp
-      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) ))
-      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &
-           &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
+      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) ))
+      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
 
       en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
@@ -351,4 +436,5 @@
       !
       END SELECT
+      ENDIF ! ln_phioc
 
       ! Bottom boundary condition on tke
@@ -536,45 +622,129 @@
       !
       CASE ( 0 )             ! Dirichlet boundary conditions
-      !
-      ! Surface value
-      zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
-      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,1) = psi(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
-      !
-      ! One level below
-      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) )))
-      zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)
-      psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,2) = 0._wp
-      z_elem_c(:,:,2) = 0._wp
-      z_elem_b(:,:,2) = 1._wp
-      ! 
-      !
+         IF( ln_phioc ) THEN   ! Wave induced mixing case 
+                               ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 
+                               ! balance between the production and the 
+                               ! dissipation terms including the wave effect 
+            IF( nn_wmix==jp_janssen ) THEN 
+               ! Surface value 
+               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
+               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,1) = psi(:,:,1)       
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! One level below 
+               zex1 = (rmm*ra_sf+rnn) 
+               zex2 = (rmm*ra_sf) 
+               zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 
+               psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 
+               z_elem_a(:,:,2) = 0._wp 
+               z_elem_c(:,:,2) = 0._wp 
+               z_elem_b(:,:,2) = 1._wp 
+               !  
+               ! 
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+               ! Surface value 
+               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
+               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,1) = psi(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! One level below 
+               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) ))) 
+               zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
+               psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,2) = 0._wp 
+               z_elem_c(:,:,2) = 0._wp 
+               z_elem_b(:,:,2) = 1._wp 
+               !  
+               ! 
+            ENDIF 
+         ELSE                  ! Wave induced mixing case with default values 
+            ! Surface value 
+            zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
+            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+            z_elem_a(:,:,1) = psi(:,:,1) 
+            z_elem_c(:,:,1) = 0._wp 
+            z_elem_b(:,:,1) = 1._wp 
+            ! 
+            ! One level below 
+            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) ))) 
+            zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
+            psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+            z_elem_a(:,:,2) = 0._wp 
+            z_elem_c(:,:,2) = 0._wp 
+            z_elem_b(:,:,2) = 1._wp 
+            !  
+            ! 
+         ENDIF
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
-      !
-      ! Surface value: Dirichlet
-      zdep(:,:)       = zhsro(:,:) * rl_sf
-      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,1) = psi(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
-      !
-      ! Neumann condition at k=2
-      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-      z_elem_a(:,:,2) = 0._wp
-      !
-      ! Set psi vertical flux at the surface:
-      zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
-      zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
-      zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
-      zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &
-             & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)
-      zflxs(:,:) = zdep(:,:) * zflxs(:,:)
-      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
-
-      !   
-      !
+         IF( ln_phioc ) THEN  ! Wave induced mixing case with forced/coupled fields 
+            IF( nn_wmix==jp_janssen ) THEN 
+               ! 
+               zdep(:,:) = rl_sf * zhsro(:,:) 
+               psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,1) = psi(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! Neumann condition at k=2 
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+               z_elem_a(:,:,2) = 0._wp 
+               ! 
+               ! Set psi vertical flux at the surface: 
+               zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 
+               zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
+                  &                   * en(:,:,1)**rmm * zdep 
+               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
+               ! 
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+               ! Surface value: Dirichlet 
+               zdep(:,:)       = zhsro(:,:) * rl_sf 
+               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,1) = psi(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! Neumann condition at k=2 
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+               z_elem_a(:,:,2) = 0._wp 
+               ! 
+               ! Set psi vertical flux at the surface: 
+               zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
+               zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
+               zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * & 
+                          (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
+               zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
+                      & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
+               zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
+               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
+               !    
+            ENDIF 
+         ELSE                 ! Wave induced mixing case with default values 
+            ! Surface value: Dirichlet 
+            zdep(:,:)       = zhsro(:,:) * rl_sf 
+            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+            z_elem_a(:,:,1) = psi(:,:,1) 
+            z_elem_c(:,:,1) = 0._wp 
+            z_elem_b(:,:,1) = 1._wp 
+            ! 
+            ! Neumann condition at k=2 
+            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+            z_elem_a(:,:,2) = 0._wp 
+            ! 
+            ! Set psi vertical flux at the surface: 
+            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
+            zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
+            zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * & 
+                       (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
+            zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
+                   & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
+            zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
+            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
+            !    
+            ! 
+         ENDIF
       END SELECT
 
@@ -867,5 +1037,5 @@
       NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
          &            rn_clim_galp, ln_sigpsi, rn_hsro,      &
-         &            rn_crban, rn_charn, rn_frac_hs,        &
+         &            rn_crban_default, rn_charn, rn_frac_hs,&
          &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
          &            nn_stab_func, nn_clos
@@ -895,5 +1065,5 @@
          WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot
          WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi
-         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
+         WRITE(numout,*) '      Craig and Banner coefficient (default)        rn_crban       = ', rn_crban_default
          WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
          WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
@@ -909,9 +1079,14 @@
 
       !                                !* Check of some namelist values
-      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' )
-      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' )
-      IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' )
-      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' )
-      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' )
+      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )  
+      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )  
+      IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )  
+      IF( nn_z0_met == 3 .AND. .NOT.ln_rough ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_rough=T' )  
+      IF( nn_z0_met .NE. 3 .AND. ln_rough ) THEN 
+         CALL ctl_warn('W A R N I N G:  ln_rough=.TRUE. but nn_z0_met is not 3 - resetting nn_z0_met to 3')  
+         nn_z0_met = 3  
+      ENDIF
+      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )  
+      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
 
       SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
@@ -1085,5 +1260,5 @@
                &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
 
-      IF ( rn_crban==0._wp ) THEN
+      IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN
          rl_sf = vkarmn
       ELSE
@@ -1124,8 +1299,8 @@
       rc03  = rc02 * rc0
       rc04  = rc03 * rc0
-      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
-      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking 
-      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
-      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer 
+      rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf 
+      rsbc_tke2         = rdt * rn_crban_default / rl_sf 
+      zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
+      rtrans_default    = 0.2_wp / zcr
       rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
       rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness 
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 12448)
@@ -53,6 +53,6 @@
       INTEGER ::   ios
       !!
-      NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,   &
-         &              ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp
+      NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,  &  
+         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp
       !!----------------------------------------------------------------------
 
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 12448)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 12448)
@@ -0,0 +1,129 @@
+MODULE zdfqiao
+   !!======================================================================
+   !!                       ***  MODULE  zdfqiao  ***
+   !! Qiao module      : vertical mixing enhancement due to surface waves 
+   !!======================================================================
+   !! History :  3.6  !  2014-10  (E. Clementi)  Original code
+   !!----------------------------------------------------------------------
+   !!   zdf_qiao        : compute Qiao parameters
+   !!----------------------------------------------------------------------
+
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE zdf_oce
+   USE sbcwave         ! wave module
+   USE dom_oce
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)  
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC zdf_qiao    ! routine called in step
+
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: qbv, qbvu, qbvv
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE zdf_qiao( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE zdf_qiao ***
+      !!
+      !! ** Purpose :Compute the Qiao term (qbv) to be added to
+      !!             vertical viscosity and diffusivity coeffs.  
+      !!
+      !! ** Method  :qbv = alpha * A * Us(0) * exp (3 * k * z)
+      !!             
+      !! ** action  :Compute the Qiao wave dependent term 
+      !!             only if ln_wave=.true.
+      !!               
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in  ) ::  kt   ! ocean time step
+      !
+      INTEGER :: jj, ji, jk   ! dummy loop indices
+      !!---------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
+         IF( .NOT. ( ln_wave .AND. ln_sdw ) )   &
+            &   CALL ctl_stop ( 'Ask for wave Qiao enhanced turbulence but ln_wave   &
+            &                    and ln_sdw have to be activated')
+         IF( zdf_qiao_alloc() /= 0 )   &
+            &   CALL ctl_stop( 'STOP', 'zdf_qiao : unable to allocate arrays' )
+      ENDIF
+
+      !
+      ! Compute the Qiao term Bv (qbv) to be added to
+      ! vertical viscosity and diffusivity
+      ! qbv = alpha * A * Us(0) * exp (3 * k * z)
+      ! alpha here is set to 1
+      !---------------------------------------------------------------------------------
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               qbv(ji,jj,jk) = 1.0 * 0.353553 * hsw(ji,jj) * tsd2d(ji,jj) *           &
+            &                  EXP(3.0 * wnum(ji,jj) *                                &                     
+            &                  (-MIN( fsdepw(ji  ,jj  ,jk), fsdepw(ji+1,jj  ,jk),     &
+            &                         fsdepw(ji  ,jj+1,jk), fsdepw(ji+1,jj+1,jk))))   &
+            &                          * wmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( qbv, 'W', 1. )   ! Lateral boundary conditions
+         
+      !
+      ! Interpolate Qiao parameter qbv into the grid_U and grid_V
+      !----------------------------------------------------------
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               qbvu(ji,jj,jk) = 0.5 * wumask(ji,jj,jk)  *              &  
+            &                  ( qbv(ji,jj,jk) + qbv(ji+1,jj  ,jk) )
+               qbvv(ji,jj,jk) = 0.5 * wvmask(ji,jj,jk)  *              &
+            &                  ( qbv(ji,jj,jk) + qbv(ji  ,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+      ! 
+      CALL lbc_lnk( qbvu, 'U', 1. ) ; CALL lbc_lnk( qbvv, 'V', 1. )   ! Lateral boundary conditions
+
+      ! Enhance vertical mixing coeff.         
+      !-------------------------------
+      !
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               avmu(ji,jj,jk) = ( avmu(ji,jj,jk) + qbvu(ji,jj,jk) ) * umask(ji,jj,jk)
+               avmv(ji,jj,jk) = ( avmv(ji,jj,jk) + qbvv(ji,jj,jk) ) * vmask(ji,jj,jk)
+               avt (ji,jj,jk) = ( avt (ji,jj,jk) + qbv (ji,jj,jk) ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+   END SUBROUTINE zdf_qiao
+
+   INTEGER FUNCTION zdf_qiao_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_qiao_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( qbv(jpi,jpj,jpk), qbvu(jpi,jpj,jpk), qbvv(jpi,jpj,jpk),   &
+         &      STAT = zdf_qiao_alloc )
+      !
+      IF( lk_mpp             )  CALL mpp_sum ( zdf_qiao_alloc )
+      IF( zdf_qiao_alloc > 0 )  CALL ctl_warn('zdf_qiao_alloc: allocation of arrays failed.')
+      !
+   END FUNCTION zdf_qiao_alloc
+      
+   !!======================================================================
+END MODULE zdfqiao
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 12448)
@@ -69,4 +69,5 @@
    USE icbini          ! handle bergs, initialisation
    USE icbstp          ! handle bergs, calving, themodynamics and transport
+   USE sbccpl
    USE cpl_oasis3      ! OASIS3 coupling
    USE c1d             ! 1D configuration
@@ -175,5 +176,17 @@
             CALL stp                         ! AGRIF: time stepping
 #else
-            CALL stp( istp )                 ! standard time stepping
+            IF (lk_oasis) CALL sbc_cpl_snd( istp )  ! Coupling to atmos   
+            CALL stp( istp )   
+            ! We don't couple on the final timestep because   
+            ! our restart file has already been written   
+            ! and contains all the necessary data for a   
+            ! restart. sbc_cpl_snd could be called here   
+            ! but it would require   
+            ! a) A test to ensure it was not performed   
+            !    on the very last time-step   
+            ! b) the presence of another call to   
+            !    sbc_cpl_snd call prior to the main DO loop   
+            ! This solution produces identical results   
+            ! with fewer lines of code.  
 #endif
             istp = istp + 1
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 12448)
@@ -25,4 +25,5 @@
    !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal
    !!                 !  2012-07  (J. Simeon, G. Madec, C. Ethe) Online coarsening of outputs
+   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves
    !!            3.7  !  2014-04  (F. Roquet, G. Madec) New equations of state
    !!----------------------------------------------------------------------
@@ -72,5 +73,5 @@
       !!              -8- Outputs and diagnostics
       !!----------------------------------------------------------------------
-      INTEGER ::   jk       ! dummy loop indice
+      INTEGER ::   ji,jj,jk ! dummy loop indice
       INTEGER ::   indic    ! error indicator if < 0
       INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
@@ -105,9 +106,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( lk_tide    )   CALL sbc_tide( kstp )
-      IF( lk_bdy     )  THEN
-         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib 
-                         CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-      ENDIF
-                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
+                         CALL sbc     ( kstp )        ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -132,4 +129,6 @@
       IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
       IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
+      IF( ln_zdfqiao )   CALL zdf_qiao( kstp )           ! Qiao vertical mixing   
+      !
       IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
       IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
@@ -399,4 +398,5 @@
                                CALL ctl_stop( 'step: indic < 0' )
                                CALL dia_wri_state( 'output.abort', kstp )
+                               CALL ctl_stop('STOP','NEMO failure in stp')
       ENDIF
       IF( ln_harm_ana_store   )   CALL harm_ana( kstp )        ! Harmonic analysis of tides 
@@ -411,5 +411,5 @@
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
+      !IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
       !
 #if defined key_iomput
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 12448)
@@ -29,4 +29,5 @@
    USE sbctide          ! Tide initialisation
    USE sbcapr           ! surface boundary condition: ssh_ib required by bdydta 
+   USE sbcwave          ! Wave intialisation
 
    USE traqsr           ! solar radiation penetration      (tra_qsr routine)
@@ -86,4 +87,5 @@
    USE zdfric           ! Richardson vertical mixing       (zdf_ric routine)
    USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine)
+   USE zdfqiao          !Qiao module wave induced mixing   (zdf_qiao routine)
 
    USE zpshde           ! partial step: hor. derivative     (zps_hde routine)
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 12448)
@@ -34,5 +34,5 @@
                       jp_fabm_r6s,   jp_fabm_r8n,  &
                       jp_fabm_r8c,   jp_fabm_r8p,  &
-                      jp_fabm_r8s,                 &
+                      jp_fabm_r8s,   jp_fabm_xeps, &
                       jp_fabm_pgrow, jp_fabm_ploss
 
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 12448)
@@ -122,11 +122,7 @@
       jp_fabm_o3ph  = fabm_diag_index( 'O3_pH' )
       jp_fabm_o3pc  = fabm_diag_index( 'O3_pCO2' )
+      jp_fabm_xeps  = fabm_diag_index( 'light_xEPS' )
       jp_fabm_pgrow = fabm_diag_index( 'p_grow_sum_result' )
       jp_fabm_ploss = fabm_diag_index( 'p_loss_sum_result' )
-      
-      MLD_MAX(:,:)   = 0.0
-      PGROW_AVG(:,:) = 0.0
-      PLOSS_AVG(:,:) = 0.0
-      PHYT_AVG(:,:)  = 0.0
 
       IF (lwp) THEN
@@ -436,4 +432,11 @@
          ENDDO
       ENDIF
+      
+      ! Initialise variables required for Hemmings et al. (2008)
+      ! nitrogen balancing for chlorophyll data assimilation
+      MLD_MAX(:,:)   = 0.0
+      PGROW_AVG(:,:) = 0.0
+      PLOSS_AVG(:,:) = 0.0
+      PHYT_AVG(:,:)  = 0.0
 
    END SUBROUTINE trc_ini_fabm
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 12448)
@@ -104,5 +104,5 @@
                   c_no(2:k_floor)=abs(w_if(2:k_floor,jn))*z2dt/ &
                                 ( 0.5_wp*(fse3t(ji,jj,2:k_floor) + &
-                                fse3t(ji,ji,1:k_floor-1)) )
+                                fse3t(ji,jj,1:k_floor-1)) )
                   cmax_no=MAXVAL(c_no(2:k_floor))
 
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90	(revision 12448)
@@ -109,9 +109,14 @@
             ii = idx%nbi(ib,igrd)
             ij = idx%nbj(ib,igrd)
-            zwgt = idx%nbw(ib,igrd)
+            zwgt = idx%nbw(ib,igrd) * rdttrc(ik) / 86400.d0  ! damping with a timescale of day
             tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) + zwgt * ( ( dta%trc(ib,ik) * dta%rn_fac)  & 
                         &  - tra(ii,ij,ik,jn) ) ) * tmask(ii,ij,ik)
          END DO
       END DO 
+! tra is overwritten at the boundary so damping doesn't work here - need neumann in addition
+! to duplicate the internal value at the boundary
+       CALL bdy_trc_nmn( jn,  idx, dta, kt )
+
+
       !
       IF( kt .eq. nit000 ) CLOSE( unit = 102 )
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 12448)
@@ -81,6 +81,5 @@
       CALL top_alloc()              ! allocate TOP arrays
 
-      l_trcdm2dc = ln_dm2dc .OR. ( ln_cpl .AND. ncpl_qsr_freq /= 1 )
-      l_trcdm2dc = l_trcdm2dc  .AND. .NOT. lk_offline
+      l_trcdm2dc = ln_dm2dc  .AND. .NOT. lk_offline
       IF( l_trcdm2dc .AND. lwp ) &
          &   CALL ctl_warn(' Coupling with passive tracers and used of diurnal cycle. &
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 12448)
@@ -24,7 +24,11 @@
    !!----------------------------------------------------------------------
    USE oce_trc
+   USE dom_oce, ONLY: fjulday
+   USE phycst, ONLY: rday
+   USE in_out_manager, ONLY: ln_rstdate
    USE trc
    USE trcnam_trp
    USE iom
+   USE ioipsl, ONLY: ju2ymds
    USE daymod
    ! +++>>> FABM
@@ -51,4 +55,7 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt       ! number of iteration
+      INTEGER             ::   iyear, imonth, iday
+      REAL (wp)           ::   zsec
+      REAL (wp)           ::   zfjulday
       !
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
@@ -80,8 +87,15 @@
       ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 2*nn_dttrc + 1)
       ! except if we write tracer restart files every tracer time step or if a tracer restart file was writen at nitend - 2*nn_dttrc + 1
-      IF( kt == nitrst - 2*nn_dttrc .OR. nstock == nn_dttrc .OR. ( kt == nitend - nn_dttrc .AND. .NOT. lrst_trc ) ) THEN
-         ! beware of the format used to write kt (default is i8.8, that should be large enough)
-         IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst
-         ELSE                        ;   WRITE(clkt,'(i8.8)') nitrst
+      IF( kt == nitrst - 2*nn_dttrc + 1 .OR. nstock == nn_dttrc .OR. ( kt == nitend - nn_dttrc .AND. .NOT. lrst_trc ) ) THEN
+         if ( ln_rstdate ) THEN
+            zfjulday = fjulday + rdttrc(1) / rday
+            IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error
+            CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )
+            WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
+         ELSE
+            ! beware of the format used to write kt (default is i8.8, that should be large enough)
+            IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst
+            ELSE                        ;   WRITE(clkt,'(i8.8)') nitrst
+            ENDIF
          ENDIF
          ! create the file
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/SETTE/sette.sh
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/SETTE/sette.sh	(revision 12447)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/SETTE/sette.sh	(revision 12448)
@@ -88,6 +88,7 @@
 #
 # Compiler among those in NEMOGCM/ARCH
-COMPILER=X64_ADA
-export BATCH_COMMAND_PAR="llsubmit"
+module load cray-netcdf-hdf5parallel  
+COMPILER=XC40_METO  
+export BATCH_COMMAND_PAR="qsub
 export BATCH_COMMAND_SEQ=$BATCH_COMMAND_PAR
 export INTERACT_FLAG="no"
