Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 2715)
@@ -6,4 +6,5 @@
    !! History :   -   ! 2003-05 (M. Vancoppenolle) UCL-ASTR first coding for LIM3-1D
    !!            3.0  ! 2005-12 (M. Vancoppenolle) adapted to the 3-D version
+   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
    !!---------------------------------------------------------------------
 #if defined key_lim3
@@ -16,9 +17,11 @@
    USE phycst           ! physical constants (ocean directory)
    USE sbc_oce          ! Surface boundary condition: ocean fields
-   USE ice              ! LIM: sea-ice variables
-   USE par_ice          ! LIM: sea-ice parameters
-   USE thd_ice          ! LIM: sea-ice thermodynamics
-   USE limvar           ! LIM: sea-ice variables
+   USE ice              ! LIM variables
+   USE par_ice          ! LIM parameters
+   USE thd_ice          ! LIM thermodynamics
+   USE limvar           ! LIM variables
+   USE wrk_nemo         ! workspace manager
    USE in_out_manager   ! I/O manager
+   USE lib_mpp         ! MPP library
 
    IMPLICIT NONE
@@ -29,5 +32,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -51,16 +54,24 @@
       INTEGER  ::   ji, jk     ! dummy loop indices 
       INTEGER  ::   zji, zjj   ! local integers
-      REAL(wp) ::   zsold, zeps, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars
+      REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars
       REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars
-      REAL(wp), DIMENSION(jpij) ::   ze_init, zhiold, zsiold   ! 1D workspace
+      !
+      REAL(wp), POINTER, DIMENSION(:) ::   ze_init, zhiold, zsiold
       !!---------------------------------------------------------------------
 
-      zeps=1.0e-06_wp
+      IF(  wrk_in_use(1, 1,2,3)  ) THEN
+         CALL ctl_stop('lim_thd_dh : requestead workspace arrays unavailable.')   ;   RETURN
+      END IF
+      ! Set-up pointers to sub-arrays of workspace arrays
+      ze_init =>  wrk_1d_1 (1:jpij)
+      zhiold  =>  wrk_1d_2 (1:jpij)
+      zsiold  =>  wrk_1d_3 (1:jpij)
 
       !------------------------------------------------------------------------------|
       ! 1) Constant salinity, constant in time                                       |
       !------------------------------------------------------------------------------|
-
+!!gm comment: if num_sal = 1 s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !!
       IF( num_sal == 1 ) THEN
+         !
          DO jk = 1, nlay_i
             DO ji = kideb, kiut
@@ -79,11 +90,5 @@
       !------------------------------------------------------------------------------|
 
-      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
-
-         !         WRITE(numout,*)
-         !         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
-         !         num_sal
-         !         WRITE(numout,*) '~~~~~~~~~~~'
-         !         WRITE(numout,*)
+      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
 
          !---------------------------------
@@ -91,58 +96,46 @@
          !---------------------------------
          DO ji = kideb, kiut
-            zhiold(ji)   =  ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) -     &
-               dh_i_surf(ji)
-         END DO ! ji
+            zhiold(ji) = ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) - dh_i_surf(ji)
+         END DO
 
          !---------------------
          ! Global heat content
          !---------------------
-
-         ze_init(:)  =  0.0
+         ze_init(:)  =  0._wp
          DO jk = 1, nlay_i
             DO ji = kideb, kiut
                ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
-            END DO ! ji
-         END DO ! jk
-
-         DO ji = kideb, kiut
-
-            !----------
+            END DO
+         END DO
+
+         DO ji = kideb, kiut
+            !
             ! Switches 
             !----------
-
-            ! iflush  : 1 if summer 
-            iflush       =  MAX( 0.0 , SIGN ( 1.0 , t_su_b(ji) - rtt ) ) 
-            ! igravdr : 1 if t_su lt t_bo
-            igravdr      =  MAX( 0.0 , SIGN ( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )
-            ! iaccrbo : 1 if bottom accretion
-            iaccrbo      =  MAX( 0.0 , SIGN ( 1.0 , dh_i_bott(ji) ) )
-            ! isnowic : 1 if snow ice formation
-            i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-2 ) )
-            isnowic      = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch
+            iflush       =         MAX( 0._wp , SIGN( 1.0 , t_su_b(ji) - rtt )        )    ! =1 if summer 
+            igravdr      =         MAX( 0._wp , SIGN( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo
+            iaccrbo      =         MAX( 0._wp , SIGN( 1.0 , dh_i_bott(ji) )           )    ! =1 if bottom accretion
+            i_ice_switch = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + 1.e-2 ) )
+            isnowic      = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - dh_snowice(ji) ) ) * i_ice_switch   ! =1 if snow ice formation
 
             !---------------------
             ! Salinity tendencies
             !---------------------
-
-            ! drainage by gravity drainage
+            !                                   ! drainage by gravity drainage
             dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice 
