Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/SBC/sbccpl.F90	(revision 14353)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/SBC/sbccpl.F90	(revision 15574)
@@ -129,4 +129,5 @@
    INTEGER, PARAMETER ::   jpr_icb    = 61
    INTEGER, PARAMETER ::   jpr_ts_ice = 62   ! Sea ice surface temp
+   !!INTEGER, PARAMETER ::   jpr_qtrice = 63   ! Transmitted solar thru sea-ice
 
    INTEGER, PARAMETER ::   jprcv      = 62   ! total number of fields received
@@ -202,4 +203,5 @@
       &             sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd
    !                                   ! Other namelist parameters
+!!   TYPE(FLD_C) ::   sn_rcv_qtrice
    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
@@ -237,5 +239,5 @@
       !!             ***  FUNCTION sbc_cpl_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(5)
+      INTEGER :: ierr(4)
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -247,9 +249,6 @@
 #endif
       ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
-#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) )
+      !
+      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 )
@@ -286,5 +285,6 @@
          &                  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   , sn_rcv_icb  , sn_rcv_isf   , sn_rcv_ts_ice
+         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_icb  , sn_rcv_isf   , sn_rcv_ts_ice,  & !!, sn_rcv_qtrice
+         &                  sn_rcv_mslp
 
       !!---------------------------------------------------------------------
@@ -327,4 +327,5 @@
          WRITE(numout,*)'      ice shelf                       = ', TRIM(sn_rcv_isf%cldes   ), ' (', TRIM(sn_rcv_isf%clcat   ), ')'
          WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
+!!       WRITE(numout,*)'      transmitted solar thru sea-ice  = ', TRIM(sn_rcv_qtrice%cldes), ' (', TRIM(sn_rcv_qtrice%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), ')'
@@ -528,6 +529,6 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '   iceshelf received from oasis '
-         CALL ctl_stop('STOP','not coded')
-      ENDIF
+      ENDIF
+      !
       !
       !                                                      ! ------------------------- !
@@ -602,7 +603,7 @@
       srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE.
       !
-      !                                                      ! ------------------------- !
-      !                                                      !  ice topmelt and botmelt  !
-      !                                                      ! ------------------------- !
+      !                                                      ! --------------------------------- !
+      !                                                      !  ice topmelt and conduction flux  !   
+      !                                                      ! --------------------------------- !
       srcv(jpr_topm )%clname = 'OTopMlt'
       srcv(jpr_botm )%clname = 'OBotMlt'
@@ -615,4 +616,16 @@
          srcv(jpr_topm:jpr_botm)%laction = .TRUE.
       ENDIF
+!!      !                                                      ! --------------------------- !
+!!      !                                                      ! transmitted solar thru ice  !   
+!!      !                                                      ! --------------------------- !
+!!      srcv(jpr_qtrice)%clname = 'OQtr'
+!!      IF( TRIM(sn_rcv_qtrice%cldes) == 'coupled' ) THEN
+!!         IF ( TRIM( sn_rcv_qtrice%clcat ) == 'yes' ) THEN
+!!            srcv(jpr_qtrice)%nct = nn_cats_cpl
+!!         ELSE
+!!           CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtrice%clcat should always be set to yes currently' )
+!!         ENDIF
+!!         srcv(jpr_qtrice)%laction = .TRUE.
+!!      ENDIF
       !                                                      ! ------------------------- !
       !                                                      !    ice skin temperature   !
@@ -888,8 +901,4 @@
       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        !
@@ -1248,5 +1257,5 @@
                   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.0_wp, frcv(jpr_oty1)%z3(:,:,1), 'V',  -1.0_wp )
+               CALL lbc_lnk( 'sbccpl', frcv(jpr_otx1)%z3(:,:,1), 'U',  -1.0_wp, frcv(jpr_oty1)%z3(:,:,1), 'V',  -1.0_wp )
             ENDIF
             llnewtx = .TRUE.
@@ -1293,5 +1302,5 @@
          IF( llnewtau ) THEN
             zcoef = 1. / ( zrhoa * zcdrag )
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
             END_2D
