Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8514)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8515)
@@ -1,4 +1,4 @@
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-!! LIM3 namelist:  
+!! ESIM namelist:  
 !!              1 - Generic parameters                 (namice_run)
 !!              2 - Ice thickness discretization       (namice_itd)
@@ -35,4 +35,5 @@
 !------------------------------------------------------------------------------
    rn_himean        =   2.0           !  expected domain-average ice thickness (m)
+   rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping
 /
 !------------------------------------------------------------------------------
@@ -46,5 +47,5 @@
       rn_perdg      =   17.0          !     ridging work divided by pot. energy change in ridging
                    ! -- ice_rdgrft -- !
-   rn_cs            =   0.5           !  fraction of shearing energy contributing to ridging
+   rn_csrdg         =   0.5           !  fraction of shearing energy contributing to ridging
               ! -- ice_rdgrft_prep -- !
    ln_partf_lin     = .false.         !  Linear ridging participation function (Thorndike et al, 1975)
@@ -100,5 +101,5 @@
 !------------------------------------------------------------------------------
    ln_icethd        = .true.            !  ice thermo   (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-                     ! -- limthd_dif -- !
+                     ! -- icethd_dif -- !
    rn_kappa_i       =   1.0             !  radiation attenuation coefficient in sea ice (m-1)
    ln_cndi_U64      = .false.           !  sea ice thermal conductivity: k = k0 + beta.S/T            (Untersteiner, 1964)
@@ -107,9 +108,9 @@
    rn_cnd_s         =   0.31            !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
                                         !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
-                      ! -- limthd_dh -- !
+                      ! -- icethd_dh -- !
    ln_icedH         = .true.            !  activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
    rn_blow_s        =   0.66            !     mesure of snow blowing into the leads
                                         !     = 1 => no snow blowing, < 1 => some snow blowing
-                      ! -- limthd_da -- !
+                      ! -- icethd_da -- !
    ln_icedA         = .true.            !  activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
       rn_beta       =   1.0             !     coef. beta for lateral melting param. Recommended range=[0.8-1.2]
@@ -120,13 +121,11 @@
                                         !      => 6  vs 8m = +40% melting at the peak (A~0.5)
                                         !         10 vs 8m = -20% melting
-                     ! -- limthd_lac -- !
+                     ! -- icethd_lac -- !
    ln_icedO         = .true.            !  activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-      rn_hinew      =   0.1             !     thickness for new ice formation in open water (m)
+      rn_hinew      =   0.1             !     thickness for new ice formation in open water (m), must be larger than rn_hnewice
       ln_frazil     = .false.           !     Frazil ice parameterization (ice collection as a function of wind)
          rn_maxfraz =   1.0             !        maximum fraction of frazil ice collecting at the ice base
          rn_vfraz   =   0.417           !        thresold drift speed for frazil ice collecting at the ice bottom (m/s)
          rn_Cfraz   =   5.0             !        squeezing coefficient for frazil ice collecting at the ice bottom
-                    ! -- limitd_th -- !
-   rn_himin         =   0.1             !  minimum ice thickness (m) used in remapping, must be smaller than rn_hnewice
                          ! -- icestp -- !
    nn_iceflx        =  -1               !  Redistribute heat flux over ice categories
@@ -143,5 +142,5 @@
 &namice_sal     !   Ice salinity
 !------------------------------------------------------------------------------
-                    ! -- limthd_sal -- !
+                    ! -- icethd_sal -- !
    ln_icedS       = .true.             !  activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
    nn_icesal      =   2                !  ice salinity option
@@ -171,8 +170,8 @@
 &namice_ini     !   Ice initialization
 !------------------------------------------------------------------------------
-                     ! -- limistate -- !
+                     ! -- iceistate -- !
    ln_iceini      = .true.             !  activate ice initialization (T) or not (F)
    ln_iceini_file = .false.            !  netcdf file provided for initialization (T) or not (F)
-   rn_thres_sst   =   2.0              !  maximum water temperature with initial ice (degC)
+   rn_thres_sst   =   2.0              !  max delta temp. above Tfreeze with initial ice = (sst - tfreeze)
    rn_hts_ini_n   =   0.3              !  initial real snow thickness (m), North
    rn_hts_ini_s   =   0.3              !        "            "             South
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8514)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8515)
@@ -36,5 +36,6 @@
    PRIVATE
 
