Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 2528)
@@ -1,51 +1,43 @@
 MODULE limthd_sal
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
-   !!----------------------------------------------------------------------
    !!======================================================================
    !!                       ***  MODULE limthd_sal ***
-   !!                   computation salinity variations in
-   !!                               the ice
+   !! LIM-3 sea-ice :  computation of salinity variations in the ice
    !!======================================================================
+   !! History :   -   ! 2003-05 (M. Vancoppenolle) UCL-ASTR first coding for LIM3-1D
+   !!            3.0  ! 2005-12 (M. Vancoppenolle) adapted to the 3-D version
+   !!---------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM-3 sea-ice model
+   !!----------------------------------------------------------------------
    !!   lim_thd_sal : salinity variations in the ice
-   !! * Modules used
+   !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
    USE phycst           ! physical constants (ocean directory)
    USE sbc_oce          ! Surface boundary condition: ocean fields
-   USE thd_ice
-   USE ice
-   USE in_out_manager
-   USE limvar
-   USE par_ice
+   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 in_out_manager   ! I/O manager
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_thd_sal        ! called by lim_thd
-   PUBLIC lim_thd_sal_init   ! called by lim_thd
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      epsi20 = 1e-20   ,  &
-      epsi13 = 1e-13   ,  &
-      zzero  = 0.e0    ,  &
-      zone   = 1.e0
-
-   !!----------------------------------------------------------------------
-   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
+   PUBLIC   lim_thd_sal        ! called by limthd module
+   PUBLIC   lim_thd_sal_init   ! called by iceini module
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE lim_thd_sal(kideb,kiut)
+   SUBROUTINE lim_thd_sal( kideb, kiut )
       !!-------------------------------------------------------------------
-      !!                ***  ROUTINE lim_thd_sal  ***       
-      !! ** Purpose :
-      !!        This routine computes new salinities in the ice
+      !!                ***  ROUTINE lim_thd_sal  ***    
+      !!   
+      !! ** Purpose :   computes new salinities in the ice
       !!
       !! ** Method  :  4 possibilities
@@ -54,67 +46,21 @@
       !!               -> num_sal = 3 -> S = S(z)   [multiyear ice]
       !!               -> num_sal = 4 -> S = S(h)   [Cox and Weeks 74]
-      !!           
-      !! ** Steps
-      !!
-      !! ** Arguments
-      !!
-      !! ** Inputs / Outputs
-      !!
-      !! ** External
-      !!
-      !! ** References
-      !!
-      !! ** History  : 
-      !!
-      !! "Je ne suis pour l'instant qu'a 80% de ma condition, mais c'est 
-      !!  les 30% qui restent qui seront les plus difficiles"
-      !!                                           E. Mpenza
-      !!
-      !!-------------------------------------------------------------------
-      !! History :
-      !!   ori  !  03-05 M. Vancoppenolle UCL-ASTR first coding for LIM-1D
-      !!   3.0  !  05-12 Routine rewritten for the 3-D version
       !!---------------------------------------------------------------------
-      !!
-      !! * Local variables
-      INTEGER, INTENT(in) :: &
-         kideb, kiut             !: thickness category index
-
-      INTEGER ::             &
-         ji, jk ,            &   !: geographic and layer index 
-         zji, zjj
-
-      REAL(wp) ::            &
-         zsold,              &   !: old salinity
-         zeps=1.0e-06   ,    &   !: very small
-         iflush         ,    &   !: flushing (1) or not (0)
-         iaccrbo        ,    &   !: bottom accretion (1) or not (0)
-         igravdr        ,    &   !: gravity drainage or not
-         isnowic        ,    &   !: gravity drainage or not
-         i_ice_switch   ,    &   !: ice thickness above a certain treshold or not
-         ztmelts        ,    &   !: freezing point of sea ice
-         zaaa           ,    &   !: dummy factor
-         zbbb           ,    &   !: dummy factor
-         zccc           ,    &   !: dummy factor
-         zdiscrim                !: dummy factor
-
-      REAL(wp), DIMENSION(jpij) ::          &
-         ze_init        ,    &   !initial total enthalpy
-         zhiold         ,    &
-         zsiold
-
+      INTEGER, INTENT(in) ::  kideb, kiut   ! thickness category index
+      !
+      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) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars
+      REAL(wp), DIMENSION(jpij) ::   ze_init, zhiold, zsiold   ! 1D workspace
       !!---------------------------------------------------------------------
 
+      zeps=1.0e-06_wp
+
       !------------------------------------------------------------------------------|
       ! 1) Constant salinity, constant in time                                       |
       !------------------------------------------------------------------------------|
 
-      IF (num_sal.eq.1) THEN
-
-         !         WRITE(numout,*)
-         !         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
-         !         num_sal
-         !         WRITE(numout,*) '~~~~~~~~~~~~'
-
+      IF( num_sal == 1 ) THEN
          DO jk = 1, nlay_i
             DO ji = kideb, kiut
@@ -122,10 +68,10 @@
             END DO ! ji
          END DO ! jk
-
+         !
          DO ji = kideb, kiut
             sm_i_b(ji)      =  bulk_sal 
          END DO ! ji
-
-      ENDIF ! num_sal .EQ. 1
+         !
+      ENDIF
 
       !------------------------------------------------------------------------------|
