Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 9270)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 9271)
@@ -475,8 +475,8 @@
                   WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl)
                   WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl)
-                  WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)
+                  WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1:nlay_s,jl)
                   WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)
                   WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl)
-                  WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl)
+                  WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1:nlay_s,jl)
                   WRITE(numout,*) ' t_i           : ', t_i(ji,jj,1:nlay_i,jl)
                   WRITE(numout,*) ' s_i           : ', s_i(ji,jj,jl)
@@ -524,5 +524,5 @@
                   WRITE(numout,*) ' h_i        : ', h_i(ji,jj,jl)              , ' h_s        : ', h_s(ji,jj,jl)
                   WRITE(numout,*) ' t_i        : ', t_i(ji,jj,1:nlay_i,jl)
-                  WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl)
+                  WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1:nlay_s,jl)
                   WRITE(numout,*) ' s_i        : ', s_i(ji,jj,jl)              , ' o_i        : ', o_i(ji,jj,jl)
                   WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' a_i_b      : ', a_i_b(ji,jj,jl)   
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90	(revision 9270)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90	(revision 9271)
@@ -38,8 +38,8 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsn , sysn , sxxsn , syysn , sxysn    ! snow thickness
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxa  , sya  , sxxa  , syya  , sxya     ! lead fraction
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow thermal content
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow layers heat content
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction
@@ -87,6 +87,7 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0smi, z0oi
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0es
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei
       !!----------------------------------------------------------------------
@@ -94,8 +95,7 @@
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
       !
-      ALLOCATE( zarea(jpi,jpj)     , z0opw(jpi,jpj, 1 ) ,                                           &
-         &      z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , z0ai(jpi,jpj,jpl) , z0es(jpi,jpj,jpl) ,   &
-         &      z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap(jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   &
-         &      z0ei (jpi,jpj,nlay_i,jpl)                                                           )
+      ALLOCATE( zarea(jpi,jpj)    , z0opw(jpi,jpj, 1 ) , z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) ,                       &
+         &      z0ai(jpi,jpj,jpl) , z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap (jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   &
+         &      z0es (jpi,jpj,nlay_s,jpl), z0ei(jpi,jpj,nlay_i,jpl) )
       !
       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !        
@@ -119,5 +119,7 @@
          z0smi(:,:,jl) = psv_i(:,:,  jl) * e1e2t(:,:)     ! Salt content
          z0oi (:,:,jl) = poa_i(:,:,  jl) * e1e2t(:,:)     ! Age content
-         z0es (:,:,jl) = pe_s (:,:,1,jl) * e1e2t(:,:)     ! Snow heat content
+         DO jk = 1, nlay_s
+            z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
+         END DO
          DO jk = 1, nlay_i
             z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
@@ -158,15 +160,15 @@
                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
-                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
-                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+               DO jk = 1, nlay_s                                                               !--- snow heat contents ---
+                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
+                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
+                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
+                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
+               END DO
                DO jk = 1, nlay_i                                                               !--- ice heat contents ---
-                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
+                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
+                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                END DO
                IF ( ln_pnd_H12 ) THEN
@@ -211,15 +213,15 @@
                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
-                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
-                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+               DO jk = 1, nlay_s                                                             !--- snow heat contents ---
+                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
+                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
+                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
+                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
+               END DO
                DO jk = 1, nlay_i                                                             !--- ice heat contents ---
-                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
+                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
+                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                END DO
                IF ( ln_pnd_H12 ) THEN
@@ -247,5 +249,7 @@
          poa_i(:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:)
          pa_i (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:)
-         pe_s (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:)
+         DO jk = 1, nlay_s
+            pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:)
+         END DO
          DO jk = 1, nlay_i
             pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:)
@@ -257,5 +261,5 @@
       END DO
       !
-      DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0es , z0smi , z0oi , z0ap , z0vp , z0ei )
+      DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0smi , z0oi , z0ap , z0vp , z0es, z0ei )
       !
       IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file
@@ -625,11 +629,14 @@
          &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
          &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
-         &      sxc0 (jpi,jpj,jpl) , syc0 (jpi,jpj,jpl) , sxxc0 (jpi,jpj,jpl) , syyc0 (jpi,jpj,jpl) , sxyc0 (jpi,jpj,jpl) ,   &
          &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
          &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
          &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
          &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
