Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 6455)
@@ -9,4 +9,5 @@
    !!             -   ! 2001-06  (M. Vancoppenolle) LIM 3.0
    !!             -   ! 2006-08  (G. Madec)  cleaning for surface module
+   !!            3.6  ! 2016-01  (C. Rousset) new parameterization for sea ice albedo
    !!----------------------------------------------------------------------
 
@@ -29,21 +30,14 @@
 
    INTEGER  ::   albd_init = 0      !: control flag for initialization
-   REAL(wp) ::   zzero     = 0.e0   ! constant values
-   REAL(wp) ::   zone      = 1.e0   !    "       "
-
-   REAL(wp) ::   c1     = 0.05    ! constants values
-   REAL(wp) ::   c2     = 0.10    !    "        "
-   REAL(wp) ::   rmue   = 0.40    !  cosine of local solar altitude
-
+  
+   REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude
+   REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
+   REAL(wp) ::   c1       = 0.05    ! snow thickness (only for nn_ice_alb=0)
+   REAL(wp) ::   c2       = 0.10    !  "        "
+   REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
+ 
    !                             !!* namelist namsbc_alb
-   REAL(wp) ::   rn_cloud         !  cloudiness effect on snow or ice albedo (Grenfell & Perovich, 1984)
-#if defined key_lim3
-   REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-#else
-   REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-#endif
-   REAL(wp) ::   rn_alphd         !  coefficients for linear interpolation used to compute
-   REAL(wp) ::   rn_alphdi        !  albedo between two extremes values (Pyane, 1972)
-   REAL(wp) ::   rn_alphc         ! 
+   INTEGER  ::   nn_ice_alb
+   REAL(wp) ::   rn_albice
 
    !!----------------------------------------------------------------------
@@ -59,11 +53,31 @@
       !!          
       !! ** Purpose :   Computation of the albedo of the snow/ice system 
-      !!                as well as the ocean one
       !!       
-      !! ** Method  : - Computation of the albedo of snow or ice (choose the 
-      !!                rignt one by a large number of tests
-      !!              - Computation of the albedo of the ocean
-      !!
-      !! References :   Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
+      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
+      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
+      !!                                                                           and Grenfell & Perovich (JGR 2004)
+      !!                Description of scheme 1:
+      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
+      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
+      !!                     0-5cm  : linear function of ice thickness
+      !!                     5-150cm: log    function of ice thickness
+      !!                     > 150cm: constant
+      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004)
+      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting
+      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004)
+      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law
+      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
+      !!
+      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
+      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
+      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
+      !!                     under melting conditions than under freezing conditions
+      !!                  3) the evolution of ice albedo as a function of ice thickness shows  
+      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
+      !!
+      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !!                Brandt et al. 2005, J. Climate, vol 18
+      !!                Grenfell & Perovich 2004, JGR, vol 109 
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
@@ -73,99 +87,133 @@
       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky
       !!
-      INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      INTEGER  ::   ijpl          ! number of ice categories (3rd dim of ice input arrays)
-      REAL(wp) ::   zalbpsnm      ! albedo of ice under clear sky when snow is melting
-      REAL(wp) ::   zalbpsnf      ! albedo of ice under clear sky when snow is freezing
-      REAL(wp) ::   zalbpsn       ! albedo of snow/ice system when ice is coverd by snow
-      REAL(wp) ::   zalbpic       ! albedo of snow/ice system when ice is free of snow
-      REAL(wp) ::   zithsn        ! = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow)
-      REAL(wp) ::   zitmlsn       ! = 1 freezinz snow (pt_ice >=rt0_snow) ; = 0 melting snow (pt_ice<rt0_snow)
-      REAL(wp) ::   zihsc1        ! = 1 hsn <= c1 ; = 0 hsn > c1
-      REAL(wp) ::   zihsc2        ! = 1 hsn >= c2 ; = 0 hsn < c2
-      !!
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalbfz    ! = rn_alphdi for freezing ice ; = rn_albice for melting ice
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zficeth   !  function of ice thickness
+      INTEGER  ::   ji, jj, jl         ! dummy loop indices
+      INTEGER  ::   ijpl               ! number of ice categories (3rd dim of ice input arrays)
+      REAL(wp)            ::   ralb_im, ralb_sf, ralb_sm, ralb_if
+      REAL(wp)            ::   zswitch, z1_c1, z1_c2
+      REAL(wp)                            ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free)
       !!---------------------------------------------------------------------
