Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/INSTALL.rst
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/INSTALL.rst	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/INSTALL.rst	(revision 11082)
@@ -104,5 +104,5 @@
    .. code:: console
 
-      $ svn co http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/branchs/xios-2.5
+      $ svn co https://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/branchs/xios-2.5
 
 Download the NEMO source code
@@ -111,5 +111,5 @@
 .. code:: console
 
-   $ svn co http://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/release-4.0
+   $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/release-4.0
 
 Description of directory tree
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/README.rst
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/README.rst	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/README.rst	(revision 11082)
@@ -1,4 +1,5 @@
-:Release: |version|
-:Date:    |today|
+:Release:  |release|
+:Date:     |today|
+:SVN rev.: |revision|
 
 NEMO_ for **Nucleus for European Modelling of the Ocean** is a state-of-the-art modelling framework for
@@ -90,8 +91,8 @@
 .. Substitutions / Links
 
-.. |NEMO manual| image:: http://zenodo.org/badge/DOI/10.5281/zenodo.1464816.svg
-.. |NEMO guide|  image:: http://zenodo.org/badge/DOI/10.5281/zenodo.1475325.svg
-.. |SI3 manual|  image:: http://zenodo.org/badge/DOI/10.5281/zenodo.1471689.svg
-.. |TOP manual|  image:: http://zenodo.org/badge/DOI/10.5281/zenodo.1471700.svg
+.. |NEMO manual| image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1464816.svg
+.. |NEMO guide|  image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1475325.svg
+.. |SI3 manual|  image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1471689.svg
+.. |TOP manual|  image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1471700.svg
 
 .. |NEMO strategy| replace:: multi-year development strategy
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml	(revision 11082)
@@ -79,9 +79,4 @@
        <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/branches/UKMO/NEMO_4.0_GO8_package/cfgs/README.rst
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/README.rst	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/README.rst	(revision 11082)
@@ -235,5 +235,5 @@
 .. literalinclude:: ../../../cfgs/GYRE_PISCES/EXPREF/namelist_ref
    :language: fortran
-   :lines: 306-333
+   :lines: 935-960
 
 Input dynamical fields for this configuration (``ORCA2_OFF_v4.0.tar``) comes from
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/field_def_nemo-ice.xml	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/field_def_nemo-ice.xml	(revision 11082)
@@ -160,9 +160,4 @@
           <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 -->
           <field id="sst_m"    unit="degC" />
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/namelist_ice_ref	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/namelist_ice_ref	(revision 11082)
@@ -154,5 +154,5 @@
 &namthd_do      !   Ice growth in open water
 !------------------------------------------------------------------------------
-   rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_hnewice
+   rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_himin
    ln_frazil        = .false.         !  Frazil ice parameterization (ice collection as a function of wind)
       rn_maxfraz    =   1.0           !     maximum fraction of frazil ice collecting at the ice base
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml	(revision 11082)
@@ -79,9 +79,4 @@
        <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/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SPITZ12/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 11082)
@@ -80,4 +80,6 @@
                      ! Sea-ice :
    nn_ice      = 2         !  SI3
+   ln_ice_embd = .false.   !  =T embedded sea-ice (pressure + mass and salt exchanges)
+      !                    !  =F levitating ice (no pressure, mass and salt exchanges)
                      ! Misc. options of sbc :
    ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr )
@@ -90,5 +92,5 @@
       !
       ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012)
-      ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015)
+      ln_Cd_L15   = .true.    !  air-ice drags = F(ice concentration) (Lupkes et al. 2015)
       !
    cn_dir = './'  !  root directory for the bulk data location
@@ -96,16 +98,4 @@
    !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !          weights filename            ! rotation ! land/sea mask !
    !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   !
-!! ERAI
-!   sn_wndi     = 'u10_era_spitz'               ,         6         , 'u10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Uwnd' , ''
-!   sn_wndj     = 'v10_era_spitz'               ,         6         , 'v10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Vwnd' , ''
-!   sn_qsr      = 'ssrd_era_spitz'              ,         6         , 'ssrd'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_qlw      = 'strd_era_spitz'              ,         6         , 'strd'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_tair     = 't2_era_spitz'                ,         6         , 't2'      ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_humi     = 'humi_era_spitz'              ,         6         , 'humi'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_prec     = 'precip_era_spitz'            ,         6         , 'precip'  ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_snow     = 'snow_era_spitz'              ,         6         , 'snow'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_tdif     = 'taudif'                      ,         6         , 'taudif'  ,   .true.    , .false. , 'yearly'  , ''             , '' , ''
-!      rn_zqt      = 10.       !  Air temperature & humidity reference height (m)
-!! MAR
    sn_wndi     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,   2 ,  'u10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Uwnd' , ''
    sn_wndj     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,   2 ,  'v10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Vwnd' , ''
@@ -129,12 +119,4 @@
    !                       ! type of penetration                        (default: NO selection)
    ln_qsr_rgb  = .true.       !  RGB light penetration (Red-Green-Blue)
-/
-!-----------------------------------------------------------------------
-&namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr =T)
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namsbc_rnf    !   runoffs                                              (ln_rnf =T)
-!-----------------------------------------------------------------------
 /
 !!======================================================================
@@ -236,19 +218,13 @@
 !-----------------------------------------------------------------------
    ln_loglayer = .true.   !  logarithmic drag: Cd = vkarmn/log(z/z0) |U|
+   ln_drgimp   = .true.   !  implicit top/bottom friction flag
 /
 !-----------------------------------------------------------------------
 &namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F)
 !-----------------------------------------------------------------------
-   rn_Cd0      =  2.5e-3   !!1.e-3    !  drag coefficient [-]
-   rn_Cdmax    =  0.01     !!0.1      !  drag value maximum [-] (logarithmic drag)
-   rn_ke0      =  0.       !!2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases)
+   rn_Cd0      =  2.5e-3   !  drag coefficient [-]
+   rn_Cdmax    =  0.1      !  drag value maximum [-] (logarithmic drag)
+   rn_ke0      =  0.       !  background kinetic energy  [m2/s2] (non-linear cases)
    rn_z0       =  3.e-3    !  roughness [m] (ln_loglayer=T)
-   ln_boost    = .false.   !  =T regional boost of Cd0 ; =F constant
-      rn_boost =  50.         !  local boost factor  [-]
-/
-!-----------------------------------------------------------------------
-&nambbc        !   bottom temperature boundary condition                (default: OFF)
-!-----------------------------------------------------------------------
-   ln_trabbc   = .false.    !  Apply a geothermal heating at the ocean bottom
 /
 !-----------------------------------------------------------------------
@@ -258,6 +234,4 @@
       nn_bbl_ldf  =  1        !  diffusive bbl (=1)   or not (=0)
       nn_bbl_adv  =  0        !  advective bbl (=1/2) or not (=0)
-      rn_ahtbbl   =  1000.    !  lateral mixing coefficient in the bbl  [m2/s]
-      rn_gambbl   =  10.      !  advective bbl coefficient                 [s]
 /
 !!======================================================================
@@ -293,5 +267,5 @@
    !                       !  Coefficients:
    nn_aht_ijk_t    = 31        !  space/time variation of eddy coefficient:
-      !                             !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing)
+   !                                !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing)
 /
 !!======================================================================
@@ -320,5 +294,5 @@
 &namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_dynvor_enT = .true. !  energy conserving scheme (T-point)
+   ln_dynvor_eeT = .true.  !  energy conserving scheme (een using e3t)
 /
 !-----------------------------------------------------------------------
@@ -352,8 +326,14 @@
 &namzdf        !   vertical physics manager                             (default: NO selection)
 !-----------------------------------------------------------------------
+   !                       ! adaptive-implicit vertical advection
+   ln_zad_Aimp = .true.      !  Courant number dependent scheme (Shchepetkin 2015)
    !                       ! type of vertical closure (required)
    ln_zdftke   = .true.       !  Turbulent Kinetic Energy closure       (T =>   fill namzdf_tke)
    ln_zdfgls   = .false.      !  Generic Length Scale closure           (T =>   fill namzdf_gls)
+   !                       ! convection
+   ln_zdfevd   = .true.       !  enhanced vertical diffusion
+   !
    ln_zdfddm   = .true.    ! double diffusive mixing
+   !
    !                       !  Coefficients
    rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
@@ -384,8 +364,4 @@
 /
 !-----------------------------------------------------------------------
-&nam_diaharm   !   Harmonic analysis of tidal constituents              ('key_diaharm')
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4")
 !-----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SPITZ12/EXPREF/namelist_ice_cfg
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SPITZ12/EXPREF/namelist_ice_cfg	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SPITZ12/EXPREF/namelist_ice_cfg	(revision 11082)
@@ -35,4 +35,5 @@
 &namdyn         !   Ice dynamics
 !------------------------------------------------------------------------------
+   ln_landfast_L16  = .true.          !  landfast: parameterization from Lemieux 2016
 /
 !------------------------------------------------------------------------------
@@ -43,4 +44,6 @@
 &namdyn_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
+   ln_rhg_EVP       = .true.          !  EVP rheology
+      ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017)
 /
 !------------------------------------------------------------------------------
@@ -67,5 +70,5 @@
 &namthd_do      !   Ice growth in open water
 !------------------------------------------------------------------------------
-   rn_hinew         =   0.02          !  thickness for new ice formation in open water (m), must be larger than rn_hnewice
+   rn_hinew         =   0.02          !  thickness for new ice formation in open water (m), must be larger than rn_himin
    ln_frazil        = .true.          !  Frazil ice parameterization (ice collection as a function of wind)
 /
@@ -77,7 +80,6 @@
 &namthd_pnd     !   Melt ponds
 !------------------------------------------------------------------------------
-   ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
-   ln_pnd_CST       = .false.         !  activate constant melt ponds
-   ln_pnd_alb       = .false.         !  melt ponds affect albedo or not
+   ln_pnd_H12       = .true.          !  activate evolutive melt ponds (from Holland et al 2012)
+   ln_pnd_alb       = .true.          !  melt ponds affect albedo or not
 /
 
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/Makefile
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/Makefile	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/Makefile	(revision 11082)
@@ -23,3 +23,3 @@
 # Browse to 127.0.0.1:8000/NEMO_guide.html
 htmllive:
-	sphinx-autobuild $(SPHINXOPTS) $(SOURCEDIR) $(BUILDDIR)/html
+	sphinx-autobuild $(SPHINXOPTS) $(SOURCEDIR) $(BUILDDIR)/htmllive
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/README.rst
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/README.rst	(revision 11082)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/README.rst	(revision 11082)
@@ -0,0 +1,54 @@
+**********************
+NEMO Quick Start Guide
+**********************
+
+| The NEMO guide is made up of several files written in `ReStructuredText <http://docutils.sourceforge.net/rst.html>`_ (`.rst` extension), a WYSIWYG markup language used in the Python community, and scattered all over the NEMO sources.
+| You can view them one by one in plain text from `./source` folder, or export all to a user-friendly guide under `./build` (only HTML format at the moment, PDF expected later).
+
+Build and export the guide in HTML
+==================================
+
+1.  Install Sphinx documentation generator, its BibTeX extension and "Read The Docs" theme thanks to `pip` packages tool
+
+    .. code-block:: console
+
+        $ pip install sphinx sphinxcontrib.bibtex sphinx_rtd_theme
+
+2.  Build the HTML export with `make` in `./build/html`
+
+    .. code-block:: console
+
+        $ make html
+
+3.  Finally browse the guide by opening `./build/html/NEMO_guide.html`
+
+
+Edit the sources and check the output in real time
+==================================================
+
+| To facilitate the update of the guide, editors can install a useful package that will automatically trigger a new build and the reload of the HTML page for every recorded change in the sources.
+| So the reviewer saves time by controlling on-line their modifications almost as it types and also by avoiding repeated interactive rebuilds.
+
+Install `sphinx-autobuild` package
+
+.. code-block:: console
+
+	$ pip install sphinx-autobuild
+
+Launch a local web server hosting a draft export of the guide (build this time in `./build/livehtml`)
+
+.. code-block:: console
+
+    $ make livehtml
+
+| Open in the same time the 2 formats of the content to review: the source file and the web page by browsing from the new guide hosted by your local server on `<http://127.0.0.1:8000/NEMO_guide.html>`_.
+| Start the update, save your changes and verify instantly the HTML export in your browser.
+
+.. warning::
+
+    | Your modifications are not taken into account?
+    | For symlink file, you will have to close it to update the HTML export. Otherwise look at the log of the Sphinx build, you probably made a typo!
+
+.. hint::
+
+    Are there broken links? Fix "Page not found" errors by running `make linkcheck`
Index: MO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/requirements.txt
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/requirements.txt	(revision 11081)
+++ 	(revision )
@@ -1,52 +1,0 @@
-alabaster==0.7.12
-Babel==2.6.0
-Brlapi==0.6.5
-certifi==2018.11.29
-chardet==3.0.4
-cryptography==1.7.1
-cupshelpers==1.0
-docutils==0.14
-httplib2==0.9.2
-idna==2.8
-imagesize==1.1.0
-Jinja2==2.10
-keyring==10.1
-keyrings.alt==1.3
-latexcodec==1.0.5
-louis==3.0.0
-Mako==1.0.6
-MarkupSafe==1.1.0
-oset==0.1.3
-packaging==18.0
-pexpect==4.2.1
-Pillow==4.0.0
-ptyprocess==0.5.1
-pyasn1==0.1.9
-pybtex==0.22.2
-pybtex-docutils==0.2.1
-pycrypto==2.6.1
-pycups==1.9.73
-pycurl==7.43.0
-Pygments==2.3.1
-pygobject==3.22.0
-pyparsing==2.3.1
-pysmbc==1.0.15.6
-python-apt==1.4.0b3
-python-debian==0.1.30
-python-debianbts==2.6.1
-pytz==2018.9
-pyxdg==0.25
-PyYAML==3.13
-reportbug==7.1.7
-reportlab==3.3.0
-requests==2.21.0
-roman==2.0.0
-SecretStorage==2.3.1
-six==1.12.0
-snowballstemmer==1.2.1
-Sphinx==1.8.3
-sphinx-rtd-theme==0.4.2
-sphinxcontrib-bibtex==0.4.2
-sphinxcontrib-websupport==1.1.0
-unattended-upgrades==0.1
-urllib3==1.24.1
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/source/_templates/layout.html
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/source/_templates/layout.html	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/source/_templates/layout.html	(revision 11082)
@@ -12,5 +12,5 @@
 	<p>Community ocean model for multifarious space and time scales</p>
 
-	<div class="version">{{ release }}</div>
+	<div class="version">{{ version }}</div>
 
 	{% include "searchbox.html" %}
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/source/conf.py
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/source/conf.py	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/doc/rst/source/conf.py	(revision 11082)
@@ -92,4 +92,5 @@
 # -- Customisation -----------------------------------------------------------
 
+# Timestamping
 import datetime
 year = datetime.date.today().year
@@ -108,2 +109,7 @@
 # Include common directives for every rst file
 rst_epilog = open('global.rst', 'r').read()
+
+# SVN revision
+import subprocess
+revision = subprocess.check_output("svnversion").decode("utf-8")
+rst_prolog = '.. |revision| replace:: %s' % revision
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icecor.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icecor.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icecor.F90	(revision 11082)
@@ -66,7 +66,6 @@
          WRITE(numout,*) '~~~~~~~'
       ENDIF
-      !
       !                             !-----------------------------------------------------
-      !                             !  ice thickness must exceed himin (for ice diff)    !
+      !                             !  ice thickness must exceed himin (for temp. diff.) !
       !                             !-----------------------------------------------------
       WHERE( a_i(:,:,:) >= epsi20 )   ;   h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
@@ -97,5 +96,4 @@
          END DO
       ENDIF
-
       !                             !-----------------------------------------------------
       !                             !  Rebin categories with thickness out of bounds     !
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icectl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icectl.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icectl.F90	(revision 11082)
@@ -69,5 +69,5 @@
       !!
       REAL(wp) ::   zv, zs, zt, zfs, zfv, zft
-      REAL(wp) ::   zvmin, zamin, zamax 
+      REAL(wp) ::   zvmin, zamin, zamax, zeimin, zesmin, zsmin
       REAL(wp) ::   zvtrp, zetrp
       REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill
@@ -141,7 +141,10 @@
          zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv
 
-         zvmin = glob_min( 'icectl', v_i )
-         zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) )
-         zamin = glob_min( 'icectl', a_i )
+         zamax  = glob_max( 'icectl', SUM( a_i, dim=3 ) )
+         zvmin  = glob_min( 'icectl', v_i )
+         zamin  = glob_min( 'icectl', a_i )
+         zsmin  = glob_min( 'icectl', sv_i )
+         zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) )
+         zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) )
 
          ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
@@ -152,13 +155,18 @@
 
          IF(lwp) THEN
-            IF ( ABS( zv   ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv
-            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 ( 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' .AND. cd_routine /= 'Hbig' ) &
-               &                           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 UMx advection scheme (see comments in icedyn_adv.F90)
+            ! check conservation issues
+            IF ( ABS( zv ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv
+            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
+            ! check maximum ice concentration
+            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
+            ! check negative values
+            IF ( zvmin  < 0. )           WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
+            IF ( zamin  < 0. )           WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+            IF ( zsmin  < 0. )           WRITE(numout,*) 'violation s_i<0               (',cd_routine,') = ',zsmin
+            IF ( zeimin < 0. )           WRITE(numout,*) 'violation e_i<0               (',cd_routine,') = ',zeimin
+            IF ( zesmin < 0. )           WRITE(numout,*) 'violation e_s<0               (',cd_routine,') = ',zesmin
+!clem: the following check fails (I think...)
 !            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN
 !                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn.F90	(revision 11082)
@@ -74,8 +74,7 @@
       INTEGER, INTENT(in) ::   kt     ! ice time step
       !!
-      INTEGER  ::   ji, jj, jl        ! dummy loop indices
+      INTEGER  ::   ji, jj        ! dummy loop indices
       REAL(wp) ::   zcoefu, zcoefv
-      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdivu_i
       !!--------------------------------------------------------------------
       !
@@ -89,48 +88,38 @@
       ENDIF
       !                      
-      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
-         tau_icebfr(:,:) = 0._wp
-         DO jl = 1, jpl
-            WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
-         END DO
-      ENDIF
-      !
-      !                                !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
-      DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), &
-                  &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), &
-                  &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), &
-                  &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) )
-               zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), &
-                  &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), &
-                  &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), &
-                  &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) )
-               zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), &
-                  &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), &
-                  &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), &
-                  &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )
-      !
-      !
-      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
+      ! retrieve thickness from volume for landfast param. and UMx advection scheme
+      WHERE( a_i(:,:,:) >= epsi20 )
+         h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
+         h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
+      ELSEWHERE
+         h_i(:,:,:) = 0._wp
+         h_s(:,:,:) = 0._wp
+      END WHERE
+      !
+      WHERE( a_ip(:,:,:) >= epsi20 )
+         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
+      ELSEWHERE
+         h_ip(:,:,:) = 0._wp
+      END WHERE
+      !
+      !
+      SELECT CASE( nice_dyn )          !-- Set which dynamics is running
 
       CASE ( np_dynALL )           !==  all dynamical processes  ==!
-         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
-         CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting 
-         CALL ice_cor       ( kt , 1 )                                             ! -- Corrections
-
+         !
+         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology  
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting 
+         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections
+         !
       CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
-         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
-         CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting)
-
+         !
+         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology  
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting)
+         CALL ice_var_zapsmall                                              ! -- zap small areas
+         !
       CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
-         ALLOCATE( zdivu_i(jpi,jpj) )
+         !
          ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
          ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
@@ -145,20 +134,8 @@
          END DO
          ! ---
-         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice
-
-         ! diagnostics: divergence at T points 
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
-         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
-
-         DEALLOCATE( zdivu_i )
-
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+         !
       CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities)
-         ALLOCATE( zdivu_i(jpi,jpj) )
+         !
          u_ice(:,:) = rn_uice * umask(:,:,1)
          v_ice(:,:) = rn_vice * vmask(:,:,1)
@@ -166,19 +143,30 @@
          !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
          ! ---
-         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice
-
-         ! diagnostics: divergence at T points 
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+
+      END SELECT
+      !
+      !
+      ! diagnostics: divergence at T points 
+      IF( iom_use('icediv') ) THEN
+         !
+         SELECT CASE( nice_dyn )
+
+         CASE ( np_dynADV1D , np_dynADV2D )
+
+            ALLOCATE( zdivu_i(jpi,jpj) )
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                     &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+               END DO
             END DO
-         END DO
-         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
-         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
-
-         DEALLOCATE( zdivu_i )
-
-      END SELECT
+            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
+            CALL iom_put( "icediv" , zdivu_i(:,:) )
+            DEALLOCATE( zdivu_i )
+
+         END SELECT
+         !
+      ENDIF
       !
       ! controls