-         &      sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe(jpi,jpj,nlay_i,jpl) , &
-         &      syye(jpi,jpj,nlay_i,jpl) , sxye(jpi,jpj,nlay_i,jpl)                            , &
+         !
+         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
+         &      syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl)                             , &
+         !
+         &      sxe  (jpi,jpj,nlay_i,jpl) , sye  (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
+         &      syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl)                             , &
          &      STAT = ierr )
       !
@@ -689,10 +696,4 @@
             CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   )
             CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   )
-            !                                                        ! snow thermal content
-            CALL iom_get( numrir, jpdom_autoglo, 'sxc0'  , sxc0   )
-            CALL iom_get( numrir, jpdom_autoglo, 'syc0'  , syc0   )
-            CALL iom_get( numrir, jpdom_autoglo, 'sxxc0' , sxxc0  )
-            CALL iom_get( numrir, jpdom_autoglo, 'syyc0' , syyc0  )
-            CALL iom_get( numrir, jpdom_autoglo, 'sxyc0' , sxyc0  )
             !                                                        ! ice salinity
             CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal  )
@@ -713,6 +714,15 @@
             CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw )
             CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw )
+            !                                                        ! snow layers heat content
+            DO jk = 1, nlay_s
+               WRITE(zchar1,'(I2.2)') jk
+               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
+               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
+               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
+            END DO
             !                                                        ! ice layers heat content
-            DO jk = 1, nlay_i 
+            DO jk = 1, nlay_i
                WRITE(zchar1,'(I2.2)') jk
                znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
@@ -744,8 +754,8 @@
             sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
             sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! lead fraction
-            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow thermal content
             sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
             sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
             sxopw = 0._wp   ;   syopw = 0._wp   ;   sxxopw = 0._wp   ;   syyopw = 0._wp   ;   sxyopw = 0._wp      ! open water in sea ice
+            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
             sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
             IF( ln_pnd_H12 ) THEN
@@ -783,10 +793,4 @@
          CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
          CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
-         !                                                           ! snow thermal content
-         CALL iom_rstput( iter, nitrst, numriw, 'sxc0'  , sxc0   )
-         CALL iom_rstput( iter, nitrst, numriw, 'syc0'  , syc0   )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxc0' , sxxc0  )
-         CALL iom_rstput( iter, nitrst, numriw, 'syyc0' , syyc0  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxyc0' , sxyc0  )
          !                                                           ! ice salinity
          CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
@@ -807,6 +811,15 @@
          CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw )
          CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw )
+         !                                                           ! snow layers heat content
+         DO jk = 1, nlay_s
+            WRITE(zchar1,'(I2.2)') jk
+            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+         END DO
          !                                                           ! ice layers heat content
-         DO jk = 1, nlay_i 
+         DO jk = 1, nlay_i
             WRITE(zchar1,'(I2.2)') jk
             znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90	(revision 9270)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90	(revision 9271)
@@ -122,5 +122,7 @@
             END DO
             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content
+            DO jk = 1, nlay_s
+               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,jl) )   ! Snow heat content
+            END DO
             IF ( ln_pnd_H12 ) THEN
                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 9270)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 9271)
@@ -570,6 +570,6 @@
       REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
       REAL(wp)                  ::   airft1, oirft1, aprft1
-      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, esrdg  ! area etc of new ridges
-      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, esrft  ! area etc of rafted ice
+      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges
+      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice
       !
       REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
@@ -577,5 +577,7 @@
       REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a
       !
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice energy of rafting ice
+      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice
+      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice  energy of rafting ice
+      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrdg     ! enth*volume of new ridges      
       REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges
       !!-------------------------------------------------------------------
@@ -633,5 +635,4 @@
                oirdg1     = oa_i_2d(ji,jl1)   * afrdg
                oirdg2(ji) = oa_i_2d(ji,jl1)   * afrdg * hi_hrdg(ji,jl1) 
-               esrdg(ji)  = ze_s_2d(ji,1,jl1) * afrdg
 
                virft(ji)  = v_i_2d (ji,jl1)   * afrft
