Index: NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/SBC/sbccpl.F90	(revision 13360)
+++ NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/SBC/sbccpl.F90	(revision 14037)
@@ -8,4 +8,5 @@
    !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
    !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi)  wave coupling updates
    !!----------------------------------------------------------------------
 
@@ -41,5 +42,5 @@
 #endif
 #if defined key_si3
-   USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE icevar         ! for CALL ice_var_snwblow
 #endif
    !
@@ -48,4 +49,8 @@
    USE lib_mpp        ! distribued memory computing library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+
+#if defined key_oasis3 
+   USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut 
+#endif 
 
    IMPLICIT NONE
@@ -102,20 +107,26 @@
    INTEGER, PARAMETER ::   jpr_fraqsr = 42   ! fraction of solar net radiation absorbed in the first ocean level
    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 
+   !**  surface wave coupling  **
+   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_tauwoc = 50   ! Stress fraction adsorbed by waves
+   INTEGER, PARAMETER ::   jpr_wstrf  = 50   ! Stress fraction adsorbed by waves
    INTEGER, PARAMETER ::   jpr_wdrag  = 51   ! Neutral surface drag coefficient
-   INTEGER, PARAMETER ::   jpr_isf    = 52
-   INTEGER, PARAMETER ::   jpr_icb    = 53
-   INTEGER, PARAMETER ::   jpr_wfreq  = 54   ! Wave peak frequency
-   INTEGER, PARAMETER ::   jpr_tauwx  = 55   ! x component of the ocean stress from waves
-   INTEGER, PARAMETER ::   jpr_tauwy  = 56   ! y component of the ocean stress from waves
-   INTEGER, PARAMETER ::   jpr_ts_ice = 57   ! Sea ice surface temp
-
-   INTEGER, PARAMETER ::   jprcv      = 57   ! total number of fields received  
+   INTEGER, PARAMETER ::   jpr_charn  = 52   ! Chranock coefficient
+   INTEGER, PARAMETER ::   jpr_twox   = 53   ! wave to ocean momentum flux
+   INTEGER, PARAMETER ::   jpr_twoy   = 54   ! wave to ocean momentum flux
+   INTEGER, PARAMETER ::   jpr_tawx   = 55   ! net wave-supported stress
+   INTEGER, PARAMETER ::   jpr_tawy   = 56   ! net wave-supported stress
+   INTEGER, PARAMETER ::   jpr_bhd    = 57   ! Bernoulli head. waves' induced surface pressure
+   INTEGER, PARAMETER ::   jpr_tusd   = 58   ! zonal stokes transport
+   INTEGER, PARAMETER ::   jpr_tvsd   = 59   ! meridional stokes tranmport
+   INTEGER, PARAMETER ::   jpr_isf    = 60
+   INTEGER, PARAMETER ::   jpr_icb    = 61
+   INTEGER, PARAMETER ::   jpr_ts_ice = 62   ! Sea ice surface temp
+
+   INTEGER, PARAMETER ::   jprcv      = 62   ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -152,5 +163,5 @@
    INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level 
    INTEGER, PARAMETER ::   jps_fice1  = 33   ! first-order ice concentration (for semi-implicit coupling of atmos-ice fluxes)
-   INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area
+   INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area fraction
    INTEGER, PARAMETER ::   jps_ht_p   = 35   ! meltpond thickness
    INTEGER, PARAMETER ::   jps_kice   = 36   ! sea ice effective conductivity
@@ -159,4 +170,12 @@
 
    INTEGER, PARAMETER ::   jpsnd      = 38   ! total number of fields sent 
+
+#if ! defined key_oasis3 
+   ! Dummy variables to enable compilation when oasis3 is not being used 
+   INTEGER                    ::   OASIS_Sent        = -1 
+   INTEGER                    ::   OASIS_SentOut     = -1 
+   INTEGER                    ::   OASIS_ToRest      = -1 
+   INTEGER                    ::   OASIS_ToRestOut   = -1 
+#endif 
 
    !                                  !!** namelist namsbc_cpl **
@@ -172,16 +191,18 @@
       &             sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr
    !                                   ! Received from the atmosphere
-   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr,  &
+   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, sn_rcv_ts_ice
    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf
-   ! Send to waves 
+   !                                   ! Send to waves 
    TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 