@@ -188,92 +176,4 @@
 
 
-   SUBROUTINE Hbig( phi_max, phs_max, phip_max )
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE Hbig  ***
-      !!
-      !! ** Purpose : Thickness correction in case advection scheme creates
-      !!              abnormally tick ice or snow
-      !!
-      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
-      !!                 (before advection) and reduce it by adapting ice concentration
-      !!              2- check whether snow thickness is larger than the surrounding 9-points
-      !!                 (before advection) and reduce it by sending the excess in the ocean
-      !!              3- check whether snow load deplets the snow-ice interface below sea level$
-      !!                 and reduce it by sending the excess in the ocean
-      !!              4- correct pond fraction to avoid a_ip > a_i
-      !!
-      !! ** input   : Max thickness of the surrounding 9-points
-      !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
-      !
-      INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra
-      !!-------------------------------------------------------------------
-      ! controls
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      !
-      CALL ice_var_zapsmall                       !-- zap small areas
-      !
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( v_i(ji,jj,jl) > 0._wp ) THEN
-                  !
-                  !                               ! -- check h_ip -- !
-                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
-                  IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN
-                     zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )
-                     IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN
-                        a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl)
-                     ENDIF
-                  ENDIF
-                  !
-                  !                               ! -- check h_i -- !
-                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
-                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
-                     a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
-                  ENDIF
-                  !
-                  !                               ! -- check h_s -- !
-                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
-                  zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
-                  IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
-                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
-                     !
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
-                     !
-                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
-                     v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl)
-                  ENDIF           
-                  !
-                  !                               ! -- check snow load -- !
-                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
-                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
-                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
-                  zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
-                  IF( zvs_excess > 0._wp ) THEN
-                     zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_s(ji,jj,jl), epsi20 )
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
-                     !
-                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
-                     v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess
-                  ENDIF
-                  
-               ENDIF
-            END DO
-         END DO
-      END DO 
-      !                                           !-- correct pond fraction to avoid a_ip > a_i
-      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
-      !
-      ! controls
-      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      !
-   END SUBROUTINE Hbig
-
-
    SUBROUTINE Hpiling
       !!-------------------------------------------------------------------
@@ -290,6 +190,4 @@
       ! controls
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      !
-      CALL ice_var_zapsmall                       !-- zap small areas
       !
       at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
@@ -337,16 +235,16 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namdyn:'
-         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL
-         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV
-         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D
-         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D
-         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')'
-         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat
-         WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16
-         WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home
-         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra
-         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr
-         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax
-         WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile
+         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL
+         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV
+         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D
+         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D
+         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')'
+         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat
+         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16
+         WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home
+         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra
+         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr
+         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax
+         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile
          WRITE(numout,*)
       ENDIF
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_adv.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_adv.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_adv.F90	(revision 11082)
@@ -64,10 +64,9 @@
       !!----------------------------------------------------------------------
       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
       !!---------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('icedyn_adv')
+      ! controls
+      IF( ln_timing    )   CALL timing_start('icedyn_adv')                                                             ! timing
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       !
       IF( kt == nit000 .AND. lwp ) THEN
@@ -76,44 +75,20 @@
          WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
-      
-      ! conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
-                     
-      !----------
-      ! Advection
-      !----------
+      !
+      !---------------!
+      !== Advection ==!
+      !---------------!
       SELECT CASE( nice_adv )
       !                                !-----------------------!
       CASE( np_advUMx )                ! ULTIMATE-MACHO scheme !
          !                             !-----------------------!
-         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
-      !                                !-----------------------!
+         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, &
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+         !                             !-----------------------!
       CASE( np_advPRA )                ! PRATHER scheme        !
          !                             !-----------------------!
-         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, &
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
       END SELECT
-
-      !----------------------------
-      ! Debug the advection schemes
-      !----------------------------
-      ! clem: At least one advection scheme above is not strictly positive => UMx
-      !       In Prather, I am not sure if the fields are bounded by 0 or not (it seems yes)
-      !       In UMx    , advected fields are not perfectly 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)
-      !
-      ! 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 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/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_adv_umx.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_adv_umx.F90	(revision 11082)
@@ -11,8 +11,8 @@
    !!   'key_si3'                                       SI3 sea-ice model
    !!----------------------------------------------------------------------
-   !!   ice_dyn_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme
+   !!   ice_dyn_adv_umx   : update the tracer fields
    !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders
-   !!   macho             : ???
-   !!   nonosc_ice        : compute monotonic tracer fluxes by a non-oscillatory algorithm 
+   !!   macho             : compute the fluxes
+   !!   nonosc_ice        : limit the fluxes using a non-oscillatory algorithm 
    !!----------------------------------------------------------------------
    USE phycst         ! physical constant
@@ -32,27 +32,22 @@
 
    PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90
-      
-   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
-   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
-   
-   ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid)
-   INTEGER ::   kn_limiter = 1
-
-   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
-   !   clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why)
-   !         but in realistic cases, it avoids having very negative ice temperature (-50) at low ice concentration 
-   LOGICAL ::   ll_neg = .TRUE.
-   
-   ! alternate directions for upstream
-   LOGICAL ::   ll_upsxy = .TRUE.
-
-   ! alternate directions for high order
-   LOGICAL ::   ll_hoxy = .TRUE.
-   
-   ! prelimiter: use it to avoid overshoot in H
-   !   clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why)
-   LOGICAL ::   ll_prelimiter_zalesak = .FALSE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D
-
-
+   !
+   INTEGER, PARAMETER ::   np_advS = 1         ! advection for S and T:    dVS/dt = -div(      uVS     ) => np_advS = 1
+   !                                                                    or dVS/dt = -div( uA * uHS / u ) => np_advS = 2
+   !                                                                    or dVS/dt = -div( uV * uS  / u ) => np_advS = 3
+   INTEGER, PARAMETER ::   np_limiter = 1      ! limiter: 1 = nonosc
+   !                                                      2 = superbee
+   !                                                      3 = h3
+   LOGICAL            ::   ll_upsxy  = .TRUE.   ! alternate directions for upstream
+   LOGICAL            ::   ll_hoxy   = .TRUE.   ! alternate directions for high order
+   LOGICAL            ::   ll_neg    = .TRUE.   ! if T interpolated at u/v points is negative or v_i < 1.e-6
+   !                                                 then interpolate T at u/v points using the upstream scheme
+   LOGICAL            ::   ll_prelim = .FALSE.  ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D
+   !
+   REAL(wp)           ::   z1_6   = 1._wp /   6._wp   ! =1/6
+   REAL(wp)           ::   z1_120 = 1._wp / 120._wp   ! =1/120
+   !
+   INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::   imsk_small, jmsk_small
+   !
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -64,5 +59,5 @@
 CONTAINS
 
-   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  &
+   SUBROUTINE ice_dyn_adv_umx( kn_umx, 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 )
       !!----------------------------------------------------------------------
@@ -79,4 +74,7 @@
       REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
       REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness
       REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
@@ -93,18 +91,43 @@
       INTEGER  ::   icycle                  ! number of sub-timestep for the advection
       REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
-      REAL(wp) ::   zdt
-      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
+      REAL(wp) ::   zdt, zvi_cen
+      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box
       REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zu_cat, zv_cat
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs 
       !!----------------------------------------------------------------------
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
       !
-      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !
-      !     When needed, the advection split is applied at the next time-step in order to avoid blocking global comm.
-      !     ...this should not affect too much the stability... Was ok on the tests we did...
+      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
+      DO jl = 1, jpl
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), &
+                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), &
+                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
+                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
+               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), &
+                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), &
+                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
+                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
+               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), &
+                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), &
+                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
+                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
+      !
+      !
+      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
+      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
+      !              this should not affect too much the stability
       zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
       zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
@@ -116,5 +139,4 @@
       ELSE                         ;   icycle = 1
       ENDIF
-      
       zdt = rdt_ice / REAL(icycle)
 
@@ -122,5 +144,11 @@
       zudy(:,:) = pu_ice(:,:) * e2u(:,:)
       zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
-
+      !
+      ! setup transport for each ice cat 
+      DO jl = 1, jpl
+         zu_cat(:,:,jl) = zudy(:,:)
+         zv_cat(:,:,jl) = zvdx(:,:)
+      END DO
+      !
       ! --- define velocity for advection: u*grad(H) --- !
       DO jj = 2, jpjm1
@@ -154,73 +182,180 @@
          END WHERE
          !
-         ! set u*a=u for advection of A only 
-         DO jl = 1, jpl
-            zua_ho(:,:,jl) = zudy(:,:)
-            zva_ho(:,:,jl) = zvdx(:,:)
-         END DO
-         
+         ! setup a mask where advection will be upstream
+         IF( ll_neg ) THEN
+            IF( .NOT. ALLOCATED(imsk_small) )   ALLOCATE( imsk_small(jpi,jpj,jpl) ) 
+            IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) ) 
+            DO jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) )
+                     IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0
+                     ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF
+                     zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) )
+                     IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0
+                     ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF
+                  END DO
+               END DO
+            END DO
+         ENDIF
+         !
+         ! ----------------------- !
+         ! ==> start advection <== !
+         ! ----------------------- !
+         !
+         !== Ice area ==!
          zamsk = 1._wp
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) !-- Ice area
-         zamsk = 0._wp
-         !
-         zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i )                !-- Ice volume
-         !
-         zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s )                !-- Snw volume
-         !
-         zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i )               !-- Salt content
-         !
-         DO jk = 1, nlay_i
-            zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) )   !-- Ice heat content
-         END DO
-         !
-         DO jk = 1, nlay_s
-            zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) )   !-- Snw heat content
-         END DO
-         !
-         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN                                                              !-- Ice Age
-            ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with
-            !       fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that
-            !       spread into the domain. Therefore we cheat and consider that ice age should be advected as ice concentration
-            !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
-            !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i )
-            ! set u*a=u for advection of ice age
-            DO jl = 1, jpl
-               zua_ho(:,:,jl) = zudy(:,:)
-               zva_ho(:,:,jl) = zvdx(:,:)
-            END DO
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, &
+            &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho )
+         !
+         !                             ! --------------------------------- !
+         IF( np_advS == 1 ) THEN       ! -- advection form: -div( uVS ) -- !
+            !                          ! --------------------------------- !
+            zamsk = 0._wp
+            !== Ice volume ==!
+            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_i, zua_ups, zva_ups )
+            !== Snw volume ==!         
+            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_s, zua_ups, zva_ups )
+            !
             zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, poa_i, poa_i )
+            !== Salt content ==!
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+               &                                      psv_i, psv_i )
+            !== Ice heat content ==!
+            DO jk = 1, nlay_i
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+                  &                                      pe_i(:,:,jk,:), pe_i(:,:,jk,:) )
+            END DO
+            !== Snw heat content ==!
+            DO jk = 1, nlay_s
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+                  &                                      pe_s(:,:,jk,:), pe_s(:,:,jk,:) )
+            END DO
+            !
+            !                          ! ------------------------------------------ !
+         ELSEIF( np_advS == 2 ) THEN   ! -- advection form: -div( uA * uHS / u ) -- !
+            !                          ! ------------------------------------------ !
             zamsk = 0._wp