@@ -1512,5 +1521,5 @@
          ! ice shelf fwf
          IF( srcv(jpr_isf)%laction )  THEN
-            fwfisf_oasis(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)
+            fwfisf_oasis(:,:) = frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf to the ocean ( > 0 = melting )
          END IF
 
@@ -1589,6 +1598,6 @@
       !! ** Action  :   return ptau_i, ptau_j, the stress over the ice
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
-      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
+      REAL(wp), INTENT(inout), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
+      REAL(wp), INTENT(inout), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
       !!
       INTEGER ::   ji, jj   ! dummy loop indices
@@ -1597,4 +1606,6 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty
       !!----------------------------------------------------------------------
+      !
+#if defined key_si3 || defined key_cice
       !
       IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1
@@ -1666,13 +1677,15 @@
                p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
             END_2D
-            CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
+            CALL lbc_lnk( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
          END SELECT
 
       ENDIF
       !
+#endif
+      !
    END SUBROUTINE sbc_cpl_ice_tau
 
 
-   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi )
+   SUBROUTINE sbc_cpl_ice_flx( kt, picefr, palbi, psst, pist, phs, phi )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
@@ -1716,4 +1729,11 @@
       !!                                                                      are provided but not included in emp here. Only runoff will
       !!                                                                      be included in emp in other parts of NEMO code
+      !!
+      !! ** Note : In case of the ice-atm coupling with conduction fluxes (such as Jules interface for the Met-Office),
+      !!              qsr_ice and qns_ice are not provided and they are not supposed to be used in the ice code.
+      !!              However, by precaution we also "fake" qns_ice and qsr_ice this way:
+      !!              qns_ice = qml_ice + qcn_ice ??
+      !!              qsr_ice = qtr_ice_top ??
+      !!
       !! ** Action  :   update at each nf_ice time step:
       !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
@@ -1724,4 +1744,5 @@
       !!                   sprecip           solid precipitation over the ocean
       !!----------------------------------------------------------------------
+      INTEGER,  INTENT(in)                                ::   kt         ! ocean model time step index (only for a_i_last_couple)
       REAL(wp), INTENT(in)   , DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
       !                                                   !!           ! optional arguments, used only in 'mixed oce-ice' case or for Met-Office coupling
@@ -1740,4 +1761,13 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
       !!----------------------------------------------------------------------
+      !
+#if defined key_si3 || defined key_cice
+      !
+      IF( kt == nit000 ) THEN
+         ! allocate ice fractions from last coupling time here and not in sbc_cpl_init because of jpl
+         IF( .NOT.ALLOCATED(a_i_last_couple) )   ALLOCATE( a_i_last_couple(jpi,jpj,jpl) )
+         ! initialize to a_i for the 1st time step
+         a_i_last_couple(:,:,:) = a_i(:,:,:)
+      ENDIF
       !
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -1768,6 +1798,4 @@
       END SELECT
 
-#if defined key_si3
-
       ! --- evaporation over ice (kg/m2/s) --- !
       IF (ln_scale_ice_flux) THEN ! typically met-office requirements
@@ -1834,6 +1862,6 @@
          rnf(:,:)    = rnf(:,:) + fwficb(:,:)
       ENDIF
-      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
-        fwfisf_oasis(:,:) = - frcv(jpr_isf)%z3(:,:,1)
+      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf > 0 mean melting)
+        fwfisf_oasis(:,:) = frcv(jpr_isf)%z3(:,:,1)
       ENDIF
 
@@ -1860,38 +1888,34 @@
       ENDIF
 
