Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 2715)
@@ -1,6 +1,3 @@
 MODULE limvar
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
-   !!----------------------------------------------------------------------
    !!======================================================================
    !!                       ***  MODULE limvar ***
@@ -32,119 +29,92 @@
    !!                        - ot_i(jpi,jpj)  !average ice age
    !!======================================================================
+   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
+   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!----------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE dom_ice
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_var_agg       : 
+   !!   lim_var_glo2eqv   :
+   !!   lim_var_eqv2glo   :
+   !!   lim_var_salprof   : 
+   !!   lim_var_salprof1d :
+   !!   lim_var_bv        :
+   !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
    USE phycst           ! physical constants (ocean directory) 
    USE sbc_oce          ! Surface boundary condition: ocean fields
-   USE thd_ice
-   USE in_out_manager
-   USE ice
-   USE par_ice
+   USE ice              ! LIM variables
+   USE par_ice          ! LIM parameters
+   USE dom_ice          ! LIM domain
+   USE thd_ice          ! LIM thermodynamics
+   USE wrk_nemo         ! workspace manager
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp         ! MPP library
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_var_agg
-   PUBLIC lim_var_glo2eqv
-   PUBLIC lim_var_eqv2glo
-   PUBLIC lim_var_salprof
-   PUBLIC lim_var_bv
-   PUBLIC lim_var_salprof1d
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      epsi20 = 1e-20   ,  &
-      epsi13 = 1e-13   ,  &
-      zzero  = 0.e0    ,  &
-      zone   = 1.e0
-
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
+   PUBLIC   lim_var_agg          !
+   PUBLIC   lim_var_glo2eqv      !
+   PUBLIC   lim_var_eqv2glo      !
+   PUBLIC   lim_var_salprof      !
+   PUBLIC   lim_var_bv           !
+   PUBLIC   lim_var_salprof1d    !
+
+   REAL(wp) ::   eps20 = 1.e-20_wp   ! module constants
+   REAL(wp) ::   eps16 = 1.e-16_wp   !    -       -
+   REAL(wp) ::   eps13 = 1.e-13_wp   !    -       -
+   REAL(wp) ::   eps10 = 1.e-10_wp   !    -       -
+   REAL(wp) ::   eps06 = 1.e-06_wp   !    -       -
+   REAL(wp) ::   zzero = 0.e0        !    -       -
+   REAL(wp) ::   zone  = 1.e0        !    -       -
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE lim_var_agg(n)
+   SUBROUTINE lim_var_agg( kn )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_var_agg  ***
-      !! ** Purpose :
-      !!        This routine aggregates ice-thickness-category variables to  
-      !!                                all-ice variables
-      !!        i.e. it turns VGLO into VAGG
+      !!
+      !! ** Purpose :   aggregates ice-thickness-category variables to all-ice variables
+      !!              i.e. it turns VGLO into VAGG
       !! ** Method  :
       !!
-      !! ** Arguments :
-      !!           kideb , kiut : Starting and ending points on which the 
-      !!                         the computation is applied
-      !!
-      !! ** Inputs / Ouputs : (global commons)
       !! ** Arguments : n = 1, at_i vt_i only
       !!                n = 2 everything
       !!
-      !! ** External : 
-      !!
-      !! ** References :
-      !!
-      !! ** History :
-      !!           (01-2006) Martin Vancoppenolle, UCL-ASTR
-      !!
       !! note : you could add an argument when you need only at_i, vt_i
       !!        and when you need everything
       !!------------------------------------------------------------------