+            !== Ice volume ==!
+            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_i, zua_ups, zva_ups )
+            !== Snw volume ==!         
+            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_s, zua_ups, zva_ups )
+            !== Salt content ==!
+            zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, psv_i, zua_ups, zva_ups )
+            !== Ice heat content ==!
+            DO jk = 1, nlay_i
+               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
+                  &                                      zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups )
+            END DO
+            !== Snw heat content ==!
+            DO jk = 1, nlay_s
+               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
+                  &                                      zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups )
+            END DO
+            !
+            !                          ! ----------------------------------------- !
+         ELSEIF( np_advS == 3 ) THEN   ! -- advection form: -div( uV * uS / u ) -- !
+            !                          ! ----------------------------------------- !
+            zamsk = 0._wp
+            !
+            ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl),  &
+               &      zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) )
+            !
+            ! inverse of Vi
+            WHERE( pv_i(:,:,:) >= epsi20 )   ;   z1_vi(:,:,:) = 1._wp / pv_i(:,:,:)
+            ELSEWHERE                        ;   z1_vi(:,:,:) = 0.
+            END WHERE
+            ! inverse of Vs
+            WHERE( pv_s(:,:,:) >= epsi20 )   ;   z1_vs(:,:,:) = 1._wp / pv_s(:,:,:)
+            ELSEWHERE                        ;   z1_vs(:,:,:) = 0.
+            END WHERE
+            !
+            ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays)
+            !
+            !== Ice volume ==!
+            zuv_ups = zua_ups
+            zvv_ups = zva_ups
+            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
+            !== Salt content ==!
+            zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, &
+               &                                      zhvar, psv_i, zuv_ups, zvv_ups )
+            !== Ice heat content ==!
+            DO jk = 1, nlay_i
+               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
+                  &                                      zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups )
+            END DO
+            !== Snow volume ==!         
+            zuv_ups = zua_ups
+            zvv_ups = zva_ups
+            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
+            !== Snw heat content ==!
+            DO jk = 1, nlay_s
+               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
+                  &                                      zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups )
+            END DO
+            !
+            DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs )
+            !
          ENDIF
          !
-         IF ( ln_pnd_H12 ) THEN                                                                                               !-- melt ponds
-            ! set u*a=u for advection of Ap only 
-            DO jl = 1, jpl
-               zua_ho(:,:,jl) = zudy(:,:)
-               zva_ho(:,:,jl) = zvdx(:,:)
-            END DO
-            !
+         !== Ice age ==!
+         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN
             zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! fraction
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+               &                                      poa_i, poa_i )
+         ENDIF
+         !
+         !== melt ponds ==!
+         IF ( ln_pnd_H12 ) THEN
+            ! fraction
+            zamsk = 1._wp
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, &
+               &                                      pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho )
+            ! volume
             zamsk = 0._wp
-            !
             zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip )                 ! volume
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_ip, zua_ups, zva_ups )
          ENDIF
          !
+         !== Open water area ==!
          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 area
+               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 
                   &                          - ( 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_umx', pato_i(:,:), 'T',  1. )
-         !
+         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. )
+         !
+         !
+         ! --- Ensure non-negative fields and in-bound thicknesses --- !
+         ! 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 )
+         !
+         ! Make sure ice thickness is not too big
+         !    (because ice thickness can be too large where ice concentration is very small)
+         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+
       END DO
       !
@@ -228,5 +363,6 @@
 
    
-   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
+   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox,  &
+      &                                            pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE adv_umx  ***
@@ -235,40 +371,51 @@
       !!                 tracers and add it to the general trend of tracer equations
       !!
-      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracer H
+      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc
       !!                 - calculate tracer H at u and v points (Ultimate)
-      !!                 - calculate the high order fluxes using alterning directions (Macho?)
+      !!                 - calculate the high order fluxes using alterning directions (Macho)
       !!                 - apply a limiter on the fluxes (nonosc_ice)
-      !!                 - convert this tracer flux to a tracer content flux (uH -> uV)
-      !!                 - calculate the high order solution for tracer content V
+      !!                 - convert this tracer flux to a "volume" flux (uH -> uV)
+      !!                 - apply a limiter a second time on the volumes fluxes (nonosc_ice)
+      !!                 - calculate the high order solution for V
       !!
-      !! ** Action : solve 2 equations => a) da/dt = -div(ua)
-      !!                                  b) dV/dt = -div(uV) using dH/dt = -u.grad(H)
-      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H.
-      !!                        - then we convert this flux to a "volume" flux this way => uH*ua/u
-      !!                             where ua is the flux from eq. a)
-      !!                        - at last we estimate dV/dt = -div(uH*ua/u)
+      !! ** Action : solve 3 equations => a) dA/dt  = -div(uA)
+      !!                                  b) dV/dt  = -div(uV)  using dH/dt = -u.grad(H)
+      !!                                  c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S)
       !!
-      !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc.
-      !!           - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now
-      !!             is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of
-      !!             the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
+      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H.
+      !!                        - then we convert this flux to a "volume" flux this way => uH * uA / u
+      !!                             where uA is the flux from eq. a)
+      !!                             this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur)
+      !!                        - at last we estimate dV/dt = -div(uH * uA / u)
+      !!
+      !!             in eq. c), one can solve the equation for  S (ln_advS=T), then dVS/dt = -div(uV * uS  / u)
+      !!                                                or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u) 
+      !!
+      !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc.
+      !!           - At the ice edge, Ultimate scheme can lead to:
+      !!                              1) negative interpolated tracers at u-v points
+      !!                              2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward
+      !!                              Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of
+      !!                                               the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
+      !!                              Solution for 2): we set it to 0 in this case
       !!           - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good.
       !!             Large values of H can appear for very small ice concentration, and when it does it messes the things up since we
-      !!             work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D.
+      !!             work on H (and not V). It is partly related to the multi-category approach
       !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice
-      !!             concentration is small).
-      !! To-do: expand the prelimiter from zalesak to make it work in 2D
-      !!----------------------------------------------------------------------
-      REAL(wp)                        , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0)
-      INTEGER                         , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2)
-      INTEGER                         , INTENT(in   )           ::   jt             ! number of sub-iteration
-      INTEGER                         , INTENT(in   )           ::   kt             ! number of iteration
-      REAL(wp)                        , INTENT(in   )           ::   pdt            ! tracer time-step
-      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
-      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
-      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
-      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt             ! tracer field
-      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc            ! tracer content field
-      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes
+      !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter
+      !!             since sv_i and e_i are still good.
+      !!----------------------------------------------------------------------
+      REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0)
+      INTEGER                         , INTENT(in   )           ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
+      INTEGER                         , INTENT(in   )           ::   jt               ! number of sub-iteration
+      INTEGER                         , INTENT(in   )           ::   kt               ! number of iteration
+      REAL(wp)                        , INTENT(in   )           ::   pdt              ! tracer time-step
+      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv        ! 2 ice velocity components => u*e2
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc       ! 2 ice velocity components => u*e2 or u*a*e2u
+      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox     ! upstream velocity
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt               ! tracer field
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc              ! tracer content field
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL ::   pua_ups, pva_ups ! upstream u*a fluxes
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho   ! high order u*a fluxes
       !
       INTEGER  ::   ji, jj, jl       ! dummy loop indices  
@@ -303,7 +450,7 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
-                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
-                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
+                  IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
+                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj)
+                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj)
                   ELSE
                      zfu_ho (ji,jj,jl) = 0._wp
@@ -311,7 +458,7 @@
                   ENDIF
                   !
-                  IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
-                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
-                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
+                  IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
+                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj)
+                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj)
                   ELSE
                      zfv_ho (ji,jj,jl) = 0._wp  
@@ -321,15 +468,36 @@
             END DO
          END DO
+
+         ! the new "volume" fluxes must also be "flux corrected"
+         ! thus we calculate the upstream solution and apply a limiter again
+         DO jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
+                  !
+                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. )
+         !
+         IF    ( np_limiter == 1 ) THEN
+            CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
+         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
+            CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho )
+            CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho )
+         ENDIF
+         !
       ENDIF
-      !                                   --ho
-      ! in case of advection of A: output u*a
-      ! -------------------------------------
+      !                                   --ho    --ups
+      ! in case of advection of A: output u*a and u*a
+      ! -----------------------------------------------
       IF( PRESENT( pua_ho ) ) THEN
          DO jl = 1, jpl
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl)
-                  pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl)
-               END DO
+                  pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
+                  pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
+              END DO
             END DO
          END DO
@@ -499,7 +667,7 @@
          END DO
          !
-         IF    ( kn_limiter == 1 ) THEN
+         IF    ( np_limiter == 1 ) THEN
             CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
-         ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
+         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
             CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
             CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
@@ -517,5 +685,5 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
 
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
@@ -538,5 +706,5 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
 
          ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
@@ -549,5 +717,5 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
@@ -570,8 +738,8 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
 
          ENDIF
-         IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
          
       ENDIF
@@ -609,7 +777,7 @@
          !
          !                                                        !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho )
+         CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )
          !                                                        !--  limiter in x --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
@@ -619,5 +787,5 @@
                      &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) &
                      &                                                                                        * pamsk           &
-                     &                             ) * pdt ) * tmask(ji,jj,1) 
+                     &                             ) * pdt ) * tmask(ji,jj,1)
                END DO
             END DO
@@ -627,10 +795,10 @@
          !                                                        !--  ultimate interpolation of pt at v-point  --!
          IF( ll_hoxy ) THEN
-            CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
+            CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
          ELSE
-            CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho )
+            CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho )
          ENDIF
          !                                                        !--  limiter in y --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
          !         
          !
@@ -638,7 +806,7 @@
          !
          !                                                        !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho )
+         CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )
          !                                                        !--  limiter in y --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
@@ -656,19 +824,19 @@
          !                                                        !--  ultimate interpolation of pt at u-point  --!
          IF( ll_hoxy ) THEN
-            CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
+            CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
          ELSE
-            CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho )
+            CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho )
          ENDIF
          !                                                        !--  limiter in x --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
          !
       ENDIF
 
-      IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+      IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
       !
    END SUBROUTINE macho
 
 
-   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho )
+   SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -680,4 +848,5 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
+      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0)
       INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
       REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
@@ -806,5 +975,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( pt_u(ji,jj,jl) < 0._wp ) THEN
+                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                      pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
                         &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
@@ -826,5 +995,5 @@
    
  
-   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho )
+   SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_y  ***
@@ -836,4 +1005,5 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
+      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0)
       INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
       REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
@@ -959,5 +1129,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( pt_v(ji,jj,jl) < 0._wp ) THEN
+                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                      pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
                         &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
@@ -1023,5 +1193,5 @@
       !                        |      |      |        |    *
       !            t_ups :       i-1     i       i+1       i+2   
-      IF( ll_prelimiter_zalesak ) THEN
+      IF( ll_prelim ) THEN
          
          DO jl = 1, jpl
@@ -1102,10 +1272,11 @@
                !
                !                                  ! up & down beta terms
-               IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
-               ELSE                    ; zbetup(ji,jj,jl) = 0._wp ! zbig
+               ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
+               IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
+               ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig
                ENDIF
                !
