Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8426)
@@ -177,5 +177,5 @@
 
    !                                     !!** ice-init namelist (namiceini) **
-                                          ! -- iceistate -- !
+                                          ! -- iceist -- !
    LOGICAL , PUBLIC ::   ln_limini        ! initialization or not
    LOGICAL , PUBLIC ::   ln_limini_file   ! Ice initialization state from 2D netcdf file
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8426)
@@ -29,5 +29,4 @@
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE timing         ! Timing
-   USE icecons        ! conservation tests
    USE icectl         ! control prints
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8426)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8426)
@@ -0,0 +1,363 @@
+MODULE icealb
+   !!======================================================================
+   !!                       ***  MODULE  icealb  ***
+   !! Ocean forcing:  bulk thermohaline forcing of the ice
+   !!=====================================================================
+   !! History :
+   !!   NEMO     4.0  ! 2017-07  (C. Rousset) Split ice and ocean albedos
+   !!----------------------------------------------------------------------
+   !!   ice_alb    : albedo for ice (clear and overcast skies)
+   !!   alb_init   : initialisation of albedo computation
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                  LIM 3.0 sea-ice model
+   !!----------------------------------------------------------------------
+   USE ice, ONLY : jpl
+   USE phycst         ! physical constants
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_alb   ! routine called icestp.F90
+
+   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
+   INTEGER  ::   albd_init = 0       ! control flag for initialization  
+   REAL(wp) ::   c1        = 0.05    ! snow thickness (only for nn_ice_alb=0)
+   REAL(wp) ::   c2        = 0.10    !  "        "
+   REAL(wp) ::   rcloud    = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
+ 
+   !                             !!* namelist namsbc_alb
+   INTEGER  ::   nn_ice_alb
+   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2010)
+   !! $Id: icealb.F90 8268 2017-07-03 15:01:04Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_alb( pt_ice, ph_ice, ph_snw, pafrac_pnd, ph_pnd, ld_pnd, pa_ice_cs, pa_ice_os )
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE ice_alb  ***
+      !!          
+      !! ** Purpose :   Computation of the albedo of the snow/ice system 
+      !!       
+      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
+      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
+      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
+      !!                                                                           and Grenfell & Perovich (JGR 2004)
+      !!                  2: fractional surface-based formulation of scheme 1 (NEW)
+      !!                Description of scheme 1:
+      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
+      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
+      !!                     0-5cm  : linear function of ice thickness
+      !!                     5-150cm: log    function of ice thickness
+      !!                     > 150cm: constant
+      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004)
+      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting
+      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004)
+      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law
+      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
+      !!
+      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
+      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
+      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
+      !!                     under melting conditions than under freezing conditions
+      !!                  3) the evolution of ice albedo as a function of ice thickness shows  
+      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
+      !!
+      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !!                Brandt et al. 2005, J. Climate, vol 18
+      !!                Grenfell & Perovich 2004, JGR, vol 109 
+      !!
+      !!----------------------------------------------------------------------
+      !!
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice              !  ice surface temperature (Kelvin)
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice              !  sea-ice thickness
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw              !  snow depth
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd          !  melt pond relative fraction (per unit ice area)
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd              !  melt pond depth
+      LOGICAL , INTENT(in   )                   ::   ld_pnd              !  melt ponds radiatively active or not
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs           !  albedo of ice under clear    sky
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os           !  albedo of ice under overcast sky
+      !
+      INTEGER  ::   ji, jj, jl                                           ! dummy loop indices
+      REAL(wp) ::   zswitch, z1_c1, z1_c2
+      REAL(wp) ::   zhref_pnd                                 
+      REAL(wp) ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
+      !
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free)
+!! MV MP
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_pnd                  ! ponded sea ice albedo
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_ice                  ! bare sea ice albedo
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_snw                  ! snow-covered sea ice albedo
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_snw                ! relative snow fraction
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_ice                ! relative ice fraction
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_pnd                ! relative ice fraction (effective)
+      !!
+      !!---------------------------------------------------------------------
+
+      IF( albd_init == 0 )   CALL alb_init      ! initialization 
+
+      !-----------------------------------------------------
+      !  Snow-free albedo (no ice thickness dependence yet)
+      !-----------------------------------------------------
+      !
+      ! Part common to nn_ice_alb = 0, 1, 2
+      !
+      IF ( .NOT. ld_pnd ) THEN   !--- No melt ponds OR radiatively inactive melt ponds
+         ! Bare ice albedo is prescribed, with implicit assumption on pond fraction and depth
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = rn_alb_imlt
+                                                       ! !!! MV I think we could replace rt0_ice by rt0 and get rid of rt0
+         ELSE WHERE                                              ;   zalb(:,:,:) = rn_alb_idry
+         END  WHERE
+      ENDIF
+
+      SELECT CASE ( nn_ice_alb )
+
+      !------------------------------------------
+      !  Shine and Henderson-Sellers (1985)
+      !------------------------------------------
+      ! NB: This parameterization is based on clear sky values
+
+      CASE( 0 )
+       
+         ! Thickness-dependent bare ice albedo
+         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
+         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              &
+            &                                                                 - 0.8608 * ph_ice * ph_ice     &
+            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
+         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
+         END WHERE
+
+         IF ( ld_pnd ) THEN
+            ! Depth-dependent ponded ice albedo
+            zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
+            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd / zhref_pnd ) 
+
+            ! Snow-free ice albedo is a function of pond fraction
+            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ; zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd ;   END WHERE
+         ENDIF 
+
+         DO jl = 1, jpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! Freezing snow
+                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
+                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
+                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
+                     &                           + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) / c1  )   &
+                     &        +         zswitch   * rn_alb_sdry  
+
+                  ! Melting snow
+                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
+                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
+                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) / c2 )   &
+                      &     +         zswitch   *   rn_alb_smlt 
+                  !
+                  ! Snow albedo
+                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+               
+                  ! Surface albedo
+                  zswitch             = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_cs(ji,jj,jl) = zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
+                  !
+               END DO
+            END DO
+         END DO
+
+         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
+
+      !------------------------------------------
+      !  New parameterization (2016)
+      !------------------------------------------
+      ! NB: This parameterization is based on overcast skies values
+      
+      CASE( 1 ) 
+
+! compilation of values from literature (reference overcast sky values)
+!        rn_alb_sdry = 0.85      ! dry snow
+!        rn_alb_smlt = 0.75      ! melting snow
+!        rn_alb_idry = 0.60      ! bare frozen ice
+!        rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015)
+!                                ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich
+! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
+!        rn_alb_sdry = 0.85      ! dry snow
+!        rn_alb_smlt = 0.72      ! melting snow
+!        rn_alb_idry = 0.65      ! bare frozen ice
+! Brandt et al 2005 (East Antarctica)
+!        rn_alb_sdry = 0.87      ! dry snow
+!        rn_alb_smlt = 0.82      ! melting snow
+!        rn_alb_idry = 0.54      ! bare frozen ice
+! 
+         ! Computation of snow-free ice albedo 
+         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+         z1_c2 = 1. / 0.05
+
+         ! Accounting for the ice-thickness dependency
+         WHERE     ( 1.5  < ph_ice                     )        ;  zalb_it = zalb
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )        ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
+            &                                                                     ( LOG(1.5) - LOG(ph_ice) )
+         ELSE WHERE                                             ;  zalb_it = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
+         END WHERE
+
+         IF ( ld_pnd ) THEN
+            ! Depth-dependent ponded ice albedo
+            zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
+            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd / zhref_pnd ) 
+
+            ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions
+            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;  zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd ;  END WHERE
+         ENDIF 
+
+         z1_c1 = 1. / 0.02
+         z1_c2 = 1. / 0.03
+         
+         ! Overcast sky surface albedo (accounting for snow, ice melt ponds)
+         DO jl = 1, jpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! Snow depth dependence of snow albedo
+                  zalb_sf = rn_alb_sdry - ( rn_alb_sdry - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
+                  zalb_sm = rn_alb_smlt - ( rn_alb_smlt - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
+
+                  ! Snow albedo (MV I guess we could use rt0 instead of rt0_snow)
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+
+                  ! Surface albedo   
+                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
+
+              END DO
+            END DO
+         END DO
+
+         ! Clear sky surface albedo
+         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
+
+      !---------------------------------------------------
+      !  Fractional surface-based parameterization (2017)
+      !---------------------------------------------------
+      CASE( 2 ) 
+ 
+      ! MV: I propose this formulation that is more elegant, and more easy to expand towards
+      !     varying pond and snow fraction.
+      !     Formulation 1 has issues to handle ponds and snow together that
+      !     can't easily be fixed. This one handles it better, I believe.
+
+          !-----------------------------------------
+          ! Snow, bare ice and ponded ice fractions 
+          !-----------------------------------------
+          ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category
+
+          !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time)
+          zafrac_pnd = 0._wp
+          IF ( ld_pnd ) THEN  
+             WHERE( ph_snw == 0._wp ) ;  zafrac_pnd = pafrac_pnd ;  END WHERE  ! Snow fully "shades" melt ponds
+          ENDIF         
+
+          !--- Specific snow fraction (for now, prescribed)
+          WHERE     ( ph_snw > 0._wp     ) ;  zafrac_snw = 1.
+          ELSE WHERE                       ;  zafrac_snw = 0.
+          END WHERE
+ 
+          !--- Specific ice fraction 
+          zafrac_ice = 1. - zafrac_snw - zafrac_pnd
+ 
+          !--------------------------------------------------
+          ! Snow-covered, pond-covered, and bare ice albedos
+          !--------------------------------------------------
+          ! Bare ice albedo
+          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+          z1_c2 = 1. / 0.05
+          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_ice = zalb
+          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_ice = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
+            &                                                                       ( LOG(1.5) - LOG(ph_ice) )
+          ELSE WHERE                                       ;  zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
+          END WHERE
+
+          ! Snow-covered ice albedo (freezing, melting cases)
+          z1_c1 = 1. / 0.02
+          z1_c2 = 1. / 0.03
+         
+          WHERE( pt_ice < rt0_snow ) ; zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 );
+          ELSE WHERE                 ; zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 );
+          END WHERE
+
+          ! Depth-dependent ponded ice albedo
+          IF ( ld_pnd ) THEN
+             zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
+             zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd / zhref_pnd ) 
+          ELSE
+             zalb_pnd  = rn_alb_dpnd
+          ENDIF
+
+          ! Surface albedo is weighted mean of snow, ponds and bare ice contributions
+          pa_ice_os = zafrac_snw * zalb_snw  +  zafrac_pnd * zalb_pnd  +  zafrac_ice * zalb_ice
+          
+          pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 )
+
+      END SELECT
+      !
+   END SUBROUTINE ice_alb
+
+   SUBROUTINE alb_init
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE alb_init  ***
+      !!
+      !! ** Purpose :   initializations for the albedo parameters
+      !!
+      !! ** Method  :   Read the namelist namicealb
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
+      NAMELIST/namicealb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
+      !!----------------------------------------------------------------------
+      !
+      albd_init = 1                     ! indicate that the initialization has been done
+      !
+      REWIND( numnam_ice_ref )              ! Namelist namicealb in reference namelist : Albedo parameters
+      READ  ( numnam_ice_ref, namicealb, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicealb in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namsbc_alb in configuration namelist : Albedo parameters
+      READ  ( numnam_ice_cfg, namicealb, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicealb in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namicealb )
+      !
+      IF(lwp) THEN                      ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'albedo : set albedo parameters'
+         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) '   Namelist namicealb : albedo '
+         WRITE(numout,*) '      choose the albedo parameterization   nn_ice_alb  = ', nn_ice_alb
+         WRITE(numout,*) '      albedo of dry snow                   rn_alb_sdry = ', rn_alb_sdry
+         WRITE(numout,*) '      albedo of melting snow               rn_alb_smlt = ', rn_alb_smlt
+         WRITE(numout,*) '      albedo of dry ice                    rn_alb_idry = ', rn_alb_idry
+         WRITE(numout,*) '      albedo of bare puddled ice           rn_alb_imlt = ', rn_alb_imlt
+         WRITE(numout,*) '      albedo of ponded ice                 rn_alb_dpnd = ', rn_alb_dpnd
+      ENDIF
+      !
+   END SUBROUTINE alb_init
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   !
+   SUBROUTINE ice_alb     ! Dummy routine
+      WRITE(*,*) 'ice_alb: You should not have seen this print! error?'
+   END SUBROUTINE ice_alb
+#endif
+
+   !!======================================================================
+END MODULE icealb
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealbedo.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealbedo.F90	(revision 8424)
+++ 	(revision )
@@ -1,363 +1,0 @@
-MODULE icealbedo
-   !!======================================================================
-   !!                       ***  MODULE  icealbedo  ***
-   !! Ocean forcing:  bulk thermohaline forcing of the ice
-   !!=====================================================================
-   !! History :
-   !!   NEMO     4.0  ! 2017-07  (C. Rousset) Split ice and ocean albedos
-   !!----------------------------------------------------------------------
-   !!   ice_albedo    : albedo for ice (clear and overcast skies)
-   !!   albedo_init   : initialisation of albedo computation
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                  LIM 3.0 sea-ice model
-   !!----------------------------------------------------------------------
-   USE ice, ONLY : jpl
-   USE phycst         ! physical constants
-   !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_albedo   ! routine called icestp.F90
-
-   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
-   INTEGER  ::   albd_init = 0       ! control flag for initialization  
-   REAL(wp) ::   c1        = 0.05    ! snow thickness (only for nn_ice_alb=0)
-   REAL(wp) ::   c2        = 0.10    !  "        "
-   REAL(wp) ::   rcloud    = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
- 
-   !                             !!* namelist namsbc_alb
-   INTEGER  ::   nn_ice_alb
-   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2010)
-   !! $Id: icealbedo.F90 8268 2017-07-03 15:01:04Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_albedo( pt_ice, ph_ice, ph_snw, pafrac_pnd, ph_pnd, ld_pnd, pa_ice_cs, pa_ice_os )
-      !!----------------------------------------------------------------------
-      !!               ***  ROUTINE ice_albedo  ***
-      !!          
-      !! ** Purpose :   Computation of the albedo of the snow/ice system 
-      !!       
-      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
-      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
-      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
-      !!                                                                           and Grenfell & Perovich (JGR 2004)
-      !!                  2: fractional surface-based formulation of scheme 1 (NEW)
-      !!                Description of scheme 1:
-      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
-      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
-      !!                     0-5cm  : linear function of ice thickness
-      !!                     5-150cm: log    function of ice thickness
-      !!                     > 150cm: constant
-      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004)
-      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting
-      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004)
-      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law
-      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
-      !!
-      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
-      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
-      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
-      !!                     under melting conditions than under freezing conditions
-      !!                  3) the evolution of ice albedo as a function of ice thickness shows  
-      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
-      !!
-      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
-      !!                Brandt et al. 2005, J. Climate, vol 18
-      !!                Grenfell & Perovich 2004, JGR, vol 109 
-      !!
-      !!----------------------------------------------------------------------
-      !!
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice              !  ice surface temperature (Kelvin)
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice              !  sea-ice thickness
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw              !  snow depth
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd          !  melt pond relative fraction (per unit ice area)
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd              !  melt pond depth
-      LOGICAL , INTENT(in   )                   ::   ld_pnd              !  melt ponds radiatively active or not
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs           !  albedo of ice under clear    sky
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os           !  albedo of ice under overcast sky
-      !
-      INTEGER  ::   ji, jj, jl                                           ! dummy loop indices
-      REAL(wp) ::   zswitch, z1_c1, z1_c2
-      REAL(wp) ::   zhref_pnd                                 
-      REAL(wp) ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
-      !
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free)
-!! MV MP
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_pnd                  ! ponded sea ice albedo
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_ice                  ! bare sea ice albedo
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_snw                  ! snow-covered sea ice albedo
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_snw                ! relative snow fraction
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_ice                ! relative ice fraction
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_pnd                ! relative ice fraction (effective)
-      !!
-      !!---------------------------------------------------------------------
-
-      IF( albd_init == 0 )   CALL albedo_init      ! initialization 
-
-      !-----------------------------------------------------
-      !  Snow-free albedo (no ice thickness dependence yet)
-      !-----------------------------------------------------
-      !
-      ! Part common to nn_ice_alb = 0, 1, 2
-      !
-      IF ( .NOT. ld_pnd ) THEN   !--- No melt ponds OR radiatively inactive melt ponds
-         ! Bare ice albedo is prescribed, with implicit assumption on pond fraction and depth
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = rn_alb_imlt
-                                                       ! !!! MV I think we could replace rt0_ice by rt0 and get rid of rt0
-         ELSE WHERE                                              ;   zalb(:,:,:) = rn_alb_idry
-         END  WHERE
-      ENDIF
-
-      SELECT CASE ( nn_ice_alb )
-
-      !------------------------------------------
-      !  Shine and Henderson-Sellers (1985)
-      !------------------------------------------
-      ! NB: This parameterization is based on clear sky values
-
-      CASE( 0 )
-       
-         ! Thickness-dependent bare ice albedo
-         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
-         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
-         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              &
-            &                                                                 - 0.8608 * ph_ice * ph_ice     &
-            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
-         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
-         END WHERE
-
-         IF ( ld_pnd ) THEN
-            ! Depth-dependent ponded ice albedo
-            zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
-            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd / zhref_pnd ) 
-
-            ! Snow-free ice albedo is a function of pond fraction
-            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ; zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd ;   END WHERE
-         ENDIF 
-
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! Freezing snow
-                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
-                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
-                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
-                     &                           + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) / c1  )   &
-                     &        +         zswitch   * rn_alb_sdry  
-
-                  ! Melting snow
-                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
-                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
-                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) / c2 )   &
-                      &     +         zswitch   *   rn_alb_smlt 
-                  !
-                  ! Snow albedo
-                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
-                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
-               
-                  ! Surface albedo
-                  zswitch             = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
-                  pa_ice_cs(ji,jj,jl) = zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
-                  !
-               END DO
-            END DO
-         END DO
-
-         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
-
-      !------------------------------------------
-      !  New parameterization (2016)
-      !------------------------------------------
-      ! NB: This parameterization is based on overcast skies values
-      
-      CASE( 1 ) 
-
-! compilation of values from literature (reference overcast sky values)
-!        rn_alb_sdry = 0.85      ! dry snow
-!        rn_alb_smlt = 0.75      ! melting snow
-!        rn_alb_idry = 0.60      ! bare frozen ice
-!        rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015)
-!                                ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich
-! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
-!        rn_alb_sdry = 0.85      ! dry snow
-!        rn_alb_smlt = 0.72      ! melting snow
-!        rn_alb_idry = 0.65      ! bare frozen ice
-! Brandt et al 2005 (East Antarctica)
-!        rn_alb_sdry = 0.87      ! dry snow
-!        rn_alb_smlt = 0.82      ! melting snow
-!        rn_alb_idry = 0.54      ! bare frozen ice
-! 
-         ! Computation of snow-free ice albedo 
-         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
-         z1_c2 = 1. / 0.05
-
-         ! Accounting for the ice-thickness dependency
-         WHERE     ( 1.5  < ph_ice                     )        ;  zalb_it = zalb
-         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )        ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
-            &                                                                     ( LOG(1.5) - LOG(ph_ice) )
-         ELSE WHERE                                             ;  zalb_it = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
-         END WHERE
-
-         IF ( ld_pnd ) THEN
-            ! Depth-dependent ponded ice albedo
-            zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
-            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd / zhref_pnd ) 
-
-            ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions
-            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;  zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd ;  END WHERE
-         ENDIF 
-
-         z1_c1 = 1. / 0.02
-         z1_c2 = 1. / 0.03
-         
-         ! Overcast sky surface albedo (accounting for snow, ice melt ponds)
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! Snow depth dependence of snow albedo
-                  zalb_sf = rn_alb_sdry - ( rn_alb_sdry - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
-                  zalb_sm = rn_alb_smlt - ( rn_alb_smlt - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
-
-                  ! Snow albedo (MV I guess we could use rt0 instead of rt0_snow)
-                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
-                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
-
-                  ! Surface albedo   
-                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
-                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
-
-              END DO
-            END DO
-         END DO
-
-         ! Clear sky surface albedo
-         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
-
-      !---------------------------------------------------
-      !  Fractional surface-based parameterization (2017)
-      !---------------------------------------------------
-      CASE( 2 ) 
- 
-      ! MV: I propose this formulation that is more elegant, and more easy to expand towards
-      !     varying pond and snow fraction.
-      !     Formulation 1 has issues to handle ponds and snow together that
-      !     can't easily be fixed. This one handles it better, I believe.
-
-          !-----------------------------------------
-          ! Snow, bare ice and ponded ice fractions 
-          !-----------------------------------------
-          ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category
-
-          !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time)
-          zafrac_pnd = 0._wp
-          IF ( ld_pnd ) THEN  
-             WHERE( ph_snw == 0._wp ) ;  zafrac_pnd = pafrac_pnd ;  END WHERE  ! Snow fully "shades" melt ponds
-          ENDIF         
-
-          !--- Specific snow fraction (for now, prescribed)
-          WHERE     ( ph_snw > 0._wp     ) ;  zafrac_snw = 1.
-          ELSE WHERE                       ;  zafrac_snw = 0.
-          END WHERE
- 
-          !--- Specific ice fraction 
-          zafrac_ice = 1. - zafrac_snw - zafrac_pnd
- 
-          !--------------------------------------------------
-          ! Snow-covered, pond-covered, and bare ice albedos
-          !--------------------------------------------------
-          ! Bare ice albedo
-          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
-          z1_c2 = 1. / 0.05
-          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_ice = zalb
-          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_ice = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
-            &                                                                       ( LOG(1.5) - LOG(ph_ice) )
-          ELSE WHERE                                       ;  zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
-          END WHERE
-
-          ! Snow-covered ice albedo (freezing, melting cases)
-          z1_c1 = 1. / 0.02
-          z1_c2 = 1. / 0.03
-         
-          WHERE( pt_ice < rt0_snow ) ; zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 );
-          ELSE WHERE                 ; zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 );
-          END WHERE
-
-          ! Depth-dependent ponded ice albedo
-          IF ( ld_pnd ) THEN
-             zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
-             zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd / zhref_pnd ) 
-          ELSE
-             zalb_pnd  = rn_alb_dpnd
-          ENDIF
-
-          ! Surface albedo is weighted mean of snow, ponds and bare ice contributions
-          pa_ice_os = zafrac_snw * zalb_snw  +  zafrac_pnd * zalb_pnd  +  zafrac_ice * zalb_ice
-          
-          pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 )
-
-      END SELECT
-      !
-   END SUBROUTINE ice_albedo
-
-   SUBROUTINE albedo_init
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE albedo_init  ***
-      !!
-      !! ** Purpose :   initializations for the albedo parameters
-      !!
-      !! ** Method  :   Read the namelist namicealb
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicealb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
-      !!----------------------------------------------------------------------
-      !
-      albd_init = 1                     ! indicate that the initialization has been done
-      !
-      REWIND( numnam_ice_ref )              ! Namelist namicealb in reference namelist : Albedo parameters
-      READ  ( numnam_ice_ref, namicealb, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicealb in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namsbc_alb in configuration namelist : Albedo parameters
-      READ  ( numnam_ice_cfg, namicealb, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicealb in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namicealb )
-      !
-      IF(lwp) THEN                      ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'albedo : set albedo parameters'
-         WRITE(numout,*) '~~~~~~~'
-         WRITE(numout,*) '   Namelist namicealb : albedo '
-         WRITE(numout,*) '      choose the albedo parameterization   nn_ice_alb  = ', nn_ice_alb
-         WRITE(numout,*) '      albedo of dry snow                   rn_alb_sdry = ', rn_alb_sdry
-         WRITE(numout,*) '      albedo of melting snow               rn_alb_smlt = ', rn_alb_smlt
-         WRITE(numout,*) '      albedo of dry ice                    rn_alb_idry = ', rn_alb_idry
-         WRITE(numout,*) '      albedo of bare puddled ice           rn_alb_imlt = ', rn_alb_imlt
-         WRITE(numout,*) '      albedo of ponded ice                 rn_alb_dpnd = ', rn_alb_dpnd
-      ENDIF
-      !
-   END SUBROUTINE albedo_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   !
-   SUBROUTINE ice_albedo     ! Dummy routine
-      WRITE(*,*) 'ice_albedo: You should not have seen this print! error?'
-   END SUBROUTINE ice_albedo
-#endif
-
-   !!======================================================================
-END MODULE icealbedo
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecons.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecons.F90	(revision 8424)
+++ 	(revision )
@@ -1,198 +1,0 @@
-MODULE icecons
-   !!======================================================================
-   !!                   ***  MODULE  icecons  ***
-   !! LIM-3 Sea Ice :   conservation check
-   !!======================================================================
-   !! History :   -   ! Original code from William H. Lipscomb, LANL
-   !!            3.0  ! 2004-06  (M. Vancoppenolle)   Energy Conservation 
-   !!            3.5  ! 2011-02  (G. Madec)  add mpp considerations
-   !!             -   ! 2014-05  (C. Rousset) add ice_cons_hsm
-   !!             -   ! 2015-03  (C. Rousset) add ice_cons_final
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM-3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!    ice_cons     :   checks whether energy, mass and salt are conserved 
-   !!----------------------------------------------------------------------
-   USE phycst         ! physical constants
-   USE ice            ! LIM-3 variables
-   USE dom_oce        ! ocean domain
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE sbc_oce , ONLY : sfx  ! Surface boundary condition: ocean fields
-   USE sbc_ice , ONLY : qevap_ice
-   
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_cons_hsm
-   PUBLIC   ice_cons_final
-
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id: icecons.F90 8409 2017-08-07 15:29:21Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
-      !!--------------------------------------------------------------------------------------------------------
-      !!                                        ***  ROUTINE ice_cons_hsm ***
-      !!
-      !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine
-      !!                     + test if ice concentration and volume are > 0
-      !!
-      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiachk=true
-      !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to
-      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
-      !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
-      !!              and a heat content of 3e5 J/kg (=latent heat of fusion) 
-      !!--------------------------------------------------------------------------------------------------------
-      INTEGER         , INTENT(in)    :: icount        ! determine wether this is the beggining of the routine (0) or the end (1)
-      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine
-      REAL(wp)        , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
-      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft
-      REAL(wp)                        :: zvmin, zamin, zamax 
-      REAL(wp)                        :: zvtrp, zetrp
-      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill
-      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
-
-      IF( icount == 0 ) THEN
-
-         ! salt flux
-         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &
-            &                ) *  e1e2t(:,:) * tmask(:,:,1) * zconv )
-
-         ! water flux
-         zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
-            &                  wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
-            &                  wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        &
-            &                ) * e1e2t(:,:) * tmask(:,:,1) * zconv )
-
-         ! heat flux
-         zft_b  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
-            &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &                ) *  e1e2t(:,:) * tmask(:,:,1) * zconv )
-
-         zvi_b  = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * tmask(:,:,1) * zconv )
-
-         zsmv_b = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * tmask(:,:,1) * zconv )
-
-         zei_b  = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
-            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    &
-                            ) * e1e2t * tmask(:,:,1) * zconv )
-
-      ELSEIF( icount == 1 ) THEN
-
-         ! salt flux
-         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
-            &              ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zfs_b
-
-         ! water flux
-         zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
-            &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
-            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        &
-            &              ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zfw_b
-
-         ! heat flux
-         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
-            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &              ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zft_b
- 
-         ! outputs
-         zvi  = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 )  &
-            &                    * e1e2t * tmask(:,:,1) * zconv ) - zvi_b ) * r1_rdtice - zfw ) * rday
-
-         zsmv = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 )  &
-            &                    * e1e2t * tmask(:,:,1) * zconv ) - zsmv_b ) * r1_rdtice + zfs ) * rday
-
-         zei  = ( glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
-            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    &
-            &                ) * e1e2t * tmask(:,:,1) * zconv ) - zei_b ) * r1_rdtice + zft
-
-         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative
-         zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t * tmask(:,:,1) * zconv ) * rday 
-         zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e1e2t * tmask(:,:,1) * zconv )
-
-         zvmin = glob_min( v_i )
-         zamax = glob_max( SUM( a_i, dim=3 ) )
-         zamin = glob_min( a_i )
-
-         ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
-         zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t * zconv ) ! in 1.e9 m2
-         zv_sill = zarea * 2.5e-5
-         zs_sill = zarea * 25.e-5
-         zh_sill = zarea * 10.e-5
-
-         IF(lwp) THEN
-            IF ( ABS( zvi  ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zvi
-            IF ( ABS( zsmv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv
-            IF ( ABS( zei  ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zei
-            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'iceadv' ) THEN
-                                         WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
-                                         WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp
-            ENDIF
-            IF (     zvmin   < -epsi10 ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
-            IF (     zamax   > MAX( rn_amax_n, rn_amax_s ) + epsi10 .AND. &
-               &                         cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' ) THEN
-                                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
-            IF (     zamax   > 1._wp   ) WRITE(numout,*) 'violation a_i>1               (',cd_routine,') = ',zamax
-            ENDIF
-            IF (      zamin  < -epsi10 ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
-         ENDIF
-
-      ENDIF
-
-   END SUBROUTINE ice_cons_hsm
-
-   SUBROUTINE ice_cons_final( cd_routine )
-      !!---------------------------------------------------------------------------------------------------------
-      !!                                   ***  ROUTINE ice_cons_final ***
-      !!
-      !! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step
-      !!
-      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiachk=true
-      !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to
-      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
-      !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
-      !!              and a heat content of 3e5 J/kg (=latent heat of fusion) 
-      !!--------------------------------------------------------------------------------------------------------
-      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine
-      REAL(wp)                        :: zhfx, zsfx, zvfx
-      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill
-      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
-
-      ! heat flux
-      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es   &
-      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check)
-         &              ) * e1e2t * tmask(:,:,1) * zconv ) 
-      ! salt flux
-      zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t * tmask(:,:,1) * zconv ) * rday
-      ! water flux
-      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t * tmask(:,:,1) * zconv ) * rday
-
-      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
-      zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t * zconv ) ! in 1.e9 m2
-      zv_sill = zarea * 2.5e-5
-      zs_sill = zarea * 25.e-5
-      zh_sill = zarea * 10.e-5
-
-      IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx    [Mt/day]       (',cd_routine,')  = ',(zvfx)
-      IF( ABS( zsfx ) > zs_sill ) WRITE(numout,*) 'violation sfx    [psu*Mt/day]   (',cd_routine,')  = ',(zsfx)
-      IF( ABS( zhfx ) > zh_sill ) WRITE(numout,*) 'violation hfx    [GW]           (',cd_routine,')  = ',(zhfx)
-
-   END SUBROUTINE ice_cons_final
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option         Empty module            NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-#endif
-   !!======================================================================
-END MODULE icecons
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8426)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8426)
@@ -0,0 +1,211 @@
+MODULE icecor
+   !!======================================================================
+   !!                     ***  MODULE  icecor  ***
+   !!   LIM-3 : Update of sea-ice global variables at the end of the time step
+   !!======================================================================
+   !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code
+   !!            3.5  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!    ice_cor   : computes update of sea-ice global variables from trend terms
+   !!----------------------------------------------------------------------
+   USE dom_oce
+   USE phycst          ! physical constants
+   USE ice
+   USE ice1D           ! LIM thermodynamic sea-ice variables
+   USE iceitd
+   USE icevar
+   USE icectl          ! control prints
+   !
+   USE in_out_manager  ! I/O manager
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lbclnk          ! lateral boundary condition - MPP link
+   USE lib_mpp         ! MPP library
+   USE timing          ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_cor
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
+   !! $Id: icecor.F90 8378 2017-07-26 13:55:59Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_cor( kt , kn )
+      !!-------------------------------------------------------------------
+      !!               ***  ROUTINE ice_cor  ***
+      !!               
+      !! ** Purpose :  Computes update of sea-ice global variables at 
+      !!               the end of the dynamics.
+      !!                
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
+      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo
+      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      REAL(wp) ::   zsal
+      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      !!-------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('icecor')
+
+      IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)' icecor ' 
+         WRITE(numout,*)' ~~~~~~ '
+      ENDIF
+
+      ! conservation test
+      IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+      !----------------------------------------------------------------------
+      ! Constrain the thickness of the smallest category above himin
+      !----------------------------------------------------------------------
+      IF( kn == 2 ) THEN
+         
+         DO jj = 1, jpj 
+            DO ji = 1, jpi
+               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes
+               ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch
+               IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
+                  a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin
+               ENDIF
+            END DO
+         END DO
+
+      ENDIF
+         
+      !----------------------------------------------------
+      ! ice concentration should not exceed amax 
+      !-----------------------------------------------------
+      at_i(:,:) = 0._wp
+      DO jl = 1, jpl
+         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
+      END DO
+
+      DO jl  = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
+                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
+               ENDIF
+            END DO
+         END DO
+      END DO
+    
+      !---------------------
+      ! Ice salinity bounds
+      !---------------------
+      IF (  nn_icesal == 2  ) THEN 
+         DO jl = 1, jpl
+            DO jj = 1, jpj 
+               DO ji = 1, jpi
+                  zsal            = smv_i(ji,jj,jl)
+                  ! salinity stays in bounds
+                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
+                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) )
+                  ! associated salt flux
+                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      !----------------------------------------------------
+      ! Rebin categories with thickness out of bounds
+      !----------------------------------------------------
+      IF ( jpl > 1 ) CALL ice_itd_reb
+
+      !-----------------
+      ! zap small values
+      !-----------------
+      CALL ice_var_zapsmall
+
+      !----------------------------------------------
+      ! Ice drift. Corrections to avoid wrong values
+      !----------------------------------------------
+      IF( kn == 2 ) THEN
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
+                  IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
+                  IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
+                  IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
+                  IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
+               ENDIF
+            END DO
+         END DO
+         !lateral boundary conditions
+         CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )
+         !mask velocities
+         u_ice(:,:) = u_ice(:,:) * umask(:,:,1)
+         v_ice(:,:) = v_ice(:,:) * vmask(:,:,1)
+      ENDIF
+      
+      ! -------------------------------------------------
+      ! Diagnostics
+      ! -------------------------------------------------
+      IF( kn == 1 ) THEN
+         DO jl  = 1, jpl
+            afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
+         END DO
+         
+         DO jj = 1, jpj
+            DO ji = 1, jpi            
+               ! heat content variation (W.m-2)
+               diag_heat(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
+                  &                   SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
+                  &                 ) * r1_rdtice
+               ! salt, volume
+               diag_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
+               diag_vice(ji,jj) = SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
+               diag_vsnw(ji,jj) = SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
+            END DO
+         END DO
+         
+      ELSEIF( kn == 2 ) THEN
+         
+         DO jl  = 1, jpl
+            oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice          ! ice natural aging
+            afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
+         END DO
+         afx_tot = afx_thd + afx_dyn
+         
+         DO jj = 1, jpj
+            DO ji = 1, jpi            
+               ! heat content variation (W.m-2)
+               diag_heat(ji,jj) = diag_heat(ji,jj) -  &
+                  &               ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
+                  &                 SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
+                  &               ) * r1_rdtice   
+               ! salt, volume
+               diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
+               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
+               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
+            END DO
+         END DO
+         
+      ENDIF
+      
+      ! conservation test
+      IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+      ! control prints
+      IF( ln_ctl )       CALL ice_prt3D( 'icecor' )
+      IF( ln_limctl .AND. kn == 2 )  &
+         &               CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )
+   
+      IF( nn_timing == 1 )  CALL timing_stop('icecor')
+
+   END SUBROUTINE ice_cor
+
+#endif
+
+END MODULE icecor
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8426)
@@ -32,4 +32,6 @@
    PRIVATE
 