-      !! * Arguments
-
-      !! * Local variables
-      INTEGER ::   ji,       &   ! spatial dummy loop index
-         jj,       &   ! spatial dummy loop index
-         jk,       &   ! vertical layering dummy loop index
-         jl            ! ice category dummy loop index
-
-      REAL ::      zeps, epsi16, zinda, epsi06
-
-      INTEGER, INTENT( in ) ::   n     ! describes what is needed
-
-      !!-- End of declarations
-      !!----------------------------------------------------------------------------------------------
-      zeps = 1.0e-13
-      epsi16 = 1.0e-16
-      epsi06 = 1.0e-6
-
-      !------------------
-      ! Zero everything
-      !------------------
-
-      vt_i(:,:)  = 0.0
-      vt_s(:,:)  = 0.0
-      at_i(:,:)  = 0.0 
-      ato_i(:,:) = 1.0 
-
-      IF ( n .GT. 1 ) THEN
-         et_s(:,:)  = 0.0
-         ot_i(:,:)  = 0.0
-         smt_i(:,:) = 0.0
-         et_i(:,:)  = 0.0
-      ENDIF
+      INTEGER, INTENT( in ) ::   kn     ! =1 at_i & vt only ; = what is needed
+      !
+      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      REAL(wp) ::   zinda
+      !!------------------------------------------------------------------
 
       !--------------------
       ! Compute variables
       !--------------------
-
+      vt_i (:,:) = 0._wp
+      vt_s (:,:) = 0._wp
+      at_i (:,:) = 0._wp
+      ato_i(:,:) = 1._wp
+      !
       DO jl = 1, jpl
          DO jj = 1, jpj
             DO ji = 1, jpi
-
+               !
                vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
                vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
                at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
-
+               !
                zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 
-               icethi(ji,jj) = vt_i(ji,jj) / MAX(at_i(ji,jj),epsi16)*zinda  
-               ! ice thickness
+               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , eps16 ) * zinda  ! ice thickness
             END DO
          END DO
@@ -153,120 +123,75 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ato_i(ji,jj) = MAX(1.0 - at_i(ji,jj), 0.0)   ! open water fraction
-         END DO
-      END DO
-
-      IF ( n .GT. 1 ) THEN
-
+            ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   ! open water fraction
+         END DO
+      END DO
+
+      IF( kn > 1 ) THEN
+         et_s (:,:) = 0._wp
+         ot_i (:,:) = 0._wp
+         smt_i(:,:) = 0._wp
+         et_i (:,:) = 0._wp
+         !
          DO jl = 1, jpl
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  et_s(ji,jj)  = et_s(ji,jj)  +     &       ! snow heat content
-                     e_s(ji,jj,1,jl)
+                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content
                   zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - 0.10 ) ) 
-                  smt_i(ji,jj) = smt_i(ji,jj) +     &       ! ice salinity
-                     smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , zeps ) * &
-                     zinda
+                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , eps13 ) * zinda   ! ice salinity
                   zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 
-                  ot_i(ji,jj)  = ot_i(ji,jj)  +     &       ! ice age
-                     oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , zeps ) * &
-                     zinda
-               END DO
-            END DO
-         END DO
-
+                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , eps13 ) * zinda   ! ice age
+               END DO
+            END DO
+         END DO
+         !
          DO jl = 1, jpl
             DO jk = 1, nlay_i
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     et_i(ji,jj) = et_i(ji,jj) + e_i(ji,jj,jk,jl) ! ice heat
-                     ! content
-                  END DO
-               END DO
-            END DO
-         END DO
-
-      ENDIF ! n .GT. 1
-
+               et_i(:,:) = et_i(:,:) + e_i(:,:,jk,jl)       ! ice heat content
+            END DO
+         END DO
+         !
+      ENDIF
+      !
    END SUBROUTINE lim_var_agg
 
-   !==============================================================================
 
    SUBROUTINE lim_var_glo2eqv
       !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_glo2eqv ***'