@@ -640,5 +641,4 @@
                oirft1     = oa_i_2d(ji,jl1)   * afrft 
                oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
-               esrft(ji)  = ze_s_2d(ji,1,jl1) * afrft
 
                IF ( ln_pnd_H12 ) THEN
@@ -663,7 +663,4 @@
                wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean
                   &                                      + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
-
-               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
-                  &                                - esrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
 
                ! Put the melt pond water into the ocean
@@ -697,9 +694,27 @@
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
                ENDIF
-               ze_s_2d(ji,1,jl1) = ze_s_2d(ji,1,jl1) - esrdg (ji) - esrft (ji)
             ENDIF
 
          END DO ! ji
 
+         ! special loop for e_s because of layers jk
+         DO jk = 1, nlay_s
+            DO ji = 1, npti
+               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+                  ! Compute ridging /rafting fractions
+                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
+                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
+                  ! Compute ridging ice and new ridges for es
+                  esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
+                  esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
+                  ! Put the snow lost by ridging into the ocean
+                  hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
+                     &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
+                  ! Update jl1
+                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
+               ENDIF
+            END DO
+         END DO
+                  
          ! special loop for e_i because of layers jk
          DO jk = 1, nlay_i
@@ -754,5 +769,5 @@
                   sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
                   v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
-                     &                                    vsrft (ji) * rn_fsnwrft * zswitch(ji) )
+                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) )
                   IF ( ln_pnd_H12 ) THEN
                      v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
@@ -761,11 +776,17 @@
                         &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
                   ENDIF
-                  ze_s_2d(ji,1,jl2) = ze_s_2d(ji,1,jl2) + ( esrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
-                     &                                      esrft (ji) * rn_fsnwrft * zswitch(ji) )
                   
                ENDIF
 
             END DO
-
+            ! for snow enthalpy
+            DO jk = 1, nlay_s
+               DO ji = 1, npti
+                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
+                     &   ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji)  +  &
+                     &                                               esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
+               END DO
+            END DO
+            ! for ice enthalpy
             DO jk = 1, nlay_i
                DO ji = 1, npti
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 9270)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 9271)
@@ -132,9 +132,12 @@
       CALL iom_rstput( iter, nitrst, numriw, 'a_ip', a_ip )
       CALL iom_rstput( iter, nitrst, numriw, 'v_ip', v_ip )
-!!gm dangerous !!!!!  ===>>>> better reading writing all snow layers !
-      ! Snow enthalpy (1st snow layer only)
-      z3d = e_s(:,:,1,:)
-      CALL iom_rstput( iter, nitrst, numriw, 'e_s_l01' , z3d )
-      ! Ice enthalpy (all ice layers)
+      ! Snow enthalpy
+      DO jk = 1, nlay_s 
+         WRITE(zchar1,'(I2.2)') jk
+         znam = 'e_s'//'_l'//zchar1
+         z3d(:,:,:) = e_s(:,:,jk,:)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+      END DO
+      ! Ice enthalpy
       DO jk = 1, nlay_i 
          WRITE(zchar1,'(I2.2)') jk
@@ -219,10 +222,13 @@
          v_ip(:,:,:) = 0._wp
       ENDIF
-!!gm dangerous !!!!!  ===>>>> better reading writing all snow layers !
-      ! Snow enthalpy (1st snow layer only)
-      CALL iom_get( numrir, jpdom_autoglo, 'e_s_l01' , z3d )
-      e_s(:,:,1,:) = z3d
-      ! Ice enthalpy (all ice layers)
-      DO jk = 1, nlay_i 
+      ! Snow enthalpy
+      DO jk = 1, nlay_s
+         WRITE(zchar1,'(I2.2)') jk
+         znam = 'e_s'//'_l'//zchar1
+         CALL iom_get( numrir, jpdom_autoglo, znam , z3d )
+         e_s(:,:,jk,:) = z3d(:,:,:)
+      END DO
+      ! Ice enthalpy
+      DO jk = 1, nlay_i
          WRITE(zchar1,'(I2.2)') jk
          znam = 'e_i'//'_l'//zchar1
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 9270)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 9271)
@@ -98,9 +98,9 @@
       REAL(wp), DIMENSION(jpij) ::   zdh_s_sub   ! snow sublimation
 