+   PUBLIC   ice_cons_hsm
+   PUBLIC   ice_cons_final
    PUBLIC   ice_ctl
    PUBLIC   ice_prt
@@ -46,4 +48,158 @@
 CONTAINS
 
+   SUBROUTINE ice_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
+      !!--------------------------------------------------------------------------------------------------------
+      !!                                        ***  ROUTINE ice_cons_hsm ***
+      !!
+      !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine
+      !!                     + test if ice concentration and volume are > 0
+      !!
+      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiachk=true
+      !!              It prints in ocean.output if there is a violation of conservation at each time-step
+      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to
+      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
+      !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
+      !!              and a heat content of 3e5 J/kg (=latent heat of fusion) 
+      !!--------------------------------------------------------------------------------------------------------
+      INTEGER         , INTENT(in)    :: icount        ! determine wether this is the beggining of the routine (0) or the end (1)
+      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine
+      REAL(wp)        , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft
+      REAL(wp)                        :: zvmin, zamin, zamax 
+      REAL(wp)                        :: zvtrp, zetrp
+      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill
+      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
+
+      IF( icount == 0 ) THEN
+
+         ! salt flux
+         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
+            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &
+            &                ) *  e1e2t(:,:) * tmask(:,:,1) * zconv )
+
+         ! water flux
+         zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
+            &                  wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
+            &                  wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        &
+            &                ) * e1e2t(:,:) * tmask(:,:,1) * zconv )
+
+         ! heat flux
+         zft_b  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
+            &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
+            &                ) *  e1e2t(:,:) * tmask(:,:,1) * zconv )
+
+         zvi_b  = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * tmask(:,:,1) * zconv )
+
+         zsmv_b = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * tmask(:,:,1) * zconv )
+
+         zei_b  = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
+            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    &
+                            ) * e1e2t * tmask(:,:,1) * zconv )
+
+      ELSEIF( icount == 1 ) THEN
+
+         ! salt flux
+         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
+            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
+            &              ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zfs_b
+
+         ! water flux
+         zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
+            &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
+            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        &
+            &              ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zfw_b
+
+         ! heat flux
+         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
+            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
+            &              ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zft_b
+ 
+         ! outputs
+         zvi  = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 )  &
+            &                    * e1e2t * tmask(:,:,1) * zconv ) - zvi_b ) * r1_rdtice - zfw ) * rday
+
+         zsmv = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 )  &
+            &                    * e1e2t * tmask(:,:,1) * zconv ) - zsmv_b ) * r1_rdtice + zfs ) * rday
+
+         zei  = ( glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
+            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    &
+            &                ) * e1e2t * tmask(:,:,1) * zconv ) - zei_b ) * r1_rdtice + zft
+
+         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative
+         zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t * tmask(:,:,1) * zconv ) * rday 
+         zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e1e2t * tmask(:,:,1) * zconv )
+
+         zvmin = glob_min( v_i )
+         zamax = glob_max( SUM( a_i, dim=3 ) )
+         zamin = glob_min( a_i )
+
+         ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
+         zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t * zconv ) ! in 1.e9 m2
+         zv_sill = zarea * 2.5e-5
+         zs_sill = zarea * 25.e-5
+         zh_sill = zarea * 10.e-5
+
+         IF(lwp) THEN
+            IF ( ABS( zvi  ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zvi
+            IF ( ABS( zsmv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv
+            IF ( ABS( zei  ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zei
+            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'iceadv' ) THEN
+                                         WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
+                                         WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp
+            ENDIF
+            IF (     zvmin   < -epsi10 ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
+            IF (     zamax   > MAX( rn_amax_n, rn_amax_s ) + epsi10 .AND. &
+               &                         cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' ) THEN
+                                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
+            IF (     zamax   > 1._wp   ) WRITE(numout,*) 'violation a_i>1               (',cd_routine,') = ',zamax
+            ENDIF
+            IF (      zamin  < -epsi10 ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+         ENDIF
+
+      ENDIF
+
+   END SUBROUTINE ice_cons_hsm
+
+
+   SUBROUTINE ice_cons_final( cd_routine )
+      !!---------------------------------------------------------------------------------------------------------
+      !!                                   ***  ROUTINE ice_cons_final ***
+      !!
+      !! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step
+      !!
+      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiachk=true
+      !!              It prints in ocean.output if there is a violation of conservation at each time-step
+      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to
+      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
+      !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
+      !!              and a heat content of 3e5 J/kg (=latent heat of fusion) 
+      !!--------------------------------------------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine
+      REAL(wp)                        :: zhfx, zsfx, zvfx
+      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill
+      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
+
+      ! heat flux
+      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es   &
+      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check)
+         &              ) * e1e2t * tmask(:,:,1) * zconv ) 
+      ! salt flux
+      zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t * tmask(:,:,1) * zconv ) * rday
+      ! water flux
+      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t * tmask(:,:,1) * zconv ) * rday
+
+      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
+      zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t * zconv ) ! in 1.e9 m2
+      zv_sill = zarea * 2.5e-5
+      zs_sill = zarea * 25.e-5
+      zh_sill = zarea * 10.e-5
+
+      IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx    [Mt/day]       (',cd_routine,')  = ',(zvfx)
+      IF( ABS( zsfx ) > zs_sill ) WRITE(numout,*) 'violation sfx    [psu*Mt/day]   (',cd_routine,')  = ',(zsfx)
+      IF( ABS( zhfx ) > zh_sill ) WRITE(numout,*) 'violation hfx    [GW]           (',cd_routine,')  = ',(zhfx)
+
+   END SUBROUTINE ice_cons_final
+
+   
    SUBROUTINE ice_ctl( kt )
       !!-----------------------------------------------------------------------
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90	(revision 8426)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90	(revision 8426)
@@ -0,0 +1,262 @@
+MODULE icedia
+   !!======================================================================
+   !!                       ***  MODULE limdia_hsb   ***
+   !!  LIM-3 sea ice model :   diagnostics of ice model 
+   !!======================================================================
+   !! History :  3.4  ! 2012-10  (C. Rousset)  original code
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                       LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_dia_hsb        : computation and output of the time evolution of keys variables
+   !!   lim_dia_hsb_init   : initialization and namelist read
+   !!----------------------------------------------------------------------
+   USE ice             ! LIM-3: sea-ice variable
+   USE dom_oce         ! ocean domain
+   USE sbc_oce, ONLY: sfx         ! surface boundary condition: ocean fields
+   USE daymod          ! model calendar
+   USE phycst          ! physical constant
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
+   USE timing          ! preformance summary
+   USE iom             ! I/O manager
+   USE lib_fortran     ! glob_sum
+   USE icerst          ! ice restart
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_dia        ! routine called by icestp.F90
+   PUBLIC   ice_dia_init   ! routine called in icestp.F90
+
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents
+   REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends
+   
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.4 , NEMO Consortium (2012)
+   !! $Id: icedia.F90 8413 2017-08-07 17:05:39Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE ice_dia( kt )
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_dia  ***
+      !!     
+      !! ** Purpose: Compute the ice global heat content, salt content and volume conservation
+      !!	
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt    ! number of iteration
+      !!
+      real(wp)   ::   zbg_ivol, zbg_svol, zbg_area, zbg_isal, zbg_item ,zbg_stem
+      REAL(wp)   ::   z_frc_voltop, z_frc_volbot, z_frc_sal, z_frc_temtop, z_frc_tembot  
+      REAL(wp)   ::   zdiff_vol, zdiff_sal, zdiff_tem  
+      !!---------------------------------------------------------------------------
+      IF( nn_timing == 1 )   CALL timing_start('ice_dia')
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'icedia'
+         WRITE(numout,*)'~~~~~~'
+      ENDIF
+
+      ! ----------------------- !
+      ! 1 -  Contents !
+      ! ----------------------- !
+      zbg_ivol = glob_sum( vt_i(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! ice volume (km3)
+      zbg_svol = glob_sum( vt_s(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! snow volume (km3)
+      zbg_area = glob_sum( at_i(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-6 )                  ! area (km2)
+      zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt content (pss*km3)
+      zbg_item = glob_sum( et_i * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J)
+      zbg_stem = glob_sum( et_s * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J)
+      
+      ! ---------------------------!
+      ! 2 - Trends due to forcing  !
+      ! ---------------------------!
+      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )  ! freshwater flux ice/snow-ocean 
+      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                     ! freshwater flux ice/snow-atm
+      z_frc_sal    = r1_rau0 * glob_sum( - sfx(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                                            ! salt fluxes ice/snow-ocean
+      z_frc_tembot =           glob_sum( hfx_out(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ocean (and below ice)
+      z_frc_temtop =           glob_sum( hfx_in (:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ice-coean
+      !
+      frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3
+      frc_volbot  = frc_volbot  + z_frc_volbot  * rdt_ice ! km3
+      frc_sal     = frc_sal     + z_frc_sal     * rdt_ice ! km3*pss
+      frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J
+      frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J
+            
+      ! ----------------------- !
+      ! 3 -  Content variations !
+      ! ----------------------- !
+      zdiff_vol = r1_rau0 * glob_sum( ( rhoic * vt_i(:,:) + rhosn * vt_s(:,:) - vol_loc_ini(:,:)  &  ! freshwater trend (km3) 
+         &                            ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) 
+      zdiff_sal = r1_rau0 * glob_sum( ( rhoic * SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:)     &  ! salt content trend (km3*pss)
+         &                            ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )
+      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:)                  &  ! heat content trend (1.e20 J)
+      !  &                            + SUM( qevap_ice * a_i_b, dim=3 ) &     !! clem: I think this line should be commented (but needs a check)
+         &                            ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )
+
+      ! ----------------------- !
+      ! 4 -  Drifts             !
+      ! ----------------------- !
+      zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot )
+      zdiff_sal = zdiff_sal - frc_sal
+      zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop )
+
+      ! ----------------------- !
+      ! 5 - Diagnostics writing !
+      ! ----------------------- !
+      !
+      IF( iom_use('ibgvolume') )  CALL iom_put( 'ibgvolume' , zdiff_vol        )   ! ice/snow volume  drift            (km3 equivalent ocean water)         
+      IF( iom_use('ibgsaltco') )  CALL iom_put( 'ibgsaltco' , zdiff_sal        )   ! ice salt content drift            (psu*km3 equivalent ocean water)
+      IF( iom_use('ibgheatco') )  CALL iom_put( 'ibgheatco' , zdiff_tem        )   ! ice/snow heat content drift       (1.e20 J)
+      IF( iom_use('ibgheatfx') )  CALL iom_put( 'ibgheatfx' , zdiff_tem /      &   ! ice/snow heat flux drift          (W/m2)
+         &                                                    glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
+
+      IF( iom_use('ibgfrcvoltop') )  CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water) 
+      IF( iom_use('ibgfrcvolbot') )  CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water) 
+      IF( iom_use('ibgfrcsal') )     CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)   
+      IF( iom_use('ibgfrctemtop') )  CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)   
+      IF( iom_use('ibgfrctembot') )  CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)   
+      IF( iom_use('ibgfrchfxtop') )  CALL iom_put( 'ibgfrchfxtop' , frc_temtop / & ! heat on top of ice/snw/ocean      (W/m2) 
+         &                                                    glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
+      IF( iom_use('ibgfrchfxbot') )  CALL iom_put( 'ibgfrchfxbot' , frc_tembot / & ! heat on top of ocean(below ice)   (W/m2) 
+         &                                                    glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
+
+      IF( iom_use('ibgvol_tot' ) )  CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                        (km3)
+      IF( iom_use('sbgvol_tot' ) )  CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                       (km3)
+      IF( iom_use('ibgarea_tot') )  CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                          (km2)
+      IF( iom_use('ibgsalt_tot') )  CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content              (pss*km3)
+      IF( iom_use('ibgheat_tot') )  CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                  (1.e20 J)
+      IF( iom_use('sbgheat_tot') )  CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                 (1.e20 J)
+      !
+      IF( lrst_ice )   CALL ice_dia_rst( 'WRITE', kt_ice )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('ice_dia')
+      !
+   END SUBROUTINE ice_dia
+
+
+   SUBROUTINE ice_dia_init
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_dia_init  ***
+      !!     
+      !! ** Purpose: Initialization for the heat salt volume budgets
+      !!	
+      !! ** Method : Compute initial heat content, salt content and volume
+      !!
+      !! ** Action : - Compute initial heat content, salt content and volume
+      !!             - Initialize forcing trends
+      !!             - Compute coefficients for conversion
+      !!---------------------------------------------------------------------------
+      INTEGER            ::   ierror   ! local integer
+      !!
+      !!NAMELIST/namicehsb/ blabla
+      !!----------------------------------------------------------------------
+      !
+      !!REWIND ( numnam_ice )              ! Read Namelist namicehsb 
+      !!READ   ( numnam_ice, namicehsb )
+      !
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_dia_init : check the heat and salt budgets'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+      !      
+      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 )  THEN
+         CALL ctl_stop( 'ice_dia: unable to allocate vol_loc_ini' )
+         RETURN
+      ENDIF
+
+      CALL ice_dia_rst( 'READ' )  !* read or initialize all required files
+      !
+   END SUBROUTINE ice_dia_init
+
+   SUBROUTINE ice_dia_rst( cdrw, kt )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE limdia_rst  ***
+     !!                     
+     !! ** Purpose :   Read or write DIA file in restart file
+     !!
+     !! ** Method  :   use of IOM library
+     !!----------------------------------------------------------------------
+     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+     INTEGER         , INTENT(in), OPTIONAL ::   kt     ! ice time-step
+     REAL(wp)                     ::   ziter
+     INTEGER                      ::   iter
+     !
+     !!----------------------------------------------------------------------
+     !
+     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+        IF( ln_rstart ) THEN                   !* Read the restart file
+           !
+           CALL iom_get( numrir, 'kt_ice' , ziter )
+           IF(lwp) WRITE(numout,*)
+           IF(lwp) WRITE(numout,*) ' ice_dia_rst read at time step = ', ziter
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
+           CALL iom_get( numrir, 'frc_voltop' , frc_voltop  )
+           CALL iom_get( numrir, 'frc_volbot' , frc_volbot  )
+           CALL iom_get( numrir, 'frc_temtop' , frc_temtop  )
+           CALL iom_get( numrir, 'frc_tembot' , frc_tembot  )
+           CALL iom_get( numrir, 'frc_sal'    , frc_sal     )
+           CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini )
+           CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini )
+           CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini )
+        ELSE
+           IF(lwp) WRITE(numout,*)
+           IF(lwp) WRITE(numout,*) ' ice_dia at initial state '
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
+           ! set trends to 0
+           frc_voltop  = 0._wp                                          
+           frc_volbot  = 0._wp                                          
+           frc_temtop  = 0._wp                                                 
+           frc_tembot  = 0._wp                                                 
+           frc_sal     = 0._wp                                                 
+           ! record initial ice volume, salt and temp
+           vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2)
+           tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J)
+           sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2)
+           
+       ENDIF
+
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+        !                                   ! -------------------
+        iter = kt + nn_fsbc - 1   ! ice restarts are written at kt == nitrst - nn_fsbc + 1
+
+        IF( iter == nitrst ) THEN
+           IF(lwp) WRITE(numout,*)
+           IF(lwp) WRITE(numout,*) ' ice_dia_rst write at time step = ', kt
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
+        ENDIF
+
+        ! Write in numriw (if iter == nitrst)
+        ! ------------------ 
+        CALL iom_rstput( iter, nitrst, numriw, 'frc_voltop' , frc_voltop  )
+        CALL iom_rstput( iter, nitrst, numriw, 'frc_volbot' , frc_volbot  )
+        CALL iom_rstput( iter, nitrst, numriw, 'frc_temtop' , frc_temtop  )
+        CALL iom_rstput( iter, nitrst, numriw, 'frc_tembot' , frc_tembot  )
+        CALL iom_rstput( iter, nitrst, numriw, 'frc_sal'    , frc_sal     )
+        CALL iom_rstput( iter, nitrst, numriw, 'vol_loc_ini', vol_loc_ini )
+        CALL iom_rstput( iter, nitrst, numriw, 'tem_loc_ini', tem_loc_ini )
+        CALL iom_rstput( iter, nitrst, numriw, 'sal_loc_ini', sal_loc_ini )
+        !
+     ENDIF
+     !
+   END SUBROUTINE ice_dia_rst
+ 
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :         Empty module          NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE ice_dia          ! Empty routine
+   END SUBROUTINE ice_dia
+#endif
+   !!======================================================================
+END MODULE icedia
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icediahsb.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icediahsb.F90	(revision 8424)
+++ 	(revision )
@@ -1,256 +1,0 @@
-MODULE icediahsb
-   !!======================================================================
-   !!                       ***  MODULE limdia_hsb   ***
-   !!  LIM-3 sea ice model :   diagnostics of ice model 
-   !!======================================================================
-   !! History :  3.4  ! 2012-10  (C. Rousset)  original code
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                       LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_dia_hsb        : computation and output of the time evolution of keys variables
-   !!   lim_dia_hsb_init   : initialization and namelist read
-   !!----------------------------------------------------------------------
-   USE ice             ! LIM-3: sea-ice variable
-   USE dom_oce         ! ocean domain
-   USE sbc_oce, ONLY: sfx         ! surface boundary condition: ocean fields
-   USE daymod          ! model calendar
-   USE phycst          ! physical constant
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! MPP library
-   USE timing          ! preformance summary
-   USE iom             ! I/O manager
-   USE lib_fortran     ! glob_sum
-   USE icerst          ! ice restart
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_diahsb        ! routine called by icestp.F90
-   PUBLIC   ice_diahsb_init   ! routine called in icestp.F90
-
-   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents
-   REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends
-   
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.4 , NEMO Consortium (2012)
-   !! $Id: icediahsb.F90 8413 2017-08-07 17:05:39Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE ice_diahsb( kt )
-      !!---------------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_diahsb  ***
-      !!     
-      !! ** Purpose: Compute the ice global heat content, salt content and volume conservation
-      !!	
-      !!---------------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt    ! number of iteration
-      !!
-      real(wp)   ::   zbg_ivol, zbg_svol, zbg_area, zbg_isal, zbg_item ,zbg_stem
-      REAL(wp)   ::   z_frc_voltop, z_frc_volbot, z_frc_sal, z_frc_temtop, z_frc_tembot  
-      REAL(wp)   ::   zdiff_vol, zdiff_sal, zdiff_tem  
-      !!---------------------------------------------------------------------------
-      IF( nn_timing == 1 )   CALL timing_start('ice_diahsb')
-
-      ! ----------------------- !
-      ! 1 -  Contents !
-      ! ----------------------- !
-      zbg_ivol = glob_sum( vt_i(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! ice volume (km3)
-      zbg_svol = glob_sum( vt_s(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! snow volume (km3)
-      zbg_area = glob_sum( at_i(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-6 )                  ! area (km2)
-      zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt content (pss*km3)
-      zbg_item = glob_sum( et_i * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J)
-      zbg_stem = glob_sum( et_s * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J)
-      
-      ! ---------------------------!
-      ! 2 - Trends due to forcing  !
-      ! ---------------------------!
-      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )  ! freshwater flux ice/snow-ocean 
-      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                     ! freshwater flux ice/snow-atm
-      z_frc_sal    = r1_rau0 * glob_sum( - sfx(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                                            ! salt fluxes ice/snow-ocean
-      z_frc_tembot =           glob_sum( hfx_out(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ocean (and below ice)
-      z_frc_temtop =           glob_sum( hfx_in (:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ice-coean
-      !
-      frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3
-      frc_volbot  = frc_volbot  + z_frc_volbot  * rdt_ice ! km3
-      frc_sal     = frc_sal     + z_frc_sal     * rdt_ice ! km3*pss
-      frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J
-      frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J
-            
-      ! ----------------------- !
-      ! 3 -  Content variations !
-      ! ----------------------- !
-      zdiff_vol = r1_rau0 * glob_sum( ( rhoic * vt_i(:,:) + rhosn * vt_s(:,:) - vol_loc_ini(:,:)  &  ! freshwater trend (km3) 
-         &                            ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) 
-      zdiff_sal = r1_rau0 * glob_sum( ( rhoic * SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:)     &  ! salt content trend (km3*pss)
-         &                            ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )
-      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:)                  &  ! heat content trend (1.e20 J)
-      !  &                            + SUM( qevap_ice * a_i_b, dim=3 ) &     !! clem: I think this line should be commented (but needs a check)
-         &                            ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )
-
-      ! ----------------------- !
-      ! 4 -  Drifts             !
-      ! ----------------------- !
-      zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot )
-      zdiff_sal = zdiff_sal - frc_sal
-      zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop )
-
-      ! ----------------------- !
-      ! 5 - Diagnostics writing !
-      ! ----------------------- !
-      !
-      IF( iom_use('ibgvolume') )  CALL iom_put( 'ibgvolume' , zdiff_vol        )   ! ice/snow volume  drift            (km3 equivalent ocean water)         
-      IF( iom_use('ibgsaltco') )  CALL iom_put( 'ibgsaltco' , zdiff_sal        )   ! ice salt content drift            (psu*km3 equivalent ocean water)
-      IF( iom_use('ibgheatco') )  CALL iom_put( 'ibgheatco' , zdiff_tem        )   ! ice/snow heat content drift       (1.e20 J)
-      IF( iom_use('ibgheatfx') )  CALL iom_put( 'ibgheatfx' , zdiff_tem /      &   ! ice/snow heat flux drift          (W/m2)
-         &                                                    glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
-
-      IF( iom_use('ibgfrcvoltop') )  CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water) 
-      IF( iom_use('ibgfrcvolbot') )  CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water) 
-      IF( iom_use('ibgfrcsal') )     CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)   
-      IF( iom_use('ibgfrctemtop') )  CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)   
-      IF( iom_use('ibgfrctembot') )  CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)   
-      IF( iom_use('ibgfrchfxtop') )  CALL iom_put( 'ibgfrchfxtop' , frc_temtop / & ! heat on top of ice/snw/ocean      (W/m2) 
-         &                                                    glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
-      IF( iom_use('ibgfrchfxbot') )  CALL iom_put( 'ibgfrchfxbot' , frc_tembot / & ! heat on top of ocean(below ice)   (W/m2) 
-         &                                                    glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
-
-      IF( iom_use('ibgvol_tot' ) )  CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                        (km3)
-      IF( iom_use('sbgvol_tot' ) )  CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                       (km3)
-      IF( iom_use('ibgarea_tot') )  CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                          (km2)
-      IF( iom_use('ibgsalt_tot') )  CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content              (pss*km3)
-      IF( iom_use('ibgheat_tot') )  CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                  (1.e20 J)
-      IF( iom_use('sbgheat_tot') )  CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                 (1.e20 J)
-      !
-      IF( lrst_ice )   CALL ice_diahsb_rst( 'WRITE', kt_ice )
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('ice_diahsb')
-      !
-   END SUBROUTINE ice_diahsb
-
-
-   SUBROUTINE ice_diahsb_init
-      !!---------------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_diahsb_init  ***
-      !!     
-      !! ** Purpose: Initialization for the heat salt volume budgets
-      !!	
-      !! ** Method : Compute initial heat content, salt content and volume
-      !!
-      !! ** Action : - Compute initial heat content, salt content and volume
-      !!             - Initialize forcing trends
-      !!             - Compute coefficients for conversion
-      !!---------------------------------------------------------------------------
-      INTEGER            ::   ierror   ! local integer
-      !!
-      !!NAMELIST/namicehsb/ blabla
-      !!----------------------------------------------------------------------
-      !
-      !!REWIND ( numnam_ice )              ! Read Namelist namicehsb 
-      !!READ   ( numnam_ice, namicehsb )
-      !
-      IF(lwp) THEN                   ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'ice_diahsb_init : check the heat and salt budgets'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-      ENDIF
-      !      
-      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror )
-      IF( ierror > 0 )  THEN
-         CALL ctl_stop( 'ice_diahsb: unable to allocate vol_loc_ini' )
-         RETURN
-      ENDIF
-
-      CALL ice_diahsb_rst( 'READ' )  !* read or initialize all required files
-      !
-   END SUBROUTINE ice_diahsb_init
-
-   SUBROUTINE ice_diahsb_rst( cdrw, kt )
-     !!---------------------------------------------------------------------
-     !!                   ***  ROUTINE limdia_rst  ***
-     !!                     
-     !! ** Purpose :   Read or write DIA file in restart file
-     !!
-     !! ** Method  :   use of IOM library
-     !!----------------------------------------------------------------------
-     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
-     INTEGER         , INTENT(in), OPTIONAL ::   kt     ! ice time-step
-     REAL(wp)                     ::   ziter
-     INTEGER                      ::   iter
-     !
-     !!----------------------------------------------------------------------
-     !
-     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
-        IF( ln_rstart ) THEN                   !* Read the restart file
-           !
-           CALL iom_get( numrir, 'kt_ice' , ziter )
-           IF(lwp) WRITE(numout,*)
-           IF(lwp) WRITE(numout,*) ' ice_diahsb_rst read at time step = ', ziter
-           IF(lwp) WRITE(numout,*) '~~~~~~~'
-           CALL iom_get( numrir, 'frc_voltop' , frc_voltop  )
-           CALL iom_get( numrir, 'frc_volbot' , frc_volbot  )
-           CALL iom_get( numrir, 'frc_temtop' , frc_temtop  )
-           CALL iom_get( numrir, 'frc_tembot' , frc_tembot  )
-           CALL iom_get( numrir, 'frc_sal'    , frc_sal     )
-           CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini )
-           CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini )
-           CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini )
-        ELSE
-           IF(lwp) WRITE(numout,*)
-           IF(lwp) WRITE(numout,*) ' ice_diahsb at initial state '
-           IF(lwp) WRITE(numout,*) '~~~~~~~'
-           ! set trends to 0
-           frc_voltop  = 0._wp                                          
-           frc_volbot  = 0._wp                                          
-           frc_temtop  = 0._wp                                                 
-           frc_tembot  = 0._wp                                                 
-           frc_sal     = 0._wp                                                 
-           ! record initial ice volume, salt and temp
-           vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2)
-           tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J)
-           sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2)
-           
-       ENDIF
-
-     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
-        !                                   ! -------------------
-        iter = kt + nn_fsbc - 1   ! ice restarts are written at kt == nitrst - nn_fsbc + 1
-
-        IF( iter == nitrst ) THEN
-           IF(lwp) WRITE(numout,*)
-           IF(lwp) WRITE(numout,*) ' ice_diahsb_rst write at time step = ', kt
-           IF(lwp) WRITE(numout,*) '~~~~~~~'
-        ENDIF
-
-        ! Write in numriw (if iter == nitrst)
-        ! ------------------ 
-        CALL iom_rstput( iter, nitrst, numriw, 'frc_voltop' , frc_voltop  )
-        CALL iom_rstput( iter, nitrst, numriw, 'frc_volbot' , frc_volbot  )
-        CALL iom_rstput( iter, nitrst, numriw, 'frc_temtop' , frc_temtop  )
-        CALL iom_rstput( iter, nitrst, numriw, 'frc_tembot' , frc_tembot  )
-        CALL iom_rstput( iter, nitrst, numriw, 'frc_sal'    , frc_sal     )
-        CALL iom_rstput( iter, nitrst, numriw, 'vol_loc_ini', vol_loc_ini )
-        CALL iom_rstput( iter, nitrst, numriw, 'tem_loc_ini', tem_loc_ini )
-        CALL iom_rstput( iter, nitrst, numriw, 'sal_loc_ini', sal_loc_ini )
-        !
-     ENDIF
-     !
-   END SUBROUTINE ice_diahsb_rst
- 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :         Empty module          NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE ice_diahsb          ! Empty routine
-   END SUBROUTINE ice_diahsb
-#endif
-   !!======================================================================
-END MODULE icediahsb
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr1.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr1.F90	(revision 8424)
+++ 	(revision )
@@ -1,153 +1,0 @@
-MODULE iceerr1
-   !!======================================================================
-   !!                     ***  MODULE  iceerr1  ***
-   !!   LIM-3 : Update of sea-ice global variables at the end of the time step
-   !!======================================================================
-   !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code
-   !!            3.5  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!    ice_err1   : computes update of sea-ice global variables from trend terms
-   !!----------------------------------------------------------------------
-   USE dom_oce
-   USE phycst          ! physical constants
-   USE ice
-   USE ice1D           ! LIM thermodynamic sea-ice variables
-   USE iceitd
-   USE icevar
-   USE icecons         ! conservation tests
-   USE icectl          ! control prints
-   !
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE in_out_manager  ! I/O manager
-   USE timing          ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_err1
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id: iceerr1.F90 8378 2017-07-26 13:55:59Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_err1( kt )
-      !!-------------------------------------------------------------------
-      !!               ***  ROUTINE ice_err1  ***
-      !!               
-      !! ** Purpose :  Computes update of sea-ice global variables at 
-      !!               the end of the dynamics.
-      !!                
-      !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! number of iteration
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zsal
-      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
-      !!-------------------------------------------------------------------
-      IF( nn_timing == 1 )  CALL timing_start('iceerr1')
-
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)'' 
-         WRITE(numout,*)' ice_err1 ' 
-         WRITE(numout,*)' ~~~~~~~~~~~ '
-      ENDIF
-
-      ! conservation test
-      IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'iceerr1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      !----------------------------------------------------
-      ! ice concentration should not exceed amax 
-      !-----------------------------------------------------
-      at_i(:,:) = 0._wp
-      DO jl = 1, jpl
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-      END DO
-
-      DO jl  = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
-                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
-               ENDIF
-            END DO
-         END DO
-      END DO
-    
-      !---------------------
-      ! Ice salinity bounds
-      !---------------------
-      IF (  nn_icesal == 2  ) THEN 
-         DO jl = 1, jpl
-            DO jj = 1, jpj 
-               DO ji = 1, jpi
-                  zsal            = smv_i(ji,jj,jl)
-                  ! salinity stays in bounds
-                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
-                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) )
-                  ! associated salt flux
-                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      !----------------------------------------------------
-      ! Rebin categories with thickness out of bounds
-      !----------------------------------------------------
-      IF ( jpl > 1 ) CALL ice_itd_reb
-
-      !-----------------
-      ! zap small values
-      !-----------------
-      CALL ice_var_zapsmall
-
-      ! -------------------------------------------------
-      ! Diagnostics
-      ! -------------------------------------------------
-      DO jl  = 1, jpl
-         afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
-      END DO
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi            
-            ! heat content variation (W.m-2)
-            diag_heat(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
-               &                   SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
-               &                 ) * r1_rdtice
-            ! salt, volume
-            diag_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
-            diag_vice(ji,jj) = SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
-            diag_vsnw(ji,jj) = SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
-         END DO
-      END DO
-
-      ! conservation test
-      IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'iceerr1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      ! control prints
-      IF( ln_ctl )       CALL ice_prt3D( 'iceerr1' )
-   
-      IF( nn_timing == 1 )  CALL timing_stop('iceerr1')
-
-   END SUBROUTINE ice_err1
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option         Empty Module               No sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE ice_err1     ! Empty routine
-   END SUBROUTINE ice_err1
-
-#endif
-
-END MODULE iceerr1
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr2.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr2.F90	(revision 8424)
+++ 	(revision )
@@ -1,192 +1,0 @@
-MODULE iceerr2
-   !!======================================================================
-   !!                     ***  MODULE  iceerr2  ***
-   !!   LIM-3 : Update of sea-ice global variables at the end of the time step
-   !!======================================================================
-   !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code
-   !!            3.5  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!    ice_err2   : computes update of sea-ice global variables from trend terms
-   !!----------------------------------------------------------------------
-   USE dom_oce
-   USE phycst          ! physical constants
-   USE ice
-   USE ice1D           ! LIM thermodynamic sea-ice variables
-   USE iceitd
-   USE icevar
-   USE icecons         ! conservation tests
-   USE icectl
-   !
-   USE lbclnk          ! lateral boundary condition - MPP exchanges
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE in_out_manager
-   USE timing          ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_err2   ! routine called by ice_step
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id: iceerr2.F90 8378 2017-07-26 13:55:59Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_err2( kt )
-      !!-------------------------------------------------------------------
-      !!               ***  ROUTINE ice_err2  ***
-      !!               
-      !! ** Purpose :  Computes update of sea-ice global variables at 
-      !!               the end of the time step.
-      !!
-      !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! number of iteration
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zsal
-      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
-      !!-------------------------------------------------------------------
-      IF( nn_timing == 1 )  CALL timing_start('iceerr2')
-
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)''
-         WRITE(numout,*)' ice_err2 '
-         WRITE(numout,*)' ~~~~~~~~~~~ '
-      ENDIF
-
-      ! conservation test
-      IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'iceerr2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      !----------------------------------------------------------------------
-      ! Constrain the thickness of the smallest category above himin
-      !----------------------------------------------------------------------
-      DO jj = 1, jpj 
-         DO ji = 1, jpi
-            rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes
-            ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch
-            IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
-               a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin
-            ENDIF
-         END DO
-      END DO
-      
-      !-----------------------------------------------------
-      ! ice concentration should not exceed amax 
-      !-----------------------------------------------------
-      at_i(:,:) = 0._wp
-      DO jl = 1, jpl
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-      END DO
-
-      DO jl  = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
-                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
-               ENDIF
-            END DO
-         END DO
-      END DO
-
-      !---------------------
-      ! Ice salinity
-      !---------------------
-      IF (  nn_icesal == 2  ) THEN 
-         DO jl = 1, jpl
-            DO jj = 1, jpj 
-               DO ji = 1, jpi
-                  zsal            = smv_i(ji,jj,jl)
-                  ! salinity stays in bounds
-                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
-                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) )
-                  ! associated salt flux
-                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      !----------------------------------------------------
-      ! Rebin categories with thickness out of bounds
-      !----------------------------------------------------
-      IF ( jpl > 1 ) CALL ice_itd_reb
-
-      !-----------------
-      ! zap small values
-      !-----------------
-      CALL ice_var_zapsmall
-
-      !------------------------------------------------------------------------------
-      ! Corrections to avoid wrong values                                        |
-      !------------------------------------------------------------------------------
-      ! Ice drift
-      !------------
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-            IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
-               IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
-               IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
-               IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
-               IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
-            ENDIF
-         END DO
-      END DO
-      !lateral boundary conditions
-      CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. )
-      !mask velocities
-      u_ice(:,:) = u_ice(:,:) * umask(:,:,1)
-      v_ice(:,:) = v_ice(:,:) * vmask(:,:,1)
- 
-      ! -------------------------------------------------
-      ! Diagnostics
-      ! -------------------------------------------------
-      DO jl  = 1, jpl
-         oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice          ! ice natural aging
-         afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
-      END DO
-      afx_tot = afx_thd + afx_dyn
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi            
-            ! heat content variation (W.m-2)
-            diag_heat(ji,jj) = diag_heat(ji,jj) -  &
-               &               ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
-               &                 SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
-               &               ) * r1_rdtice   
-            ! salt, volume
-            diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
-            diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
-            diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
-         END DO
-      END DO
-
-      ! conservation test
-      IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'iceerr2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      ! control prints
-      IF( ln_limctl )    CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )
-      IF( ln_ctl )       CALL ice_prt3D( 'iceerr2' )
-   
-      IF( nn_timing == 1 )  CALL timing_stop('iceerr2')
-
-   END SUBROUTINE ice_err2
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option         Empty Module               No sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE ice_err2     ! Empty routine
-   END SUBROUTINE ice_err2
-
-#endif
-
-END MODULE iceerr2
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8426)
@@ -19,5 +19,5 @@
    USE sbcblk          ! Surface boundary condition: bulk
    USE sbccpl          ! Surface boundary condition: coupled interface