-      !! ** Purpose :
-      !!        This routine computes equivalent variables as function of    
-      !!                              global variables 
-      !!        i.e. it turns VGLO into VEQV
-      !! ** Method  :
-      !!
-      !! ** Arguments :
-      !!           kideb , kiut : Starting and ending points on which the 
-      !!                         the computation is applied
-      !!
-      !! ** Inputs / Ouputs : 
-      !!
-      !! ** External : 
-      !!
-      !! ** References :
-      !!
-      !! ** History :
-      !!           (01-2006) Martin Vancoppenolle, UCL-ASTR
-      !!
-      !!------------------------------------------------------------------
-
-      !! * Local variables
-      INTEGER ::   ji,       &   ! spatial dummy loop index
-         jj,       &   ! spatial dummy loop index
-         jk,       &   ! vertical layering dummy loop index
-         jl            ! ice category dummy loop index
-
-      REAL :: zq_i, zaaa, zbbb, zccc, zdiscrim, &
-         ztmelts, zindb, zq_s, zfac1, zfac2
-
-      REAL :: zeps, epsi06
-
-      zeps    = 1.0e-10
-      epsi06  = 1.0e-06
-
-      !!-- End of declarations
-      !!------------------------------------------------------------------------------
+      !!                ***  ROUTINE lim_var_glo2eqv ***
+      !!
+      !! ** Purpose :   computes equivalent variables as function of global variables 
+      !!              i.e. it turns VGLO into VEQV
+      !!------------------------------------------------------------------
+      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   !   -      -
+      !!------------------------------------------------------------------
 
       !-------------------------------------------------------
       ! Ice thickness, snow thickness, ice salinity, ice age
       !-------------------------------------------------------
-!CDIR NOVERRCHK
       DO jl = 1, jpl
-!CDIR NOVERRCHK
          DO jj = 1, jpj
-!CDIR NOVERRCHK
             DO ji = 1, jpi
-               zindb          = 1.0-MAX(0.0,SIGN(1.0,- a_i(ji,jj,jl))) !0 if no ice and 1 if yes
-               ht_i(ji,jj,jl) = v_i(ji,jj,jl)   / MAX( a_i(ji,jj,jl) , zeps ) * zindb
-               ht_s(ji,jj,jl) = v_s(ji,jj,jl)   / MAX( a_i(ji,jj,jl) , zeps ) * zindb
-               o_i(ji,jj,jl)  = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , zeps ) * zindb
-            END DO
-         END DO
-      END DO
-
-      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) )THEN
-
-!CDIR NOVERRCHK
+               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) ) )   !0 if no ice and 1 if yes
+               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps10 ) * zindb
+               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps10 ) * zindb
+               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps10 ) * zindb
+            END DO
+         END DO
+      END DO
+
+      IF(  num_sal == 2  .OR.  num_sal == 4  )THEN
          DO jl = 1, jpl
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  zindb          = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes
-                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX(v_i(ji,jj,jl),zeps) * zindb
-               END DO
-            END DO
-         END DO
-
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) ) )   !0 if no ice and 1 if yes
+                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , eps10 ) * zindb
+               END DO
+            END DO
+         END DO
       ENDIF
 
-      ! salinity profile
-      CALL lim_var_salprof
+      CALL lim_var_salprof      ! salinity profile
 
       !-------------------
@@ -281,22 +206,16 @@
 !CDIR NOVERRCHK
                DO ji = 1, jpi
-                  !Energy of melting q(S,T) [J.m-3]
-                  zq_i       = e_i(ji,jj,jk,jl) / area(ji,jj) / &
-                     MAX( v_i(ji,jj,jl) , epsi06 ) * nlay_i 
-                  ! zindb = 0 if no ice and 1 if yes
-                  zindb      = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) )
-                  !convert units ! very important that this line is here
-                  zq_i       = zq_i * unit_fac * zindb
-                  !Ice layer melt temperature
-                  ztmelts    =  -tmut*s_i(ji,jj,jk,jl) + rtt
-                  !Conversion q(S,T) -> T (second order equation)
-                  zaaa       =  cpic
-                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + &
-                     zq_i / rhoic - lfus
+                  !                                                              ! Energy of melting q(S,T) [J.m-3]
+                  zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , eps06 ) * REAL(nlay_i,wp) 
+                  zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes
+                  zq_i    = zq_i * unit_fac * zindb                              !convert units
+                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature
+                  !
+                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
+                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i / rhoic - lfus
                   zccc       =  lfus * (ztmelts-rtt)
-                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) )
-                  t_i(ji,jj,jk,jl) = rtt + zindb *( - zbbb - zdiscrim ) / & 
-                     ( 2.0 *zaaa )
-                  t_i(ji,jj,jk,jl) = MIN( rtt, MAX(173.15, t_i(ji,jj,jk,jl) ) )
+                  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) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rtt < t_i < rtt
                END DO
             END DO
