Index: NEMO/releases/release-4.0.1/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/releases/release-4.0.1/src/ICE/icedyn_adv_pra.F90	(revision 11697)
+++ NEMO/releases/release-4.0.1/src/ICE/icedyn_adv_pra.F90	(revision 11730)
@@ -39,5 +39,4 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow layers heat content
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content
@@ -81,10 +80,11 @@
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
       !
-      INTEGER  ::   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), 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)        ::   zarea
-      REAL(wp), DIMENSION(jpi,jpj,1)          ::   z0opw
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp
@@ -109,130 +109,138 @@
       zdt = rdt_ice / REAL(icycle)
       
-      !-------------------------
-      ! transported fields                                        
-      !-------------------------
-      z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:)              ! Open water area 
-      DO jl = 1, jpl
-         zarea(:,:,jl) = e1e2t(:,:)
-         z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume
-         z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume
-         z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area
-         z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content
-         z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content
-         DO jk = 1, nlay_s
-            z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
-         END DO
-         DO jk = 1, nlay_i
-            z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
-         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
-         ENDIF
-      END DO
-
-      !                                                    !--------------------------------------------!
-      IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
-         !                                                 !--------------------------------------------!
-         DO jt = 1, icycle
-            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water
-            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw )
-            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
-            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
-            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
-            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
-            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
-            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
-            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
-            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
-            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
-            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
+      ! --- transport --- !
+      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
+      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
+
+      DO jt = 1, icycle
+
+         ! record at_i before advection (for open water)
+         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
+         
+         ! --- transported fields --- !                                        
+         DO jl = 1, jpl
+            zarea(:,:,jl) = e1e2t(:,:)
+            z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume
+            z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume
+            z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area
+            z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content
+            z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content
+            DO jk = 1, nlay_s
+               z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
+            END DO
+            DO jk = 1, nlay_i
+               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
+            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
+            ENDIF
+         END DO
+         !
+         !                                                                  !--------------------------------------------!
+         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
+            !                                                               !--------------------------------------------!
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
             !
-            DO jk = 1, nlay_s                                                                             !--- snow heat content
-               CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
-                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
-               CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
-                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
-            END DO
-            DO jk = 1, nlay_i                                                                             !--- ice heat content
-               CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
-                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
-               CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
-                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+            DO jk = 1, nlay_s                                                                           !--- snow heat content
+               CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+            END DO
+            DO jk = 1, nlay_i                                                                           !--- ice heat content
+               CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
             END DO
             !
             IF ( ln_pnd_H12 ) THEN
-               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
-               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
-               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
-               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
+               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 ) 
             ENDIF
-         END DO
-      !                                                    !--------------------------------------------!
-      ELSE                                                 !== even ice time step:  adv_y then adv_x  ==!
-         !                                                 !--------------------------------------------!
-         DO jt = 1, icycle
-            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water
-            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw )
-            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
-            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
-            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
-            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
-            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
-            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
-            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
-            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
-            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
-            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
-            DO jk = 1, nlay_s                                                                             !--- snow heat content
-               CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
-                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
-               CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
-                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
-            END DO
-            DO jk = 1, nlay_i                                                                             !--- ice heat content
-               CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
-                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
-               CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
-                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+            !                                                               !--------------------------------------------!
+         ELSE                                                               !== even ice time step:  adv_y then adv_x  ==!
+            !                                                               !--------------------------------------------!
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
+            DO jk = 1, nlay_s                                                                           !--- snow heat content
+               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+               CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+            END DO
+            DO jk = 1, nlay_i                                                                           !--- ice heat content
+               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+               CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
             END DO
             IF ( ln_pnd_H12 ) THEN
-               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
-               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
-               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
-               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
+               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 )
             ENDIF
-         END DO
-      ENDIF
-
-      !-------------------------------------------
-      ! Recover the properties from their contents
-      !-------------------------------------------
-      pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1)
-      DO jl = 1, jpl
-         pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         DO jk = 1, nlay_s
-            pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         END DO
-         DO jk = 1, nlay_i
-            pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         END DO
-         IF ( ln_pnd_H12 ) THEN
-            pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-            pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            !
          ENDIF
+
+         ! --- Recover the properties from their contents --- !
+         DO jl = 1, jpl
+            pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            DO jk = 1, nlay_s
+               pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            END DO
+            DO jk = 1, nlay_i
+               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            END DO
+            IF ( ln_pnd_H12 ) THEN
+               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            ENDIF
+         END DO
+         !
+         ! derive open water from ice concentration
+         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water
+                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. )
+         !
+         ! --- 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 )
+         !
+         ! --- Ensure snow load is not too big --- !
+         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
+         !
       END DO
