Index: NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icedyn_adv_pra.F90	(revision 13360)
+++ NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icedyn_adv_pra.F90	(revision 14037)
@@ -44,4 +44,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvl , syvl , sxxvl , syyvl , sxyvl    ! melt pond lid volume
 
    !! * Substitutions
@@ -55,5 +56,5 @@
 
    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 )
+      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                **  routine ice_dyn_adv_pra  **
@@ -81,41 +82,55 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid thickness
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
       !
-      INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
       INTEGER  ::   icycle                  ! number of sub-timestep for the advection
-      REAL(wp) ::   zdt                     !   -      -
+      REAL(wp) ::   zdt, z1_dt              !   -      -
       REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication
       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)        ::   zhi_max, zhs_max, zhip_max, zs_i, zsi_max
+      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   ze_i, zei_max
+      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   ze_s, zes_max
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
-      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0vl
       REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es
       REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei
+      !! diagnostics
+      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat      
       !!----------------------------------------------------------------------
       !
       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_2D( 0, 0, 0, 0 )
-            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_2D
+      ! --- Record max of the surrounding 9-pts (for call Hbig) --- !
+      ! thickness and salinity
+      WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:)
+      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp
+      END WHERE
+      CALL icemax3D( ph_i , zhi_max )
+      CALL icemax3D( ph_s , zhs_max )
+      CALL icemax3D( ph_ip, zhip_max)
+      CALL icemax3D( zs_i , zsi_max )
+      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp )
+      !
+      ! enthalpies
+      DO jk = 1, nlay_i
+         WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:)
+         ELSEWHERE                      ; ze_i(:,:,jk,:) = 0._wp
+         END WHERE
       END DO
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp )
+      DO jk = 1, nlay_s
+         WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:)
+         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp
+         END WHERE
+      END DO   
+      CALL icemax4D( ze_i , zei_max )
+      CALL icemax4D( ze_s , zes_max )
+      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp )
+      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp )
+      !
       !
       ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
@@ -132,4 +147,5 @@
       ENDIF
       zdt = rDt_ice / REAL(icycle)
+      z1_dt = 1._wp / zdt
       
       ! --- transport --- !
@@ -138,4 +154,11 @@
 
       DO jt = 1, icycle
+
+         ! diagnostics
+         zdiag_adv_mass(:,:) =   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
+            &                  + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow
+         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi
+         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
+            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 )
 
          ! record at_i before advection (for open water)
@@ -156,7 +179,10 @@
                z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
             END DO
-            IF ( ln_pnd_H12 ) THEN
-               z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
-               z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
+            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
+               z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction
+               z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume
+               IF ( ln_pnd_lids ) THEN
+                  z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:)   ! Melt pond lid volume
+               ENDIF
             ENDIF
          END DO
@@ -189,9 +215,13 @@
             END DO
             !
-            IF ( ln_pnd_H12 ) THEN
+            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
                CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
                CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
                CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
+               IF ( ln_pnd_lids ) THEN
+                  CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
+                  CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 
+               ENDIF
             ENDIF
             !                                                               !--------------------------------------------!
@@ -220,11 +250,42 @@
                   &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
             END DO
-            IF ( ln_pnd_H12 ) THEN
+            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
                CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
                CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
                CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
+               IF ( ln_pnd_lids ) THEN
+                  CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
+                  CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 
+               ENDIF
             ENDIF
             !
+         ENDIF
+         
+         ! --- Lateral boundary conditions --- !
+         !     caution: for gradients (sx and sy) the sign changes
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp  & ! ice volume
+            &                                , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  &
+            &                                , z0snw , 'T', 1._wp, sxsn  , 'T', -1._wp, sysn  , 'T', -1._wp  & ! snw volume
+            &                                , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  )
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity
+            &                                , sxxsal, 'T', 1._wp, syysal, 'T',  1._wp, sxysal, 'T',  1._wp  &
+            &                                , z0ai  , 'T', 1._wp, sxa   , 'T', -1._wp, sya   , 'T', -1._wp  & ! ice concentration
+            &                                , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  )
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age
+            &                                , sxxage, 'T', 1._wp, syyage, 'T',  1._wp, sxyage, 'T',  1._wp  )
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es  , 'T', 1._wp, sxc0  , 'T', -1._wp, syc0  , 'T', -1._wp  & ! snw enthalpy
+            &                                , sxxc0 , 'T', 1._wp, syyc0 , 'T',  1._wp, sxyc0 , 'T',  1._wp  ) 
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy
+            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  )
+         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
+            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
+               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
+               &                                , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume
+               &                                , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  ) 
+            IF ( ln_pnd_lids ) THEN
+               CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume
+                  &                                , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  ) 
+            ENDIF
          ENDIF
 
