Index: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4763)
+++ branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5208)
@@ -12,4 +12,5 @@
    !!            3.4  ! 2011-01  (A Porter)  dynamical allocation
    !!             -   ! 2012-10  (C. Rousset)  add lim_diahsb
+   !!            3.6  ! 2014-07  (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -59,5 +60,4 @@
    USE prtctl          ! Print control
    USE lib_fortran     ! 
-   USE cpl_oasis3, ONLY : lk_cpl
 
 #if defined key_bdy 
@@ -80,33 +80,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-   FUNCTION fice_cell_ave ( ptab)
-      !!--------------------------------------------------------------------------
-      !! * Compute average over categories, for grid cell (ice covered and free ocean)
-      !!--------------------------------------------------------------------------
-      REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
-      REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
-      INTEGER :: jl ! Dummy loop index
-      
-      fice_cell_ave (:,:) = 0.0_wp
-      
-      DO jl = 1, jpl
-         fice_cell_ave (:,:) = fice_cell_ave (:,:) &
-            &                  + a_i (:,:,jl) * ptab (:,:,jl)
-      END DO
-      
-   END FUNCTION fice_cell_ave
-   
-   FUNCTION fice_ice_ave ( ptab)
-      !!--------------------------------------------------------------------------
-      !! * Compute average over categories, for ice covered part of grid cell
-      !!--------------------------------------------------------------------------
-      REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
-      REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
-
-      fice_ice_ave (:,:) = 0.0_wp
-      WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
-
-   END FUNCTION fice_ice_ave
 
    !!======================================================================
@@ -133,33 +104,13 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt      ! ocean time step
-      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE)
+      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE, =5 COUPLED)
       !!
-      INTEGER  ::   ji, jj, jl, jk      ! dummy loop index
+      INTEGER  ::   jl      ! dummy loop index
       REAL(wp) ::   zcoef   ! local scalar
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice (for coupled)
-
-      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all    ! Mean albedo over all categories
-      REAL(wp), POINTER, DIMENSION(:,:) :: ztem_ice_all    ! Mean temperature over all categories
-      
-      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_ice_all   ! Mean solar heat flux over all categories
-      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_ice_all   ! Mean non solar heat flux over all categories
-      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_ice_all   ! Mean latent heat flux over all categories
-      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories
-      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories
-      REAL(wp) ::   ztmelts           ! clem 2014: for HC diags
-      REAL(wp) ::   epsi20 = 1.e-20   !
+      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)
       !!----------------------------------------------------------------------
 
-      !- O.M. : why do we allocate all these arrays even when MOD( kt-1, nn_fsbc ) /= 0 ?????
-
       IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim')
-
-      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice )
-
-      IF( lk_cpl ) THEN
-         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
-            &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
-      ENDIF
 
       IF( kt == nit000 ) THEN
@@ -182,28 +133,30 @@
          !                                           !----------------!
          !
-         u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
-         v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean)
-
-         ! masked sea surface freezing temperature [Kelvin]
-         t_bo(:,:) = ( tfreez( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )
-
-         CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo
-
+         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)                     ! mean surface ocean current at ice velocity point
+         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)                    ! (C-grid dynamics :  U- & V-points as the ocean)
+         !
+         t_bo(:,:) = ( eos_fzp( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )  ! masked sea surface freezing temperature [Kelvin]
+         !                                                                                  ! (set to rt0 over land)
+         !                                           ! Ice albedo
+         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )      
+
+         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos
+
+         SELECT CASE( kblk )
+         CASE( jp_core , jp_cpl )   ! CORE and COUPLED bulk formulations
+
+            ! albedo depends on cloud fraction because of non-linear spectral effects
+            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 
+            ! (zalb_ice) is computed within the bulk routine
+            
+         END SELECT
+         
+         !                                           ! Mask sea ice surface temperature
          DO jl = 1, jpl
             t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) )
          END DO
-
-         IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
-         
-         IF( lk_cpl ) THEN
-            IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
-               !
-               ! Compute mean albedo and temperature
-               zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) ) 
-               ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) ) 
-               !
-            ENDIF
-         ENDIF
-                                               ! Bulk formulea - provides the following fields:
+     
+         ! Bulk formulae  - provides the following fields:
          ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2]
          ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2]
