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
