Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/dom_ice_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 1855)
@@ -2,13 +2,11 @@
    !!======================================================================
    !!                   ***  MODULE  dom_ice  ***
-   !! LIM 2.0 Sea Ice :   Domain  variables
+   !! LIM-2 Sea Ice :   Domain  variables
    !!======================================================================
-   !! History :   2.0  !  03-08  (C. Ethe)  Free form and module
+   !! History :   2.0  !  2003-08  (C. Ethe)  Free form and module
    !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
-   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
    USE par_ice_2
@@ -18,17 +16,24 @@
 
    LOGICAL, PUBLIC ::   l_jeq     = .TRUE.     !: Equator inside the domain flag
+   INTEGER, PUBLIC ::   njeq , njeqm1          !: j-index of the equator if it is inside the domain
 
-   INTEGER, PUBLIC ::   njeq , njeqm1          !: j-index of the equator if it is inside the domain
-      !                                        !  (otherwise = jpj+10 (SH) or -10 (SH) )
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor, fcor     !: coriolis factor and coeficient
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   covrai           !: sine of geographic latitude
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   area             !: surface of grid cell 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tms   , tmu      !: temperature and velocity points masks
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght             !: weight of the 4 neighbours to compute averages
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   akappa, bkappa   !: first and third group of metric coefficients
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd           !: second group of metric coefficients
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor,   &  !: coriolis factor and coeficient
-      &                                              covrai ,         &  !: sine of geographic latitude
-      &                                              area   ,         &  !: surface of grid cell 
-      &                                              tms    , tmu        !: temperature and velocity points masks
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght   ,         &  !: weight of the 4 neighbours to compute averages
-      &                                              akappa , bkappa     !: first and third group of metric coefficients
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd   !: second group of metric coefficients
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
 
+   !!-------- -------------------------------------------------------------
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
-#endif
 END MODULE dom_ice_2
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90	(revision 1855)
@@ -4,5 +4,6 @@
    !! Sea Ice physics:  diagnostics variables of ice defined in memory
    !!=====================================================================
-   !! History :  2.0  !  03-08  (C. Ethe)  F90: Free form and module
+   !! History :  2.0  !  2003-08 (C. Ethe)  F90: Free form and module
+   !!             -   !  2010-05 (G. Madec) suppression of thd_ice_2
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -10,9 +11,4 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006)
-   !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_ice_2          ! LIM sea-ice parameters
 
@@ -20,72 +16,126 @@
    PRIVATE
 
-   !!* Share parameters namelist (namicerun read in iceini) *
-   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in  = "restart_ice_in"   !: suffix of ice restart name (input)
-   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output)
-   LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F)
-   LOGICAL               , PUBLIC ::   ln_limdmp     = .FALSE.            !: Ice damping
-   REAL(wp)              , PUBLIC ::   hsndif        = 0.e0               !: computation of temp. in snow (0) or not (9999)
-   REAL(wp)              , PUBLIC ::   hicdif        = 0.e0               !: computation of temp. in ice (0) or not (9999)
-   REAL(wp), DIMENSION(2), PUBLIC ::   acrit = (/ 1.e-06 , 1.e-06 /)      !: minimum fraction for leads in 
-   !                                                                      !: north and south hemisphere
-   !!* ice-dynamic namelist (namicedyn) *
-   INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation
-   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation
-   REAL(wp), PUBLIC ::   rdt_ice            !: ice time step
-   REAL(wp), PUBLIC ::   epsd   = 1.0e-20   !: tolerance parameter for dynamic
-   REAL(wp), PUBLIC ::   alpha  = 0.5       !: coefficient for semi-implicit coriolis
-   REAL(wp), PUBLIC ::   dm     = 0.6e+03   !: diffusion constant for dynamics
-   REAL(wp), PUBLIC ::   om     = 0.5       !: relaxation constant
-   REAL(wp), PUBLIC ::   resl   = 5.0e-05   !: maximum value for the residual of relaxation
-   REAL(wp), PUBLIC ::   cw     = 5.0e-03   !: drag coefficient for oceanic stress
-   REAL(wp), PUBLIC ::   angvg  = 0.e0      !: turning angle for oceanic stress
-   REAL(wp), PUBLIC ::   pstar  = 1.0e+04   !: first bulk-rheology parameter
-   REAL(wp), PUBLIC ::   c_rhg  = 20.e0     !: second bulk-rhelogy parameter
-   REAL(wp), PUBLIC ::   etamn  = 0.e+07    !: minimun value for viscosity
-   REAL(wp), PUBLIC ::   creepl = 2.e-08    !: creep limit
-   REAL(wp), PUBLIC ::   ecc    = 2.e0      !: eccentricity of the elliptical yield curve
-   REAL(wp), PUBLIC ::   ahi0   = 350.e0    !: sea-ice hor. eddy diffusivity coeff. (m2/s)
+   !---------------------------------------------------!
+   !                 Sea-Ice namelist                  !
+   !---------------------------------------------------!
+   !                                                                !!! **  Namelist namicerun  **
+   CHARACTER(len=32), PUBLIC ::   cn_icerst_in  = "restart_ice_in"   !: suffix of ice restart name (input)
+   CHARACTER(len=32), PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output)
+   LOGICAL          , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F)
+   LOGICAL          , PUBLIC ::   ln_limdmp     = .FALSE.            !: Ice damping
+   REAL(wp)         , PUBLIC ::   hsndif        = 0.e0               !: computation of temp. in snow (0) or not (9999)
+   REAL(wp)         , PUBLIC ::   hicdif        = 0.e0               !: computation of temp. in ice (0) or not (9999)
+   REAL(wp)         , PUBLIC ::   acrit(2) = (/ 1.e-06 , 1.e-06 /)   !: min lead fraction in north and south hemispheres
+   
+   !                                         !!! **  Namelist  namicedyn  **
+   INTEGER , PUBLIC ::   nbiter =   1         !: number of sub-time steps for relaxation
+   INTEGER , PUBLIC ::   nbitdr = 250         !: maximum number of iterations for relaxation
+   REAL(wp), PUBLIC ::   rdt_ice              !: ice time step
+   REAL(wp), PUBLIC ::   epsd   =   1.0e-20   !: tolerance parameter for dynamic
+   REAL(wp), PUBLIC ::   alpha  =   0.5       !: coefficient for semi-implicit coriolis
+   REAL(wp), PUBLIC ::   dm     =   0.6e+03   !: diffusion constant for dynamics
+   REAL(wp), PUBLIC ::   om     =   0.5       !: relaxation constant
+   REAL(wp), PUBLIC ::   resl   =   5.0e-05   !: maximum value for the residual of relaxation
+   REAL(wp), PUBLIC ::   cw     =   5.0e-03   !: drag coefficient for oceanic stress
+   REAL(wp), PUBLIC ::   angvg  =   0.e0      !: turning angle for oceanic stress
+   REAL(wp), PUBLIC ::   pstar  =   1.0e+04   !: first bulk-rheology parameter
+   REAL(wp), PUBLIC ::   c_rhg  =  20.e0      !: second bulk-rhelogy parameter
+   REAL(wp), PUBLIC ::   etamn  =   0.e+07    !: minimun value for viscosity
+   REAL(wp), PUBLIC ::   creepl =   2.e-08    !: creep limit
+   REAL(wp), PUBLIC ::   ecc    =   2.e0      !: eccentricity of the elliptical yield curve
+   REAL(wp), PUBLIC ::   ahi0   = 350.e0      !: sea-ice hor. eddy diffusivity coeff. (m2/s)
 
-   REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc )
-   REAL(wp), PUBLIC ::   rhoco              !: = rau0 * cw
-   REAL(wp), PUBLIC ::   sangvg, cangvg     !: sin and cos of the turning angle for ocean stress
-   REAL(wp), PUBLIC ::   pstarh             !: pstar / 2.0
+   !                                              !!! ** Namelist  namicethd  **
+   REAL(wp), PUBLIC ::   hmelt      = -0.15        !: maximum melting at the bottom
+   REAL(wp), PUBLIC ::   hicmin     =  0.2         !: ice th. corr. to max. ener. in brine pocket
+   REAL(wp), PUBLIC ::   hiclim     =  0.05        !: minimum ice thickness
+   REAL(wp), PUBLIC ::   amax       =  0.999       !: maximum lead fraction
+   REAL(wp), PUBLIC ::   swiqst     =  1.0         !: energy stored in brine pocket (1) or not (0)
+   REAL(wp), PUBLIC ::   sbeta      =  1.0         !: numerical scheme for diffusion in ice 
+   REAL(wp), PUBLIC ::   parlat     =  0.0         !: percent. of energy used for lateral ablation
+   REAL(wp), PUBLIC ::   hakspl     =  0.5         !: slope of distr. for Hakkinen-Mellro's lat. melt
+   REAL(wp), PUBLIC ::   hibspl     =  0.5         !: slope of distribution for Hibler's lat. melt
+   REAL(wp), PUBLIC ::   exld       =  2.0         !: exponent for leads-closure rate
+   REAL(wp), PUBLIC ::   hakdif     =  1.0         !: coefficient for diffusions of ice and snow
+   REAL(wp), PUBLIC ::   thth       =  0.2         !: thick. for comp. of eq. thermal conduct
+   REAL(wp), PUBLIC ::   hnzst      =  0.1         !: thick. of the surf. layer in temp. comp.
+   REAL(wp), PUBLIC ::   parsub     =  1.0         !: switch for snow sublimation or not
+   REAL(wp), PUBLIC ::   alphs      =  1.0         !: coef. for snow density when snow-ice formation
+   REAL(wp), PUBLIC ::   hiccrit(2) = (/0.3,0.3/)  !: ice th. for lateral accretion in the NH (SH) (m)
+
+   !---------------------------------------------------!
+   !                Sea-Ice variables                  !
+   !---------------------------------------------------!
+
+   REAL(wp), PUBLIC ::   usecc2           !:  = 1.0 / ( ecc * ecc )
+   REAL(wp), PUBLIC ::   rhoco            !: = rau0 * cw
+   REAL(wp), PUBLIC ::   sangvg, cangvg   !: sin and cos of the turning angle for ocean stress
+   REAL(wp), PUBLIC ::   pstarh           !: pstar / 2.0
+   REAL(wp), PUBLIC ::   uscomi           !: inverse of minimum lead fraction      !!gm to be suppress
+   REAL(wp), PUBLIC ::   cnscg            !: ratio  rcpsn/rcpic                    !!gm to be suppress
 
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv   !: hor. diffusivity coeff. at ocean U- and V-points (m2/s)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv   !: ice hor. eddy diffusivity coef. at ocean U- and V-points
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm   !: mean snow and ice thicknesses
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s                 !: friction velocity
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s         !: friction velocity
 
-   !!* diagnostic quantities
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tfu           !: Freezing/Melting point temperature of sea water at SSS
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicif         !: Ice thickness
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnif         !: Snow thickness
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicifp        !: Ice production/melting
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   frld          !: Leads fraction = 1-a/totalarea
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   phicif        !: ice thickness  at previous time 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pfrld         !: Leads fraction at previous time  
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qstoif        !: Energy stored in the brine pockets
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fbif          !: Heat flux at the ice base
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmsnif       !: Variation of snow mass
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdtcn         !: net downward heat flux from the ice to the ocean
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   thcm          !: part of the solar energy used in the lead heat budget
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fstric        !: Solar flux transmitted trough the ice
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max heat contained in brine pockets (?)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fscmbq        !: Linked with the solar flux below the ice (?)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fsbbq         !: Also linked with the solar flux below the ice (?)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy in case of toral lateral ablation (?)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   dmgwi         !: Variation of the mass of snow ice
+   INTEGER , PUBLIC, DIMENSION(jpij) ::   npb         !: number of points where computations has to be done
+   INTEGER , PUBLIC, DIMENSION(jpij) ::   npac        !: correspondance between the points
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice   !: two components of the ice   velocity at I-point (m/s)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce   !: two components of the ocean velocity at I-point (m/s)
+   !                                                         !!* sea-ice variables *
+   REAL(wp), PUBLIC ::   sist   (jpi,jpj), sist_1d   (jpij)   !: Sea-Ice Surface Temperature (Kelvin)
+   REAL(wp), PUBLIC ::   tfu    (jpi,jpj), tfu_1d    (jpij)   !: Freezing/Melting point temperature of sea water at SSS
+   REAL(wp), PUBLIC ::   hicif  (jpi,jpj), h_ice_1d  (jpij)   !: ice thickness
+   REAL(wp), PUBLIC ::   phicif (jpi,jpj)                     !: ice thickness  at previous time 
+   REAL(wp), PUBLIC ::   hsnif  (jpi,jpj), h_snow_1d (jpij)   !: Snow thickness
+   REAL(wp), PUBLIC ::   hicifp (jpi,jpj)                     !: Ice production/melting
+   REAL(wp), PUBLIC ::   frld   (jpi,jpj), frld_1d   (jpij)   !: Leads fraction = 1-a/totalarea
+   REAL(wp), PUBLIC ::   pfrld  (jpi,jpj)                     !: Leads fraction at previous time  
+   REAL(wp), PUBLIC ::   qstoif (jpi,jpj), qstbif_1d (jpij)   !: Energy stored in the brine pockets         !!gm 
+   REAL(wp), PUBLIC ::   fbif   (jpi,jpj), fbif_1d   (jpij)   !: Heat flux at the ice base
+   REAL(wp), PUBLIC ::   rdmsnif(jpi,jpj), rdmsnif_1d(jpij)   !: Variation of snow mass
+   REAL(wp), PUBLIC ::   rdmicif(jpi,jpj), rdmicif_1d(jpij)   !: Variation of ice mass
+   REAL(wp), PUBLIC ::   qldif  (jpi,jpj), qldif_1d  (jpij)   !: heat balance of the lead (or of the open ocean)
+   REAL(wp), PUBLIC ::   qcmif  (jpi,jpj), qcmif_1d  (jpij)   !: Energy needed to freeze the ocean surface layer 
+   REAL(wp), PUBLIC ::   fdtcn  (jpi,jpj)                     !: net downward heat flux from the ice to the ocean
+   REAL(wp), PUBLIC ::   qdtcn  (jpi,jpj)                     !: energy from the ice to the ocean point (at a factor 2)
+   REAL(wp), PUBLIC ::   thcm   (jpi,jpj), thcm_1d   (jpij)   !: part of the solar energy used in the lead heat budget
+   REAL(wp), PUBLIC ::   fstric (jpi,jpj), fstric_1d (jpij)   !: Solar flux transmitted trough the ice
+   REAL(wp), PUBLIC ::   ffltbif(jpi,jpj), fltbif_1d (jpij)   !: Array linked with the max heat contained in brine pockets (?)
+   REAL(wp), PUBLIC ::   fscmbq (jpi,jpj), fscbq_1d  (jpij)   !: Linked with the solar flux below the ice (?)
+   REAL(wp), PUBLIC ::   fsbbq  (jpi,jpj)                     !: Also linked with the solar flux below the ice (?)
+   REAL(wp), PUBLIC ::   qfvbq  (jpi,jpj), qfvbq_1d  (jpij)   !: Array used to store energy in case of total lateral ablation (?)
+   REAL(wp), PUBLIC ::   dmgwi(jpi,jpj), dmgwi_1d(jpij)       !: Variation of the mass of snow ice
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   tbif  !: Temperature inside the ice/snow layer
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice     !: two components of the ice   velocity at I-point (m/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce     !: two components of the ocean velocity at I-point (m/s)
 
-   !!* moment used in the advection scheme
+   REAL(wp), PUBLIC ::   tbif(jpi,jpj,jplayersp1), tbif_1d(jpij,jplayersp1)  !: Temperature inside the ice+snow layer
+
+   ! Surface forcing transford into  1D field  (2D field defined in /OPA_SRC/SBC/sbc_ice.F90)
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   qns_ice_1d   !:   qns_ice   
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   qsr_ice_1d   !:   qsr_ice
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   qla_ice_1d   !:   qla_ice
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   dqla_ice_1d  !:   dqla_ice
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   dqns_ice_1d  !:   dqns_ice
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   fr1_i0_1d    !:   fr1_i0
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   fr2_i0_1d    !:   fr2_i0
+
+   ! Surface forcing transformed into 1D field  (2D field defined in /OPA_SRC/SBC/sbc_oce.F90)
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   sprecip_1d   !:   sprecip
+   
+   ! local variable  transformed into 1D field  (2D field defined in /LIM_SRC_2/limthd_2.F90)
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   qlbbq_1d     !:   zqlbsbq
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   dvsbq_1d     !:   zdvosif
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   dvbbq_1d     !:   zdvobif
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   dvlbq_1d     !:   zdvolif
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   dvnbq_1d     !:   zdvonif
+   
+   
+   REAL(wp), PUBLIC, DIMENSION(jpij) ::   rdvomif_1d    !:   rdvomif       <<<== !gm should be suppressed !
+
+   !---------------------------------------------------!
+   !   Prather' advection scheme : 1st & 2nd moments   !
+   !---------------------------------------------------!
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxice, syice, sxxice, syyice, sxyice   !: for ice  volume
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: for snow volume
@@ -102,4 +152,8 @@
 #endif
 
+   !!-------- -------------------------------------------------------------
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE ice_2
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/iceini_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/iceini_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/iceini_2.F90	(revision 1855)
@@ -4,6 +4,6 @@
    !!   Sea-ice model : LIM 2.0 Sea ice model Initialization
    !!======================================================================
-   !! History :   1.0  !  02-08  (G. Madec)  F90: Free form and modules
-   !!             2.0  !  03-08  (C. Ethe)  add ice_run
+   !! History :  1.0  !  2002-08  (G. Madec)  F90: Free form and modules
+   !!            2.0  !  2003-08  (C. Ethe)  add ice_run
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -11,28 +11,26 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
+   !!   ice_init_2      : sea-ice model initialization
    !!----------------------------------------------------------------------
-   !!   ice_init_2       : sea-ice model initialization
-   !!   ice_run_2        : Definition some run parameter for ice model
-   !!----------------------------------------------------------------------
-   USE dom_oce
-   USE dom_ice_2
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE sbc_ice         ! surface boundary condition: ice
-   USE phycst          ! Define parameters for the routines
-   USE ice_2
-   USE limmsh_2
-   USE limistate_2
-   USE limrst_2   
-   USE in_out_manager
+   USE dom_oce          ! ocean domain
+   USE dom_ice_2        ! LIM-2 ice domain
+   USE sbc_oce          ! surface boundary condition: ocean
+   USE sbc_ice          ! surface boundary condition: sea ice
+   USE phycst           ! physical constant
+   USE ice_2            ! LIM-2 ice varibles
+   USE limmsh_2         ! LIM-2 mesh
+   USE limistate_2      ! LIM-2 initial state
+   USE limrst_2         ! LIM-2 restart
+   USE in_out_manager   ! I/O manager
       
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   ice_init_2               ! called by sbcice_lim_2.F90
+   PUBLIC   ice_init_2   ! called by sbcice_lim_2.F90
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -45,45 +43,14 @@
       !! ** purpose :   
       !!----------------------------------------------------------------------
-      !
-      ! Open the namelist file 
-      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )      
-      CALL ice_run_2                    !  read in namelist some run parameters
-                 
-      ! Louvain la Neuve Ice model
-      rdt_ice = nn_fsbc * rdttra(1)
-
-      CALL lim_msh_2                  ! ice mesh initialization
-     
-      ! Initial sea-ice state
-      IF( .NOT.ln_rstart ) THEN
-         CALL lim_istate_2            ! start from rest: sea-ice deduced from sst
-      ELSE
-         CALL lim_rst_read_2          ! start from a restart file
-      ENDIF
-      
-      tn_ice(:,:,1) = sist(:,:)         ! initialisation of ice temperature   
-      fr_i  (:,:) = 1.0 - frld(:,:)   ! initialisation of sea-ice fraction    
-      !
-   END SUBROUTINE ice_init_2
-
-
-   SUBROUTINE ice_run_2
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_run_2 ***
-      !!                 
-      !! ** Purpose :   Definition some run parameter for ice model
-      !!
-      !! ** Method  :   Read the namicerun namelist and check the parameter 
-      !!       values called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namelist namicerun
-      !!-------------------------------------------------------------------
       NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, ln_limdmp, acrit, hsndif, hicdif
       !!-------------------------------------------------------------------
