Index: /NEMO/trunk/src/ICE/ice.F90
===================================================================
--- /NEMO/trunk/src/ICE/ice.F90	(revision 10414)
+++ /NEMO/trunk/src/ICE/ice.F90	(revision 10415)
@@ -332,10 +332,10 @@
    !! * Old values of global variables
    !!----------------------------------------------------------------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b    !: snow and ice volumes/thickness
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b         !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                         !: snow heat content
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                         !: ice temperatures
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b              !: ice velocity
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                        !: ice concentration (total)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b, h_ip_b    !: snow and ice volumes/thickness
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b                 !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                                 !: snow heat content
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                                 !: ice temperatures
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b                      !: ice velocity
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                                !: ice concentration (total)
             
    !!----------------------------------------------------------------------
@@ -441,6 +441,6 @@
       ! * Old values of global variables
       ii = ii + 1
-      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl)        ,   &
-         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,   &
+      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl), h_ip_b(jpi,jpj,jpl),  &
+         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,               &
          &      oa_i_b(jpi,jpj,jpl)                                                   , STAT=ierr(ii) )
 
Index: /NEMO/trunk/src/ICE/icectl.F90
===================================================================
--- /NEMO/trunk/src/ICE/icectl.F90	(revision 10414)
+++ /NEMO/trunk/src/ICE/icectl.F90	(revision 10415)
@@ -190,7 +190,5 @@
 
       ! heat flux
-      zhfx  = glob_sum( ( qt_atm_oi - qt_oce_ai - 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 ) * zconv
+      zhfx  = glob_sum( ( qt_atm_oi - qt_oce_ai - diag_heat ) * e1e2t ) * zconv
 
       ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
Index: /NEMO/trunk/src/ICE/icedyn.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn.F90	(revision 10414)
+++ /NEMO/trunk/src/ICE/icedyn.F90	(revision 10415)
@@ -76,5 +76,5 @@
       INTEGER  ::   ji, jj, jl        ! dummy loop indices
       REAL(wp) ::   zcoefu, zcoefv
-      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max
+      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
       !!--------------------------------------------------------------------
@@ -101,16 +101,20 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               zhi_max(ji,jj,jl) = MAX( epsi20, h_i_b(ji,jj,jl), h_i_b(ji+1,jj  ,jl), h_i_b(ji  ,jj+1,jl), &
-                  &                                              h_i_b(ji-1,jj  ,jl), h_i_b(ji  ,jj-1,jl), &
-                  &                                              h_i_b(ji+1,jj+1,jl), h_i_b(ji-1,jj-1,jl), &
-                  &                                              h_i_b(ji+1,jj-1,jl), h_i_b(ji-1,jj+1,jl) )
-               zhs_max(ji,jj,jl) = MAX( epsi20, h_s_b(ji,jj,jl), h_s_b(ji+1,jj  ,jl), h_s_b(ji  ,jj+1,jl), &
-                  &                                              h_s_b(ji-1,jj  ,jl), h_s_b(ji  ,jj-1,jl), &
-                  &                                              h_s_b(ji+1,jj+1,jl), h_s_b(ji-1,jj-1,jl), &
-                  &                                              h_s_b(ji+1,jj-1,jl), h_s_b(ji-1,jj+1,jl) )
+               zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), &
+                  &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), &
+                  &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), &
+                  &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) )
+               zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), &
+                  &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), &
+                  &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), &
+                  &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) )
+               zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), &
+                  &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), &
+                  &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), &
+                  &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) )
             END DO
          END DO
       END DO
-      CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1. )
+      CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )
       !
       !
@@ -118,13 +122,13 @@
 
       CASE ( np_dynALL )           !==  all dynamical processes  ==!
-         CALL ice_dyn_rhg   ( kt )                                       ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness
-         CALL ice_dyn_rdgrft( kt )                                       ! -- ridging/rafting 
-         CALL ice_cor       ( kt , 1 )                                   ! -- Corrections
+         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
+         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
+         CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting 
+         CALL ice_cor       ( kt , 1 )                                             ! -- Corrections
 
       CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
-         CALL ice_dyn_rhg   ( kt )                                       ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness
-         CALL Hpiling                                                    ! -- simple pile-up (replaces ridging/rafting)
+         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
+         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
+         CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting)
 
       CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
@@ -163,5 +167,5 @@
          !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
          ! ---
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness
+         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
 
          ! diagnostics: divergence at T points 
@@ -186,5 +190,5 @@
 
 
-   SUBROUTINE Hbig( phi_max, phs_max )
+   SUBROUTINE Hbig( phi_max, phs_max, phip_max )
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE Hbig  ***
@@ -203,8 +207,8 @@
       !! ** input   : Max thickness of the surrounding 9-points
       !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max   ! max ice thick from surrounding 9-pts
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
       !
       INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zhi, zhs, zvs_excess, zfra
+      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra
       !!-------------------------------------------------------------------
       !
@@ -216,10 +220,16 @@
                IF ( v_i(ji,jj,jl) > 0._wp ) THEN
                   !
+                  !                               ! -- check h_ip -- !
+                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
+                  IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN
+                     zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )
+                     IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN
+                        a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl)
+                     ENDIF
+                  ENDIF
+                  !
                   !                               ! -- check h_i -- !
                   ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
                   zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-!!clem                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)  
-!!clem                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &
-!!clem                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN
                   IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
                      a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