@@ -242,7 +303,10 @@
                pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
             END DO
-            IF ( ln_pnd_H12 ) THEN
+            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
                pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+               IF ( ln_pnd_lids ) THEN
+                  pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+               ENDIF
             ENDIF
          END DO
@@ -256,12 +320,23 @@
          CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1.0_wp )
          !
+         ! --- diagnostics --- !
+         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
+            &                                        + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow &
+            &                                        - zdiag_adv_mass(:,:) ) * z1_dt
+         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi &
+            &                                        - zdiag_adv_salt(:,:) ) * z1_dt
+         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
+            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) &
+            &                                        - zdiag_adv_heat(:,:) ) * z1_dt
+         !
          ! --- Ensure non-negative fields --- !
          !     Remove negative values (conservation is ensured)
          !     (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 )
+         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, 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 )
+         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, &
+            &            pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
          !
          ! --- Ensure snow load is not too big --- !
@@ -292,11 +367,18 @@
       !! 
       INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
+      INTEGER  ::   jj0                                  ! dummy loop indices
       REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars
       REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
+      REAL(wp) ::   zpsm, zps0
+      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy
       REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
       REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
       !-----------------------------------------------------------------------
+      ! in order to avoid lbc_lnk (communications):
+      !    jj loop must be 1:jpj   if adv_x is called first
+      !                and 2:jpj-1 if adv_x is called second
+      jj0 = NINT(pcrh)
       !
       jcat = SIZE( ps0 , 3 )   ! size of input arrays
@@ -305,128 +387,159 @@
          !
          ! Limitation of moments.                                           
-         DO_2D( 0, 0, 1, 1 )
-            !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
-            psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
-            !
-            zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
-            zs1max  = 1.5 * zslpmax
-            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) )
-            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
-               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  )
-            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
-
-            ps0 (ji,jj,jl) = zslpmax  
-            psx (ji,jj,jl) = zs1new         * rswitch
-            psxx(ji,jj,jl) = zs2new         * rswitch
-            psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch
-            psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch
-            psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
-         END_2D
-
-         !  Calculate fluxes and moments between boxes i<-->i+1              
-         DO_2D( 0, 0, 1, 1 )
-            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
-            zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
-            zalfq        =  zalf * zalf
-            zalf1        =  1.0 - zalf
-            zalf1q       =  zalf1 * zalf1
-            !
-            zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl)
-            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) )
-            zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) )
-            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq
-            zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
-            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl)
-            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl)
-
-            !  Readjust moments remaining in the box.
-            psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
-            ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
-            psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) )
-            psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl)
-            psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj)
-            psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj)
-            psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
-         END_2D
-
-         DO_2D( 0, 0, 1, 0 )
-            zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
-            zalg  (ji,jj) = zalf
-            zalfq         = zalf * zalf
-            zalf1         = 1.0 - zalf
-            zalg1 (ji,jj) = zalf1
-            zalf1q        = zalf1 * zalf1
-            zalg1q(ji,jj) = zalf1q
-            !
-            zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
-            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
-               &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
-            zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
-            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
-            zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
-            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
-            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
-         END_2D
-
-         DO_2D( 0, 0, 0, 0 )
-            zbt  =       zbet(ji-1,jj)
-            zbt1 = 1.0 - zbet(ji-1,jj)
-            !
-            psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) )
-            ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) )
-            psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) )
-            psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl)
-            psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) )
-            psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) )
-            psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl)
-         END_2D
-
-         !   Put the temporary moments into appropriate neighboring boxes.    
-         DO_2D( 0, 0, 0, 0 )
-            zbt  =       zbet(ji-1,jj)
-            zbt1 = 1.0 - zbet(ji-1,jj)
-            psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)
-            zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)
-            zalf1         = 1.0 - zalf
-            ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj)
-            !
-            ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)
-            psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)
-            psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             &
-               &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) &
-               &            + zbt1 * psxx(ji,jj,jl)
-            psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             &
-               &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   &
-               &            + zbt1 * psxy(ji,jj,jl)
-            psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)
-            psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)
-         END_2D
-
-         DO_2D( 0, 0, 0, 0 )
-            zbt  =       zbet(ji,jj)
-            zbt1 = 1.0 - zbet(ji,jj)
-            psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
-            zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
-            zalf1         = 1.0 - zalf
-            ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
-            !
-            ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
-            psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp )
-            psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
-               &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    &
-               &                                           + ( zalf1 - zalf ) * ztemp ) )
-            psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
-               &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) )
-            psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) )
-            psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) )
-         END_2D
-
+         DO jj = Njs0 - jj0, Nje0 + jj0
+            
+            DO ji = Nis0 - 1, Nie0 + 1
+
+               zpsm  = psm (ji,jj,jl) ! optimization
+               zps0  = ps0 (ji,jj,jl)
+               zpsx  = psx (ji,jj,jl)
+               zpsxx = psxx(ji,jj,jl)
+               zpsy  = psy (ji,jj,jl)
+               zpsyy = psyy(ji,jj,jl)
+               zpsxy = psxy(ji,jj,jl)
+
+               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
+               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
+               !
+               zslpmax = MAX( 0._wp, zps0 )
+               zs1max  = 1.5 * zslpmax
+               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) )
+               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
+               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+
+               zps0  = zslpmax  
+               zpsx  = zs1new  * rswitch
+               zpsxx = zs2new  * rswitch
+               zpsy  = zpsy    * rswitch
+               zpsyy = zpsyy   * rswitch
+               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
+
+               !  Calculate fluxes and moments between boxes i<-->i+1              
+               !                                !  Flux from i to i+1 WHEN u GT 0 
+               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
+               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
+               zalfq        =  zalf * zalf
+               zalf1        =  1.0 - zalf
+               zalf1q       =  zalf1 * zalf1
+               !
+               zfm (ji,jj)  =  zalf  *   zpsm 
+               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
+               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx )
+               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq
+               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy )
+               zfxy(ji,jj)  =  zalfq *   zpsxy
+               zfyy(ji,jj)  =  zalf  *   zpsyy
+
+               !                                !  Readjust moments remaining in the box.
+               zpsm  =  zpsm  - zfm(ji,jj)
+               zps0  =  zps0  - zf0(ji,jj)
+               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
+               zpsxx =  zalf1  * zalf1q * zpsxx
+               zpsy  =  zpsy  - zfy (ji,jj)
+               zpsyy =  zpsyy - zfyy(ji,jj)
+               zpsxy =  zalf1q * zpsxy
+               !
+               psm (ji,jj,jl) = zpsm ! optimization
+               ps0 (ji,jj,jl) = zps0 
+               psx (ji,jj,jl) = zpsx 
+               psxx(ji,jj,jl) = zpsxx
+               psy (ji,jj,jl) = zpsy 
+               psyy(ji,jj,jl) = zpsyy
+               psxy(ji,jj,jl) = zpsxy
+               !
+            END DO
+            
+            DO ji = Nis0 - 1, Nie0
+               !                                !  Flux from i+1 to i when u LT 0.
+               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
+               zalg  (ji,jj) = zalf
+               zalfq         = zalf * zalf
+               zalf1         = 1.0 - zalf
+               zalg1 (ji,jj) = zalf1
+               zalf1q        = zalf1 * zalf1
+               zalg1q(ji,jj) = zalf1q
+               !
+               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
+               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
+                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
+               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
+               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
+               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
+               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
+               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
+            END DO
+
+            DO ji = Nis0, Nie0
+               !
+               zpsm  = psm (ji,jj,jl) ! optimization
+               zps0  = ps0 (ji,jj,jl)
+               zpsx  = psx (ji,jj,jl)
+               zpsxx = psxx(ji,jj,jl)
+               zpsy  = psy (ji,jj,jl)
+               zpsyy = psyy(ji,jj,jl)
+               zpsxy = psxy(ji,jj,jl)
+               !                                !  Readjust moments remaining in the box.
+               zbt  =       zbet(ji-1,jj)
+               zbt1 = 1.0 - zbet(ji-1,jj)
+               !
+               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
+               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
+               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
+               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
+               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  - zfy (ji-1,jj) )
+               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
+               zpsxy = zalg1q(ji-1,jj) * zpsxy
+
+               !   Put the temporary moments into appropriate neighboring boxes.    
+               !                                !   Flux from i to i+1 IF u GT 0.
+               zbt   =       zbet(ji-1,jj)
+               zbt1  = 1.0 - zbet(ji-1,jj)
+               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm
+               zalf  = zbt * zfm(ji-1,jj) / zpsm
+               zalf1 = 1.0 - zalf
+               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
+               !
+               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0
+               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx
+               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            &
+                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &
+                  &            + zbt1 * zpsxx
+               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            &
+                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  &
+                  &            + zbt1 * zpsxy
+               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy 
+               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy
+
+               !                                !  Flux from i+1 to i IF u LT 0.
+               zbt   =       zbet(ji,jj)
+               zbt1  = 1.0 - zbet(ji,jj)
+               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
+               zalf  = zbt1 * zfm(ji,jj) / zpsm
+               zalf1 = 1.0 - zalf
+               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
+               !
+               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) )
+               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp )
+               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx &
+                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    &
+                  &                         + ( zalf1 - zalf ) * ztemp ) )
+               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  &
+                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) )
+               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) )
+               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) )
+               !
+               psm (ji,jj,jl) = zpsm  ! optimization
+               ps0 (ji,jj,jl) = zps0 
+               psx (ji,jj,jl) = zpsx 
+               psxx(ji,jj,jl) = zpsxx
+               psy (ji,jj,jl) = zpsy 
+               psyy(ji,jj,jl) = zpsyy
+               psxy(ji,jj,jl) = zpsxy
+            END DO
+            !
+         END DO
+         !
       END DO