@@ -214,6 +167,6 @@
          !
          SELECT CASE( kblk )
-         CASE( 3 )                                       ! CLIO bulk formulation
-            CALL blk_ice_clio( t_su , zalb_ice_cs, zalb_ice_os,                           &
+         CASE( jp_clio )                                       ! CLIO bulk formulation
+            CALL blk_ice_clio( t_su , zalb_cs    , zalb_os    , zalb_ice  ,               &
                &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   &
                &                      qla_ice    , dqns_ice   , dqla_ice  ,               &
@@ -221,13 +174,8 @@
                &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  )
             !         
-         CASE( 4 )                                       ! CORE bulk formulation
-            ! MV 2014
-            ! We must account for cloud fraction in the computation of the albedo
-            ! The present ref just uses the clear sky value
-            ! The overcast sky value is 0.06 higher, and polar skies are mostly overcast
-            ! CORE has no cloud fraction, hence we must prescribe it
-            ! Mean summer cloud fraction computed from CLIO = 0.81
-            zalb_ice(:,:,:) = 0.19 * zalb_ice_cs(:,:,:) + 0.81 * zalb_ice_os(:,:,:)
-            ! Following line, we replace zalb_ice_cs by simply zalb_ice
+            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
+               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
+
+         CASE( jp_core )                                       ! CORE bulk formulation
             CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               &
                &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   &
@@ -235,51 +183,25 @@
                &                      tprecip   , sprecip   ,                            &
                &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
+               !
+            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
+               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
             !
-         CASE ( 5 )
-            zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
+         CASE ( jp_cpl )
             
             CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
 
-            CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=tn_ice    )
-
-            ! Latent heat flux is forced to 0 in coupled :
-            !  it is included in qns (non-solar heat flux)
-            qla_ice  (:,:,:) = 0.0e0_wp
-            dqla_ice (:,:,:) = 0.0e0_wp
+            ! MV -> seb 
+!           CALL sbc_cpl_ice_flx( p_frld=ato_i, 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, qla_ice, dqla_ice, nn_limflx )
+!           ! Latent heat flux is forced to 0 in coupled :
+!           !  it is included in qns (non-solar heat flux)
+!           qla_ice  (:,:,:) = 0._wp
+!           dqla_ice (:,:,:) = 0._wp
+            ! END MV -> seb
             !
          END SELECT
-
-         ! Average over all categories
-         IF( lk_cpl ) THEN
-         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
-
-            z_qns_ice_all  (:,:) = fice_ice_ave ( qns_ice  (:,:,:) )
-            z_qsr_ice_all  (:,:) = fice_ice_ave ( qsr_ice  (:,:,:) )
-            z_dqns_ice_all (:,:) = fice_ice_ave ( dqns_ice (:,:,:) )
-            z_qla_ice_all  (:,:) = fice_ice_ave ( qla_ice  (:,:,:) )
-            z_dqla_ice_all (:,:) = fice_ice_ave ( dqla_ice (:,:,:) )
-
-            DO jl = 1, jpl
-               dqns_ice (:,:,jl) = z_dqns_ice_all (:,:)
-               dqla_ice (:,:,jl) = z_dqla_ice_all (:,:)
-            END DO
-            !
-            IF ( ln_iceflx_ave ) THEN
-               DO jl = 1, jpl
-                  qns_ice  (:,:,jl) = z_qns_ice_all  (:,:)
-                  qsr_ice  (:,:,jl) = z_qsr_ice_all  (:,:)
-                  qla_ice  (:,:,jl) = z_qla_ice_all  (:,:)
-               END DO
-            END IF
-            !
-            IF ( ln_iceflx_linear ) THEN
-               DO jl = 1, jpl
-                  qns_ice  (:,:,jl) = z_qns_ice_all(:,:) + z_dqns_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:))
-                  qla_ice  (:,:,jl) = z_qla_ice_all(:,:) + z_dqla_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:))
-                  qsr_ice  (:,:,jl) = (1.0e0_wp-zalb_ice(:,:,jl)) / (1.0e0_wp-zalb_ice_all(:,:)) * z_qsr_ice_all(:,:)
-               END DO
-            END IF
-         END IF
-         ENDIF
+         
          !                                           !----------------------!
          !                                           ! LIM-3  time-stepping !