@@ -307,24 +226,17 @@
       ! Snow temperatures
       !--------------------
-      zfac1 = 1. / ( rhosn * cpic )
+      zfac1 = 1._wp / ( rhosn * cpic )
       zfac2 = lfus / cpic  
-!CDIR NOVERRCHK
       DO jl = 1, jpl
-!CDIR NOVERRCHK
          DO jk = 1, nlay_s
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
+            DO jj = 1, jpj
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  zq_s       = e_s(ji,jj,jk,jl) / area(ji,jj) / &
-                     MAX( v_s(ji,jj,jl) , epsi06 ) * nlay_s 
-                  ! zindb = 0 if no ice and 1 if yes
-                  zindb      = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) ) )
-                  !convert units ! very important that this line is here
-                  zq_s       = zq_s * unit_fac * zindb
+                  zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , eps06 ) ) * REAL(nlay_s,wp)
+                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes
+                  zq_s  = zq_s * unit_fac * zindb                                    ! convert units
+                  !
                   t_s(ji,jj,jk,jl) = rtt + zindb * ( - zfac1 * zq_s + zfac2 )
-                  t_s(ji,jj,jk,jl) = MIN( rtt, MAX(173.15, t_s(ji,jj,jk,jl) ) )
-
+                  t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rtt < t_i < rtt
                END DO
             END DO
@@ -335,66 +247,45 @@
       ! Mean temperature
       !-------------------
-      tm_i(:,:) = 0.0
-!CDIR NOVERRCHK
+      tm_i(:,:) = 0._wp
       DO jl = 1, jpl
-!CDIR NOVERRCHK
          DO jk = 1, nlay_i
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  zindb          = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)))
-                  zindb          = zindb*1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl)))
-                  tm_i(ji,jj) = tm_i(ji,jj) + t_i(ji,jj,jk,jl)*v_i(ji,jj,jl) / &
-                     REAL(nlay_i) / MAX( vt_i(ji,jj) , zeps )
-               END DO
-            END DO
-         END DO
-      END DO
-
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , -a_i(ji,jj,jl) ) )  )   &
+                     &  * (  1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) )  )
+                  tm_i(ji,jj) = tm_i(ji,jj) + t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
+                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , eps10 )  )
+               END DO
+            END DO
+         END DO
+      END DO
+      !
    END SUBROUTINE lim_var_glo2eqv
 
-   !===============================================================================
 
    SUBROUTINE lim_var_eqv2glo
       !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_eqv2glo ***'
-      !! ** Purpose :
-      !!        This routine computes global     variables as function of    
-      !!                              equivalent variables
-      !!        i.e. it turns VEQV into VGLO
+      !!                ***  ROUTINE lim_var_eqv2glo ***
+      !!
+      !! ** Purpose :   computes global variables as function of equivalent variables
+      !!                i.e. it turns VEQV into VGLO
       !! ** Method  :
       !!
-      !! ** Arguments :
-      !!
-      !! ** Inputs / Ouputs : (global commons)
-      !!
-      !! ** External : 
-      !!
-      !! ** References :
-      !!
-      !! ** History :
-      !!           (01-2006) Martin Vancoppenolle, UCL-ASTR
-      !!                     Take it easy man
-      !!                     Life is just a simple game, between 
-      !!                     ups / and downs \ :@)
-      !!
-      !!------------------------------------------------------------------
-
+      !! ** History :  (01-2006) Martin Vancoppenolle, UCL-ASTR
+      !!------------------------------------------------------------------
+      !
       v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
       v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
       smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
       oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:)
-
+      !
    END SUBROUTINE lim_var_eqv2glo
 
-   !===============================================================================
 
    SUBROUTINE lim_var_salprof
       !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_salprof ***'
-      !! ** Purpose :
-      !!        This routine computes salinity profile in function of
-      !!        bulk salinity     
+      !!                ***  ROUTINE lim_var_salprof ***
+      !!
+      !! ** Purpose :   computes salinity profile in function of bulk salinity     
       !!
       !! ** Method  : If bulk salinity greater than s_i_1, 
@@ -406,50 +297,23 @@
       !!
       !! ** References : Vancoppenolle et al., 2007 (in preparation)