-   PUBLIC   ice_istate      ! called by icestp.F90
+   PUBLIC   ice_istate        ! called by icestp.F90
+   PUBLIC   ice_istate_init   ! called by icestp.F90
 
    INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
@@ -102,6 +103,6 @@
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'ice_istate : sea-ice initialization '
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~ '
+      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
 
       !--------------------------------------------------------------------
@@ -109,6 +110,4 @@
       !--------------------------------------------------------------------
       !
-      CALL ice_istate_init
-
       ! init surface temperature
       DO jl = 1, jpl
@@ -154,23 +153,19 @@
             !-----------------------------
             ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( ff_t(ji,jj) >= 0._wp ) THEN
-                     zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj)
-                     zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj)
-                     zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj)
-                     zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
-                     zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj)
-                     ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
-                  ELSE
-                     zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj)
-                     zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj)
-                     zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj)
-                     zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
-                     zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj)
-                     ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
-                  ENDIF
-               END DO
-            END DO
+            WHERE( ff_t(:,:) >= 0._wp )
+               zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
+               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
+               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
+               zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
+               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
+               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
+            ELSEWHERE
+               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
+               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
+               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
+               zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
+               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
+               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
+            END WHERE
             !
          ENDIF ! ln_iceini_file
@@ -554,26 +549,24 @@
       slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
       slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
-
-      ! Define the initial parameters
-      ! -------------------------
-
-      IF(lwp) THEN
+      !
+      !
+      IF(lwp) THEN                          ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_ini'
-         WRITE(numout,*) '      initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
-         WRITE(numout,*) '      ice initialization from a netcdf file      ln_iceini_file  = ', ln_iceini_file
-         WRITE(numout,*) '      threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
-         WRITE(numout,*) '      initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
-         WRITE(numout,*) '      initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
-         WRITE(numout,*) '      initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
-         WRITE(numout,*) '      initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
-         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
-         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
-         WRITE(numout,*) '      initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
-         WRITE(numout,*) '      initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
-         WRITE(numout,*) '      initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
-         WRITE(numout,*) '      initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
+         WRITE(numout,*) '      initialization with ice (T) or not (F)                 ln_iceini     = ', ln_iceini
+         WRITE(numout,*) '      ice initialization from a netcdf file                ln_iceini_file  = ', ln_iceini_file
+         WRITE(numout,*) '      max delta ocean temp. above Tfreeze with initial ice   rn_thres_sst  = ', rn_thres_sst
+         WRITE(numout,*) '      initial snow thickness in the north                    rn_hts_ini_n  = ', rn_hts_ini_n
+         WRITE(numout,*) '      initial snow thickness in the south                    rn_hts_ini_s  = ', rn_hts_ini_s 
+         WRITE(numout,*) '      initial ice thickness  in the north                    rn_hti_ini_n  = ', rn_hti_ini_n
+         WRITE(numout,*) '      initial ice thickness  in the south                    rn_hti_ini_s  = ', rn_hti_ini_s
+         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_n  = ', rn_ati_ini_n
+         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_s  = ', rn_ati_ini_s
+         WRITE(numout,*) '      initial  ice salinity  in the north                    rn_smi_ini_n  = ', rn_smi_ini_n
+         WRITE(numout,*) '      initial  ice salinity  in the south                    rn_smi_ini_s  = ', rn_smi_ini_s
+         WRITE(numout,*) '      initial  ice/snw temp  in the north                    rn_tmi_ini_n  = ', rn_tmi_ini_n
+         WRITE(numout,*) '      initial  ice/snw temp  in the south                    rn_tmi_ini_s  = ', rn_tmi_ini_s
       ENDIF
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8514)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8515)
@@ -650,5 +650,5 @@
       REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
       !!
-      NAMELIST/namice_itd/ rn_himean
+      NAMELIST/namice_itd/ rn_himean, rn_himin
       !!------------------------------------------------------------------
       !
@@ -668,4 +668,5 @@
          WRITE(numout,*) '   Namelist namice_itd : '
          WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean
+         WRITE(numout,*) '      minimum ice thickness                          rn_himin  = ', rn_himin 
       ENDIF
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8514)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8515)
@@ -52,5 +52,5 @@
    !
    ! ** namelist (namice_rdgrft) **