-      
+
       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
-
-      CALL wrk_alloc( jpi,jpj,ijpl, zalbfz, zficeth )
+      
+      CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
 
       IF( albd_init == 0 )   CALL albedo_init      ! initialization 
 
-      !---------------------------
-      !  Computation of  zficeth
-      !---------------------------
-      ! ice free of snow and melts
-      WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalbfz(:,:,:) = rn_albice
-      ELSE WHERE                                              ;   zalbfz(:,:,:) = rn_alphdi
-      END  WHERE
-
-      WHERE     ( 1.5  < ph_ice                     )  ;  zficeth = zalbfz
-      ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 )
-      ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              &
-         &                                                                 - 0.8608 * ph_ice * ph_ice     &
-         &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
-      ELSE WHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice
-      END WHERE
-
-!!gm old code
-!      DO jl = 1, ijpl
-!         DO jj = 1, jpj
-!            DO ji = 1, jpi
-!               IF( ph_ice(ji,jj,jl) > 1.5 ) THEN
-!                  zficeth(ji,jj,jl) = zalbfz(ji,jj,jl)
-!               ELSEIF( ph_ice(ji,jj,jl) > 1.0  .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN
-!                  zficeth(ji,jj,jl) = 0.472 + 2.0 * ( zalbfz(ji,jj,jl) - 0.472 ) * ( ph_ice(ji,jj,jl) - 1.0 )
-!               ELSEIF( ph_ice(ji,jj,jl) > 0.05 .AND. ph_ice(ji,jj,jl) <= 1.0 ) THEN
-!                  zficeth(ji,jj,jl) = 0.2467 + 0.7049 * ph_ice(ji,jj,jl)                               &
-!                     &                    - 0.8608 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl)                 &
-!                     &                    + 0.3812 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl) * ph_ice (ji,jj,jl)
-!               ELSE
-!                  zficeth(ji,jj,jl) = 0.1 + 3.6 * ph_ice(ji,jj,jl) 
-!               ENDIF
-!            END DO
-!         END DO
-!      END DO
-!!gm end old code
-      
-      !----------------------------------------------- 
-      !    Computation of the snow/ice albedo system 
-      !-------------------------- ---------------------
-      
-      !    Albedo of snow-ice for clear sky.
-      !-----------------------------------------------    
-      DO jl = 1, ijpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !  Case of ice covered by snow.             
-               !                                        !  freezing snow        
-               zihsc1   = 1.0 - MAX( zzero , SIGN( zone , - ( ph_snw(ji,jj,jl) - c1 ) ) )
-               zalbpsnf = ( 1.0 - zihsc1 ) * (  zficeth(ji,jj,jl)                                             &
-                  &                           + ph_snw(ji,jj,jl) * ( rn_alphd - zficeth(ji,jj,jl) ) / c1  )   &
-                  &     +         zihsc1   * rn_alphd  
-               !                                        !  melting snow                
-               zihsc2   = MAX( zzero , SIGN( zone , ph_snw(ji,jj,jl) - c2 ) )
-               zalbpsnm = ( 1.0 - zihsc2 ) * ( rn_albice + ph_snw(ji,jj,jl) * ( rn_alphc - rn_albice ) / c2 )   &
-                  &     +         zihsc2   *   rn_alphc 
-               !
-               zitmlsn  =  MAX( zzero , SIGN( zone , pt_ice(ji,jj,jl) - rt0_snow ) )   
-               zalbpsn  =  zitmlsn * zalbpsnm + ( 1.0 - zitmlsn ) * zalbpsnf
-            
-               !  Case of ice free of snow.
-               zalbpic  = zficeth(ji,jj,jl) 
-            
-               ! albedo of the system   
-               zithsn   = 1.0 - MAX( zzero , SIGN( zone , - ph_snw(ji,jj,jl) ) )
-               pa_ice_cs(ji,jj,jl) =  zithsn * zalbpsn + ( 1.0 - zithsn ) *  zalbpic
+      
+      SELECT CASE ( nn_ice_alb )
+
+      !------------------------------------------
+      !  Shine and Henderson-Sellers (1985)
+      !------------------------------------------
+      CASE( 0 )
+       
+         ralb_sf = 0.80       ! dry snow
+         ralb_sm = 0.65       ! melting snow
+         ralb_if = 0.72       ! bare frozen ice
+         ralb_im = rn_albice  ! bare puddled ice 
+         
+         !  Computation of ice albedo (free of snow)
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
+         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
+         END  WHERE
+      
+         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
+         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              &
+            &                                                                 - 0.8608 * ph_ice * ph_ice     &
+            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
+         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
+         END WHERE
+     
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! freezing snow
+                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
+                  !                                        !  freezing snow        
+                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
+                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
+                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
+                     &        +         zswitch   * ralb_sf  
+
+                  ! melting snow
+                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
+                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
+                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
+                      &     +         zswitch   *   ralb_sm 
+                  !
+                  ! snow albedo
+                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+               
+                  ! Ice/snow albedo
+                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_cs(ji,jj,jl) =  zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
+                  !
+               END DO
             END DO
          END DO