-      !                    
-      REWIND ( numnam_ice )                       ! Read Namelist namicerun 
+      !
+      !                                      ! Open the ice namelist file 
+      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )    
+        
+      REWIND ( numnam_ice )                  ! Read namicerun namelist
       READ   ( numnam_ice , namicerun )
-
-      IF(lwp) THEN
+      !
+      IF(lwp) THEN                           ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
@@ -95,6 +62,19 @@
          WRITE(numout,*) '   computation of temp. in ice  (=0) or not (=9999) hicdif = ', hicdif
       ENDIF
+                       
+      rdt_ice = nn_fsbc * rdttra(1)          ! ice time step 
+
+      CALL lim_msh_2                         ! ice mesh initialization
+     
+      !                                      ! Initial sea-ice state
+      IF( .NOT.ln_rstart ) THEN   ;   CALL lim_istate_2        ! start from rest: sea-ice deduced from sst
+      ELSE                        ;   CALL lim_rst_read_2      ! restart from a file
+      ENDIF
+
+      !                                      ! ice variables see by the ocean      
+      tn_ice(:,:,1) = sist(:,:)                 ! surface ice temperature   
+      fr_i  (:,:)   = 1.0 - frld(:,:)           ! sea-ice fraction    
       !
-   END SUBROUTINE ice_run_2
+   END SUBROUTINE ice_init_2
 
 #else
@@ -102,7 +82,4 @@
    !!   Default option :        Empty module       NO LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE ice_init_2      ! Empty routine
-   END SUBROUTINE ice_init_2
 #endif
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limadv_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limadv_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limadv_2.F90	(revision 1855)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                       ***  MODULE limadv_2   ***
-   !! LIM 2.0 sea-ice model : sea-ice advection
+   !! LIM-2 sea-ice model : sea-ice advection
    !!======================================================================
    !! History :  OPA  ! 2000-01 (LIM)  Original code
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdia_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdia_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdia_2.F90	(revision 1855)
@@ -12,5 +12,4 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
    !!   lim_dia_2      : computation of the time evolution of keys var.
    !!   lim_dia_init_2 : initialization and namelist read
@@ -28,25 +27,23 @@
    PRIVATE
 
-   PUBLIC               lim_dia_2          ! called by sbc_ice_lim_2
-   INTEGER, PUBLIC ::   ntmoy   = 1 ,   &  !: instantaneous values of ice evolution or averaging ntmoy
-      &                 ninfo   = 1        !: frequency of ouputs on file ice_evolu in case of averaging
-
-   INTEGER, PARAMETER ::   &  ! Parameters for outputs to files "evolu"
-      jpinfmx = 100         ,    &  ! maximum number of key variables
-      jpchinf = 5           ,    &  ! ???
-      jpchsep = jpchinf + 2         ! ???
-
-   INTEGER ::   &
-      nfrinf  = 4 ,     &  ! number of variables written in one line 
-      nferme ,          &  ! last time step at which the var. are written on file
-      nvinfo ,          &  ! number of total variables 
-      nbvt   ,          &  ! number of time variables
-      naveg                ! number of step for accumulation before averaging
-
-   CHARACTER(len= 8) ::   fmtinf  = '1PE13.5 ' ! format of the output values  
-   CHARACTER(len=30) ::   fmtw  ,           &  ! formats
-      &                   fmtr  ,           &  ! ???
-      &                   fmtitr               ! ???
-   CHARACTER(len=jpchsep), DIMENSION(jpinfmx) ::   titvar               ! title of key variables
+   PUBLIC   lim_dia_2   ! called by sbc_ice_lim_2
+   
+   INTEGER, PUBLIC ::   ntmoy   = 1   !: instantaneous values of ice evolution or averaging ntmoy
+   INTEGER, PUBLIC ::   ninfo   = 1   !: frequency of ouputs on file ice_evolu in case of averaging
+
+   !                                              !!! Parameters for outputs to files "evolu"
+   INTEGER, PARAMETER ::   jpinfmx = 100           ! maximum number of key variables
+   INTEGER, PARAMETER ::   jpchinf = 5             ! ???
+   INTEGER, PARAMETER ::   jpchsep = jpchinf + 2   ! ???
+
+   INTEGER ::   nfrinf  = 4   ! number of variables written in one line 
+   INTEGER ::   nferme        ! last time step at which the var. are written on file
+   INTEGER ::   nvinfo        ! number of total variables 
+   INTEGER ::   nbvt          ! number of time variables
+   INTEGER ::   naveg         ! number of step for accumulation before averaging
+
+   CHARACTER(len= 8) ::   fmtinf  = '1PE13.5 '   ! format of the output values  
+   CHARACTER(len=30) ::   fmtw, fmtr, fmtitr     ! various formats
+   CHARACTER(len=jpchsep), DIMENSION(jpinfmx) ::   titvar   ! title of key variables
  
    REAL(wp)                     ::   epsi06 = 1.e-06      ! ???
@@ -57,5 +54,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -71,12 +68,12 @@
       !!      the temporal evolution of some key variables
       !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! number of iteration
+      INTEGER, INTENT(in) ::   kt   ! number of iteration
       !!
       INTEGER  ::   jv,ji, jj   ! dummy loop indices
       INTEGER  ::   nv          ! indice of variable 
-      REAL(wp) ::   zarea    , zldarea  ,    &  ! sea-ice and leads area
-         &          zextent15, zextent85,    &  ! sea-ice extent (15% and 85%)
-         &          zicevol  , zsnwvol  ,    &  ! sea-ice and snow volume volume
-         &          zicespd                     ! sea-ice velocity
+      REAL(wp) ::   zarea    , zldarea     ! sea-ice and leads area
+      REAL(wp) ::   zextent15, zextent85   ! sea-ice extent (15% and 85%)
+      REAL(wp) ::   zicevol  , zsnwvol     ! sea-ice and snow volume volume
+      REAL(wp) ::   zicespd                     ! sea-ice velocity
       REAL(wp), DIMENSION(jpinfmx) ::   vinfor  ! temporary working space 
       !!-------------------------------------------------------------------
@@ -120,5 +117,4 @@
       vinfor(nv+13) = SQRT( vinfor(nv+13) / MAX( vinfor(nv+9) , epsi06 ) )
 
-
       ! variables in southern Hemis
        nv = nv + 1
@@ -177,13 +173,12 @@
        INTEGER  ::   nv            ! indice of variable 
        REAL(wp) ::   zxx0, zxx1    ! temporary scalars
-
+       !!
        NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy
        !!-------------------------------------------------------------------
 
-       ! Read Namelist namicedia
-       REWIND ( numnam_ice )
+       REWIND ( numnam_ice )                 ! Read Namelist namicedia
        READ   ( numnam_ice  , namicedia )
        
-       IF(lwp) THEN
+       IF(lwp) THEN                          ! control print
           WRITE(numout,*)
           WRITE(numout,*) 'lim_dia_init_2 : ice parameters for ice diagnostics '
@@ -274,7 +269,4 @@
    !!   Default option :                           NO LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_dia_2         ! Empty routine
-   END SUBROUTINE lim_dia_2
 #endif
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdmp_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 1855)
@@ -4,5 +4,5 @@
    !!  Ice model : restoring Ice thickness and Fraction leads
    !!======================================================================
-   !! History :   2.0  !  04-04 (S. Theetten) Original code
+   !! History :   2.0  ! 2004-04 (S. Theetten) Original code
    !!----------------------------------------------------------------------
 #if defined key_lim2   &&   defined key_tradmp
@@ -11,14 +11,13 @@
    !!   'key_tradmp'                                                Damping
    !!----------------------------------------------------------------------
+   !!   lim_dmp_2       : ice model damping
    !!----------------------------------------------------------------------
-   !!   lim_dmp_2      : ice model damping
-   !!----------------------------------------------------------------------
-   USE in_out_manager  ! I/O manager
-   USE phycst          ! physical constants
-   USE ice_2
-   USE tradmp
-   USE dom_oce
-   USE oce
-   USE iom
+   USE oce              ! ocean variables
+   USE dom_oce          ! ocean domain
+   USE phycst           ! physical constants
+   USE ice_2            ! LIM-2 variables
+   USE tradmp           ! traceur damping
+   USE in_out_manager   ! I/O manager
+   USE iom              !
    
    IMPLICIT NONE
@@ -27,15 +26,15 @@
    PUBLIC   lim_dmp_2     ! called by ice_step_2
    
-   INTEGER                        ::   nice1, nice2,  &  ! first and second record used
-      &                                inumice_dmp       ! logical unit for ice variables (damping)
-   REAL(wp), DIMENSION(jpi,jpj)   ::   hicif_dta  ,   &  ! ice thickness at a given time
-      &                                frld_dta          ! fraction lead at a given time
-   REAL(wp), DIMENSION(jpi,jpj,2) ::   hicif_data ,   &  ! ice thickness data at two consecutive times
-      &                                frld_data         ! fraction lead data at two consecutive times
+   INTEGER                        ::   nice1, nice2   ! first and second record used
+   INTEGER                        ::   inumice_dmp    ! logical unit for ice variables (damping)
+   REAL(wp), DIMENSION(jpi,jpj)   ::   hicif_dta      ! ice thickness at a given time
+   REAL(wp), DIMENSION(jpi,jpj)   ::   frld_dta       ! fraction lead at a given time
+   REAL(wp), DIMENSION(jpi,jpj,2) ::   hicif_data     ! ice thickness data at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj,2) ::   frld_data      ! fraction lead data at two consecutive times
 
    !! * Substitution
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0 , UCL-LOCEAN-IPSL  (2006)
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -44,5 +43,5 @@
 CONTAINS
 
-   SUBROUTINE lim_dmp_2(kt)
+   SUBROUTINE lim_dmp_2( kt )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE lim_dmp_2 ***
@@ -53,15 +52,15 @@
       !! ** method  : the key_tradmp must be used to compute resto(:,:) coef.
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! ocean time-step
-      !
-      INTEGER             ::   ji, jj         ! dummy loop indices
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step
+      !!
+      INTEGER ::   ji, jj   ! dummy loop indices
       !!---------------------------------------------------------------------
       !
       CALL dta_lim_2( kt )
-
+      !
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
             hicif(ji,jj) = hicif(ji,jj) - rdt_ice * resto(ji,jj,1) * ( hicif(ji,jj) - hicif_dta(ji,jj) )
-            frld(ji,jj)  = frld (ji,jj) - rdt_ice * resto(ji,jj,1) * ( frld(ji,jj)  - frld_dta (ji,jj) )  
+            frld(ji,jj)  = frld (ji,jj) - rdt_ice * resto(ji,jj,1) * ( frld (ji,jj) - frld_dta (ji,jj) )  
          END DO
       END DO
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdyn_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 1855)
@@ -17,10 +17,10 @@
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
-   USE sbc_oce        !
-   USE phycst         !
-   USE ice_2          !
-   USE dom_ice_2      !
-   USE limistate_2    !
-   USE limrhg_2       ! ice rheology
+   USE sbc_oce        ! surface boundary condition : ocean
+   USE phycst         ! physical constants
+   USE ice_2          ! LIM-2 sea-ice variables
+   USE dom_ice_2      ! LIM-2 domain
+   USE limistate_2    ! LIM-2 initial state
+   USE limrhg_2       ! LIM-2 rheology
 
    USE lbclnk         !
@@ -34,10 +34,9 @@
    PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim
 
-   !! * Module variables
    REAL(wp)  ::  rone    = 1.e0   ! constant value
 
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -59,9 +58,9 @@
       !!              - treatment of the case if no ice dynamic
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! number of iteration
-      !!
-      INTEGER  ::   ji, jj             ! dummy loop indices
-      INTEGER  ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology
-      REAL(wp) ::   zcoef              ! temporary scalar
+      INTEGER, INTENT(in) ::   kt   ! number of iteration
+      !!
+      INTEGER  ::   ji, jj        ! dummy loop indices
+      INTEGER  ::   i_j1, i_jpj   ! Starting/ending j-indices for rheology
+      REAL(wp) ::   zcoef         ! temporary scalar
       REAL(wp), DIMENSION(jpj)     ::   zind           ! i-averaged indicator of sea-ice
       REAL(wp), DIMENSION(jpj)     ::   zmsk           ! i-averaged of tmask
@@ -144,5 +143,5 @@
          ENDIF
 
-         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
+         IF(ln_ctl)   CALL prt_ctl( tab2d_1=u_ice, clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :' )
          
          ! computation of friction velocity
@@ -179,5 +178,5 @@
       CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
       !
-      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ust2s, clinfo1=' lim_dyn  : ust2s :' )
 
    END SUBROUTINE lim_dyn_2
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limhdf_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 1855)
@@ -4,4 +4,8 @@
    !! LIM 2.0 ice model : horizontal diffusion of sea-ice quantities
    !!======================================================================
+   !! History :  LIM  ! 2000-01 (UCL)  Original code
+   !!            1.0  ! 2001-05 (G. Madec, R. Hordoir) opa norm
+   !!            2.0  ! 2002-08 (C. Ethe)  F90, free form
+   !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
@@ -10,29 +14,26 @@
    !!   lim_hdf_2  : diffusion trend on sea-ice variable
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE dom_oce
-   USE in_out_manager
-   USE ice_2
-   USE lbclnk
-   USE lib_mpp
-   USE prtctl          ! Print control
+   USE dom_oce          ! ocean domain
+   USE ice_2            ! LIM-2 variables
+   USE prtctl           ! Print control
+   USE lbclnk           !
+   USE lib_mpp          !
+   USE in_out_manager   ! I/O manager
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_hdf_2    ! called by lim_tra_2
+   PUBLIC   lim_hdf_2   ! called by lim_tra_2
 
-   !! * Module variables
-   LOGICAL  ::   linit = .TRUE.              ! ???
+   LOGICAL  ::   linit = .TRUE.              ! indictor of initialisation
    REAL(wp) ::   epsi04 = 1e-04              ! constant
-   REAL(wp), DIMENSION(jpi,jpj) ::   zfact   ! ???
+   REAL(wp), DIMENSION(jpi,jpj) ::   zfact   ! metric term
 
    !! * Substitution 
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -43,33 +44,19 @@
       !!                  ***  ROUTINE lim_hdf_2  ***
       !!
-      !! ** purpose :   Compute and add the diffusive trend on sea-ice
-      !!      variables
+      !! ** purpose :   Compute and add the diffusive trend on sea-ice variables
       !!
       !! ** method  :   Second order diffusive operator evaluated using a
-      !!      Cranck-Nicholson time Scheme.
+      !!              Cranck-Nicholson time Scheme.
       !!
-      !! ** Action  :    update ptab with the diffusive contribution
+      !! ** Action  :   update ptab with the diffusive contribution
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   ptab   ! Field on which the diffusion is applied  
       !!
-      !! History :
-      !!        !  00-01 (LIM) Original code
-      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
-      !!        !  02-08 (C. Ethe)  F90, free form
-      !!-------------------------------------------------------------------
-      ! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   &
-         ptab                 ! Field on which the diffusion is applied  
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         ptab0                ! ???
-
-      ! * Local variables
       INTEGER ::  ji, jj      ! dummy loop indices
-      INTEGER ::  &
-         its, iter            ! temporary integers
+      INTEGER ::  its, iter   ! temporary integers
       CHARACTER (len=55) :: charout
-      REAL(wp) ::  &
-         zalfa, zrlxint, zconv, zeps   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::  & 
-         zrlx, zflu, zflv, &  ! temporary workspaces
-         zdiv0, zdiv          !    "           "
+      REAL(wp) ::   zalfa, zrlxint, zconv, zeps   ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zrlx, zflu, zflv, zdiv0, zdiv  ! temporary workspaces
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztab0                ! ???
       !!-------------------------------------------------------------------
 