-      !!
-      !! ** History :
-      !!           (08-2006) Martin Vancoppenolle, UCL-ASTR
-      !!                     Take it easy man
-      !!                     Life is just a simple game, between ups 
-      !!                     / and downs \ :@)
-      !!
-      !!------------------------------------------------------------------
-      !! * Arguments
-
-      !! * Local variables
-      INTEGER ::             &
-         ji            ,     & !: spatial dummy loop index
-         jj            ,     & !: spatial dummy loop index
-         jk            ,     & !: vertical layering dummy loop index
-         jl                    !: ice category dummy loop index
-
-      REAL(wp) ::            &
-         dummy_fac0    ,     & !: dummy factor used in computations
-         dummy_fac1    ,     & !: dummy factor used in computations
-         dummy_fac     ,     & !: dummy factor used in computations
-         zind0         ,     & !: switch, = 1 if sm_i lt s_i_0
-         zind01        ,     & !: switch, = 1 if sm_i between s_i_0 and s_i_1
-         zindbal       ,     & !: switch, = 1, if 2*sm_i gt sss_m
-         zargtemp              !: dummy factor
-
-      REAL(wp), DIMENSION(nlay_i) ::      &
-         zs_zero               !: linear salinity profile for salinities under
-      !: s_i_0
-
-      REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
-         z_slope_s     ,     & !: slope of the salinity profile
-         zalpha                !: weight factor for s between s_i_0 and s_i_1
-
-      !!-- End of declarations
-      !!------------------------------------------------------------------------------
+      !!------------------------------------------------------------------
+      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), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha   ! 3D pointer
+      !!------------------------------------------------------------------
+
+      IF( wrk_in_use( 2, 1,2 ) ) THEN
+         CALL ctl_stop( 'lim_var_salprof: requested workspace arrays unavailable' )   ;   RETURN
+      END IF
+
+      z_slope_s => wrk_3d_1(:,:,1:jpl)   ! slope of the salinity profile
+      zalpha    => wrk_3d_2(:,:,1:jpl)   ! weight factor for s between s_i_0 and s_i_1
 
       !---------------------------------------
       ! Vertically constant, constant in time
       !---------------------------------------
-
-      IF ( num_sal .EQ. 1 ) THEN
-
-         s_i(:,:,:,:) = bulk_sal
-
-      ENDIF
+      IF( num_sal == 1 )   s_i(:,:,:,:) = bulk_sal
 
       !-----------------------------------
@@ -457,45 +321,32 @@
       !-----------------------------------
 
-      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) )THEN
-
+      IF(   num_sal == 2  .OR.   num_sal == 4   ) THEN
+         !
          DO jk = 1, nlay_i
             s_i(:,:,jk,:)  = sm_i(:,:,:)
-         END DO ! jk
-
-         ! Slope of the linear profile zs_zero
-         !-------------------------------------
+         END DO
+         !
+         DO jl = 1, jpl                               ! Slope of the linear profile 
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( 0.01 , ht_i(ji,jj,jl) )
+               END DO
+            END DO
+         END DO
+         !
+         dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 )       ! Weighting factor between zs_zero and zs_inf
+         dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 )
+
+         zalpha(:,:,:) = 0._wp
          DO jl = 1, jpl
             DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z_slope_s(ji,jj,jl) = 2.0 * sm_i(ji,jj,jl) / MAX( 0.01      &
-                     , ht_i(ji,jj,jl) )
-               END DO ! ji
-            END DO ! jj
-         END DO ! jl
-
-         ! Weighting factor between zs_zero and zs_inf
-         !---------------------------------------------
-         dummy_fac0 = 1. / ( ( s_i_0 - s_i_1 ) )
-         dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 )
-
-         zalpha(:,:,:) = 0.0
-
-!CDIR NOVERRCHK
-         DO jl = 1, jpl
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
                DO ji = 1, jpi
                   ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise
                   zind0  = MAX( 0.0   , SIGN( 1.0  , 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.0 - zind0 ) *                                  &
-                     MAX( 0.0   , SIGN( 1.0  , s_i_1 - sm_i(ji,jj,jl) ) ) 
+                  zind01 = ( 1.0 - zind0 ) * MAX( 0.0   , SIGN( 1.0  , s_i_1 - sm_i(ji,jj,jl) ) ) 
                   ! If 2.sm_i GE sss_m then zindbal = 1
-                  zindbal = MAX( 0.0 , SIGN( 1.0 , 2. * sm_i(ji,jj,jl) -      &
-                     sss_m(ji,jj) ) )
-                  zalpha(ji,jj,jl) = zind0  * 1.0                             &
-                     + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + &
-                     dummy_fac1 )
+                  zindbal = MAX( 0.0 , SIGN( 1.0 , 2. * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
+                  zalpha(ji,jj,jl) = zind0  * 1.0 + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )
                   zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1.0 - zindbal )
                END DO