-      END DO
-      
-      !    Albedo of snow-ice for overcast sky.
-      !----------------------------------------------  
-      pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rn_cloud       ! Oberhuber correction
-      !
-      CALL wrk_dealloc( jpi,jpj,ijpl, zalbfz, zficeth )
+
+         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
+
+      !------------------------------------------
+      !  New parameterization (2016)
+      !------------------------------------------
+      CASE( 1 ) 
+
+         ralb_im = rn_albice  ! bare puddled ice
+! compilation of values from literature
+         ralb_sf = 0.85      ! dry snow
+         ralb_sm = 0.75      ! melting snow
+         ralb_if = 0.60      ! bare frozen ice
+! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
+!         ralb_sf = 0.85       ! dry snow
+!         ralb_sm = 0.72       ! melting snow
+!         ralb_if = 0.65       ! bare frozen ice
+! Brandt et al 2005 (East Antarctica)
+!         ralb_sf = 0.87      ! dry snow
+!         ralb_sm = 0.82      ! melting snow
+!         ralb_if = 0.54      ! bare frozen ice
+! 
+         !  Computation of ice albedo (free of snow)
+         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+         z1_c2 = 1. / 0.05
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im
+         ELSE WHERE                                              ;   zalb = ralb_if
+         END  WHERE
+         
+         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
+            &                                                                     ( LOG(1.5) - LOG(ph_ice) )
+         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice
+         END WHERE
+
+         z1_c1 = 1. / 0.02
+         z1_c2 = 1. / 0.03
+         !  Computation of the snow/ice albedo
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
+                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
+
+                   ! snow albedo
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+
+                  ! Ice/snow albedo   
+                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
+
+              END DO
+            END DO
+         END DO
+         ! Effect of the clouds (2d order polynomial)
+         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
+
+      END SELECT
+      
+      CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
       !
    END SUBROUTINE albedo_ice
@@ -181,10 +229,10 @@
       REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
       !!
-      REAL(wp) ::   zcoef   ! local scalar
-      !!----------------------------------------------------------------------
-      !
-      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )      ! Parameterization of Briegled and Ramanathan, 1982 
-      pa_oce_cs(:,:) = zcoef               
-      pa_oce_os(:,:)  = 0.06                         ! Parameterization of Kondratyev, 1969 and Payne, 1972
+      REAL(wp) :: zcoef 
+      !!----------------------------------------------------------------------
+      !
+      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
+      pa_oce_cs(:,:) = zcoef 
+      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
       !
    END SUBROUTINE albedo_oce
@@ -200,5 +248,5 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ rn_cloud, rn_albice, rn_alphd, rn_alphdi, rn_alphc
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
       !!----------------------------------------------------------------------
       !