-
-      !-- Lateral boundary conditions
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   &
-         &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes
-         &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp )
-      !
+      !      
    END SUBROUTINE adv_x
 
@@ -449,11 +562,18 @@
       !!
       INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
+      INTEGER  ::   ji0                                  ! dummy loop indices
       REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars
       REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
+      REAL(wp) ::   zpsm, zps0
+      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy
       REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
       REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
       !---------------------------------------------------------------------
+      ! in order to avoid lbc_lnk (communications):
+      !    ji loop must be 1:jpi   if adv_y is called first
+      !                and 2:jpi-1 if adv_y is called second
+      ji0 = NINT(pcrh)
       !
       jcat = SIZE( ps0 , 3 )   ! size of input arrays
@@ -462,50 +582,66 @@
          !
          ! Limitation of moments.
-         DO_2D( 1, 1, 0, 0 )
-            !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
-            psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  )
-            !
-            zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
+         DO_2D( 1, 1, ji0, ji0 )
+            !
+            zpsm  = psm (ji,jj,jl) ! optimization
+            zps0  = ps0 (ji,jj,jl)
+            zpsx  = psx (ji,jj,jl)
+            zpsxx = psxx(ji,jj,jl)
+            zpsy  = psy (ji,jj,jl)
+            zpsyy = psyy(ji,jj,jl)
+            zpsxy = psxy(ji,jj,jl)
+            !
+            !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise)
+            zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  )
+            !
+            zslpmax = MAX( 0._wp, zps0 )
             zs1max  = 1.5 * zslpmax
-            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) )
-            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
-               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  )
+            zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) )
+            zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) )
             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
             !
-            ps0 (ji,jj,jl) = zslpmax  
-            psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch
-            psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch
-            psy (ji,jj,jl) = zs1new         * rswitch
-            psyy(ji,jj,jl) = zs2new         * rswitch
-            psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
-         END_2D
- 
-         !  Calculate fluxes and moments between boxes j<-->j+1              
-         DO_2D( 1, 1, 0, 0 )
+            zps0  = zslpmax  
+            zpsx  = zpsx  * rswitch
+            zpsxx = zpsxx * rswitch
+            zpsy  = zs1new         * rswitch
+            zpsyy = zs2new         * rswitch
+            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
+
+            !  Calculate fluxes and moments between boxes j<-->j+1              
+            !                                !  Flux from j to j+1 WHEN v GT 0   
             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
-            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
+            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm
             zalfq        =  zalf * zalf
             zalf1        =  1.0 - zalf
             zalf1q       =  zalf1 * zalf1
             !