-   REAL(wp) ::   rn_cs            ! fraction of shearing energy contributing to ridging            
+   REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging            
    LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
    REAL(wp) ::   rn_gstar         !    fractional area of young ice contributing to ridging
@@ -170,5 +170,5 @@
             !  (thick, newly ridged ice).
 
-            closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
+            closing_net(ji,jj) = rn_csrdg * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
 
             ! 2.2 divu_adv
@@ -867,5 +867,5 @@
       NAMELIST/namice_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
          &                    ln_str_R75, rn_perdg,          &
-         &                    rn_cs     ,                    &
+         &                    rn_csrdg  ,                    &
          &                    ln_partf_lin, rn_gstar,        &
          &                    ln_partf_exp, rn_astar,        & 
@@ -893,5 +893,5 @@
          WRITE(numout,*) '      ice strength parameterization Rothrock (1975)            ln_str_R75   = ', ln_str_R75 
          WRITE(numout,*) '            Ratio of ridging work to PotEner change in ridging rn_perdg     = ', rn_perdg 
-         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_cs        = ', rn_cs 
+         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg 
          WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
          WRITE(numout,*) '            Fraction of ice coverage contributing to ridging   rn_gstar     = ', rn_gstar
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 8514)
+++ /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
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8514)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8515)
@@ -285,4 +285,5 @@
       !                                ! Initial sea-ice state
       IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
+         CALL ice_istate_init
          CALL ice_istate
       ELSE                                    ! start from a restart file
@@ -301,11 +302,7 @@
       tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
       !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( gphit(ji,jj) > 0._wp ) THEN   ;   rn_amax_2d(ji,jj) = rn_amax_n  ! NH
-            ELSE                              ;   rn_amax_2d(ji,jj) = rn_amax_s  ! SH
-            ENDIF
-         END DO
-      END DO
+      WHERE( gphit(:,:) > 0._wp )   ;   rn_amax_2d(:,:) = rn_amax_n  ! NH
+      ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH
+      END WHERE
       !
    END SUBROUTINE ice_init
@@ -381,43 +378,23 @@
       !!----------------------------------------------------------------------
       !
-      DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               a_i_b  (ji,jj,jl)   = a_i  (ji,jj,jl)     ! ice area
-               v_i_b  (ji,jj,jl)   = v_i  (ji,jj,jl)     ! ice volume
-               v_s_b  (ji,jj,jl)   = v_s  (ji,jj,jl)     ! snow volume
-               smv_i_b(ji,jj,jl)   = smv_i(ji,jj,jl)     ! salt content
-               oa_i_b (ji,jj,jl)   = oa_i (ji,jj,jl)     ! areal age content
-               e_s_b  (ji,jj,:,jl) = e_s  (ji,jj,:,jl)   ! snow thermal energy
-               e_i_b  (ji,jj,:,jl) = e_i  (ji,jj,:,jl)   ! ice thermal energy
-               !                                         ! ice thickness
-               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes
-               ht_i_b(ji,jj,jl) = v_i_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
-               ht_s_b(ji,jj,jl) = v_s_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
-            END DO
-         END DO   
-      END DO
-      CALL lbc_lnk_multi(  a_i_b, 'T', 1., v_i_b , 'T', 1., ht_i_b, 'T', 1., smv_i_b, 'T', 1.,   &
-         &                oa_i_b, 'T', 1., v_s_b , 'T', 1., ht_s_b, 'T', 1. )
-      CALL lbc_lnk( e_i_b, 'T', 1. )
-      CALL lbc_lnk( e_s_b, 'T', 1. )
+      a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area
+      v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume
+      v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume
+      smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content
+      oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content
+      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
+      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
+      WHERE( a_i_b(:,:,:) >= epsi20 )
+         ht_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:)   ! ice thickness
+         ht_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:)   ! snw thickness
+      ELSEWHERE
+         ht_i_b(:,:,:) = 0._wp
+         ht_s_b(:,:,:) = 0._wp
+      END WHERE
       
