Index: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4763)
+++ branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5208)
@@ -29,4 +29,5 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! work arrays
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE limthd_ent
@@ -36,7 +37,4 @@
 
    PUBLIC lim_thd_lac     ! called by lim_thd
-
-   REAL(wp) ::   epsi10 = 1.e-10_wp   !
-   REAL(wp) ::   epsi20 = 1.e-20_wp   !
 
    !!----------------------------------------------------------------------
@@ -71,10 +69,10 @@
       !!             - Computation of variation of ice volume and mass
       !!             - Computation of frldb after lateral accretion and 
-      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)      
+      !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)      
       !!------------------------------------------------------------------------
-      INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices
-      INTEGER ::   layer, nbpac     ! local integers 
-      INTEGER ::   ii, ij, iter   !   -       -
-      REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zindb, zinda, zde  ! local scalars
+      INTEGER ::   ji,jj,jk,jl      ! dummy loop indices
+      INTEGER ::   nbpac            ! local integers 
+      INTEGER ::   ii, ij, iter     !   -       -
+      REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zde                         ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
@@ -89,5 +87,5 @@
       REAL(wp) ::   zv_newfra
   
-      INTEGER , POINTER, DIMENSION(:) ::   jcat      ! indexes of categories where new ice grows
+      INTEGER , POINTER, DIMENSION(:) ::   jcat        ! indexes of categories where new ice grows
       REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not
 
@@ -101,14 +99,13 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zda_res     ! residual area in case of excessive heat budget
       REAL(wp), POINTER, DIMENSION(:) ::   zat_i_1d    ! total ice fraction    
-      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_lev   ! total ice fraction for level ice only (type 1)   
-      REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb   ! accretion of frazil ice at the ice bottom
+      REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
       REAL(wp), POINTER, DIMENSION(:) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
 
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_old      ! old volume of ice in category jl
-      REAL(wp), POINTER, DIMENSION(:,:) ::   za_old      ! old area of ice in category jl
-      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d     ! 1-D version of a_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d     ! 1-D version of v_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d    ! 1-D version of oa_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d   ! 1-D version of smv_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_b      ! old volume of ice in category jl
+      REAL(wp), POINTER, DIMENSION(:,:) ::   za_b      ! old area of ice in category jl
+      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d   ! 1-D version of a_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d   ! 1-D version of v_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d  ! 1-D version of oa_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d ! 1-D version of smv_i
 
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d   !: 1-D version of e_i
@@ -119,7 +116,7 @@
       CALL wrk_alloc( jpij, jcat )   ! integer
       CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
-      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d )
-      CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
-      CALL wrk_alloc( jpij,jkmax,jpl, ze_i_1d )
+      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
+      CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
+      CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d )
       CALL wrk_alloc( jpi,jpj, zvrel )
 
@@ -132,6 +129,7 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes
-                  e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i )
+                  rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes
+                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) &
+                      &   / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i, wp )
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac
                END DO
@@ -171,7 +169,6 @@
          zgamafr = 0.03
 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
+         DO jj = 2, jpj
+            DO ji = 2, jpi
                IF ( qlead(ji,jj) < 0._wp ) THEN
                   !-------------
@@ -189,7 +186,7 @@
                   ! Frazil ice velocity
                   !---------------------
-                  zindb = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
-                  zvfrx = zindb * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
-                  zvfry = zindb * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
+                  rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
+                  zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
+                  zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
 
                   !-------------------
@@ -197,9 +194,7 @@
                   !-------------------
                   ! C-grid ice velocity
-                  zindb = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  )
-                  zvgx  = zindb * (  u_ice(ji-1,jj  ) * tmu(ji-1,jj  )    &
-                     &             + u_ice(ji,jj    ) * tmu(ji  ,jj  )  ) * 0.5_wp
-                  zvgy  = zindb * (  v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)    &
-                     &             + v_ice(ji,jj    ) * tmv(ji  ,jj  )  ) * 0.5_wp
+                  rswitch = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  )
+                  zvgx    = rswitch * ( u_ice(ji-1,jj  ) * tmu(ji-1,jj  )  + u_ice(ji,jj) * tmu(ji,jj) ) * 0.5_wp
+                  zvgy    = rswitch * ( v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)  + v_ice(ji,jj) * tmv(ji,jj) ) * 0.5_wp
 
                   !-----------------------------------
@@ -243,4 +238,7 @@
             END DO ! loop on ji ends
          END DO ! loop on jj ends