@@ -219,9 +267,6 @@
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*) '   Namelist namsbc_alb : albedo '
-         WRITE(numout,*) '      correction for snow and ice albedo                  rn_cloud  = ', rn_cloud
-         WRITE(numout,*) '      albedo of melting ice in the arctic and antarctic   rn_albice = ', rn_albice
-         WRITE(numout,*) '      coefficients for linear                             rn_alphd  = ', rn_alphd
-         WRITE(numout,*) '      interpolation used to compute albedo                rn_alphdi = ', rn_alphdi
-         WRITE(numout,*) '      between two extremes values (Pyane, 1972)           rn_alphc  = ', rn_alphc
+         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb
+         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice
       ENDIF
       !
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 6455)
@@ -80,5 +80,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_oce       !: heat flux of precip and evap over ocean     [W/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_ice       !: heat flux of precip and evap over ice       [W/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: heat flux of precip over ice                [J/m3]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qevap_ice      !: heat flux of evap over ice                  [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: enthalpy of precip over ice                 [J/m3]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s]
 #endif
@@ -144,7 +145,7 @@
 #endif
 #if defined key_lim3
-         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,  &
-         &      qemp_ice(jpi,jpj)     , qemp_oce(jpi,jpj)      ,                       &
-         &      qns_oce (jpi,jpj)     , qsr_oce (jpi,jpj)      , emp_oce (jpi,jpj)  ,  &
+         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   &
+         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   &
+         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce (jpi,jpj)  ,   &
 #endif
          &      emp_ice(jpi,jpj)      ,  STAT= ierr(1) )
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 6455)
@@ -684,4 +684,10 @@
       qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
 
+      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
+      DO jl = 1, jpl
+         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) - lfus )
+                                   ! but then qemp_ice should also include sublimation 
+      END DO
+
       CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
 #endif
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 6455)
@@ -612,7 +612,7 @@
       ! --- evaporation --- !
       z1_lsub = 1._wp / Lsub
-      evap_ice (:,:,:) = qla_ice (:,:,:) * z1_lsub ! sublimation
-      devap_ice(:,:,:) = dqla_ice(:,:,:) * z1_lsub
-      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean
+      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation
+      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT
+      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean
 
       ! --- evaporation minus precipitation --- !
@@ -637,4 +637,10 @@
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
       qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+
+      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
+      DO jl = 1, jpl
+         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
+                                   ! But we do not have Tice => consider it at 0°C => evap=0 
+      END DO
 
       CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6455)
@@ -1378,14 +1378,16 @@
       !
       INTEGER ::   jl         ! dummy loop index
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk, zsnw
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
       !
-      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
-      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
+      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
+      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
+      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
 
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -1423,8 +1425,57 @@
       END SELECT
 
-      IF( iom_use('subl_ai_cea') )   &
-         CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
-      !   
-      !                                                           ! runoffs and calving (put in emp_tot)
+#if defined key_lim3
+      ! zsnw = snow percentage over ice after wind blowing
+      zsnw(:,:) = 0._wp
+      CALL lim_thd_snwblow( p_frld, zsnw )
+      
+      ! --- evaporation (kg/m2/s) --- !
+      zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1)
+      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
+      ! therefore, sublimation is not redistributed over the ice categories in case no subgrid scale fluxes are provided by atm.
+      zdevap_ice(:,:) = 0._wp
+      
+      ! --- evaporation minus precipitation corrected for the effect of wind blowing on snow --- !
+      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) - zsprecip * (1._wp - zsnw)
+      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip * (1._wp - zsnw)          
+
+      ! Sublimation over sea-ice (cell average)
+      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:) )
+      ! runoffs and calving (put in emp_tot)
+      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+      IF( srcv(jpr_cal)%laction ) THEN 
+         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
+         CALL iom_put( 'calving_cea', frcv(jpr_cal)%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(:,:)
+         emp_oce(:,:) = emp_oce(:,:) * xcplmask(:,:,0) + zemp_oce(:,:) * zmsk(:,:)
+         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
+         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
+         DO jl=1,jpl
+            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:)
+            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:)
+         ENDDO
+      ELSE
+         emp_tot(:,:) =         zemp_tot(:,:)
+         emp_ice(:,:) =         zemp_ice(:,:)
+         emp_oce(:,:) =         zemp_oce(:,:)     
+         sprecip(:,:) =         zsprecip(:,:)
+         tprecip(:,:) =         ztprecip(:,:)
+         DO jl=1,jpl
+            evap_ice (:,:,jl) = zevap_ice (:,:)
+            devap_ice(:,:,jl) = zdevap_ice(:,:)
+         ENDDO
+      ENDIF
+
+                                     CALL iom_put( 'snowpre'    , sprecip                         )  ! Snow
+      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw ) )  ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw   )  ! Snow over sea-ice         (cell average)    
+#else
+      ! Sublimation over sea-ice (cell average)
+      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )
+      ! runoffs and calving (put in emp_tot)
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
       IF( srcv(jpr_cal)%laction ) THEN 