-               IF( zneg > 0._wp ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
-               ELSE                    ; zbetdo(ji,jj,jl) = 0._wp ! zbig
+               IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
+               ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig
                ENDIF
                !
@@ -1149,18 +1320,4 @@
          END DO
 
-         ! clem test
-!!         DO jj = 2, jpjm1
-!!            DO ji = 2, fs_jpim1   ! vector opt.
-!!               zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
-!!                  &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
-!!                  &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
-!!                  &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
-!!                  &         ) * tmask(ji,jj,1)
-!!               IF( zzt < -epsi20 ) THEN
-!!                  WRITE(numout,*) 'T<0 nonosc_ice',zzt
-!!               ENDIF
-!!            END DO
-!!         END DO
-
       END DO
       !
@@ -1203,5 +1360,5 @@
                Rjp = zslpx(ji+1,jj,jl)
 
-               IF( kn_limiter == 3 ) THEN
+               IF( np_limiter == 3 ) THEN
 
                   IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
@@ -1219,5 +1376,5 @@
                   pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
 
-               ELSEIF( kn_limiter == 2 ) THEN
+               ELSEIF( np_limiter == 2 ) THEN
                   IF( Rj /= 0. ) THEN
                      IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
@@ -1298,5 +1455,5 @@
                Rjp = zslpy(ji,jj+1,jl)
 
-               IF( kn_limiter == 3 ) THEN
+               IF( np_limiter == 3 ) THEN
 
                   IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
@@ -1314,5 +1471,5 @@
                   pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
 
-               ELSEIF( kn_limiter == 2 ) THEN
+               ELSEIF( np_limiter == 2 ) THEN
 
                   IF( Rj /= 0. ) THEN
@@ -1358,4 +1515,94 @@
    END SUBROUTINE limiter_y
 
+
+   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE Hbig  ***
+      !!
+      !! ** Purpose : Thickness correction in case advection scheme creates
+      !!              abnormally tick ice or snow
+      !!
+      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
+      !!                 (before advection) and reduce it by adapting ice concentration
+      !!              2- check whether snow thickness is larger than the surrounding 9-points
+      !!                 (before advection) and reduce it by sending the excess in the ocean
+      !!              3- check whether snow load deplets the snow-ice interface below sea level$
+      !!                 and reduce it by sending the excess in the ocean
+      !!              4- correct pond fraction to avoid a_ip > a_i
+      !!
+      !! ** 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, psv_i, poa_i, pa_i, pa_ip, pv_ip
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i
+      !
+      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
+      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zvs_excess, zfra
+      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
+      !!-------------------------------------------------------------------
+      !
+      z1_dt = 1._wp / pdt
+      !
+      DO jl = 1, jpl
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+                  !
+                  !                               ! -- check h_ip -- !
+                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
+                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
+                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
+                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
+                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
+                     ENDIF
+                  ENDIF
+                  !
+                  !                               ! -- check h_i -- !
+                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
+                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
+                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                     pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
+                  ENDIF
+                  !
+                  !                               ! -- check h_s -- !
+                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
+                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
+                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
+                     !
+                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
+                     !
+                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
+                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
+                  ENDIF           
+                  !
+                  !                               ! -- check snow load -- !
+                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
+                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
+                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
+                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
+                  IF( zvs_excess > 0._wp ) THEN
+                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
+                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
+                     !
+                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
+                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
+                  ENDIF
+                  
+               ENDIF
+            END DO
+         END DO
+      END DO 
+      !                                           !-- correct pond fraction to avoid a_ip > a_i
+      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
+      !
+      !
+   END SUBROUTINE Hbig
+   
 #else
    !!----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_rdgrft.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_rdgrft.F90	(revision 11082)
@@ -142,8 +142,4 @@
       INTEGER, PARAMETER ::   jp_itermax = 20    
       !!-------------------------------------------------------------------
-      ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem)
-      !       likely due to truncation error ( i.e. 1. - 1. /= 0 )
-      !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)
-      
       ! controls
       IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing
@@ -156,14 +152,14 @@
       ENDIF      
 
-      CALL ice_var_zapsmall   ! Zero out categories with very small areas
-
       !--------------------------------
       ! 0) Identify grid cells with ice
       !--------------------------------
+      at_i(:,:) = SUM( a_i, dim=3 )
+      !
       npti = 0   ;   nptidx(:) = 0
       ipti = 0   ;   iptidx(:) = 0
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF ( at_i(ji,jj) > 0._wp ) THEN
+            IF ( at_i(ji,jj) > epsi10 ) THEN
                npti           = npti + 1
                nptidx( npti ) = (jj - 1) * jpi + ji
@@ -178,10 +174,10 @@
         
          ! 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(:,:) )
+         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_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(:,:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_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   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i )
 
          DO ji = 1, npti
@@ -310,6 +306,6 @@
 
       !                       ! Ice thickness needed for rafting
-      WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
-      ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp
+      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
+      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
       END WHERE
 
@@ -329,6 +325,6 @@
       zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
       !
-      WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
-      ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp
+      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
+      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp
       END WHERE
       !
@@ -455,6 +451,6 @@
       ! 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) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
-      ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp
+      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
+      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
       END WHERE
       
@@ -466,5 +462,5 @@
          DO ji = 1, npti
             zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
-            IF( zfac > pa_i(ji,jl) ) THEN
+            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
                closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice
             ENDIF
@@ -510,4 +506,5 @@
       REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
       REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a
+      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i)
       !
       REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice
@@ -518,5 +515,7 @@
       INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation
       !!-------------------------------------------------------------------
-
+      !
+      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume
+      !
       ! 1) Change in open water area due to closing and opening
       !--------------------------------------------------------
@@ -535,6 +534,8 @@
             IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
 
-               z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
-
+               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
+               ELSE                                 ;   z1_ai(ji) = 0._wp
+               ENDIF
+               
                ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
                airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
@@ -549,10 +550,13 @@
 
                ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
-               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg
+               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg
+               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0
+               ELSE                           ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0
+               ENDIF
                ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
 
                ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
                virdg1     = v_i_2d (ji,jl1)   * afrdg
-               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg )
+               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw
                vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg
                sirdg1     = sv_i_2d(ji,jl1)   * afrdg
@@ -726,7 +730,8 @@
       END DO ! jl1
       !
+      ! roundoff errors
+      !----------------
       ! In case ridging/rafting lead to very small negative values (sometimes it happens)
-      WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp
-      WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp
+      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )
       !
    END SUBROUTINE rdgrft_shift
@@ -854,5 +859,5 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
-
+         !
          !                 !---------------------!
       CASE( 2 )            !==  from 1D to 2D  ==!
@@ -945,5 +950,19 @@
          CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
       ENDIF
-      !                              ! allocate tke arrays
+      !
+      IF( .NOT. ln_icethd ) THEN
+         rn_porordg = 0._wp
+         rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp
+         rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp
+         IF( lwp ) THEN
+            WRITE(numout,*) '      ==> only ice dynamics is activated, thus some parameters must be changed'
+            WRITE(numout,*) '            rn_porordg   = ', rn_porordg
+            WRITE(numout,*) '            rn_fsnwrdg   = ', rn_fsnwrdg 
+            WRITE(numout,*) '            rn_fpndrdg   = ', rn_fpndrdg 
+            WRITE(numout,*) '            rn_fsnwrft   = ', rn_fsnwrft 
+            WRITE(numout,*) '            rn_fpndrft   = ', rn_fpndrft 
+         ENDIF
+      ENDIF
+      !                              ! allocate arrays
       IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
       !
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_rhg.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_rhg.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_rhg.F90	(revision 11082)
@@ -58,4 +58,6 @@
       !!--------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! ice time step
+      !
+      INTEGER  ::   jl   ! dummy loop indices
       !!--------------------------------------------------------------------
       ! controls
@@ -68,8 +70,15 @@
          WRITE(numout,*)'~~~~~~~~~~~'
       ENDIF
-
-      ! --------
-      ! Rheology
-      ! --------   
+      !
+      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
+         tau_icebfr(:,:) = 0._wp
+         DO jl = 1, jpl
+            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
+         END DO
+      ENDIF
+      !
+      !--------------!
+      !== Rheology ==!
+      !--------------!   
       SELECT CASE( nice_rhg )
       !                                !------------------------!
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_rhg_evp.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icedyn_rhg_evp.F90	(revision 11082)
@@ -112,5 +112,5 @@
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
       !!
-      LOGICAL, PARAMETER ::   ll_bdy_substep = .FALSE. ! temporary option to call bdy at each sub-time step (T)
+      LOGICAL, PARAMETER ::   ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T)
       !                                                                              or only at the main time step (F)
       INTEGER ::   ji, jj       ! dummy loop indices
@@ -160,5 +160,6 @@
 
       REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
-      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
+      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small
+      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
       !! --- diags
       REAL(wp), DIMENSION(jpi,jpj) ::   zswi
@@ -319,6 +320,8 @@
 
             ! switches
-            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
-            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
+            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zswitchU(ji,jj) = 0._wp
+            ELSE                                                   ;   zswitchU(ji,jj) = 1._wp   ;   ENDIF
+            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zswitchV(ji,jj) = 0._wp
+            ELSE                                                   ;   zswitchV(ji,jj) = 1._wp   ;   ENDIF
 
          END DO
@@ -519,5 +522,5 @@
                      &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
                      &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
-                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           ) * zmaskV(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
@@ -526,5 +529,5 @@
                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
                      &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
+                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &            ) * zmaskV(ji,jj)
                   ENDIF
@@ -567,5 +570,5 @@
                      &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
-                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin 
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
                      &            ) * zmaskU(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
@@ -574,5 +577,5 @@
                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
                      &            ) * zmaskU(ji,jj)
                   ENDIF
@@ -617,5 +620,5 @@
                      &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
-                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin 
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
                      &            ) * zmaskU(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
@@ -624,5 +627,5 @@
                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &            ) * zmaskU(ji,jj)
                   ENDIF
@@ -665,5 +668,5 @@
                      &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
                      &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
-                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           ) * zmaskV(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
@@ -672,5 +675,5 @@
                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
                      &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
+                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &            ) * zmaskV(ji,jj)
                   ENDIF
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/iceitd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/iceitd.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/iceitd.F90	(revision 11082)
@@ -21,4 +21,5 @@
    USE ice1D          ! sea-ice: thermodynamic variables
    USE ice            ! sea-ice: variables
+   USE icevar         ! sea-ice: operations
    USE icectl         ! sea-ice: conservation tests
    USE icetab         ! sea-ice: convert 1D<=>2D
@@ -91,4 +92,6 @@
       !  1) Identify grid cells with ice
       !-----------------------------------------------------------------------------------------------
+      at_i(:,:) = SUM( a_i, dim=3 )
+      !
       npti = 0   ;   nptidx(:) = 0
       DO jj = 1, jpj
