Index: NEMO/branches/2020/dev_12905_xios_ancil/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2020/dev_12905_xios_ancil/src/OCE/SBC/sbccpl.F90	(revision 13003)
+++ NEMO/branches/2020/dev_12905_xios_ancil/src/OCE/SBC/sbccpl.F90	(revision 13766)
@@ -41,5 +41,5 @@
 #endif
 #if defined key_si3
-   USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE icevar         ! for CALL ice_var_snwblow
 #endif
    !
@@ -48,4 +48,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
@@ -152,5 +156,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 +163,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 **
@@ -184,4 +196,6 @@
    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 +205,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] 
@@ -199,4 +216,5 @@
    !! Substitution
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -210,5 +228,5 @@
       !!             ***  FUNCTION sbc_cpl_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(4)
+      INTEGER :: ierr(5)
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -220,6 +238,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 )
@@ -248,14 +269,14 @@
       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_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_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_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal   ,  &
+         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                &
+         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 &
          &                  sn_rcv_ts_ice
-
       !!---------------------------------------------------------------------
       !
@@ -277,4 +298,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  ), ')'
@@ -325,7 +350,4 @@
          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
 
@@ -364,5 +386,7 @@
       ! 
       ! 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%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.
       
@@ -695,5 +719,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))
@@ -819,4 +843,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        ! 
@@ -1036,6 +1064,6 @@
          xcplmask(:,:,:) = 0.
          CALL iom_open( 'cplmask', inum )
-         CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:nlci,1:nlcj,1:nn_cplmodel),   &
-            &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nn_cplmodel /) )
+         CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:jpi,1:jpj,1:nn_cplmodel),   &
+            &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ jpi,jpj,nn_cplmodel /) )
          CALL iom_close( inum )
       ELSE
@@ -1107,5 +1135,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
       !!----------------------------------------------------------------------
       !
@@ -1167,9 +1195,9 @@
             !                              
             IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
-               DO_2D_00_00
+               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) )
                END_2D
-               CALL lbc_lnk_multi( 'sbccpl', frcv(jpr_otx1)%z3(:,:,1), 'U',  -1., frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
+               CALL lbc_lnk_multi( 'sbccpl', frcv(jpr_otx1)%z3(:,:,1), 'U',  -1.0_wp, frcv(jpr_oty1)%z3(:,:,1), 'V',  -1.0_wp )
             ENDIF
             llnewtx = .TRUE.
@@ -1191,10 +1219,10 @@
          ! => need to be done only when otx1 was changed
          IF( llnewtx ) THEN
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
                zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
                frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
             END_2D
-            CALL lbc_lnk( 'sbccpl', frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
+            CALL lbc_lnk( 'sbccpl', frcv(jpr_taum)%z3(:,:,1), 'T', 1.0_wp )
             llnewtau = .TRUE.
          ELSE
@@ -1216,10 +1244,22 @@
          IF( llnewtau ) THEN 
             zcoef = 1. / ( zrhoa * zcdrag ) 
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
             END_2D
          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      
@@ -1481,4 +1521,5 @@
       INTEGER ::   ji, jj   ! dummy loop indices
       INTEGER ::   itx      ! index of taux over ice
+      REAL(wp)                     ::   zztmp1, zztmp2
       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty 
       !!----------------------------------------------------------------------
@@ -1544,23 +1585,14 @@
             p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
             p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
-         CASE( 'F' )
-            DO_2D_00_00
-               p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
-               p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
+         CASE( 'T' )
+            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) )
+               zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) )
+               p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
+               p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
             END_2D
-         CASE( 'T' )
-            DO_2D_00_00
-               p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
-               p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
-            END_2D
-         CASE( 'I' )
-            DO_2D_00_00
-               p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
-               p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
-            END_2D
+            CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
          END SELECT
-         IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN 
-            CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
-         ENDIF
          
       ENDIF
@@ -1628,9 +1660,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
       !!----------------------------------------------------------------------
       !
@@ -1652,5 +1685,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)
@@ -1664,6 +1696,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) --- !
@@ -1672,11 +1743,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
@@ -1756,15 +1821,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
       !
@@ -1774,4 +1839,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)
@@ -1791,7 +1859,6 @@
             ENDDO
          ELSE
-            qns_tot(:,:) = qns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
+            zqns_tot(:,:) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
             DO jl = 1, jpl
-               zqns_tot(:,:   ) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
                zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
             END DO
@@ -1804,5 +1871,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
@@ -1810,5 +1877,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
@@ -1916,4 +1983,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)
@@ -1934,7 +2004,6 @@
             END DO
          ELSE
-            qsr_tot(:,:   ) = qsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
+            zqsr_tot(:,:) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
             DO jl = 1, jpl
-               zqsr_tot(:,:   ) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
             END DO
@@ -2002,4 +2071,6 @@
             ENDDO
          ENDIF
+      CASE( 'none' ) 
+         zdqns_ice(:,:,:) = 0._wp
       END SELECT
       
@@ -2017,9 +2088,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(:,:,:)
@@ -2032,19 +2106,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
          !
@@ -2103,4 +2190,5 @@
       !
       isec = ( kt - nit000 ) * NINT( rn_Dt )        ! date of exchanges
+      info = OASIS_idle
 
       zfr_l(:,:) = 1.- fr_i(:,:)
@@ -2241,4 +2329,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 )
@@ -2304,5 +2404,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' )  
@@ -2310,6 +2410,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' )  
@@ -2372,10 +2472,10 @@
             SELECT CASE( TRIM( sn_snd_crt%cldes ) )
             CASE( 'oce only'             )      ! C-grid ==> T
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) )
                   zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji  ,jj-1,1,Kmm) ) 
                END_2D
             CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T  
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)  
                   zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)
@@ -2383,7 +2483,7 @@
                   zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  )     + v_ice(ji  ,jj-1  )     ) *  fr_i(ji,jj)
                END_2D
-               CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1., zity1, 'T', -1. )
+               CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1.0_wp, zity1, 'T', -1.0_wp )
             CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T
-               DO_2D_00_00
+               DO_2D( 0, 0, 0, 0 )
                   zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   &
                      &         + 0.5 * ( u_ice(ji,jj  )     + u_ice(ji-1,jj    )     ) *  fr_i(ji,jj)
@@ -2392,5 +2492,5 @@
                END_2D
             END SELECT
-            CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.,  zoty1, ssnd(jps_ocy1)%clgrid, -1. )
+            CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp,  zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
             !
          ENDIF
@@ -2449,10 +2549,10 @@
           SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
           CASE( 'oce only'             )      ! C-grid ==> T 
-             DO_2D_00_00
+             DO_2D( 0, 0, 0, 0 )
                 zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) ) 
                 zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji , jj-1,1,Kmm) )  
              END_2D
           CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T   
-             DO_2D_00_00
+             DO_2D( 0, 0, 0, 0 )
                 zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   
                 zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj) 
@@ -2460,7 +2560,7 @@
                 zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
              END_2D
-             CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1.,  zity1, 'T', -1. ) 
+             CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1.0_wp,  zity1, 'T', -1.0_wp ) 
           CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T  
-             DO_2D_00_00
+             DO_2D( 0, 0, 0, 0 )
                 zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   & 
                    &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
@@ -2469,5 +2569,5 @@
              END_2D
           END SELECT
-         CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocxw)%clgrid, -1., zoty1, ssnd(jps_ocyw)%clgrid, -1. ) 
+         CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocxw)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocyw)%clgrid, -1.0_wp ) 
          ! 
          ! 