-            zfm (ji,jj)  =  zalf  * psm(ji,jj,jl)
-            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 
-            zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) )
-            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl)
-            zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
-            zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl)
-            zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl)
-            !
-            !  Readjust moments remaining in the box.
-            psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
-            ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
-            psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) )
-            psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl)
-            psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj)
-            psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj)
-            psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
+            zfm (ji,jj)  =  zalf  * zpsm
+            zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) ) 
+            zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy )
+            zfyy(ji,jj)  =  zalf  * zalfq * zpsyy
+            zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy )
+            zfxy(ji,jj)  =  zalfq * zpsxy
+            zfxx(ji,jj)  =  zalf  * zpsxx
+            !
+            !                                !  Readjust moments remaining in the box.
+            zpsm   =  zpsm  - zfm(ji,jj)
+            zps0   =  zps0  - zf0(ji,jj)
+            zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy )
+            zpsyy  =  zalf1 * zalf1q * zpsyy
+            zpsx   =  zpsx  - zfx(ji,jj)
+            zpsxx  =  zpsxx - zfxx(ji,jj)
+            zpsxy  =  zalf1q * zpsxy
+            !
+            psm (ji,jj,jl) = zpsm ! optimization
+            ps0 (ji,jj,jl) = zps0 
+            psx (ji,jj,jl) = zpsx 
+            psxx(ji,jj,jl) = zpsxx
+            psy (ji,jj,jl) = zpsy 
+            psyy(ji,jj,jl) = zpsyy
+            psxy(ji,jj,jl) = zpsxy
          END_2D
          !
-         DO_2D( 1, 0, 0, 0 )
+         DO_2D( 1, 0, ji0, ji0 )
+            !                                !  Flux from j+1 to j when v LT 0.
             zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 
             zalg  (ji,jj) = zalf
@@ -526,70 +662,80 @@
          END_2D
 
-         !  Readjust moments remaining in the box. 
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, ji0, ji0 )
+            !                                !  Readjust moments remaining in the box.
             zbt  =         zbet(ji,jj-1)
             zbt1 = ( 1.0 - zbet(ji,jj-1) )
             !
-            psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) )
-            ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) )
-            psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) )
-            psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl)
-            psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) )
-            psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) )
-            psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl)
+            zpsm  = psm (ji,jj,jl) ! optimization
+            zps0  = ps0 (ji,jj,jl)
+            zpsx  = psx (ji,jj,jl)
+            zpsxx = psxx(ji,jj,jl)
+            zpsy  = psy (ji,jj,jl)
+            zpsyy = psyy(ji,jj,jl)
+            zpsxy = psxy(ji,jj,jl)
+            !
+            zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
+            zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
+            zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy )
+            zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy
+            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) )
+            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
+            zpsxy = zalg1q(ji,jj-1) * zpsxy
+
+            !   Put the temporary moments into appropriate neighboring boxes.    
+            !                                !   Flux from j to j+1 IF v GT 0.
+            zbt   =       zbet(ji,jj-1)
+            zbt1  = 1.0 - zbet(ji,jj-1)
+            zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm 
+            zalf  = zbt * zfm(ji,jj-1) / zpsm 
+            zalf1 = 1.0 - zalf
+            ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
+            !
+            zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0
+            zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  &
+               &             + zbt1 * zpsy  
+            zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           &
+               &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
+               &             + zbt1 * zpsyy
+            zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             &
+               &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  &
+               &             + zbt1 * zpsxy
+            zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx 
+            zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx
+
+            !                                !  Flux from j+1 to j IF v LT 0.
+            zbt   =       zbet(ji,jj)
+            zbt1  = 1.0 - zbet(ji,jj)
+            zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
+            zalf  = zbt1 * zfm(ji,jj) / zpsm
+            zalf1 = 1.0 - zalf
+            ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
+            !
+            zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) )
+            zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp )
+            zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy &
+               &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     &
+               &                         + ( zalf1 - zalf ) * ztemp ) )
+            zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  &
+               &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) )
+            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) )
+            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) )
+            !
+            psm (ji,jj,jl) = zpsm ! optimization
+            ps0 (ji,jj,jl) = zps0 
+            psx (ji,jj,jl) = zpsx 
+            psxx(ji,jj,jl) = zpsxx
+            psy (ji,jj,jl) = zpsy 
+            psyy(ji,jj,jl) = zpsyy
+            psxy(ji,jj,jl) = zpsxy
          END_2D