@@ -249,5 +252,5 @@
             ! --- g(h) for each thickness category --- !  
             CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in
-               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR   (1:npti,jl)   )   ! out
+               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL    (1:npti,jl), hR    (1:npti,jl)   )   ! out
             !
          END DO
@@ -389,9 +392,11 @@
       REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary
       !
-      INTEGER  ::   ji, jj, jl, jk     ! dummy loop indices
-      INTEGER  ::   ii, ij, jl2, jl1   ! local integers
+      INTEGER  ::   ji, jl, jk         ! dummy loop indices
+      INTEGER  ::   jl2, jl1           ! local integers
       REAL(wp) ::   ztrans             ! ice/snow transferred
-      REAL(wp), DIMENSION(jpij)     ::   zworka, zworkv   ! workspace
-      REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn           !  -    -
+      REAL(wp), DIMENSION(jpij)            ::   zworka, zworkv   ! workspace
+      REAL(wp), DIMENSION(jpij,jpl)        ::   zaTsfn           !  -    -
+      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d
+      REAL(wp), DIMENSION(jpij,nlay_s,jpl) ::   ze_s_2d
       !!------------------------------------------------------------------
          
@@ -405,4 +410,14 @@
       CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
       CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
+      DO jl = 1, jpl
+         DO jk = 1, nlay_s
+            CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
+         END DO
+         DO jk = 1, nlay_i
+            CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+         END DO
+      END DO
+      ! to correct roundoff errors on a_i
+      CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d )
 
       !----------------------------------------------------------------------------------------------
@@ -435,8 +450,4 @@
                ELSE                                  ;   zworka(ji) = 0._wp
                ENDIF
-               !
-               ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
-               !       because of truncation error ( i.e. 1. - 1. /= 0 )
-               !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)
                !
                a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
@@ -476,8 +487,5 @@
          !
          DO jk = 1, nlay_s         !--- Snow heat content
-            !
             DO ji = 1, npti
-               ii = MOD( nptidx(ji) - 1, jpi ) + 1
-               ij = ( nptidx(ji) - 1 ) / jpi + 1
                !
                jl1 = kdonor(ji,jl)
@@ -487,8 +495,7 @@
                   ELSE                ;  jl2 = jl
                   ENDIF
-                  !
-                  ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
-                  e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
-                  e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
+                  ztrans             = ze_s_2d(ji,jk,jl1) * zworkv(ji)
+                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans
+                  ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans
                ENDIF
             END DO
@@ -497,6 +504,4 @@
          DO jk = 1, nlay_i         !--- Ice heat content
             DO ji = 1, npti
-               ii = MOD( nptidx(ji) - 1, jpi ) + 1
-               ij = ( nptidx(ji) - 1 ) / jpi + 1
                !
                jl1 = kdonor(ji,jl)
@@ -506,8 +511,7 @@
                   ELSE                ;  jl2 = jl
                   ENDIF
-                  !
-                  ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
-                  e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
-                  e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
+                  ztrans             = ze_i_2d(ji,jk,jl1) * zworkv(ji)
+                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans
+                  ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans
                ENDIF
             END DO
@@ -515,7 +519,21 @@
          !
       END DO                   ! boundaries, 1 to jpl-1
+
+      !-------------------
+      ! 3) roundoff errors
+      !-------------------
+      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
+      !       because of truncation error ( i.e. 1. - 1. /= 0 )
+      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )
+
+      ! at_i must be <= rn_amax
+      zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
+      DO jl  = 1, jpl
+         WHERE( zworka(1:npti) > rn_amax_1d(1:npti) )   &
+            &   a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti)
+      END DO
       
       !-------------------------------------------------------------------------------
-      ! 3) Update ice thickness and temperature
+      ! 4) Update ice thickness and temperature
       !-------------------------------------------------------------------------------
       WHERE( a_i_2d(1:npti,:) >= epsi20 )
@@ -536,4 +554,12 @@
       CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
       CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
+      DO jl = 1, jpl
+         DO jk = 1, nlay_s
+            CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
+         END DO
+         DO jk = 1, nlay_i
+            CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+         END DO
+      END DO
       !
    END SUBROUTINE itd_shiftice
@@ -558,4 +584,6 @@
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
+      !
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
       !
       jdonor(:,:) = 0
@@ -635,4 +663,6 @@
       END DO
       !
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+      !
    END SUBROUTINE ice_itd_reb
 
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icestp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icestp.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icestp.F90	(revision 11082)
@@ -189,5 +189,5 @@
          IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics      
          !
-         IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections
+                                        CALL ice_cor( kt , 2 )        ! -- Corrections
          !
                                         CALL ice_var_glo2eqv          ! necessary calls (at least for coupling)
@@ -323,4 +323,7 @@
          WRITE(numout,*) '         maximum ice concentration for SH                              = ', rn_amax_s
       ENDIF
+      !                                        !--- change max ice concentration for roundoff errors
+      rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 )
+      rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 )
       !                                        !--- check consistency
       IF ( jpl > 1 .AND. ln_virtual_itd ) THEN
@@ -431,7 +434,9 @@
       t_si       (:,:,:) = rt0   ! temp at the ice-snow interface
 
-      tau_icebfr(:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here)
-      cnd_ice   (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T)
-      qtr_ice_bot(:,:,:) = 0._wp  ! initialization: part of solar radiation transmitted through the ice needed at least for outputs
+      tau_icebfr (:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here)
+      cnd_ice    (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T)
+      qcn_ice    (:,:,:) = 0._wp   ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T)
+      qtr_ice_bot(:,:,:) = 0._wp   ! initialization: part of solar radiation transmitted through the ice needed at least for outputs
+      qsb_ice_bot(:,:)   = 0._wp   ! (needed if ln_icethd=F)
       !
       ! for control checks (ln_icediachk)
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icethd.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icethd.F90	(revision 11082)
@@ -102,6 +102,4 @@
       ENDIF
       
-      CALL ice_var_glo2eqv
-
       !---------------------------------------------!
       ! computation of friction velocity at T points
@@ -162,6 +160,7 @@
             qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
 
-            ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 
-            IF( zqld > 0._wp ) THEN
+            ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 
+            ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget
+            IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN
                fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
                qlead(ji,jj) = 0._wp
@@ -178,7 +177,6 @@
       ! In case we bypass open-water ice formation
       IF( .NOT. ln_icedO )  qlead(:,:) = 0._wp
-      ! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
+      ! In case we bypass growing/melting from top and bottom
       IF( .NOT. ln_icedH ) THEN
-         qt_atm_oi  (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
          qsb_ice_bot(:,:) = 0._wp
          fhld       (:,:) = 0._wp
@@ -221,36 +219,32 @@
             dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp
             dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
-            !
-            IF( ln_icedH ) THEN                                     ! --- growing/melting --- !
-                              CALL ice_thd_zdf                             ! Ice/Snow Temperature profile
-                              CALL ice_thd_dh                              ! Ice/Snow thickness   
-                              CALL ice_thd_pnd                             ! Melt ponds formation
-                              CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping
+            !                                      
+                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- !
+            !
+            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- !
+                              CALL ice_thd_dh                           ! Ice-Snow thickness   
+                              CALL ice_thd_pnd                          ! Melt ponds formation
+                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping
             ENDIF
-            !
                               CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !    
             !
-                              CALL ice_thd_temp                     ! --- temperature update --- !
+                              CALL ice_thd_temp                     ! --- Temperature update --- !
             !
             IF( ln_icedH .AND. ln_virtual_itd ) &
-               &              CALL ice_thd_mono                     ! --- extra lateral melting if virtual_itd --- !
-            !
-            IF( ln_icedA )    CALL ice_thd_da                       ! --- lateral melting --- !
+               &              CALL ice_thd_mono                     ! --- Extra lateral melting if virtual_itd --- !
+            !
+            IF( ln_icedA )    CALL ice_thd_da                       ! --- Lateral melting --- !
             !
                               CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
             !                                                       ! --- & Move to 2D arrays --- !
-            !
          ENDIF
          !
       END DO
-
+      !
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
-      !
-                           CALL ice_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- !
-      !                                                    !     & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0
       !                   
-      IF( jpl > 1      )   CALL ice_itd_rem( kt )          ! --- Transport ice between thickness categories --- !
-      !
-      IF( ln_icedO     )   CALL ice_thd_do                 ! --- frazil ice growing in leads --- !
+      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- !
+      !
+      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- !
       !
       ! controls
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icethd_do.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icethd_do.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icethd_do.F90	(revision 11082)
@@ -114,7 +114,5 @@
       IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
 
-      CALL ice_var_agg(1)
-      CALL ice_var_glo2eqv
-
+      at_i(:,:) = SUM( a_i, dim=3 )
       !------------------------------------------------------------------------------!
       ! 1) Collection thickness of ice formed in leads and polynyas
@@ -319,12 +317,12 @@
 
          ! --- lateral ice growth --- !
-         ! If lateral ice growth gives an ice concentration gt 1, then
+         ! If lateral ice growth gives an ice concentration > amax, then
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, npti
-            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) )
+            IF ( za_newice(ji) >  MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
+               zda_res(ji)   = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
-               za_newice(ji) = za_newice(ji) - zda_res  (ji)
-               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji)
+               za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res  (ji) )
+               zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res  (ji) )
             ELSE
                zda_res(ji) = 0._wp
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icethd_zdf_bl99.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icethd_zdf_bl99.F90	(revision 11082)
@@ -206,11 +206,11 @@
       !
       l_T_converged(:) = .FALSE.
-      !                                                          !============================!
       ! Convergence calculated until all sub-domain grid points have converged
       ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation)
       ! but values are not taken into account (results independant of MPI partitioning)
       !
+      !                                                                            !============================!
       DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max )   ! Iterative procedure begins !
-         !                                                       !============================!
+         !                                                                         !============================!
          iconv = iconv + 1
          !
@@ -742,5 +742,5 @@
                   zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
                   ! t_i
-                  DO jk = 0, nlay_i
+                  DO jk = 1, nlay_i
                      ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0 
                      t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
@@ -856,4 +856,8 @@
          t_i_1d    (1:npti,:) = ztiold        (1:npti,:)
          qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti)
+
+         !!clem
+         ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input
+         !clem
       ENDIF
       !
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icevar.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icevar.F90	(revision 11082)
@@ -44,5 +44,6 @@
    !!   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_zapneg    : remove negative ice fields
+   !!   ice_var_roundoff  : remove negative values arising from roundoff erros
    !!   ice_var_itd       : convert 1-cat to jpl-cat
    !!   ice_var_itd2      : convert N-cat to jpl-cat
@@ -71,4 +72,5 @@
    PUBLIC   ice_var_zapsmall
    PUBLIC   ice_var_zapneg
+   PUBLIC   ice_var_roundoff
    PUBLIC   ice_var_itd
    PUBLIC   ice_var_itd2
@@ -229,5 +231,5 @@
                   IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
                      !