-#else
-      zsnw(:,:) = picefr(:,:)
-      ! --- Continental fluxes --- !
-      IF( srcv(jpr_rnf)%laction ) THEN   ! runoffs (included in emp later on)
-         rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
-      ENDIF
-      IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot)
-         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
-      ENDIF
-      IF( srcv(jpr_icb)%laction ) THEN   ! iceberg added to runoffs
-         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
-         rnf(:,:)    = rnf(:,:) + fwficb(:,:)
-      ENDIF
-      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
-        fwfisf_oasis(:,:) = - frcv(jpr_isf)%z3(:,:,1)
-      ENDIF
-      !
-      IF( ln_mixcpl ) THEN
-         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
-         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
-         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
-         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
-      ELSE
-         emp_tot(:,:) =                                  zemp_tot(:,:)
-         emp_ice(:,:) =                                  zemp_ice(:,:)
-         sprecip(:,:) =                                  zsprecip(:,:)
-         tprecip(:,:) =                                  ztprecip(:,:)
-      ENDIF
-      !
-#endif
-
+!! for CICE ??
+!!$      zsnw(:,:) = picefr(:,:)
+!!$      ! --- Continental fluxes --- !
+!!$      IF( srcv(jpr_rnf)%laction ) THEN   ! runoffs (included in emp later on)
+!!$         rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+!!$      ENDIF
+!!$      IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot)
+!!$         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
+!!$      ENDIF
+!!$      IF( srcv(jpr_icb)%laction ) THEN   ! iceberg added to runoffs
+!!$         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
+!!$         rnf(:,:)    = rnf(:,:) + fwficb(:,:)
+!!$      ENDIF
+!!$      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf >0 mean melting)
+!!$        fwfisf_oasis(:,:) = frcv(jpr_isf)%z3(:,:,1)
+!!$      ENDIF
+!!$      !
+!!$      IF( ln_mixcpl ) THEN
+!!$         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
+!!$         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
+!!$         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
+!!$         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
+!!$      ELSE
+!!$         emp_tot(:,:) =                                  zemp_tot(:,:)
+!!$         emp_ice(:,:) =                                  zemp_ice(:,:)
+!!$         sprecip(:,:) =                                  zsprecip(:,:)
+!!$         tprecip(:,:) =                                  ztprecip(:,:)
+!!$      ENDIF
+      !
       ! outputs
-!!      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
@@ -1901,9 +1925,29 @@
       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('rain_ao_cea') )   CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * ziceld(:,:)         )  ! liquid precipitation over ocean (cell average)
+      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , zevap_ice_total(:,:) * 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)
+         &                                                         - zevap_ice_total(:,:) * 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
+!!      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
+      !
+      !                                                      ! ========================= !
+      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
+      !                                                      ! ========================= !
+      CASE ('coupled')
+         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(:,:,:)
+            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
+         ENDIF
+      END SELECT
       !
       !                                                      ! ========================= !
@@ -1911,8 +1955,13 @@
       !                                                      ! ========================= !
       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
+         ! Get the sea ice non solar heat flux from conductive, melting and sublimation fluxes
+         IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
+            zqns_ice(:,:,:) = qml_ice(:,:,:) + qcn_ice(:,:,:)
+         ELSE
+            zqns_ice(:,:,:) = 0._wp
+         ENDIF
+         ! Calculate the total non solar heat flux. The ocean only non solar heat flux (zqns_oce) will be recalculated after this CASE
+         ! statement to be consistent with other coupling methods even though .zqns_oce = frcv(jpr_qnsoce)%z3(:,:,1)
+         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1) + SUM( zqns_ice(:,:,:) * a_i(:,:,:), dim=3 )
       CASE( 'conservative' )     ! the required fields are directly provided
          zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
