Index: /trunk/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /trunk/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 5145)
+++ /trunk/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 5146)
@@ -94,4 +94,5 @@
                                    !     3: activate G(he) only              --- temporary option
                                    !     4: activate lateral melting only    --- temporary option
+  ln_it_qnsice = .true.            !  iterate the surface non-solar flux with surface temperature (T) or not (F)
 /
 !------------------------------------------------------------------------------
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5145)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5146)
@@ -198,4 +198,5 @@
    INTEGER , PUBLIC ::   nn_ice_thcon        !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007)
    INTEGER , PUBLIC ::   nn_monocat          !: virtual ITD mono-category parameterizations (1) or not (0)
+   LOGICAL , PUBLIC ::   ln_it_qnsice        !: iterate surface flux with changing surface temperature or not (F)
 
    !                                     !!** ice-mechanical redistribution namelist (namiceitdme)
@@ -285,4 +286,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations 
@@ -433,5 +435,6 @@
          &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,                        &
          &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,    &
-         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) , hfx_err_rem(jpi,jpj), &
+         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,     & 
+         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) ,                                   &
          &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,                           &
          &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) ,    &
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5145)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5146)
@@ -146,4 +146,8 @@
             hfx_out(ji,jj) = hfx_out(ji,jj) + zf_mass + zfcm1
 
+            ! Add the residual from heat diffusion equation (W.m-2)
+            !-------------------------------------------------------
+            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj)
+
             ! New qsr and qns used to compute the oceanic heat flux at the next time step
             !---------------------------------------------------
@@ -164,6 +168,7 @@
             !  computing freshwater exchanges at the ice/ocean interface
             IF( lk_cpl ) THEN 
-               zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !
-                  &   + wfx_snw(ji,jj)
+                zemp =   emp_tot(ji,jj)                                    &   ! net mass flux over grid cell
+                   &   - emp_ice(ji,jj) * ( 1._wp - pfrld(ji,jj) )         &   ! minus the mass flux intercepted by sea ice
+                   &   + sprecip(ji,jj) * ( pfrld(ji,jj) - pfrld(ji,jj)**rn_betas )   !
             ELSE
                zemp =   emp(ji,jj)     *           pfrld(ji,jj)            &   ! evaporation over oceanic fraction
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5145)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5146)
@@ -601,4 +601,5 @@
          CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err         , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res         , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, hfx_err_dif_1d (1:nbpb), hfx_err_dif , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
 
@@ -651,4 +652,5 @@
          CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj )
          CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj )
+         CALL tab_1d_2d( nbpb, hfx_err_dif   , npb, hfx_err_dif_1d(1:nbpb), jpi, jpj )
          !
          CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
@@ -675,5 +677,5 @@
       NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       &
          &                rn_himin, parsub, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &
-         &                nn_monocat
+         &                nn_monocat, ln_it_qnsice
       !!-------------------------------------------------------------------
       !
@@ -718,4 +720,5 @@
          WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i
          WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat
+         WRITE(numout,*)'      iterate the surface non-solar flux (T) or not (F)       ln_it_qnsice = ', ln_it_qnsice
       ENDIF
       !
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5145)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5146)
@@ -120,4 +120,5 @@
       REAL(wp), POINTER, DIMENSION(:)     ::   zh_s        ! snow layer thickness
       REAL(wp), POINTER, DIMENSION(:)     ::   zfsw        ! solar radiation absorbed at the surface
+      REAL(wp), POINTER, DIMENSION(:)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
       REAL(wp), POINTER, DIMENSION(:)     ::   zf          ! surface flux function
       REAL(wp), POINTER, DIMENSION(:)     ::   dzf         ! derivative of the surface flux function
@@ -168,5 +169,5 @@
       CALL wrk_alloc( jpij, numeqmin, numeqmax )
       CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
-      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe )
+      CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe )
       CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)
       CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0)
@@ -242,7 +243,8 @@
       !-------------------------------------------------------
       DO ji = kideb , kiut
-         zfsw   (ji) =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
-         zftrice(ji) =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
-         dzf    (ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
+         zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
+         zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
+         dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
+         zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value
       END DO
 
@@ -452,5 +454,5 @@
          !------------------------------------------------------------------------------|
          !
-         IF( .NOT. lk_cpl ) THEN   !--- forced atmosphere case
+         IF ( ln_it_qnsice ) THEN 
             DO ji = kideb , kiut
                ! update of the non solar flux according to the update in T_su
@@ -757,4 +759,6 @@
       CALL lim_thd_enmelt( kideb, kiut )
 
+      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
+      IF ( ln_it_qnsice ) hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - ( qns_ice_1d(:)  - zqns_ice_b(:) ) * a_i_1d(:) 
 
       ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
@@ -770,19 +774,5 @@
       END DO 
 
-      ! diagnose external surface (forced case) or bottom (forced case) from heat conservation
-      IF( .NOT. lk_cpl ) THEN   ! --- forced case: qns_ice and fc_su are diagnosed
-         !
-         DO ji = kideb, kiut
-            qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err(ji)
-            fc_su     (ji) = fc_su(ji)      - zhfx_err(ji)
-         END DO
-         !
-      ELSE                      ! --- coupled case: ocean turbulent heat flux is diagnosed
-         !
-         DO ji = kideb, kiut
-            fhtur_1d  (ji) = fhtur_1d(ji)   - zhfx_err(ji)
-         END DO
-         !
-      ENDIF
+      hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - zhfx_err(:) * a_i_1d(:)
 
       !-----------------------------------------
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5145)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5146)
@@ -55,4 +55,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_dif_1d
 
    ! heat flux associated with ice-atmosphere mass exchange
@@ -144,10 +145,10 @@
          &      qsr_ice_1d (jpij) ,     &
          &      fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qns_ice_1d(jpij) ,     &
-         &      t_bo_1d   (jpij) ,                                          &
-         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,     & 
-         &      hfx_dif_1d(jpij) ,hfx_opw_1d(jpij) , &
-         &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , &
-         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) , &
-         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij),       STAT=ierr(1) )
+         &      t_bo_1d   (jpij) ,                                         &
+         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,    & 
+         &      hfx_dif_1d(jpij) , hfx_opw_1d(jpij) ,                      &
+         &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) ,                      &
+         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) ,   &
+         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , STAT=ierr(1) )
       !
       ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_1d     (jpij) ,     &
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5145)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5146)
@@ -1199,6 +1199,6 @@
             ENDDO
          ELSE
+            qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
             DO jl=1,jpl
-               qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
                qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
             ENDDO
@@ -1258,6 +1258,6 @@
             ENDDO
          ELSE
+            qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
             DO jl=1,jpl
-               qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
                qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
             ENDDO
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5145)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5146)
@@ -602,4 +602,5 @@
       hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
       hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
+      hfx_err_dif(:,:) = 0._wp
 
       afx_tot(:,:) = 0._wp   ;