+      ! 
+      CALL lbc_lnk( zvrel(:,:), 'T', 1. )
+      CALL lbc_lnk( hicol(:,:), 'T', 1. )
 
       ENDIF ! End of computation of frazil ice collection thickness
@@ -255,4 +253,6 @@
       ! This occurs if open water energy budget is negative
       nbpac = 0
+      npac(:) = 0
+      !
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -298,9 +298,8 @@
 
          CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, t_bo_1d   (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac)     , sfx_opw, jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, hicol_1d  (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
 
@@ -315,7 +314,6 @@
          ! Keep old ice areas and volume in memory
          !-----------------------------------------
-         zv_old(:,:) = zv_i_1d(:,:) 
-         za_old(:,:) = za_i_1d(:,:)
-
+         zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 
+         za_b(1:nbpac,:) = za_i_1d(1:nbpac,:)
          !----------------------
          ! Thickness of new ice
@@ -324,5 +322,5 @@
             zh_newice(ji) = hiccrit
          END DO
-         IF( fraz_swi == 1 ) zh_newice(:) = hicol_b(:)
+         IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
 
          !----------------------
@@ -331,5 +329,5 @@
          SELECT CASE ( num_sal )
          CASE ( 1 )                    ! Sice = constant 
-            zs_newice(:) = bulk_sal
+            zs_newice(1:nbpac) = bulk_sal
          CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
             DO ji = 1, nbpac
@@ -339,5 +337,5 @@
             END DO
          CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
-            zs_newice(:) =   2.3
+            zs_newice(1:nbpac) =   2.3
          END SELECT
 
@@ -348,6 +346,6 @@
          DO ji = 1, nbpac
             ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
-            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
-               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_b(ji) - rtt, -epsi10 ) )   &
+            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                             &
+               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) )   &
                &                       - rcp  *         ( ztmelts - rtt )  )
          END DO ! ji
@@ -367,5 +365,5 @@
             zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg]
 
-            zEw           = rcp * ( t_bo_b(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_b [J/kg]
+            zEw           = rcp * ( t_bo_1d(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_1d [J/kg]
                                                                    ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 
                                                                    
@@ -388,6 +386,6 @@
 
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            zinda         = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
-            zfrazb        = zinda * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
+            rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
+            zfrazb        = rswitch * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
             zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
             zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
@@ -438,5 +436,5 @@
          DO ji = 1, nbpac
             jl = jcat(ji)                                                    ! categroy in which new ice is put
-            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) ) )   ! 0 if old ice
+            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
          END DO
 
@@ -444,8 +442,8 @@
             DO ji = 1, nbpac
                jl = jcat(ji)
-               zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
+               rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
                ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                      &
-                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_old(ji,jl) )  &
-                  &        * zinda / MAX( zv_i_1d(ji,jl), epsi20 )
+                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  &
+                  &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )
             END DO
          END DO
@@ -468,19 +466,14 @@
             ! new volumes including lateral/bottom accretion + residual
             DO ji = 1, nbpac
-               zinda          = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
-               zv_newfra      = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
-               za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl)               
+               rswitch        = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
+               zv_newfra      = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
+               za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)               
                zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
-
                ! for remapping
                h_i_old (ji,nlay_i+1) = zv_newfra
                qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
             ENDDO
-
             ! --- Ice enthalpy remapping --- !
-            IF( zv_newfra > 0._wp ) THEN
-               CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
-            ENDIF
-
+            CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
          ENDDO
 
@@ -490,6 +483,6 @@
          DO jl = 1, jpl
             DO ji = 1, nbpac
-               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes
-               zoa_i_1d(ji,jl)  = za_old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb   
+               rswitch          = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes
+               zoa_i_1d(ji,jl)  = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch   
             END DO 
          END DO   
@@ -500,5 +493,5 @@
          DO jl = 1, jpl
             DO ji = 1, nbpac
-               zdv   = zv_i_1d(ji,jl) - zv_old(ji,jl)
+               zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl)
                zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
             END DO
@@ -519,5 +512,4 @@
          CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj )
          CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )
-         CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )
 
          CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj )
@@ -534,5 +526,5 @@
                DO ji = 1, jpi
                   ! heat content in Joules
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) * unit_fac ) 
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac ) 
                END DO
             END DO
@@ -543,7 +535,7 @@
       CALL wrk_dealloc( jpij, jcat )   ! integer
       CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
-      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d )
-      CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
-      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_1d )
+      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
+      CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
+      CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d )
       CALL wrk_dealloc( jpi,jpj, zvrel )
       !