@@ -1961,5 +2010,4 @@
       IF( srcv(jpr_icb)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus  ! remove latent heat of iceberg melting
 
-#if defined key_si3
       ! --- non solar flux over ocean --- !
       !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
@@ -2014,26 +2062,26 @@
       ENDIF
 
-#else
-      zcptsnw (:,:) = zcptn(:,:)
-      zcptrain(:,:) = zcptn(:,:)
-
-      ! clem: this formulation is certainly wrong... but better than it was...
-      zqns_tot(:,:) = zqns_tot(:,:)                             &          ! zqns_tot update over free ocean with:
-         &          - (  ziceld(:,:) * zsprecip(:,:) * rLfus )  &          ! remove the latent heat flux of solid precip. melting
-         &          - (  zemp_tot(:,:)                          &          ! remove the heat content of mass flux (assumed to be at SST)
-         &             - zemp_ice(:,:) ) * zcptn(:,:)
-
-     IF( ln_mixcpl ) THEN
-         qns_tot(:,:) = qns(:,:) * ziceld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
-         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
-         DO jl=1,jpl
-            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
-         ENDDO
-      ELSE
-         qns_tot(:,:  ) = zqns_tot(:,:  )
-         qns_ice(:,:,:) = zqns_ice(:,:,:)
-      ENDIF
-
-#endif
+!! for CICE ??
+!!$      ! --- non solar flux over ocean --- !
+!!$      zcptsnw (:,:) = zcptn(:,:)
+!!$      zcptrain(:,:) = zcptn(:,:)
+!!$
+!!$      ! clem: this formulation is certainly wrong... but better than it was...
+!!$      zqns_tot(:,:) = zqns_tot(:,:)                             &          ! zqns_tot update over free ocean with:
+!!$         &          - (  ziceld(:,:) * zsprecip(:,:) * rLfus )  &          ! remove the latent heat flux of solid precip. melting
+!!$         &          - (  zemp_tot(:,:)                          &          ! remove the heat content of mass flux (assumed to be at SST)
+!!$         &             - zemp_ice(:,:) ) * zcptn(:,:)
+!!$
+!!$     IF( ln_mixcpl ) THEN
+!!$         qns_tot(:,:) = qns(:,:) * ziceld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
+!!$         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
+!!$         DO jl=1,jpl
+!!$            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
+!!$         ENDDO
+!!$      ELSE
+!!$         qns_tot(:,:  ) = zqns_tot(:,:  )
+!!$         qns_ice(:,:,:) = zqns_ice(:,:,:)
+!!$      ENDIF
+
       ! outputs
       IF ( srcv(jpr_cal)%laction ) CALL iom_put('hflx_cal_cea' , - frcv(jpr_cal)%z3(:,:,1) * rLfus ) ! latent heat from calving
@@ -2053,6 +2101,31 @@
       IF (        iom_use('hflx_snow_ai_cea') ) &                                                    ! heat flux from snow (over ice)
          &   CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) *  zsnw(:,:) )
+      IF(         iom_use('hflx_subl_cea') )    &                                                    ! heat flux from sublimation
+         &   CALL iom_put('hflx_subl_cea' ,   SUM( qevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) * tmask(:,:,1) )
       ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
       !