@@ -289,28 +211,17 @@
          !
          !                                           ! Store previous ice values
-!!gm : remark   old_...   should becomes ...b  as tn versus tb  
-         old_a_i  (:,:,:)   = a_i  (:,:,:)     ! ice area
-         old_e_i  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
-         old_v_i  (:,:,:)   = v_i  (:,:,:)     ! ice volume
-         old_v_s  (:,:,:)   = v_s  (:,:,:)     ! snow volume 
-         old_e_s  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
-         old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content
-         old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content
-         old_u_ice(:,:)     = u_ice(:,:)
-         old_v_ice(:,:)     = v_ice(:,:)
-
-         ! trends    !!gm is it truly necessary ???
-         d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp
-         d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp
-         d_e_i_thd  (:,:,:,:) = 0._wp   ;   d_e_i_trp  (:,:,:,:) = 0._wp
-         d_v_s_thd  (:,:,:)   = 0._wp   ;   d_v_s_trp  (:,:,:)   = 0._wp
-         d_e_s_thd  (:,:,:,:) = 0._wp   ;   d_e_s_trp  (:,:,:,:) = 0._wp
-         d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp
-         d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp
-         d_u_ice_dyn(:,:)     = 0._wp   ;   d_v_ice_dyn(:,:)     = 0._wp
+         a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area
+         e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
+         v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume
+         v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume 
+         e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
+         smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content
+         oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content
+         u_ice_b(:,:)     = u_ice(:,:)
+         v_ice_b(:,:)     = v_ice(:,:)
 
          ! salt, heat and mass fluxes
          sfx    (:,:) = 0._wp   ;
-         sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
+         sfx_bri(:,:) = 0._wp   ; 
          sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
          sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
@@ -333,17 +244,4 @@
          hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
          hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
-
-         !
-         fhld  (:,:)    = 0._wp 
-         fmmflx(:,:)    = 0._wp     
-         ! part of solar radiation transmitted through the ice
-         ftr_ice(:,:,:) = 0._wp
-
-         ! diags
-         diag_trp_vi  (:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp
-         diag_heat_dhc(:,:) = 0._wp  
-
-         ! dynamical invariants
-         delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp
 
                           CALL lim_rst_opn( kt )     ! Open Ice restart file
@@ -371,13 +269,13 @@
          ENDIF
 !                         !- Change old values for new values
-                          old_u_ice(:,:)   = u_ice (:,:)
-                          old_v_ice(:,:)   = v_ice (:,:)
-                          old_a_i(:,:,:)   = a_i (:,:,:)
-                          old_v_s(:,:,:)   = v_s (:,:,:)
-                          old_v_i(:,:,:)   = v_i (:,:,:)
-                          old_e_s(:,:,:,:) = e_s (:,:,:,:)
-                          old_e_i(:,:,:,:) = e_i (:,:,:,:)
-                          old_oa_i(:,:,:)  = oa_i(:,:,:)
-                          old_smv_i(:,:,:) = smv_i (:,:,:)
+                          u_ice_b(:,:)     = u_ice(:,:)
+                          v_ice_b(:,:)     = v_ice(:,:)
+                          a_i_b  (:,:,:)   = a_i  (:,:,:)
+                          v_s_b  (:,:,:)   = v_s  (:,:,:)
+                          v_i_b  (:,:,:)   = v_i  (:,:,:)
+                          e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
+                          e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
+                          oa_i_b (:,:,:)   = oa_i (:,:,:)
+                          smv_i_b(:,:,:)   = smv_i(:,:,:)
  
          ! ----------------------------------------------
@@ -389,4 +287,17 @@
                           pfrld(:,:)   = 1._wp - at_i(:,:)
                           phicif(:,:)  = vt_i(:,:)
+
+                          ! MV -> seb
+                          SELECT CASE( kblk )
+                             CASE ( jp_cpl )
+                             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, qla_ice, dqla_ice, nn_limflx )
+                           ! Latent heat flux is forced to 0 in coupled :
+                           !  it is included in qns (non-solar heat flux)
+                             qla_ice  (:,:,:) = 0._wp
+                             dqla_ice (:,:,:) = 0._wp
+                          END SELECT
+                          ! END MV -> seb
                           !
                           CALL lim_var_bv                 ! bulk brine volume (diag)