-!!gm Question:  here , a_i_b, u_ice and v_ice  are defined over the whole domain, 
-!!gm            so why not just a copy over the whole domain and no lbc_lnk ?
-!!gm            that is some thing like:
-!            at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 )
-!            u_ice_b(:,:) = u_ice(:,:)
-!            v_ice_b(:,:) = v_ice(:,:)
-!    idem for the loop above
-!!gm
       ! ice velocities & total concentration
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-            at_i_b(ji,jj)  = SUM( a_i_b(ji,jj,:) )
-            u_ice_b(ji,jj) = u_ice(ji,jj)
-            v_ice_b(ji,jj) = v_ice(ji,jj)
-         END DO
-      END DO
-      CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. )
+      at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 )
+      u_ice_b(:,:) = u_ice(:,:)
+      v_ice_b(:,:) = v_ice(:,:)
       !
    END SUBROUTINE store_fields
@@ -433,47 +410,43 @@
       INTEGER  ::   ji, jj      ! dummy loop index
       !!----------------------------------------------------------------------
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            sfx    (ji,jj) = 0._wp   ;
-            sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp
-            sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp
-            sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp
-            sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp
-            sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp
-            !
-            wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp
-            wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp
-            wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp
-            wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp
-            wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp
-            wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp  
-            wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
-            wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
-            wfx_snw_sni(ji,jj) = 0._wp 
-            ! MV MP 2016
-            wfx_pnd(ji,jj) = 0._wp
-            ! END MV MP 2016
-            
-            hfx_thd(ji,jj) = 0._wp   ;
-            hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp
-            hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp
-            hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp
-            hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp
-            hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp
-            hfx_err(ji,jj) = 0._wp   ;   hfx_err_rem(ji,jj) = 0._wp
-            hfx_err_dif(ji,jj) = 0._wp
-            wfx_err_sub(ji,jj) = 0._wp
-            !
-            afx_tot(ji,jj) = 0._wp   ;
-            !
-            diag_heat(ji,jj) = 0._wp ;   diag_smvi(ji,jj) = 0._wp
-            diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp
-            
-            ! SIMIP diagnostics
-            diag_fc_bo(ji,jj)    = 0._wp ; diag_fc_su(ji,jj)    = 0._wp
-            
-            tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
-         END DO
-      END DO
+      sfx    (:,:) = 0._wp   ;
+      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp
+      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
+      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
+      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
+      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
+      !
+      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
+      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
+      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
+      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
+      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
+      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp  
+      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
+      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
+      wfx_snw_sni(:,:) = 0._wp 
+      ! MV MP 2016
+      wfx_pnd(:,:) = 0._wp
+      ! END MV MP 2016
+
+      hfx_thd(:,:) = 0._wp   ;
+      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
+      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
+      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
+      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
+      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
+      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
+      hfx_err_dif(:,:) = 0._wp
+      wfx_err_sub(:,:) = 0._wp
+      !
+      afx_tot(:,:) = 0._wp   ;
+      !
+      diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp
+      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
+
+      ! SIMIP diagnostics
+      diag_fc_bo(:,:)    = 0._wp ; diag_fc_su(:,:)    = 0._wp
+
+      tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
       
    END SUBROUTINE ice_diag0
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8514)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8515)
@@ -539,9 +539,9 @@
       INTEGER  ::   ios   ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_thd/ ln_icethd, rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s,          &
-         &                ln_icedH, rn_blow_s,                                                             &
-         &                ln_icedA, rn_beta, rn_dmin,                                                     &
-         &                ln_icedO, rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz, rn_himin,   &
-         &                nn_iceflx
+      NAMELIST/namice_thd/ ln_icethd, rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s,   &
+         &                 ln_icedH, rn_blow_s,                                                    &
+         &                 ln_icedA, rn_beta, rn_dmin,                                             &
+         &                 ln_icedO, rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz,          &
+         &                 nn_iceflx
       !!-------------------------------------------------------------------
       !
@@ -581,6 +581,4 @@
          WRITE(numout,*) '      Threshold relative drift speed for collection of frazil rn_vfraz     = ', rn_vfraz
          WRITE(numout,*) '      Squeezing coefficient for collection of frazil          rn_Cfraz     = ', rn_Cfraz
-         WRITE(numout,*) '   -- iceitd --'
-         WRITE(numout,*) '      minimum ice thickness                                   rn_himin     = ', rn_himin 
          WRITE(numout,*) '   -- icestp --'
          WRITE(numout,*) '      Multicategory heat flux formulation                     nn_iceflx    = ', nn_iceflx