-   USE icealbedo       ! ice albedo
+   USE icealb          ! ice albedo
    !
    USE iom             ! I/O manager library
@@ -64,4 +64,10 @@
       IF( nn_timing == 1 )   CALL timing_start('ice_forcing_tau')
 
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'ice_forcing_tau'
+         WRITE(numout,*)'~~~~~~~~~~~~~~~'
+      ENDIF
+
       SELECT CASE( ksbc )
          CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
@@ -78,6 +84,5 @@
             END DO
          END DO
-         CALL lbc_lnk( utau_ice, 'U', -1. )
-         CALL lbc_lnk( vtau_ice, 'V', -1. )
+         CALL lbc_lnk_multi( utau_ice, 'U', -1., vtau_ice, 'V', -1. )
       ENDIF
 
@@ -119,6 +124,12 @@
       IF( nn_timing == 1 )   CALL timing_start('ice_forcing_flx')
 
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'ice_forcing_flx'
+         WRITE(numout,*)'~~~~~~~~~~~~~~~'
+      ENDIF
+
       ! --- cloud-sky and overcast-sky ice albedos --- !
-      CALL ice_albedo( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os )
+      CALL ice_alb( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os )
 
       ! albedo depends on cloud fraction because of non-linear spectral effects
@@ -140,9 +151,9 @@
                                    CALL blk_ice_flx( t_su, alb_ice )
             IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
-            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            IF( nn_limflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
 
          CASE ( jp_purecpl )                                     ! coupled formulation
                                    CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
-            IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            IF( nn_limflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
       END SELECT
 
@@ -162,7 +173,7 @@
 
 
-   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_lim_flx  ***
+   SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_flx_dist  ***
       !!
       !! ** Purpose :   update the ice surface boundary condition by averaging and / or
@@ -195,5 +206,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('ice_lim_flx')
+      IF( nn_timing == 1 )  CALL timing_start('ice_flx_dist')
       !
       SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==!
@@ -231,7 +242,7 @@
       END SELECT
       !
-      IF( nn_timing == 1 )  CALL timing_stop('ice_lim_flx')
-      !
-   END SUBROUTINE ice_lim_flx
+      IF( nn_timing == 1 )  CALL timing_stop('ice_flx_dist')
+      !
+   END SUBROUTINE ice_flx_dist
 
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90	(revision 8426)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90	(revision 8426)
@@ -0,0 +1,599 @@
+MODULE iceist
+   !!======================================================================
+   !!                     ***  MODULE  iceist  ***
+   !!              Initialisation of diagnostics ice variables
+   !!======================================================================
+   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
+   !!            3.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!             -   ! 2014    (C. Rousset) add N/S initializations
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                    LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_ist      :  Initialisation of diagnostics ice variables
+   !!   ice_ist_init :  initialization of ice state and namelist read
+   !!----------------------------------------------------------------------
+   USE phycst         ! physical constant
+   USE oce            ! dynamics and tracers variables
+   USE dom_oce        ! ocean domain
+   USE sbc_oce , ONLY : sst_m, sss_m 
+   USE sbc_ice , ONLY : tn_ice
+   USE eosbn2         ! equation of state
+   USE ice            ! sea-ice variables
+   USE par_oce        ! ocean parameters
+   USE icevar         ! ice_var_salprof
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE fldread        ! read input fields
+   USE iom
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_ist      ! routine called by lim_init.F90
+
+   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
+   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness (m)    at T-point
+   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snow thicknes (m)    at T-point
+   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point
+   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point
+   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point
+   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
+   !!----------------------------------------------------------------------
+   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
+   !! $Id: iceist.F90 8378 2017-07-26 13:55:59Z clem $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_ist
+      !!-------------------------------------------------------------------
+      !!                    ***  ROUTINE ice_ist  ***
+      !!
+      !! ** Purpose :   defined the sea-ice initial state
+      !!
+      !! ** Method  :   This routine will put some ice where ocean
+      !!                is at the freezing point, then fill in ice 
+      !!                state variables using prescribed initial 
+      !!                values in the namelist            
+      !!
+      !! ** Steps   :   1) Read namelist
+      !!                2) Basal temperature; ice and hemisphere masks
+      !!                3) Fill in the ice thickness distribution using gaussian
+      !!                4) Fill in space-dependent arrays for state variables
+      !!                5) Diagnostic arrays
+      !!                6) Lateral boundary conditions
+      !!
+      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
+      !!              where there is no ice (clem: I do not know why, is it mandatory?) 
+      !!
+      !! History :
+      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
+      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
+      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
+      !!--------------------------------------------------------------------
+      !! * Local variables
+      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
+      REAL(wp)   :: ztmelts, zdh
+      INTEGER    :: i_hemis, i_fill, jl0
+      REAL(wp)   :: zarg, zV, zconv, zdv
+      REAL(wp), DIMENSION(jpi,jpj)     :: zswitch    ! ice indicator
+      REAL(wp), DIMENSION(jpi,jpj)     :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
+      REAL(wp), DIMENSION(jpi,jpj)     :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini, za_i_ini                         !data by cattegories to fill
+      INTEGER , DIMENSION(4)           :: itest
+      !--------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'ice_ist : sea-ice initialization '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~ '
+
+      !--------------------------------------------------------------------
+      ! 1) Read namelist
+      !--------------------------------------------------------------------
+      !
+      CALL ice_ist_init
+
+      ! init surface temperature
+      DO jl = 1, jpl
+         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
+         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
+      END DO
+
+      ! init basal temperature (considered at freezing point)
+      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
+      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
+
+
+      !--------------------------------------------------------------------
+      ! 2) Initialization of sea ice state variables
+      !--------------------------------------------------------------------
+      IF( ln_limini ) THEN
+         !
+         IF( ln_limini_file )THEN
+         !
+            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1)
+            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1)
+            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1)
+            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1)
+            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1)
+            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1)
+            !
+            WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 
+            ELSEWHERE                       ; zswitch(:,:) = 0._wp
+            END WHERE
+            !
+         ELSE ! ln_limini_file = F
+
+            !--------------------------------------------------------------------
+            ! 3) Basal temperature, ice mask
+            !--------------------------------------------------------------------
+            ! no ice if sst <= t-freez + ttest
+            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 
+            ELSEWHERE                                                                  ; zswitch(:,:) = tmask(:,:,1)
+            END WHERE
+
+            !-----------------------------
+            ! 3.1) Hemisphere-dependent arrays
+            !-----------------------------
+            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( ff_t(ji,jj) >= 0._wp ) THEN
+                     zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj)
+                     zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj)
+                     zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj)
+                     zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
+                     zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj)
+                     ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
+                  ELSE
+                     zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj)
+                     zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj)
+                     zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj)
+                     zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
+                     zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj)
+                     ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
+                  ENDIF
+               END DO
+            END DO
+            !
+         ENDIF ! ln_limini_file
+         
+         zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)   ! ice volume
+         !---------------------------------------------------------------------
+         ! 3.2) Distribute ice concentration and thickness into the categories
+         !---------------------------------------------------------------------
+         ! a gaussian distribution for ice concentration is used
+         ! then we check whether the distribution fullfills
+         ! volume and area conservation, positivity and ice categories bounds
+         zh_i_ini(:,:,:) = 0._wp 
+         za_i_ini(:,:,:) = 0._wp
+         !
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
+
+                  !--- jl0: most likely index where cc will be maximum
+                  jl0 = jpl
+                  DO jl = 1, jpl
+                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
+                        jl0 = jl
+                        CYCLE
+                     ENDIF
+                  END DO
+                  !
+                  ! initialisation of tests
+                  itest(:)  = 0
+                  
+                  i_fill = jpl + 1                                             !====================================
+                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
+                     ! iteration                                               !====================================
+                     i_fill = i_fill - 1
+
+                     ! initialisation of ice variables for each try
+                     zh_i_ini(ji,jj,:) = 0._wp 
+                     za_i_ini(ji,jj,:) = 0._wp
+                     itest(:) = 0
+                     !
+                     ! *** case very thin ice: fill only category 1
+                     IF ( i_fill == 1 ) THEN
+                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
+                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
+
+                     ! *** case ice is thicker: fill categories >1
+                     ELSE
+
+                        ! Fill ice thicknesses in the (i_fill-1) cat by hmean 
+                        DO jl = 1, i_fill-1
+                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
+                        END DO
+                        !
+                        !--- Concentrations
+                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
+                        DO jl = 1, i_fill - 1
+                           IF( jl /= jl0 )THEN
+                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
+                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
+                           ENDIF
+                        END DO
+                        !
+                        ! Concentration in the last (i_fill) category
+                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
+
+                        ! Ice thickness in the last (i_fill) category
+                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
+                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
+
+                        ! clem: correction if concentration of upper cat is greater than lower cat
+                        !       (it should be a gaussian around jl0 but sometimes it is not)
+                        IF ( jl0 /= jpl ) THEN
+                           DO jl = jpl, jl0+1, -1
+                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
+                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
+                                 zh_i_ini(ji,jj,jl    ) = 0._wp
+                                 za_i_ini(ji,jj,jl    ) = 0._wp
+                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
+                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
+                              END IF
+                           ENDDO
+                        ENDIF
+                        !
+                     ENDIF ! case ice is thick or thin
+
+                     !---------------------
+                     ! Compatibility tests
+                     !---------------------
+                     ! Test 1: area conservation
+                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )
+                     IF ( zconv < epsi06 ) itest(1) = 1
+                     
+                     ! Test 2: volume conservation
+                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &
+                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
+                     IF ( zconv < epsi06 ) itest(2) = 1
+                     
+                     ! Test 3: thickness of the last category is in-bounds ?
+                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
+                     
+                     ! Test 4: positivity of ice concentrations
+                     itest(4) = 1
+                     DO jl = 1, i_fill
+                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0
+                     END DO
+                     !                                      !============================
+                  END DO                                    ! end iteration on categories
+                  !                                         !============================
+                  !
+                  IF( lwp .AND. SUM(itest) /= 4 ) THEN 
+                     WRITE(numout,*)
+                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
+                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
+                     WRITE(numout,*)
+                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
+                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
+                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
+                  ENDIF
+               
+               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp
+               !
+            END DO   
+         END DO   
+
+         !---------------------------------------------------------------------
+         ! 3.3) Space-dependent arrays for ice state variables
+         !---------------------------------------------------------------------
+
+         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
+         DO jl = 1, jpl ! loop over categories
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
+                  ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
+                  sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
+                  o_i(ji,jj,jl)   = 0._wp                                                     ! age (0 day)
+                  t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
+
+                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
+                    ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
+                  ELSE
+                    ht_s(ji,jj,jl)= 0._wp
+                  ENDIF
+
+                  ! This case below should not be used if (ht_s/ht_i) is ok in namelist
+                  ! In case snow load is in excess that would lead to transformation from snow to ice
+                  ! Then, transfer the snow excess into the ice (different from icethd_dh)
+                  zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
+                  ! recompute ht_i, ht_s avoiding out of bounds values
+                  ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
+                  ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
+
+                  ! ice volume, salt content, age content
+                  v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
+                  v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
+                  smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
+                  oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
+               END DO
+            END DO
+         END DO
+
+         ! for constant salinity in time
+         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
+            CALL ice_var_salprof
+            smv_i = sm_i * v_i
+         ENDIF
+            
+         ! Snow temperature and heat content
+         DO jk = 1, nlay_s
+            DO jl = 1, jpl ! loop over categories
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
+                     ! Snow energy of melting
+                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
+
+                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
+                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
+                  END DO
+               END DO
+            END DO
+         END DO
+
+         ! Ice salinity, temperature and heat content
+         DO jk = 1, nlay_i
+            DO jl = 1, jpl ! loop over categories
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
+                     s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
+                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
+
+                     ! heat content per unit volume
+                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
+                        -   rcp     * ( ztmelts - rt0 ) )
+
+                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
+                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
+                  END DO
+               END DO
+            END DO
+         END DO
+
+         tn_ice (:,:,:) = t_su (:,:,:)
+
+         ! MV MP 2016
+         ! Melt pond volume and fraction
+         IF ( ln_pnd ) THEN
+
+            DO jl = 1, jpl
+
+               a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:)
+               h_ip     (:,:,jl) = 0.05 * zswitch(:,:)
+               a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
+               v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl)
+
+            END DO
+
+         ELSE
+
+            a_ip(:,:,:)      = 0._wp
+            v_ip(:,:,:)      = 0._wp
+            a_ip_frac(:,:,:) = 0._wp
+            h_ip     (:,:,:) = 0._wp
+
+         ENDIF
+         ! END MV MP 2016
+
+      ELSE ! if ln_limini=false
+         a_i  (:,:,:) = 0._wp
+         v_i  (:,:,:) = 0._wp
+         v_s  (:,:,:) = 0._wp
+         smv_i(:,:,:) = 0._wp
+         oa_i (:,:,:) = 0._wp
+         ht_i (:,:,:) = 0._wp
+         ht_s (:,:,:) = 0._wp
+         sm_i (:,:,:) = 0._wp
+         o_i  (:,:,:) = 0._wp
+
+         e_i(:,:,:,:) = 0._wp
+         e_s(:,:,:,:) = 0._wp
+
+         DO jl = 1, jpl
+            DO jk = 1, nlay_i
+               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
+            END DO
+            DO jk = 1, nlay_s
+               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
+            END DO
+         END DO
+
+         a_ip(:,:,:)      = 0._wp
+         v_ip(:,:,:)      = 0._wp
+         a_ip_frac(:,:,:) = 0._wp
+         h_ip     (:,:,:) = 0._wp
+
+      ENDIF ! ln_limini
+      
+      at_i (:,:) = 0.0_wp
+      DO jl = 1, jpl
+         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
+      END DO
+
+      !--------------------------------------------------------------------
+      ! 4) Global ice variables for output diagnostics                    | 
+      !--------------------------------------------------------------------
+      u_ice (:,:)     = 0._wp
+      v_ice (:,:)     = 0._wp
+      stress1_i(:,:)  = 0._wp
+      stress2_i(:,:)  = 0._wp
+      stress12_i(:,:) = 0._wp
+
+      !--------------------------------------------------------------------
+      ! 5) Moments for advection
+      !--------------------------------------------------------------------
+
+      sxopw (:,:) = 0._wp 
+      syopw (:,:) = 0._wp 
+      sxxopw(:,:) = 0._wp 
+      syyopw(:,:) = 0._wp 
+      sxyopw(:,:) = 0._wp
+
+      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
+      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
+      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
+      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
+      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
+
+      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
+      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
+      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
+      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
+      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
+
+      sxsal  (:,:,:)  = 0._wp
+      sysal  (:,:,:)  = 0._wp
+      sxxsal (:,:,:)  = 0._wp
+      syysal (:,:,:)  = 0._wp
+      sxysal (:,:,:)  = 0._wp
+
+      sxage  (:,:,:)  = 0._wp
+      syage  (:,:,:)  = 0._wp
+      sxxage (:,:,:)  = 0._wp
+      syyage (:,:,:)  = 0._wp
+      sxyage (:,:,:)  = 0._wp
+
+      !------------------------------------
+      ! 6) store fields at before time-step
+      !------------------------------------
+      ! it is only necessary for the 1st interpolation by Agrif
+      a_i_b  (:,:,:)   = a_i  (:,:,:)
+      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
+      v_i_b  (:,:,:)   = v_i  (:,:,:)
+      v_s_b  (:,:,:)   = v_s  (:,:,:)
+      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
+      smv_i_b(:,:,:)   = smv_i(:,:,:)
+      oa_i_b (:,:,:)   = oa_i (:,:,:)
+      u_ice_b(:,:)     = u_ice(:,:)
+      v_ice_b(:,:)     = v_ice(:,:)
+
+      ! MV MP 2016
+      IF ( nn_pnd_scheme > 0 ) THEN
+         sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
+         syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
+         sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp 
+         syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp 
+         sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp
+      ENDIF
+      ! END MV MP 2016
+
+!!!clem
+!!      ! Output the initial state and forcings
+!!      CALL dia_wri_state( 'output.init', nit000 )
+!!!      
+
+   END SUBROUTINE ice_ist
+
+   SUBROUTINE ice_ist_init
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_ist_init  ***
+      !!        
+      !! ** Purpose : Definition of initial state of the ice 
+      !!
+      !! ** Method : Read the namiceini namelist and check the parameter 
+      !!       values called at the first timestep (nit000)
+      !!
+      !! ** input : 
+      !!        Namelist namiceini
+      !!
+      !! history :
+      !!  8.5  ! 03-08 (C. Ethe) original code 
+      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
+      !!-----------------------------------------------------------------------------
+      !
+      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
+      INTEGER :: ji,jj
+      INTEGER :: ifpr, ierror
+      !
+      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
+      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
+      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
+      !
+      NAMELIST/namiceini/ ln_limini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
+         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
+         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
+         &                sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
+      !!-----------------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
+      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
+      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namiceini )
+
+      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
+      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
+      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
+
+      ! Define the initial parameters
+      ! -------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_ist_init : ice parameters inititialisation '
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_limini     = ', ln_limini
+         WRITE(numout,*) '   ice initialization from a netcdf file      ln_limini_file  = ', ln_limini_file
+         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
+         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
+         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
+         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
+         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
+         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
+         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
+         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
+         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
+         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
+         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
+      ENDIF
+
+      IF( ln_limini_file ) THEN                      ! Ice initialization using input file
+         !
+         ! set si structure
+         ALLOCATE( si(jpfldi), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'Ice_ini in iceist: unable to allocate si structure' )   ;   RETURN
+         ENDIF
+
+         DO ifpr = 1, jpfldi
+            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
+            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
+         END DO
+
+         ! fill si with slf_i and control print
+         CALL fld_fill( si, slf_i, cn_dir, 'ice_ist', 'lim istate ini', 'numnam_ice' )
+
+         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
+
+      ENDIF
+
+   END SUBROUTINE ice_ist_init
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :         Empty module          NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE ice_ist          ! Empty routine
+   END SUBROUTINE ice_ist
+#endif
+
+   !!======================================================================
+END MODULE iceist
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8424)
+++ 	(revision )
@@ -1,599 +1,0 @@
-MODULE iceistate
-   !!======================================================================
-   !!                     ***  MODULE  iceistate  ***
-   !!              Initialisation of diagnostics ice variables
-   !!======================================================================
-   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
-   !!            3.0  ! 2011-02 (G. Madec) dynamical allocation
-   !!             -   ! 2014    (C. Rousset) add N/S initializations
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                    LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   ice_istate      :  Initialisation of diagnostics ice variables
-   !!   ice_istate_init :  initialization of ice state and namelist read
-   !!----------------------------------------------------------------------
-   USE phycst         ! physical constant
-   USE oce            ! dynamics and tracers variables
-   USE dom_oce        ! ocean domain
-   USE sbc_oce , ONLY : sst_m, sss_m 
-   USE sbc_ice , ONLY : tn_ice
-   USE eosbn2         ! equation of state
-   USE ice            ! sea-ice variables
-   USE par_oce        ! ocean parameters
-   USE icevar         ! ice_var_salprof
-   !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE fldread        ! read input fields
-   USE iom
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_istate      ! routine called by lim_init.F90
-
-   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
-   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness (m)    at T-point
-   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snow thicknes (m)    at T-point
-   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point
-   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point
-   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point
-   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
-   !!----------------------------------------------------------------------
-   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
-   !! $Id: iceistate.F90 8378 2017-07-26 13:55:59Z clem $
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_istate
-      !!-------------------------------------------------------------------
-      !!                    ***  ROUTINE ice_istate  ***
-      !!
-      !! ** Purpose :   defined the sea-ice initial state
-      !!
-      !! ** Method  :   This routine will put some ice where ocean
-      !!                is at the freezing point, then fill in ice 
-      !!                state variables using prescribed initial 
-      !!                values in the namelist            
-      !!
-      !! ** Steps   :   1) Read namelist
-      !!                2) Basal temperature; ice and hemisphere masks
-      !!                3) Fill in the ice thickness distribution using gaussian
-      !!                4) Fill in space-dependent arrays for state variables
-      !!                5) Diagnostic arrays
-      !!                6) Lateral boundary conditions
-      !!
-      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
-      !!              where there is no ice (clem: I do not know why, is it mandatory?) 
-      !!
-      !! History :
-      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
-      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
-      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
-      !!--------------------------------------------------------------------
-      !! * Local variables
-      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
-      REAL(wp)   :: ztmelts, zdh
-      INTEGER    :: i_hemis, i_fill, jl0
-      REAL(wp)   :: zarg, zV, zconv, zdv
-      REAL(wp), DIMENSION(jpi,jpj)     :: zswitch    ! ice indicator
-      REAL(wp), DIMENSION(jpi,jpj)     :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj)     :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini, za_i_ini                         !data by cattegories to fill
-      INTEGER , DIMENSION(4)           :: itest
-      !--------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'ice_istate : sea-ice initialization '
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~ '
-
-      !--------------------------------------------------------------------
-      ! 1) Read namelist
-      !--------------------------------------------------------------------
-      !
-      CALL ice_istate_init
-
-      ! init surface temperature
-      DO jl = 1, jpl
-         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
-         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
-      END DO
-
-      ! init basal temperature (considered at freezing point)
-      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
-      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
-
-
-      !--------------------------------------------------------------------
-      ! 2) Initialization of sea ice state variables
-      !--------------------------------------------------------------------
-      IF( ln_limini ) THEN
-         !
-         IF( ln_limini_file )THEN
-         !
-            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1)
-            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1)
-            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1)
-            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1)
-            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1)
-            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1)
-            !
-            WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 
-            ELSEWHERE                       ; zswitch(:,:) = 0._wp
-            END WHERE
-            !
-         ELSE ! ln_limini_file = F
-
-            !--------------------------------------------------------------------
-            ! 3) Basal temperature, ice mask
-            !--------------------------------------------------------------------
-            ! no ice if sst <= t-freez + ttest
-            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 
-            ELSEWHERE                                                                  ; zswitch(:,:) = tmask(:,:,1)
-            END WHERE
-
-            !-----------------------------
-            ! 3.1) Hemisphere-dependent arrays
-            !-----------------------------
-            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( ff_t(ji,jj) >= 0._wp ) THEN
-                     zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj)
-                     zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj)
-                     zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj)
-                     zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
-                     zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj)
-                     ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
-                  ELSE
-                     zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj)
-                     zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj)
-                     zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj)
-                     zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
-                     zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj)
-                     ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
-                  ENDIF
-               END DO
-            END DO
-            !
-         ENDIF ! ln_limini_file
-         
-         zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)   ! ice volume
-         !---------------------------------------------------------------------
-         ! 3.2) Distribute ice concentration and thickness into the categories
-         !---------------------------------------------------------------------
-         ! a gaussian distribution for ice concentration is used
-         ! then we check whether the distribution fullfills
-         ! volume and area conservation, positivity and ice categories bounds
-         zh_i_ini(:,:,:) = 0._wp 
-         za_i_ini(:,:,:) = 0._wp
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
-
-                  !--- jl0: most likely index where cc will be maximum
-                  jl0 = jpl
-                  DO jl = 1, jpl
-                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
-                        jl0 = jl
-                        CYCLE
-                     ENDIF
-                  END DO
-                  !
-                  ! initialisation of tests
-                  itest(:)  = 0
-                  
-                  i_fill = jpl + 1                                             !====================================
-                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
-                     ! iteration                                               !====================================
-                     i_fill = i_fill - 1
-
-                     ! initialisation of ice variables for each try
-                     zh_i_ini(ji,jj,:) = 0._wp 
-                     za_i_ini(ji,jj,:) = 0._wp
-                     itest(:) = 0
-                     !
-                     ! *** case very thin ice: fill only category 1
-                     IF ( i_fill == 1 ) THEN
-                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
-                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
-
-                     ! *** case ice is thicker: fill categories >1
-                     ELSE
-
-                        ! Fill ice thicknesses in the (i_fill-1) cat by hmean 
-                        DO jl = 1, i_fill-1
-                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
-                        END DO
-                        !
-                        !--- Concentrations
-                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
-                        DO jl = 1, i_fill - 1
-                           IF( jl /= jl0 )THEN
-                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
-                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
-                           ENDIF
-                        END DO
-                        !
-                        ! Concentration in the last (i_fill) category
-                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
-
-                        ! Ice thickness in the last (i_fill) category
-                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
-                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
-
-                        ! clem: correction if concentration of upper cat is greater than lower cat
-                        !       (it should be a gaussian around jl0 but sometimes it is not)
-                        IF ( jl0 /= jpl ) THEN
-                           DO jl = jpl, jl0+1, -1
-                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
-                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
-                                 zh_i_ini(ji,jj,jl    ) = 0._wp
-                                 za_i_ini(ji,jj,jl    ) = 0._wp
-                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
-                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
-                              END IF
-                           ENDDO
-                        ENDIF
-                        !
-                     ENDIF ! case ice is thick or thin
-
-                     !---------------------
-                     ! Compatibility tests
-                     !---------------------
-                     ! Test 1: area conservation
-                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )
-                     IF ( zconv < epsi06 ) itest(1) = 1
-                     
-                     ! Test 2: volume conservation
-                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &
-                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
-                     IF ( zconv < epsi06 ) itest(2) = 1
-                     
-                     ! Test 3: thickness of the last category is in-bounds ?
-                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
-                     
-                     ! Test 4: positivity of ice concentrations
-                     itest(4) = 1
-                     DO jl = 1, i_fill
-                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0
-                     END DO
-                     !                                      !============================
-                  END DO                                    ! end iteration on categories
-                  !                                         !============================
-                  !
-                  IF( lwp .AND. SUM(itest) /= 4 ) THEN 
-                     WRITE(numout,*)
-                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
-                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
-                     WRITE(numout,*)
-                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
-                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
-                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
-                  ENDIF
-               
-               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp
-               !
-            END DO   
-         END DO   
-
-         !---------------------------------------------------------------------
-         ! 3.3) Space-dependent arrays for ice state variables
-         !---------------------------------------------------------------------
-
-         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
-         DO jl = 1, jpl ! loop over categories
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
-                  ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
-                  sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
-                  o_i(ji,jj,jl)   = 0._wp                                                     ! age (0 day)
-                  t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
-
-                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
-                    ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
-                  ELSE
-                    ht_s(ji,jj,jl)= 0._wp
-                  ENDIF
-
-                  ! This case below should not be used if (ht_s/ht_i) is ok in namelist
-                  ! In case snow load is in excess that would lead to transformation from snow to ice
-                  ! Then, transfer the snow excess into the ice (different from icethd_dh)
-                  zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
-                  ! recompute ht_i, ht_s avoiding out of bounds values
-                  ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
-                  ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
-
-                  ! ice volume, salt content, age content
-                  v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
-                  v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
-                  smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
-                  oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
-               END DO
-            END DO
-         END DO
-
-         ! for constant salinity in time
-         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
-            CALL ice_var_salprof
-            smv_i = sm_i * v_i
-         ENDIF
-            
-         ! Snow temperature and heat content
-         DO jk = 1, nlay_s
-            DO jl = 1, jpl ! loop over categories
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
-                     ! Snow energy of melting
-                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
-
-                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
-                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
-                  END DO
-               END DO
-            END DO
-         END DO
-
-         ! Ice salinity, temperature and heat content
-         DO jk = 1, nlay_i
-            DO jl = 1, jpl ! loop over categories
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
-                     s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
-                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
-
-                     ! heat content per unit volume
-                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
-                        -   rcp     * ( ztmelts - rt0 ) )
-
-                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
-                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
-                  END DO
-               END DO
-            END DO
-         END DO
-
-         tn_ice (:,:,:) = t_su (:,:,:)
-
-         ! MV MP 2016
-         ! Melt pond volume and fraction
-         IF ( ln_pnd ) THEN
-
-            DO jl = 1, jpl
-
-               a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:)
-               h_ip     (:,:,jl) = 0.05 * zswitch(:,:)
-               a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
-               v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl)
-
-            END DO
-
-         ELSE
-
-            a_ip(:,:,:)      = 0._wp
-            v_ip(:,:,:)      = 0._wp
-            a_ip_frac(:,:,:) = 0._wp
-            h_ip     (:,:,:) = 0._wp
-
-         ENDIF
-         ! END MV MP 2016
-
-      ELSE ! if ln_limini=false
-         a_i  (:,:,:) = 0._wp
-         v_i  (:,:,:) = 0._wp
-         v_s  (:,:,:) = 0._wp
-         smv_i(:,:,:) = 0._wp
-         oa_i (:,:,:) = 0._wp
-         ht_i (:,:,:) = 0._wp
-         ht_s (:,:,:) = 0._wp
-         sm_i (:,:,:) = 0._wp
-         o_i  (:,:,:) = 0._wp
-
-         e_i(:,:,:,:) = 0._wp
-         e_s(:,:,:,:) = 0._wp
-
-         DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
-            END DO
-            DO jk = 1, nlay_s
-               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
-            END DO
-         END DO
-
-         a_ip(:,:,:)      = 0._wp
-         v_ip(:,:,:)      = 0._wp
-         a_ip_frac(:,:,:) = 0._wp
-         h_ip     (:,:,:) = 0._wp
-
-      ENDIF ! ln_limini
-      
-      at_i (:,:) = 0.0_wp
-      DO jl = 1, jpl
-         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
-      END DO
-
-      !--------------------------------------------------------------------
-      ! 4) Global ice variables for output diagnostics                    | 
-      !--------------------------------------------------------------------
-      u_ice (:,:)     = 0._wp
-      v_ice (:,:)     = 0._wp
-      stress1_i(:,:)  = 0._wp
-      stress2_i(:,:)  = 0._wp
-      stress12_i(:,:) = 0._wp
-
-      !--------------------------------------------------------------------
-      ! 5) Moments for advection
-      !--------------------------------------------------------------------
-
-      sxopw (:,:) = 0._wp 
-      syopw (:,:) = 0._wp 
-      sxxopw(:,:) = 0._wp 
-      syyopw(:,:) = 0._wp 
-      sxyopw(:,:) = 0._wp
-
-      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
-      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
-      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
-      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
-      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
-
-      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
-      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
-      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
-      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
-      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
-
-      sxsal  (:,:,:)  = 0._wp
-      sysal  (:,:,:)  = 0._wp
-      sxxsal (:,:,:)  = 0._wp
-      syysal (:,:,:)  = 0._wp
-      sxysal (:,:,:)  = 0._wp
-
-      sxage  (:,:,:)  = 0._wp
-      syage  (:,:,:)  = 0._wp
-      sxxage (:,:,:)  = 0._wp
-      syyage (:,:,:)  = 0._wp
-      sxyage (:,:,:)  = 0._wp
-
-      !------------------------------------
-      ! 6) store fields at before time-step
-      !------------------------------------
-      ! it is only necessary for the 1st interpolation by Agrif
-      a_i_b  (:,:,:)   = a_i  (:,:,:)
-      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
-      v_i_b  (:,:,:)   = v_i  (:,:,:)
-      v_s_b  (:,:,:)   = v_s  (:,:,:)
-      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
-      smv_i_b(:,:,:)   = smv_i(:,:,:)
-      oa_i_b (:,:,:)   = oa_i (:,:,:)
-      u_ice_b(:,:)     = u_ice(:,:)
-      v_ice_b(:,:)     = v_ice(:,:)
-
-      ! MV MP 2016
-      IF ( nn_pnd_scheme > 0 ) THEN
-         sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
-         syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
-         sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp 
-         syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp 
-         sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp
-      ENDIF
-      ! END MV MP 2016
-
-!!!clem
-!!      ! Output the initial state and forcings
-!!      CALL dia_wri_state( 'output.init', nit000 )
-!!!      
-
-   END SUBROUTINE ice_istate
-
-   SUBROUTINE ice_istate_init
-      !!-------------------------------------------------------------------
-      !!                   ***  ROUTINE ice_istate_init  ***
-      !!        
-      !! ** Purpose : Definition of initial state of the ice 
-      !!
-      !! ** Method : Read the namiceini namelist and check the parameter 
-      !!       values called at the first timestep (nit000)
-      !!
-      !! ** input : 
-      !!        Namelist namiceini
-      !!
-      !! history :
-      !!  8.5  ! 03-08 (C. Ethe) original code 
-      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
-      !!-----------------------------------------------------------------------------
-      !
-      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
-      INTEGER :: ji,jj
-      INTEGER :: ifpr, ierror
-      !
-      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
-      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
-      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
-      !
-      NAMELIST/namiceini/ ln_limini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
-         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
-         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
-         &                sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
-      !!-----------------------------------------------------------------------------
-      !
-      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
-      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
-      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namiceini )
-
-      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
-      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
-      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
-
-      ! Define the initial parameters
-      ! -------------------------
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'ice_istate_init : ice parameters inititialisation '
-         WRITE(numout,*) '~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_limini     = ', ln_limini
-         WRITE(numout,*) '   ice initialization from a netcdf file      ln_limini_file  = ', ln_limini_file
-         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
-         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
-         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
-         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
-         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
-         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
-         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
-         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
-         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
-         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
-         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
-      ENDIF
-
-      IF( ln_limini_file ) THEN                      ! Ice initialization using input file
-         !
-         ! set si structure
-         ALLOCATE( si(jpfldi), STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'Ice_ini in iceistate: unable to allocate si structure' )   ;   RETURN
-         ENDIF
-
-         DO ifpr = 1, jpfldi
-            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
-            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
-         END DO
-
-         ! fill si with slf_i and control print
-         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'lim istate ini', 'numnam_ice' )
-
-         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
-
-      ENDIF
-
-   END SUBROUTINE ice_istate_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :         Empty module          NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE ice_istate          ! Empty routine
-   END SUBROUTINE ice_istate
-#endif
-
-   !!======================================================================
-END MODULE iceistate
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8426)
@@ -23,5 +23,5 @@
    USE ice1D            ! LIM-3 thermodynamic variables
    USE ice              ! LIM-3 variables