-
-         !   Put the temporary moments into appropriate neighboring boxes.    
-         DO_2D( 0, 0, 0, 0 )
-            zbt  =       zbet(ji,jj-1)
-            zbt1 = 1.0 - zbet(ji,jj-1)
-            psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 
-            zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 
-            zalf1         = 1.0 - zalf
-            ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1)
-            !
-            ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)
-            psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  &
-               &             + zbt1 * psy(ji,jj,jl)  
-            psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           &
-               &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
-               &             + zbt1 * psyy(ji,jj,jl)
-            psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            &
-               &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  &
-               &             + zbt1 * psxy(ji,jj,jl)
-            psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)
-            psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)
-         END_2D
-
-         DO_2D( 0, 0, 0, 0 )
-            zbt  =       zbet(ji,jj)
-            zbt1 = 1.0 - zbet(ji,jj)
-            psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
-            zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
-            zalf1         = 1.0 - zalf
-            ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
-            !
-            ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) )
-            psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )
-            psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) &
-               &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    &
-               &                                            + ( zalf1 - zalf ) * ztemp ) )
-            psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
-               &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) )
-            psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) )
-            psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) )
-         END_2D
-
+         !
       END DO
-
-      !-- Lateral boundary conditions
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   &
-         &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes
-         &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp )
       !
    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 )
+   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, &
+      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE Hbig  ***
@@ -605,11 +751,14 @@
       !! ** 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)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
-      !
-      INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i
+      !
+      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
+      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra
       !!-------------------------------------------------------------------
       !
@@ -617,5 +766,4 @@
       !
       DO jl = 1, jpl
-
          DO_2D( 1, 1, 1, 1 )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
@@ -623,5 +771,5 @@
                !                               ! -- 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
+               IF( ln_pnd_LEV .OR. ln_pnd_TOPO .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
@@ -650,7 +798,45 @@
                ENDIF           
                !                  
+               !                               ! -- check s_i -- !
+               ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean
+               zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl)
+               IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                  zfra = psi_max(ji,jj,jl) / zsi
+                  sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt
+                  psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra
+               ENDIF
+               !
             ENDIF
          END_2D
       END DO 
+      !
+      !                                           ! -- check e_i/v_i -- !
+      DO jl = 1, jpl
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
+            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+               ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
+               zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl)
+               IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                  zfra = pei_max(ji,jj,jk,jl) / zei
+                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
+                  pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra
+               ENDIF
+            ENDIF
+         END_3D
+      END DO
+      !                                           ! -- check e_s/v_s -- !
+      DO jl = 1, jpl
+         DO_3D( 1, 1, 1, 1, 1, nlay_s )
+            IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
+               ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
+               zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl)
+               IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                  zfra = pes_max(ji,jj,jk,jl) / zes
+                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
+                  pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra
+               ENDIF
+            ENDIF
+         END_3D
+      END DO
       !
    END SUBROUTINE Hbig
@@ -724,6 +910,7 @@
          &      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) ,   &
+         &      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) ,   &
+         &      sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) ,   &
          !
          &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
@@ -772,30 +959,30 @@
             !
             !                                                        ! ice thickness
-            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice  )
-            CALL iom_get( numrir, jpdom_auto, 'syice' , syice  )
+            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice )
             CALL iom_get( numrir, jpdom_auto, 'syyice', syyice )
             CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice )
             !                                                        ! snow thickness
-            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn   )
-            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn   )
+            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn  )
             CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn  )
             CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn  )
             !                                                        ! ice concentration
-            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa    )
-            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya    )
+            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa   , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya   , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxa'  , sxxa   )
             CALL iom_get( numrir, jpdom_auto, 'syya'  , syya   )
             CALL iom_get( numrir, jpdom_auto, 'sxya'  , sxya   )
             !                                                        ! ice salinity
-            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal  )
-            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal  )
+            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal )
             CALL iom_get( numrir, jpdom_auto, 'syysal', syysal )
             CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal )
             !                                                        ! ice age
-            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage  )
-            CALL iom_get( numrir, jpdom_auto, 'syage' , syage  )
+            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage )
             CALL iom_get( numrir, jpdom_auto, 'syyage', syyage )
@@ -804,32 +991,59 @@
             DO jk = 1, nlay_s
                WRITE(zchar1,'(I2.2)') jk
-               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
-               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
-               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
-               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
-               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxc0'//'_l'//zchar1  
+               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
+               znam = 'syc0'//'_l'//zchar1  
+               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxxc0'//'_l'//zchar1 
+               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
+               znam = 'syyc0'//'_l'//zchar1 
+               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxyc0'//'_l'//zchar1 
+               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
             END DO
             !                                                        ! ice layers heat content
             DO jk = 1, nlay_i
                WRITE(zchar1,'(I2.2)') jk
