Index: branches/dev_003_CPL/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/dev_003_CPL/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 893)
+++ branches/dev_003_CPL/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 990)
@@ -24,5 +24,5 @@
    USE par_ice
    USE lib_mpp
-      
+
    IMPLICIT NONE
    PRIVATE
@@ -46,139 +46,134 @@
 
    SUBROUTINE lim_thd_dh(kideb,kiut,jl)
-       !!------------------------------------------------------------------
-       !!                ***  ROUTINE lim_thd_dh  ***
-       !!------------------------------------------------------------------
-       !! ** Purpose :
-       !!           This routine determines variations of ice and snow thicknesses.
-       !! ** Method  :
-       !!           Ice/Snow surface melting arises from imbalance in surface fluxes
-       !!           Bottom accretion/ablation arises from flux budget
-       !!           Snow thickness can increase by precipitation and decrease by 
-       !!              sublimation
-       !!           If snow load excesses Archmiede limit, snow-ice is formed by
-       !!              the flooding of sea-water in the snow
-       !! ** Steps  
-       !!           1) Compute available flux of heat for surface ablation
-       !!           2) Compute snow and sea ice enthalpies
-       !!           3) Surface ablation and sublimation
-       !!           4) Bottom accretion/ablation
-       !!           5) Case of Total ablation
-       !!           6) Snow ice formation
-       !!
-       !! ** Arguments
-       !!
-       !! ** Inputs / Outputs
-       !!
-       !! ** External
-       !!
-       !! ** References : Bitz and Lipscomb, JGR 99
-       !!                 Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
-       !!                 Vancoppenolle, Fichefet and Bitz, GRL 2005
-       !!                 Vancoppenolle et al., OM08
-       !!
-       !! ** History  : 
-       !!   original code    01-04 (LIM)
-       !!   original routine
-       !!               (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium
-       !!               (06/07-2005) 3D version
-       !!               (03-2008)    Clean code
-       !!
-       !!------------------------------------------------------------------
-       !! * Arguments
-       INTEGER , INTENT (IN) ::  &
-          kideb     ,         &  !: Start point on which the  the computation is applied
-          kiut      ,         &  !: End point on which the  the computation is applied
-          jl                     !: Thickness cateogry number
-
-       !! * Local variables
-       INTEGER ::             & 
-          ji        ,         &  !: space index 
-          jk        ,         &  !: ice layer index
-          isnow     ,         &  !: switch for presence (1) or absence (0) of snow
-          zji       ,         &  !: 2D corresponding indices to ji
-          zjj       ,         &
-          isnowic   ,         &  !: snow ice formation not
-          i_ice_switch   ,    &  !: ice thickness above a certain treshold or not
-          iter
-
-       REAL(wp) ::            &
-          zhsnew    ,         &  !: new snow thickness
-          zihgnew   ,         &  !: switch for total ablation
-          ztmelts   ,         &  !: melting point
-          zhn       ,         &
-          zdhcf     ,         &
-          zdhbf     ,         &
-          zhni      ,         &
-          zhnfi     ,         &
-          zihg      ,         &
-          zdhnm     ,         &
-          zhnnew    ,         &
-          zeps = 1.0e-13,     &
-          zhisn     ,         &
-          zfracs    ,         &  !: fractionation coefficient for bottom salt
-                                 !: entrapment
-          zds       ,         &  !: increment of bottom ice salinity
-          zcoeff    ,         &  !: dummy argument for snowfall partitioning
-                                 !: over ice and leads
-          zsm_snowice,        &  !: snow-ice salinity
-          zswi1     ,         &  !: switch for computation of bottom salinity
-          zswi12    ,         &  !: switch for computation of bottom salinity
-          zswi2     ,         &  !: switch for computation of bottom salinity
-          zgrr      ,         &  !: bottom growth rate
-          zihic     ,         &  !: 
-          ztform                 !: bottom formation temperature
-
-       REAL(wp) , DIMENSION(jpij) ::  &
-          zh_i      ,         &  ! ice layer thickness
-          zh_s      ,         &  ! snow layer thickness
-          ztfs      ,         &  ! melting point
-          zhsold    ,         &  ! old snow thickness
-          zqprec    ,         &  !: energy of fallen snow
-          zqfont_su ,         &  ! incoming, remaining surface energy
-          zqfont_bo              ! incoming, bottom energy
-
-       REAL(wp) , DIMENSION(jpij) :: &
-          z_f_surf,           &  ! surface heat for ablation
-          zhgnew                 ! new ice thickness
-
-       REAL(wp), DIMENSION(jpij) :: &
-          zdh_s_mel  ,        &  ! snow melt 
-          zdh_s_pre  ,        &  ! snow precipitation 
-          zdh_s_sub  ,        &  ! snow sublimation
-          zfsalt_melt            ! salt flux due to ice melt
-
-       REAL(wp) , DIMENSION(jpij,jkmax) :: &
-          zdeltah
-
-       ! Pathological cases
-       REAL(wp), DIMENSION(jpij) :: &
-          zfdt_init  ,        &  !: total incoming heat for ice melt
-          zfdt_final ,        &  !: total remaing heat for ice melt
-          zqt_i      ,        &  !: total ice heat content
-          zqt_s      ,        &  !: total snow heat content
-          zqt_dummy              !: dummy heat content
-
-       REAL(wp), DIMENSION(jpij,jkmax) :: &
-          zqt_i_lay              !: total ice heat content
-
-       ! Heat conservation 
-       REAL(wp), DIMENSION(jpij) :: &
-          zfbase,             &
-          zdq_i
-
-       INTEGER, DIMENSION(jpij) ::  &
-          innermelt
-
-       REAL(wp) :: &
-          meance_dh
-
-       INTEGER ::                   &
-          num_iter_max,       &
-          numce_dh
-
-!!-----------------------------------------------------------------------------
-
-      WRITE(numout,*) 'lim_thd_dh : computation of vertical snow/ice accretion/ablation'
-      WRITE(numout,*) '~~~~~~~~~'
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE lim_thd_dh  ***
+      !!------------------------------------------------------------------
+      !! ** Purpose :
+      !!           This routine determines variations of ice and snow thicknesses.
+      !! ** Method  :
+      !!           Ice/Snow surface melting arises from imbalance in surface fluxes
+      !!           Bottom accretion/ablation arises from flux budget
+      !!           Snow thickness can increase by precipitation and decrease by 
+      !!              sublimation
+      !!           If snow load excesses Archmiede limit, snow-ice is formed by
+      !!              the flooding of sea-water in the snow
+      !! ** Steps  
+      !!           1) Compute available flux of heat for surface ablation
+      !!           2) Compute snow and sea ice enthalpies
+      !!           3) Surface ablation and sublimation
+      !!           4) Bottom accretion/ablation
+      !!           5) Case of Total ablation
+      !!           6) Snow ice formation
+      !!
+      !! ** Arguments
+      !!
+      !! ** Inputs / Outputs
+      !!
+      !! ** External
+      !!
+      !! ** References : Bitz and Lipscomb, JGR 99
+      !!                 Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
+      !!                 Vancoppenolle, Fichefet and Bitz, GRL 2005
+      !!                 Vancoppenolle et al., OM08
+      !!
+      !! ** History  : 
+      !!   original code    01-04 (LIM)
+      !!   original routine
+      !!               (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium
+      !!               (06/07-2005) 3D version
+      !!               (03-2008)    Clean code
+      !!
+      !!------------------------------------------------------------------
+      !! * Arguments
+      INTEGER , INTENT (IN) ::  &
+         kideb     ,         &  !: Start point on which the  the computation is applied
+         kiut      ,         &  !: End point on which the  the computation is applied
+         jl                     !: Thickness cateogry number
+
+      !! * Local variables
+      INTEGER ::             & 
+         ji        ,         &  !: space index 
+         jk        ,         &  !: ice layer index
+         isnow     ,         &  !: switch for presence (1) or absence (0) of snow
+         zji       ,         &  !: 2D corresponding indices to ji
+         zjj       ,         &
+         isnowic   ,         &  !: snow ice formation not
+         i_ice_switch   ,    &  !: ice thickness above a certain treshold or not
+         iter
+
+      REAL(wp) ::            &
+         zhsnew    ,         &  !: new snow thickness
+         zihgnew   ,         &  !: switch for total ablation
+         ztmelts   ,         &  !: melting point
+         zhn       ,         &
+         zdhcf     ,         &
+         zdhbf     ,         &
+         zhni      ,         &
+         zhnfi     ,         &
+         zihg      ,         &
+         zdhnm     ,         &
+         zhnnew    ,         &
+         zeps = 1.0e-13,     &
+         zhisn     ,         &
+         zfracs    ,         &  !: fractionation coefficient for bottom salt
+                                !: entrapment
+         zds       ,         &  !: increment of bottom ice salinity
+         zcoeff    ,         &  !: dummy argument for snowfall partitioning
+                                !: over ice and leads
+         zsm_snowice,        &  !: snow-ice salinity
+         zswi1     ,         &  !: switch for computation of bottom salinity
+         zswi12    ,         &  !: switch for computation of bottom salinity
+         zswi2     ,         &  !: switch for computation of bottom salinity
+         zgrr      ,         &  !: bottom growth rate
+         zihic     ,         &  !: 
+         ztform                 !: bottom formation temperature
+
+      REAL(wp) , DIMENSION(jpij) ::  &
+         zh_i      ,         &  ! ice layer thickness
+         zh_s      ,         &  ! snow layer thickness
+         ztfs      ,         &  ! melting point
+         zhsold    ,         &  ! old snow thickness
+         zqprec    ,         &  !: energy of fallen snow
+         zqfont_su ,         &  ! incoming, remaining surface energy
+         zqfont_bo              ! incoming, bottom energy
+
+      REAL(wp) , DIMENSION(jpij) :: &
+         z_f_surf,           &  ! surface heat for ablation
+         zhgnew                 ! new ice thickness
+
+      REAL(wp), DIMENSION(jpij) :: &
+         zdh_s_mel  ,        &  ! snow melt 
+         zdh_s_pre  ,        &  ! snow precipitation 
+         zdh_s_sub  ,        &  ! snow sublimation
+         zfsalt_melt            ! salt flux due to ice melt
+
+      REAL(wp) , DIMENSION(jpij,jkmax) :: &
+         zdeltah
+
+      ! Pathological cases
+      REAL(wp), DIMENSION(jpij) :: &
+         zfdt_init  ,        &  !: total incoming heat for ice melt
+         zfdt_final ,        &  !: total remaing heat for ice melt
+         zqt_i      ,        &  !: total ice heat content
+         zqt_s      ,        &  !: total snow heat content
+         zqt_dummy              !: dummy heat content
+
+      REAL(wp), DIMENSION(jpij,jkmax) :: &
+         zqt_i_lay              !: total ice heat content
+
+      ! Heat conservation 
+      REAL(wp), DIMENSION(jpij) :: &
+         zfbase,             &
+         zdq_i
+
+      INTEGER, DIMENSION(jpij) ::  &
+         innermelt
+
+      REAL(wp) :: &
+         meance_dh
+
+      INTEGER ::                   &
+         num_iter_max,       &
+         numce_dh
 
       zfsalt_melt(:)  = 0.0
@@ -191,19 +186,19 @@
          old_ht_s_b(ji) = ht_s_b(ji)
       END DO
-!
-!------------------------------------------------------------------------------!
-!  1) Calculate available heat for surface ablation                            !
-!------------------------------------------------------------------------------!
-!
+      !
+      !------------------------------------------------------------------------------!
+      !  1) Calculate available heat for surface ablation                            !
+      !------------------------------------------------------------------------------!
+      !
       DO ji = kideb,kiut
          isnow         = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) )
          ztfs(ji)      = isnow * rtt + ( 1.0 - isnow ) * rtt
          z_f_surf(ji)  = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) *                 & 
-                         qsr_ice_1d(ji) - fc_su(ji)
+            qsr_ice_1d(ji) - fc_su(ji)
          z_f_surf(ji)  = MAX( zzero , z_f_surf(ji) ) *                        &
-                         MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) )
+            MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) )
          zfdt_init(ji) = ( z_f_surf(ji) + &
-                       MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) )  &
-                       * rdt_ice
+            MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) )  &
+            * rdt_ice
       END DO ! ji
 
@@ -212,9 +207,9 @@
       dsm_i_se_1d(:) = 0.0     
       dsm_i_si_1d(:) = 0.0     
-!
-!------------------------------------------------------------------------------!
-!  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
-!------------------------------------------------------------------------------!
-!
+      !
+      !------------------------------------------------------------------------------!
+      !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
+      !------------------------------------------------------------------------------!
+      !
       ! Layer thickness
       DO ji = kideb,kiut
@@ -239,9 +234,9 @@
          END DO
       END DO
-!
-!------------------------------------------------------------------------------|
-!  3) Surface ablation and sublimation                                         |
-!------------------------------------------------------------------------------|
-!
+      !
+      !------------------------------------------------------------------------------|
+      !  3) Surface ablation and sublimation                                         |
+      !------------------------------------------------------------------------------|
+      !
       !-------------------------
       ! 3.1 Snow precips / melt
@@ -272,5 +267,5 @@
          zdeltah(ji,1)  =  MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
          zqfont_su(ji)  =  MAX( 0.0 , - zdh_s_pre(ji) - zdeltah(ji,1) )      * &
-                           zqprec(ji)
+            zqprec(ji)
          zdeltah(ji,1)  =  MAX( - zdh_s_pre(ji) , zdeltah(ji,1) )
          zdh_s_mel(ji)  =  zdh_s_mel(ji) + zdeltah(ji,1)
@@ -289,5 +284,5 @@
             zdeltah(ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk)
             zqfont_su(ji)  =  MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * &
-                              q_s_b(ji,jk) 
+               q_s_b(ji,jk) 
             zdeltah(ji,jk) =  MAX( zdeltah(ji,jk) , - zh_s(ji) )
             zdh_s_mel(ji)  =  zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow    
@@ -306,5 +301,5 @@
          ! Volume and mass variations of snow
          dvsbq_1d(ji)   =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji)              &
-                        - zdh_s_mel(ji) )
+            - zdh_s_mel(ji) )
          dvsbq_1d(ji)   =  MIN( zzero, dvsbq_1d(ji) )
          rdmsnif_1d(ji) =  rhosn*dvsbq_1d(ji)
@@ -327,5 +322,5 @@
             ! recompute heat available
             zqfont_su(ji)       = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) *  &
-                                  q_i_b(ji,jk) 
+               q_i_b(ji,jk) 
             ! melt of layer jk cannot be higher than its thickness
             zdeltah(ji,jk)      = MAX( zdeltah(ji,jk) , - zh_i(ji) )
@@ -334,12 +329,12 @@
             ! for energy conservation
             zdq_i(ji)           = zdq_i(ji) + zdeltah(ji,jk) *                &
-                                        q_i_b(ji,jk) / rdt_ice
+               q_i_b(ji,jk) / rdt_ice
             ! contribution to ice-ocean salt flux 
             zji                 = MOD( npb(ji) - 1, jpi ) + 1
             zjj                 = ( npb(ji) - 1 ) / jpi + 1
             zfsalt_melt(ji)     = zfsalt_melt(ji) +                           &
-                                  ( sss_m(zji,zjj) - sm_i_b(ji)   ) *         &
-                                  a_i_b(ji) *                                 &
-                                  MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
+               ( sss_m(zji,zjj) - sm_i_b(ji)   ) *         &
+               a_i_b(ji) *                                 &
+               MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
          END DO ! ji
       END DO ! jk
@@ -349,34 +344,34 @@
       !-------------------
       IF ( con_i ) THEN
-      numce_dh = 0
-      meance_dh = 0.0
-      DO ji = kideb, kiut
-
-         IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
-            numce_dh  = numce_dh + 1
-            meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)
-         ENDIF
-
-         IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
-            WRITE(numout,*) ' ALERTE heat loss for surface melt '
-            WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
-            WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
-            WRITE(numout,*) ' z_f_surf  : ', z_f_surf(ji)
-            WRITE(numout,*) ' zdq_i   : ', zdq_i(ji)
-            WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
-            WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
-            WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
-            WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
-            WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
-            WRITE(numout,*) ' sss_m   : ', sss_m(zji,zjj)
-         ENDIF
-
-      END DO ! ji
-
-      IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
-      WRITE(numout,*) ' Error report - Category : ', jl
-      WRITE(numout,*) ' ~~~~~~~~~~~~ '
-      WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh
-      WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh
+         numce_dh = 0
+         meance_dh = 0.0
+         DO ji = kideb, kiut
+
+            IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
+               numce_dh  = numce_dh + 1
+               meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)
+            ENDIF
+
+            IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
+               WRITE(numout,*) ' ALERTE heat loss for surface melt '
+               WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
+               WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
+               WRITE(numout,*) ' z_f_surf  : ', z_f_surf(ji)
+               WRITE(numout,*) ' zdq_i   : ', zdq_i(ji)
+               WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
+               WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
+               WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
+               WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
+               WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
+               WRITE(numout,*) ' sss_m   : ', sss_m(zji,zjj)
+            ENDIF
+
+         END DO ! ji
+
+         IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
+         WRITE(numout,*) ' Error report - Category : ', jl
+         WRITE(numout,*) ' ~~~~~~~~~~~~ '
+         WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh
+         WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh
 
       ENDIF ! con_i