-   USE icecons          ! conservation tests
+   USE icectl           ! conservation tests
    USE icetab
    !
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8426)
@@ -19,5 +19,4 @@
    USE ice              ! LIM variables
    USE icevar           ! LIM
-   USE icecons          ! conservation tests
    USE icectl           ! control prints
    !
@@ -78,5 +77,5 @@
 
 
-   SUBROUTINE ice_rdgrft
+   SUBROUTINE ice_rdgrft( kt )
       !!---------------------------------------------------------------------!
       !!                ***  ROUTINE ice_rdgrft ***
@@ -104,4 +103,6 @@
       !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged
       !!--------------------------------------------------------------------!
+      INTEGER, INTENT(in) ::   kt     ! number of iteration
+      !!
       INTEGER  ::   ji, jj, jk, jl        ! dummy loop index
       INTEGER  ::   niter                 ! local integer 
@@ -120,4 +121,10 @@
       !!-----------------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('icerdgrft')
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'icerdgrft'
+         WRITE(numout,*)'~~~~~~~~~'
+      ENDIF
 
       ! conservation test
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8426)
@@ -19,5 +19,4 @@
    USE ice              ! LIM-3 variables
    USE icerhg_evp       ! EVP rheology
-   USE icecons          ! conservation tests
    USE icectl           ! control prints
    USE icevar
