Index: branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5279)
+++ branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5568)
@@ -22,5 +22,4 @@
    USE thd_ice        ! LIM thermodynamics
    USE dom_ice        ! LIM domain
-   USE par_ice        ! LIM parameters
    USE ice            ! LIM variables
    USE limtab         ! LIM 2D <==> 1D
@@ -32,4 +31,5 @@
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE limthd_ent
+   USE limvar
 
    IMPLICIT NONE
@@ -106,5 +106,4 @@
       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
 
@@ -112,4 +111,6 @@
 
       REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity
+
+      REAL(wp) :: zcai = 1.4e-3_wp
       !!-----------------------------------------------------------------------!
 
@@ -117,8 +118,10 @@
       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, 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( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d )
+      CALL wrk_alloc( jpij,nlay_i,jpl, ze_i_1d )
       CALL wrk_alloc( jpi,jpj, zvrel )
 
+      CALL lim_var_agg(1)
+      CALL lim_var_glo2eqv
       !------------------------------------------------------------------------------|
       ! 2) Convert units for ice internal energy
@@ -129,8 +132,6 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  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
+                  rswitch          = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  )   !0 if no ice
+                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp )
                END DO
             END DO
@@ -155,7 +156,7 @@
 
       ! Default new ice thickness 
-      hicol(:,:) = hiccrit
-
-      IF( fraz_swi == 1 ) THEN
+      hicol(:,:) = rn_hnewice
+
+      IF( ln_frazil ) THEN
 
          !--------------------
@@ -166,5 +167,5 @@
          zhicrit = 0.04 ! frazil ice thickness
          ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
-         zsqcd   = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
+         zsqcd   = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag)
          zgamafr = 0.03
 
@@ -176,8 +177,8 @@
                   !-------------
                   ! C-grid wind stress components
-                  ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  )   &
-                     &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) * 0.5_wp
-                  ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1)   &
-                     &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp
+                  ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
+                     &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
+                  ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   &
+                     &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
                   ! Square root of wind stress
                   ztenagm       =  SQRT( SQRT( ztaux**2 + ztauy**2 ) )
@@ -195,6 +196,6 @@
                   ! C-grid ice velocity
                   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
+                  zvgx    = rswitch * ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
+                  zvgy    = rswitch * ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
 
                   !-----------------------------------
@@ -222,5 +223,5 @@
                   iterate_frazil = .true.
 
-                  DO WHILE ( iter .LT. 100 .AND. iterate_frazil ) 
+                  DO WHILE ( iter < 100 .AND. iterate_frazil ) 
                      zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) &
                         - hicol(ji,jj) * zhicrit * ztwogp * zvrel2
@@ -266,7 +267,7 @@
       ! debug point to follow
       jiindex_1d = 0
-      IF( ln_nicep ) THEN
-         DO ji = mi0(jiindx), mi1(jiindx)
-            DO jj = mj0(jjindx), mj1(jjindx)
+      IF( ln_icectl ) THEN
+         DO ji = mi0(iiceprt), mi1(iiceprt)
+            DO jj = mj0(jiceprt), mj1(jiceprt)
                IF ( qlead(ji,jj)  <  0._wp ) THEN
                   jiindex_1d = (jj - 1) * jpi + ji
@@ -276,5 +277,5 @@
       ENDIF
    
-      IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
+      IF( ln_icectl ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
 
       !------------------------------
@@ -290,10 +291,9 @@
             CALL tab_2d_1d( nbpac, za_i_1d  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
             CALL tab_2d_1d( nbpac, zv_i_1d  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
             CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
             DO jk = 1, nlay_i
                CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
-            END DO ! jk
-         END DO ! jl
+            END DO
+         END DO
 
          CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) )
@@ -320,19 +320,19 @@
          !----------------------
          DO ji = 1, nbpac
-            zh_newice(ji) = hiccrit
-         END DO
-         IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
+            zh_newice(ji) = rn_hnewice
+         END DO
+         IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
 
          !----------------------
          ! Salinity of new ice 
          !----------------------
-         SELECT CASE ( num_sal )
+         SELECT CASE ( nn_icesal )
          CASE ( 1 )                    ! Sice = constant 
-            zs_newice(1:nbpac) = bulk_sal
+            zs_newice(1:nbpac) = rn_icesal
          CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
             DO ji = 1, nbpac
                ii =   MOD( npac(ji) - 1 , jpi ) + 1
                ij =      ( npac(ji) - 1 ) / jpi + 1
-               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  )
+               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij)  )
             END DO
          CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
@@ -345,9 +345,9 @@
          ! We assume that new ice is formed at the seawater freezing point
          DO ji = 1, nbpac
-            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
+            ztmelts       = - tmut * zs_newice(ji) + rt0                  ! Melting point (K)
             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
+               &                       + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
+               &                       - rcp  *         ( ztmelts - rt0 )  )
+         END DO
 
          !----------------
@@ -356,5 +356,5 @@
          DO ji = 1, nbpac
             zo_newice(ji) = 0._wp
-         END DO ! ji
+         END DO
 
          !-------------------
@@ -363,7 +363,7 @@
          DO ji = 1, nbpac
 
-            zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg]
-
-            zEw           = rcp * ( t_bo_1d(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_1d [J/kg]
+            zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [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) 
                                                                    
@@ -372,5 +372,5 @@
             zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0) 
                                                                    ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point   
-            zv_newice(ji) = - zfmdt / rhoic
+            zv_newice(ji) = - zfmdt * r1_rhoic
 
             zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux  
@@ -387,5 +387,5 @@
             ! 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 ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
+            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)
@@ -409,6 +409,6 @@
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, nbpac
-            IF ( za_newice(ji) >  ( amax - zat_i_1d(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_1d(ji) )
+            IF ( za_newice(ji) >  ( rn_amax - zat_i_1d(ji) ) ) THEN
+               zda_res(ji)   = za_newice(ji) - ( rn_amax - zat_i_1d(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
                za_newice(ji) = za_newice(ji) - zda_res  (ji)
@@ -459,5 +459,5 @@
             DO jk = 1, nlay_i
                DO ji = 1, nbpac
-                  h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i )
+                  h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i
                   qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
                END DO
@@ -478,14 +478,4 @@
          ENDDO
 
-         !------------
-         ! Update age 
-         !------------
-         DO jl = 1, jpl
-            DO ji = 1, nbpac
-               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   
-
          !-----------------
          ! Update salinity
@@ -504,5 +494,4 @@
             CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj )
             CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj )
-            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj )
             CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj )
             DO jk = 1, nlay_i
@@ -525,6 +514,6 @@
             DO jj = 1, jpj
                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 ,wp ) * unit_fac ) 
+                  ! heat content in J/m2
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
                END DO
             END DO
@@ -536,6 +525,6 @@
       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, 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( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d )
+      CALL wrk_dealloc( jpij,nlay_i,jpl, ze_i_1d )
       CALL wrk_dealloc( jpi,jpj, zvrel )
       !