@@ -409,17 +404,17 @@
       DO jk = 1, nlay_s !n 
          DO ji = kideb, kiut !n
-         ! In case of disparition of the snow, we have to update the snow 
-         ! temperatures
+            ! In case of disparition of the snow, we have to update the snow 
+            ! temperatures
             zhisn  =  MAX( zzero , SIGN( zone, - ht_s_b(ji) ) )
             t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
             q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
          END DO
-      END DO 
-
-!
-!------------------------------------------------------------------------------!
-! 4) Basal growth / melt                                                       !
-!------------------------------------------------------------------------------!
-!
+      END DO
+
+      !
+      !------------------------------------------------------------------------------!
+      ! 4) Basal growth / melt                                                       !
+      !------------------------------------------------------------------------------!
+      !
       ! Ice basal growth / melt is given by the ratio of heat budget over basal
       ! ice heat content.  Basal heat budget is given by the difference between
@@ -439,13 +434,13 @@
                ! New ice heat content (Bitz and Lipscomb, 1999)
                ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the
-                                                        ! Baltic
+               ! Baltic
                q_i_b(ji,nlay_i+1)  =  rhoic * &
-                                      ( cpic * ( ztmelts - ztform     )       &
-                                      + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
-                                                 ( ztform     - rtt ) )       &
-                                      - rcp * ( ztmelts-rtt ) )
+                  ( cpic * ( ztmelts - ztform     )       &
+                  + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
+                  ( ztform     - rtt ) )       &
+                  - rcp * ( ztmelts-rtt ) )
                ! Basal growth rate = - F*dt / q
                dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
-                                      qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
+                  qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
             ENDIF ! heat budget
          END DO ! ji
@@ -476,11 +471,11 @@
                   ! New ice heat content (Bitz and Lipscomb, 1999)
                   q_i_b(ji,nlay_i+1)  =  rhoic * &
-                                      ( cpic * ( ztmelts - t_bo_b(ji) )       &
-                                      + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
-                                                 ( t_bo_b(ji) - rtt ) )       &
-                                      - rcp * ( ztmelts-rtt ) )
+                     ( cpic * ( ztmelts - t_bo_b(ji) )       &
+                     + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
+                     ( t_bo_b(ji) - rtt ) )       &
+                     - rcp * ( ztmelts-rtt ) )
                   ! Bottom growth rate = - F*dt / q
                   dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) &
-                                      + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
+                     + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
                   ! New ice salinity ( Cox and Weeks, JGR, 1988 )
                   ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7
@@ -492,7 +487,7 @@
                   zswi1  = 1. - zswi2 * zswi12 
                   zfracs = zswi1  * 0.12 +  &
-                           zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + &
-                           zswi2  * 0.26 /  &
-                           ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
+                     zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + &
+                     zswi2  * 0.26 /  &
+                     ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
                   zds         = zfracs*sss_m(zji,zjj) - s_i_new(ji)
                   s_i_new(ji) = zfracs * sss_m(zji,zjj)
@@ -510,17 +505,17 @@
                ! New ice heat content (Bitz and Lipscomb, 1999)
                q_i_b(ji,nlay_i+1)  =  rhoic *                              &