@@ -82,7 +69,5 @@
 
       ! Arrays initialization
-      ptab0 (:, : ) = ptab(:,:)
-!bug  zflu (:,jpj) = 0.e0
-!bug  zflv (:,jpj) = 0.e0
+      ztab0(:, : ) = ptab(:,:)
       zdiv0(:, 1 ) = 0.e0
       zdiv0(:,jpj) = 0.e0
@@ -98,6 +83,5 @@
          DO jj = 2, jpjm1  
             DO ji = fs_2 , fs_jpim1   ! vector opt.
-               zfact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj  ) + e1v(ji,jj) + e1v(ji,jj-1) ) &
-                  &          / ( e1t(ji,jj) * e2t(ji,jj) )
+               zfact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) )
             END DO
          END DO
@@ -110,12 +94,10 @@
       iter  = 0
 
-      !                                                   !===================
-      DO WHILE ( ( zconv > zeps ) .AND. (iter <= its) )   ! Sub-time step loop
-         !                                                !===================
-         ! incrementation of the sub-time step number
-         iter = iter + 1
+      !                                                   !======================!
+      DO WHILE ( ( zconv > zeps ) .AND. (iter <= its) )   !  Sub-time step loop  !
+         !                                                !======================!
+         iter = iter + 1               ! incrementation of the sub-time step number
 
-         ! diffusive fluxes in U- and V- direction
-         DO jj = 1, jpjm1
+         DO jj = 1, jpjm1              ! diffusive fluxes in U- and V- direction
             DO ji = 1 , fs_jpim1   ! vector opt.
                zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
@@ -123,7 +105,5 @@
             END DO
          END DO
-
-         ! diffusive trend : divergence of the fluxes
-         DO jj= 2, jpjm1
+         DO jj= 2, jpjm1               ! diffusive trend : divergence of the fluxes
             DO ji = fs_2 , fs_jpim1   ! vector opt. 
                zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   &
@@ -131,12 +111,10 @@
             END DO
          END DO
-
-         ! save the first evaluation of the diffusive trend in zdiv0
+         !                             ! save the first evaluation of the diffusive trend in zdiv0
          IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)       
 
-         ! XXXX iterative evaluation?????
-         DO jj = 2, jpjm1
+         DO jj = 2, jpjm1              ! XXXX iterative evaluation?????
             DO ji = fs_2 , fs_jpim1   ! vector opt.
-               zrlxint = (   ptab0(ji,jj)    &
+               zrlxint = (   ztab0(ji,jj)    &
                   &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + zfact(ji,jj) * ptab(ji,jj) )   &
                   &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             & 
@@ -145,10 +123,7 @@
             END DO
          END DO
-
-         ! lateral boundary condition on ptab
          CALL lbc_lnk( zrlx, 'T', 1. )
 
-         ! convergence test
-         zconv = 0.e0
+         zconv = 0.e0                  ! convergence test
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -157,17 +132,16 @@
          END DO
          IF( lk_mpp )   CALL mpp_max( zconv )   ! max over the global domain
-
-         ptab(:,:) = zrlx(:,:)
-
-         !                                         !==========================
-      END DO                                       ! end of sub-time step loop
-      !                                            !==========================
+         !
+         ptab(:,:) = zrlx(:,:)         ! update value
+         !                                         !=============================!
+      END DO                                       !  end of sub-time step loop  !
+      !                                            !=============================!
 
       IF(ln_ctl)   THEN
-         zrlx(:,:) = ptab(:,:) - ptab0(:,:)
+         zrlx(:,:) = ptab(:,:) - ztab0(:,:)
          WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
-         CALL prt_ctl(tab2d_1=zrlx, clinfo1=charout)
+         CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout )
       ENDIF
-
+      !
    END SUBROUTINE lim_hdf_2
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limistate_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limistate_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limistate_2.F90	(revision 1855)
@@ -4,9 +4,9 @@
    !!              Initialisation of diagnostics ice variables
    !!======================================================================
-   !! History :   1.0  !  01-04  (C. Ethe, G. Madec)  Original code
-   !!             2.0  !  03-08  (G. Madec)  add lim_istate_init
-   !!                  !  04-04  (S. Theetten) initialization from a file
-   !!                  !  06-07  (S. Masson)  IOM to read the restart
-   !!                  !  07-10  (G. Madec)  surface module
+   !! History :   1.0  !  2001-04  (C. Ethe, G. Madec)  Original code
+   !!             2.0  !  2003-08  (G. Madec)  add lim_istate_init
+   !!              -   !  2004-04  (S. Theetten) initialization from a file
+   !!              -   !  2006-07  (S. Masson)  IOM to read the restart
+   !!              -   !  2007-10  (G. Madec)  surface module
    !!--------------------------------------------------------------------
 #if defined key_lim2
@@ -14,25 +14,24 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
    !!   lim_istate_2      :  Initialisation of diagnostics ice variables
    !!   lim_istate_init_2 :  initialization of ice state and namelist read
    !!----------------------------------------------------------------------
-   USE phycst
-   USE par_ice_2       ! ice parameters
-   USE dom_ice_2
-   USE eosbn2          ! equation of state
-   USE lbclnk
-   USE oce
-   USE ice_2
-   USE iom
-   USE in_out_manager
+   USE oce              ! ocean variables
+   USE ice_2            ! LIM-2 variables
+   USE par_ice_2        ! LIM-2 ice parameters
+   USE dom_ice_2        ! LIM-2 domain
+   USE phycst           ! physical constants
+   USE eosbn2           ! equation of state
+   USE lbclnk           ! 
+   USE iom              !
+   USE in_out_manager   !
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC lim_istate_2      ! routine called by lim_init_2.F90
-
-!!! ** init namelist (namiceini) **
-   LOGICAL  ::   ln_limini = .FALSE.  !: Ice initialization state
+   PUBLIC   lim_istate_2   ! routine called by lim_init_2.F90
+
+   !                                 !!! ** init namelist (namiceini) **
+   LOGICAL  ::   ln_limini = .FALSE.  ! Ice initialization state
    REAL(wp) ::   ttest     = 2.0      ! threshold water temperature for initial sea ice
    REAL(wp) ::   hninn     = 0.5      ! initial snow thickness in the north
@@ -45,6 +44,7 @@
    REAL(wp) ::   zero      = 0.e0     ! constant value = 0
    REAL(wp) ::   zone      = 1.e0     ! constant value = 1
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -62,6 +62,6 @@
       !!                or from arbitrary sea-ice conditions
       !!--------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk                ! dummy loop indices
-      REAL(wp) ::   zidto                     ! temporary scalar
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zidto        ! temporary scalar
       !--------------------------------------------------------------------
  
@@ -69,7 +69,7 @@
 
       IF( .NOT. ln_limini ) THEN  
-         
+         !
          tfu(:,:) = tfreez( sn(:,:,1) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
-
+         !
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -126,6 +126,5 @@
 
       !-- lateral boundary conditions
-      CALL lbc_lnk( hicif, 'T', 1. )
-      CALL lbc_lnk( frld , 'T', 1. )
+      CALL lbc_lnk( hicif, 'T', 1. )   ;   CALL lbc_lnk( frld , 'T', 1. )
 
       ! C A U T I O N  frld = 1 over land and lbc_lnk put zero along 
@@ -140,5 +139,5 @@
       CALL lbc_lnk( fsbbq  , 'T', 1. )
       CALL lbc_lnk( qstoif , 'T', 1. )
-
+      !
    END SUBROUTINE lim_istate_2
 
@@ -151,13 +150,13 @@
       !!
       !! ** Method  :   Read the namiceini namelist and check the parameter 
-      !!       values called at the first timestep (nit000)
+      !!              values called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist namiceini
       !!-------------------------------------------------------------------
-      INTEGER :: inum_ice
-      INTEGER :: ji,jj
-
+      INTEGER ::   ji,jj      ! dummy loop indices
+      INTEGER ::   inum_ice   ! temporary integer
+      !!
       NAMELIST/namiceini/ ln_limini, ttest, hninn, hginn, alinn, &
-         &                hnins, hgins, alins
+         &                                  hnins, hgins, alins
       !!-------------------------------------------------------------------
       !
@@ -165,5 +164,5 @@
       READ   ( numnam_ice , namiceini )
       !
-      IF(lwp) THEN
+      IF(lwp) THEN                        ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'lim_istate_init_2 : ice parameters inititialisation '
@@ -179,5 +178,5 @@
       ENDIF
 
-      IF( ln_limini ) THEN                      ! Ice initialization using input file
+      IF( ln_limini ) THEN                ! Ice initialization using input file
          !
          CALL iom_open( 'Ice_initialization.nc', inum_ice )
@@ -186,5 +185,5 @@
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '                  ice state initialization with : Ice_initialization.nc'
-            
+            !
             CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif )      
             CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif )      
@@ -192,9 +191,9 @@
             CALL iom_get( inum_ice, jpdom_data, 'ts'   , sist  )
             CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif(1:nlci,1:nlcj,:),   &
-                 &        kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) )
+               &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) )
             ! put some values in the extra-halo...
             DO jj = nlcj+1, jpj   ;   tbif(1:nlci,jj,:) = tbif(1:nlci,nlej,:)   ;   END DO
             DO ji = nlci+1, jpi   ;   tbif(ji    ,: ,:) = tbif(nlei  ,:   ,:)   ;   END DO
-
+            !
             CALL iom_close( inum_ice)
             !
@@ -208,7 +207,4 @@
    !!   Default option :         Empty module      NO LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_istate_2        ! Empty routine
-   END SUBROUTINE lim_istate_2
 #endif
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limmsh_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 1855)
@@ -4,4 +4,7 @@
    !! LIM 2.0 ice model :   definition of the ice mesh parameters
    !!======================================================================
+   !! ** History :  LIM  ! 2001-04 (UCL)  original code
+   !!               2.0  ! 2002-08 (C. Ethe, G. Madec) F90, module
+   !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
@@ -10,5 +13,4 @@
    !!   lim_msh_2   : definition of the ice mesh
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE phycst
    USE dom_oce
@@ -20,11 +22,10 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC lim_msh_2      ! routine called by ice_ini_2.F90
-
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   PUBLIC   lim_msh_2   ! routine called by ice_ini_2.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -42,20 +43,10 @@
       !!              - Initialization of the ice masks (tmsk, umsk)
       !! 
-      !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10 
-      !!
-      !! ** History :
-      !!         original    : 01-04 (LIM)
-      !!         addition    : 02-08 (C. Ethe, G. Madec)
+      !! Reference  : Deleersnijder et al. Ocean Modelling 100, 7-10 
       !!--------------------------------------------------------------------- 
-      !! * Local variables
-      INTEGER :: ji, jj      ! dummy loop indices
-
-      REAL(wp), DIMENSION(jpi,jpj) ::  &
-         zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction
-         !                   ! (resp. y direction) (defined at the center)
-      REAL(wp) ::         &
-         zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid
-         zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid
-         zusden, zusden2     ! temporary scalars
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   zh1p, zd1d2p, zusden    ! temporary scalars
+      REAL(wp) ::   zh2p, zd2d1p, zusden2   !    -         -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zd2d1 , zd1d2   ! 2D workspace
       !!---------------------------------------------------------------------
 
@@ -112,8 +103,8 @@
       !-------------------
 !!ibug ???
-      akappa(:,:,:,:) = 0.e0
-      wght(:,:,:,:) = 0.e0
+      akappa(:,:,:,:)     = 0.e0
+      wght  (:,:,:,:)     = 0.e0
       alambd(:,:,:,:,:,:) = 0.e0
-      tmu(:,:) = 0.e0
+      tmu   (:,:)         = 0.e0
 !!i
       
@@ -238,10 +229,7 @@
          END DO
       END DO
-      
-      !--lateral boundary conditions    
-      CALL lbc_lnk( tmu(:,:), 'I', 1. )
-      
-      ! unmasked and masked area of T-grid cell
-      area(:,:) = e1t(:,:) * e2t(:,:)
+      CALL lbc_lnk( tmu(:,:), 'I', 1. )      !--lateral boundary conditions    
+      
+      area(:,:) = e1t(:,:) * e2t(:,:)        ! unmasked and masked area of T-grid cell
       
    END SUBROUTINE lim_msh_2
@@ -251,7 +239,4 @@
    !!   Default option            Dummy Module         NO LIM sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_msh_2           ! Dummy routine
-   END SUBROUTINE lim_msh_2
 #endif
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limrhg_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 1855)
@@ -4,9 +4,9 @@
    !!   Ice rheology :  performs sea ice rheology
    !!======================================================================
-   !! History :  0.0  !  93-12  (M.A. Morales Maqueda.)  Original code
-   !!            1.0  !  94-12  (H. Goosse) 
-   !!            2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp
-   !!            " "  !  06-08  (G. Madec)  surface module, ice-stress at I-point
-   !!            " "  !  09-09  (G. Madec)  Huge verctor optimisation
+   !! History :  0.0  !  1993-12  (M.A. Morales Maqueda.)  Original code
+   !!            1.0  !  1994-12  (H. Goosse) 
+   !!            2.0  !  1903-12  (C. Ethe, G. Madec)  F90, mpp
+   !!             -   !  2006-08  (G. Madec)  surface module, ice-stress at I-point
+   !!             -   !  2009-09  (G. Madec)  Huge verctor optimisation
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -14,6 +14,5 @@
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   lim_rhg_2   : computes ice velocities
+   !!   lim_rhg_2     : computes ice velocities
    !!----------------------------------------------------------------------
    USE par_oce        ! ocean parameter
@@ -34,11 +33,11 @@
    PUBLIC   lim_rhg_2 ! routine called by lim_dyn
 
-   REAL(wp) ::   rzero   = 0.e0   ! constant value: zero
-   REAL(wp) ::   rone    = 1.e0   !            and  one
+   REAL(wp) ::   rzero = 0.e0   ! constant value: zero
+   REAL(wp) ::   rone  = 1.e0   !            and  one
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -56,6 +55,5 @@
       !!  viscous-plastic law including shear strength and a bulk rheology.
       !!
-      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity defined
-      !!              at I-point
+      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity defined at I-point
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
@@ -582,8 +580,4 @@
    !!   Default option          Dummy module       NO 2.0 LIM sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_rhg_2( k1 , k2 )       ! Dummy routine
-      WRITE(*,*) 'lim_rhg_2: You should not have seen this print! error?', k1, k2
-   END SUBROUTINE lim_rhg_2
 #endif
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limrst_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limrst_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limrst_2.F90	(revision 1855)
@@ -2,13 +2,12 @@
    !!======================================================================
    !!                     ***  MODULE  limrst_2  ***
-   !! Ice restart :  write the ice restart file
+   !! LIM-2 ice model : open/read/write the restart file
    !!======================================================================
-   !! History :  2.0  !  01-04  (C. Ethe, G. Madec)  Original code
-   !!                 !  06-07  (S. Masson)  use IOM for restart read/write
+   !! History :  2.0  ! 2001-04  (C. Ethe, G. Madec)  Original code
+   !!                 ! 2006-07  (S. Masson)  use IOM for restart read/write
    !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   lim_rst_opn_2   : open ice restart file
@@ -20,5 +19,4 @@
    USE sbc_oce
    USE sbc_ice
-
    USE in_out_manager
    USE iom
@@ -35,5 +33,5 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -84,4 +82,5 @@
    END SUBROUTINE lim_rst_opn_2
 
+
    SUBROUTINE lim_rst_write_2( kt )
       !!----------------------------------------------------------------------
@@ -91,5 +90,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! number of iteration
-      !
+      !!
       INTEGER ::   iter   ! kt + nn_fsbc -1
       !!----------------------------------------------------------------------
@@ -181,5 +180,5 @@
       ENDIF
 
-      IF ( jprstlib == jprstdimg ) THEN
+      IF( jprstlib == jprstdimg ) THEN
         ! eventually read netcdf file (monobloc)  for restarting on different number of processors
         ! if {cn_icerst_in}.nc exists, then set jlibalt to jpnf90
@@ -188,5 +187,5 @@
       ENDIF
 
-      CALL iom_open ( cn_icerst_in, numrir, kiolib = jlibalt )
+      CALL iom_open( cn_icerst_in, numrir, kiolib = jlibalt )
 
       CALL iom_get( numrir, 'kt_ice' , ziter )
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 1855)
@@ -2,14 +2,13 @@
    !!======================================================================
    !!                       ***  MODULE limsbc_2   ***
-   !!           computation of the flux at the sea ice/ocean interface
+   !!   LIM 2 ice model :   heat, salt, mass and momentum fluxes at the sea ice/ocean interface
    !!======================================================================
-   !! History : 00-01 (H. Goosse) Original code
-   !!           02-07 (C. Ethe, G. Madec) re-writing F90
-   !!           06-07 (G. Madec) surface module
+   !! History :  LIM  ! 2000-01 (H. Goosse) Original code
+   !!            2.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90
+   !!             -   ! 2006-07 (G. Madec) surface module
    !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   lim_sbc_2  : flux at the ice / ocean interface
@@ -20,5 +19,5 @@
    USE sbc_oce          ! surface boundary condition
    USE phycst           ! physical constants
-   USE ice_2            ! LIM sea-ice variables
+   USE ice_2            ! LIM 2 sea-ice variables
 
    USE lbclnk           ! ocean lateral boundary condition
@@ -33,5 +32,5 @@
    PRIVATE
 
-   PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2
+   PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2
 
    REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
@@ -44,5 +43,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -109,8 +108,5 @@
          sice_r(:,:) = sice
          !
-         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN
-            !                                        ! =======================
-            !                                        !  ORCA_R2 configuration
-            !                                        ! =======================
+         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN             !  ORCA_R2 configuration
             ii0 = 145   ;   ii1 = 180        ! Baltic Sea
             ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0
@@ -194,7 +190,7 @@
       END DO
 
-      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )      
-      CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )      
-      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) )
+      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:)                         )      
+      CALL iom_put( 'qns_io_cea'  , qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )      
+      CALL iom_put( 'qsr_io_cea'  , fstric(:,:) * (1. - pfrld(:,:))      )
 
       !------------------------------------------!
@@ -302,8 +298,7 @@
       !-----------------------------------------------!
 
-      IF ( lk_cpl ) THEN           
-         ! Ice surface temperature 
+      IF( lk_cpl ) THEN           
          tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
-         ! Computation of snow/ice and ocean albedo
+         !                                  ! snow/ice and ocean albedos
          CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb )
          alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys)
@@ -312,12 +307,12 @@
 
       IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
-         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
-         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
-            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
-         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
+         CALL prt_ctl( tab2d_1=qsr , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns  , clinfo2=' qns     : ' )
+         CALL prt_ctl( tab2d_1=emp , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps , clinfo2=' emps    : ' )
+         CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
+            &          tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
+         CALL prt_ctl( tab2d_1=fr_i, clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ' )
       ENDIF 
-   
-    END SUBROUTINE lim_sbc_2
+      !
+   END SUBROUTINE lim_sbc_2
 
 #else
@@ -325,7 +320,4 @@
    !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_sbc_2         ! Dummy routine
-   END SUBROUTINE lim_sbc_2
 #endif 
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limtab_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limtab_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limtab_2.F90	(revision 1855)
@@ -2,12 +2,13 @@
    !!======================================================================
    !!                       ***  MODULE limtab_2   ***
-   !!              transform 1D (2D) array to a 2D (1D) table
+   !!   LIM 2 ice model : transform 1D (2D) array to a 2D (1D) array
    !!======================================================================
 #if defined key_lim2
    !!----------------------------------------------------------------------
-   !!   tab_2d_1d  : 2-D to 1-D
-   !!   tab_1d_2d  : 1-D to 2-D
+   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !! * Modules used
+   !!   tab_2d_1d_2  : 2-D to 1-D
+   !!   tab_1d_2d_2  : 1-D to 2-D
+   !!----------------------------------------------------------------------
    USE par_kind
 
@@ -15,10 +16,9 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC tab_2d_1d_2  ! called by lim_ther
-   PUBLIC tab_1d_2d_2  ! called by lim_ther
+   PUBLIC   tab_2d_1d_2   ! called by lim_thd_2
+   PUBLIC   tab_1d_2d_2   ! called by lim_thd_2
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -26,55 +26,49 @@
 CONTAINS
 
-   SUBROUTINE tab_2d_1d_2 ( ndim1d, tab1d, tab2d, ndim2d_x, ndim2d_y, tab_ind )
-
-      INTEGER, INTENT(in) :: &
-         ndim1d, ndim2d_x, ndim2d_y
-
-      REAL(wp), DIMENSION (ndim2d_x, ndim2d_y), INTENT(in) ::  &
-         tab2d
-
-      INTEGER, DIMENSION ( ndim1d), INTENT ( in) :: &
-         tab_ind
-
-      REAL(wp), DIMENSION(ndim1d), INTENT ( out) ::  & 
-         tab1d
-
-      INTEGER ::  &
-         jn , jid, jjd
-        
-      DO jn = 1, ndim1d
-         jid        = MOD( tab_ind(jn) - 1, ndim2d_x ) + 1
-         jjd        = ( tab_ind(jn) - 1 ) / ndim2d_x + 1
-         tab1d( jn) = tab2d( jid, jjd)
+   SUBROUTINE tab_2d_1d_2( kdim1d, ptab1d, ptab2d, kdim2d_i, kdim2d_j, ptab_ind )
+      !!-------------------------------------------------------------------
+      INTEGER , INTENT(in   )                                ::   kdim1d
+      INTEGER , INTENT(in   )                                ::   kdim2d_i, kdim2d_j
+      REAL(wp), INTENT(in   ), DIMENSION(kdim2d_i, kdim2d_j) ::   ptab2d
+      INTEGER , INTENT(in   ), DIMENSION(kdim1d)             ::   ptab_ind
+      REAL(wp), INTENT(  out), DIMENSION(kdim1d)             ::   ptab1d
+      !!
+      INTEGER ::   jn , jid, jjd   ! dummy loop indices
+      !!-------------------------------------------------------------------
+      ! 
+      DO jn = 1, kdim1d
+         jid = MOD( ptab_ind(jn) - 1, kdim2d_i ) + 1
+         jjd = ( ptab_ind(jn) - 1 ) / kdim2d_i + 1
+         ptab1d(jn) = ptab2d(jid,jjd)
       END DO 
-
+      !
    END SUBROUTINE tab_2d_1d_2
 
 
-   SUBROUTINE tab_1d_2d_2 ( ndim1d, tab2d, tab_ind, tab1d, ndim2d_x, ndim2d_y )
-
-      INTEGER, INTENT ( in) :: &
-         ndim1d, ndim2d_x, ndim2d_y
-
-      INTEGER, DIMENSION (ndim1d) , INTENT (in) :: &
-         tab_ind
-
-      REAL(wp), DIMENSION(ndim1d), INTENT (in) ::  &
-         tab1d  
-
-      REAL(wp), DIMENSION (ndim2d_x, ndim2d_y), INTENT ( out) :: &
-         tab2d
-
-      INTEGER :: &
-         jn, jid, jjd
-
-      DO jn = 1, ndim1d
-         jid             = MOD( tab_ind(jn) - 1, ndim2d_x) + 1
-         jjd             =    ( tab_ind(jn) - 1 ) / ndim2d_x  + 1
-         tab2d(jid, jjd) = tab1d( jn)
+   SUBROUTINE tab_1d_2d_2( kdim1d, ptab2d, ptab_ind, ptab1d, kdim2d_i, kdim2d_j )
+      !!-------------------------------------------------------------------
+      INTEGER , INTENT(in   )                                ::   kdim1d
+      INTEGER , INTENT(in   )                                ::   kdim2d_i, kdim2d_j
+      INTEGER , INTENT(in   ), DIMENSION(kdim1d)             ::   ptab_ind
+      REAL(wp), INTENT(in   ), DIMENSION(kdim1d)             ::   ptab1d
+      REAL(wp), INTENT(  out), DIMENSION(kdim2d_i, kdim2d_j) ::   ptab2d
+      !!
+      INTEGER ::   jn, jid, jjd   ! dummy loop indices
+      !!-------------------------------------------------------------------
+      !
+      DO jn = 1, kdim1d
+         jid = MOD( ptab_ind(jn) - 1, kdim2d_i) + 1
+         jjd =    ( ptab_ind(jn) - 1 ) / kdim2d_i  + 1
+         ptab2d(jid,jjd) = ptab1d(jn)
       END DO
-
+      !
    END SUBROUTINE tab_1d_2d_2
 
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
 #endif
+
+   !!======================================================================
 END MODULE limtab_2
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90	(revision 1855)
@@ -2,10 +2,10 @@
    !!======================================================================
    !!                  ***  MODULE limthd_2   ***
-   !!              LIM thermo ice model : ice thermodynamic
+   !!   LIM 2 ice model : thermodynamics
    !!======================================================================
-   !! History :  1.0  ! 2000-01 (LIM)
+   !! History :  LIM  ! 2000-01 (UCL) Original code
    !!            2.0  ! 2002-07 (C. Ethe, G. Madec) F90
-   !!            2.0  ! 2003-08 (C. Ethe)  add lim_thd_init
-   !!             -   ! 2008-2008  (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface
+   !!             -   ! 2003-08 (C. Ethe)  add lim_thd_init
+   !!             -   ! 2008-08  (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface
    !!---------------------------------------------------------------------
 #if defined key_lim2
@@ -13,5 +13,5 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!   lim_thd_2      : thermodynamic of sea ice
+   !!   lim_thd_2      : thermodynamics of sea ice
    !!   lim_thd_init_2 : initialisation of sea-ice thermodynamic
    !!----------------------------------------------------------------------
@@ -23,9 +23,8 @@
    USE lib_mpp
    USE iom             ! IOM library
-   USE ice_2           ! LIM sea-ice variables
+   USE ice_2           ! LIM 2 sea-ice variables
    USE sbc_oce         ! 
    USE sbc_ice         ! 
-   USE thd_ice_2       ! LIM thermodynamic sea-ice variables
-   USE dom_ice_2       ! LIM sea-ice domain
+   USE dom_ice_2       ! LIM 2 sea-ice domain
    USE limthd_zdf_2
    USE limthd_lac_2
@@ -49,6 +48,6 @@
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
-   !!-------- -------------------------------------------------------------
-   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2009) 
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -566,7 +565,4 @@
    !!   Default option          Dummy module       NO LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_thd_2         ! Dummy routine
-   END SUBROUTINE lim_thd_2
 #endif
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_lac_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 1855)
@@ -1,34 +1,35 @@
 MODULE limthd_lac_2
-#if defined key_lim2
    !!======================================================================
    !!                       ***  MODULE limthd_lac_2   ***
-   !!                lateral thermodynamic growth of the ice 
+   !!   LIM 2 ice model :   thermodynamics -- lateral accretion of the ice 
    !!======================================================================
-
-   !!----------------------------------------------------------------------
-   !!   lim_lat_acr_2    : lateral accretion of ice
-   !! * Modules used
+   !! History :  LIM  !  2001-04 (UCL)  original code
+   !!            2.0  !  2002-08 (C. Ethe, G. Madec)  F90, mpp
+   !!----------------------------------------------------------------------
+#if defined key_lim2
+   !!----------------------------------------------------------------------
+   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_lat_acr_2   : lateral accretion of ice
+   !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
-   USE phycst
-   USE thd_ice_2
-   USE ice_2
-   USE limistate_2 
+   USE phycst           ! physical constants
+   USE ice_2            ! LIM 2 sea-ice variables
+   USE limistate_2      ! LIM 2 initial state
      
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_thd_lac_2   ! called by lim_thd_2
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      epsi20 = 1.e-20  ,  &
-      epsi13 = 1.e-13  ,  &
-      zzero  = 0.e0    ,  &
-      zone   = 1.e0
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   PUBLIC   lim_thd_lac_2   ! called by lim_thd_2
+
+   REAL(wp) ::   epsi20 = 1.e-20   ! constant values
+   REAL(wp) ::   epsi13 = 1.e-13   !
+   REAL(wp) ::   rzero  = 0.e0     !
+   REAL(wp) ::   rone   = 1.e0     !
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -60,50 +61,37 @@
       !!               update h_snow_1d, h_ice_1d and tbif_1d(:,:)      
       !! 
-      !! ** References :
-      !!      M. Maqueda, 1995, PhD Thesis, Univesidad Complutense Madrid
-      !!      Fichefet T. and M. Maqueda 1997, J. Geo. Res., 102(C6), 
-      !!                                                12609 -12646   
-      !! History :
-      !!   1.0  !  01-04 (LIM)  original code
-      !!   2.0  !  02-08 (C. Ethe, G. Madec)  F90, mpp
+      !! References :   M. Maqueda, 1995, PhD Thesis, Univesidad Complutense Madrid
+      !!                Fichefet T. and M. Maqueda 1997, J. Geo. Res., 102(C6), 12609 -12646   
       !!-------------------------------------------------------------------
-      !! * Arguments
-      INTEGER , INTENT(IN)::  &
-         kideb          ,   &  ! start point on which the the computation is applied
-         kiut                  ! end point on which the the computation is applied
-
-      !! * Local variables
-      INTEGER ::            &
-         ji             ,   &  !  dummy loop indices
-         iicefr         ,   &  !  1 = existing ice ; 0 = no ice
-         iiceform       ,   &  !  1 = ice formed   ; 0 = no ice formed
-         ihemis                !  dummy indice
-      REAL(wp), DIMENSION(jpij) :: &
-         zqbgow           ,  &  !  heat budget of the open water (negative)
-         zfrl_old         ,  &  !  previous sea/ice fraction
-         zhice_old        ,  &  !  previous ice thickness
-         zhice0           ,  &  !  thickness of newly formed ice in leads
-         zfrlmin          ,  &  !  minimum fraction for leads
-         zdhicbot               !  part of thickness of newly formed ice in leads which 
-                                !  has been already used in transport for example
-      REAL(wp)  ::  &
-         zhemis           ,  &  !  hemisphere (0 = North, 1 = South)
-         zhicenew         ,  &  !  new ice thickness
-         zholds2          ,  &  !  ratio of previous ice thickness and 2 
-         zhnews2          ,  &  !  ratio of new ice thickness and 2 
-         zfrlnew          ,  &  !  new sea/ice fraction
-         zfrld            ,  &  !  ratio of sea/ice fraction and minimum fraction for leads
-         zfrrate          ,  &  !  leads-closure rate
-         zdfrl                  !  sea-ice fraction increment
-      REAL(wp)  ::  &
-         zdh1 , zdh2 , zdh3 , zdh4, zdh5   , &   ! tempory scalars
-         ztint , zta1 , zta2 , zta3 , zta4 , &
-         zah, zalpha , zbeta
+      INTEGER, INTENT(in) ::   kideb   ! start point on which the the computation is applied
+      INTEGER, INTENT(in) ::   kiut    ! end point on which the the computation is applied
+      !!
+      INTEGER ::   ji         !  dummy loop indices
+      INTEGER ::   iicefr     !  1 = existing ice ; 0 = no ice
+      INTEGER ::   iiceform   !  1 = ice formed   ; 0 = no ice formed
+      INTEGER ::   ihemis     !  hemisphere indicator 
+      REAL(wp), DIMENSION(jpij) ::   zqbgow      !  heat budget of the open water (negative)
+      REAL(wp), DIMENSION(jpij) ::   zfrl_old    !  previous sea/ice fraction
+      REAL(wp), DIMENSION(jpij) ::   zhice_old   !  previous ice thickness
+      REAL(wp), DIMENSION(jpij) ::   zhice0      !  thickness of newly formed ice in leads
+      REAL(wp), DIMENSION(jpij) ::   zfrlmin     !  minimum fraction for leads
+      REAL(wp), DIMENSION(jpij) ::   zdhicbot    !  part of thickness of newly formed ice in leads which 
+      !                                          !  has been already used in transport for example
+      REAL(wp) ::   zhemis     !  hemisphere (0 = North, 1 = South)
+      REAL(wp) ::   zhicenew   !  new ice thickness
+      REAL(wp) ::   zholds2    !  ratio of previous ice thickness and 2 
+      REAL(wp) ::   zhnews2    !  ratio of new ice thickness and 2 
+      REAL(wp) ::   zfrlnew    !  new sea/ice fraction
+      REAL(wp) ::   zfrld      !  ratio of sea/ice fraction and minimum fraction for leads
+      REAL(wp) ::   zfrrate    !  leads-closure rate
+      REAL(wp) ::   zdfrl      !  sea-ice fraction increment
+      REAL(wp) ::   zdh1 , zdh2, zdh3, zdh4, zdh5   ! tempory scalars
+      REAL(wp) ::   ztint, zta1, zta2, zta3, zta4   !    -       -
+      REAL(wp) ::   zah, zalpha, zbeta              !    -       -
       !!---------------------------------------------------------------------      
                    
       !--------------------------------------------------------------
       !   Computation of the heat budget of the open water (negative)
-      !--------------------------------------------------------------
-      
+      !-------------------------------------------------------------- 
       DO ji = kideb , kiut      
          zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)
@@ -131,5 +119,4 @@
       !     Morales Maqueda, 1995 - Fichefet and Morales Maqueda, 1997
       !---------------------------------------------------------------------
-      
 !CDIR NOVERRCHK
       DO ji = kideb , kiut
@@ -145,5 +132,5 @@
          !--computation of the remaining part of ice thickness which has been already used
          zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) ) & 
-                      -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic ) 
+            &         -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic ) 
       END DO
  
@@ -155,5 +142,4 @@
       !         (1-Anew) * hsnew = (1-Aold) * hsold            
       !--------------------------------------------------------------------------------------------
-      
       DO ji = kideb , kiut
          iicefr       = 1 - MAX( 0, INT( SIGN( 1.5 * zone , zfrl_old(ji) - 1.0 + epsi13 ) ) )
@@ -169,5 +155,4 @@
       !   Ajustement of ice internal temperatures
       !-------------------------------------------------------
-      
       DO ji = kideb , kiut
          iicefr      = 1 - MAX( 0, INT( SIGN( 1.5 * zone , zfrl_old(ji) - 1.0 + epsi13 ) ) )
@@ -213,19 +198,17 @@
       !           dV = Vnew - Vold
       !-------------------------------------------------------------
-      
       DO ji = kideb , kiut
          dvlbq_1d  (ji) = ( 1. - frld_1d(ji) ) * h_ice_1d(ji) - ( 1. - zfrl_old(ji) ) * zhice_old(ji)
          rdmicif_1d(ji) = rdmicif_1d(ji) + rhoic * dvlbq_1d(ji)
       END DO
-      
+      !
    END SUBROUTINE lim_thd_lac_2
+   
 #else
+   !!----------------------------------------------------------------------
+   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
    !!======================================================================
-   !!                       ***  MODULE limthd_lac_2   ***
-   !!                           no sea ice model
-   !!======================================================================
-CONTAINS
-   SUBROUTINE lim_thd_lac_2           ! Empty routine
-   END SUBROUTINE lim_thd_lac_2
-#endif
 END MODULE limthd_lac_2
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 1855)
@@ -4,6 +4,6 @@
    !!                thermodynamic growth and decay of the ice 
    !!======================================================================
-   !! History :  1.0  !  01-04 (LIM) Original code
-   !!            2.0  !  02-08 (C. Ethe, G. Madec) F90
+   !! History :  1.0  !  2001-04 (LIM) Original code
+   !!            2.0  !  2002-08 (C. Ethe, G. Madec) F90
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -11,14 +11,12 @@
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
    !!   lim_thd_zdf_2 : vertical accr./abl. and lateral ablation of sea ice
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_oce          ! ocean parameters
-   USE phycst           ! ???
-   USE thd_ice_2
-   USE ice_2
-   USE limistate_2
-   USE in_out_manager
+   USE phycst           ! physical constants
+   USE ice_2            ! LIM 2 sea-ice variables
+   USE limistate_2      ! LIM 2 initial state
+   USE in_out_manager   ! I/O manager
+   !!
    USE cpl_oasis3, ONLY : lk_cpl
       