-                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
+                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i             ! Energy of melting e(S,T) [J.m-3]
                      ztmelts          = - sz_i(ji,jj,jk,jl) * rTmlt                                 ! Ice layer melt temperature [C]
                      ! Conversion q(S,T) -> T (second order equation)
@@ -236,5 +238,5 @@
                      t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
                      !
-                  ELSE                                !--- no ice
+                  ELSE                                   !--- no ice
                      t_i(ji,jj,jk,jl) = rt0
                   ENDIF
@@ -537,5 +539,5 @@
 
 
-   SUBROUTINE ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_zapneg ***
@@ -543,6 +545,5 @@
       !! ** Purpose :   Remove negative sea ice fields and correct fluxes
       !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
-      !
+      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
       REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
@@ -555,9 +556,16 @@
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
-      !!-------------------------------------------------------------------
-      !
+      !
+      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
+      REAL(wp) ::   z1_dt
+      !!-------------------------------------------------------------------
+      !
+      z1_dt = 1._wp / pdt
       !
       DO jl = 1, jpl       !==  loop over the categories  ==!
          !
+         ! make sure a_i=0 where v_i<=0
+         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp
+
          !----------------------------------------
          ! zap ice energy and send it to the ocean
@@ -566,6 +574,6 @@
             DO jj = 1 , jpj
                DO ji = 1 , jpi
-                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
-                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 >0
+                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
                      pe_i(ji,jj,jk,jl) = 0._wp
                   ENDIF
@@ -577,6 +585,6 @@
             DO jj = 1 , jpj
                DO ji = 1 , jpi
-                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
-                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
                      pe_s(ji,jj,jk,jl) = 0._wp
                   ENDIF
@@ -590,14 +598,14 @@
          DO jj = 1 , jpj
             DO ji = 1 , jpi
-               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
-                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice
+               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
                   pv_i   (ji,jj,jl) = 0._wp
                ENDIF
-               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
-                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice
+               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
                   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
-                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice
+               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_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
                ENDIF
@@ -616,4 +624,36 @@
    END SUBROUTINE ice_var_zapneg
 
+
+   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_var_roundoff ***
+      !!
+      !! ** Purpose :   Remove negative sea ice values arising from roundoff errors
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
+      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_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( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
+      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
+      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
+      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
+      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
+      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
+      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
+      IF ( ln_pnd_H12 ) THEN
+         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
+         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
+      ENDIF
+      !
+   END SUBROUTINE ice_var_roundoff
+   
    
    SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
@@ -650,6 +690,6 @@
       INTEGER  :: idim, i_fill, jl0  
       REAL(wp) :: zarg, zV, zconv, zdh, zdv
-      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
+      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
       INTEGER , DIMENSION(4)                  ::   itest
       !!-------------------------------------------------------------------
@@ -780,5 +820,5 @@
       !!
       !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
-      !!                   by removing 25% ice area from jlmin and jlmax (resp.) 
+       !!                   by removing 25% ice area from jlmin and jlmax (resp.) 
       !!              
       !!               3) Expand the filling to the empty cat between jlmin and jlmax 
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icewri.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icewri.F90	(revision 11082)
@@ -145,5 +145,5 @@
 
         ! record presence of fast ice
-         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp
+         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp
          ELSEWHERE                                                ; zfast(:,:) = 0._wp
          END WHERE
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/NST/agrif_top_update.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/NST/agrif_top_update.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/NST/agrif_top_update.F90	(revision 11082)
@@ -138,5 +138,5 @@
                      DO ji=i1,i2
                         IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
-                           trb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
+                           trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn) & 
                                  & + atfp * ( tabres_child(ji,jj,jk,jn) &
                                  &          - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdy_oce.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdy_oce.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdy_oce.F90	(revision 11082)
@@ -85,5 +85,5 @@
    !
    INTEGER                    ::   nb_bdy                   !: number of open boundary sets
-   INTEGER                    ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run)
+   INTEGER, DIMENSION(jp_bdy) ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run)
    INTEGER, DIMENSION(jp_bdy) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme
    INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P 
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdydta.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdydta.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdydta.F90	(revision 11082)
@@ -243,5 +243,5 @@
                         IF( ln_full_vel_array(jbdy) ) THEN
                            CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
-                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy,   &
+                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy(jbdy),   &
                                      & fvl=ln_full_vel_array(jbdy)  )
                         ELSE
@@ -313,5 +313,5 @@
                      jend = jstart + dta%nread(1) - 1
                      CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), &
-                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy,   &
+                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy(jbdy),   &
                                   & fvl=ln_full_vel_array(jbdy) )
                   ENDIF
@@ -446,5 +446,5 @@
       NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s
 #endif
-      NAMELIST/nambdy_dta/ ln_full_vel, nb_jpk_bdy
+      NAMELIST/nambdy_dta/ ln_full_vel
       !!---------------------------------------------------------------------------
       !
@@ -508,9 +508,9 @@
       ! Read namelists
       ! --------------
-      REWIND(numnam_ref)
       REWIND(numnam_cfg)
       jfld = 0 
       DO jbdy = 1, nb_bdy         
          IF( nn_dta(jbdy) == 1 ) THEN
+            REWIND(numnam_ref)
             READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
 901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdydyn2d.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdydyn2d.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdydyn2d.F90	(revision 11082)
@@ -187,5 +187,5 @@
          ! Use characteristics method instead
          zflag = ABS(flagu)
-         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
+         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(ii+NINT(flagu),ij)
          pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
       END DO
@@ -205,5 +205,5 @@
          ! Use characteristics method instead
          zflag = ABS(flagv)
-         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
+         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ij+NINT(flagv))
          pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
       END DO
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdyice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdyice.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdyice.F90	(revision 11082)
@@ -57,6 +57,7 @@
       INTEGER ::   jbdy   ! BDY set index
       !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('bdy_ice_thd')
+      ! controls
+      IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing
+      IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       !
       CALL ice_var_glo2eqv
@@ -78,6 +79,8 @@
       CALL ice_var_agg(1)
       !
-      IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
-      IF( ln_timing )   CALL timing_stop('bdy_ice_thd')
+      ! controls
+      IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints
+      IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing
       !
    END SUBROUTINE bdy_ice
@@ -148,8 +151,8 @@
             jpbound = 0   ;   ib = ji   ;   jb = jj
             !
-            IF( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji+1 ; jb = jj
-            IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji-1 ; jb = jj
-            IF( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. )   jpbound = 1 ; ib = ji   ; jb = jj+1
-            IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. )   jpbound = 1 ; ib = ji   ; jb = jj-1
+            IF( u_ice(ji  ,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji+1
+            IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji  ,jj  ,1) == 0. )   jpbound = 1 ; ib = ji-1
+            IF( v_ice(ji  ,jj  ) < 0. .AND. vmask(ji  ,jj-1,1) == 0. )   jpbound = 1 ; jb = jj+1
+            IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj  ,1) == 0. )   jpbound = 1 ; jb = jj-1
             !
             IF( nn_ice_dta(jbdy) == 0 )   jpbound = 0 ; ib = ji ; jb = jj   ! case ice boundaries = initial conditions
@@ -306,5 +309,5 @@
                      ! one of the two zmsk is always 0 (because of zflag)
                      zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
-                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
+                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) )   ! 0 if no ice
                      !  
                      ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice)
@@ -329,5 +332,5 @@
                      ! one of the two zmsk is always 0 (because of zflag)
                      zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
-                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
+                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) )   ! 0 if no ice
                      !  
                      ! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice)
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdyini.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdyini.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/BDY/bdyini.F90	(revision 11082)
@@ -140,5 +140,5 @@
       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
       CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid
-      INTEGER :: com_east, com_west, com_south, com_north          ! Flags for boundaries sending
+      INTEGER :: com_east, com_west, com_south, com_north, jpk_max  ! Flags for boundaries sending
       INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving
       INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates
@@ -397,5 +397,5 @@
           IF(lwp) WRITE(numout,*)
         ENDIF
-        IF( nb_jpk_bdy > 0 ) THEN
+        IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN
            IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***'
         ELSE
@@ -516,25 +516,16 @@
          ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
             &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
-
-         IF( nb_jpk_bdy>0 ) THEN
-            ALLOCATE( dta_global(jpbdtau, 1, nb_jpk_bdy) )
-            ALLOCATE( dta_global_z(jpbdtau, 1, nb_jpk_bdy) )
-            ALLOCATE( dta_global_dz(jpbdtau, 1, nb_jpk_bdy) )
-         ELSE
-            ALLOCATE( dta_global(jpbdtau, 1, jpk) )
-            ALLOCATE( dta_global_z(jpbdtau, 1, jpk) ) ! needed ?? TODO
-            ALLOCATE( dta_global_dz(jpbdtau, 1, jpk) )! needed ?? TODO
-         ENDIF
+         
+         jpk_max = MAXVAL(nb_jpk_bdy)
+         jpk_max = MAX(jpk_max, jpk)
+
+         ALLOCATE( dta_global(jpbdtau, 1, jpk_max) )
+         ALLOCATE( dta_global_z(jpbdtau, 1, jpk_max) ) ! needed ?? TODO
+         ALLOCATE( dta_global_dz(jpbdtau, 1, jpk_max) )! needed ?? TODO
 
          IF ( icount>0 ) THEN
-            IF( nb_jpk_bdy>0 ) THEN
-               ALLOCATE( dta_global2(jpbdtas, nrimmax, nb_jpk_bdy) )
-               ALLOCATE( dta_global2_z(jpbdtas, nrimmax, nb_jpk_bdy) )
-               ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, nb_jpk_bdy) )
-            ELSE
-               ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
-               ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk) ) ! needed ?? TODO
-               ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk) )! needed ?? TODO  
-            ENDIF
+            ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk_max) )
+            ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk_max) ) ! needed ?? TODO
+            ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk_max) )! needed ?? TODO  
          ENDIF
          ! 
@@ -960,5 +951,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
-                       if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
+                       if( ii == (nlcit(nowe+1)-1) ) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
@@ -974,5 +965,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
-                       if((com_east_b .ne. 1) .and. (ii == 2)) then
+                       if( ii == 2 ) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
@@ -989,5 +980,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
-                       if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
+                       if( ii == (nlcit(nowe+1)-1) ) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
@@ -1004,5 +995,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
-                       if((com_east_b .ne. 1) .and. (ii == 2)) then
+                       if( ii == 2 ) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DIA/diacfl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DIA/diacfl.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DIA/diacfl.F90	(revision 11082)
@@ -29,9 +29,4 @@
    REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number 
 
-!!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc !
-!!gm          I don't understand why.
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
-!!gm end
-
    PUBLIC   dia_cfl       ! routine called by step.F90
    PUBLIC   dia_cfl_init  ! routine called by nemogcm