@@ -1450,4 +1501,5 @@
       IF( iom_use('snow_ai_cea') )   &
          CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average)
+#endif
 
       !                                                      ! ========================= !
@@ -1505,18 +1557,7 @@
       IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
 
-#if defined key_lim3
-      CALL wrk_alloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
-
+#if defined key_lim3      
       ! --- evaporation --- !
-      ! clem: evap_ice is set to 0 for LIM3 since we still do not know what to do with sublimation
-      ! the problem is: the atm. imposes both mass evaporation and heat removed from the snow/ice
-      !                 but it is incoherent WITH the ice model  
-      DO jl=1,jpl
-         evap_ice(:,:,jl) = 0._wp  ! should be: frcv(jpr_ievp)%z3(:,:,1)
-      ENDDO
       zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean
-
-      ! --- evaporation minus precipitation --- !
-      emp_oce(:,:) = emp_tot(:,:) - emp_ice(:,:)
 
       ! --- non solar flux over ocean --- !
@@ -1525,18 +1566,23 @@
       WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
 
-      ! --- heat flux associated with emp --- !
-      zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing
+      ! --- heat flux associated with emp (W/m2) --- !
       zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
          &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip
          &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean
-      qemp_ice(:,:)  = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
-         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
-
+!      zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
+!         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
+      zqemp_ice(:,:) =      zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice (only)
+                                                                                                       ! qevap_ice=0 since we consider Tice=0°C
+      
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
       zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
 
-      ! --- total non solar flux --- !
-      zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:)
+      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
+      DO jl = 1, jpl
+         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but we do not have Tice, so we consider Tice=0°C
+      END DO
+
+      ! --- total non solar flux (including evap/precip) --- !
+      zqns_tot(:,:) = zqns_tot(:,:) + zqemp_ice(:,:) + zqemp_oce(:,:)
 
       ! --- in case both coupled/forced are active, we must mix values --- ! 
@@ -1545,20 +1591,20 @@
          qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
          DO jl=1,jpl
-            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
+            qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:)
+            qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:)
          ENDDO
          qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
          qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
-!!clem         evap_ice(:,:) = evap_ice(:,:) * xcplmask(:,:,0)
+         qemp_ice (:,:) =  qemp_ice(:,:) * xcplmask(:,:,0) +  zqemp_ice(:,:)* zmsk(:,:)
       ELSE
          qns_tot  (:,:  ) = zqns_tot  (:,:  )
          qns_oce  (:,:  ) = zqns_oce  (:,:  )
          qns_ice  (:,:,:) = zqns_ice  (:,:,:)
-         qprec_ice(:,:)   = zqprec_ice(:,:)
-         qemp_oce (:,:)   = zqemp_oce (:,:)
-      ENDIF
-
-      CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
+         qevap_ice(:,:,:) = zqevap_ice(:,:,:)
+         qprec_ice(:,:  ) = zqprec_ice(:,:  )
+         qemp_oce (:,:  ) = zqemp_oce (:,:  )
+         qemp_ice (:,:  ) = zqemp_ice (:,:  )
+      ENDIF
 #else