+      REAL(wp), DIMENSION(jpij,nlay_s) ::   zh_s      ! snw layer thickness
+      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
       REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
       INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting 
 
-      REAL(wp), DIMENSION(jpij) ::   zeh_i       ! total ice heat content  (J.m-2)
       REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing
 
@@ -161,31 +161,40 @@
       END DO
 
+      !--------------------------------------------!
+      !  2) Computing layer thicknesses            !
+      !--------------------------------------------!
+      DO jk = 1, nlay_i
+         DO ji = 1, npti
+            zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
+         END DO
+      END DO
+
+      DO jk = 1, nlay_s
+         DO ji = 1, npti
+            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
+         END DO
+      END DO
+
       !------------------------------------------------------------------------------!
       ! If snow temperature is above freezing point, then snow melts 
       ! (should not happen but sometimes it does)
       !------------------------------------------------------------------------------!
-      DO ji = 1, npti
-         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
-            ! Contribution to heat flux to the ocean [W.m-2], < 0  
-            hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,1) * h_s_1d(ji) * a_i_1d(ji) * r1_rdtice
-            ! Contribution to mass flux
-            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * h_s_1d(ji) * a_i_1d(ji) * r1_rdtice
-            dh_s_mlt(ji)       = dh_s_mlt(ji) - h_s_1d(ji)
-            ! updates
-            h_s_1d(ji)    = 0._wp
-            e_s_1d (ji,1) = 0._wp
-            t_s_1d (ji,1) = rt0
-         END IF
-      END DO
-
-      !------------------------------------------------------------!
-      !  2) Computing layer thicknesses and enthalpies.            !
-      !------------------------------------------------------------!
-      DO jk = 1, nlay_i
-         DO ji = 1, npti
-            zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
-            zeh_i(ji)   = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk)
-         END DO
-      END DO
+      DO jk = 1, nlay_s
+         DO ji = 1, npti
+            IF( t_s_1d(ji,jk) > rt0 ) THEN !!! Internal melting
+               ! Contribution to heat flux to the ocean [W.m-2], < 0  
+               hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice
+               ! Contribution to mass flux
+               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice
+               dh_s_mlt(ji)       = dh_s_mlt(ji) - zh_s(ji,jk)
+               ! updates
+               h_s_1d(ji)    = h_s_1d(ji) - zh_s(ji,jk)
+               zh_s  (ji,jk) = 0._wp
+               e_s_1d(ji,jk) = 0._wp
+               t_s_1d(ji,jk) = rt0
+            END IF
+         END DO
+      END DO
+         
 
       !------------------------------------------------------------------------------|
@@ -244,4 +253,11 @@
       END DO
 
+      ! recalculate snow layers
+      DO jk = 1, nlay_s
+         DO ji = 1, npti
+            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
+         END DO
+      END DO
+
       ! If heat still available (zq_su > 0), then melt more snow
       zdeltah(1:npti,:) = 0._wp
@@ -250,8 +266,8 @@
          DO ji = 1, npti
             ! thickness change
-            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) 
+            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zh_s(ji,jk) ) ) 
             rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) ) 
             zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 )
-            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - h_s_1d(ji) ) ! bound melting
+            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting
             zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)    
             ! heat used to melt snow(W.m-2, >0)
@@ -263,4 +279,5 @@
             zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
             h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
+            zh_s(ji,jk)= MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
          END DO
       END DO
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 9270)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 9271)
@@ -474,4 +474,15 @@
          END DO
          !
+         DO jk = 1, nlay_s
+            DO jj = 1 , jpj
+               DO ji = 1 , jpi
+                  ! update exchanges with ocean
+                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
+                  t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
+               END DO
+            END DO
+         END DO
+         !
          DO jj = 1 , jpj
             DO ji = 1 , jpi
@@ -480,10 +491,4 @@
                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_rdtice
                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
-               hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0
-               !-----------------------------------------------------------------
-               ! Zap snow energy 
-               !-----------------------------------------------------------------
-               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
-               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
                !
                !-----------------------------------------------------------------