@@ -26,12 +24,13 @@
    PRIVATE
 
-   PUBLIC   lim_thd_zdf_2        ! called by lim_thd_2
-
-   REAL(wp) ::   epsi20 = 1.e-20  ,  &  ! constant values
-      &          epsi13 = 1.e-13  ,  &
-      &          zzero  = 0.e0    ,  &
-      &          zone   = 1.e0
+   PUBLIC   lim_thd_zdf_2   ! called by lim_thd_2
+
+   REAL(wp) ::   epsi20 = 1.e-20   ! constant values
+   REAL(wp) ::   epsi13 = 1.e-13   !
+   REAL(wp) ::   rzero  = 0.e0     !
+   REAL(wp) ::   rone   = 1.e0     !
+
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -73,328 +72,322 @@
       !!
       INTEGER ::   ji       ! dummy loop indices
-      REAL(wp), DIMENSION(jpij,2) ::   zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 )
-      REAL(wp), DIMENSION(jpij) ::  &
-         ztsmlt      &    ! snow/ice surface melting temperature
-         ,ztbif      &    ! int. temp. at the mid-point of the 1st layer of the snow/ice sys. 
-         ,zksn       &    ! effective conductivity of snow
-         ,zkic       &    ! effective conductivity of ice
-         ,zksndh     &    ! thermal cond. at the mid-point of the 1st layer of the snow/ice sys. 
-         , zfcsu     &    ! conductive heat flux at the surface of the snow/ice system 
-         , zfcsudt   &    ! = zfcsu * dt
-         , zi0       &    ! frac. of the net SW rad. which is not absorbed at the surface
-         , z1mi0     &    ! fraction of the net SW radiation absorbed at the surface
-         , zqmax     &    ! maximum energy stored in brine pockets
-         , zrcpdt    &    ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer)
-         , zts_old   &    ! previous surface temperature
-         , zidsn , z1midsn , zidsnic ! tempory variables
-      REAL(wp), DIMENSION(jpij) ::   &
-          zfnet       &  ! net heat flux at the top surface( incl. conductive heat flux)
-          , zsprecip  &    ! snow accumulation
-          , zhsnw_old &    ! previous snow thickness
-          , zdhictop  &    ! change in ice thickness due to top surf ablation/accretion
-          , zdhicbot  &    ! change in ice thickness due to bottom surf abl/acc
-          , zqsup     &    ! energy transmitted to ocean (coming from top surf abl/acc)
-          , zqocea    &    ! energy transmitted to ocean (coming from bottom sur abl/acc)
-          , zfrl_old  &    ! previous sea/ice fraction
-          , zfrld_1d    &    ! new sea/ice fraction
-          , zep            ! internal temperature of the 2nd layer of the snow/ice system
-       REAL(wp), DIMENSION(3) :: & 
-          zplediag  &    ! principle diagonal, subdiag. and supdiag. of the 
-          , zsubdiag  &    ! tri-diagonal matrix coming from the computation
-          , zsupdiag  &    ! of the temperatures inside the snow-ice system
-          , zsmbr          ! second member
-       REAL(wp) :: & 
-          zhsu     &     ! thickness of surface layer
-          , zhe      &     ! effective thickness for compu. of equ. thermal conductivity
-          , zheshth  &     ! = zhe / thth
-          , zghe     &     ! correction factor of the thermal conductivity
-          , zumsb    &     ! parameter for numerical method to solve heat-diffusion eq.
-          , zkhsn    &     ! conductivity at the snow layer
-          , zkhic    &     ! conductivity at the ice layers
-          , zkint    &     ! equivalent conductivity at the snow-ice interface
-          , zkhsnint &     ! = zkint*dt / (hsn*rhosn*cpsn)  
-          , zkhicint &     ! = 2*zkint*dt / (hic*rhoic*cpic)
-          , zpiv1 , zpiv2  &       ! tempory scalars used to solve the tri-diagonal system
-          , zb2 , zd2 , zb3 , zd3 &
-          , ztint          ! equivalent temperature at the snow-ice interface
-       REAL(wp) :: & 
-          zexp      &     ! exponential function of the ice thickness
-          , zfsab     &     ! part of solar radiation stored in brine pockets
-          , zfts      &     ! value of energy balance function when the temp. equal surf. temp.
-          , zdfts     &     ! value of derivative of ztfs when the temp. equal surf. temp.
-          , zdts      &     ! surface temperature increment
-          , zqsnw_mlt &     ! energy needed to melt snow
-          , zdhsmlt   &     ! change in snow thickness due to melt
-          , zhsn      &     ! snow thickness (previous+accumulation-melt)
-          , zqsn_mlt_rem &  ! remaining heat coming from snow melting
-          , zqice_top_mlt & ! energy used to melt ice at top surface
-          , zdhssub      &  ! change in snow thick. due to sublimation or evaporation
-          , zdhisub      &  ! change in ice thick. due to sublimation or evaporation    
-          , zdhsn        &  ! snow ice thickness increment
-          , zdtsn        &  ! snow internal temp. increment
-          , zdtic        &  ! ice internal temp. increment
-          , zqnes          ! conductive energy due to ice melting in the first ice layer
-       REAL(wp) :: & 
-          ztbot     &      ! temperature at the bottom surface
-          , zfcbot    &      ! conductive heat flux at bottom surface
-          , zqice_bot &      ! energy used for bottom melting/growing
-          , zqice_bot_mlt &  ! energy used for bottom melting
-          , zqstbif_bot  &  ! part of energy stored in brine pockets used for bottom melting
-          , zqstbif_old  &  ! tempory var. for zqstbif_bot
-          , zdhicmlt      &  ! change in ice thickness due to bottom melting
-          , zdhicm        &  ! change in ice thickness var. 
-          , zdhsnm        &  ! change in snow thickness var. 
-          , zhsnfi        &  ! snow thickness var. 
-          , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables
-          , ztb2, ztb3
-       REAL(wp) :: & 
-          zdrmh         &   ! change in snow/ice thick. after snow-ice formation
-          , zhicnew       &   ! new ice thickness
-          , zhsnnew       &   ! new snow thickness
-          , zquot , ztneq &   ! tempory temp. variables
-          , zqice, zqicetot & ! total heat inside the snow/ice system
-          , zdfrl         &   ! change in ice concentration
-          , zdvsnvol      &   ! change in snow volume
-          , zdrfrl1, zdrfrl2 &  ! tempory scalars
-          , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind
-       !!----------------------------------------------------------------------
-
-       !-----------------------------------------------------------------------
-       !  1. Boundaries conditions for snow/ice system internal temperature
-       !       - If tbif_1d(ji,1) > rt0_snow, tbif_1d(ji,1) = rt0_snow 
-       !       - If tbif_1d(ji,2/3) > rt0_ice, tbif_1d(ji,2/3) = rt0_ice 
-       !     Computation of energies due to surface and bottom melting 
-       !-----------------------------------------------------------------------
-       
-       DO ji = kideb , kiut
-          zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) )
-          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) )
-          !--computation of energy due to surface melting
-          zqcmlt(ji,1) = ( MAX ( zzero ,  &
-             &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
-          !--computation of energy due to bottom melting
-          zqcmlt(ji,2) = ( MAX( zzero , &
-             &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
-             &           + MAX( zzero , &
-             &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
-             &           ) * ( 1.0 - zihic  )
-          !--limitation of  snow/ice system internal temperature
-          tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) )
-          tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) )
-          tbif_1d(ji,3)   = MIN( rt0_ice , tbif_1d(ji,3) )
-       END DO
-
-       !-------------------------------------------
-       !  2. Calculate some intermediate variables.  
-       !-------------------------------------------
-       
-       ! initialisation of the thickness of surface layer
-       zhsu = hnzst  
-
-       DO ji = kideb , kiut
-          zind   = MAX( zzero , SIGN( zone , zhsu - h_snow_1d(ji) ) )
-          zihsn  = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) )
-          zihsn  = MAX( zihsn , zind )
-          zihic  = MAX( zzero , sign( zone , hicdif - h_ice_1d(ji) ) )
-          !     2.1. Computation of surface melting temperature
-          !----------------------------------------------------
-          zind  = MAX( zzero , SIGN( zone , -h_snow_1d(ji) ) )
-          ztsmlt(ji) = ( 1.0 - zind ) * rt0_snow + zind * rt0_ice
-          !
-          !     2.2. Effective conductivity of snow and ice
-          !-----------------------------------------------
-
-          !---computation of the correction factor on the thermal conductivity
-          !-- (Morales Maqueda, 1995 ; Fichefet and Morales Maqueda, 1997)
-          zhe      =  ( rcdsn / ( rcdsn + rcdic ) ) * h_ice_1d(ji)   &
-             &     + ( rcdic / ( rcdsn + rcdic ) ) * h_snow_1d(ji) 
-          zihe     = MAX( zzero , SIGN( zone , 2.0 * zhe - thth ) )
-          zheshth  = zhe / thth
-          zghe     = ( 1.0 - zihe ) * zheshth * ( 2.0 - zheshth )   &
-             &     +         zihe   * 0.5 * ( 1.5 + LOG( 2.0 * zheshth ) )
-
-          !---effective conductivities 
-          zksn(ji)  = zghe * rcdsn  
-          zkic(ji)  = zghe * rcdic
-
-          !
-          !     2.3. Computation of the conductive heat flux from the snow/ice 
-          !          system interior toward the top surface
-          !------------------------------------------------------------------
-
-          !---Thermal conductivity at the mid-point of the first snow/ice system layer
-          zksndh(ji) =   ( ( 1.0 - zihsn ) * 2.0 * zksn(ji) + zihsn * 4.0 * zkic(ji) )   &
-             &         / ( ( 1.0 - zihsn ) *  h_snow_1d(ji)                              &
-             &           +        zihsn   *  ( ( 1.0 + 3.0 * zihic ) * h_ice_1d(ji)      &
-             &           + 4.0 * zkic(ji)/zksn(ji) * h_snow_1d(ji) ) )
-
-          !---internal temperature at the mid-point of the first snow/ice system layer
-          ztbif(ji)  = ( 1.0 - zihsn ) * tbif_1d(ji,1)                       &
-             &       +         zihsn   * ( ( 1.0 - zihic ) * tbif_1d(ji,2)   &
-             &       +         zihic   * tfu_1d(ji)   )
-          !---conductive heat flux 
-          zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
-
-       END DO
-
-       !--------------------------------------------------------------------
-       !  3. Calculate : 
-       !     - fstbif_1d, part of solar radiation absorbing inside the ice
-       !       assuming an exponential absorption (Grenfell and Maykut, 1977)
-       !     - zqmax,  maximum energy stored in brine pockets
-       !     - qstbif_1d, total energy stored in brine pockets (updating)
-       !-------------------------------------------------------------------
-
-       DO ji = kideb , kiut
-          zihsn  = MAX( zzero , SIGN (zone , -h_snow_1d(ji) ) )
-          zihic  = MAX( zzero , 1.0 - ( h_ice_1d(ji) / zhsu ) )     
-          zind   = MAX( zzero , SIGN (zone , hicdif - h_ice_1d(ji) ) )
-          !--Computation of the fraction of the net shortwave radiation which
-          !--penetrates inside the ice cover ( See Forcat)
-          zi0(ji)  = zihsn * ( fr1_i0_1d(ji) + zihic * fr2_i0_1d(ji) )
-          zexp     = MIN( zone , EXP( -1.5 * ( h_ice_1d(ji) - zhsu ) ) )
-          fstbif_1d(ji) = zi0(ji) * qsr_ice_1d(ji) * zexp
-          !--Computation of maximum energy stored in brine pockets zqmax and update
-          !--the total energy stored in brine pockets, if less than zqmax
-          zqmax(ji) = MAX( zzero , 0.5 * xlic * ( h_ice_1d(ji) - hicmin ) )
-          zfsab   = zi0(ji) * qsr_ice_1d(ji) * ( 1.0 - zexp )
-          zihq    = ( 1.0 - zind ) * MAX(zzero, SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) ) &
-             &    +         zind   * zone
-          qstbif_1d(ji) = ( qstbif_1d(ji) + ( 1.0 - zihq ) * zfsab * rdt_ice ) * swiqst
-          !--fraction of shortwave radiation absorbed at surface
-          ziexp = zihq * zexp + ( 1.0 - zihq ) * ( swiqst + ( 1.0 - swiqst ) * zexp )
-          z1mi0(ji) = 1.0 - zi0(ji) * ziexp
-       END DO
-
-       !--------------------------------------------------------------------------------
-       !  4. Computation of the surface temperature : determined by considering the 
-       !     budget of a thin layer of thick. zhsu at the top surface (H. Grenier, 1995) 
-       !     and based on a surface energy balance : 
-       !     hsu * rcp * dT/dt = Fsr + Fnsr(T) + Fcs(T),
-       !     where - Fsr is the net absorbed solar radiation, 
-       !           - Fnsr is the total non solar radiation (incoming and outgoing long-wave,
-       !             sensible and latent heat fluxes)
-       !           - Fcs the conductive heat flux at the top of surface
-       !------------------------------------------------------------------------------
-
-       !     4.1. Computation of intermediate values
-       !---------------------------------------------
-       DO ji = kideb, kiut
-          zrcpdt(ji) = ( rcpsn * MIN( h_snow_1d(ji) , zhsu )    &
-             &       + rcpic * MAX( zhsu - h_snow_1d(ji) , zzero ) ) / rdt_ice
-          zts_old(ji) =  sist_1d(ji)
-       END DO
-
-       !     4.2. Computation of surface temperature by expanding the eq. of energy balance
-       !          with Ts = Tp + DT. One obtain , F(Tp) + DT * DF(Tp) = 0
-       !          where  - F(Tp) = Fsr + Fnsr(Tp) + Fcs(Tp) 
-       !                 - DF(Tp)= (dFnsr(Tp)/dT) + (dFcs(Tp)/dT) - hsu*rcp/dt
-       !---------------------------------------------------------------------------------
-
-       DO ji = kideb, kiut
-          !---computation of the derivative of energy balance function 
-          zdfts    =  zksndh(ji)   & ! contribution of the conductive heat flux
-             &      + zrcpdt(ji)   & ! contribution of hsu * rcp / dt
-             &      - dqns_ice_1d (ji)     ! contribution of the total non solar radiation 
-          !---computation of the energy balance function 
-          zfts    = - z1mi0 (ji) * qsr_ice_1d(ji)   & ! net absorbed solar radiation
-             &      - qns_ice_1d(ji)                & ! total non solar radiation
-             &      - zfcsu (ji)                      ! conductive heat flux from the surface
-          !---computation of surface temperature increment  
-          zdts    = -zfts / zdfts
-          !---computation of the new surface temperature 
-          sist_1d(ji) = sist_1d(ji) + zdts
-       END DO
-
-       !----------------------------------------------------------------------------
-       !  5. Boundary condition at the top surface
-       !--    IF Tsb < Tmelt, Fnet = Fcs (the net heat flux equal the conductive heat flux)
-       !      Otherwise Tsb = Tmelt and Qnet(Tmelt) > 0 
-       !      Fnet(Tmelt) is therefore the net surface flux needed for melting
-       !----------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpij,2) ::   zqcmlt   ! energy due to surface( /1 ) and bottom melting( /2 )
+      REAL(wp), DIMENSION(jpij) ::   ztsmlt     ! snow/ice surface melting temperature
+      REAL(wp), DIMENSION(jpij) ::   ztbif      ! int. temp. at the mid-point of the 1st layer of the snow/ice sys. 
+      REAL(wp), DIMENSION(jpij) ::   zksn       ! effective conductivity of snow
+      REAL(wp), DIMENSION(jpij) ::   zkic       ! effective conductivity of ice
+      REAL(wp), DIMENSION(jpij) ::   zksndh     ! thermal cond. at the mid-point of the 1st layer of the snow/ice sys. 
+      REAL(wp), DIMENSION(jpij) ::   zfcsu      ! conductive heat flux at the surface of the snow/ice system 
+      REAL(wp), DIMENSION(jpij) ::   zfcsudt    ! = zfcsu * dt
+      REAL(wp), DIMENSION(jpij) ::   zi0        ! frac. of the net SW rad. which is not absorbed at the surface
+      REAL(wp), DIMENSION(jpij) ::   z1mi0      ! fraction of the net SW radiation absorbed at the surface
+      REAL(wp), DIMENSION(jpij) ::   zqmax      ! maximum energy stored in brine pockets
+      REAL(wp), DIMENSION(jpij) ::   zrcpdt     ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer)
+      REAL(wp), DIMENSION(jpij) ::   zts_old    ! previous surface temperature
+      REAL(wp), DIMENSION(jpij) ::   zidsn , z1midsn , zidsnic ! tempory variables
+      !!
+      REAL(wp), DIMENSION(jpij) ::   zfnet       ! net heat flux at the top surface( incl. conductive heat flux)
+      REAL(wp), DIMENSION(jpij) ::   zsprecip    ! snow accumulation
+      REAL(wp), DIMENSION(jpij) ::   zhsnw_old   ! previous snow thickness
+      REAL(wp), DIMENSION(jpij) ::   zdhictop    ! change in ice thickness due to top surf ablation/accretion
+      REAL(wp), DIMENSION(jpij) ::   zdhicbot    ! change in ice thickness due to bottom surf abl/acc
+      REAL(wp), DIMENSION(jpij) ::   zqsup       ! energy transmitted to ocean (coming from top surf abl/acc)
+      REAL(wp), DIMENSION(jpij) ::   zqocea      ! energy transmitted to ocean (coming from bottom sur abl/acc)
+      REAL(wp), DIMENSION(jpij) ::   zfrl_old    ! previous sea/ice fraction
+      REAL(wp), DIMENSION(jpij) ::   zfrld_1d    ! new sea/ice fraction
+      REAL(wp), DIMENSION(jpij) ::   zep         ! internal temperature of the 2nd layer of the snow/ice system
+      !!
+      REAL(wp), DIMENSION(3) ::   zplediag       ! principle diagonal, subdiag. and supdiag. of the 
+      REAL(wp), DIMENSION(3) ::   zsubdiag       ! tri-diagonal matrix coming from the computation
+      REAL(wp), DIMENSION(3) ::   zsupdiag       ! of the temperatures inside the snow-ice system
+      REAL(wp), DIMENSION(3) ::   zsmbr          ! second member
+      REAL(wp) ::   zhsu        ! thickness of surface layer
+      REAL(wp) ::   zhe         ! effective thickness for compu. of equ. thermal conductivity
+      REAL(wp) ::   zheshth     ! = zhe / thth
+      REAL(wp) ::   zghe        ! correction factor of the thermal conductivity
+      REAL(wp) ::   zumsb       ! parameter for numerical method to solve heat-diffusion eq.
+      REAL(wp) ::   zkhsn       ! conductivity at the snow layer
+      REAL(wp) ::   zkhic       ! conductivity at the ice layers
+      REAL(wp) ::   zkint       ! equivalent conductivity at the snow-ice interface
+      REAL(wp) ::   zkhsnint    ! = zkint*dt / (hsn*rhosn*cpsn)  
+      REAL(wp) ::   zkhicint    ! = 2*zkint*dt / (hic*rhoic*cpic)
+      REAL(wp) ::   zpiv1 , zpiv2        ! tempory scalars used to solve the tri-diagonal system
+      REAL(wp) ::   zb2, zd2, zb3, zd3
+      REAL(wp) ::   ztint          ! equivalent temperature at the snow-ice interface
+      REAL(wp) ::   zexp         ! exponential function of the ice thickness
+      REAL(wp) ::   zfsab        ! part of solar radiation stored in brine pockets
+      REAL(wp) ::   zfts         ! value of energy balance function when the temp. equal surf. temp.
+      REAL(wp) ::   zdfts        ! value of derivative of ztfs when the temp. equal surf. temp.
+      REAL(wp) ::   zdts         ! surface temperature increment
+      REAL(wp) ::   zqsnw_mlt    ! energy needed to melt snow
+      REAL(wp) ::   zdhsmlt      ! change in snow thickness due to melt
+      REAL(wp) ::   zhsn         ! snow thickness (previous+accumulation-melt)
+      REAL(wp) ::   zqsn_mlt_rem   ! remaining heat coming from snow melting
+      REAL(wp) ::   zqice_top_mlt   ! energy used to melt ice at top surface
+      REAL(wp) ::   zdhssub        ! change in snow thick. due to sublimation or evaporation
+      REAL(wp) ::   zdhisub        ! change in ice thick. due to sublimation or evaporation    
+      REAL(wp) ::   zdhsn          ! snow ice thickness increment
+      REAL(wp) ::   zdtsn          ! snow internal temp. increment
+      REAL(wp) ::   zdtic          ! ice internal temp. increment
+      REAL(wp) ::   zqnes          ! conductive energy due to ice melting in the first ice layer
+      !!
+      REAL(wp) ::   ztbot           ! temperature at the bottom surface
+      REAL(wp) ::   zfcbot          ! conductive heat flux at bottom surface
+      REAL(wp) ::   zqice_bot       ! energy used for bottom melting/growing
+      REAL(wp) ::   zqice_bot_mlt   ! energy used for bottom melting
+      REAL(wp) ::   zqstbif_bot    ! part of energy stored in brine pockets used for bottom melting
+      REAL(wp) ::   zqstbif_old    ! tempory var. for zqstbif_bot
+      REAL(wp) ::   zdhicmlt        ! change in ice thickness due to bottom melting
+      REAL(wp) ::   zdhicm          ! change in ice thickness var. 
+      REAL(wp) ::   zdhsnm          ! change in snow thickness var. 
+      REAL(wp) ::   zhsnfi          ! snow thickness var. 
+      REAL(wp) ::   zc1, zpc1, zc2, zpc2, zp1, zp2   ! temporary scalars
+      REAL(wp) ::   ztb2, ztb3
+      !!
+      REAL(wp) ::   zdrmh           ! change in snow/ice thick. after snow-ice formation
+      REAL(wp) ::   zhicnew         ! new ice thickness
+      REAL(wp) ::   zhsnnew         ! new snow thickness
+      REAL(wp) ::   zquot , ztneq   ! tempory temp. variables
+      REAL(wp) ::   zqice, zqicetot & ! total heat inside the snow/ice system
+      REAL(wp) ::   zdfrl           ! change in ice concentration
+      REAL(wp) ::   zdvsnvol        ! change in snow volume
+      REAL(wp) ::   zdrfrl1, zdrfrl2   ! tempory scalars
+      REAL(wp) ::   zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind
+      !!----------------------------------------------------------------------
+
+      !-----------------------------------------------------------------------
+      !  1. Boundaries conditions for snow/ice system internal temperature
+      !       - If tbif_1d(ji,1) > rt0_snow, tbif_1d(ji,1) = rt0_snow 
+      !       - If tbif_1d(ji,2/3) > rt0_ice, tbif_1d(ji,2/3) = rt0_ice 
+      !     Computation of energies due to surface and bottom melting 
+      !-----------------------------------------------------------------------
+      DO ji = kideb , kiut
+         zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) )
+         zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) )
+         !--computation of energy due to surface melting
+         zqcmlt(ji,1) = ( MAX ( zzero ,  &
+            &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
+         !--computation of energy due to bottom melting
+         zqcmlt(ji,2) = ( MAX( zzero , &
+            &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
+            &           + MAX( zzero , &
+            &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
+            &           ) * ( 1.0 - zihic  )
+         !--limitation of  snow/ice system internal temperature
+         tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) )
+         tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) )
+         tbif_1d(ji,3)   = MIN( rt0_ice , tbif_1d(ji,3) )
+      END DO
+
+      !-------------------------------------------
+      !  2. Calculate some intermediate variables.  
+      !-------------------------------------------
+      zhsu = hnzst         ! initialisation of the thickness of surface layer
+      !
+      DO ji = kideb , kiut
+         zind   = MAX(  zzero , SIGN( zone , zhsu   - h_snow_1d(ji) )  )
+         zihsn  = MAX(  zzero , SIGN( zone , hsndif - h_snow_1d(ji) )  )
+         zihsn  = MAX(  zihsn , zind                                   )
+         zihic  = MAX(  zzero , sign( zone , hicdif - h_ice_1d (ji) )  )
+         !
+         !     2.1. Computation of surface melting temperature
+         !----------------------------------------------------
+         zind  = MAX( zzero , SIGN( zone , -h_snow_1d(ji) ) )
+         ztsmlt(ji) = ( 1.0 - zind ) * rt0_snow + zind * rt0_ice
+         !
+         !     2.2. Effective conductivity of snow and ice
+         !-----------------------------------------------
+
+         !---computation of the correction factor on the thermal conductivity
+         !-- (Morales Maqueda, 1995 ; Fichefet and Morales Maqueda, 1997)
+         zhe      =  ( rcdsn / ( rcdsn + rcdic ) ) * h_ice_1d(ji)   &
+            &     + ( rcdic / ( rcdsn + rcdic ) ) * h_snow_1d(ji) 
+         zihe     = MAX( zzero , SIGN( zone , 2.0 * zhe - thth ) )
+         zheshth  = zhe / thth
+         zghe     = ( 1.0 - zihe ) * zheshth * ( 2.0 - zheshth )   &
+            &     +         zihe   * 0.5 * ( 1.5 + LOG( 2.0 * zheshth ) )
+
+         !---effective conductivities 
+         zksn(ji)  = zghe * rcdsn  
+         zkic(ji)  = zghe * rcdic
+         !
+         !     2.3. Computation of the conductive heat flux from the snow/ice 
+         !          system interior toward the top surface
+         !------------------------------------------------------------------
+
+         !---Thermal conductivity at the mid-point of the first snow/ice system layer
+         zksndh(ji) =   ( ( 1.0 - zihsn ) * 2.0 * zksn(ji) + zihsn * 4.0 * zkic(ji) )   &
+            &         / ( ( 1.0 - zihsn ) *  h_snow_1d(ji)                              &
+            &           +        zihsn   *  ( ( 1.0 + 3.0 * zihic ) * h_ice_1d(ji)      &
+            &           + 4.0 * zkic(ji)/zksn(ji) * h_snow_1d(ji) ) )
+
+         !---internal temperature at the mid-point of the first snow/ice system layer
+         ztbif(ji)  = ( 1.0 - zihsn ) * tbif_1d(ji,1)                       &
+            &       +         zihsn   * ( ( 1.0 - zihic ) * tbif_1d(ji,2)   &
+            &       +         zihic   * tfu_1d(ji)   )
+         !---conductive heat flux 
+         zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
+
+      END DO
+
+      !--------------------------------------------------------------------
+      !  3. Calculate : 
+      !     - fstbif_1d, part of solar radiation absorbing inside the ice
+      !       assuming an exponential absorption (Grenfell and Maykut, 1977)
+      !     - zqmax,  maximum energy stored in brine pockets
+      !     - qstbif_1d, total energy stored in brine pockets (updating)
+      !-------------------------------------------------------------------
+
+      DO ji = kideb , kiut
+         zihsn  = MAX( zzero , SIGN (zone , -h_snow_1d(ji) ) )
+         zihic  = MAX( zzero , 1.0 - ( h_ice_1d(ji) / zhsu ) )     
+         zind   = MAX( zzero , SIGN (zone , hicdif - h_ice_1d(ji) ) )
+         !--Computation of the fraction of the net shortwave radiation which
+         !--penetrates inside the ice cover ( See Forcat)
+         zi0(ji)  = zihsn * ( fr1_i0_1d(ji) + zihic * fr2_i0_1d(ji) )
+         zexp     = MIN( zone , EXP( -1.5 * ( h_ice_1d(ji) - zhsu ) ) )
+         fstbif_1d(ji) = zi0(ji) * qsr_ice_1d(ji) * zexp
+         !--Computation of maximum energy stored in brine pockets zqmax and update
+         !--the total energy stored in brine pockets, if less than zqmax
+         zqmax(ji) = MAX( zzero , 0.5 * xlic * ( h_ice_1d(ji) - hicmin ) )
+         zfsab   = zi0(ji) * qsr_ice_1d(ji) * ( 1.0 - zexp )
+         zihq    = ( 1.0 - zind ) * MAX(zzero, SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) ) &
+            &    +         zind   * zone
+         qstbif_1d(ji) = ( qstbif_1d(ji) + ( 1.0 - zihq ) * zfsab * rdt_ice ) * swiqst
+         !--fraction of shortwave radiation absorbed at surface
+         ziexp = zihq * zexp + ( 1.0 - zihq ) * ( swiqst + ( 1.0 - swiqst ) * zexp )
+         z1mi0(ji) = 1.0 - zi0(ji) * ziexp
+      END DO
+
+      !--------------------------------------------------------------------------------
+      !  4. Computation of the surface temperature : determined by considering the 
+      !     budget of a thin layer of thick. zhsu at the top surface (H. Grenier, 1995) 
+      !     and based on a surface energy balance : 
+      !     hsu * rcp * dT/dt = Fsr + Fnsr(T) + Fcs(T),
+      !     where - Fsr is the net absorbed solar radiation, 
+      !           - Fnsr is the total non solar radiation (incoming and outgoing long-wave,
+      !             sensible and latent heat fluxes)
+      !           - Fcs the conductive heat flux at the top of surface
+      !------------------------------------------------------------------------------
+
+      !     4.1. Computation of intermediate values
+      !---------------------------------------------
+      DO ji = kideb, kiut
+         zrcpdt(ji) = ( rcpsn * MIN( h_snow_1d(ji) , zhsu )    &
+            &       + rcpic * MAX( zhsu - h_snow_1d(ji) , zzero ) ) / rdt_ice
+         zts_old(ji) =  sist_1d(ji)
+      END DO
+
+      !     4.2. Computation of surface temperature by expanding the eq. of energy balance
+      !          with Ts = Tp + DT. One obtain , F(Tp) + DT * DF(Tp) = 0
+      !          where  - F(Tp) = Fsr + Fnsr(Tp) + Fcs(Tp) 
+      !                 - DF(Tp)= (dFnsr(Tp)/dT) + (dFcs(Tp)/dT) - hsu*rcp/dt
+      !---------------------------------------------------------------------------------
+
+      DO ji = kideb, kiut
+         !---computation of the derivative of energy balance function 
+         zdfts    =  zksndh(ji)   & ! contribution of the conductive heat flux
+            &      + zrcpdt(ji)   & ! contribution of hsu * rcp / dt
+            &      - dqns_ice_1d (ji)     ! contribution of the total non solar radiation 
+         !---computation of the energy balance function 
+         zfts    = - z1mi0 (ji) * qsr_ice_1d(ji)   & ! net absorbed solar radiation
+            &      - qns_ice_1d(ji)                & ! total non solar radiation
+            &      - zfcsu (ji)                      ! conductive heat flux from the surface
+         !---computation of surface temperature increment  
+         zdts    = -zfts / zdfts
+         !---computation of the new surface temperature 
+         sist_1d(ji) = sist_1d(ji) + zdts
+      END DO
+
+      !----------------------------------------------------------------------------
+      !  5. Boundary condition at the top surface
+      !--    IF Tsb < Tmelt, Fnet = Fcs (the net heat flux equal the conductive heat flux)
+      !      Otherwise Tsb = Tmelt and Qnet(Tmelt) > 0 
+      !      Fnet(Tmelt) is therefore the net surface flux needed for melting
+      !----------------------------------------------------------------------------
        
        
-       !     5.1.  Limitation of surface temperature and update total non solar fluxes,
-       !          latent heat flux and conductive flux at the top surface 
-       !----------------------------------------------------------------------  
+      !     5.1.  Limitation of surface temperature and update total non solar fluxes,
+      !          latent heat flux and conductive flux at the top surface 
+      !----------------------------------------------------------------------  
                      
