Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4733)
@@ -230,5 +230,4 @@
    ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core)
    ln_blk_mfs  = .false.   !  MFS bulk formulation                      (T => fill namsbc_mfs )
-   ln_cpl      = .false.   !  Coupled formulation                       (T => fill namsbc_cpl )
    ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
    nn_ice      = 2         !  =0 no ice boundary condition   ,
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4733)
@@ -29,5 +29,4 @@
    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE wrk_nemo         ! work arrays
-   USE cpl_oasis3, ONLY : lk_cpl
 
    IMPLICIT NONE
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4733)
@@ -98,5 +98,5 @@
       !!              - fr_i    : ice fraction
       !!              - tn_ice  : sea-ice surface temperature
-      !!              - alb_ice : sea-ice alberdo (lk_cpl=T)
+      !!              - alb_ice : sea-ice albedo (lk_cpl=T)
       !!
       !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
@@ -136,8 +136,7 @@
             ! Solar heat flux reaching the ocean = zfcm1 (W.m-2) 
             !---------------------------------------------------
-            IF( lk_cpl ) THEN ! be careful: not been tested yet
-               ! original line
+            IF( lk_cpl ) THEN 
+               !!! LIM2 version zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
                zfcm1 = qsr_tot(ji,jj)
-               !!! LIM2 version zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
                DO jl = 1, jpl
                   zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * old_a_i(ji,jj,jl)
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4733)
@@ -68,5 +68,5 @@
       !!                ***  ROUTINE lim_thd  ***       
       !!  
-      !! ** Purpose : This routine manages the ice thermodynamic.
+      !! ** Purpose : This routine manages ice thermodynamics
       !!         
       !! ** Action : - Initialisation of some variables
@@ -74,11 +74,12 @@
       !!               at the ice base, snow acc.,heat budget of the leads)
       !!             - selection of the icy points and put them in an array
-      !!             - call lim_vert_ther for vert ice thermodynamic
-      !!             - back to the geographic grid
-      !!             - selection of points for lateral accretion
-      !!             - call lim_lat_acc  for the ice accretion
+      !!             - call lim_thd_dif  for vertical heat diffusion
+      !!             - call lim_thd_dh   for vertical ice growth and melt
+      !!             - call lim_thd_ent  for enthalpy remapping
+      !!             - call lim_thd_sal  for ice desalination
+      !!             - call lim_thd_temp to  retrieve temperature from ice enthalpy
       !!             - back to the geographic grid
       !!     
-      !! ** References : H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
+      !! ** References : 
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt    ! number of iteration
@@ -93,5 +94,9 @@
       !
       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      !
+      REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns
       !!-------------------------------------------------------------------
+      CALL wrk_alloc( jpi, jpj, zqsr, zqns )
+
       IF( nn_timing == 1 )  CALL timing_start('limthd')
 
@@ -137,4 +142,24 @@
       !-----------------------------------------------------------------------------!
 
+      !--- Ocean solar and non solar fluxes to be used in zqld
+      IF ( .NOT. lk_cpl ) THEN   ! --- forced case, fluxes to the lead are the same as over the ocean
+         !
+         zqsr(:,:) = qsr(:,:)      ; zqns(:,:) = qns(:,:)
+         !
+      ELSE                       ! --- coupled case, fluxes to the lead are total - intercepted
+         !
+         zqsr(:,:) = qsr_tot(:,:)  ; zqns(:,:) = qns_tot(:,:)
+         !
+         DO jl = 1, jpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * old_a_i(ji,jj,jl)
+                  zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * old_a_i(ji,jj,jl)
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+
 !CDIR NOVERRCHK
       DO jj = 1, jpj
@@ -149,11 +174,13 @@
             !           !  temperature and turbulent mixing (McPhee, 1992)
             !
+
             ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
-            zqld =  tms(ji,jj) * rdt_ice *                                       &
-               &  ( pfrld(ji,jj)         * ( qsr(ji,jj) * oatte(ji,jj)           &   ! solar heat + clem modif
-               &                           + qns(ji,jj) )                        &   ! non solar heat
-               ! latent heat of precip (note that precip is included in qns but not in qns_ice)
-               &    + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  &
+            zqld =  tms(ji,jj) * rdt_ice *                                          &
+               &  ( pfrld(ji,jj)         * ( zqsr(ji,jj) * oatte(ji,jj) + zqns(ji,jj) ) 
+               &    + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) *     & ! heat content of precip
+               &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  &
                &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) )
+               ! REMARK valid at least in forced mode from clem
+               ! precip is included in qns but not in qns_ice
 
             !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
@@ -482,8 +509,13 @@
       ENDIF
       !
+      !
+      CALL wrk_dealloc( jpi, jpj, zqsr, zqns )
+
+      !
       ! conservation test
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
       IF( nn_timing == 1 )  CALL timing_stop('limthd')
+
    END SUBROUTINE lim_thd 
 
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4733)
@@ -166,9 +166,9 @@
       !
       DO ji = kideb, kiut
-         zinda         = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )
-         ztmelts       = zinda * rtt + ( 1._wp - zinda ) * rtt
-
-         zfdum     = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
-         zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
+         zinda      = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )
+         ztmelts    = zinda * rtt + ( 1._wp - zinda ) * rtt
+
+         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
+         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
 
          zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) )
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4733)
@@ -146,6 +146,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid   ! tridiagonal system terms
       ! diag errors on heat
