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 9123)
+++ 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 )