@@ -61,4 +60,10 @@
 
       IF( nn_timing == 1 )  CALL timing_start('icerhg')
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'icerhg'
+         WRITE(numout,*)'~~~~~~'
+      ENDIF
 
       CALL ice_var_agg(1)   ! aggregate ice categories
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8426)
@@ -11,5 +11,5 @@
    !!            3.3  ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
    !!            3.4  ! 2011-01  (A Porter)  dynamical allocation
-   !!             -   ! 2012-10  (C. Rousset)  add ice_diahsb
+   !!             -   ! 2012-10  (C. Rousset)  add ice_dia
    !!            3.6  ! 2014-07  (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface
    !!            4.0  ! 2016-06  (L. Brodeau) new unified bulk routine (based on AeroBulk)
@@ -37,9 +37,8 @@
    USE icerdgrft       ! Ice ridging/rafting
    USE iceupdate       ! sea surface boundary condition
-   USE icediahsb       ! Ice budget diagnostics
+   USE icedia          ! Ice budget diagnostics
    USE icewri          ! Ice outputs
    USE icerst          ! Ice restarts
-   USE iceerr1         ! Ice corrections after dynamics
-   USE iceerr2         ! Ice corrections after thermo
+   USE icecor          ! Ice corrections
    USE icevar          ! Ice variables switch
    USE icectl          !
@@ -47,5 +46,5 @@
    USE limmp
    ! END MV MP 2016
-   USE iceistate       ! LIM initial state
+   USE iceist          ! LIM initial state
    USE icethd_sal      ! LIM ice thermodynamics: salinity
    !
@@ -130,9 +129,4 @@
          CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
          t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
-
-         ! Mask sea ice surface temperature (set to rt0 over land)
-         DO jl = 1, jpl
-            t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
-         END DO
          !
                                       CALL ice_bef         ! Store previous ice values
@@ -156,10 +150,10 @@
                                       CALL ice_adv( kt )       ! -- advection
             IF( nn_limdyn == 2 .AND. nn_monocat /= 2 )  &      ! -- ridging/rafting
-               &                      CALL ice_rdgrft         
-            IF( nn_limdyn == 2 )      CALL ice_err1( kt )      ! -- Corrections
+               &                      CALL ice_rdgrft( kt )         
+            IF( nn_limdyn == 2 )      CALL ice_cor( kt , 1 )   ! -- Corrections
             !
          ENDIF
-
          ! ---
+         
 #if defined key_agrif
          IF( .NOT. Agrif_Root() )     CALL agrif_interp_lim3('T')
@@ -196,5 +190,5 @@
          ! END MV MP 2016
 
-         IF( ln_limthd )              CALL ice_err2( kt )       ! -- Corrections
+         IF( ln_limthd )              CALL ice_cor( kt , 2 )    ! -- Corrections
          ! ---
 # if defined key_agrif
@@ -213,5 +207,5 @@
 !!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
 !!# endif
-         IF( ln_limdiahsb )           CALL ice_diahsb( kt )     ! -- Diagnostics and outputs 
+         IF( ln_limdiahsb )           CALL ice_dia( kt )     ! -- Diagnostics and outputs 
          !
                                       CALL ice_wri( 1 )         ! -- Ice outputs 
@@ -282,5 +276,5 @@
       !                                ! Initial sea-ice state
       IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
-         CALL ice_istate
+         CALL ice_ist
       ELSE                                    ! start from a restart file
          CALL ice_rst_read
@@ -291,5 +285,5 @@
       CALL ice_update_init                 ! ice surface boundary condition
       !
-      IF( ln_limdiahsb) CALL ice_diahsb_init  ! initialization for diags
+      IF( ln_limdiahsb) CALL ice_dia_init  ! initialization for diags
       !
       fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8426)