-                                   ( cpic * ( ztmelts - t_bo_b(ji) )       &
-                                   + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
-                                              ( t_bo_b(ji) - rtt ) )       &
-                                   - rcp * ( ztmelts-rtt ) )
+                  ( cpic * ( ztmelts - t_bo_b(ji) )       &
+                  + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
+                  ( t_bo_b(ji) - rtt ) )       &
+                  - rcp * ( ztmelts-rtt ) )
                ! Basal growth rate = - F*dt / q
                dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
-                                      qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
+                  qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
                ! Salinity update
                ! entrapment during bottom growth
                dsm_i_se_1d(ji) =  ( s_i_new(ji)*dh_i_bott(ji) +              & 
-                                   sm_i_b(ji)*ht_i_b(ji) ) /                 & 
-                                   MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps )   &
-                                   - sm_i_b(ji)
+                  sm_i_b(ji)*ht_i_b(ji) ) /                 & 
+                  MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps )   &
+                  - sm_i_b(ji)
             ENDIF ! heat budget
          END DO ! ji
@@ -537,5 +532,5 @@
          ! heat convergence at the surface > 0
          IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
-                 
+
             s_i_new(ji)   =  s_i_b(ji,nlay_i)
             zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
@@ -559,16 +554,16 @@
                   zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk)
                   zqfont_bo(ji)   = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * &
-                                    q_i_b(ji,jk)
+                     q_i_b(ji,jk)
                   zdeltah(ji,jk)  = MAX(zdeltah(ji,jk), - zh_i(ji) )
                   dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
                   zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * &
-                                    q_i_b(ji,jk) / rdt_ice
-               ! contribution to salt flux
+                     q_i_b(ji,jk) / rdt_ice
+                  ! contribution to salt flux
                   zji             = MOD( npb(ji) - 1, jpi ) + 1
                   zjj             = ( npb(ji) - 1 ) / jpi + 1
                   zfsalt_melt(ji) = zfsalt_melt(ji) +                         &
-                                   ( sss_m(zji,zjj) - sm_i_b(ji)   ) *        &
-                                   a_i_b(ji) * &
-                                   MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
+                     ( sss_m(zji,zjj) - sm_i_b(ji)   ) *        &
+                     a_i_b(ji) * &
+                     MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
                ENDIF
             ENDIF
@@ -580,40 +575,40 @@
       !-------------------
       IF ( con_i ) THEN
-      DO ji = kideb, kiut
-         IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
-            IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
-               numce_dh = numce_dh + 1
-               meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
-            ENDIF
-            IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
-                WRITE(numout,*) ' ALERTE heat loss for basal  melt '
-                WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
-                WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
-                WRITE(numout,*) ' zfbase  : ', zfbase(ji)
-                WRITE(numout,*) ' zdq_i   : ', zdq_i(ji)
-                WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
-                WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
-                WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
-                WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
-                WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
-                WRITE(numout,*) ' sss_m   : ', sss_m(zji,zjj)
-                WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
-                WRITE(numout,*) ' innermelt : ', innermelt(ji)
-            ENDIF
-         ENDIF ! heat convergence at the surface
-      END DO ! ji
-
-      IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
-      WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh
-      WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh
-      WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)
+         DO ji = kideb, kiut
+            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
+               IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
+                  numce_dh = numce_dh + 1
+                  meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
+               ENDIF
+               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
+                  WRITE(numout,*) ' ALERTE heat loss for basal  melt '
+                  WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
+                  WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
+                  WRITE(numout,*) ' zfbase  : ', zfbase(ji)
+                  WRITE(numout,*) ' zdq_i   : ', zdq_i(ji)
+                  WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
+                  WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
+                  WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
+                  WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
+                  WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
+                  WRITE(numout,*) ' sss_m   : ', sss_m(zji,zjj)
+                  WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
+                  WRITE(numout,*) ' innermelt : ', innermelt(ji)
+               ENDIF
+            ENDIF ! heat convergence at the surface
+         END DO ! ji
+
+         IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
+         WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh
+         WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh
+         WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)
 
       ENDIF ! con_i
 
