Index: branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 6439)
+++ branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 6440)
@@ -75,7 +75,6 @@
       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) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
-      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness
       CHARACTER (len = 15) :: fieldid
 
@@ -108,9 +107,9 @@
       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
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity
-
-      REAL(wp) :: zcai = 1.4e-3_wp
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d !: 1-D version of e_i
+
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel     ! relative ice / frazil velocity
+
+      REAL(wp) :: zcai = 1.4e-3_wp                     ! ice-air drag (clem: should be dependent on coupling/forcing used)
       !!-----------------------------------------------------------------------!
 
@@ -143,6 +142,5 @@
       !------------------------------------------------------------------------------!    
       ! hicol is the thickness of new ice formed in open water
-      ! hicol can be either prescribed (frazswi = 0)
-      ! or computed (frazswi = 1)
+      ! hicol can be either prescribed (frazswi = 0) or computed (frazswi = 1)
       ! Frazil ice forms in open water, is transported by wind
       ! accumulates at the edge of the consolidated ice edge
@@ -155,6 +153,8 @@
       zvrel(:,:) = 0._wp
 
-      ! Default new ice thickness 
-      hicol(:,:) = rn_hnewice
+      ! Default new ice thickness
+      WHERE( qlead(:,:) < 0._wp ) ; hicol = rn_hnewice
+      ELSEWHERE                   ; hicol = 0._wp
+      END WHERE
 
       IF( ln_frazil ) THEN
@@ -182,5 +182,5 @@
                      &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
                   ! Square root of wind stress
-                  ztenagm       =  SQRT( SQRT( ztaux**2 + ztauy**2 ) )
+                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
 
                   !---------------------
@@ -205,41 +205,29 @@
                   zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
                      &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 )
-                  zvrel(ji,jj)  = SQRT( zvrel2 )
+                  zvrel(ji,jj) = SQRT( zvrel2 )
 
                   !---------------------
                   ! Iterative procedure
                   !---------------------
-                  hicol(ji,jj) = zhicrit + 0.1 
-                  hicol(ji,jj) = zhicrit +   hicol(ji,jj)    &
-                     &                   / ( hicol(ji,jj) * hicol(ji,jj) -  zhicrit * zhicrit ) * ztwogp * zvrel2
-
-!!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1)
-!!gm                                                   = zhicrit**2 + 0.2*zhicrit +0.01
-!!gm                therefore the 2 lines with hicol can be replaced by 1 line:
-!!gm              hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2
-!!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop
+                  hicol(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
+                     &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
 
                   iter = 1
-                  iterate_frazil = .true.
-
-                  DO WHILE ( iter < 100 .AND. iterate_frazil ) 
-                     zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) &
-                        - hicol(ji,jj) * zhicrit * ztwogp * zvrel2
-                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) &
-                        - zhicrit * ztwogp * zvrel2
-                     zhicol_new = hicol(ji,jj) - zf/zfp
-                     hicol(ji,jj)   = zhicol_new
-
+                  DO WHILE ( iter < 20 ) 
+                     zf  = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) -   &
+                        &    hicol(ji,jj) * zhicrit * ztwogp * zvrel2
+                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
+
+                     hicol(ji,jj) = hicol(ji,jj) - zf/zfp
                      iter = iter + 1
-
-                  END DO ! do while
+                  END DO
 
                ENDIF ! end of selection of pixels where ice forms
 
-            END DO ! loop on ji ends
-         END DO ! loop on jj ends
-      ! 
-      CALL lbc_lnk( zvrel(:,:), 'T', 1. )
-      CALL lbc_lnk( hicol(:,:), 'T', 1. )
+            END DO 
+         END DO 
+         ! 
+         CALL lbc_lnk( zvrel(:,:), 'T', 1. )
+         CALL lbc_lnk( hicol(:,:), 'T', 1. )
 
       ENDIF ! End of computation of frazil ice collection thickness
@@ -282,6 +270,5 @@
       ! Move from 2-D to 1-D vectors
       !------------------------------
-      ! If ocean gains heat do nothing 
-      ! 0therwise compute new ice formation
+      ! If ocean gains heat do nothing. Otherwise compute new ice formation
 
       IF ( nbpac > 0 ) THEN
@@ -297,13 +284,14 @@
          END DO
 
-         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , 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, 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) )
-
-         CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac)     , hfx_thd, jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac)     , hfx_opw, jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead     , 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, 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) )
+
+         CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac)     , hfx_thd   , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac)     , hfx_opw   , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, rn_amax_1d(1:nbpac)     , rn_amax_2d, jpi, jpj, npac(1:nbpac) )
 
          !------------------------------------------------------------------------------!
@@ -316,11 +304,9 @@
          zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 
          za_b(1:nbpac,:) = za_i_1d(1:nbpac,:)
+
          !----------------------
          ! Thickness of new ice
          !----------------------
-         DO ji = 1, nbpac
-            zh_newice(ji) = rn_hnewice
-         END DO
-         IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
+         zh_newice(1:nbpac) = hicol_1d(1:nbpac)
 
          !----------------------
@@ -384,12 +370,17 @@
             ! salt flux
             sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice
-
+         END DO
+         
+         zv_frazb(:) = 0._wp
+         IF( ln_frazil ) THEN
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
-            zfrazb        = rswitch * ( TANH ( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb
-            zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
-            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
-         END DO
-
+            DO ji = 1, nbpac
+               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
+               zfrazb        = rswitch * ( TANH( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb
+               zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
+               zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
+            END DO
+         END IF
+         
          !-----------------
          ! Area of new ice
@@ -409,6 +400,6 @@
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, nbpac
-            IF ( za_newice(ji) >  ( rn_amax - zat_i_1d(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( rn_amax - zat_i_1d(ji) )
+            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN
+               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
                za_newice(ji) = za_newice(ji) - zda_res  (ji)
@@ -443,5 +434,5 @@
                jl = jcat(ji)
                rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
-               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                      &
+               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_b(ji,jl) )  &
                   &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )
