Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8512)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8515)
@@ -52,5 +52,5 @@
 CONTAINS
 
-   SUBROUTINE ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
+   SUBROUTINE ice_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i )
       !!-------------------------------------------------------------------
       !!                 ***  SUBROUTINE ice_rhg_evp  ***
@@ -108,6 +108,6 @@
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt      ! time step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   stress1_i, stress2_i, stress12_i
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   u_ice, v_ice, shear_i, divu_i, delta_i 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i 
       !!
       INTEGER ::   ji, jj       ! dummy loop indices
@@ -247,7 +247,7 @@
 
       ! Initialise stress tensor 
-      zs1 (:,:) = stress1_i (:,:) 
-      zs2 (:,:) = stress2_i (:,:)
-      zs12(:,:) = stress12_i(:,:)
+      zs1 (:,:) = pstress1_i (:,:) 
+      zs2 (:,:) = pstress2_i (:,:)
+      zs12(:,:) = pstress12_i(:,:)
 
       ! Ice strength
@@ -340,6 +340,6 @@
          IF(ln_ctl) THEN   ! Convergence test
             DO jj = 1, jpjm1
-               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
-               zv_ice(:,jj) = v_ice(:,jj)
+               zu_ice(:,jj) = pu_ice(:,jj) ! velocity at previous time step
+               zv_ice(:,jj) = pv_ice(:,jj)
             END DO
          ENDIF
@@ -350,6 +350,6 @@
 
                ! 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)   &
+               zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+                  &         + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
                   &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
 
@@ -367,12 +367,12 @@
               
                ! 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)   &
+               zdiv  = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj)   &
+                  &    + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1)   &
                   &    ) * r1_e1e2t(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)   &
+               zdt  = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &   - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
                   &   ) * r1_e1e2t(ji,jj)
                zdt2 = zdt * zdt
@@ -426,10 +426,10 @@
                   !
                   !                !--- u_ice at V point
-               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
-                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+               u_iceV(ji,jj) = 0.5_wp * ( ( pu_ice(ji,jj  ) + pu_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
+                  &                     + ( pu_ice(ji,jj+1) + pu_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
                   !
                   !                !--- v_ice at U point
-               v_iceU(ji,jj) = 0.5_wp * ( ( 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) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+               v_iceU(ji,jj) = 0.5_wp * ( ( pv_ice(ji  ,jj) + pv_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
+                  &                     + ( pv_ice(ji+1,jj) + pv_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
                !
             END DO
@@ -444,17 +444,17 @@
                DO ji = fs_2, fs_jpim1
                      !                 !--- tau_io/(v_oce - v_ice)
-                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) )  &
                      &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
                      !                 !--- Ocean-to-Ice stress
-                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) )
                      !
                      !                 !--- tau_bottom/v_ice
-                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
+                  zvel  = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
                   zTauB = - tau_icebfr(ji,jj) / zvel
                      !
                      !                 !--- Coriolis at V-points (energy conserving formulation)
                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
-                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
-                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * pu_ice(ji,jj  ) + e2u(ji-1,jj  ) * pu_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) )
                      !
                      !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
@@ -465,13 +465,13 @@
                      !
                      !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
-                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
-                     &           ) * zmaskV(ji,jj)
+                  pv_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * pv_ice(ji,jj)                            &  ! previous velocity
+                     &                                      + zTauE + zTauO * pv_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                      ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
+                     &              + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin
+                     &            ) * zmaskV(ji,jj)
                      !
                   ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
+                  !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) )                  &
                   !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
@@ -479,5 +479,5 @@
                END DO
             END DO
-            CALL lbc_lnk( v_ice, 'V', -1. )
+            CALL lbc_lnk( pv_ice, 'V', -1. )
             !
 #if defined key_agrif
@@ -491,18 +491,18 @@
                                
                   ! tau_io/(u_oce - u_ice)
-                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) )  &
                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
 
                   ! Ocean-to-Ice stress
-                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) )
 
                   ! tau_bottom/u_ice
-                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
+                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) )
                   zTauB = - tau_icebfr(ji,jj) / zvel
 
                   ! Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
-                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
-                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * pv_ice(ji  ,jj) + e1v(ji  ,jj-1) * pv_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) )
                   
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
@@ -513,17 +513,17 @@
 
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
-                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
-                     &           ) * zmaskU(ji,jj)
+                  pu_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * pu_ice(ji,jj)                            &  ! previous velocity
+                     &                                      + zTauE + zTauO * pu_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                      ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
+                     &              + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin 
+                     &            ) * zmaskU(ji,jj)
                   ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
+                  !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) )                  &
                   !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
                END DO
             END DO
-            CALL lbc_lnk( u_ice, 'U', -1. )
+            CALL lbc_lnk( pu_ice, 'U', -1. )
             !
 #if defined key_agrif
@@ -539,18 +539,18 @@
 
                   ! tau_io/(u_oce - u_ice)