@@ -420,4 +331,6 @@
          IF( ln_nicep )   CALL lim_ctl( kt )              ! alerts in case of model crash
          !
+         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+         !
       ENDIF                                    ! End sea-ice time step only
 
@@ -429,18 +342,92 @@
       !                                                   ! otherwise the atm.-ocean stresses are used everywhere
       IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
-      
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
-      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice )
-
-      IF( lk_cpl ) THEN
-         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
-            &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
-      ENDIF
+
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim')
       !
    END SUBROUTINE sbc_ice_lim
-
-
+   
+   
+      SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   &
+         &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_ice_lim  ***
+      !!                   
+      !! ** Purpose :   update the ice surface boundary condition by averaging and / or
+      !!                redistributing fluxes on ice categories                   
+      !!
+      !! ** Method  :   average then redistribute 
+      !!
+      !! ** Action  :   
+      !!---------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   k_limflx   ! =-1 do nothing; =0 average ; 
+                                                                ! =1 average and redistribute ; =2 redistribute
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature 
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqla_ice   ! latent heat flux
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdql_ice   ! latent heat flux sensitivity
+      !
+      INTEGER  ::   jl      ! dummy loop index
+      !
+      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_m    ! Mean albedo over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: ztem_m    ! Mean temperature over all categories
+      !
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m   ! Mean solar heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m   ! Mean non solar heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_m   ! Mean latent heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m   ! Mean d(qns)/dT over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_dql_m   ! Mean d(qla)/dT over all categories
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('ice_lim_flx')
+      !
+      !
+      SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==!
+      CASE( 0 , 1 )
+         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m)
+         !
+         z_qns_m(:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
+         z_qsr_m(:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
+         z_dqn_m(:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
+         z_qla_m(:,:) = fice_ice_ave ( pqla_ice (:,:,:) )
+         z_dql_m(:,:) = fice_ice_ave ( pdql_ice (:,:,:) )
+         DO jl = 1, jpl
+            pdqn_ice(:,:,jl) = z_dqn_m(:,:)
+            pdql_ice(:,:,jl) = z_dql_m(:,:)
+         END DO
+         !
+         DO jl = 1, jpl
+            pqns_ice(:,:,jl) = z_qns_m(:,:)
+            pqsr_ice(:,:,jl) = z_qsr_m(:,:)
+            pqla_ice(:,:,jl) = z_qla_m(:,:)
+         END DO
+         !
+         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m)
+      END SELECT
+
+      SELECT CASE( k_limflx )                              !==  redistribution on all ice categories  ==!
+      CASE( 1 , 2 )
+         CALL wrk_alloc( jpi,jpj, zalb_m, ztem_m )
+         !
+         zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) ) 
+         ztem_m(:,:) = fice_ice_ave ( ptn_ice  (:,:,:) ) 
+         DO jl = 1, jpl
+            pqns_ice(:,:,jl) = pqns_ice(:,:,jl) + pdqn_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))
+            pqla_ice(:,:,jl) = pqla_ice(:,:,jl) + pdql_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))
+            pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) ) 
+         END DO
+         !
+         CALL wrk_dealloc( jpi,jpj, zalb_m, ztem_m )
+      END SELECT
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('ice_lim_flx')
+      !
+   END SUBROUTINE ice_lim_flx
+   
+   
    SUBROUTINE lim_ctl( kt )
       !!-----------------------------------------------------------------------
@@ -472,6 +459,6 @@
                   !WRITE(numout,*) ' at_i     ', at_i(ji,jj)
                   !WRITE(numout,*) ' Point - category', ji, jj, jl
-                  !WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl)
-                  !WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl)
+                  !WRITE(numout,*) ' a_i *** a_i_b   ', a_i      (ji,jj,jl), a_i_b  (ji,jj,jl)
+                  !WRITE(numout,*) ' v_i *** v_i_b   ', v_i      (ji,jj,jl), v_i_b  (ji,jj,jl)
                   !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
                   !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