-   ! Received from waves 
-   TYPE(FLD_C) ::   sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, sn_rcv_tauwoc, &
-                    sn_rcv_wdrag, sn_rcv_wfreq
+   !                                   ! Received from waves 
+   TYPE(FLD_C) ::   sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, &
+      &             sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd
    !                                   ! Other namelist parameters
    INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
    LOGICAL     ::   ln_usecplmask         !  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)
+   LOGICAL     ::   ln_scale_ice_flux     !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration) 
+
    TYPE ::   DYNARR     
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3   
@@ -191,4 +212,7 @@
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
+#if defined key_si3 || defined key_cice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple !: Ice fractional area at last coupling time
+#endif
 
    REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2] 
@@ -211,5 +235,5 @@
       !!             ***  FUNCTION sbc_cpl_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(4)
+      INTEGER :: ierr(5)
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -221,6 +245,9 @@
 #endif
       ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
-      !
-      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) ) 
+#if defined key_si3 || defined key_cice
+      ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(4) )
+#endif
+      !
+      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(5) )
 
       sbc_cpl_alloc = MAXVAL( ierr )
@@ -249,13 +276,13 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos
       !!
-      NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   & 
+      NAMELIST/namsbc_cpl/  nn_cplmodel  , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux,             &
+         &                  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2   ,  & 
          &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  & 
-         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   & 
-         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr  ,   & 
-         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_tauwoc,  &
+         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc ,  & 
+         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr   ,  & 
+         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_wstrf ,  &
+         &                  sn_rcv_charn , sn_rcv_taw   , sn_rcv_bhd  , sn_rcv_tusd  , sn_rcv_tvsd,    &
          &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal  ,   &
-         &                  sn_rcv_iceflx, sn_rcv_co2   , nn_cplmodel , ln_usecplmask, sn_rcv_mslp ,   &
-         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq , sn_rcv_tauw, nn_cats_cpl  ,   &
-         &                  sn_rcv_ts_ice
+         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_icb  , sn_rcv_isf   , sn_rcv_ts_ice 
 
       !!---------------------------------------------------------------------
@@ -278,4 +305,8 @@
       ENDIF
       IF( lwp .AND. ln_cpl ) THEN                        ! control print
+         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
+         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
+         WRITE(numout,*)'  ln_scale_ice_flux                   = ', ln_scale_ice_flux
+         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
          WRITE(numout,*)'  received fields (mutiple ice categogies)'
          WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
@@ -295,4 +326,6 @@
          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,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')'
+         WRITE(numout,*)'      surface waves:'
          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 ), ')' 
@@ -301,9 +334,7 @@
          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_tauwoc%cldes), ' (', TRIM(sn_rcv_tauwoc%clcat ), ')' 
-         WRITE(numout,*)'      Stress components by waves      = ', TRIM(sn_rcv_tauw%cldes  ), ' (', TRIM(sn_rcv_tauw%clcat  ), ')'
+         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'
          WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
-         WRITE(numout,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')' 
+         WRITE(numout,*)'      Charnock coefficient            = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')'
          WRITE(numout,*)'  sent fields (multiple ice categories)'
          WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
@@ -326,9 +357,25 @@
          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,*)'  nn_cats_cpl                         = ', nn_cats_cpl