-                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) )  &
                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
 
                   ! Ocean-to-Ice stress
-                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) )
 
                   ! tau_bottom/u_ice
-                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
+                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) )
                   zTauB = - tau_icebfr(ji,jj) / zvel
 
                   ! Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
-                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
-                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * pv_ice(ji  ,jj) + e1v(ji  ,jj-1) * pv_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) )
 
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
@@ -561,17 +561,17 @@
 
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
-                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
-                     &           ) * zmaskU(ji,jj)
+                  pu_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * pu_ice(ji,jj)                            &  ! previous velocity
+                     &                                      + zTauE + zTauO * pu_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                      ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
+                     &              + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin
+                     &            ) * zmaskU(ji,jj)
                   ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
+                  !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) )                  &
                   !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
                END DO
             END DO
-            CALL lbc_lnk( u_ice, 'U', -1. )
+            CALL lbc_lnk( pu_ice, 'U', -1. )
             !
 #if defined key_agrif
@@ -585,18 +585,18 @@
          
                   ! tau_io/(v_oce - v_ice)
-                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) )  &
                      &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
 
                   ! Ocean-to-Ice stress
-                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) )
 
                   ! tau_bottom/v_ice
-                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
+                  zvel  = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
                   ztauB = - tau_icebfr(ji,jj) / zvel
                   
                   ! Coriolis at V-points (energy conserving formulation)
                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
-                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
-                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * pu_ice(ji,jj  ) + e2u(ji-1,jj  ) * pu_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) )
 
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
@@ -607,17 +607,17 @@
                   
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
-                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
-                     &           ) * zmaskV(ji,jj)
+                  pv_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * pv_ice(ji,jj)                            &  ! previous velocity
+                     &                                      + zTauE + zTauO * pv_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                      ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
+                     &              + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin
+                     &            ) * zmaskV(ji,jj)
                   ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
+                  !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) )                  &
                   !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
                END DO
             END DO
-            CALL lbc_lnk( v_ice, 'V', -1. )
+            CALL lbc_lnk( pv_ice, 'V', -1. )
             !
 #if defined key_agrif
@@ -631,5 +631,5 @@
          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) ) )
+               zresr(:,jj) = MAX( ABS( pu_ice(:,jj) - zu_ice(:,jj) ), ABS( pv_ice(:,jj) - zv_ice(:,jj) ) )
             END DO
             zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
@@ -648,6 +648,6 @@
 
             ! 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)   &
+            zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+               &         + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
                &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
 
@@ -660,6 +660,6 @@
             
             ! 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)   &
+            zdt  = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+               &   - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
                &   ) * r1_e1e2t(ji,jj)
             zdt2 = zdt * zdt
@@ -671,24 +671,24 @@
             
             ! shear at T points
-            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
+            pshear_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)   &
-               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
-               &            ) * r1_e1e2t(ji,jj)
+            pdivu_i(ji,jj) = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj)   &
+               &             + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1)   &
+               &             ) * r1_e1e2t(ji,jj)
             
             ! delta at T points
-            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
+            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
             rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
-            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
+            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
 
          END DO
       END DO
-      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
+      CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
       
       ! --- Store the stress tensor for the next time step --- !
-      stress1_i (:,:) = zs1 (:,:)
-      stress2_i (:,:) = zs2 (:,:)
-      stress12_i(:,:) = zs12(:,:)
+      pstress1_i (:,:) = zs1 (:,:)
+      pstress2_i (:,:) = zs2 (:,:)
+      pstress12_i(:,:) = zs12(:,:)
       !
 
@@ -703,6 +703,6 @@
 
       ! --- divergence, shear and strength --- !
-      IF( iom_use('idive')  )   CALL iom_put( "idive"  , divu_i(:,:)   * zswi(:,:) )   ! divergence
-      IF( iom_use('ishear') )   CALL iom_put( "ishear" , shear_i(:,:)  * zswi(:,:) )   ! shear
+      IF( iom_use('idive')  )   CALL iom_put( "idive"  , pdivu_i(:,:)   * zswi(:,:) )   ! divergence
+      IF( iom_use('ishear') )   CALL iom_put( "ishear" , pshear_i(:,:)  * zswi(:,:) )   ! shear
       IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
 
@@ -714,19 +714,19 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
-                  &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  &
-                  &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
-
-               zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
+               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
+                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
+                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
+
+               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
 
                zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
 
-!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
-!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
-!!               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
+!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
+!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
+!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
 !!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
-               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
+               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
                zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                    ! shear stress
-               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
+               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
             END DO
          END DO
@@ -773,6 +773,6 @@
                
                ! 2D ice mass, snow mass, area transport arrays (X, Y)
-               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
-               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
+               zfac_x = 0.5 * pu_ice(ji,jj) * e2u(ji,jj) * rswitch
+               zfac_y = 0.5 * pv_ice(ji,jj) * e1v(ji,jj) * rswitch
                
                zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
