Index: branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4766)
+++ branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5581)
@@ -25,13 +25,10 @@
    USE par_oce          ! ocean parameters
    USE phycst           ! physical constants
-   USE par_ice          ! ice parameters
    USE dom_oce          ! ocean domain
-   USE dom_ice,    ONLY : tms, area
    USE ice              ! LIM sea-ice variables
    USE sbc_ice          ! Surface boundary condition: sea-ice fields
    USE sbc_oce          ! Surface boundary condition: ocean fields
    USE sbccpl
-   USE cpl_oasis3, ONLY : lk_cpl
-   USE oce       , ONLY : iatte, oatte, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass
+   USE oce       , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass
    USE albedo           ! albedo parameters
    USE lbclnk           ! ocean lateral boundary condition - MPP exchanges
@@ -41,17 +38,16 @@
    USE prtctl           ! Print control
    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget
+   USE traqsr           ! add penetration of solar flux in the calculation of heat budget
    USE iom
    USE domvvl           ! Variable volume
+   USE limctl
+   USE limcons
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   lim_sbc_init   ! called by ice_init
+   PUBLIC   lim_sbc_init   ! called by sbc_lim_init
    PUBLIC   lim_sbc_flx    ! called by sbc_ice_lim
    PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim
-
-   REAL(wp)  ::   epsi10 = 1.e-10   !
-   REAL(wp)  ::   epsi20 = 1.e-20   !
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2]
@@ -98,28 +94,30 @@
       !!              - fr_i    : ice fraction
       !!              - tn_ice  : sea-ice surface temperature
-      !!              - alb_ice : sea-ice alberdo (lk_cpl=T)
+      !!              - alb_ice : sea-ice albedo (only useful in coupled mode)
       !!
       !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
       !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
+      !!              These refs are now obsolete since everything has been revised
+      !!              The ref should be Rousset et al., 2015
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! number of iteration
-      !
-      INTEGER  ::   ji, jj, jl, jk           ! dummy loop indices
-      REAL(wp) ::   zinda, zemp      ! local scalars
-      REAL(wp) ::   zf_mass         ! Heat flux associated with mass exchange ice->ocean (W.m-2)
-      REAL(wp) ::   zfcm1           ! New solar flux received by the ocean
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace
+      INTEGER, INTENT(in) ::   kt                                  ! number of iteration
+      INTEGER  ::   ji, jj, jl, jk                                 ! dummy loop indices
+      REAL(wp) ::   zqmass                                         ! Heat flux associated with mass exchange ice->ocean (W.m-2)
+      REAL(wp) ::   zqsr                                           ! New solar flux received by the ocean
+      !
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 2D/3D workspace
       !!---------------------------------------------------------------------
-      
-      IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp )
 
       ! make calls for heat fluxes before it is modified
-      CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) )   !     solar flux at ocean surface
-      CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) )   ! non-solar flux at ocean surface
-      CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux at ice surface
-      CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  ! non-solar flux at ice surface
-      CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux transmitted thru ice
-      CALL iom_put( "qt_oce"  , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) )  
-      CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * old_a_i(:,:,:), dim=3 ) )
+      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface
+      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface
+      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
+      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface
+      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
+      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )  
+      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
+         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
+      IF( iom_use('qemp_oce' ) )   CALL iom_put( "qemp_oce"  , qemp_oce(:,:) )  
+      IF( iom_use('qemp_ice' ) )   CALL iom_put( "qemp_ice"  , qemp_ice(:,:) )  
 
       ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
@@ -130,33 +128,24 @@
             !      heat flux at the ocean surface      !
             !------------------------------------------!
-            zinda   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( 1._wp - pfrld(ji,jj) ) ) ) ! 1 if ice
-
-            ! Solar heat flux reaching the ocean = zfcm1 (W.m-2) 
+            ! Solar heat flux reaching the ocean = zqsr (W.m-2) 
             !---------------------------------------------------
-            IF( lk_cpl ) THEN ! be carfeful: not been tested yet
-               ! original line
-               zfcm1 = qsr_tot(ji,jj)
-               !!!zfcm1 = qsr_tot(ji,jj) + ftr_ice(ji,jj) * ( 1._wp - pfrld(ji,jj) ) / ( 1._wp - zinda + zinda * iatte(ji,jj) )
-               DO jl = 1, jpl
-                  zfcm1 = zfcm1 - ( qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) * old_a_i(ji,jj,jl)
-               END DO
-            ELSE
-               !!!zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + &
-               !!!     &    ( 1._wp - pfrld(ji,jj) ) * ftr_ice(ji,jj) / ( 1._wp - zinda + zinda * iatte(ji,jj) )
-               zfcm1   = pfrld(ji,jj) * qsr(ji,jj)
-               DO jl = 1, jpl
-                  zfcm1   = zfcm1 + old_a_i(ji,jj,jl) * ftr_ice(ji,jj,jl)
-               END DO
-            ENDIF
+            zqsr = qsr_tot(ji,jj)
+            DO jl = 1, jpl
+               zqsr = zqsr - a_i_b(ji,jj,jl) * (  qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) 
+            END DO
 
             ! Total heat flux reaching the ocean = hfx_out (W.m-2) 
             !---------------------------------------------------