@@ -233,8 +243,8 @@
                      !
                      wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
                      !
-                     e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra
-                     v_s(ji,jj,jl)   = a_i(ji,jj,jl) * phs_max(ji,jj,jl)
+                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
+                     v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl)
                   ENDIF           
                   !
@@ -242,13 +252,13 @@
                   ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
                   !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
-                  !    this imposed mini can artificially make the snow thickness very high (if concentration decreases drastically)
+                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
                   zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
                   IF( zvs_excess > 0._wp ) THEN
                      zfra = zvs_excess / MAX( v_s(ji,jj,jl), epsi20 )
                      wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
                      !
-                     e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra
-                     v_s(ji,jj,jl)   = v_s(ji,jj,jl) - zvs_excess
+                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
+                     v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess
                   ENDIF
                   
Index: /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 10414)
+++ /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 10415)
@@ -147,7 +147,7 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
-      REAL(wp), DIMENSION(jpi,jpj) ::   zssh_lead_m                     ! array used for the calculation of ice surface slope:
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
       !                                                                 !    ocean surface (ssh_m) if ice is not embedded
-      !                                                                 !    ice top surface if ice is embedded   
+      !                                                                 !    ice bottom surface if ice is embedded   
       REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
@@ -268,8 +268,8 @@
       ! 2) Wind / ocean stress, mass terms, coriolis terms
       !------------------------------------------------------------------------------!
-
-      !== embedded sea ice: compute representative ice top surface      ==!
-      !== non-embedded sea ice: use ocean surface for slope calculation ==!
-      zssh_lead_m(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
+      ! sea surface height
+      !    embedded sea ice: compute representative ice top surface
+      !    non-embedded sea ice: use ocean surface for slope calculation
+      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
 
       DO jj = 2, jpjm1
@@ -309,6 +309,6 @@
 
             ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
-            zspgU(ji,jj)    = - zmassU * grav * ( zssh_lead_m(ji+1,jj) - zssh_lead_m(ji,jj) ) * r1_e1u(ji,jj)
-            zspgV(ji,jj)    = - zmassV * grav * ( zssh_lead_m(ji,jj+1) - zssh_lead_m(ji,jj) ) * r1_e2v(ji,jj)
+            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
+            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
 
             ! masks
Index: /NEMO/trunk/src/ICE/icestp.F90
===================================================================
--- /NEMO/trunk/src/ICE/icestp.F90	(revision 10414)
+++ /NEMO/trunk/src/ICE/icestp.F90	(revision 10415)
@@ -194,13 +194,6 @@
                                         CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling)
          !
-!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
-!!       moreover it should only be called at the update frequency (as in agrif_ice_update.F90)
-!# if defined key_agrif
-!         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid()
-!# endif
                                         CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
-!# if defined key_agrif
-!         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid()
-!# endif
+         !
          IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs 
          !
@@ -368,9 +361,15 @@
       e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy
       WHERE( a_i_b(:,:,:) >= epsi20 )
-         h_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:)   ! ice thickness
-         h_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:)   ! snw thickness
+         h_i_b(:,:,:) = v_i_b(:,:,:) / a_i_b(:,:,:)   ! ice thickness
+         h_s_b(:,:,:) = v_s_b(:,:,:) / a_i_b(:,:,:)   ! snw thickness
       ELSEWHERE
          h_i_b(:,:,:) = 0._wp
          h_s_b(:,:,:) = 0._wp
+      END WHERE
+      
+      WHERE( a_ip(:,:,:) >= epsi20 )
+         h_ip_b(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)   ! ice pond thickness
+      ELSEWHERE
+         h_ip_b(:,:,:) = 0._wp
       END WHERE
       !
Index: /NEMO/trunk/src/ICE/icevar.F90
===================================================================
--- /NEMO/trunk/src/ICE/icevar.F90	(revision 10414)
+++ /NEMO/trunk/src/ICE/icevar.F90	(revision 10415)
@@ -954,5 +954,5 @@
    FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
       !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE rhg_evp_rst  ***
+      !!                   ***  ROUTINE ice_var_sshdyn  ***
       !!                     
       !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
@@ -963,5 +963,4 @@
       !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*, 
       !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
-      !!
       !!----------------------------------------------------------------------
       !
Index: /NEMO/trunk/src/OCE/stpctl.F90
===================================================================
--- /NEMO/trunk/src/OCE/stpctl.F90	(revision 10414)
+++ /NEMO/trunk/src/OCE/stpctl.F90	(revision 10415)
@@ -50,5 +50,5 @@
       !!              - Print it each 50 time steps
       !!              - Stop the run IF problem encountered by setting indic=-3
-      !!                Problems checked: |ssh| maximum larger than 10 m
+      !!                Problems checked: |ssh| maximum larger than 20 m
       !!                                  |U|   maximum larger than 10 m/s 
       !!                                  negative sea surface salinity
@@ -136,5 +136,5 @@
       ENDIF
       !
-      IF (  zmax(1) >   15._wp .OR.   &                    ! too large sea surface height ( > 15 m )
+      IF (  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m )
          &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s)
          &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity
@@ -159,5 +159,5 @@
          IF(lwp) THEN
             WRITE(numout,cform_err)
-            WRITE(numout,*) ' stp_ctl: |ssh| > 10 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
+            WRITE(numout,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
             WRITE(numout,*) ' ======= '
             WRITE(numout,9100) kt,   zmax(1), iih , ijh
