Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 6762)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 6763)
@@ -90,5 +90,6 @@
    rn_gamma       =    0.2          !     (ln_landfast=T)  fraction of ocean depth that ice must reach to initiate landfast
                                     !                      recommended range: [0.1 ; 0.25]
-   rn_icebfr      =    1.e-05       !     (ln_landfast=T)  maximum bottom stress per unit area of contact                  
+   rn_icebfr      =    5.6e-06      !     (ln_landfast=T)  maximum bottom stress per unit area of contact                  
+                                    !                      a very large value ensures ice velocity=0 even with a small contact area
                                     !                      recommended range: ??
 /
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6762)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6763)
@@ -336,4 +336,6 @@
    rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity
                            !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds)
+   ln_Cd_L12   = .false.   !  Modify the drag ice-atm and oce-atm depending on ice concentration
+                           !  This parameterization is from Lupkes et al. (JGR 2012)
 /
 !-----------------------------------------------------------------------
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6762)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6763)
@@ -422,4 +422,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                      !: ice temperatures
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b           !: ice velocity
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                     !: ice concentration (total)
             
    !!--------------------------------------------------------------------------
@@ -528,5 +529,5 @@
          &      oa_i_b (jpi,jpj,jpl)                                                    , STAT=ierr(ii) )
       ii = ii + 1
-      ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) )
+      ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) )
       
       ! * Ice thickness distribution variables
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6762)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6763)
@@ -116,5 +116,5 @@
       REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2  ! alpha and beta from Bouillon 2009 and 2013
       REAL(wp) ::   zm1, zm2, zm3                              ! ice/snow mass
-      REAL(wp) ::   delta, zp_delf, zds2
+      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2
       REAL(wp) ::   zTauO, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel  ! temporary scalars
 
@@ -134,5 +134,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2               ! internal stresses
       
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt, zds                  ! tension and shear
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                       ! shear
       REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12            ! stress tensor components
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr     ! check convergence
@@ -149,5 +149,5 @@
       CALL wrk_alloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 )
       CALL wrk_alloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 )
-      CALL wrk_alloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
+      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
       CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2 )
 
@@ -257,5 +257,4 @@
          END DO
       END DO
-
       !
       !------------------------------------------------------------------------------!
@@ -266,6 +265,5 @@
       DO jter = 1 , nn_nevp                           !    loop over jter    !
          !                                            !----------------------!        
-         ! Convergence test
-         IF(ln_ctl) THEN
+         IF(ln_ctl) THEN   ! Convergence test
             DO jj = 1, jpjm1
                zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
@@ -275,16 +273,6 @@
 
          ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               
-               ! divergence at T points
-               divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
-                  &            ) * r1_e12t(ji,jj)
-
-               ! tension at T points
-               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
-                  &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
-                  &         ) * r1_e12t(ji,jj)
+         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
+            DO ji = 1, fs_jpim1
 
                ! shear at F points
@@ -292,4 +280,76 @@
                   &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
                   &         ) * r1_e12f(ji,jj)
+
+            END DO
+         END DO
+         CALL lbc_lnk( zds, 'F', 1. )
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 ! no vector loop
+
+               ! shear**2 at T points (doc eq. A16)
+               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e12t(ji,jj)
+              
+               ! divergence at T points
+               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+                  &    ) * r1_e12t(ji,jj)
+               zdiv2 = zdiv * zdiv
+               
+               ! tension at T points
+               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &   ) * r1_e12t(ji,jj)
+               zdt2 = zdt * zdt
+               
+               ! delta at T points
+               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )  
+
+               ! P/delta at T points
+               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl )
+               
+               ! stress at T points
+               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
+             
+            END DO
+         END DO
+         CALL lbc_lnk( zp_delt, 'T', 1. )
+
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
+               ! P/delta at F points
+               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
+               
+               ! stress at F points
+               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
+
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+ 
+
+         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1               
+
+               ! U points
+               zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
+                  &                    ) * r1_e2u(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
+                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
+                  &                  ) * r1_e12u(ji,jj)
+
+               ! V points
+               zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
+                  &                    ) * r1_e1v(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
+                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
+                  &                  ) * r1_e12v(ji,jj)
 
                ! u_ice at V point
