Index: /NEMO/trunk/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml
===================================================================
--- /NEMO/trunk/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml	(revision 9942)
+++ /NEMO/trunk/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml	(revision 9943)
@@ -79,4 +79,9 @@
        <field field_ref="vfxsnw"           name="vfxsnw" />
        
+       <!-- diag error for negative ice volume after advection -->
+       <field field_ref="iceneg_pres"      name="sineg_pres" />
+       <field field_ref="iceneg_volu"      name="sineg_volu" />
+       <field field_ref="iceneg_hfx"       name="sineg_hfx"  />
+
        <!-- categories -->
        <field field_ref="icemask_cat"      name="simskcat"/>
Index: /NEMO/trunk/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/trunk/cfgs/SHARED/field_def_nemo-ice.xml	(revision 9942)
+++ /NEMO/trunk/cfgs/SHARED/field_def_nemo-ice.xml	(revision 9943)
@@ -158,4 +158,9 @@
           <!-- diags -->
           <field id="hfxdhc"       long_name="Heat content variation in snow and ice (neg = ice cooling)"   unit="W/m2" />
+
+	  <!-- diagnostics of the negative values resulting from the advection scheme -->
+	  <field id="iceneg_pres"  long_name="Fraction of time steps with negative sea ice volume"                    unit=""     />
+	  <field id="iceneg_volu"  long_name="Negative sea ice volume per area arising from advection"                unit="m"    />
+          <field id="iceneg_hfx"   long_name="Negative sea ice heat content (eq. heat flux) arising from advection"   unit="W/m2" />
 
 	  <!-- sbcssm variables -->
Index: /NEMO/trunk/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml
===================================================================
--- /NEMO/trunk/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml	(revision 9942)
+++ /NEMO/trunk/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml	(revision 9943)
@@ -78,4 +78,9 @@
        <field field_ref="vfxice"           name="vfxice" />
        <field field_ref="vfxsnw"           name="vfxsnw" />
+
+       <!-- diag error for negative ice volume after advection -->
+       <field field_ref="iceneg_pres"      name="sineg_pres" />
+       <field field_ref="iceneg_volu"      name="sineg_volu" />
+       <field field_ref="iceneg_hfx"       name="sineg_hfx"  />
        
        <!-- categories -->