-      !
-      ! --- 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 )
-      !
-      ! --- Ensure snow load is not too big --- !
-      CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
       !
       IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file
@@ -296,5 +304,5 @@
             DO ji = 1, jpi
                zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
-               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji,jj,jl)
+               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
                zalfq        =  zalf * zalf
                zalf1        =  1.0 - zalf
@@ -322,5 +330,5 @@
          DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0.
             DO ji = 1, fs_jpim1
-               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji+1,jj,jl) 
+               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
                zalg  (ji,jj) = zalf
                zalfq         = zalf * zalf
@@ -465,5 +473,5 @@
             DO ji = fs_2, fs_jpim1
                zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
-               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt * e1v(ji,jj) / psm(ji,jj,jl)
+               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
                zalfq        =  zalf * zalf
                zalf1        =  1.0 - zalf
@@ -491,5 +499,5 @@
          DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
             DO ji = fs_2, fs_jpim1
-               zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * pdt * e1v(ji,jj) ) / psm(ji,jj+1,jl) 
+               zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 
                zalg  (ji,jj) = zalf
                zalfq         = zalf * zalf
@@ -645,6 +653,5 @@
       !
       !                             !* allocate prather fields
-      ALLOCATE( sxopw(jpi,jpj,1)   , syopw(jpi,jpj,1)   , sxxopw(jpi,jpj,1)   , syyopw(jpi,jpj,1)   , sxyopw(jpi,jpj,1)   ,   &
-         &      sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
+      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
          &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
          &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
@@ -692,5 +699,5 @@
          !                                   !==========================!
          !
-         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )    ! file exist: id1>0
+         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. )    ! file exist: id1>0
          ELSE                   ;   id1 = 0                                                  ! no restart: id1=0
          ENDIF
@@ -728,10 +735,4 @@
             CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
             CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
-            !                                                        ! open water in sea ice
-            CALL iom_get( numrir, jpdom_autoglo, 'sxopw' , sxopw  )
-            CALL iom_get( numrir, jpdom_autoglo, 'syopw' , syopw  )
-            CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw )
-            CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw )
-            CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw )
             !                                                        ! snow layers heat content
             DO jk = 1, nlay_s
@@ -776,5 +777,4 @@
             sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
             sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
-            sxopw = 0._wp   ;   syopw = 0._wp   ;   sxxopw = 0._wp   ;   syyopw = 0._wp   ;   sxyopw = 0._wp      ! open water in sea ice
             sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
             sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
@@ -825,10 +825,4 @@
          CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
          CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
-         !                                                           ! open water in sea ice
-         CALL iom_rstput( iter, nitrst, numriw, 'sxopw' , sxopw  )
-         CALL iom_rstput( iter, nitrst, numriw, 'syopw' , syopw  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw )
-         CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw )
          !                                                           ! snow layers heat content
          DO jk = 1, nlay_s
Index: NEMO/releases/release-4.0.1/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/releases/release-4.0.1/src/ICE/icedyn_rdgrft.F90	(revision 11697)
+++ NEMO/releases/release-4.0.1/src/ICE/icedyn_rdgrft.F90	(revision 11730)
@@ -86,7 +86,7 @@
       !!                ***  ROUTINE ice_dyn_rdgrft_alloc ***
       !!-------------------------------------------------------------------
-      ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   &
-         &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),   &
-         &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  &
+      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,               &
+         &      apartf(jpij,0:jpl) , hrmin  (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), &
+         &      hrmax (jpij,jpl)   , hi_hrdg(jpij,jpl) , araft(jpij,jpl) ,                   &
          &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
 
@@ -137,5 +137,4 @@
       REAL(wp) ::   zfac                       ! local scalar
       INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not
-      REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s)
       REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
       !
@@ -175,7 +174,7 @@
         
          ! just needed here
-         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i  )
          CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i )
          ! needed here and in the iteration loop
+         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i) ! zdivu is used as a work array here (no change in divu_i)
          CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   )
          CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   )
@@ -187,16 +186,8 @@
             closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
             !