@@ -503,18 +354,13 @@
          END DO
 
-         ! Computation of the profile
-         !----------------------------
-         dummy_fac = 1. / nlay_i
-
+         dummy_fac = 1._wp / nlay_i                   ! Computation of the profile
          DO jl = 1, jpl
             DO jk = 1, nlay_i
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     ! linear profile with 0 at the surface
-                     zs_zero(jk)      = z_slope_s(ji,jj,jl) * ( jk - 1./2. ) * &
-                        ht_i(ji,jj,jl) * dummy_fac
-                     ! weighting the profile
-                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero(jk) +       &
-                        ( 1.0 - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
+                     !                                      ! linear profile with 0 at the surface
+                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * dummy_fac
+                     !                                      ! weighting the profile
+                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
                   END DO ! ji
                END DO ! jj
@@ -527,137 +373,83 @@
       ! Vertically varying salinity profile, constant in time
       !-------------------------------------------------------
-      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
-
-      IF ( num_sal .EQ. 3 ) THEN
-
-         sm_i(:,:,:) = 2.30
-
-!CDIR NOVERRCHK
+
+      IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
+         !
+         sm_i(:,:,:) = 2.30_wp
+         !
          DO jl = 1, jpl
 !CDIR NOVERRCHK
             DO jk = 1, nlay_i
-!CDIR NOVERRCHK
-               DO jj = 1, jpj
-!CDIR NOVERRCHK
-                  DO ji = 1, jpi
-                     zargtemp  = ( jk - 0.5 ) / nlay_i
-                     s_i(ji,jj,jk,jl) =  1.6 - 1.6 * COS( 3.14169265 * &
-                        ( zargtemp**(0.407/           &
-                        ( 0.573 + zargtemp ) ) ) )
-                  END DO ! ji
-               END DO ! jj
-            END DO ! jk
-         END DO ! jl
+               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
+               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
+               s_i(:,:,jk,jl) =  zsal
+            END DO
+         END DO
 
       ENDIF ! num_sal
-
+      !
+      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('lim_var_salprof: failed to release workspace arrays.')
+      !
    END SUBROUTINE lim_var_salprof
 
-   !===============================================================================
 
    SUBROUTINE lim_var_bv
       !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_bv ***'
-      !! ** Purpose :
-      !!        This routine computes mean brine volume (%) in sea ice
+      !!                ***  ROUTINE lim_var_bv ***
+      !!
+      !! ** Purpose :   computes mean brine volume (%) in sea ice
       !!
       !! ** Method  : e = - 0.054 * S (ppt) / T (C)
       !!
-      !! ** Arguments :
-      !!
-      !! ** Inputs / Ouputs : (global commons)
-      !!
-      !! ** External : 
-      !!
-      !! ** References : Vancoppenolle et al., JGR, 2007
-      !!
-      !! ** History :
-      !!           (08-2006) Martin Vancoppenolle, UCL-ASTR
-      !!
-      !!------------------------------------------------------------------
-      !! * Arguments
-
-      !! * Local variables
-      INTEGER ::   ji,       &   ! spatial dummy loop index
-         jj,       &   ! spatial dummy loop index
-         jk,       &   ! vertical layering dummy loop index
-         jl            ! ice category dummy loop index
-
-      REAL :: zbvi,          &   ! brine volume for a single ice category
-         zeps,          &   ! very small value
-         zindb              ! is there ice or not
-
-      !!-- End of declarations
-      !!------------------------------------------------------------------------------
-
-      zeps = 1.0e-13
-      bv_i(:,:) = 0.0
-!CDIR NOVERRCHK
+      !! References : Vancoppenolle et al., JGR, 2007
+      !!------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      REAL(wp) ::   zbvi, zindb      ! local scalars
+      !!------------------------------------------------------------------
+      !
+      bv_i(:,:) = 0._wp
       DO jl = 1, jpl
-!CDIR NOVERRCHK
          DO jk = 1, nlay_i
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  zindb          = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes
-                  zbvi = - zindb * tmut *s_i(ji,jj,jk,jl) /             & 
-                     MIN( t_i(ji,jj,jk,jl) - 273.15 , zeps )         &
-                     * v_i(ji,jj,jl) / REAL(nlay_i)
-                  bv_i(ji,jj) = bv_i(ji,jj) + zbvi  &
-                     / MAX( vt_i(ji,jj) , zeps )
-               END DO
-            END DO
-         END DO
-      END DO
-
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes
+                  zbvi  = - zindb * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - 273.15 , eps13 )   &
+                     &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp)
+                  bv_i(ji,jj) = bv_i(ji,jj) + zbvi  / MAX( vt_i(ji,jj) , eps13 )
+               END DO
+            END DO
+         END DO
+      END DO
+      !
    END SUBROUTINE lim_var_bv
 