-
-            ! drainage by flushing  
-            dsm_i_fl_1d(ji)  = - iflush * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
+            !                                   ! drainage by flushing  
+            dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
 
             !-----------------
             ! Update salinity   
             !-----------------
-
             ! only drainage terms ( gravity drainage and flushing )
-            ! snow ice / bottom sources are added in lim_thd_ent
-            ! to conserve energy
+            ! snow ice / bottom sources are added in lim_thd_ent to conserve energy
             zsiold(ji) = sm_i_b(ji)
             sm_i_b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)
 
-            ! if no ice, salinity eq 0.1
+            ! if no ice, salinity = 0.1
             i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
-            sm_i_b(ji)   = i_ice_switch*sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch )
+            sm_i_b(ji)   = i_ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch )
          END DO ! ji
 
@@ -155,4 +148,5 @@
 
          DO ji = kideb, kiut
+!!gm useless
             ! iflush  : 1 if summer 
             iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) ) 
@@ -161,4 +155,5 @@
             ! iaccrbo : 1 if bottom accretion
             iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) )
+!!gm end useless
             !
             fhbri_1d(ji) = 0._wp
@@ -186,8 +181,8 @@
                zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
                zccc         =  lfus * ( ztmelts - rtt )
-               zdiscrim     =  SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) )
+               zdiscrim     =  SQRT(  MAX( zbbb*zbbb - 4.0*zaaa*zccc, 0._wp )  )
                t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa )
-            END DO !ji
-         END DO !jk
+            END DO
+         END DO
          !
       ENDIF ! num_sal .EQ. 2
@@ -197,14 +192,5 @@
       !------------------------------------------------------------------------------|
 
-      IF( num_sal .EQ. 3 ) THEN
-
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
-            num_sal
-         WRITE(numout,*) '~~~~~~~~~~~~'
-
-         CALL lim_var_salprof1d(kideb,kiut)
-
-      ENDIF ! num_sal .EQ. 3
+      IF( num_sal == 3 )   CALL lim_var_salprof1d( kideb, kiut )
 
       !------------------------------------------------------------------------------|
@@ -212,33 +198,22 @@
       !------------------------------------------------------------------------------|
 
-      ! Cox and Weeks, 1974
-      IF (num_sal.eq.5) THEN
-
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
-            num_sal
-         WRITE(numout,*) '~~~~~~~~~~~~'
-
-         DO ji = kideb, kiut
-
+      IF( num_sal == 5 ) THEN      ! Cox and Weeks, 1974
+         !
+         DO ji = kideb, kiut
             zsold = sm_i_b(ji)
-
-            IF (ht_i_b(ji).lt.0.4) THEN
-               sm_i_b(ji)    = 14.24 - 19.39*ht_i_b(ji) 
+            IF( ht_i_b(ji) < 0.4 ) THEN
+               sm_i_b(ji) = 14.24 - 19.39 * ht_i_b(ji) 
             ELSE
-               sm_i_b(ji)    =  7.88 - 1.59*ht_i_b(ji)
-               sm_i_b(ji)    = MIN(sm_i_b(ji),zsold)  
+               sm_i_b(ji) =  7.88 - 1.59 * ht_i_b(ji)
+               sm_i_b(ji) = MIN( sm_i_b(ji) , zsold )  
             ENDIF
-
-            IF ( ht_i_b(ji) .GT. 3.06918239 ) THEN 
-               sm_i_b(ji)     = 3.0
+            IF( ht_i_b(ji) > 3.06918239 ) THEN 
+               sm_i_b(ji) = 3._wp
             ENDIF
-
             DO jk = 1, nlay_i
                s_i_b(ji,jk)   = sm_i_b(ji)
             END DO
-
-         END DO ! ji
-
+         END DO
+         !
       ENDIF ! num_sal
 
@@ -247,8 +222,8 @@
       !------------------------------------------------------------------------------|
 
-      IF ( num_sal .EQ. 4 ) THEN
-         DO ji = kideb, kiut
-            zji                 = MOD( npb(ji) - 1, jpi ) + 1
-            zjj                 = ( npb(ji) - 1 ) / jpi + 1
+      IF ( num_sal == 4 ) THEN
+         DO ji = kideb, kiut
+            zji = MOD( npb(ji) - 1 , jpi ) + 1
+            zjj =    ( npb(ji) - 1 ) / jpi + 1
             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal    )               &
                &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
@@ -256,10 +231,12 @@
       ELSE
          DO ji = kideb, kiut
-            zji                 = MOD( npb(ji) - 1, jpi ) + 1
-            zjj                 = ( npb(ji) - 1 ) / jpi + 1
+            zji = MOD( npb(ji) - 1 , jpi ) + 1
+            zjj =    ( npb(ji) - 1 ) / jpi + 1
             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - s_i_new(ji) )               &
                &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
-         END DO ! ji
+         END DO
       ENDIF
+      !
+      IF( wrk_not_released(1, 1,2,3) )   CALL ctl_stop( 'lim_thd_lac : failed to release workspace arrays' )
       !
    END SUBROUTINE lim_thd_sal
