Index: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icectl.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icectl.F90	(revision 14589)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icectl.F90	(revision 14590)
@@ -801,6 +801,6 @@
       !!   DO jj = 1, jpj
       !!      DO ji = 1, jpi
-      !!         zdiag_mass2D(ji,jj) =   wfx_ice(ji,jj)   + wfx_snw(ji,jj)   + wfx_spr(ji,jj) + wfx_sub(ji,jj) &
-      !!            &                  + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj)
+      !!         zdiag_mass2D(ji,jj) =   wfx_ice(ji,jj)   + wfx_snw(ji,jj)   + wfx_spr(ji,jj)   + wfx_sub(ji,jj) + wfx_pnd(ji,jj) &
+      !!            &                  + diag_vice(ji,jj) + diag_vsnw(ji,jj) + diag_vpnd(ji,jj) - diag_adv_mass(ji,jj)
       !!         zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj)
       !!         zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj)
@@ -815,6 +815,6 @@
 
       ! -- mass diag -- !
-      zdiag_mass     = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub &
-         &                                  + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice
+      zdiag_mass     = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub + wfx_pnd &
+         &                                  + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) * rdt_ice
       zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rdt_ice
 
Index: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icesbc.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icesbc.F90	(revision 14589)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icesbc.F90	(revision 14590)
@@ -140,5 +140,5 @@
       CASE( jp_blk )              !--- bulk formulation
                                   CALL blk_ice_flx    ( t_su, h_s, h_i, alb_ice )    ! 
-         IF( ln_mixcpl        )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
+         IF( ln_mixcpl        )   CALL sbc_cpl_ice_flx( kt, picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
          IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
          !                        !    compute conduction flux and surface temperature (as in Jules surface module)
@@ -146,5 +146,5 @@
             &                     CALL blk_ice_qcn    ( ln_virtual_itd, t_su, t_bo, h_s, h_i )
       CASE ( jp_purecpl )         !--- coupled formulation
-                                  CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
+                                  CALL sbc_cpl_ice_flx( kt, picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
          IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
       END SELECT
Index: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceupdate.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceupdate.F90	(revision 14589)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceupdate.F90	(revision 14590)
@@ -105,5 +105,10 @@
       ! Net heat flux on top of the ice-ocean (W.m-2)
       !----------------------------------------------
-      qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:) 
+      IF( ln_cndflx ) THEN   ! ice-atm interface = conduction (and melting) fluxes
+         qt_atm_oi(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:) + &
+            &             SUM( a_i_b * ( qcn_ice + qml_ice + qtr_ice_top ), dim=3 ) + qemp_ice(:,:)
+      ELSE                   ! ice-atm interface = solar and non-solar fluxes
+         qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:) 
+      ENDIF
 
       ! --- case we bypass ice thermodynamics --- !
@@ -121,13 +126,19 @@
             ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2) 
             !---------------------------------------------------
-            zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) )
-
+            IF( ln_cndflx ) THEN   ! ice-atm interface = conduction (and melting) fluxes
+               zqsr = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) + SUM( a_i_b (ji,jj,:) * qtr_ice_bot(ji,jj,:) )
+            ELSE                   ! ice-atm interface = solar and non-solar fluxes
+               zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) )
+            ENDIF
+         
             ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2) 
             !---------------------------------------------------
-            qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) &
-               &                                - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) &
-               &                                + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) &
-               &                                + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj)                 
-            
+            IF( ln_icethd ) THEN
+               qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) &
+                  &                                - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) &
+                  &                                + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) &
+                  &                                + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj)
+            ENDIF
+         
             ! New qsr and qns used to compute the oceanic heat flux at the next time step
             !----------------------------------------------------------------------------
@@ -282,7 +293,7 @@
       IF( iom_use('hfxcndtop'  ) )   CALL iom_put( 'hfxcndtop'  , SUM( qcn_ice_top * a_i_b, dim=3 ) )   ! Surface conduction flux
 !!    IF( iom_use('hfxmelt'    ) )   CALL iom_put( 'hfxmelt'    , SUM( qml_ice     * a_i_b, dim=3 ) )   ! Surface melt flux
-!!    IF( iom_use('hfxldmelt'  ) )   CALL iom_put( 'hfxldmelt'  ,      fhld        * at_i_b         )   ! Heat in lead for ice melting
+!!    IF( iom_use('hfxldmelt'  ) )   CALL iom_put( 'hfxldmelt'  ,      fhld        * at_i_b         )   ! Heat in lead for ice melting 
 !!    IF( iom_use('hfxldgrow'  ) )   CALL iom_put( 'hfxldgrow'  ,      qlead       * r1_rdtice      )   ! Heat in lead for ice growth
-      
+
       ! controls
       !---------
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbccpl.F90	(revision 14589)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbccpl.F90	(revision 14590)
@@ -120,4 +120,5 @@
    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 ::   jpr_qtrice = 58   ! Transmitted solar thru sea-ice
 
    INTEGER, PARAMETER ::   jprcv      = 57   ! total number of fields received  
@@ -193,4 +194,6 @@
                     sn_rcv_wdrag, sn_rcv_wfreq
    !                                   ! Other namelist parameters
+!!   TYPE(FLD_C) ::   sn_rcv_qtrice
+!!   !                                   ! 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
@@ -227,5 +230,5 @@
       !!             ***  FUNCTION sbc_cpl_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(5)
+      INTEGER :: ierr(4)
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -237,9 +240,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 )
@@ -277,5 +277,5 @@
          &                  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
+         &                  sn_rcv_ts_ice !!, sn_rcv_qtrice
       !!---------------------------------------------------------------------
       !
@@ -319,4 +319,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,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')' 
@@ -575,7 +576,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'
@@ -588,4 +589,17 @@
          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   !   
@@ -844,8 +858,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        ! 
@@ -1516,6 +1526,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
@@ -1524,4 +1534,6 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty 
       !!----------------------------------------------------------------------
+      !
+#if defined key_si3 || defined key_cice
       !
       IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
@@ -1600,8 +1612,10 @@
       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  ***
@@ -1645,4 +1659,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
@@ -1653,4 +1674,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
@@ -1669,4 +1691,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)
@@ -1696,6 +1727,4 @@
          CALL ctl_stop('STOP', 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl')
       END SELECT
-
-#if defined key_si3
 
       ! --- evaporation over ice (kg/m2/s) --- !
@@ -1789,38 +1818,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(:,:) = - 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(:,:) = - 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
@@ -1833,15 +1858,39 @@
       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)  &
-         &                                                         - zevap_ice_total(:,:) * 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 conductive flux  !
+      !                                                      ! ================================= !
+      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
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
       !                                                      ! ========================= !
       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)
@@ -1890,5 +1939,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
@@ -1943,26 +1991,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
@@ -1985,10 +2033,34 @@
       !
       !                                                      ! ========================= !
+      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)
       !                                                      ! ========================= !
       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' )
@@ -2043,68 +2115,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]
@@ -2138,17 +2146,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)
-         zqtr_ice_top(:,:,:) = 0._wp
+!!         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
+!!         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.   !
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcice_cice.F90	(revision 14589)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcice_cice.F90	(revision 14590)
@@ -133,5 +133,5 @@
             CALL cice_sbc_force(kt)
          ELSE IF ( ksbc == jp_purecpl ) THEN
-            CALL sbc_cpl_ice_flx( fr_i )
+            CALL sbc_cpl_ice_flx( kt, fr_i )
          ENDIF
 