-            ! divergence given by the advection scheme
-            !   (which may not be equal to divu as computed from the velocity field)
-            IF    ( ln_adv_Pra ) THEN
-               zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
-            ELSEIF( ln_adv_UMx ) THEN
-               zdivu_adv(ji) = zdivu(ji)
-            ENDIF
-            !
-            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough
-            !                                                                                        ! to give asum = 1.0 after ridging
+            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough
+            !                                                                                ! to give asum = 1.0 after ridging
             ! Opening rate (non-negative) that will give asum = 1.0 after ridging.
-            opning(ji) = closing_net(ji) + zdivu_adv(ji)
+            opning(ji) = closing_net(ji) + zdivu(ji)
          END DO
          !
@@ -215,5 +206,5 @@
                ato_i_1d   (ipti)   = ato_i_1d   (ji)
                closing_net(ipti)   = closing_net(ji)
-               zdivu_adv  (ipti)   = zdivu_adv  (ji)
+               zdivu      (ipti)   = zdivu      (ji)
                opning     (ipti)   = opning     (ji)
             ENDIF
@@ -259,7 +250,7 @@
                ELSE
                   iterate_ridging  = 1
-                  zdivu_adv  (ji) = zfac * r1_rdtice
-                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) )
-                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) )
+                  zdivu      (ji) = zfac * r1_rdtice
+                  closing_net(ji) = MAX( 0._wp, -zdivu(ji) )
+                  opning     (ji) = MAX( 0._wp,  zdivu(ji) )
                ENDIF
             END DO
@@ -309,5 +300,5 @@
 
       !                       ! Ice thickness needed for rafting
-      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
+      WHERE( pa_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
       ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
       END WHERE
@@ -328,5 +319,5 @@
       zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
       !
-      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
+      WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
       ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp
       END WHERE
@@ -454,5 +445,5 @@
       ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.  
       ! NOTE: 0 < aksum <= 1
-      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
+      WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
       ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
       END WHERE
@@ -537,5 +528,5 @@
             IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
 
-               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
+               IF( a_i_2d(ji,jl1) > epsi10 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
                ELSE                                 ;   z1_ai(ji) = 0._wp
                ENDIF
@@ -595,5 +586,5 @@
                ! virtual salt flux to keep salinity constant
                IF( nn_icesal /= 2 )  THEN
-                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
+                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )       ! ridge salinity = s_i
                   sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice  &  ! put back sss_m into the ocean
                      &                            - s_i_1d(ji) * vsw * rhoi * r1_rdtice     ! and get  s_i  from the ocean 
Index: NEMO/releases/release-4.0.1/src/ICE/iceitd.F90
===================================================================
--- NEMO/releases/release-4.0.1/src/ICE/iceitd.F90	(revision 11697)
+++ NEMO/releases/release-4.0.1/src/ICE/iceitd.F90	(revision 11730)
@@ -211,5 +211,5 @@
                CALL itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in
                   &              g0  (1:npti,1), g1  (1:npti,1), hL     (1:npti,1), hR    (1:npti,1)   )   ! out
-                  !
+               !
                ! Area lost due to melting of thin ice
                DO ji = 1, npti
@@ -218,5 +218,5 @@
                      !
                      zdh0 =  h_i_1d(ji) - h_ib_1d(ji)                
-                     IF( zdh0 < 0.0 ) THEN      !remove area from category 1
+                     IF( zdh0 < 0.0 ) THEN      ! remove area from category 1
                         zdh0 = MIN( -zdh0, hi_max(1) )
                         !Integrate g(1) from 0 to dh0 to estimate area melted
@@ -226,8 +226,7 @@
                            zx1    = zetamax
                            zx2    = 0.5 * zetamax * zetamax 
-                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed
+                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                ! ice area removed
                            zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i                
-                           zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting 
-                           !     of thin ice (zdamax > 0)
+                           zda0   = MIN( zda0, zdamax )                            ! ice area lost due to melting of thin ice (zdamax > 0)
                            ! Remove area, conserving volume
                            h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
@@ -349,5 +348,5 @@
       DO ji = 1, npti
          !
-         IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN
+         IF( paice(ji) > epsi10  .AND. phice(ji) > epsi10 )  THEN
             !
             ! Initialize hL and hR
Index: NEMO/releases/release-4.0.1/src/ICE/icevar.F90
===================================================================
--- NEMO/releases/release-4.0.1/src/ICE/icevar.F90	(revision 11697)
+++ NEMO/releases/release-4.0.1/src/ICE/icevar.F90	(revision 11730)
@@ -622,5 +622,5 @@
                   pv_s   (ji,jj,jl) = 0._wp
                ENDIF
-               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN
                   sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
                   psv_i  (ji,jj,jl) = 0._wp