-   !===============================================================================
-
-   SUBROUTINE lim_var_salprof1d(kideb,kiut)
+
+   SUBROUTINE lim_var_salprof1d( kideb, kiut )
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE lim_thd_salprof1d  ***
       !!
       !! ** Purpose :   1d computation of the sea ice salinity profile
-      !!                Works with 1d vectors and is used by thermodynamic
-      !!                modules
-      !!
-      !! history :
-      !!   3.0  !  May  2007 M. Vancoppenolle  Original code
+      !!                Works with 1d vectors and is used by thermodynamic modules
       !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) :: &
-         kideb, kiut             ! thickness category index
-
-      INTEGER ::             &
-         ji, jk,             &   ! geographic and layer index 
-         zji, zjj
-
-      REAL(wp) ::            &
-         dummy_fac0,         &   ! dummy factors
-         dummy_fac1,         &
-         dummy_fac2,         &
-         zalpha    ,         &   ! weighting factor
-         zind0     ,         &   ! switches as in limvar
-         zind01    ,         &   ! switch
-         zindbal   ,         &   ! switch if in freshwater area
-         zargtemp
-
-      REAL(wp), DIMENSION(jpij) ::            &
-         z_slope_s
-
-      REAL(wp), DIMENSION(jpij,jkmax) ::      &
-         zs_zero
-      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
+      !
+      INTEGER  ::   ji, jk    ! dummy loop indices
+      INTEGER  ::   zji, zjj  ! local integers
+      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars
+      REAL(wp) ::   zalpha, zind0, zind01, zindbal, zs_zero              !   -      -
+      !
+      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
+      !!---------------------------------------------------------------------
+
+      IF(  wrk_in_use(1, 1)  ) THEN
+         CALL ctl_stop('lim_var_salprof1d : requestead workspace arrays unavailable.')   ;   RETURN
+      END IF
+      ! Set-up pointers to sub-arrays of workspace arrays
+      z_slope_s  =>  wrk_1d_1 (1:jpij)
 
       !---------------------------------------
       ! Vertically constant, constant in time
       !---------------------------------------
-
-      IF ( num_sal .EQ. 1 ) THEN
-
-         s_i_b(:,:) = bulk_sal
-
-      ENDIF
+      IF( num_sal == 1 )   s_i_b(:,:) = bulk_sal
 
       !------------------------------------------------------
@@ -665,19 +457,15 @@
       !------------------------------------------------------
 
