Index: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4763)
+++ branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5208)
@@ -66,6 +66,4 @@
    PUBLIC   lim_var_salprof1d    !
 
-   REAL(wp) ::   epsi10 = 1.e-10_wp   !    -       -
-
    !!----------------------------------------------------------------------
    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
@@ -92,5 +90,4 @@
       !
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zinda, zindb
       !!------------------------------------------------------------------
 
@@ -111,6 +108,6 @@
                at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
                !
-               zinda = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
-               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda  ! ice thickness
+               rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
+               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch  ! ice thickness
             END DO
          END DO
@@ -132,9 +129,9 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zinda = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
-                  zindb = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
                   et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content
-                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda   ! ice salinity
-                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi10 ) * zindb   ! ice age
+                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
+                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch   ! ice salinity
+                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
+                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi10 ) * rswitch   ! ice age
                END DO
             END DO
@@ -161,5 +158,5 @@
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       REAL(wp) ::   zq_i, zaaa, zbbb, zccc, zdiscrim     ! local scalars
-      REAL(wp) ::   ztmelts, zindb, zq_s, zfac1, zfac2   !   -      -
+      REAL(wp) ::   ztmelts, zq_s, zfac1, zfac2   !   -      -
       !!------------------------------------------------------------------
 
@@ -170,8 +167,8 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes
-               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
-               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
-               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
+               rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes
+               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
+               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
+               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
             END DO
          END DO
@@ -182,6 +179,6 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes
-                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * zindb
+                  rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes
+                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * rswitch
                END DO
             END DO
@@ -203,6 +200,6 @@
                DO ji = 1, jpi
                   !                                                              ! Energy of melting q(S,T) [J.m-3]
-                  zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes
-                  zq_i    = zindb * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 
+                  rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes
+                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 
                   zq_i    = zq_i * unit_fac                             !convert units
                   ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature
@@ -212,5 +209,5 @@
                   zccc       =  lfus * (ztmelts-rtt)
                   zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) )
-                  t_i(ji,jj,jk,jl) = rtt + zindb *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )
+                  t_i(ji,jj,jk,jl) = rtt + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )
                   t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rtt < t_i < rtt
                END DO
@@ -229,9 +226,9 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes
-                  zq_s  = zindb * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)
+                  rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes
+                  zq_s  = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)
                   zq_s  = zq_s * unit_fac                                    ! convert units
                   !
-                  t_s(ji,jj,jk,jl) = rtt + zindb * ( - zfac1 * zq_s + zfac2 )
+                  t_s(ji,jj,jk,jl) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 )
                   t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rtt < t_i < rtt
                END DO
@@ -248,6 +245,6 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
-                  tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
+                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
+                  tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
                      &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  )
                END DO
@@ -295,5 +292,5 @@
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
       REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac, zsal      ! local scalar
-      REAL(wp) ::   zind0, zind01, zindbal, zargtemp , zs_zero   !   -      -
+      REAL(wp) ::   zswi0, zswi01, zswibal, zargtemp , zs_zero   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha   ! 3D pointer
       !!------------------------------------------------------------------
@@ -330,13 +327,13 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise
-                  zind0  = MAX( 0._wp   , SIGN( 1._wp  , s_i_0 - sm_i(ji,jj,jl) ) ) 
-                  ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws 
-                  zind01 = ( 1._wp - zind0 ) * MAX( 0._wp   , SIGN( 1._wp  , s_i_1 - sm_i(ji,jj,jl) ) ) 
-                  ! If 2.sm_i GE sss_m then zindbal = 1
+                  ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise
+                  zswi0  = MAX( 0._wp   , SIGN( 1._wp  , s_i_0 - sm_i(ji,jj,jl) ) ) 
+                  ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws 
+                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , s_i_1 - sm_i(ji,jj,jl) ) ) 
+                  ! If 2.sm_i GE sss_m then zswibal = 1
                   ! this is to force a constant salinity profile in the Baltic Sea
-                  zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
-                  zalpha(ji,jj,jl) = zind0  + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )
-                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zindbal )
+                  zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
+                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )
+                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal )
                END DO
             END DO
@@ -390,5 +387,4 @@
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zindb   !   -      -
       !!------------------------------------------------------------------
 
@@ -399,6 +395,6 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
-                  tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
+                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
+                  tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
                      &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  )
                END DO
@@ -421,5 +417,5 @@
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zbvi, zinda, zindb      ! local scalars
+      REAL(wp) ::   zbvi             ! local scalars
       !!------------------------------------------------------------------
       !
@@ -429,9 +425,9 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zinda = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  )
-                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
-                  zbvi  = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   &
+                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  )
+                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   &
                      &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp)
-                  bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi  / MAX( vt_i(ji,jj) , epsi10 )
+                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
+                  bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi10 )
                END DO
             END DO
@@ -454,5 +450,5 @@
       INTEGER  ::   ii, ij  ! local integers
       REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars
-      REAL(wp) ::   zalpha, zind0, zind01, zindbal, zs_zero              !   -      -
+      REAL(wp) ::   zalpha, zswi0, zswi01, zswibal, zs_zero              !   -      -
       !
       REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
@@ -464,5 +460,5 @@
       ! Vertically constant, constant in time
       !---------------------------------------
-      IF( num_sal == 1 )   s_i_b(:,:) = bulk_sal
+      IF( num_sal == 1 )   s_i_1d(:,:) = bulk_sal
 
       !------------------------------------------------------
@@ -473,5 +469,5 @@
          !
          DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
-            z_slope_s(ji) = 2._wp * sm_i_b(ji) / MAX( epsi10 , ht_i_b(ji) )
+            z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) )
          END DO
 
@@ -488,17 +484,17 @@
                ii =  MOD( npb(ji) - 1 , jpi ) + 1
                ij =     ( npb(ji) - 1 ) / jpi + 1
-               ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise
-               zind0  = MAX( 0._wp , SIGN( 1._wp  , s_i_0 - sm_i_b(ji) ) ) 
-               ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws 
-               zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) ) 
-               ! if 2.sm_i GE sss_m then zindbal = 1
+               ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise
+               zswi0  = MAX( 0._wp , SIGN( 1._wp  , s_i_0 - sm_i_1d(ji) ) ) 
+               ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws 
+               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_1d(ji) ) ) 
+               ! if 2.sm_i GE sss_m then zswibal = 1
                ! this is to force a constant salinity profile in the Baltic Sea
-               zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) )
+               zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
                !
-               zalpha = (  zind0 + zind01 * ( sm_i_b(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zindbal )
+               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zswibal )
                !
-               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_b(ji) * dummy_fac2
+               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2
                ! weighting the profile
-               s_i_b(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_b(ji)
+               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
             END DO ! ji
          END DO ! jk
@@ -512,5 +508,5 @@
       IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
          !
-         sm_i_b(:) = 2.30_wp
+         sm_i_1d(:) = 2.30_wp
          !
 !CDIR NOVERRCHK
@@ -519,5 +515,5 @@
             zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
             DO ji = kideb, kiut
-               s_i_b(ji,jk) = zsal
+               s_i_1d(ji,jk) = zsal
             END DO
          END DO