-       IF ( .NOT. lk_cpl ) THEN   ! duplicate the loop for performances issues
-          DO ji = kideb, kiut
-             sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) )
-             qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
-             qla_ice_1d(ji) = qla_ice_1d(ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
-             zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
-          END DO
-       ELSE
-          DO ji = kideb, kiut
-             sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) )
-             zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
-          END DO
-       ENDIF
-
-       !     5.2. Calculate available heat for surface ablation. 
-       !---------------------------------------------------------------------
-
-       DO ji = kideb, kiut
-          zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)          
-          zfnet(ji) = MAX( zzero , zfnet(ji) )
-          zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) )
-       END DO
-
-       !-------------------------------------------------------------------------
-       !  6. Calculate changes in internal temperature due to vertical diffusion   
-       !     processes. The evolution of this temperature is governed by the one-
-       !     dimensionnal heat-diffusion equation. 
-       !     Given the temperature tbif(1/2/3), at time m we solve a set
-       !     of finite difference equations to obtain new tempe. Each tempe is coupled
-       !     to the temp. immediatly above and below by heat conduction terms. Thus 
-       !     we have a set of equations of the form A * T = B, where A is a tridiagonal
-       !     matrix, T a vector whose components are the unknown new temp.
-       !-------------------------------------------------------------------------
+      IF ( .NOT. lk_cpl ) THEN   ! duplicate the loop for performances issues
+         DO ji = kideb, kiut
+            sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) )
+            qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
+            qla_ice_1d(ji) = qla_ice_1d(ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
+            zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
+         END DO
+      ELSE
+         DO ji = kideb, kiut
+            sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) )
+            zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
+         END DO
+      ENDIF
+
+      !     5.2. Calculate available heat for surface ablation. 
+      !---------------------------------------------------------------------
+
+      DO ji = kideb, kiut
+         zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)          
+         zfnet(ji) = MAX( zzero , zfnet(ji) )
+         zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) )
+      END DO
+
+      !-------------------------------------------------------------------------
+      !  6. Calculate changes in internal temperature due to vertical diffusion   
+      !     processes. The evolution of this temperature is governed by the one-
+      !     dimensionnal heat-diffusion equation. 
+      !     Given the temperature tbif(1/2/3), at time m we solve a set
+      !     of finite difference equations to obtain new tempe. Each tempe is coupled
+      !     to the temp. immediatly above and below by heat conduction terms. Thus 
+      !     we have a set of equations of the form A * T = B, where A is a tridiagonal
+      !     matrix, T a vector whose components are the unknown new temp.
+      !-------------------------------------------------------------------------
        
-       !--parameter for the numerical methode use to solve the heat-diffusion equation
-       !- implicit, explicit or Crank-Nicholson
-       zumsb = 1.0 - sbeta  
-       DO ji = kideb, kiut
-          zidsn(ji)   = MAX ( zzero, SIGN( zone, hsndif - h_snow_1d(ji) ) ) 
-          z1midsn(ji) = 1.0 - zidsn(ji)
-          zihic       = MAX ( zzero, SIGN( zone, hicdif - h_ice_1d(ji) ) ) 
-          zidsnic(ji) = zidsn(ji) *  zihic 
-          zfcsudt(ji) = zfcsu(ji) * rdt_ice 
-       END DO 
+      !--parameter for the numerical methode use to solve the heat-diffusion equation
+      !- implicit, explicit or Crank-Nicholson
+      zumsb = 1.0 - sbeta  
+      DO ji = kideb, kiut
+         zidsn(ji)   = MAX ( zzero, SIGN( zone, hsndif - h_snow_1d(ji) ) ) 
+         z1midsn(ji) = 1.0 - zidsn(ji)
+         zihic       = MAX ( zzero, SIGN( zone, hicdif - h_ice_1d(ji) ) ) 
+         zidsnic(ji) = zidsn(ji) *  zihic 
+         zfcsudt(ji) = zfcsu(ji) * rdt_ice 
+      END DO 
    