@@ -36,5 +36,4 @@
    USE icetab         ! 1D <==> 2D transformation
    USE icevar         !
-   USE icecons        ! conservation tests
    USE icectl         ! control print
    !
@@ -93,7 +92,7 @@
 
       IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)'' 
-         WRITE(numout,*)' ice_thd '
-         WRITE(numout,*)' ~~~~~~~~'
+         WRITE(numout,*)
+         WRITE(numout,*)' icethd '
+         WRITE(numout,*)' ~~~~~~~'
       ENDIF
       
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90	(revision 8426)
@@ -23,5 +23,5 @@
    USE ice            ! LIM variables
    USE icetab         ! LIM 2D <==> 1D
-   USE icecons        ! LIM conservation
+   USE icectl         ! LIM conservation
    USE icethd_ent
    USE icevar
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8426)
@@ -32,9 +32,8 @@
    USE sbc_oce , ONLY : nn_fsbc, ln_ice_embd, sfx, fr_i, qsr_tot, qns, qsr, fmmflx, emp, taum, utau, vtau
    USE sbccpl         ! Surface boundary condition: coupled interface
-   USE icealbedo      ! albedo parameters
+   USE icealb         ! albedo parameters
    USE traqsr         ! add penetration of solar flux in the calculation of heat budget
    USE domvvl         ! Variable volume
    USE icectl         !
