Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv.F90	(revision 12236)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv.F90	(revision 12252)
@@ -88,5 +88,5 @@
       CASE( np_advPRA )                ! PRATHER scheme        !
          !                             !-----------------------!
-         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, &
+         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, h_i, h_s, h_ip, &
             &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
       END SELECT
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_pra.F90	(revision 12236)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_pra.F90	(revision 12252)
@@ -54,5 +54,5 @@
 CONTAINS
 
-   SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  &
+   SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  &
       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!----------------------------------------------------------------------
@@ -70,4 +70,7 @@
       REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
       REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness
       REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
@@ -87,4 +90,5 @@
       REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2
       REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
@@ -95,4 +99,25 @@
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
+      !
+      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
+      DO jl = 1, jpl
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), &
+                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), &
+                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
+                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
+               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), &
+                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), &
+                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
+                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
+               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), &
+                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), &
+                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
+                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
       !
       ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
@@ -239,4 +264,8 @@
          !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+         !
+         ! --- Make sure ice thickness is not too big --- !
+         !     (because ice thickness can be too large where ice concentration is very small)
+         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s )
          !
          ! --- Ensure snow load is not too big --- !
@@ -588,4 +617,72 @@
       !
    END SUBROUTINE adv_y
+
+
+   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE Hbig  ***
+      !!
+      !! ** Purpose : Thickness correction in case advection scheme creates
+      !!              abnormally tick ice or snow
+      !!
+      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
+      !!                 (before advection) and reduce it by adapting ice concentration
+      !!              2- check whether snow thickness is larger than the surrounding 9-points
+      !!                 (before advection) and reduce it by sending the excess in the ocean
+      !!
+      !! ** input   : Max thickness of the surrounding 9-points
+      !!-------------------------------------------------------------------
+      REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
+      !
+      INTEGER  ::   ji, jj, jl         ! dummy loop indices
+      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra
+      !!-------------------------------------------------------------------
+      !
+      z1_dt = 1._wp / pdt
+      !
+      DO jl = 1, jpl
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( pv_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. pv_ip(ji,jj,jl) > 0._wp ) THEN
+                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
+                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
+                        pa_ip(ji,jj,jl) = pv_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 = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
+                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                     pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
+                  ENDIF
+                  !
+                  !                               ! -- check h_s -- !
+                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
+                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
+                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
+                     !
+                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
+                     !
+                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
+                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
+                  ENDIF           
+                  !                  
+               ENDIF
+            END DO
+         END DO
+      END DO 
+      !
+   END SUBROUTINE Hbig
 
 
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90	(revision 12236)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90	(revision 12252)
@@ -352,8 +352,11 @@
          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
          !
-         ! Make sure ice thickness is not too big
-         !    (because ice thickness can be too large where ice concentration is very small)
-         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
-
+         ! --- Make sure ice thickness is not too big --- !
+         !     (because ice thickness can be too large where ice concentration is very small)
+         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s )
+         !
+         ! --- Ensure snow load is not too big --- !
+         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
+         !
       END DO
       !
@@ -1514,5 +1517,5 @@
 
 
-   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s )
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE Hbig  ***
@@ -1525,7 +1528,4 @@
       !!              2- check whether snow thickness is larger than the surrounding 9-points
       !!                 (before advection) and reduce it by sending the excess in the ocean
-      !!              3- check whether snow load deplets the snow-ice interface below sea level$
-      !!                 and reduce it by sending the excess in the ocean
-      !!              4- correct pond concentration to avoid a_ip > a_i
       !!
       !! ** input   : Max thickness of the surrounding 9-points
@@ -1533,11 +1533,9 @@
       REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step
       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
-      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i
-      !
-      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
-      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zvs_excess, zfra
-      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
+      !
+      INTEGER  ::   ji, jj, jl         ! dummy loop indices
+      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra
       !!-------------------------------------------------------------------
       !
@@ -1578,29 +1576,67 @@
                      pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
                   ENDIF           
+                  !                  
+               ENDIF
+            END DO
+         END DO
+      END DO 
+      !
+   END SUBROUTINE Hbig
+
+
+   SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE Hsnow  ***
+      !!
+      !! ** Purpose : 1- Check snow load after advection
+      !!              2- Correct pond concentration to avoid a_ip > a_i
+      !!
+      !! ** Method :  If snow load makes snow-ice interface to deplet below the ocean surface
+      !!              then put the snow excess in the ocean
+      !!
+      !! ** Notes :   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 very thick (if concentration decreases drastically)
+      !!              This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
+      !!-------------------------------------------------------------------
+      REAL(wp)                    , INTENT(in   ) ::   pdt   ! tracer time-step
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
+      !
+      INTEGER  ::   ji, jj, jl   ! dummy loop indices
+      REAL(wp) ::   z1_dt, zvs_excess, zfra
+      !!-------------------------------------------------------------------
+      !
+      z1_dt = 1._wp / pdt
+      !
+      ! -- check snow load -- !
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                   !
-                  !                               ! -- check snow load -- !
-                  ! 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 very thick (if concentration decreases drastically)
                   zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
-                  IF( zvs_excess > 0._wp ) THEN
+                  !
+                  IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
+                     ! put snow excess in the ocean
                      zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
                      wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
                      hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
-                     !
+                     ! correct snow volume and heat content
                      pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
                      pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
                   ENDIF
-                  
+                  !
                ENDIF
             END DO
          END DO
-      END DO 
-      !                                           !-- correct pond concentration to avoid a_ip > a_i
+      END DO
+      !
+      !-- correct pond concentration to avoid a_ip > a_i -- !
       WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
       !
-      !
-   END SUBROUTINE Hbig
-   
+   END SUBROUTINE Hsnow
+
+
 #else
    !!----------------------------------------------------------------------