@@ -303,60 +363,4 @@
             END DO
          END DO
-         CALL lbc_lnk( zds, 'F', 1. )
-
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1 ! no vector loop
-
-               ! shear**2 at T points (doc eq. A16)
-               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
-                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
-                  &   ) * 0.25_wp * r1_e12t(ji,jj)
-              
-               ! delta at T points
-               delta          = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zds2 ) * usecc2 )  
-               delta_i(ji,jj) = delta + rn_creepl
-
-               ! P/delta at T points
-               zp_delt(ji,jj) = zpresh(ji,jj) / delta_i(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( zp_delt, 'T', 1. )
-
-         ! --- Stress tensor --- !
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1 ! no vector loop
-               
-               ! P/delta at F points
-               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
-               
-               ! stress tensor at T points
-               zs1(ji,jj) = ( zs1 (ji,jj) * zalph1 + zp_delt(ji,jj) * ( divu_i(ji,jj) - (delta_i(ji,jj)-rn_creepl) ) ) * z1_alph1
-               zs2(ji,jj) = ( zs2 (ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt(ji,jj) * z1_ecc2  )                      ) * z1_alph2
-               ! F points
-               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf        * ( zds(ji,jj) * z1_ecc2  ) * 0.5_wp             ) * z1_alph2
-            END DO
-         END DO
-         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
- 
-         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1               
-               ! U points
-               zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
-                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
-                  &                    ) * r1_e2u(ji,jj)                                                                      &
-                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
-                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
-                  &                  ) * r1_e12u(ji,jj)
-
-               ! V points
-               zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
-                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
-                  &                    ) * r1_e1v(ji,jj)                                                                      &
-                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
-                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
-                  &                  ) * r1_e12v(ji,jj)
-            END DO
-         END DO
          !
          ! --- Computation of ice velocity --- !
@@ -368,10 +372,10 @@
                DO ji = fs_2, fs_jpim1
 
-                  zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
-                     &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) )
-
+                  zvel    = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) )
+                  
                   zTauA   =   za2(ji,jj) * vtau_ice(ji,jj)
-                  zTauO   =   za2(ji,jj) * rhoco * zvel
-                  ztauB   = - tau_icebfr(ji,jj) / zvel
+                  zTauO   =   za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                                 + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) )
+                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
                   zCor    =   zcor2(ji,jj)  * u_ice2(ji,jj)
                   zmt     =   zmass2(ji,jj) * z1_dtevp
@@ -403,10 +407,10 @@
                      &               + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
                   
-                  zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
-                     &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) )
-
+                  zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  
                   zTauA   =   za1(ji,jj) * utau_ice(ji,jj)
-                  zTauO   =   za1(ji,jj) * rhoco * zvel
-                  ztauB   = - tau_icebfr(ji,jj) / zvel
+                  zTauO   =   za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                                 + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) )
+                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
                   zCor    =   zcor1(ji,jj)  * zv_ice1
                   zmt     =   zmass1(ji,jj) * z1_dtevp
@@ -435,10 +439,10 @@
                DO ji = fs_2, fs_jpim1
 
-                  zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
-                     &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) )
+                  zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
 
                   zTauA   =   za1(ji,jj) * utau_ice(ji,jj)
-                  zTauO   =   za1(ji,jj) * rhoco * zvel
-                  ztauB   = - tau_icebfr(ji,jj) / zvel
+                  zTauO   =   za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                                 + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) )
+                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
                   zCor    =   zcor1(ji,jj)  * v_ice1(ji,jj)
                   zmt     =   zmass1(ji,jj) * z1_dtevp
@@ -469,10 +473,10 @@
                      &               + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
 
-                  zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
-                     &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) )
+                  zvel    = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) )
          
                   zTauA   =   za2(ji,jj) * vtau_ice(ji,jj)
