Index: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4763)
+++ branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 5208)
@@ -50,6 +50,4 @@
    PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2)
 
-   REAL(wp) ::   epsi10 = 1.e-10_wp   !
-      
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -119,5 +117,6 @@
       CHARACTER (len=50) ::   charout
       REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
-      REAL(wp) ::   za, zstms, zsang, zmask   ! local scalars
+      REAL(wp) ::   za, zstms, zmask   ! local scalars
+      REAL(wp) ::   zc1, zc2, zc3             ! ice mass
 
       REAL(wp) ::   dtevp              ! time step for subcycling
@@ -125,10 +124,10 @@
       REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars
       REAL(wp) ::   zu_ice2, zv_ice1   !
-      REAL(wp) ::   zddc, zdtc, zzdst   ! delta on corners and on centre
+      REAL(wp) ::   zddc, zdtc         ! delta on corners and on centre
+      REAL(wp) ::   zdst               ! shear at the center of the grid point
       REAL(wp) ::   zdsshx, zdsshy     ! term for the gradient of ocean surface
       REAL(wp) ::   sigma1, sigma2     ! internal ice stress
 
       REAL(wp) ::   zresm         ! Maximal error on ice velocity
-      REAL(wp) ::   zindb         ! ice (1) or not (0)      
       REAL(wp) ::   zdummy        ! dummy argument
       REAL(wp) ::   zintb, zintn  ! dummy argument
@@ -140,6 +139,4 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points
       REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zc1              ! ice mass
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zusw             ! temporary weight for ice strength computation
       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce1, v_oce1   ! ocean u/v component on U points                           
       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2, v_oce2   ! ocean u/v component on V points
@@ -147,8 +144,6 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
       
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zdd   , zdt      ! Divergence and tension at centre of grid cells
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells
       REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zdst             ! Shear on centre of grid cells
-      REAL(wp), POINTER, DIMENSION(:,:) ::   deltat, deltac   ! Delta at centre and corners of grid cells
       REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2 
       REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12
@@ -160,7 +155,7 @@
 
       CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_alloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                )
-      CALL wrk_alloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds  , zdst  )
-      CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
+      CALL wrk_alloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                )
+      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
+      CALL wrk_alloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
 
 #if  defined key_lim2 && ! defined key_lim2_vp
@@ -179,12 +174,14 @@
       !
       !------------------------------------------------------------------------------!
-      ! 1) Ice-Snow mass (zc1), ice strength (zpresh)                                !
+      ! 1) Ice strength (zpresh)                                !
       !------------------------------------------------------------------------------!
       !
       ! Put every vector to 0
-      zpresh (:,:) = 0._wp   ;   zc1   (:,:) = 0._wp
+      delta_i(:,:) = 0._wp   ;
+      zpresh (:,:) = 0._wp   ;  
       zpreshc(:,:) = 0._wp
       u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp
-      zdd    (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
+      divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
+      shear_i(:,:) = 0._wp
 
 #if defined key_lim3
@@ -196,5 +193,4 @@
 !CDIR NOVERRCHK
          DO ji = 1 , jpi
-            zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
 #if defined key_lim3
             zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj)
@@ -218,10 +214,9 @@
                &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &
                &              tms(ji,jj)     * wght(ji+1,jj+1,1,1)
-            zusw(ji,jj)    = 1.0 / MAX( zstms, epsd )
             zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
                &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
                &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
                &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   &
-               &             ) * zusw(ji,jj)
+               &             ) / MAX( zstms, epsd )
          END DO
       END DO
@@ -265,4 +260,8 @@
          DO ji = fs_2, fs_jpim1
 
+            zc1 = tms(ji  ,jj  ) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
+            zc2 = tms(ji+1,jj  ) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
+            zc3 = tms(ji  ,jj+1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
+
             zt11 = tms(ji  ,jj) * e1t(ji  ,jj)
             zt12 = tms(ji+1,jj) * e1t(ji+1,jj)
@@ -275,6 +274,6 @@
 
             ! Mass, coriolis coeff. and currents
-            zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd)
-            zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd)
+            zmass1(ji,jj) = ( zt12*zc1 + zt11*zc2 ) / (zt11+zt12+epsd)
+            zmass2(ji,jj) = ( zt22*zc1 + zt21*zc3 ) / (zt21+zt22+epsd)
             zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) )   &
                &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd )
@@ -344,5 +343,5 @@
                !  
                !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
-               !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells
+               !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells
                !- zds(:,:): shear on northeast corner of grid cells
                !
@@ -353,5 +352,5 @@
                !                      bugs (Martin, for Miguel).
                !
-               !- ALSO: arrays zdd, zdt, zds and delta could 
+               !- ALSO: arrays zdt, zds and delta could 
                !  be removed in the future to minimise memory demand.
                !
@@ -361,10 +360,10 @@
                !
                !
-               zdd(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)                  &
-                  &          )                                             &
-                  &         / area(ji,jj)
+               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)                  &
+                  &             )                                             &
+                  &            / area(ji,jj)
 
                zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
@@ -408,5 +407,5 @@
 
                !- Calculate Delta at centre of grid cells
-               zzdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          &
+               zdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          &
                   &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          &
                   &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          &
@@ -415,26 +414,11 @@
                   &         / area(ji,jj)
 
-               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 )  
-               ! MV rewriting
-               ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl )
-               !!gm faster to replace the line above with simply:
-               !!                deltat(ji,jj) = MAX( delta, creepl )
-               !!gm end  
-               deltat(ji,jj) = delta + creepl
-               ! END MV
+               delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
+               delta_i(ji,jj) = delta + creepl
                !-Calculate stress tensor components zs1 and zs2 
                !-at centre of grid cells (see section 3.5 of CICE user's guide).
-               !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
-               !   &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
-               !   &          / ( 1._wp + alphaevp * dtotel )
-
-               !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
-               !              zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
-               !   &          / ( 1._wp + alphaevp * ecc2 * dtotel )
-
-               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
-               zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
+               zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) )  &
                   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
-               zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
+               zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) )  &
                   &         / ( 1._wp + dtotel )
 
@@ -468,14 +452,9 @@
                   &        / ( e1f(ji,jj) * e2f(ji,jj) )
 
-               deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
+               zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
 
                !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-               !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
-               !   &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
-               !   &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
-
-               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
                zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
-                  &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+                  &          ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) )  &
                   &          / ( 1.0 + dtotel ) 
 
@@ -513,5 +492,4 @@
                DO ji = fs_2, fs_jpim1
                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
-                  zsang        = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg
                   z0           = zmass1(ji,jj)/dtevp
 
@@ -523,7 +501,7 @@
                      (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj))
                   zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
-                     za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
-                  zcca         = z0+za*cangvg
-                  zccb         = zcorl1(ji,jj)+za*zsang
+                     za*(u_oce1(ji,jj))
+                  zcca         = z0+za
+                  zccb         = zcorl1(ji,jj)
                   u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
 
@@ -536,5 +514,4 @@
 #endif
 #if defined key_bdy
-         ! clem: change u_ice and v_ice at the boundary for each iteration
          CALL bdy_ice_lim_dyn( 'U' )
 #endif         
@@ -546,5 +523,4 @@
 
                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)
-                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
                   z0           = zmass2(ji,jj)/dtevp
                   ! SB modif because ocean has no slip boundary condition
@@ -555,7 +531,7 @@
                      (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
                   zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
-                     za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
-                  zcca         = z0+za*cangvg
-                  zccb         = zcorl2(ji,jj)+za*zsang
+                     za2ct(ji,jj) + za*(v_oce2(ji,jj))
+                  zcca         = z0+za
+                  zccb         = zcorl2(ji,jj)
                   v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
 
@@ -568,5 +544,4 @@
 #endif
 #if defined key_bdy
-         ! clem: change u_ice and v_ice at the boundary for each iteration
          CALL bdy_ice_lim_dyn( 'V' )
 #endif         
@@ -578,5 +553,4 @@
                DO ji = fs_2, fs_jpim1
                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)
-                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
                   z0           = zmass2(ji,jj)/dtevp
                   ! SB modif because ocean has no slip boundary condition