+      !                                                      ! ========================= !
+      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
+      !                                                      ! ========================= !
+      CASE ('coupled')
+         IF( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
+            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
+         ELSE
+            ! Set all category values equal for the moment
+            DO jl=1,jpl
+               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
+            ENDDO
+         ENDIF
+      CASE( 'none' )
+         zdqns_ice(:,:,:) = 0._wp
+      END SELECT
+
+      IF( ln_mixcpl ) THEN
+         DO jl=1,jpl
+            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
+         ENDDO
+      ELSE
+         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
+      ENDIF
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
@@ -2060,6 +2133,6 @@
       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.
+         ! For the Met Office the only sea ice solar flux is the transmitted qsr which is added onto zqsr_ice
+         ! further down. Therefore start zqsr_ice off at zero.
          zqsr_ice(:,:,:) = 0._wp
       CASE( 'conservative' )
@@ -2114,68 +2187,4 @@
          END DO
       ENDIF
-
-#if defined key_si3
-      ! --- solar flux over ocean --- !
-      !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
-      zqsr_oce = 0._wp
-      WHERE( ziceld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / ziceld(:,:)
-
-      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
-      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
-#endif
-
-      IF( ln_mixcpl ) THEN
-         qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
-         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
-         DO jl = 1, jpl
-            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
-         END DO
-      ELSE
-         qsr_tot(:,:  ) = zqsr_tot(:,:  )
-         qsr_ice(:,:,:) = zqsr_ice(:,:,:)
-      ENDIF
-
-      !                                                      ! ========================= !
-      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
-      !                                                      ! ========================= !
-      CASE ('coupled')
-         IF( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
-            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
-         ELSE
-            ! Set all category values equal for the moment
-            DO jl=1,jpl
-               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
-            ENDDO
-         ENDIF
-      CASE( 'none' )
-         zdqns_ice(:,:,:) = 0._wp
-      END SELECT
-
-      IF( ln_mixcpl ) THEN
-         DO jl=1,jpl
-            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
-         ENDDO
-      ELSE
-         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
-      ENDIF
-
-#if defined key_si3
-      !                                                      ! ========================= !
-      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
-      !                                                      ! ========================= !
-      CASE ('coupled')
-         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(:,:,:)
-            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
-         ENDIF
-      END SELECT
       !                                                      ! ========================= !
       !                                                      !      Transmitted Qsr      !   [W/m2]
@@ -2209,17 +2218,50 @@
       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)
+!!         SELECT CASE( TRIM( sn_rcv_qtrice%cldes ) )
+!!            !
+!!            !      ! ===> here we receive the qtr_ice_top array from the coupler
+!!         CASE ('coupled')
+!!            IF (ln_scale_ice_flux) THEN
+!!               WHERE( a_i(:,:,:) > 1.e-10_wp )
+!!                  zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+!!               ELSEWHERE
+!!                  zqtr_ice_top(:,:,:) = 0.0_wp
+!!               ENDWHERE
+!!            ELSE
+!!               zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:)
+!!            ENDIF
+!!           
+!!            ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
+!!            zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
+!!            zqsr_tot(:,:)   = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
+!!            
+!!            !      if we are not getting this data from the coupler then assume zero (fully opaque ice)
+!!         CASE ('none')
          zqtr_ice_top(:,:,:) = 0._wp
-         !
-      ENDIF
-      !
+!!         END SELECT
+            !
+      ENDIF
+
       IF( ln_mixcpl ) THEN
-         DO jl=1,jpl
+         qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
+         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) + zqsr_tot(:,:) * zmsk(:,:)
+         DO jl = 1, jpl
+            qsr_ice    (:,:,jl) = qsr_ice    (:,:,jl) * xcplmask(:,:,0) + zqsr_ice    (:,:,jl) * zmsk(:,:)
             qtr_ice_top(:,:,jl) = qtr_ice_top(:,:,jl) * xcplmask(:,:,0) + zqtr_ice_top(:,:,jl) * zmsk(:,:)
-         ENDDO
+         END DO
       ELSE
+         qsr_tot    (:,:  ) = zqsr_tot    (:,:  )
+         qsr_ice    (:,:,:) = zqsr_ice    (:,:,:)
          qtr_ice_top(:,:,:) = zqtr_ice_top(:,:,:)
       ENDIF
+      
+      ! --- solar flux over ocean --- !
+      !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
+      zqsr_oce = 0._wp
+      WHERE( ziceld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / ziceld(:,:)
+
+      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
+      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
+
       !                                                      ! ================== !
       !                                                      !   ice skin temp.   !
@@ -2560,5 +2602,5 @@
                   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.0_wp, zity1, 'T', -1.0_wp )
+               CALL lbc_lnk( 'sbccpl', zitx1, 'T', -1.0_wp, zity1, 'T', -1.0_wp )
             CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T
                DO_2D( 0, 0, 0, 0 )
@@ -2569,5 +2611,5 @@
                END_2D
             END SELECT
-            CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp,  zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
+            CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp,  zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
             !
          ENDIF
@@ -2637,5 +2679,5 @@
                 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.0_wp,  zity1, 'T', -1.0_wp )
+             CALL lbc_lnk( 'sbccpl', zitx1, 'T', -1.0_wp,  zity1, 'T', -1.0_wp )
           CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T
              DO_2D( 0, 0, 0, 0 )
@@ -2646,5 +2688,5 @@
              END_2D
           END SELECT
-         CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocxw)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocyw)%clgrid, -1.0_wp )
+         CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocxw)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocyw)%clgrid, -1.0_wp )
          !
          !