-      ENDIF
-
+      ENDIF
+      IF( lwp .AND. ln_wave) THEN                        ! control print
+      WRITE(numout,*)'      surface waves:'
+         WRITE(numout,*)'      Significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'
+         WRITE(numout,*)'      Wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'
+         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'
+         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')'
+         WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')'
+         WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')'
+         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'
+         WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')'
+         WRITE(numout,*)'      Charnock coefficient            = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')'
+         WRITE(numout,*)' Transport associated to Stokes drift grid u = ', TRIM(sn_rcv_tusd%cldes ), ' (', TRIM(sn_rcv_tusd%clcat ), ')'
+         WRITE(numout,*)' Transport associated to Stokes drift grid v = ', TRIM(sn_rcv_tvsd%cldes ), ' (', TRIM(sn_rcv_tvsd%clcat ), ')'
+         WRITE(numout,*)'      Bernouilli pressure head        = ', TRIM(sn_rcv_bhd%cldes   ), ' (', TRIM(sn_rcv_bhd%clcat  ), ')'
+         WRITE(numout,*)'Wave to ocean momentum flux and Net wave-supported stress = ', TRIM(sn_rcv_taw%cldes ), ' (', TRIM(sn_rcv_taw%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
+      ENDIF
       !                                   ! allocate sbccpl arrays
       IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
@@ -367,5 +414,5 @@
       IF(       TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice'  &
            .OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-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.
       
@@ -608,9 +655,4 @@
          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
@@ -618,15 +660,8 @@
          cpl_wnum = .TRUE.
       ENDIF
-      srcv(jpr_tauwoc)%clname = 'O_TauOce'   ! stress fraction adsorbed by the wave
-      IF( TRIM(sn_rcv_tauwoc%cldes ) == 'coupled' )  THEN
-         srcv(jpr_tauwoc)%laction = .TRUE.
-         cpl_tauwoc = .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.
+      srcv(jpr_wstrf)%clname = 'O_WStrf'     ! stress fraction adsorbed by the wave
+      IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' )  THEN
+         srcv(jpr_wstrf)%laction = .TRUE.
+         cpl_wstrf = .TRUE.
       ENDIF
       srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient
@@ -635,7 +670,36 @@
          cpl_wdrag = .TRUE.
       ENDIF
-      IF( srcv(jpr_tauwoc)%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_tauwoc=coupled and sn_rcv_tauw=coupled)' )
+      srcv(jpr_charn)%clname = 'O_Charn'     ! Chranock coefficient
+      IF( TRIM(sn_rcv_charn%cldes ) == 'coupled' )  THEN
+         srcv(jpr_charn)%laction = .TRUE.
+         cpl_charn = .TRUE.
+      ENDIF
+      srcv(jpr_bhd)%clname = 'O_Bhd'     ! Bernoulli head. waves' induced surface pressure
+      IF( TRIM(sn_rcv_bhd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_bhd)%laction = .TRUE.
+         cpl_bhd = .TRUE.
+      ENDIF
+      srcv(jpr_tusd)%clname = 'O_Tusd'     ! zonal stokes transport
+      IF( TRIM(sn_rcv_tusd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_tusd)%laction = .TRUE.
+         cpl_tusd = .TRUE.
+      ENDIF
+      srcv(jpr_tvsd)%clname = 'O_Tvsd'     ! meridional stokes tranmport
+      IF( TRIM(sn_rcv_tvsd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_tvsd)%laction = .TRUE.
+         cpl_tvsd = .TRUE.
+      ENDIF
+
+      srcv(jpr_twox)%clname = 'O_Twox'     ! wave to ocean momentum flux in the u direction
+      srcv(jpr_twoy)%clname = 'O_Twoy'     ! wave to ocean momentum flux in the v direction
+      srcv(jpr_tawx)%clname = 'O_Tawx'     ! Net wave-supported stress in the u direction
+      srcv(jpr_tawy)%clname = 'O_Tawy'     ! Net wave-supported stress in the v direction
+      IF( TRIM(sn_rcv_taw%cldes ) == 'coupled' )  THEN
+         srcv(jpr_twox)%laction = .TRUE.
+         srcv(jpr_twoy)%laction = .TRUE.
+         srcv(jpr_tawx)%laction = .TRUE.
+         srcv(jpr_tawy)%laction = .TRUE.
+         cpl_taw = .TRUE.
+      ENDIF
       !
       !                                                      ! ------------------------------- !
@@ -698,5 +762,5 @@
          ! Change first letter to couple with atmosphere if already coupled OPA
          ! this is nedeed as each variable name used in the namcouple must be unique:
-         ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere
+         ! for example O_Runoff received by OPA from SAS and therefore S_Runoff received by SAS from the Atmosphere
          DO jn = 1, jprcv
             IF( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
@@ -822,4 +886,8 @@
       END SELECT
 
+      ! Initialise ice fractions from last coupling time to zero (needed by Met-Office)
+#if defined key_si3 || defined key_cice
+       a_i_last_couple(:,:,:) = 0._wp
+#endif
       !                                                      ! ------------------------- ! 
       !                                                      !      Ice Meltponds        ! 
@@ -1033,5 +1101,4 @@
       !   initialisation of the coupler  !
       ! ================================ !
-
       CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
       
@@ -1046,4 +1113,5 @@
       ENDIF
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
+      !
       !
    END SUBROUTINE sbc_cpl_init
@@ -1110,5 +1178,5 @@
       REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
       REAL(wp) ::   zzx, zzy               ! temporary variables
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra
       !!----------------------------------------------------------------------
       !
@@ -1121,4 +1189,8 @@
          IF( ncpl_qsr_freq /= 0) ncpl_qsr_freq = 86400 / ncpl_qsr_freq ! used by top
          
+         IF ( ln_wave .AND. nn_components == 0 ) THEN
+            ncpl_qsr_freq = 1;
+            WRITE(numout,*) 'ncpl_qsr_freq is set to 1 when coupling NEMO with wave (without SAS) '
+         ENDIF
       ENDIF
       !
@@ -1170,5 +1242,5 @@
             !                              
             IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 0, 0, 0, 0 )                                        ! T ==> (U,V)
                   frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
                   frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
@@ -1224,5 +1296,17 @@
          ENDIF
       ENDIF
-
+!!$      !                                                      ! ========================= !
+!!$      SELECT CASE( TRIM( sn_rcv_clouds%cldes ) )             !       cloud fraction      !
+!!$      !                                                      ! ========================= !
+!!$      cloud_fra(:,:) = frcv(jpr_clfra)*z3(:,:,1)
+!!$      END SELECT
+!!$
+      zcloud_fra(:,:) = pp_cldf   ! should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
+      IF( ln_mixcpl ) THEN
+         cloud_fra(:,:) = cloud_fra(:,:) * xcplmask(:,:,0) + zcloud_fra(:,:)* zmsk(:,:)
+      ELSE
+         cloud_fra(:,:) = zcloud_fra(:,:)
+      ENDIF
+      !                                                      ! ========================= !
       ! u(v)tau and taum will be modified by ice model
       ! -> need to be reset before each call of the ice/fsbc      
@@ -1283,9 +1367,4 @@
          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   !
@@ -1294,6 +1373,6 @@
 
          ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode
-         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &
-                                      .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) THEN
+         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. &
+             srcv(jpr_wper)%laction .OR. srcv(jpr_hsig)%laction )   THEN
             CALL sbc_stokes( Kmm )
          ENDIF
@@ -1302,19 +1381,48 @@
       !                                                      ! Stress adsorbed by waves  !
       !                                                      ! ========================= ! 
-      IF( srcv(jpr_tauwoc)%laction .AND. ln_tauwoc ) tauoc_wave(:,:) = frcv(jpr_tauwoc)%z3(:,:,1)
-
-      !                                                      ! ========================= !  
-      !                                                      ! 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)
-      ENDIF
-
+      IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc )  tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1)
+      !
       !                                                      ! ========================= ! 
       !                                                      !   Wave drag coefficient   !
       !                                                      ! ========================= ! 
       IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw )   cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1)