-               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
-               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:)
-               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
-               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
-               znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
-            END DO
-            !
-            IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction
-               CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap  )
-               CALL iom_get( numrir, jpdom_auto, 'syap' , syap  )
-               CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap )
-               CALL iom_get( numrir, jpdom_auto, 'syyap', syyap )
-               CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap )
-               !                                                     ! melt pond volume
-               CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp  )
-               CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp  )
-               CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp )
-               CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp )
-               CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp )
+               znam = 'sxe'//'_l'//zchar1   
+               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
+               znam = 'sye'//'_l'//zchar1   
+               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sye (:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxxe'//'_l'//zchar1  
+               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
+               znam = 'syye'//'_l'//zchar1  
+               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxye'//'_l'//zchar1  
+               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
+            END DO
+            !
+            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                    ! melt pond fraction
+               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
+                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp )
+                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp )
+                  CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap )
+                  CALL iom_get( numrir, jpdom_auto, 'syyap', syyap )
+                  CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap )
+                  !                                                     ! melt pond volume
+                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp )
+                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp )
+                  CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp )
+                  CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp )
+                  CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp )
+               ELSE
+                  sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp   ! melt pond fraction
+                  sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp   ! melt pond volume
+               ENDIF
+                  !
+               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume
+                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN
+                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp )
+                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp )
+                     CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl )
+                     CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl )
+                     CALL iom_get( numrir, jpdom_auto, 'sxyvl', sxyvl )
+                  ELSE
+                     sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp   ! melt pond lid volume
+                  ENDIF
+               ENDIF
             ENDIF
             !
@@ -845,7 +1059,10 @@
             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
-               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
-               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
+            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
+               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction
+               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume
+               IF ( ln_pnd_lids ) THEN
+                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume
+               ENDIF
             ENDIF
          ENDIF
@@ -862,53 +1079,63 @@
          !
          !                                                           ! ice thickness
-         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
-         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
-         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice)
+         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice)
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice)
+         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice)
+         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice)
          !                                                           ! snow thickness
-         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
-         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
-         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn )
+         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn )
+         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn )
          !                                                           ! ice concentration
-         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
-         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
-         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa  )
+         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya  )
          !                                                           ! ice salinity
-         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
-         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal)
+         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal)
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal)
+         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal)
+         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal)
          !                                                           ! ice age
-         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
-         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
-         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage)
+         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage)
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage)
+         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage)
+         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage)
          !                                                           ! 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 )
+            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
             WRITE(zchar1,'(I2.2)') jk
-            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
-            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
-            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
-            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
-            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z3d)
+            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z3d)
+            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z3d)
+            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z3d)
+            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z3d)
          END DO
          !
-         IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction
+         IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                       ! melt pond fraction
             CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
             CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
@@ -922,4 +1149,12 @@
             CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
             CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
+            !
+            IF ( ln_pnd_lids ) THEN                                  ! melt pond lid volume
+               CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  )
+               CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  )
+               CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl )
+               CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl )
+               CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl )
+            ENDIF
          ENDIF
          !
@@ -927,4 +1162,55 @@
       !
    END SUBROUTINE adv_pra_rst
+
+   SUBROUTINE icemax3D( pice , pmax )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE icemax3D ***                     
+      !! ** Purpose :  compute the max of the 9 points around
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output
+      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array
+      INTEGER  ::   ji, jj, jl   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      DO jl = 1, jpl
+         DO jj = Njs0-1, Nje0+1    
+            DO ji = Nis0, Nie0
+               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) )
+            END DO
+         END DO
+         DO jj = Njs0, Nje0    
+            DO ji = Nis0, Nie0
+               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
+            END DO
+         END DO
+      END DO
+   END SUBROUTINE icemax3D
+
+   SUBROUTINE icemax4D( pice , pmax )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE icemax4D ***                     
+      !! ** Purpose :  compute the max of the 9 points around
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input
+      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output
+      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array
+      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      jlay = SIZE( pice , 3 )   ! size of input arrays
+      DO jl = 1, jpl
+         DO jk = 1, jlay
+            DO jj = Njs0-1, Nje0+1    
+               DO ji = Nis0, Nie0
+                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) )
+               END DO
+            END DO
+            DO jj = Njs0, Nje0    
+               DO ji = Nis0, Nie0
+                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
+               END DO
+            END DO
+         END DO
+      END DO
+   END SUBROUTINE icemax4D
 
 #else