-       DO ji = kideb, kiut
-
-          !     6.1 Calculate intermediate variables.
-          !----------------------------------------
-
-          !--conductivity at the snow surface
-          zkhsn = 2.0 * zksn(ji) * rdt_ice / rcpsn
-          !--conductivity at the ice surface
-          zkhic = 4.0 * zkic(ji) * rdt_ice / MAX( h_ice_1d(ji) * h_ice_1d(ji) * rcpic , epsi20 )
-          !--conductivity at the snow/ice interface 
-          zkint = 4.0 * zksn(ji) * zkic(ji)  &
-             &        / ( zksn(ji) * h_ice_1d(ji) + 2.0 * zkic(ji) * h_snow_1d(ji) * z1midsn(ji)) 
-          zkhsnint = zkint * rdt_ice / rcpsn
-          zkhicint = zkint * 2.0 * rdt_ice / MAX( h_ice_1d(ji) * rcpic , epsi20 )
-          
-          !     6.2. Fulfill the linear system matrix.
-          !-----------------------------------------
+      DO ji = kideb, kiut
+
+         !     6.1 Calculate intermediate variables.
+         !----------------------------------------
+
+         !--conductivity at the snow surface
+         zkhsn = 2.0 * zksn(ji) * rdt_ice / rcpsn
+         !--conductivity at the ice surface
+         zkhic = 4.0 * zkic(ji) * rdt_ice / MAX( h_ice_1d(ji) * h_ice_1d(ji) * rcpic , epsi20 )
+         !--conductivity at the snow/ice interface 
+         zkint = 4.0 * zksn(ji) * zkic(ji)  &
+            &        / ( zksn(ji) * h_ice_1d(ji) + 2.0 * zkic(ji) * h_snow_1d(ji) * z1midsn(ji)) 
+         zkhsnint = zkint * rdt_ice / rcpsn
+         zkhicint = zkint * 2.0 * rdt_ice / MAX( h_ice_1d(ji) * rcpic , epsi20 )
+          
+         !     6.2. Fulfill the linear system matrix.
+         !-----------------------------------------
 !$$$          zplediag(1) = 1 + sbeta * z1midsn(ji) * ( zkhsn + zkhsnint )       
-          zplediag(1) =   zidsn(ji) + z1midsn(ji) * h_snow_1d(ji)   &
-             &          + sbeta * z1midsn(ji) * zkhsnint 
-          zplediag(2) = 1 + sbeta * ( z1midsn(ji) * zkhicint + zkhic ) 
-          zplediag(3) = 1 + 3.0 * sbeta * zkhic   
-
-          zsubdiag(1) =  0.e0              
-          zsubdiag(2) = -1.e0 * z1midsn(ji) * sbeta * zkhicint
-          zsubdiag(3) = -1.e0 * sbeta * zkhic 
-
-          zsupdiag(1) = -1.e0 * z1midsn(ji) * sbeta * zkhsnint 
-          zsupdiag(2) = zsubdiag(3)
-          zsupdiag(3) =  0.e0
-          
-          !     6.3. Fulfill the idependent term vector.
-          !-------------------------------------------
+         zplediag(1) =   zidsn(ji) + z1midsn(ji) * h_snow_1d(ji)   &
+            &          + sbeta * z1midsn(ji) * zkhsnint 
+         zplediag(2) = 1 + sbeta * ( z1midsn(ji) * zkhicint + zkhic ) 
+         zplediag(3) = 1 + 3.0 * sbeta * zkhic   
+
+         zsubdiag(1) =  0.e0              
+         zsubdiag(2) = -1.e0 * z1midsn(ji) * sbeta * zkhicint
+         zsubdiag(3) = -1.e0 * sbeta * zkhic 
+
+         zsupdiag(1) = -1.e0 * z1midsn(ji) * sbeta * zkhsnint 
+         zsupdiag(2) = zsubdiag(3)
+         zsupdiag(3) =  0.e0
+          
+         !     6.3. Fulfill the idependent term vector.
+         !-------------------------------------------
           
 !$$$          zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *   &
@@ -402,58 +395,58 @@
 !$$$             &         - zumsb * ( zkhsn * tbif_1d(ji,1)
 !$$$             &                   + zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) ) )
-          zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *    &
-             &       ( h_snow_1d(ji) * tbif_1d(ji,1) - ( zfcsudt(ji) / rcpsn )  &
-             &       - zumsb * zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) )
-
-          zsmbr(2) =  tbif_1d(ji,2)  &
-             &      - zidsn(ji) * ( 1.0 - zidsnic(ji) ) &
-             &        * ( zfcsudt(ji) / MAX( h_ice_1d(ji) * rcpic , epsi20 ) ) &
-             &      + zumsb * ( zkhicint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) &
-             &                   - zkhic * ( tbif_1d(ji,2) - tbif_1d(ji,3) )  )
-
-          zsmbr(3) =  tbif_1d(ji,3)  &
-             &      + zkhic * ( 2.0 * tfu_1d(ji) &
-             &                + zumsb * ( tbif_1d(ji,2) - 3.0 * tbif_1d(ji,3) ) )
-          
-          !     6.4. Solve the system (Gauss elimination method).
-          !----------------------------------------------------
-          
-          zpiv1 = zsubdiag(2) / zplediag(1) 
-          zb2   = zplediag(2) - zpiv1 * zsupdiag(1)
-          zd2   = zsmbr(2) - zpiv1 * zsmbr(1)
-
-          zpiv2 = zsubdiag(3) / zb2
-          zb3   = zplediag(3) - zpiv2 * zsupdiag(2)
-          zd3   = zsmbr(3) - zpiv2 * zd2
-
-          tbif_1d(ji,3) = zd3 / zb3
-          tbif_1d(ji,2) = ( zd2 - zsupdiag(2) * tbif_1d(ji,3) ) / zb2
-          tbif_1d(ji,1) = ( zsmbr(1) - zsupdiag(1) * tbif_1d(ji,2) ) / zplediag(1)            
-
-          !--- taking into account the particular case of  zidsnic(ji) = 1
-          ztint =  (  zkic(ji) * h_snow_1d(ji) * tfu_1d (ji)    &
-             &      + zksn(ji) * h_ice_1d(ji) * sist_1d(ji) )   &
-             &   / ( zkic(ji) * h_snow_1d(ji) + zksn(ji) * h_ice_1d(ji) ) 
-
-          tbif_1d(ji,1) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,1)   &
-             &                + zidsnic(ji)   * ( ztint + sist_1d(ji) ) / 2.0
-          tbif_1d(ji,2) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,2)   &
-             &                + zidsnic(ji)   * ( 3.0 * ztint + tfu_1d(ji) ) / 4.0
-          tbif_1d(ji,3) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,3)   &
-             &                + zidsnic(ji)   * ( ztint + 3.0 * tfu_1d(ji) ) / 4.0     
-       END DO
+         zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *    &
+            &       ( h_snow_1d(ji) * tbif_1d(ji,1) - ( zfcsudt(ji) / rcpsn )  &
+            &       - zumsb * zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) )
+
+         zsmbr(2) =  tbif_1d(ji,2)  &
+            &      - zidsn(ji) * ( 1.0 - zidsnic(ji) ) &
+            &        * ( zfcsudt(ji) / MAX( h_ice_1d(ji) * rcpic , epsi20 ) ) &
+            &      + zumsb * ( zkhicint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) &
+            &                   - zkhic * ( tbif_1d(ji,2) - tbif_1d(ji,3) )  )
+
+         zsmbr(3) =  tbif_1d(ji,3)  &
+            &      + zkhic * ( 2.0 * tfu_1d(ji) &
+            &                + zumsb * ( tbif_1d(ji,2) - 3.0 * tbif_1d(ji,3) ) )
+          
+         !     6.4. Solve the system (Gauss elimination method).
+         !----------------------------------------------------
+          
+         zpiv1 = zsubdiag(2) / zplediag(1) 
+         zb2   = zplediag(2) - zpiv1 * zsupdiag(1)
+         zd2   = zsmbr(2) - zpiv1 * zsmbr(1)
+
+         zpiv2 = zsubdiag(3) / zb2
+         zb3   = zplediag(3) - zpiv2 * zsupdiag(2)
+         zd3   = zsmbr(3) - zpiv2 * zd2
+
+         tbif_1d(ji,3) = zd3 / zb3
+         tbif_1d(ji,2) = ( zd2 - zsupdiag(2) * tbif_1d(ji,3) ) / zb2
+         tbif_1d(ji,1) = ( zsmbr(1) - zsupdiag(1) * tbif_1d(ji,2) ) / zplediag(1)            
+
+         !--- taking into account the particular case of  zidsnic(ji) = 1
+         ztint =  (  zkic(ji) * h_snow_1d(ji) * tfu_1d (ji)    &
+            &      + zksn(ji) * h_ice_1d(ji) * sist_1d(ji) )   &
+            &   / ( zkic(ji) * h_snow_1d(ji) + zksn(ji) * h_ice_1d(ji) ) 
+
+         tbif_1d(ji,1) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,1)   &
+            &                + zidsnic(ji)   * ( ztint + sist_1d(ji) ) / 2.0
+         tbif_1d(ji,2) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,2)   &
+            &                + zidsnic(ji)   * ( 3.0 * ztint + tfu_1d(ji) ) / 4.0
+         tbif_1d(ji,3) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,3)   &
+            &                + zidsnic(ji)   * ( ztint + 3.0 * tfu_1d(ji) ) / 4.0     
+      END DO
  
-       !----------------------------------------------------------------------
-       !  9. Take into account surface ablation and bottom accretion-ablation.|
-       !----------------------------------------------------------------------
+      !----------------------------------------------------------------------
+      !  9. Take into account surface ablation and bottom accretion-ablation.|
+      !----------------------------------------------------------------------
        
-       !---Snow accumulation in one thermodynamic time step
-       zsprecip(kideb:kiut) = sprecip_1d(kideb:kiut) * rdt_ice / rhosn
-
-
-       DO ji = kideb, kiut
-          
-          !      9.1. Surface ablation and update of snow thickness and qstbif_1d
-          !--------------------------------------------------------------------
+      !---Snow accumulation in one thermodynamic time step
+      zsprecip(kideb:kiut) = sprecip_1d(kideb:kiut) * rdt_ice / rhosn
+
+
+      DO ji = kideb, kiut
+          
+         !      9.1. Surface ablation and update of snow thickness and qstbif_1d
+         !--------------------------------------------------------------------
           
           !--------------------------------------------------------------------------
@@ -759,8 +752,8 @@
           qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji)
           frld_1d(ji)    = zfrld_1d(ji)
-          !
-       END DO
-       ! 
-    END SUBROUTINE lim_thd_zdf_2
+         !
+      END DO
+      ! 
+   END SUBROUTINE lim_thd_zdf_2
 
 #else
@@ -768,7 +761,4 @@
    !!   Default Option                                     NO sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_thd_zdf_2          ! Empty routine
-   END SUBROUTINE lim_thd_zdf_2
 #endif
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limtrp_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 1855)
@@ -4,4 +4,8 @@
    !! LIM 2.0 transport ice model : sea-ice advection/diffusion
    !!======================================================================
+   !! History :  LIM  ! 2000-01 (LIM)  Original code
+   !!            1.0  ! 2001-05 (G. Madec, R. Hordoir) opa norm
+   !!            2.0  ! 2004-01 (G. Madec, C. Ethe)  F90, mpp
+   !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
@@ -11,5 +15,4 @@
    !!   lim_trp_init_2 : initialization and namelist read
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE phycst
    USE dom_oce
@@ -26,25 +29,20 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_trp_2     ! called by sbc_ice_lim_2
-
-   !! * Shared module variables
-   REAL(wp), PUBLIC  ::   &  !:
-      bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip)
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      epsi06 = 1.e-06  ,  &
-      epsi03 = 1.e-03  ,  &
-      epsi16 = 1.e-16  ,  &
-      rzero  = 0.e0    ,  &
-      rone   = 1.e0
+   PUBLIC   lim_trp_2   ! called by sbc_ice_lim_2
+
+   REAL(wp), PUBLIC  ::   bound  = 0.e0   !: boundary condit. (0.0 no-slip, 1.0 free-slip)
+
+   REAL(wp) ::   epsi06 = 1.e-06   ! constant values
+   REAL(wp) ::   epsi03 = 1.e-03   !
+   REAL(wp) ::   epsi16 = 1.e-16   !
+   REAL(wp) ::   rzero  = 0.e0     !
+   REAL(wp) ::   rone   = 1.e0     !
 
    !! * Substitution
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -62,34 +60,13 @@
       !!
       !! ** action :
-      !!
-      !! History :
-      !!   1.0  !  00-01 (LIM)  Original code
-      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
-      !!   2.0  !  04-01 (G. Madec, C. Ethe)  F90, mpp
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
-
-      INTEGER  ::   ji, jj, jk,   &  ! dummy loop indices
-         &          initad           ! number of sub-timestep for the advection
-
-      REAL(wp) ::  &                              
-         zindb  ,  &
-         zacrith, &
-         zindsn , &
-         zindic , &
-         zusvosn, &
-         zusvoic, &
-         zignm  , &
-         zindhe , &
-         zvbord , &
-         zcfl   , &
-         zusnit , &
-         zrtt, ztsn, ztic1, ztic2
-
-      REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace
-         zui_u , zvi_v , zsm   ,         &
-         zs0ice, zs0sn , zs0a  ,         &
-         zs0c0 , zs0c1 , zs0c2 ,         &
-         zs0st
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   initad       ! number of sub-timestep for the advection
+      REAL(wp) ::   zindb , zacrith, zindsn , zignm , zvbord, zrtt, ztic1, zusnit
+      REAL(wp) ::   zindic, zusvosn, zusvoic, zindhe, zcfl  , ztsn, ztic2
+      REAL(wp), DIMENSION(jpi,jpj) ::   zui_u , zs0sn, zs0c0, zs0a , zsm      ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zvi_v , zs0st, zs0c1, zs0c2, zs0ice   !  -      -
       !---------------------------------------------------------------------
 
@@ -98,8 +75,7 @@
       zsm(:,:) = area(:,:)
       
-      IF( ln_limdyn ) THEN
-         !-------------------------------------!
-         !   Advection of sea ice properties   !
-         !-------------------------------------!
+      !                             !-------------------------------------!
+      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
+         !                          !-------------------------------------!
 
          ! ice velocities at ocean U- and V-points (zui_u,zvi_v)
@@ -113,7 +89,5 @@
             END DO
          END DO
-         ! Lateral boundary conditions on zui_u, zvi_v
-         CALL lbc_lnk( zui_u, 'U', -1. )
-         CALL lbc_lnk( zvi_v, 'V', -1. )
+         CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )         ! Lateral boundary conditions
 
          ! CFL test for stability
@@ -123,7 +97,7 @@
          zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
 
-         IF (lk_mpp ) CALL mpp_max(zcfl)
-
-         IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
+         IF( lk_mpp ) CALL mpp_max( zcfl )
+
+         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : CFL violation at the ',nday,'th day, cfl = ',zcfl
 
          ! content of properties
@@ -144,5 +118,5 @@
          zusnit = 1.0 / REAL( initad ) 
          
-         IF ( MOD( nday , 2 ) == 0) THEN
+         IF( MOD( nday , 2 ) == 0) THEN
             DO jk = 1,initad
                CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
@@ -228,12 +202,8 @@
          !   Up-dating and limitation of sea ice properties after transport   !
          ! -------------------------------------------------------------------!
-
-         ! Up-dating and limitation of sea ice properties after transport.
          DO jj = 1, jpj
-!!!iii      zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) )              !ibug mpp  !!bugmpp  njeq!
             zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
             DO ji = 1, jpi
-
-               ! Recover mean values over the grid squares.
+               !                                                        ! Recover mean values over the grid squares.
                zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) )
                zs0ice(ji,jj) = MAX( rzero, zs0ice(ji,jj)/area(ji,jj) )
@@ -243,6 +213,5 @@
                zs0c2 (ji,jj) = MAX( rzero, zs0c2 (ji,jj)/area(ji,jj) )
                zs0st (ji,jj) = MAX( rzero, zs0st (ji,jj)/area(ji,jj) )
-
-               ! Recover in situ values.
+               !                                                        ! Recover in situ values.
                zindb         = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) )
                zacrith       = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )
@@ -262,8 +231,8 @@
                zrtt          = 173.15 * rone 
                ztsn          =          zignm   * tbif(ji,jj,1)  &
-                              + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) ) 
+                  &           + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) ) 
                ztic1          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) )
                ztic2          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) )
- 
+               !
                tbif(ji,jj,1) = zindsn * ztsn  + ( 1.0 - zindsn ) * tfu(ji,jj)               
                tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)
@@ -272,7 +241,7 @@
             END DO
          END DO
-         
+         !
       ENDIF
-      
+      !
    END SUBROUTINE lim_trp_2
 
@@ -288,15 +257,11 @@
       !!
       !! ** input   :   Namelist namicetrp
-      !!
-      !! history :
-      !!   2.0  !  03-08 (C. Ethe)  Original code
       !!-------------------------------------------------------------------
       NAMELIST/namicetrp/ bound
       !!-------------------------------------------------------------------
-
-      ! Read Namelist namicetrp
-      REWIND ( numnam_ice )
+      !
+      REWIND ( numnam_ice )                  ! Read Namelist namicetrp
       READ   ( numnam_ice  , namicetrp )
-      IF(lwp) THEN
+      IF(lwp) THEN                           ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'lim_trp_init_2 : Ice parameters for advection '
@@ -304,5 +269,5 @@
          WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound
       ENDIF
-            
+      ! 
    END SUBROUTINE lim_trp_init_2
 
@@ -311,7 +276,4 @@
    !!   Default option         Empty Module                No sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_trp_2        ! Empty routine
-   END SUBROUTINE lim_trp_2
 #endif
 
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limwri_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limwri_2.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limwri_2.F90	(revision 1855)
@@ -12,10 +12,7 @@
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   lim_wri_2      : write of the diagnostics variables in ouput file 
-   !!   lim_wri_init_2 : initialization and namelist read
-   !!   lim_wri_state_2 : write for initial state or/and abandon:
-   !!                     > output.init.nc (if ninist = 1 in namelist)
-   !!                     > output.abort.nc
+   !!   lim_wri_2       : write of the diagnostics variables in ouput file 
+   !!   lim_wri_init_2  : initialization and namelist read
+   !!   lim_wri_state_2 : write for initial state (output.init.nc if ninist=1) or/and in the abort file
    !!----------------------------------------------------------------------
    USE phycst