-                  zTauO   =   za2(ji,jj) * rhoco * zvel
-                  ztauB   = - tau_icebfr(ji,jj) / zvel
+                  zTauO   =   za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                                 + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) )
+                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
                   zCor    =   zcor2(ji,jj)  * zu_ice2
                   zmt     =   zmass2(ji,jj) * z1_dtevp
@@ -485,5 +489,4 @@
                   !!   &           + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  &
                   !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
-
                END DO
             END DO
@@ -499,6 +502,5 @@
          ENDIF
          
-         ! Convergence test
-         IF(ln_ctl) THEN
+         IF(ln_ctl) THEN   ! Convergence test
             DO jj = 2 , jpjm1
                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
@@ -539,9 +541,33 @@
       ! 5) Recompute delta, shear and div (inputs for mechanical redistribution) 
       !------------------------------------------------------------------------------!
-
-      ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1
+
+            ! shear at F points
+            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               &         ) * r1_e12f(ji,jj)
+
+         END DO
+      END DO           
+      CALL lbc_lnk( zds, 'F', 1. )
+      
       DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            
+         DO ji = 2, jpim1 ! no vector loop
+            
+            ! tension**2 at T points
+            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+               &   ) * r1_e12t(ji,jj)
+            zdt2 = zdt * zdt
+            
+            ! shear**2 at T points (doc eq. A16)
+            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
+               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
+               &   ) * 0.25_wp * r1_e12t(ji,jj)
+            
+            ! shear at T points
+            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
+
             ! divergence at T points
             divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
@@ -549,35 +575,11 @@
                &            ) * r1_e12t(ji,jj)
             
-            ! tension at T points
-            zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
-               &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
-               &         ) * r1_e12t(ji,jj)
-            
-            ! shear at F points
-            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
-               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-               &         ) * r1_e12f(ji,jj)
-            
+            ! delta at T points
+            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )  
+            delta_i(ji,jj) = zdelta + rn_creepl
+
          END DO
       END DO
-      CALL lbc_lnk_multi( divu_i, 'T', 1., zds, 'F', 1. )
-      
-
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1 ! no vector loop
-            
-            ! shear**2 at T points (doc eq. A16)
-            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
-               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
-               &   ) * 0.25_wp * r1_e12t(ji,jj)
-            
-            ! delta at T points
-            delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zds2 ) * usecc2 )  
-            delta_i(ji,jj) = delta + rn_creepl
-            
-            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zds2 )
-         END DO
-      END DO
-      CALL lbc_lnk_multi( delta_i, 'T', 1.,  shear_i, 'T', 1. )
+      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
       
       ! --- Store the stress tensor for the next time step --- !
@@ -626,5 +628,5 @@
       CALL wrk_dealloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 )
       CALL wrk_dealloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 )
-      CALL wrk_dealloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
+      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
       CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2 )
 
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6762)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6763)
@@ -137,10 +137,10 @@
 
       zalb(:,:) = 0._wp
-      WHERE     ( SUM( a_i_b, dim=3 ) <= epsi06 )  ;  zalb(:,:) = 0.066_wp
-      ELSEWHERE                                    ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
+      WHERE     ( at_i_b <= epsi06 )  ;  zalb(:,:) = 0.066_wp
+      ELSEWHERE                       ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
       END WHERE
       IF( iom_use('alb_ice' ) )  CALL iom_put( "alb_ice"  , zalb(:,:) )           ! ice albedo output
 
-      zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - SUM( a_i_b, dim=3 ) )      
+      zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - at_i_b )      
       IF( iom_use('albedo'  ) )  CALL iom_put( "albedo"  , zalb(:,:) )           ! ice albedo output
 
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6762)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6763)
@@ -185,5 +185,5 @@
       CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !  
       CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !  
-      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base 
+      CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base 
       CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice 
       CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 6762)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 6763)
@@ -16,4 +16,5 @@
    !!            3.4  !  2011-11  (C. Harris) Fill arrays required by CICE
    !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk
+   !!            4.0  !  2016-06  (C. Rousset) Add new param of drags with sea-ice (Lupkes at al 2012)
    !!----------------------------------------------------------------------
 