-      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
-
-         ! Slope of the linear profile zs_zero
-         !-------------------------------------
-!CDIR NOVERRCHK
-         DO ji = kideb, kiut 
-            z_slope_s(ji) = 2.0 * sm_i_b(ji) / MAX( 0.01      &
-               , ht_i_b(ji) )
-         END DO ! ji
+      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
+         !
+         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
+            z_slope_s(ji) = 2._wp * sm_i_b(ji) / MAX( 0.01 , ht_i_b(ji) )
+         END DO
 
          ! Weighting factor between zs_zero and zs_inf
          !---------------------------------------------
-         dummy_fac0 = 1. / ( ( s_i_0 - s_i_1 ) )
+         dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 )
          dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 )
-         dummy_fac2 = 1. / nlay_i
+         dummy_fac2 = 1._wp / REAL(nlay_i,wp)
 
 !CDIR NOVERRCHK
@@ -685,26 +473,18 @@
 !CDIR NOVERRCHK
             DO ji = kideb, kiut
-               zji    =  MOD( npb(ji) - 1, jpi ) + 1
-               zjj    =  ( npb(ji) - 1 ) / jpi + 1
-               zalpha = 0.0
+               zji =  MOD( npb(ji) - 1 , jpi ) + 1
+               zjj =     ( npb(ji) - 1 ) / jpi + 1
                ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise
-               zind0  = MAX( 0.0   , SIGN( 1.0  , s_i_0 - sm_i_b(ji) ) ) 
+               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.0 - zind0 ) *                                  &
-                  MAX( 0.0   , SIGN( 1.0  , s_i_1 - sm_i_b(ji) ) ) 
+               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
-               zindbal = MAX( 0.0 , SIGN( 1.0 , 2. * sm_i_b(ji) -      &
-                  sss_m(zji,zjj) ) )
-
-               zalpha = zind0  * 1.0                                       &
-                  + zind01 * ( sm_i_b(ji) * dummy_fac0 +           &
-                  dummy_fac1 )
-               zalpha = zalpha * ( 1.0 - zindbal )
-
-               zs_zero(ji,jk) = z_slope_s(ji) * ( jk - 1./2. ) * &
-                  ht_i_b(ji) * dummy_fac2
+               zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(zji,zjj) ) )
+               !
+               zalpha = (  zind0 + zind01 * ( sm_i_b(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zindbal )
+               !
+               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_b(ji) * dummy_fac2
                ! weighting the profile
-               s_i_b(ji,jk) = zalpha * zs_zero(ji,jk) +       &
-                  ( 1.0 - zalpha ) * sm_i_b(ji)
+               s_i_b(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_b(ji)
             END DO ! ji
          END DO ! jk
@@ -715,31 +495,28 @@
       ! Vertically varying salinity profile, constant in time
       !-------------------------------------------------------
-      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
-
-      IF ( num_sal .EQ. 3 ) THEN
-
-         sm_i_b(:) = 2.30
-
-!CDIR NOVERRCHK
-         DO ji = kideb, kiut
-!CDIR NOVERRCHK
-            DO jk = 1, nlay_i
-               zargtemp  = ( jk - 0.5 ) / nlay_i
-               s_i_b(ji,jk)  =  1.6 - 1.6*cos(3.14169265*(zargtemp**(0.407/ &
-                  (0.573+zargtemp))))
-            END DO ! jk
-         END DO ! ji
-
-      ENDIF ! num_sal
-
+
+      IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
+         !
+         sm_i_b(:) = 2.30_wp
+         !
+!CDIR NOVERRCHK
+         DO jk = 1, nlay_i
+            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
+            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
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+      IF( wrk_not_released(1, 1) )   CALL ctl_stop( 'lim_var_salprof1d : failed to release workspace arrays' )
+      !
    END SUBROUTINE lim_var_salprof1d
 
-   !===============================================================================
-
 #else
-   !!======================================================================
-   !!                       ***  MODULE limvar  ***
-   !!                          no sea ice model
-   !!======================================================================
+   !!----------------------------------------------------------------------
+   !!   Default option         Dummy module          NO  LIM3 sea-ice model
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE lim_var_agg          ! Empty routines
@@ -755,5 +532,6 @@
    SUBROUTINE lim_var_salprof1d    ! Emtpy routines
    END SUBROUTINE lim_var_salprof1d
-
 #endif
+
+   !!======================================================================
 END MODULE limvar