Index: /NEMO/trunk/src/ICE/icectl.F90
===================================================================
--- /NEMO/trunk/src/ICE/icectl.F90	(revision 9942)
+++ /NEMO/trunk/src/ICE/icectl.F90	(revision 9943)
@@ -148,12 +148,13 @@
             IF ( ABS( zs   ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs
             IF ( ABS( zt   ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt
-            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN
-                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
-                                           WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp
-            ENDIF
             IF ( zvmin < -epsi10 )         WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
             IF ( zamax > MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' )  &
                &                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
             IF ( zamin < -epsi10 )         WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+!clem: the following check fails when using UM3-5 advection scheme (see comments in icedyn_adv.F90)
+!            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN
+!                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
+!                                           WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp
+!            ENDIF
          ENDIF
          !
Index: /NEMO/trunk/src/ICE/icedyn_adv.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv.F90	(revision 9942)
+++ /NEMO/trunk/src/ICE/icedyn_adv.F90	(revision 9943)
@@ -66,4 +66,7 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! number of iteration
+      !
+      INTEGER ::   jl   ! dummy loop indice
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmask  ! fraction of time step with v_i < 0
       !!---------------------------------------------------------------------
       !
@@ -98,16 +101,23 @@
       ! Debug the advection schemes
       !----------------------------
-      ! clem: The 2 advection schemes above are not strictly positive.
-      !       In Prather, advected fields are bounded by 0 (not anymore?) in the routine with a MAX(0,field) ==> likely conservation issues
-      !       In UMx    , advected fields are not bounded and negative values can appear.
+      ! clem: At least one advection scheme above is not strictly positive => UM from 3d to 5th order
+      !       In Prather, I am not sure if the fields are bounded by 0 or not (it seems not)
+      !       In UM3-5  , advected fields are not bounded and negative values can appear.
       !                   These values are usually very small but in some occasions they can also be non-negligible
       !                   Therefore one needs to bound the advected fields by 0 (though this is not a clean fix)
-      ! ==> 1) remove negative ice areas and volumes (conservation is ensure)
-      CALL ice_var_zapsmall 
-      ! ==> 2) remove remaining negative advected fields (conservation is not preserved) => conservation issue
-      WHERE( v_s (:,:,:)   < 0._wp )   v_s (:,:,:)   = 0._wp
-      WHERE( sv_i(:,:,:)   < 0._wp )   sv_i(:,:,:)   = 0._wp
-      WHERE( e_i (:,:,:,:) < 0._wp )   e_i (:,:,:,:) = 0._wp
-      WHERE( e_s (:,:,:,:) < 0._wp )   e_s (:,:,:,:) = 0._wp
+      !
+      ! record the negative values resulting from UMx
+      zmask(:,:) = 0._wp ! keep the init to 0 here
+      DO jl = 1, jpl
+         WHERE( v_i(:,:,jl) < 0._wp )   zmask(:,:) = 1._wp
+      END DO
+      IF( iom_use('iceneg_pres') )   CALL iom_put("iceneg_pres", zmask                                      )  ! fraction of time step with v_i < 0
+      IF( iom_use('iceneg_volu') )   CALL iom_put("iceneg_volu", SUM(MIN( v_i, 0. ), dim=3 )                )  ! negative ice volume (only)
+      IF( iom_use('iceneg_hfx' ) )   CALL iom_put("iceneg_hfx" , ( SUM(SUM( MIN( e_i(:,:,1:nlay_i,:), 0. )  &  ! negative ice heat content (only)
+         &                                                                  , dim=4 ), dim=3 ) )* r1_rdtice )  ! -- eq. heat flux [W/m2]
+      !
+      ! ==> conservation is ensured by calling this routine below,
+      !     however the global ice volume is then changed by advection (but errors are very small) 
+      CALL ice_var_zapneg( ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
 
       !------------
Index: /NEMO/trunk/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_pnd.F90	(revision 9942)
+++ /NEMO/trunk/src/ICE/icethd_pnd.F90	(revision 9943)
@@ -163,5 +163,4 @@
             !--- Pond gowth ---!
             ! v_ip should never be negative, otherwise code crashes
-            ! MV: as far as I saw, UM5 can create very small negative v_ip values (not Prather)
             v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt )
             !
Index: /NEMO/trunk/src/ICE/icevar.F90
===================================================================
--- /NEMO/trunk/src/ICE/icevar.F90	(revision 9942)
+++ /NEMO/trunk/src/ICE/icevar.F90	(revision 9943)
@@ -44,4 +44,5 @@
    !!   ice_var_salprof1d : salinity profile in the ice 1D
    !!   ice_var_zapsmall  : remove very small area and volume
+   !!   ice_var_zapneg    : remove negative ice fields (to debug the advection scheme UM3-5)
    !!   ice_var_itd       : convert 1-cat to jpl-cat
    !!   ice_var_itd2      : convert N-cat to jpl-cat
@@ -68,4 +69,5 @@
    PUBLIC   ice_var_salprof1d    
    PUBLIC   ice_var_zapsmall
+   PUBLIC   ice_var_zapneg
    PUBLIC   ice_var_itd
    PUBLIC   ice_var_itd2
@@ -462,15 +464,15 @@
       DO jl = 1, jpl       !==  loop over the categories  ==!
          !
+         WHERE( a_i(:,:,jl) > epsi10 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
+         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
+         END WHERE
+         !
+         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
+         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
+         END WHERE
+         !
          !-----------------------------------------------------------------
          ! Zap ice energy and use ocean heat to melt ice
          !-----------------------------------------------------------------
-         WHERE( a_i(:,:,jl) > epsi10 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
-         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
-         END WHERE
-         !
-         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
-         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
-         END WHERE
-         !
          DO jk = 1, nlay_i
             DO jj = 1 , jpj
@@ -495,4 +497,7 @@
          END DO
          !
+         !-----------------------------------------------------------------
+         ! zap ice and snow volume, add water and salt to ocean
+         !-----------------------------------------------------------------
          DO jj = 1 , jpj
             DO ji = 1 , jpi
@@ -502,7 +507,4 @@
                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_rdtice
                !
-               !-----------------------------------------------------------------
-               ! zap ice and snow volume, add water and salt to ocean
-               !-----------------------------------------------------------------
                a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
                v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
@@ -533,4 +535,83 @@
 
 
+   SUBROUTINE ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_var_zapneg ***
+      !!
+      !! ** Purpose :   Remove negative sea ice fields and correct fluxes
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
+      !
+      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
+      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) ::   pe_s       ! snw heat content
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
+      !!-------------------------------------------------------------------
+      !
+      WHERE( pato_i(:,:)   < 0._wp )   pato_i(:,:)   = 0._wp
+      WHERE( poa_i (:,:,:) < 0._wp )   poa_i (:,:,:) = 0._wp
+      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp
+      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
+      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
+      !                                                        but it does not change conservation, so keep it this way is ok
+      !
+      DO jl = 1, jpl       !==  loop over the categories  ==!
+         !
+         !----------------------------------------
+         ! zap ice energy and send it to the ocean
+         !----------------------------------------
+         DO jk = 1, nlay_i
+            DO jj = 1 , jpj
+               DO ji = 1 , jpi
+                  IF( pe_i(ji,jj,jk,jl) < 0._wp ) THEN
+                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+                     pe_i(ji,jj,jk,jl) = 0._wp
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         !
+         DO jk = 1, nlay_s
+            DO jj = 1 , jpj
+               DO ji = 1 , jpi
+                  IF( pe_s(ji,jj,jk,jl) < 0._wp ) THEN
+                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+                     pe_s(ji,jj,jk,jl) = 0._wp
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         !
+         !-----------------------------------------------------
+         ! zap ice and snow volume, add water and salt to ocean
+         !-----------------------------------------------------
+         DO jj = 1 , jpj
+            DO ji = 1 , jpi
+              IF( pv_i(ji,jj,jl) < 0._wp ) THEN
+                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice
+                  pv_i   (ji,jj,jl) = 0._wp
+               ENDIF
+               IF( pv_s(ji,jj,jl) < 0._wp ) THEN
+                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice
+                  pv_s   (ji,jj,jl) = 0._wp
+               ENDIF
+               IF( psv_i(ji,jj,jl) < 0._wp ) THEN
+                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice
+                  psv_i  (ji,jj,jl) = 0._wp
+               ENDIF
+            END DO
+         END DO
+         !
+      END DO 
+      !
+   END SUBROUTINE ice_var_zapneg
+
+   
    SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
       !!-------------------------------------------------------------------