@@ -45,5 +46,5 @@
    USE lib_fortran     ! to use key_nosignedzero
 #if defined key_lim3
-   USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b
+   USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b
    USE limthd_dh       ! for CALL lim_thd_snwblow
 #elif defined key_lim2
@@ -60,5 +61,5 @@
    PUBLIC   blk_ice_core_flx     ! routine called in sbc_ice_lim module
 #endif
-   PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module
+   PUBLIC   turb_core_2z         ! routine called in sbcblk_mfs module
 
    INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read
@@ -76,11 +77,11 @@
 
    !                                             !!! CORE bulk parameters
-   REAL(wp), PARAMETER ::   rhoa =    1.22        ! air density
-   REAL(wp), PARAMETER ::   cpa  = 1000.5         ! specific heat of air
-   REAL(wp), PARAMETER ::   Lv   =    2.5e6       ! latent heat of vaporization
-   REAL(wp), PARAMETER ::   Ls   =    2.839e6     ! latent heat of sublimation
-   REAL(wp), PARAMETER ::   Stef =    5.67e-8     ! Stefan Boltzmann constant
-   REAL(wp), PARAMETER ::   Cice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice
-   REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant
+   REAL(wp), PARAMETER ::   rhoa   =    1.22        ! air density
+   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air
+   REAL(wp), PARAMETER ::   Lv     =    2.5e6       ! latent heat of vaporization
+   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation
+   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant
+   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice
+   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant
 
    !                                  !!* Namelist namsbc_core : CORE bulk parameters
@@ -91,5 +92,9 @@
    REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
-
+   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012)
+
+   !
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem)
+   
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -102,4 +107,15 @@
 CONTAINS
 
+   INTEGER FUNCTION sbc_blk_core_alloc()
+      !!-------------------------------------------------------------------
+      !!             ***  ROUTINE sbc_blk_core_alloc (clem) ***
+      !!-------------------------------------------------------------------
+      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_core_alloc )
+      !
+      IF( lk_mpp                  )   CALL mpp_sum( sbc_blk_core_alloc )
+      IF( sbc_blk_core_alloc /= 0 )   CALL ctl_warn('sbc_blk_core_alloc: failed to allocate arrays')
+   END FUNCTION sbc_blk_core_alloc
+
+   
    SUBROUTINE sbc_blk_core( kt )
       !!---------------------------------------------------------------------
@@ -151,5 +167,5 @@
          &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
          &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
-         &                  sn_tdif, rn_zqt,  rn_zu
+         &                  sn_tdif, rn_zqt,  rn_zu    , ln_Cd_L12
       !!---------------------------------------------------------------------
       !
@@ -157,4 +173,7 @@
       IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
          !                                      ! ====================== !
+         !
+         !                                      ! allocate sbc_blk_core array (clem)
+         IF( sbc_blk_core_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_core : unable to allocate standard arrays' )
          !
          REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters
@@ -319,4 +338,15 @@
          &               Cd, Ch, Ce, zt_zu, zq_zu )
     
+      ! Make ocean-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem)
+#if defined key_lim3
+      IF( ln_Cd_L12 ) THEN
+
+         Cd_oce(:,:) = Cd(:,:)    ! record value of pure ocean-atm. drag
+
+         CALL Cdn10_Lupkes2012( Cd )  ! calculate new drag from Lupkes(2012) equations
+
+      ENDIF
+#endif
+      
       ! ... tau module, i and j component
       DO jj = 1, jpj
@@ -437,14 +467,23 @@
       !!---------------------------------------------------------------------
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2
       REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point
       REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point
+      REAL(wp), DIMENSION(:,:), POINTER ::   Cd               ! transfer coefficient for momentum      (tau)
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau')
       !