@@ -55,8 +50,8 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER                ::   ji, jj, jk                            ! dummy loop indices
-      REAL(wp)               ::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars
-      INTEGER , DIMENSION(3) ::   iloc_u , iloc_v , iloc_w , iloc       ! workspace
-!!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
+      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices
+      REAL(wp)                         ::   z2dt, zCu_max, zCv_max, zCw_max  ! local scalars
+      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
       !!----------------------------------------------------------------------
       !
@@ -71,5 +66,5 @@
       DO jk = 1, jpk       ! calculate Courant numbers
          DO jj = 1, jpj
-            DO ji = 1, fs_jpim1   ! vector opt.
+            DO ji = 1, jpi
                zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction
                zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction
@@ -111,5 +106,5 @@
       !                    ! write out to file
       IF( lwp ) THEN
-         WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
+         WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
          WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
          WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
@@ -172,8 +167,4 @@
       rCw_max = 0._wp
       !
-!!gm required to work
-      ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) )
-!!gm end
-      !      
    END SUBROUTINE dia_cfl_init
 
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DYN/dynkeg.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DYN/dynkeg.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DYN/dynkeg.F90	(revision 11082)
@@ -74,10 +74,10 @@
       INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme 
       !
-      INTEGER  ::   ji, jj, jk, jb    ! dummy loop indices
-      INTEGER  ::   ii, ifu, ib_bdy   ! local integers
-      INTEGER  ::   ij, ifv, igrd     !   -       -
-      REAL(wp) ::   zu, zv            ! local scalars
+      INTEGER  ::   ji, jj, jk, jb           ! dummy loop indices
+      INTEGER  ::   ifu, ifv, igrd, ib_bdy   ! local integers
+      REAL(wp) ::   zu, zv                   ! local scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
+      REAL(wp)  :: zweightu, zweightv
       !!----------------------------------------------------------------------
       !
@@ -97,32 +97,4 @@
       
       zhke(:,:,jpk) = 0._wp
-      
-      IF (ln_bdy) THEN
-         ! Maria Luneva & Fred Wobus: July-2016
-         ! compensate for lack of turbulent kinetic energy on liquid bdy points
-         DO ib_bdy = 1, nb_bdy
-            IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
-               igrd = 2           ! Copying normal velocity into points outside bdy
-               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-                  DO jk = 1, jpkm1
-                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
-                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
-                     ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
-                     un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
-                  END DO
-               END DO
-               !
-               igrd = 3           ! Copying normal velocity into points outside bdy
-               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-                  DO jk = 1, jpkm1
-                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
-                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
-                     ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
-                     vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
-                  END DO
-               END DO
-            ENDIF
-         ENDDO  
-      ENDIF 
 
       SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
@@ -140,4 +112,28 @@
             END DO
          END DO
+         !
+         IF (ln_bdy) THEN
+            ! Maria Luneva & Fred Wobus: July-2016
+            ! compensate for lack of turbulent kinetic energy on liquid bdy points
+            DO ib_bdy = 1, nb_bdy
+               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
+                  igrd = 1           ! compensating null velocity on the bdy
+                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1
+                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1
+                     DO jk = 1, jpkm1
+                        zhke(ji,jj,jk) = 0._wp
+                        zweightu = umask(ji-1,jj  ,jk) + umask(ji,jj,jk)
+                        zweightv = vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)
+                        zu = un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)  +  un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
+                        zv = vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  +  vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
+                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zu / (2._wp * zweightu) 
+                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zv / (2._wp * zweightv) 
+                     END DO
+                  END DO
+               END IF
+               CALL lbc_bdy_lnk( 'dynkeg', zhke, 'T', 1., ib_bdy )   ! send 2 and recv jpi, jpj used in the computation of the speed tendencies
+            END DO
+         END IF
          !
       CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
@@ -158,14 +154,37 @@
             END DO
          END DO
+         IF (ln_bdy) THEN
+            ! Maria Luneva & Fred Wobus: July-2016
+            ! compensate for lack of turbulent kinetic energy on liquid bdy points
+            DO ib_bdy = 1, nb_bdy
+               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
+                  igrd = 1           ! compensation null velocity on land at the bdy
+                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1
+                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1
+                     DO jk = 1, jpkm1
+                        zhke(ji,jj,jk) = 0._wp
+                        zweightu = 8._wp * ( umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) ) &
+                             &   + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji  ,jj-1,jk) + umask(ji  ,jj+1,jk) ) 
+                        zweightv = 8._wp * ( vmask(ji  ,jj-1,jk) + vmask(ji  ,jj-1,jk) ) &
+                             &   + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj  ,jk) + vmask(ji+1,jj  ,jk) )
+                        zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
+                           &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
+                           &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
+                           &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
+                        zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
+                           &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
+                           &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
+                           &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
+                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zu / ( 2._wp * zweightu )
+                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zv / ( 2._wp * zweightv )
+                     END DO
+                  END DO
+               END IF
+            END DO
+         END IF
          CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
          !
-      END SELECT
-
-      IF (ln_bdy) THEN
-         ! restore velocity masks at points outside boundary
-         un(:,:,:) = un(:,:,:) * umask(:,:,:)
-         vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
-      ENDIF      
-
+      END SELECT 
       !
       DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DYN/sshwzv.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DYN/sshwzv.F90	(revision 11082)
@@ -297,7 +297,6 @@
          IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
          IF(lwp) WRITE(numout,*) '~~~~~ '
-         !
-         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all)
-      ENDIF
+      ENDIF
+      !
       !
       DO jk = 1, jpkm1            ! calculate Courant numbers
@@ -305,27 +304,27 @@
             DO ji = 2, fs_jpim1   ! vector opt.
                z1_e3w = 1._wp / e3w_n(ji,jj,jk)
-               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &
-                  &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
-                  &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
-                  &                        * r1_e1e2t(ji,jj)                                                  &
-                  &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
-                  &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
-                  &                        * r1_e1e2t(ji,jj)                                                  &
-                  &                      ) * z1_e3w
+               Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &  ! 2*rdt and not r2dt (for restartability)
+                  &                             + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
+                  &                                 MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
+                  &                               * r1_e1e2t(ji,jj)                                                  &
+                  &                             + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
+                  &                                 MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
+                  &                               * r1_e1e2t(ji,jj)                                                  &
+                  &                             ) * z1_e3w
             END DO
          END DO
       END DO
+      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
       !
       CALL iom_put("Courant",Cu_adv)
       !
-      wi(:,:,:) = 0._wp                                 ! Includes top and bottom values set to zero
       IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
          DO jk = 1, jpkm1                               ! or scan Courant criterion and partition
-            DO jj = 2, jpjm1                            ! w where necessary
-               DO ji = 2, fs_jpim1   ! vector opt.
+            DO jj = 1, jpj                              ! w where necessary
+               DO ji = 1, jpi
                   !
                   zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
                   !
-                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit
+                  IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
                      zcff = 0._wp
                   ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
@@ -346,5 +345,6 @@
       ELSE
          ! Fully explicit everywhere
-         Cu_adv = 0.0_wp                                ! Reuse array to output coefficient
+         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient
+         wi    (:,:,:) = 0._wp
       ENDIF
       CALL iom_put("wimp",wi) 
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/LBC/lib_mpp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/LBC/lib_mpp.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/LBC/lib_mpp.F90	(revision 11082)
@@ -1480,6 +1480,7 @@
       LOGICAL         , OPTIONAL, INTENT(in   ) ::   ld_lbc, ld_glb, ld_dlg
       !!
+      CHARACTER(len=128)                        ::   ccountname  ! name of a subroutine to count communications
       LOGICAL ::   ll_lbc, ll_glb, ll_dlg
-      INTEGER ::    ji,  jj,  jk,  jh, jf   ! dummy loop indices
+      INTEGER ::    ji,  jj,  jk,  jh, jf, jcount   ! dummy loop indices
       !!----------------------------------------------------------------------
       !
@@ -1538,13 +1539,21 @@
          WRITE(numcom,*) ' '
          WRITE(numcom,*) ' lbc_lnk called'
-         jj = 1
-         DO ji = 2, n_sequence_lbc
-            IF( crname_lbc(ji-1) /= crname_lbc(ji) ) THEN
-               WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_lbc(ji-1))
-               jj = 0
+         DO ji = 1, n_sequence_lbc - 1
+            IF ( crname_lbc(ji) /= 'already counted' ) THEN
+               ccountname = crname_lbc(ji)
+               crname_lbc(ji) = 'already counted'
+               jcount = 1
+               DO jj = ji + 1, n_sequence_lbc
+                  IF ( ccountname ==  crname_lbc(jj) ) THEN
+                     jcount = jcount + 1
+                     crname_lbc(jj) = 'already counted'
+                  END IF
+               END DO
+               WRITE(numcom,'(A, I4, A, A)') ' - ', jcount,' times by subroutine ', TRIM(ccountname)
             END IF
-            jj = jj + 1 
          END DO
-         WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_lbc(n_sequence_lbc))
+         IF ( crname_lbc(n_sequence_lbc) /= 'already counted' ) THEN
+            WRITE(numcom,'(A, I4, A, A)') ' - ', 1,' times by subroutine ', TRIM(crname_lbc(ncom_rec_max))
+         END IF
          WRITE(numcom,*) ' '
          IF ( n_sequence_glb > 0 ) THEN
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/ZDF/zdfphy.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/ZDF/zdfphy.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/ZDF/zdfphy.F90	(revision 11082)
@@ -132,6 +132,7 @@
       IF( ln_zad_Aimp ) THEN
          IF( zdf_phy_alloc() /= 0 )   &
-        &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' )
-         wi(:,:,:) = 0._wp
+            &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' )
+         Cu_adv(:,:,:) = 0._wp
+         wi    (:,:,:) = 0._wp
       ENDIF
       !                          !==  Background eddy viscosity and diffusivity  ==!
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/TOP/PISCES/P4Z/p4zpoc.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 11082)
@@ -166,5 +166,5 @@
                        &   + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) / tgfunc(ji,jj,jk-1) * ( 1.           &
                        &   - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) * zsizek ) + prodgoc(ji,jj,jk) &
-                       &   / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn) 
+                       &   / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn) * alphan(jn) 
                        alphat = alphat + alphag(ji,jj,jk,jn)
                        remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90	(revision 11081)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90	(revision 11082)
@@ -96,8 +96,8 @@
       !                                ! ==>>>  set by hand non-zero value on first/last columns & rows 
       DO ji = mi0(1), mi1(1)              ! first row of global domain only
-         zhu(ji,2) = zht(1,2)
-      END DO
-       DO ji = mi0(jpi), mi1(jpi)         ! last  row of global domain only
-         zhu(ji,2) = zht(jpi,2)
+         zhu(ji,2) = zht(ji,2)
+      END DO
+       DO ji = mi0(jpiglo), mi1(jpiglo)   ! last  row of global domain only
+         zhu(ji,2) = zht(ji,2)
       END DO
       zhu(:,1) = zhu(:,2)