-!
-!------------------------------------------------------------------------------!
-!  5) Pathological cases                                                       !
-!------------------------------------------------------------------------------!
-!
+      !
+      !------------------------------------------------------------------------------!
+      !  5) Pathological cases                                                       !
+      !------------------------------------------------------------------------------!
+      !
       !----------------------------------------------
       ! 5.1 Excessive ablation in a 1-category model
@@ -626,5 +621,5 @@
          ! excessive energy is sent to lateral ablation
          fsup(ji)            =  rhoic*lfus * at_i_b(ji) / MAX( ( 1.0 - at_i_b(ji) ),epsi13) &
-                             *  ( zdhbf - dh_i_bott(ji) ) / rdt_ice
+            *  ( zdhbf - dh_i_bott(ji) ) / rdt_ice
 
          dh_i_bott(ji)  = zdhbf
@@ -638,9 +633,9 @@
          zjj                 = ( npb(ji) - 1 ) / jpi + 1
          diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) &
-                              / rdt_ice
+            / rdt_ice
          diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) &
-                              / rdt_ice
+            / rdt_ice
          diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) &
-                              / rdt_ice
+            / rdt_ice
       END DO
 
@@ -667,5 +662,5 @@
          zqt_s(ji)  =  ( 1. - zihg) * zqt_s(ji) / MAX( zhni, zeps )
          zdhnm      =  - ( 1. - zihg ) * ( 1. - zihgnew ) * ( zfdt_final(ji) /  &
-                           MAX( zqt_s(ji) , zeps ) )
+            MAX( zqt_s(ji) , zeps ) )
          zhnfi          =  zhni + zdhnm
          zfdt_final(ji) =  MAX ( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 )
@@ -676,8 +671,8 @@
          !---------------------------------
          rdmicif_1d(ji) =  rdmicif_1d(ji) + a_i_b(ji) *                        &
-                           (zhgnew(ji)-ht_i_b(ji))*rhoic ! good
+            (zhgnew(ji)-ht_i_b(ji))*rhoic ! good
 
          rdmsnif_1d(ji) =  rdmsnif_1d(ji) + a_i_b(ji) * &
-                           (ht_s_b(ji)-zhni)*rhosn ! good too
+            (ht_s_b(ji)-zhni)*rhosn ! good too
 
          ! Remaining heat to the ocean 
@@ -700,18 +695,18 @@
          zjj           = ( npb(ji) - 1 ) / jpi + 1
          IF ( num_sal .NE. 4 ) &
-         fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
-                          (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
-                          ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice
+            fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
+            (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
+            ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice
          ! new lines
          IF ( num_sal .EQ. 4 ) &
-         fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
-                          (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
-                          ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice
+            fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
+            (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
+            ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice
          ! Heat flux
          ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
          ! excessive total ablation energy (focea) sent to the ocean
          qfvbq_1d(ji)  = qfvbq_1d(ji) + &
-                         fsup(ji) + ( 1.0 - zihgnew ) *        & 
-                         focea(ji) * a_i_b(ji) * rdt_ice
+            fsup(ji) + ( 1.0 - zihgnew ) *        & 
+            focea(ji) * a_i_b(ji) * rdt_ice
 
          zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) )
@@ -719,7 +714,7 @@
          fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
          qldif_1d(ji)  = qldif_1d(ji)                                         &
-                       + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & 
-                         * rdt_ice                               &
-                       + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice
+            + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & 
+            * rdt_ice                               &
+            + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice
       END DO  ! ji
 
@@ -743,9 +738,9 @@
          ht_i_b(ji) = zhgnew(ji)
       END DO  ! ji
-!
-!------------------------------------------------------------------------------|
-!  6) Snow-Ice formation                                                       |
-!------------------------------------------------------------------------------|
-!
+      !
+      !------------------------------------------------------------------------------|
+      !  6) Snow-Ice formation                                                       |
+      !------------------------------------------------------------------------------|
+      !
       ! When snow load excesses Archimede's limit, snow-ice interface goes down
       ! under sea-level, flooding of seawater transforms snow into ice