-   USE icecons        !
    USE bdy_oce  , ONLY: ln_bdy
    !
@@ -44,4 +43,5 @@
    USE lib_mpp        ! MPP library
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -108,4 +108,11 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace
       !!---------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('ice_update_flx')
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'ice_update_flx'
+         WRITE(numout,*)'~~~~~~~~~~~~~~'
+      ENDIF
 
       ! --- case we bypass ice thermodynamics --- !
@@ -201,5 +208,5 @@
       !    Snow/ice albedo (only if sent to coupler, useless in forced mode)   !
       !------------------------------------------------------------------------!
-      CALL ice_albedo( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
+      CALL ice_alb( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
 
       alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
@@ -211,4 +218,6 @@
       IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt, 3, ' - Final state ice_update - ' )
       IF( ln_ctl )      CALL ice_prt3D( 'iceupdate' )
+
+      IF( nn_timing == 1 )  CALL timing_stop('ice_update_flx')
 
    END SUBROUTINE ice_update_flx
@@ -247,4 +256,13 @@
       REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
       !!---------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('ice_update_tau')
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'ice_update_tau'
+         WRITE(numout,*)'~~~~~~~~~~~~~~'
+      ENDIF
+
       zrhoco = rau0 * rn_cio
       !
@@ -285,4 +303,5 @@
       CALL lbc_lnk_multi( utau, 'U', -1., vtau, 'V', -1. )   ! lateral boundary condition
       !
+      IF( nn_timing == 1 )  CALL timing_stop('ice_update_tau')
       !  
    END SUBROUTINE ice_update_tau
@@ -304,5 +323,5 @@
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'ice_update_init : LIM-3 sea-ice - surface boundary condition'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~   '
 
       !                                      ! allocate ice_update array
@@ -368,11 +387,4 @@
    END SUBROUTINE ice_update_init
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :        Dummy module       NO LIM 3.0 sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE ice_update           ! Dummy routine
-   END SUBROUTINE ice_update
 #endif 
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8426)
@@ -612,5 +612,5 @@
       !!                  ice thickness distribution follows a gaussian law
       !!               around the concentration of the most likely ice thickness
-      !!                           (similar as iceistate.F90)
+      !!                           (similar as iceist.F90)
       !!
       !! ** Method:   Iterative procedure