-
       ! clem: this formulation is certainly wrong... but better than it was...
       zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with:
@@ -1577,5 +1623,4 @@
          qns_ice(:,:,:) = zqns_ice(:,:,:)
       ENDIF
-
 #endif
 
@@ -1628,5 +1673,4 @@
 
 #if defined key_lim3
-      CALL wrk_alloc( jpi,jpj, zqsr_oce ) 
       ! --- solar flux over ocean --- !
       !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
@@ -1636,6 +1680,4 @@
       IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
       ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
-
-      CALL wrk_dealloc( jpi,jpj, zqsr_oce ) 
 #endif
 
@@ -1688,6 +1730,8 @@
       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
 
-      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
-      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
+      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
+      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
+      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6455)
@@ -110,5 +110,4 @@
       INTEGER  ::   jl                 ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled)
       REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice 
       !!----------------------------------------------------------------------
@@ -197,5 +196,5 @@
          ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%]
          !----------------------------------------------------------------------------------------
-         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs )
          CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
 
@@ -203,24 +202,21 @@
          CASE( jp_clio )                                       ! CLIO bulk formulation
             ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 
-            ! (zalb_ice) is computed within the bulk routine
-            CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, zalb_ice )
-            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
-            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            ! (alb_ice) is computed within the bulk routine
+                                 CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice )
+            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          CASE( jp_core )                                       ! CORE bulk formulation
             ! albedo depends on cloud fraction because of non-linear spectral effects
-            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-            CALL blk_ice_core_flx( t_su, zalb_ice )
-            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
-            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+                                 CALL blk_ice_core_flx( t_su, alb_ice )
+            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          CASE ( jp_purecpl )
             ! albedo depends on cloud fraction because of non-linear spectral effects
-            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
-            ! clem: evap_ice is forced to 0 in coupled mode for now 
-            !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models
-            evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp
-            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          END SELECT
-         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs )
 
          !----------------------------!
@@ -265,4 +261,5 @@
       !!----------------------------------------------------------------------
       INTEGER :: ierr
+      INTEGER :: ji, jj
       !!----------------------------------------------------------------------
       IF(lwp) WRITE(numout,*)
@@ -321,4 +318,12 @@
       tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
       !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( gphit(ji,jj) > 0._wp ) THEN  ;  rn_amax_2d(ji,jj) = rn_amax_n  ! NH
+            ELSE                             ;  rn_amax_2d(ji,jj) = rn_amax_s  ! SH
+            ENDIF
+        ENDDO
+      ENDDO 
+      !
       nstart = numit  + nn_fsbc      
       nitrun = nitend - nit000 + 1 
@@ -343,5 +348,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir,  &
-         &                ln_limdyn, rn_amax, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt  
+         &                ln_limdyn, rn_amax_n, rn_amax_s, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt  
       !!-------------------------------------------------------------------
       !                    
@@ -364,5 +369,6 @@
          WRITE(numout,*) '   number of snow layers                                   = ', nlay_s
          WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
-         WRITE(numout,*) '   maximum ice concentration                               = ', rn_amax 
+         WRITE(numout,*) '   maximum ice concentration for NH                        = ', rn_amax_n 
+         WRITE(numout,*) '   maximum ice concentration for SH                        = ', rn_amax_s
          WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb
          WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout
@@ -579,5 +585,5 @@
       sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
       sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
-      sfx_res(:,:) = 0._wp
+      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
       
       wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
@@ -595,6 +601,7 @@
       hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
       hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
-      hfx_err_dif(:,:) = 0._wp   ;
-
+      hfx_err_dif(:,:) = 0._wp
+      wfx_err_sub(:,:) = 0._wp
+      
       afx_tot(:,:) = 0._wp   ;
       afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6455)
@@ -456,5 +456,6 @@
       !                                                ! ---------------------------------------- !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-         CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux
+         CALL iom_put( "empmr"  , emp    - rnf )                ! upward water flux
+         CALL iom_put( "empbmr" , emp_b  - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline )
          CALL iom_put( "saltflx", sfx  )                        ! downward salt flux  
                                                                 ! (includes virtual salt flux beneath ice 