-
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Chranock coefficient    !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_charn)%laction .AND. ln_charn )  charn(:,:) = frcv(jpr_charn)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      ! net wave-supported stress !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_tawx)%laction .AND. ln_taw )     tawx(:,:) = frcv(jpr_tawx)%z3(:,:,1)
+      IF( srcv(jpr_tawy)%laction .AND. ln_taw )     tawy(:,:) = frcv(jpr_tawy)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !wave to ocean momentum flux!
+      !                                                      ! ========================= !
+      IF( srcv(jpr_twox)%laction .AND. ln_taw )     twox(:,:) = frcv(jpr_twox)%z3(:,:,1)
+      IF( srcv(jpr_twoy)%laction .AND. ln_taw )     twoy(:,:) = frcv(jpr_twoy)%z3(:,:,1)
+      !                                                      
+      !                                                      ! ========================= !
+      !                                                      !    wave TKE flux at sfc   !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_phioc)%laction .AND. ln_phioc )     phioc(:,:) = frcv(jpr_phioc)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !      Bernoulli head       !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_bhd)%laction .AND. ln_bern_srfc )   bhd_wave(:,:) = frcv(jpr_bhd)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Stokes transport u dir  !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_tusd)%laction .AND. ln_breivikFV_2016 )    tusd(:,:) = frcv(jpr_tusd)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Stokes transport v dir  !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_tvsd)%laction .AND. ln_breivikFV_2016 )     tvsd(:,:) = frcv(jpr_tvsd)%z3(:,:,1)
+      !
       !  Fields received by SAS when OASIS coupling
       !  (arrays no more filled at sbcssm stage)
@@ -1549,5 +1657,5 @@
             p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
          CASE( 'T' )
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )                    ! T ==> (U,V)
                ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology 
                zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
@@ -1623,9 +1731,10 @@
       !
       INTEGER  ::   ji, jj, jl   ! dummy loop index