@@ -174,6 +120,5 @@
             ! 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
+            isnowic      = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch
 
             !---------------------
@@ -182,12 +127,8 @@
 
             ! drainage by gravity drainage
-            dsm_i_gd_1d(ji) = - igravdr *                                     & 
-               MAX( sm_i_b(ji) - sal_G , 0.0 ) /             &
-               time_G * rdt_ice 
+            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.0 ) /             & 
-               time_F * rdt_ice
+            dsm_i_fl_1d(ji)  = - iflush * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
 
             !-----------------
@@ -199,15 +140,13 @@
             ! 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)  !                 &
+            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
-            i_ice_switch  = 1.0 - MAX ( 0.0, SIGN (1.0 , - ht_i_b(ji) ) )
-            sm_i_b(ji)    = i_ice_switch*sm_i_b(ji) + s_i_min * ( 1.0 -       &
-               i_ice_switch )
+            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 )
          END DO ! ji
 
          ! Salinity profile
-         CALL lim_var_salprof1d(kideb,kiut)
+         CALL lim_var_salprof1d( kideb, kiut )
 
          !----------------------------
@@ -217,11 +156,11 @@
          DO ji = kideb, kiut
             ! iflush  : 1 if summer 
-            iflush       =  MAX( 0.0 , SIGN ( 1.0 , t_su_b(ji) - rtt ) ) 
+            iflush  =  MAX( 0._wp , SIGN ( 1._wp , 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) ) ) 
+            igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) ) 
             ! iaccrbo : 1 if bottom accretion
-            iaccrbo      =  MAX( 0.0 , SIGN ( 1.0 , dh_i_bott(ji) ) )
-
-            fhbri_1d(ji) = 0.0
+            iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) )
+            !
+            fhbri_1d(ji) = 0._wp
          END DO ! ji
 
@@ -230,11 +169,8 @@
          !----------------------------
          DO ji = kideb, kiut
-            i_ice_switch  = 1.0 - MAX ( 0.0, SIGN (1.0 , - ht_i_b(ji) ) )
-            fsbri_1d(ji) = fsbri_1d(ji) -  &
-               i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) *  &
-               ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), &
-               sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice
-            IF ( num_sal .EQ. 4 ) fsbri_1d(ji) = 0.0
-
+            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
+            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         &
+               &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice
+            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp
          END DO ! ji
 
@@ -244,20 +180,15 @@
          !--------------------
          DO jk = 1, nlay_i
-
             DO ji = kideb, kiut
-
                ztmelts    =  -tmut*s_i_b(ji,jk) + rtt
                !Conversion q(S,T) -> T (second order equation)
                zaaa         =  cpic
-               zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + &
-                  q_i_b(ji,jk) / rhoic - lfus
+               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) )
-               t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / & 
-                  ( 2.0 *zaaa )
+               t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa )
             END DO !ji
-
          END DO !jk
-
+         !
       ENDIF ! num_sal .EQ. 2
 
@@ -266,5 +197,5 @@
       !------------------------------------------------------------------------------|
 
-      IF ( num_sal .EQ. 3 ) THEN
+      IF( num_sal .EQ. 3 ) THEN
 
          WRITE(numout,*)
@@ -320,8 +251,6 @@
             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
+            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
          END DO
       ELSE
@@ -329,14 +258,11 @@
             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
+            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
       ENDIF
-
-      !-- End of salinity computations
+      !
    END SUBROUTINE lim_thd_sal
-   !==============================================================================
+
 
    SUBROUTINE lim_thd_sal_init
@@ -346,20 +272,17 @@
       !! ** Purpose :   initialization of ice salinity parameters
       !!
-      !! ** Method  : Read the namicesal namelist and check the parameter
-      !!       values called at the first timestep (nit000)
+      !! ** Method  :   Read the namicesal namelist and check the parameter
+      !!              values called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist namicesal
-      !!
-      !! history :
-      !!   3.0  !  July 2005 M. Vancoppenolle  Original code
       !!-------------------------------------------------------------------
-      NAMELIST/namicesal/ num_sal, bulk_sal, sal_G, time_G, sal_F, time_F, &
-         s_i_max, s_i_min, s_i_0, s_i_1
+      NAMELIST/namicesal/ num_sal, bulk_sal, sal_G, time_G, sal_F, time_F,   &
+         &                s_i_max, s_i_min, s_i_0, s_i_1
       !!-------------------------------------------------------------------
-
-      ! Read Namelist namicesal
-      REWIND ( numnam_ice )
-      READ   ( numnam_ice  , namicesal )
-      IF(lwp) THEN
+      !
+      REWIND( numnam_ice )                   ! Read Namelist namicesal
+      READ  ( numnam_ice  , namicesal )
+      !
+      IF(lwp) THEN                           ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity '
@@ -376,14 +299,11 @@
          WRITE(numout,*) ' 2nd salinity for salinity profile  : ', s_i_1
       ENDIF
-
+      !
    END SUBROUTINE lim_thd_sal_init
 
 #else
    !!----------------------------------------------------------------------
-   !!   Default option         Empty Module                No sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_thd_sal        ! Empty routine
-   END SUBROUTINE lim_thd_sal
+   !!   Default option         Dummy Module          No LIM-3 sea-ice model
+   !!----------------------------------------------------------------------
 #endif
    !!======================================================================