-      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini
-      REAL(wp)                        :: zhfx_err
+      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err
       !!------------------------------------------------------------------     
       ! 
@@ -158,5 +157,5 @@
       CALL wrk_alloc( jpij, jkmax+2, 3, ztrid )
 
-      CALL wrk_alloc( jpij, zdq, zq_ini )
+      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err )
 
       ! --- diag error on heat diffusion - PART 1 --- !
@@ -407,11 +406,16 @@
          !------------------------------------------------------------------------------|
          !
-         DO ji = kideb , kiut
-            ! update of the non solar flux according to the update in T_su
-            qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) )
-
+         IF( .NOT. lk_cpl ) THEN   !--- forced atmosphere case
+            DO ji = kideb , kiut
+               ! update of the non solar flux according to the update in T_su
+               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) )
+            END DO
+         ENDIF
+
+         ! Update incoming flux
+         DO ji = kideb , kiut
             ! update incoming flux
             zf(ji)    =   zfsw(ji)              & ! net absorbed solar radiation
-               + qns_ice_1d(ji)                  ! non solar total flux 
+               + qns_ice_1d(ji)                   ! non solar total flux 
             ! (LWup, LWdw, SH, LH)
          END DO
@@ -737,14 +741,30 @@
       CALL lim_thd_enmelt( kideb, kiut )
 
-      ! --- diag error on heat diffusion - PART 2 --- !
+      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
       DO ji = kideb, kiut
          zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  &
             &                              SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )
-         zhfx_err    = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 
-         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji)
-         ! --- correction of qns_ice and surface conduction flux --- !
-         qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err 
-         fc_su     (ji) = fc_su     (ji) - zhfx_err 
-         ! --- Heat flux at the ice surface in W.m-2 --- !
+         zhfx_err(ji)   = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 
+         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_b(ji)
+      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
+
+      ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m2)
+      DO ji = kideb, kiut
          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
          hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
@@ -759,5 +779,5 @@
       CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )
       CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid )
-      CALL wrk_dealloc( jpij, zdq, zq_ini )
+      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
 
    END SUBROUTINE lim_thd_dif
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 4733)
@@ -35,5 +35,4 @@
    LOGICAL , PUBLIC ::   ln_blk_core    !: CORE bulk formulation
    LOGICAL , PUBLIC ::   ln_blk_mfs     !: MFS  bulk formulation
-   LOGICAL , PUBLIC ::   ln_cpl         !: coupled   formulation (overwritten by key_sbc_coupled )
    LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr)
    LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4733)
@@ -518,5 +518,5 @@
       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pui      ! ice surface velocity (i- and i- components      [m/s]
       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid)
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb     ! ice albedo (clear sky) (alb_ice_cs)               [%]
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb     ! ice albedo (all skies)                            [%]
       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2]
       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid)
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 4733)
@@ -1125,6 +1125,6 @@
       REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
       ! optional arguments, used only in 'mixed oce-ice' case
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo 
-      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature     [Celcius]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! all skies ice albedo 
+      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature     [Celsius]
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature     [Kelvin]
       !
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4733)
@@ -82,5 +82,5 @@
       INTEGER ::   icpt   ! local integer
       !!
-      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   &
+      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core,           &
          &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
          &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, nn_lsm, nn_limflx
@@ -123,5 +123,4 @@
          WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core = ', ln_blk_core
          WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
-         WRITE(numout,*) '              coupled    formulation (T if key_sbc_cpl)  ln_cpl      = ', ln_cpl
          WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx
          WRITE(numout,*) '           Misc. options of sbc : '
@@ -237,5 +236,5 @@
       IF( ln_blk_core     ) THEN   ;   nsbc =  4   ; icpt = icpt + 1   ;   ENDIF       ! CORE bulk       formulation
       IF( ln_blk_mfs      ) THEN   ;   nsbc =  6   ; icpt = icpt + 1   ;   ENDIF       ! MFS  bulk       formulation
-      IF( ln_cpl          ) THEN   ;   nsbc =  5   ; icpt = icpt + 1   ;   ENDIF       ! Coupled         formulation
+      IF( lk_cpl          ) THEN   ;   nsbc =  5   ; icpt = icpt + 1   ;   ENDIF       ! Coupled         formulation
       IF( cp_cfg == 'gyre') THEN   ;   nsbc =  0                       ;   ENDIF       ! GYRE analytical formulation
       IF( lk_esopa        )            nsbc = -1                                       ! esopa test, ALL formulations
Index: /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/SAS_SRC/sbcssm.F90
===================================================================
--- /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/SAS_SRC/sbcssm.F90	(revision 4732)
+++ /branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/SAS_SRC/sbcssm.F90	(revision 4733)
@@ -166,8 +166,4 @@
       !! note that we need sbc_ssm called first in sbc
       !
-      IF( ln_cpl ) THEN
-         IF( lwp ) WRITE(numout,*) 'Coupled mode not sensible with StandAlone Surface scheme'
-         ln_cpl = .FALSE.
-      ENDIF
       IF( ln_apr_dyn ) THEN
          IF( lwp ) WRITE(numout,*) 'No atmospheric gradient needed with StandAlone Surface scheme'