-      ! local scalars ( place there for vector optimisation purposes)
-      zcoef_wnorm  = rhoa * Cice
-      zcoef_wnorm2 = rhoa * Cice * 0.5
+      CALL wrk_alloc( jpi,jpj, Cd )
+
+      Cd(:,:) = Cd_ice
+      
+      ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem)
+#if defined key_lim3
+      IF( ln_Cd_L12 ) THEN
+
+         CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations
+
+      ENDIF
+#endif
 
 !!gm brutal....
@@ -467,5 +506,5 @@
                zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   &
                   &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj)
-               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
+               zwnorm_f = rhoa * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
                ! ... ice stress at I-point
                utau_ice(ji,jj) = zwnorm_f * zwndi_f
@@ -476,5 +515,5 @@
                zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   &
                   &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  )
-               wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
+               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
             END DO
          END DO
@@ -493,7 +532,7 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vect. opt.
-               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          &
+               utau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          &
                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
-               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          &
+               vtau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          &
                   &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
             END DO
@@ -510,4 +549,6 @@
       ENDIF
 
+      CALL wrk_dealloc( jpi,jpj, Cd )
+
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau')
       
@@ -548,6 +589,6 @@
       ! local scalars ( place there for vector optimisation purposes)
       zcoef_dqlw   = 4.0 * 0.95 * Stef
-      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8)
-      zcoef_dqsb   = rhoa * cpa * Cice
+      zcoef_dqla   = -Ls * Cd_ice * 11637800. * (-5897.8)
+      zcoef_dqsb   = rhoa * cpa * Cd_ice
 
       zztmp = 1. / ( 1. - albo )
@@ -575,7 +616,7 @@
                ! ... turbulent heat fluxes
                ! Sensible Heat
-               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
+               z_qsb(ji,jj,jl) = rhoa * cpa * Cd_ice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
                ! Latent Heat
-               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                           
+               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cd_ice * wndm_ice(ji,jj)   &                           
                   &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
               ! Latent heat sensitivity for ice (Dqla/Dt)
@@ -820,5 +861,5 @@
          !
       END DO
-
+      
       CALL wrk_dealloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd )
       CALL wrk_dealloc( jpi,jpj, zeta_u, stab )
@@ -903,4 +944,50 @@
    END FUNCTION psi_h
 
+
+   SUBROUTINE Cdn10_Lupkes2012( Cd )
+      !!----------------------------------------------------------------------
+      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
+      !!
+      !! ** Purpose :    Recompute the ice-atm and ocean-atm drags at 10m height to make
+      !!                 them dependent on edges at leads, melt ponds and flows.
+      !!                 After some approximations, this can be resumed to a dependency
+      !!                 on ice concentration.
+      !!                
+      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
+      !!                 with the highest level of approximation: level4, eq.(59)
+      !!                 The drag can be re-written as follows:
+      !!
+      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu
+      !!
+      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59)
+      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
+      !!                    A is the concentration of ice minus melt ponds (if any)
+      !!
+      !!                 This new drag has a parabolic shape (as a function of A) starting at
+      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 
+      !!                 and going down to Cdi(say 1.4e-3) for A=1
+      !!
+      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
+      !!                 => see Lupkes et al (2013)
+      !!
+      !! ** References : Lupkes et al. JGR 2012 (theory)
+      !!                 Lupkes et al. GRL 2013 (application to GCM)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
+      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
+      REAL(wp), PARAMETER ::   znu   = 1._wp
+      REAL(wp), PARAMETER ::   zmu   = 1._wp
+      REAL(wp), PARAMETER ::   zbeta = 1._wp
+      REAL(wp)            ::   zcoef
+      !!----------------------------------------------------------------------
+      zcoef = znu + 1._wp / ( 10._wp * zbeta )
+
+      Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
+         &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
+         &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
+      
+   END SUBROUTINE Cdn10_Lupkes2012
+      
    !!======================================================================
 END MODULE sbcblk_core
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6762)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6763)
@@ -617,4 +617,5 @@
       u_ice_b(:,:)     = u_ice(:,:)
       v_ice_b(:,:)     = v_ice(:,:)
+      at_i_b (:,:)     = SUM( a_i_b(:,:,:), dim=3 )
       
    END SUBROUTINE sbc_lim_bef