@@ -588,7 +562,7 @@
                      (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
                   zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
-                     za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
-                  zcca         = z0+za*cangvg
-                  zccb         = zcorl2(ji,jj)+za*zsang
+                     za2ct(ji,jj) + za*(v_oce2(ji,jj))
+                  zcca         = z0+za
+                  zccb         = zcorl2(ji,jj)
                   v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
 
@@ -601,5 +575,4 @@
 #endif
 #if defined key_bdy
-         ! clem: change u_ice and v_ice at the boundary for each iteration
          CALL bdy_ice_lim_dyn( 'V' )
 #endif         
@@ -610,11 +583,5 @@
                DO ji = fs_2, fs_jpim1
                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
-                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
                   z0           = zmass1(ji,jj)/dtevp
-                  ! SB modif because ocean has no slip boundary condition
-                  ! GG Bug
-                  !                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)      &
-                  !                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   &
-                  !                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      &
                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
@@ -624,7 +591,7 @@
                      (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj))
                   zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
-                     za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
-                  zcca         = z0+za*cangvg
-                  zccb         = zcorl1(ji,jj)+za*zsang
+                     za*(u_oce1(ji,jj))
+                  zcca         = z0+za
+                  zccb         = zcorl1(ji,jj)
                   u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
                END DO ! ji
@@ -636,5 +603,4 @@
 #endif
 #if defined key_bdy
-         ! clem: change u_ice and v_ice at the boundary for each iteration
          CALL bdy_ice_lim_dyn( 'U' )
 #endif         
@@ -666,6 +632,4 @@
 !CDIR NOVERRCHK
          DO ji = fs_2, fs_jpim1
-            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 
-            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )
             zdummy = vt_i(ji,jj)
             IF ( zdummy .LE. hminrhg ) THEN
@@ -683,5 +647,4 @@
 #endif
 #if defined key_bdy
-      ! clem: change u_ice and v_ice at the boundary
       CALL bdy_ice_lim_dyn( 'U' )
       CALL bdy_ice_lim_dyn( 'V' )
@@ -690,6 +653,4 @@
       DO jj = k_j1+1, k_jpj-1 
          DO ji = fs_2, fs_jpim1
-            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 
-            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )
             zdummy = vt_i(ji,jj)
             IF ( zdummy .LE. hminrhg ) THEN
@@ -713,17 +674,15 @@
 !CDIR NOVERRCHK
          DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
-            !- zdd(:,:), zdt(:,:): divergence and tension at centre 
+            !- divu_i(:,:), zdt(:,:): divergence and tension at centre 
             !- zds(:,:): shear on northeast corner of grid cells
-            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 
-            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )
             zdummy = vt_i(ji,jj)
             IF ( zdummy .LE. hminrhg ) THEN
 
-               zdd(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)                  &
-                  &         )                                              &
-                  &         / area(ji,jj)
+               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)                  &
+                  &            )                                              &
+                  &            / area(ji,jj)
 
                zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
@@ -747,16 +706,11 @@
                   &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
 
-               zdst(ji,jj) = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    &
+               zdst = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    &
                   &           - e2u( ji-1, jj   ) * v_ice1(ji-1,jj  )    &
                   &           + e1v( ji  , jj   ) * u_ice2(ji  ,jj  )    &
                   &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj)
 
-!              deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   & 
-!                  &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 
-!                  &                          ) + creepl
-               ! MV rewriting
-               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 )  
-               deltat(ji,jj) = delta + creepl
-               ! END MV
+               delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
+               delta_i(ji,jj) = delta + creepl
             
             ENDIF ! zdummy
@@ -773,12 +727,10 @@
       DO jj = k_j1+1, k_jpj-1
          DO ji = fs_2, fs_jpim1
-            divu_i (ji,jj) = zdd   (ji,jj)
-            delta_i(ji,jj) = deltat(ji,jj)
             ! begin TECLIM change 
-            zdst(ji,jj)= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &   
+            zdst= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &   
                &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &   
                &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &   
                &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 
-            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst(ji,jj) * zdst(ji,jj) )
+            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
             ! end TECLIM change
          END DO
@@ -834,7 +786,7 @@
       !
       CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_dealloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                )
-      CALL wrk_dealloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds  , zdst  )
-      CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
+      CALL wrk_dealloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                )
+      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
+      CALL wrk_dealloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
 
    END SUBROUTINE lim_rhg
