Index: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4763)
+++ branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5208)
@@ -32,6 +32,5 @@
    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 : fraqsr_1lev, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass
    USE albedo           ! albedo parameters
    USE lbclnk           ! ocean lateral boundary condition - MPP exchanges
@@ -51,7 +50,4 @@
    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 (lk_cpl=T)
       !!
       !! 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) ::   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_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(:,:) * pfrld(:,:) )   !     solar flux at ocean surface
+      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) )   ! 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 ) ) ! 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(:,:) + qns(:,:) ) * pfrld(:,:) )  
+      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * a_i_b(:,:,:), dim=3 ) )
 
       ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
@@ -130,21 +128,17 @@
             !      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) 
             !---------------------------------------------------
-            IF( lk_cpl ) THEN ! be carfeful: 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)
-               !!!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)
+                  zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * a_i_b(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) )
+               !!! LIM2 version zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(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)
+                  zfcm1   = zfcm1 + a_i_b(ji,jj,jl) * ftr_ice(ji,jj,jl)
                END DO
             ENDIF
@@ -182,9 +176,10 @@
 
             ! 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) + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj)
+            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
+                           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj)
 
             ! 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) * r1_rdtice                    ! F/M mass flux save at least for biogeochemical model
+            emp(ji,jj)    = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj)       ! mass flux + F/M mass flux (always ice/ocean mass exchange)
             
          END DO
@@ -194,5 +189,6 @@
       !      salt flux at the ocean surface      !
       !------------------------------------------!
-      sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:)
+      sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:)   &
+         &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:)
 
       !-------------------------------------------------------------!
@@ -215,9 +211,16 @@
 
       !------------------------------------------------!
-      !    Computation of snow/ice and ocean albedo    !
+      !    Snow/ice albedo (only if sent to coupler)   !
       !------------------------------------------------!
       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)
+
+            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 )
+
       ENDIF
 
@@ -229,7 +232,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
 
@@ -344,6 +345,5 @@
       ! clem modif
       IF( .NOT. ln_rstart ) THEN
-         iatte(:,:) = 1._wp
-         oatte(:,:) = 1._wp
+         fraqsr_1lev(:,:) = 1._wp
       ENDIF
       !