@@ -42,6 +39,6 @@
    INTEGER, PARAMETER                       ::   jpnoumax = 40   ! maximum number of variable for ice output
    INTEGER                                  ::   noumef          ! number of fields
-   REAL(wp)           , DIMENSION(jpnoumax) ::   cmulti ,     &  ! multiplicative constant
-      &                                          cadd            ! additive constant
+   REAL(wp)           , DIMENSION(jpnoumax) ::   cmulti         ! multiplicative constant
+   REAL(wp)           , DIMENSION(jpnoumax) ::   cadd            ! additive constant
    CHARACTER(len = 35), DIMENSION(jpnoumax) ::   titn            ! title of the field
    CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   nam             ! name of the field
@@ -52,13 +49,12 @@
    INTEGER , DIMENSION( jpij ) ::   ndex51              ! ????
 
-   REAL(wp)  ::            &  ! constant values
-      epsi16 = 1.e-16   ,  &
-      zzero  = 0.e0     ,  &
-      zone   = 1.e0
+   REAL(wp) ::   epsi16 = 1.e-16   ! constant values
+   REAL(wp) ::   rzero  = 0.e0     !
+   REAL(wp) ::   rone   = 1.e0     !
 
    !! * Substitutions
 #   include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006)
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -81,14 +77,13 @@
       !!
       !! ** Method  :   computes the average of some variables and write
-      !!      it in the NetCDF ouput files
-      !!      CAUTION: the sea-ice time-step must be an integer fraction
-      !!      of a day
+      !!              it in the NetCDF ouput files
+      !! CAUTION: the sea-ice time-step must be an integer fraction of a day
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
       !!
-      INTEGER  ::   ji, jj, jf                      ! dummy loop indices
+      INTEGER  ::   ji, jj, jf   ! dummy loop indices
       CHARACTER(len = 40)  ::   clhstnam, clop
-      REAL(wp) ::   zsto, zjulian, zout,   &  ! temporary scalars
-         &          zindh, zinda, zindb, ztmu
+      REAL(wp) ::   zsto, zjulian, zout          ! temporary scalars
+      REAL(wp) ::   zindh, zinda, zindb, ztmu
       REAL(wp), DIMENSION(1)                ::   zdept
       REAL(wp), DIMENSION(jpi,jpj)          ::   zfield
@@ -132,8 +127,8 @@
       DO jj = 2 , jpjm1
          DO ji = 1 , jpim1   ! NO vector opt.
-            zindh  = MAX( zzero , SIGN( zone , hicif(ji,jj) * (1.0 - frld(ji,jj) ) - 0.10 ) )
-            zinda  = MAX( zzero , SIGN( zone , ( 1.0 - frld(ji,jj) ) - 0.10 ) )
+            zindh  = MAX( rzero , SIGN( rone , hicif(ji,jj) * (1.0 - frld(ji,jj) ) - 0.10 ) )
+            zinda  = MAX( rzero , SIGN( rone , ( 1.0 - frld(ji,jj) ) - 0.10 ) )
             zindb  = zindh * zinda
-            ztmu   = MAX( 0.5 * zone , ( tmu(ji,jj) + tmu(ji+1,jj) + tmu(ji,jj+1) + tmu(ji+1,jj+1) ) ) 
+            ztmu   = MAX( 0.5 * rone , ( tmu(ji,jj) + tmu(ji+1,jj) + tmu(ji,jj+1) + tmu(ji+1,jj+1) ) ) 
             zcmo(ji,jj,1)  = hsnif (ji,jj)
             zcmo(ji,jj,2)  = hicif (ji,jj)
@@ -143,10 +138,10 @@
             zcmo(ji,jj,6)  = fbif  (ji,jj)
             zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
-                                      + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
-                                  / ztmu 
+               &                      + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
+               &                  / ztmu 
 
             zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
-                                      + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
-                                  / ztmu
+               &                      + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
+               &                  / ztmu
             zcmo(ji,jj,9)  = sst_m(ji,jj)
             zcmo(ji,jj,10) = sss_m(ji,jj)
@@ -200,5 +195,5 @@
       !! ** input   :   Namelist namicewri
       !!-------------------------------------------------------------------
-      INTEGER ::   nf      ! ???
+      INTEGER ::   jf      ! dummy loop indices
       TYPE FIELD 
          CHARACTER(len = 35) :: ztitle 
@@ -209,16 +204,13 @@
          REAL                :: zcadd        
       END TYPE FIELD
-      TYPE(FIELD) ::  &
-         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   &
-         field_7 , field_8 , field_9 , field_10, field_11, field_12,   &
-         field_13, field_14, field_15, field_16, field_17, field_18,   &
-         field_19
+      TYPE(FIELD) ::   field_1 , field_2 , field_3 , field_4 , field_5 , field_6  
+      TYPE(FIELD) ::   field_7 , field_8 , field_9 , field_10, field_11, field_12
+      TYPE(FIELD) ::   field_13, field_14, field_15, field_16, field_17, field_18, field_19
       TYPE(FIELD) , DIMENSION(jpnoumax) :: zfield
-
-      NAMELIST/namiceout/ noumef, &
-         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   &
-         field_7 , field_8 , field_9 , field_10, field_11, field_12,   &
-         field_13, field_14, field_15, field_16, field_17, field_18,   &
-         field_19
+      !!
+      NAMELIST/namiceout/ noumef,                                          &
+         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,       &
+         field_7 , field_8 , field_9 , field_10, field_11, field_12,       &
+         field_13, field_14, field_15, field_16, field_17, field_18, field_19
       !!-------------------------------------------------------------------
 
@@ -246,11 +238,11 @@
       zfield(19) = field_19
       
-      DO nf = 1, noumef
-         titn  (nf) = zfield(nf)%ztitle
-         nam   (nf) = zfield(nf)%zname
-         uni   (nf) = zfield(nf)%zunit
-         nc    (nf) = zfield(nf)%znc
-         cmulti(nf) = zfield(nf)%zcmulti
-         cadd  (nf) = zfield(nf)%zcadd
+      DO jf = 1, noumef
+         titn  (jf) = zfield(jf)%ztitle
+         nam   (jf) = zfield(jf)%zname
+         uni   (jf) = zfield(jf)%zunit
+         nc    (jf) = zfield(jf)%znc
+         cmulti(jf) = zfield(jf)%zcmulti
+         cadd  (jf) = zfield(jf)%zcadd
       END DO
 
@@ -262,7 +254,7 @@
          WRITE(numout,*) '           title                            name     unit      Saving (1/0) ',   &
             &            '    multiplicative constant       additive constant '
-         DO nf = 1 , noumef         
-            WRITE(numout,*) '   ', titn(nf), '   ', nam(nf),'      ', uni(nf),'  ', nc(nf),'        ', cmulti(nf),   &
-               &       '        ', cadd(nf)
+         DO jf = 1 , noumef         
+            WRITE(numout,*) '   ', titn(jf), '   ', nam(jf),'      ', uni(jf),'  ', nc(jf),'        ', cmulti(jf),   &
+               &       '        ', cadd(jf)
          END DO
       ENDIF
@@ -281,7 +273,4 @@
       !!        Used to find errors in the initial state or save the last
       !!      ocean state in case of abnormal end of a simulation
-      !!
-      !! History :
-      !!   2.0  !  2009-06  (B. Lemaire)
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt               ! ocean time-step index)
Index: anches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 1854)
+++ 	(revision )
@@ -1,87 +1,0 @@
-MODULE thd_ice_2
-#if defined key_lim2
-   !!======================================================================
-   !!                       ***  MODULE thd_ice_2  ***
-   !! LIM 2.0 sea-ice :   Ice thermodynamics in 1D
-   !!=====================================================================
-   !! History :
-   !!   2.0  !  02-11  (C. Ethe)  F90: Free form and module
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE par_ice_2
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Share Module variables
-   REAL(wp) , PUBLIC ::   & !!! ** ice-thermo namelist (namicethd) **
-      hmelt   = -0.15  ,  &  !: maximum melting at the bottom
-      hicmin  = 0.2    ,  &  !: ice th. corr. to max. ener. in brine pocket
-      hiclim  = 0.05   ,  &  !: minimum ice thickness
-      amax    = 0.999  ,  &  !: maximum lead fraction
-      swiqst  = 1.0    ,  &  !: energy stored in brine pocket (1) or not (0)
-      sbeta   = 1.0    ,  &  !: numerical scheme for diffusion in ice 
-      parlat  = 0.0    ,  &  !: percent. of energy used for lateral ablation
-      hakspl  = 0.5    ,  &  !: slope of distr. for Hakkinen-Mellro's lat. melt
-      hibspl  = 0.5    ,  &  !: slope of distribution for Hibler's lat. melt
-      exld    = 2.0    ,  &  !: exponent for leads-closure rate
-      hakdif  = 1.0    ,  &  !: coefficient for diffusions of ice and snow
-      thth    = 0.2    ,  &  !: thick. for comp. of eq. thermal conduct
-      hnzst   = 0.1    ,  &  !: thick. of the surf. layer in temp. comp.
-      parsub  = 1.0    ,  &  !: switch for snow sublimation or not
-      alphs   = 1.0          !: coef. for snow density when snow-ice formation
-
-   REAL(wp), PUBLIC, DIMENSION(2)  ::  &  !:   
-      hiccrit = (/0.3,0.3/)  !: ice th. for lateral accretion in the NH (SH) (m)
-
-   REAL(wp) , PUBLIC ::   &  !:
-      uscomi,             &  !: inverse of minimum lead fraction
-      cnscg                  !: ratio  rcpsn/rcpic
-
-   INTEGER , PUBLIC, DIMENSION(jpij) ::   &  !:
-      npb     ,   &   !: number of points where computations has to be done
-      npac            !: correspondance between the points
-
-   REAL(wp), PUBLIC, DIMENSION(jpij) ::   &  !: 
-      qldif_1d    ,     &  !: corresponding to the 2D var  qldif
-      qcmif_1d    ,     &  !: corresponding to the 2D var  qcmif
-      thcm_1d     ,     &  !:    "                  "      thcm
-      fstbif_1d   ,     &  !:    "                  "      fstric
-      fltbif_1d   ,     &  !:    "                  "      ffltbif
-      fscbq_1d    ,     &  !:    "                  "      fscmcbq
-      qsr_ice_1d  ,     &  !:    "                  "      qsr_ice
-      fr1_i0_1d   ,     &  !:    "                  "      fr1_i0
-      fr2_i0_1d   ,     &  !:    "                  "      fr2_i0
-      qns_ice_1d  ,     &  !:    "                  "      qns_ice
-      qfvbq_1d    ,     &  !:    "                  "      qfvbq
-      sist_1d     ,     &  !:    "                  "      sist
-      tfu_1d      ,     &  !:    "                  "      tfu
-      sprecip_1d  ,     &  !:    "                  "      sprecip
-      h_snow_1d   ,     &  !:    "                  "      h_snow
-      h_ice_1d    ,     &  !:    "                  "      h_ice
-      frld_1d     ,     &  !:    "                  "      frld
-      qstbif_1d   ,     &  !:    "                  "      qstoif
-      fbif_1d     ,     &  !:    "                  "      fbif
-      rdmicif_1d  ,     &  !:    "                  "      rdmicif
-      rdmsnif_1d  ,     &  !:    "                  "      rdmsnif
-      qlbbq_1d    ,     &  !:    "                  "      qlbsbq
-      dmgwi_1d    ,     &  !:    "                  "      dmgwi
-      dvsbq_1d    ,     &  !:    "                  "      rdvosif
-      rdvomif_1d  ,     &  !:    "                  "      rdvomif
-      dvbbq_1d    ,     &  !:    "                  "      rdvobif
-      dvlbq_1d    ,     &  !:    "                  "      rdvolif
-      dvnbq_1d    ,     &  !:    "                  "      rdvolif
-      dqns_ice_1d ,     &  !:    "                  "      dqns_ice
-      qla_ice_1d  ,     &  !:    "                  "      qla_ice
-      dqla_ice_1d          !:    "                  "      dqla_ice
-
-   REAL(wp), PUBLIC, DIMENSION(jpij,jplayersp1) ::   &  !:
-      tbif_1d              !: corresponding to the 2D var  tbif
-
-   !!======================================================================
-#endif
-END MODULE thd_ice_2
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/limtab.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/limtab.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/limtab.F90	(revision 1855)
@@ -2,14 +2,13 @@
    !!======================================================================
    !!                       ***  MODULE limtab   ***
-   !!              transform 1D (2D) array to a 2D (1D) table
+   !!   LIM-3 ice model : transform 1D (2D) array to a 2D (1D) array
    !!======================================================================
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!   'key_lim3' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
    !!   tab_2d_1d  : 2-D to 1-D
    !!   tab_1d_2d  : 1-D to 2-D
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_kind
 
@@ -17,10 +16,9 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC tab_2d_1d  ! called by lim_ther
-   PUBLIC tab_1d_2d  ! called by lim_ther
+   PUBLIC   tab_2d_1d   ! called by lim_thd
+   PUBLIC   tab_1d_2d   ! called by lim_thd
 
    !!----------------------------------------------------------------------
-   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008) 
+   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -28,55 +26,49 @@
 CONTAINS
 
-   SUBROUTINE tab_2d_1d ( ndim1d, tab1d, tab2d, ndim2d_x, ndim2d_y, tab_ind )
-
-      INTEGER, INTENT(in) :: &
-         ndim1d, ndim2d_x, ndim2d_y
-
-      REAL(wp), DIMENSION (ndim2d_x, ndim2d_y), INTENT(in) ::  &
-         tab2d
-
-      INTEGER, DIMENSION ( ndim1d), INTENT ( in) :: &
-         tab_ind
-
-      REAL(wp), DIMENSION(ndim1d), INTENT ( out) ::  & 
-         tab1d
-
-      INTEGER ::  &
-         jn , jid, jjd
-
-      DO jn = 1, ndim1d
-         jid        = MOD( tab_ind(jn) - 1, ndim2d_x ) + 1
-         jjd        = ( tab_ind(jn) - 1 ) / ndim2d_x + 1
-         tab1d( jn) = tab2d( jid, jjd)
-      END DO
-
+   SUBROUTINE tab_2d_1d( kdim1d, ptab1d, ptab2d, kdim2d_i, kdim2d_j, ptab_ind )
+      !!-------------------------------------------------------------------
+      INTEGER , INTENT(in   )                                ::   kdim1d
+      INTEGER , INTENT(in   )                                ::   kdim2d_i, kdim2d_j
+      REAL(wp), INTENT(in   ), DIMENSION(kdim2d_i, kdim2d_j) ::   ptab2d
+      INTEGER , INTENT(in   ), DIMENSION(kdim1d)             ::   ptab_ind
+      REAL(wp), INTENT(  out), DIMENSION(kdim1d)             ::   ptab1d
+      !!
+      INTEGER ::   jn , jid, jjd   ! dummy loop indices
+      !!-------------------------------------------------------------------
+      ! 
+      DO jn = 1, kdim1d
+         jid = MOD( ptab_ind(jn) - 1, kdim2d_i ) + 1
+         jjd = ( ptab_ind(jn) - 1 ) / kdim2d_i + 1
+         ptab1d(jn) = ptab2d(jid,jjd)
+      END DO 
+      !
    END SUBROUTINE tab_2d_1d
 
 
-   SUBROUTINE tab_1d_2d ( ndim1d, tab2d, tab_ind, tab1d, ndim2d_x, ndim2d_y )
-
-      INTEGER, INTENT ( in) :: &
-         ndim1d, ndim2d_x, ndim2d_y
-
-      INTEGER, DIMENSION (ndim1d) , INTENT (in) :: &
-         tab_ind
-
-      REAL(wp), DIMENSION(ndim1d), INTENT (in) ::  &
-         tab1d  
-
-      REAL(wp), DIMENSION (ndim2d_x, ndim2d_y), INTENT ( out) :: &
-         tab2d
-
-      INTEGER :: &
-         jn, jid, jjd
-
-      DO jn = 1, ndim1d
-         jid             = MOD( tab_ind(jn) - 1, ndim2d_x) + 1
-         jjd             =    ( tab_ind(jn) - 1 ) / ndim2d_x  + 1
-         tab2d(jid, jjd) = tab1d( jn)
+   SUBROUTINE tab_1d_2d( kdim1d, ptab2d, ptab_ind, ptab1d, kdim2d_i, kdim2d_j )
+      !!-------------------------------------------------------------------
+      INTEGER , INTENT(in   )                                ::   kdim1d
+      INTEGER , INTENT(in   )                                ::   kdim2d_i, kdim2d_j
+      INTEGER , INTENT(in   ), DIMENSION(kdim1d)             ::   ptab_ind
+      REAL(wp), INTENT(in   ), DIMENSION(kdim1d)             ::   ptab1d
+      REAL(wp), INTENT(  out), DIMENSION(kdim2d_i, kdim2d_j) ::   ptab2d
+      !!
+      INTEGER ::   jn, jid, jjd   ! dummy loop indices
+      !!-------------------------------------------------------------------
+      !
+      DO jn = 1, kdim1d
+         jid = MOD( ptab_ind(jn) - 1, kdim2d_i) + 1
+         jjd =    ( ptab_ind(jn) - 1 ) / kdim2d_i  + 1
+         ptab2d(jid,jjd) = ptab1d(jn)
       END DO
-
+      !
    END SUBROUTINE tab_1d_2d
 
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
 #endif
+
+   !!======================================================================
 END MODULE limtab
Index: /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 1854)
+++ /branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 1855)
@@ -41,9 +41,9 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   alb_ice   !: albedo of ice
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau_ice    !: u-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau_ice    !: v-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr1_i0      !: 1st fraction of sol. rad. which penetrate inside the ice cover
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr2_i0      !: 2nd fraction of sol. rad. which penetrate inside the ice cover
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_ice     !: solid freshwater budget over ice: sublivation - snow
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau_ice    !: u-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau_ice    !: v-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2]
 
 # if defined key_lim3