-      REAL(wp) ::   ztri         ! local scalar
       REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
       REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice
       REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap_ice_total
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
       !!----------------------------------------------------------------------
       !
@@ -1647,5 +1756,4 @@
          ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
          zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
-         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * picefr(:,:)
       CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
@@ -1659,6 +1767,45 @@
 
 #if defined key_si3
+
+      ! --- evaporation over ice (kg/m2/s) --- !
+      IF (ln_scale_ice_flux) THEN ! typically met-office requirements
+         IF (sn_rcv_emp%clcat == 'yes') THEN
+            WHERE( a_i(:,:,:) > 1.e-10 )  ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+            ELSEWHERE                     ; zevap_ice(:,:,:) = 0._wp
+            END WHERE
+            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
+            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
+            END WHERE
+         ELSE
+            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * SUM( a_i_last_couple, dim=3 ) / picefr(:,:)
+            ELSEWHERE                     ; zevap_ice(:,:,1) = 0._wp
+            END WHERE
+            zevap_ice_total(:,:) = zevap_ice(:,:,1)
+            DO jl = 2, jpl
+               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
+            ENDDO
+         ENDIF
+      ELSE
+         IF (sn_rcv_emp%clcat == 'yes') THEN
+            zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl)
+            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
+            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
+            END WHERE
+         ELSE
+            zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1)
+            zevap_ice_total(:,:) = zevap_ice(:,:,1)
+            DO jl = 2, jpl
+               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
+            ENDDO
+         ENDIF
+      ENDIF
+
+      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
+         ! For conservative case zemp_ice has not been defined yet. Do it now.
+         zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:)
+      ENDIF
+
       ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing)
-      zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw )
+      zsnw(:,:) = 0._wp   ;   CALL ice_var_snwblow( ziceld, zsnw )
       
       ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
@@ -1667,11 +1814,5 @@
 
       ! --- evaporation over ocean (used later for qemp) --- !
-      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:)
-
-      ! --- evaporation over ice (kg/m2/s) --- !
-      DO jl=1,jpl
-         IF(sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl)
-         ELSE                                  ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,1 )   ;   ENDIF
-      ENDDO
+      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
 
       ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
@@ -1751,15 +1892,15 @@
 !!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
 !!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) * tmask(:,:,1)                         )  ! iceshelf
-      IF( srcv(jpr_cal)%laction )   CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving
-      IF( srcv(jpr_icb)%laction )   CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs
-      IF( iom_use('snowpre') )      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
-      IF( iom_use('precip') )       CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
-      IF( iom_use('rain') )         CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation 
-      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
-      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
-      IF( iom_use('rain_ao_cea') )  CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average)
-      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average)
-      IF( iom_use('evap_ao_cea') )  CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
-         &                                                        - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average)
+      IF( srcv(jpr_cal)%laction )    CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving
+      IF( srcv(jpr_icb)%laction )    CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs
+      IF( iom_use('snowpre') )       CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
+      IF( iom_use('precip') )        CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
+      IF( iom_use('rain') )          CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation 
+      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
+      IF( iom_use('rain_ao_cea') )   CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average)
+      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )     ! Sublimation over sea-ice (cell average)
+      IF( iom_use('evap_ao_cea') )   CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
+         &                                                         - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
       ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
       !
@@ -1769,4 +1910,7 @@
       CASE( 'oce only' )         ! the required field is directly provided
          zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
+         ! For Met Office sea ice non-solar fluxes are already delt with by JULES so setting to zero
+         ! here so the only flux is the ocean only one.
+         zqns_ice(:,:,:) = 0._wp 
       CASE( 'conservative' )     ! the required fields are directly provided
          zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
@@ -1798,5 +1942,5 @@
                zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    &
                   &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
-                  &                                                                + pist(:,:,jl) * picefr(:,:) ) )
+                  &                                             + pist(:,:,jl) * picefr(:,:) ) )
             END DO
          ELSE
@@ -1804,5 +1948,5 @@
                zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    &
                   &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
-                  &                                                                + pist(:,:,jl) * picefr(:,:) ) )
+                  &                                             + pist(:,:,jl) * picefr(:,:) ) )
             END DO
          ENDIF
@@ -1910,4 +2054,7 @@
       CASE( 'oce only' )
          zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
+         ! For Met Office sea ice solar fluxes are already delt with by JULES so setting to zero
+         ! here so the only flux is the ocean only one.
+         zqsr_ice(:,:,:) = 0._wp
       CASE( 'conservative' )
          zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