@@ -754,51 +749,51 @@
 
          dh_snowice(ji) = MAX(zzero,(rhosn*ht_s_b(ji)+(rhoic-rau0) &
-                            * ht_i_b(ji))/(rhosn+rau0-rhoic))
+            * ht_i_b(ji))/(rhosn+rau0-rhoic))
          zhgnew(ji)     = MAX(zhgnew(ji),zhgnew(ji)+dh_snowice(ji))
          zhnnew            = MIN(ht_s_b(ji),ht_s_b(ji)-dh_snowice(ji))
 
-      !  Changes in ice volume and ice mass.
+         !  Changes in ice volume and ice mass.
          dvnbq_1d(ji)      = a_i_b(ji) * (zhgnew(ji)-ht_i_b(ji))
          dmgwi_1d(ji)      = dmgwi_1d(ji) + a_i_b(ji) &
-                             *(ht_s_b(ji)-zhnnew)*rhosn
+            *(ht_s_b(ji)-zhnnew)*rhosn
 
          rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & 
-                                         * ( zhgnew(ji) - ht_i_b(ji) )*rhoic 
+            * ( zhgnew(ji) - ht_i_b(ji) )*rhoic 
          rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) & 
-                                         * ( zhnnew       - ht_s_b(ji) )*rhosn
-
-!        Equivalent salt flux (1) Snow-ice formation component
-!        -----------------------------------------------------
+            * ( zhnnew       - ht_s_b(ji) )*rhosn
+
+         !        Equivalent salt flux (1) Snow-ice formation component
+         !        -----------------------------------------------------
          zji                 = MOD( npb(ji) - 1, jpi ) + 1
          zjj                 = ( npb(ji) - 1 ) / jpi + 1
 
          zsm_snowice  = ( rhoic - rhosn ) / rhoic *            &
-                        sss_m(zji,zjj) 
+            sss_m(zji,zjj) 
 
          IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji)
 
          IF ( num_sal .NE. 4 ) &
-         fseqv_1d(ji)   = fseqv_1d(ji)   + &
-                          ( sss_m(zji,zjj) - zsm_snowice ) * &
-                            a_i_b(ji)   * &
-                          ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
+            fseqv_1d(ji)   = fseqv_1d(ji)   + &
+            ( sss_m(zji,zjj) - zsm_snowice ) * &
+            a_i_b(ji)   * &
+            ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
          ! new lines
          IF ( num_sal .EQ. 4 ) &
-         fseqv_1d(ji)   = fseqv_1d(ji)   + &
-                          ( sss_m(zji,zjj) - bulk_sal    ) * &
-                            a_i_b(ji)   * &
-                          ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
+            fseqv_1d(ji)   = fseqv_1d(ji)   + &
+            ( sss_m(zji,zjj) - bulk_sal    ) * &
+            a_i_b(ji)   * &
+            ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
 
          ! entrapment during snow ice formation
          i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-6 ) )
          isnowic      = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * &
-                        i_ice_switch
+            i_ice_switch
          IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
-             dsm_i_si_1d(ji)  = ( zsm_snowice*dh_snowice(ji) &
-                            + sm_i_b(ji) * ht_i_b(ji)                          & 
-                            / MAX( ht_i_b(ji) + dh_snowice(ji), zeps)          &
-                            - sm_i_b(ji) ) * isnowic     
-
-!  Actualize new snow and ice thickness.
+            dsm_i_si_1d(ji)  = ( zsm_snowice*dh_snowice(ji) &
+            + sm_i_b(ji) * ht_i_b(ji)                          & 
+            / MAX( ht_i_b(ji) + dh_snowice(ji), zeps)          &
+            - sm_i_b(ji) ) * isnowic     
+
+         !  Actualize new snow and ice thickness.
          ht_s_b(ji)  = zhnnew
          ht_i_b(ji)  = zhgnew(ji)
@@ -811,9 +806,9 @@
          zjj                 = ( npb(ji) - 1 ) / jpi + 1
          diag_sni_gr(zji,zjj)  = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / &
-                                 rdt_ice
+            rdt_ice
 
       END DO !ji
 
-    END SUBROUTINE lim_thd_dh
+   END SUBROUTINE lim_thd_dh
 #else
    !!======================================================================
@@ -825,3 +820,3 @@
    END SUBROUTINE lim_thd_dh
 #endif
- END MODULE limthd_dh
+END MODULE limthd_dh