-            zf_mass        = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
-            hfx_out(ji,jj) = hfx_out(ji,jj) + zf_mass + zfcm1
+            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
+            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr
+
+            ! 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
             !---------------------------------------------------
-            qsr(ji,jj) = zfcm1                                      
-            qns(ji,jj) = hfx_out(ji,jj) - zfcm1              
+            qsr(ji,jj) = zqsr                                      
+            qns(ji,jj) = hfx_out(ji,jj) - zqsr              
 
             !------------------------------------------!
@@ -171,14 +160,4 @@
             !                     Even if i see Ice melting as a FW and SALT flux
             !        
-            !  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)
-            ELSE
-               zemp =   emp(ji,jj)     *           pfrld(ji,jj)            &   ! evaporation over oceanic fraction
-                  &   - tprecip(ji,jj) * ( 1._wp - pfrld(ji,jj) )          &   ! all precipitation reach the ocean
-                  &   + sprecip(ji,jj) * ( 1._wp - pfrld(ji,jj)**betas )       ! except solid precip intercepted by sea-ice
-            ENDIF
-
             ! mass flux from ice/ocean
             wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
@@ -186,6 +165,6 @@
 
             ! mass flux at the ocean/ice interface
-            fmmflx(ji,jj) = - wfx_ice(ji,jj) * rdt_ice                   ! F/M mass flux save at least for biogeochemical model
-            emp(ji,jj)    = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
+            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) ) * r1_rdtice  ! F/M mass flux save at least for biogeochemical model
+            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
             
          END DO
@@ -205,5 +184,5 @@
          snwice_mass_b(:,:) = snwice_mass(:,:)                  
          ! new mass per unit area
-         snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
+         snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
          ! time evolution of snow+ice mass
          snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice
@@ -216,12 +195,17 @@
       tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                      
 
-      !------------------------------------------------!
-      !    Computation of snow/ice and ocean albedo    !
-      !------------------------------------------------!
-      IF( lk_cpl ) THEN          ! coupled case
-         CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )                  ! snow/ice albedo
-         alb_ice(:,:,:) =  0.5_wp * zalbp(:,:,:) + 0.5_wp * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys)
-      ENDIF
-
+      !------------------------------------------------------------------------!
+      !    Snow/ice albedo (only if sent to coupler, useless in forced mode)   !
+      !------------------------------------------------------------------------!
+      CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os )    
+      CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos
+      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+      CALL wrk_dealloc( jpi, jpj, jpl, zalb_cs, zalb_os )
+
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_final( 'limsbc' )
+
+      ! control prints
+      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 3, ' - Final state lim_sbc - ' )
 
       IF(ln_ctl) THEN
@@ -231,7 +215,5 @@
          CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
       ENDIF
-      !
-      IF( lk_cpl )   CALL wrk_dealloc( jpi, jpj, jpl, zalb, zalbp )
-      ! 
+
    END SUBROUTINE lim_sbc_flx
 
@@ -271,7 +253,5 @@
       !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
-!CDIR NOVERRCHK
          DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
-!CDIR NOVERRCHK
             DO ji = fs_2, fs_jpim1
                !                                               ! 2*(U_ice-U_oce) at T-point
@@ -323,6 +303,4 @@
       !! ** input   : Namelist namicedia
       !!-------------------------------------------------------------------
-      REAL(wp) :: zsum, zarea
-      !
       INTEGER  ::   ji, jj, jk                       ! dummy loop indices
       REAL(wp) ::   zcoefu, zcoefv, zcoeff          ! local scalar
@@ -344,15 +322,9 @@
          END WHERE
       ENDIF
-      ! clem modif
-      IF( .NOT. ln_rstart ) THEN
-         iatte(:,:) = 1._wp
-         oatte(:,:) = 1._wp
-      ENDIF
-      !
-      ! clem: snwice_mass in the restart file now
+      !
       IF( .NOT. ln_rstart ) THEN
          !                                      ! embedded sea ice
          IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
-            snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )
+            snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )
             snwice_mass_b(:,:) = snwice_mass(:,:)
          ELSE