@@ -1995,4 +2142,6 @@
             ENDDO
          ENDIF
+      CASE( 'none' ) 
+         zdqns_ice(:,:,:) = 0._wp
       END SELECT
       
@@ -2010,9 +2159,12 @@
       !                                                      ! ========================= !
       CASE ('coupled')
-         IF( ln_mixcpl ) THEN
-            DO jl=1,jpl
-               qml_ice(:,:,jl) = qml_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_topm)%z3(:,:,jl) * zmsk(:,:)
-               qcn_ice(:,:,jl) = qcn_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_botm)%z3(:,:,jl) * zmsk(:,:)
-            ENDDO
+         IF (ln_scale_ice_flux) THEN
+            WHERE( a_i(:,:,:) > 1.e-10_wp )
+               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+            ELSEWHERE
+               qml_ice(:,:,:) = 0.0_wp
+               qcn_ice(:,:,:) = 0.0_wp
+            END WHERE
          ELSE
             qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
@@ -2025,19 +2177,32 @@
       IF( .NOT.ln_cndflx ) THEN                              !==  No conduction flux as surface forcing  ==!
          !
-         !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission when hi>10cm (Grenfell Maykut 77)
-         !
-         WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
-            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( ztri + ( 1._wp - ztri ) * ( 1._wp - phi(:,:,:) * 10._wp ) )
-         ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (ztri) when hi>10cm
-            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ztri
-         ELSEWHERE                                                         ! zero when hs>0
-            zqtr_ice_top(:,:,:) = 0._wp
-         END WHERE
+         IF( nn_qtrice == 0 ) THEN
+            ! formulation derived from Grenfell and Maykut (1977), where transmission rate
+            !    1) depends on cloudiness
+            !       ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
+            !       !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
+            !    2) is 0 when there is any snow
+            !    3) tends to 1 for thin ice
+            ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm
+            DO jl = 1, jpl
+               WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
+               ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )       ! constant (ztri) when hi>10cm
+                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ztri(:,:)
+               ELSEWHERE                                                           ! zero when hs>0
+                  zqtr_ice_top(:,:,jl) = 0._wp 
+               END WHERE
+            ENDDO
+         ELSEIF( nn_qtrice == 1 ) THEN
+            ! formulation is derived from the thesis of M. Lebrun (2019).
+            !    It represents the best fit using several sets of observations
+            !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
+            zqtr_ice_top(:,:,:) = 0.3_wp * zqsr_ice(:,:,:)
+         ENDIF
          !     
       ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
          !
-         !                    ! ===> here we must receive the qtr_ice_top array from the coupler
-         !                           for now just assume zero (fully opaque ice)
+         !          ! ===> here we must receive the qtr_ice_top array from the coupler
+         !                 for now just assume zero (fully opaque ice)
          zqtr_ice_top(:,:,:) = 0._wp
          !
@@ -2096,4 +2261,5 @@
       !
       isec = ( kt - nit000 ) * NINT( rn_Dt )        ! date of exchanges
+      info = OASIS_idle
 
       zfr_l(:,:) = 1.- fr_i(:,:)
@@ -2234,4 +2400,16 @@
       ENDIF
 
+#if defined key_si3 || defined key_cice
+      ! If this coupling was successful then save ice fraction for use between coupling points. 
+      ! This is needed for some calculations where the ice fraction at the last coupling point 
+      ! is needed. 
+      IF(  info == OASIS_Sent    .OR. info == OASIS_ToRest .OR. & 
+         & info == OASIS_SentOut .OR. info == OASIS_ToRestOut ) THEN 
+         IF ( sn_snd_thick%clcat == 'yes' ) THEN 
+           a_i_last_couple(:,:,1:jpl) = a_i(:,:,1:jpl)
+         ENDIF
+      ENDIF
+#endif
+
       IF( ssnd(jps_fice1)%laction ) THEN
          SELECT CASE( sn_snd_thick1%clcat )
@@ -2297,5 +2475,5 @@
             SELECT CASE( sn_snd_mpnd%clcat )  
             CASE( 'yes' )  
-               ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl)
+               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl)
                ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)  
             CASE( 'no' )  
@@ -2303,6 +2481,6 @@
                ztmp4(:,:,:) = 0.0  
                DO jl=1,jpl  
-                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)  
-                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl) 
+                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)
+                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl)
                ENDDO  
             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )  