@@ -583,7 +570,7 @@
                !DO jl = 1, jpl
                   !WRITE(numout,*) ' Category no: ', jl
-                  !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)   
+                  !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' a_i_b      : ', a_i_b  (ji,jj,jl)   
                   !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
-                  !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)   
+                  !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' v_i_b      : ', v_i_b  (ji,jj,jl)   
                   !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
                   !WRITE(numout,*) ' '
@@ -674,5 +661,5 @@
       !!                n : number of the option
       !!-------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt      ! ocean time step
+      INTEGER         , INTENT(in) ::   kt            ! ocean time step
       INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices
       CHARACTER(len=*), INTENT(in) ::   cd1           !
@@ -761,5 +748,5 @@
                WRITE(numout,*) ' strength      : ', strength(ji,jj)
                WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj)
-               WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ji,jj)  , ' old_v_ice     : ', old_v_ice(ji,jj)  
+               WRITE(numout,*) ' u_ice_b       : ', u_ice_b(ji,jj)    , ' v_ice_b       : ', v_ice_b(ji,jj)  
                WRITE(numout,*)
                
@@ -771,19 +758,19 @@
                   WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl)
                   WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl)
-                  WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' old_a_i    : ', old_a_i(ji,jj,jl)   
+                  WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' a_i_b      : ', a_i_b(ji,jj,jl)   
                   WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
-                  WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' old_v_i    : ', old_v_i(ji,jj,jl)   
+                  WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' v_i_b      : ', v_i_b(ji,jj,jl)   
                   WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
-                  WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' old_v_s    : ', old_v_s(ji,jj,jl)  
+                  WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' v_s_b      : ', v_s_b(ji,jj,jl)  
                   WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl)
-                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ji,jj,1,jl)/1.0e9 
+                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' ei1        : ', e_i_b(ji,jj,1,jl)/1.0e9 
                   WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ji,jj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ji,jj,1,jl)/1.0e9
-                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ji,jj,2,jl)/1.0e9  
+                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' ei2_b      : ', e_i_b(ji,jj,2,jl)/1.0e9  
                   WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ji,jj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ji,jj,2,jl)/1.0e9
-                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' old_e_snow : ', old_e_s(ji,jj,1,jl) 
+                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl) 
                   WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ji,jj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ji,jj,1,jl)
-                  WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' old_smv_i  : ', old_smv_i(ji,jj,jl)   
+                  WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' smv_i_b    : ', smv_i_b(ji,jj,jl)   
                   WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl) 
-                  WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' old_oa_i   : ', old_oa_i(ji,jj,jl)
+                  WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' oa_i_b     : ', oa_i_b(ji,jj,jl)
                   WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl)
                END DO !jl
@@ -793,6 +780,6 @@
                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
                WRITE(numout,*) ' - Heat fluxes in and out the ice ***'
-               WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( old_a_i(ji,jj,:) * qsr_ice(ji,jj,:) )
-               WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( old_a_i(ji,jj,:) * qns_ice(ji,jj,:) )
+               WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) )
+               WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) )
                WRITE(numout,*)
                WRITE(numout,*) 
@@ -852,6 +839,38 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE lim_prt_state
+   
+     
+   FUNCTION fice_cell_ave ( ptab )
+      !!--------------------------------------------------------------------------
+      !! * Compute average over categories, for grid cell (ice covered and free ocean)
+      !!--------------------------------------------------------------------------
+      REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
+      REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
+      INTEGER :: jl ! Dummy loop index
+      
+      fice_cell_ave (:,:) = 0.0_wp
+      
+      DO jl = 1, jpl
+         fice_cell_ave (:,:) = fice_cell_ave (:,:) &
+            &                  + a_i (:,:,jl) * ptab (:,:,jl)
+      END DO
+      
+   END FUNCTION fice_cell_ave
+   
+   
+   FUNCTION fice_ice_ave ( ptab )
+      !!--------------------------------------------------------------------------
+      !! * Compute average over categories, for ice covered part of grid cell
+      !!--------------------------------------------------------------------------
+      REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
+      REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
+
+      fice_ice_ave (:,:) = 0.0_wp
+      WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
+
+   END FUNCTION fice_ice_ave
+
 
 #else
