Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12812)
@@ -49,4 +49,6 @@
           <field id="icehpnd"      long_name="melt pond depth"                                         standard_name="sea_ice_meltpond_depth"                    unit="m" /> 
           <field id="icevpnd"      long_name="melt pond volume"                                        standard_name="sea_ice_meltpond_volume"                   unit="m" /> 
+          <field id="icehlid"      long_name="melt pond lid depth"                                     standard_name="sea_ice_meltpondlid_depth"                 unit="m" /> 
+          <field id="icevlid"      long_name="melt pond lid volume"                                    standard_name="sea_ice_meltpondlid_volume"                unit="m" /> 
 	  
 	  <!-- heat -->
@@ -171,4 +173,7 @@
           <field id="frq_m"    unit="-"    />
 
+          <!-- rheology convergence tests -->
+          <field id="uice_cvg"   long_name="sea ice velocity convergence"   standard_name="sea_ice_velocity_convergence"   unit="m/s" />
+
 	  <!-- ================= -->
           <!-- Add-ons for SIMIP -->
@@ -252,5 +257,7 @@
           <field id="iceapnd_cat"  long_name="Ice melt pond concentration per category"          unit=""        /> 
           <field id="icehpnd_cat"  long_name="Ice melt pond thickness per category"              unit="m"       detect_missing_value="true" /> 
+          <field id="icehlid_cat"  long_name="Ice melt pond lid thickness per category"          unit="m"       detect_missing_value="true" /> 
           <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"               unit=""        /> 
+          <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category"     unit=""        /> 
           <field id="icemask_cat"  long_name="Fraction of time step with sea ice (per category)" unit=""        />
           <field id="iceage_cat"   long_name="Ice age per category"                              unit="days"    detect_missing_value="true" />
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12812)
@@ -27,4 +27,13 @@
         <field id="soce"         long_name="salinity"                            standard_name="sea_water_practical_salinity"      unit="1e-3"     grid_ref="grid_T_3D"/>
         <field id="soce_e3t"     long_name="salinity    (thickness weighted)"                                                      unit="1e-3"     grid_ref="grid_T_3D" > soce * e3t </field >
+
+	<!--- additions to diawri.F90 --->
+        <field id="socegrad"    long_name="module of salinity gradient"              unit="psu/m"   grid_ref="grid_T_3D"/>
+        <field id="socegrad2"   long_name="square of module of salinity gradient"    unit="psu2/m2" grid_ref="grid_T_3D"/>
+        <field id="eken_int"    long_name="vertical integration of kinetic energy"   unit="m3/s2"   />
+        <field id="relvor"      long_name="relative vorticity"                       unit="s-1"    grid_ref="grid_T_3D"/>
+        <field id="absvor"      long_name="absolute vorticity"                       unit="s-1"    grid_ref="grid_T_3D"/>
+        <field id="potvor"      long_name="potential vorticity"                      unit="s-1"    grid_ref="grid_T_3D"/>
+        <field id="salt2c"      long_name="Salt content vertically integrated"       unit="1e-3*kg/m2" />
 
         <!-- t-eddy viscosity coefficients (ldfdyn) -->
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/namelist_ice_ref	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/namelist_ice_ref	(revision 12812)
@@ -97,4 +97,5 @@
       rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 
                                       !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300)
+   ln_rhg_chkcvg    = .false.         !  check convergence of rheology (outputs: file ice_cvg.nc & variable uice_cvg)
 /
 !------------------------------------------------------------------------------
@@ -126,4 +127,6 @@
    ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F)
    ln_icedS         = .true.          !  activate brine drainage (T) or not (F)
+   !
+   ln_leadhfx       = .true.          !  heat in the leads is used to melt sea-ice before warming the ocean
 /
 !------------------------------------------------------------------------------
@@ -176,5 +179,9 @@
 !------------------------------------------------------------------------------
    ln_pnd           = .false.         !  activate melt ponds or not
-     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
+     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
+       ln_pnd_lids  = .true.          !  ponds with frozen lids
+       ln_pnd_flush = .true.          !  ponds flushing trhu the ice  
+       rn_apnd_min  =   0.15          !  minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ??
+       rn_apnd_max  =   0.85          !  maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ??
      ln_pnd_CST     = .false.         !  activate constant  melt ponds
        rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC
@@ -186,5 +193,7 @@
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file   = .false.         !  netcdf file provided for initialization (T) or not (F)
+   nn_iceini_file   =   0             !     0 = Initialise sea ice based on SSTs
+                                      !     1 = Initialise sea ice from single category netcdf file
+                                      !     2 = Initialise sea ice from multi category restart file
    rn_thres_sst     =   2.0           !  max temp. above Tfreeze with initial ice = (sst - tfreeze)
    rn_hti_ini_n     =   3.0           !  initial ice thickness       (m), North
@@ -206,5 +215,7 @@
    rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North
    rn_hpd_ini_s     =   0.05          !        "            "             South
-   ! -- for ln_iceini_file = T
+   rn_hld_ini_n     =   0.0           !  initial pond lid depth      (m), North
+   rn_hld_ini_s     =   0.0           !        "            "             South
+   ! -- for nn_iceini_file = 1
    sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    sn_hts = 'Ice_initialization'    , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', ''
@@ -217,4 +228,5 @@
    sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_hld = 'NOT USED'              , -12 ,'hld'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    cn_dir='./'
 /
@@ -238,5 +250,5 @@
    ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F)
    ln_icectl        = .false.         !  ice points output for debug (T or F)
-   iiceprt          =  10             !  i-index for debug
-   jiceprt          =  10             !  j-index for debug
-/
+      iiceprt       =  10             !     i-index for debug
+      jiceprt       =  10             !     j-index for debug
+/
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/namelist_ref	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/cfgs/SHARED/namelist_ref	(revision 12812)
@@ -281,4 +281,5 @@
    sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
+   sn_cc       = 'NOT USED'                   ,   24         , 'CC'      ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_tdif     = 'taudif_core'                ,   24         , 'taudif'  ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
 /
@@ -286,8 +287,9 @@
 &namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_oasis3")
 !-----------------------------------------------------------------------
-   nn_cplmodel   =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data
-   ln_usecplmask = .false. !  use a coupling mask file to merge data received from several models
-   !                       !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
-   nn_cats_cpl   =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1)
+   nn_cplmodel       =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data
+   ln_usecplmask     = .false. !  use a coupling mask file to merge data received from several models
+   !                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   ln_scale_ice_flux = .false. !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration)
+   nn_cats_cpl       =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1)
 
    !_____________!__________________________!____________!_____________!______________________!________!
@@ -641,4 +643,5 @@
    bn_aip      = 'NOT USED'              ,         24.       , 'siapnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    bn_hip      = 'NOT USED'              ,         24.       , 'sihpnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
+   bn_hil      = 'NOT USED'              ,         24.       , 'sihlid'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    ! if bn_t_i etc are "not used", then define arbitrary temperatures and salinity and ponds
    rn_ice_tem  = 270.         !  arbitrary temperature               of incoming sea ice
@@ -647,4 +650,5 @@
    rn_ice_apnd = 0.2          !       --   pond fraction = a_ip/a_i            --
    rn_ice_hpnd = 0.05         !       --   pond depth                          --
+   rn_ice_hlid = 0.0          !       --   pond lid depth                      --
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/nambdy_dta
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/nambdy_dta	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/nambdy_dta	(revision 12812)
@@ -29,4 +29,5 @@
    bn_aip      = 'NOT USED'              ,         24.       , 'siapnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    bn_hip      = 'NOT USED'              ,         24.       , 'sihpnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
+   bn_hil      = 'NOT USED'              ,         24.       , 'sihlid'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    ! if bn_t_i etc are "not used", then define arbitrary temperatures and salinity and ponds
    rn_ice_tem  = 270.         !  arbitrary temperature               of incoming sea ice
@@ -35,3 +36,4 @@
    rn_ice_apnd = 0.2          !       --   pond fraction = a_ip/a_i            --
    rn_ice_hpnd = 0.05         !       --   pond depth                          --
+   rn_ice_hlid = 0.0          !       --   pond lid depth                      --
 /
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namdia
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namdia	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namdia	(revision 12812)
@@ -8,5 +8,5 @@
    ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F)
    ln_icectl        = .false.         !  ice points output for debug (T or F)
-   iiceprt          =  10             !  i-index for debug
-   jiceprt          =  10             !  j-index for debug
+      iiceprt       =  10             !     i-index for debug
+      jiceprt       =  10             !     j-index for debug
 /
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namdyn_rhg
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namdyn_rhg	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namdyn_rhg	(revision 12812)
@@ -9,3 +9,4 @@
       rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 
                                       !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300)
+   ln_rhg_chkcvg    = .false.         !  check convergence of rheology (outputs: file ice_cvg.nc & variable uice_cvg)
 /
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namini
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namini	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namini	(revision 12812)
@@ -3,5 +3,7 @@
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file   = .false.         !  netcdf file provided for initialization (T) or not (F)
+   nn_iceini_file   =   0             !     0 = Initialise sea ice based on SSTs
+                                      !     1 = Initialise sea ice from single category netcdf file
+                                      !     2 = Initialise sea ice from multi category restart file
    rn_thres_sst     =   2.0           !  max temp. above Tfreeze with initial ice = (sst - tfreeze)
    rn_hti_ini_n     =   3.0           !  initial ice thickness       (m), North
@@ -23,5 +25,7 @@
    rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North
    rn_hpd_ini_s     =   0.05          !        "            "             South
-   ! -- for ln_iceini_file = T
+   rn_hld_ini_n     =   0.0           !  initial pond lid depth      (m), North
+   rn_hld_ini_s     =   0.0           !        "            "             South
+   ! -- for nn_iceini_file = 1
    sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    sn_hts = 'Ice_initialization'    , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', ''
@@ -34,4 +38,5 @@
    sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_hld = 'NOT USED'              , -12 ,'hld'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    cn_dir='./'
 /
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namsbc_blk
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namsbc_blk	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namsbc_blk	(revision 12812)
@@ -31,4 +31,5 @@
    sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
+   sn_cc       = 'NOT USED'                   ,   24         , 'CC'      ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_tdif     = 'taudif_core'                ,   24         , 'taudif'  ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
 /
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namsbc_cpl
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namsbc_cpl	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namsbc_cpl	(revision 12812)
@@ -2,9 +2,9 @@
 &namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_oasis3")
 !-----------------------------------------------------------------------
-   nn_cplmodel   =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data
-   ln_usecplmask = .false. !  use a coupling mask file to merge data received from several models
-   !                       !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
-   nn_cats_cpl   =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1)
-
+   nn_cplmodel       =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data
+   ln_usecplmask     = .false. !  use a coupling mask file to merge data received from several models
+   !                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   ln_scale_ice_flux = .false. !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration)
+   nn_cats_cpl       =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1)
    !_____________!__________________________!____________!_____________!______________________!________!
    !             !        description       !  multiple  !    vector   !       vector         ! vector !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namthd
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namthd	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namthd	(revision 12812)
@@ -6,3 +6,5 @@
    ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F)
    ln_icedS         = .true.          !  activate brine drainage (T) or not (F)
+   !                                   
+   ln_leadhfx       = .true.          !  heat in the leads is used to melt sea-ice before warming the ocean
 /
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namthd_pnd
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namthd_pnd	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/doc/namelists/namthd_pnd	(revision 12812)
@@ -3,5 +3,9 @@
 !------------------------------------------------------------------------------
    ln_pnd           = .false.         !  activate melt ponds or not
-     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
+     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
+       ln_pnd_lids  = .true.          !  ponds with frozen lids
+       ln_pnd_flush = .true.          !  ponds flushing trhu the ice  
+       rn_apnd_min  =   0.15          !  minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ??
+       rn_apnd_max  =   0.85          !  maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ??
      ln_pnd_CST     = .false.         !  activate constant  melt ponds
        rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/ice.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/ice.F90	(revision 12812)
@@ -70,4 +70,5 @@
    !! a_ip        |      -      |    Ice pond concentration       |       |
    !! v_ip        |      -      |    Ice pond volume per unit area| m     |
+   !! v_il        |    v_il_1d  |    Ice pond lid volume per area | m     |
    !!                                                                     |
    !!-------------|-------------|---------------------------------|-------|
@@ -85,4 +86,5 @@
    !! t_su        ! t_su_1d     |    Sea ice surface temperature  ! K     |
    !! h_ip        | h_ip_1d     |    Ice pond thickness           | m     |
+   !! h_il        | h_il_1d     |    Ice pond lid thickness       | m     |
    !!                                                                     |
    !! notes: the ice model only sees a bulk (i.e., vertically averaged)   |
@@ -112,4 +114,6 @@
    !! hm_ip       |      -      |    Mean ice pond depth          | m     |
    !! vt_ip       |      -      |    Total ice pond vol. per unit area| m |
+   !! hm_il       |      -      |    Mean ice pond lid depth      | m     |
+   !! vt_il       |      -      |    Total ice pond lid vol. per area | m |
    !!=====================================================================
 
@@ -151,4 +155,5 @@
    INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling
    REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
+   LOGICAL , PUBLIC ::   ln_rhg_chkcvg    !: check ice rheology convergence 
    !
    !                                     !!** ice-advection namelist (namdyn_adv) **
@@ -190,5 +195,9 @@
    !                                     !!** ice-ponds namelist (namthd_pnd)
    LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F)
-   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012
+   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010)
+   LOGICAL,  PUBLIC ::   ln_pnd_lids      !: Allow ponds to have frozen lids
+   LOGICAL,  PUBLIC ::   ln_pnd_flush     !: Allow ponds to flush thru the ice
+   REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum ice fraction that contributes to melt ponds
+   REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum ice fraction that contributes to melt ponds
    LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth
    REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1)
@@ -271,5 +280,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt                             [W.m-2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping => must be 0   [W.m-2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qt_atm_oi   !: heat flux at the interface atm-[oce+ice]            [W.m-2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qt_oce_ai   !: heat flux at the interface oce-[atm+ice]            [W.m-2]
@@ -331,19 +339,24 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area      [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond fraction (a_ip/a_i)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_eff   !: melt pond effective fraction (not covered up by lid) (a_ip/a_i)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond depth                          [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_il       !: melt pond lid volume                     [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_il       !: melt pond lid thickness                  [m]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond concentration
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hm_ip      !: mean melt pond depth                     [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per gridcell area [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hm_il      !: mean melt pond lid depth                     [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_il      !: total melt pond lid volume per gridcell area [m]
 
    !!----------------------------------------------------------------------
    !! * Old values of global variables
    !!----------------------------------------------------------------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b, h_ip_b    !: snow and ice volumes/thickness
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b                 !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                                 !: snow heat content
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                                 !: ice temperatures
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b                      !: ice velocity
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                                !: ice concentration (total)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b     !: snow and ice volumes/thickness
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b                  !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                          !: snow heat content
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                          !: ice temperatures
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b               !: ice velocity
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                         !: ice concentration (total)
             
    !!----------------------------------------------------------------------
@@ -386,4 +399,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice_top   !: Surface conduction flux (W/m2)
 
+   !!----------------------------------------------------------------------
+   !! * Only for atmospheric coupling
+   !!----------------------------------------------------------------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple !: Ice fractional area at last coupling time
    !
    !!----------------------------------------------------------------------
@@ -400,5 +417,5 @@
       INTEGER :: ice_alloc
       !
-      INTEGER :: ierr(16), ii
+      INTEGER :: ierr(17), ii
       !!-----------------------------------------------------------------
       ierr(:) = 0
@@ -424,5 +441,5 @@
          &      hfx_sum    (jpi,jpj) , hfx_bom   (jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) ,     &
          &      hfx_opw    (jpi,jpj) , hfx_thd   (jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) ,     &
-         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj)             ,  STAT=ierr(ii) )
+         &      hfx_err_dif(jpi,jpj) , wfx_err_sub(jpi,jpj)                   , STAT=ierr(ii) )
 
       ! * Ice global state variables
@@ -448,14 +465,15 @@
 
       ii = ii + 1
-      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) )
-
-      ii = ii + 1
-      ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) )
+      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl),  &
+         &      v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , STAT = ierr(ii) )
+
+      ii = ii + 1
+      ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , hm_il(jpi,jpj) , vt_il(jpi,jpj) , STAT = ierr(ii) )
 
       ! * Old values of global variables
       ii = ii + 1
-      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl), h_ip_b(jpi,jpj,jpl),  &
-         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,               &
-         &      oa_i_b(jpi,jpj,jpl)                                                   , STAT=ierr(ii) )
+      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl),         &
+         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , &
+         &      STAT=ierr(ii) )
 
       ii = ii + 1
@@ -481,7 +499,12 @@
       ALLOCATE( t_si(jpi,jpj,jpl) , tm_si(jpi,jpj) , qcn_ice_bot(jpi,jpj,jpl) , qcn_ice_top(jpi,jpj,jpl) , STAT = ierr(ii) )
 
+      ! * For atmospheric coupling
+      ii = ii + 1
+      ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(ii) )
+
       ice_alloc = MAXVAL( ierr(:) )
       IF( ice_alloc /= 0 )   CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' )
       !
+
    END FUNCTION ice_alloc
 
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/ice1d.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/ice1d.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/ice1d.F90	(revision 12812)
@@ -51,5 +51,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_dyn_1d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_dif_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qt_oce_ai_1d
@@ -124,8 +123,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oa_i_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   o_i_1d        !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !: ice ponds
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_ip_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ip_1d       !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_il_1d       !: Ice pond lid
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_il_1d       !:
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s
@@ -157,4 +157,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_il_2d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_2d
@@ -189,5 +190,5 @@
          &      hfx_thd_1d(jpij) , hfx_spr_1d    (jpij) ,                      &
          &      hfx_snw_1d(jpij) , hfx_sub_1d    (jpij) ,                      &
-         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , qt_oce_ai_1d(jpij), STAT=ierr(ii) )
+         &      hfx_res_1d(jpij) , hfx_err_dif_1d(jpij) , qt_oce_ai_1d(jpij), STAT=ierr(ii) )
       !
       ii = ii + 1
@@ -208,6 +209,6 @@
          &      dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm  (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) ,  &    
          &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  &
-         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) ,                   &
-         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) ,                                                   &
+         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , v_il_1d (jpij) ,  &
+         &      h_il_1d (jpij) , h_ip_1d (jpij) ,                                                       &
          &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) )
       !
@@ -226,5 +227,5 @@
       ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,  &
          &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,  &
-         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) ,                      &
+         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , v_il_2d(jpij,jpl) ,  &
          &      STAT=ierr(ii) )
 
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icealb.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icealb.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icealb.F90	(revision 12812)
@@ -45,5 +45,5 @@
 CONTAINS
 
-   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, palb_cs, palb_os )
+   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, pcloud_fra, palb_ice )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE ice_alb  ***
@@ -97,6 +97,6 @@
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area)
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_os      !  albedo of ice under overcast sky
+      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   pcloud_fra   !  cloud fraction
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_ice     !  albedo of ice
       !
       INTEGER  ::   ji, jj, jl                ! dummy loop indices
@@ -106,4 +106,5 @@
       REAL(wp) ::   zalb_ice, zafrac_ice      ! bare sea ice albedo & relative ice fraction
       REAL(wp) ::   zalb_snw, zafrac_snw      ! snow-covered sea ice albedo & relative snow fraction
+      REAL(wp) ::   zalb_cs, zalb_os          ! albedo of ice under clear/overcast sky
       !!---------------------------------------------------------------------
       !
@@ -119,6 +120,8 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               !                       !--- Specific snow, ice and pond fractions (for now, we prevent melt ponds and snow at the same time)
-               IF( ph_snw(ji,jj,jl) == 0._wp ) THEN
+               !---------------------------------------------!
+               !--- Specific snow, ice and pond fractions ---!
+               !---------------------------------------------!               
+               IF( ph_snw(ji,jj,jl) == 0._wp ) THEN   !--- no snow : we prevent melt ponds and snow at the same time (for now)
                   zafrac_snw = 0._wp
                   IF( ld_pnd_alb ) THEN
@@ -129,9 +132,12 @@
                   zafrac_ice = 1._wp - zafrac_pnd
                ELSE
-                  zafrac_snw = 1._wp      ! Snow fully "shades" melt ponds and ice
+                  zafrac_snw = 1._wp                  !--- snow : fully "shades" melt ponds and ice
                   zafrac_pnd = 0._wp
                   zafrac_ice = 0._wp
                ENDIF
                !
+               !---------------!
+               !--- Albedos ---!
+               !---------------!               
                !                       !--- Bare ice albedo (for hi > 150cm)
                IF( ld_pnd_alb ) THEN
@@ -161,10 +167,12 @@
                ENDIF
                !                       !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions
-               palb_os(ji,jj,jl) = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1)
-               !
-               palb_cs(ji,jj,jl) = palb_os(ji,jj,jl)  &
-                  &                - ( - 0.1010 * palb_os(ji,jj,jl) * palb_os(ji,jj,jl)  &
-                  &                    + 0.1933 * palb_os(ji,jj,jl) - 0.0148 ) * tmask(ji,jj,1)
-               !
+               zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1)
+               !
+               zalb_cs = zalb_os - ( - 0.1010 * zalb_os * zalb_os  &
+                  &                  + 0.1933 * zalb_os - 0.0148 ) * tmask(ji,jj,1)
+               !
+               ! albedo depends on cloud fraction because of non-linear spectral effects
+               palb_ice(ji,jj,jl) = ( 1._wp - pcloud_fra(ji,jj) ) * zalb_cs + pcloud_fra(ji,jj) * zalb_os
+
             END DO
          END DO
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icectl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icectl.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icectl.F90	(revision 12812)
@@ -625,5 +625,4 @@
                   WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl) 
                   WRITE(numout,*) ' sv_i       : ', sv_i(ji,jj,jl)             , ' sv_i_b     : ', sv_i_b(ji,jj,jl)   
-                  WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' oa_i_b     : ', oa_i_b(ji,jj,jl)
                END DO !jl
                
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn.F90	(revision 12812)
@@ -99,6 +99,8 @@
       WHERE( a_ip(:,:,:) >= epsi20 )
          h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
+         h_il(:,:,:) = v_il(:,:,:) / a_ip(:,:,:)
       ELSEWHERE
          h_ip(:,:,:) = 0._wp
+         h_il(:,:,:) = 0._wp
       END WHERE
       !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_adv.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_adv.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_adv.F90	(revision 12812)
@@ -84,10 +84,10 @@
          !                             !-----------------------!
          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 )
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, v_il, 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 )
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, v_il, e_s, e_i )
       END SELECT
 
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_adv_pra.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_adv_pra.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_adv_pra.F90	(revision 12812)
@@ -44,4 +44,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvl , syvl , sxxvl , syyvl , sxyvl    ! melt pond lid volume
 
    !! * Substitutions
@@ -55,5 +56,5 @@
 
    SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  &
-      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                **  routine ice_dyn_adv_pra  **
@@ -78,4 +79,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il     ! melt pond lid thickness
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -88,5 +90,5 @@
       REAL(wp), DIMENSION(jpi,jpj,1)          ::   z0opw
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
-      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0vl
       REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es
       REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei
@@ -129,4 +131,7 @@
             z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
             z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
+            IF ( ln_pnd_lids ) THEN
+               z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:)   ! Melt pond lid volume
+            ENDIF
          ENDIF
       END DO
@@ -167,4 +172,8 @@
                CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
                CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
+               IF ( ln_pnd_lids ) THEN
+                  CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
+                  CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 
+               ENDIF
             ENDIF
          END DO
@@ -202,4 +211,8 @@
                CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
                CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
+               IF ( ln_pnd_lids ) THEN
+                  CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
+                  CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 
+               ENDIF
             ENDIF
          END DO
@@ -225,4 +238,7 @@
             pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
             pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            IF ( ln_pnd_lids ) THEN
+               pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            ENDIF
          ENDIF
       END DO
@@ -651,6 +667,7 @@
          &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
          &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
-         &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
-         &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
+         &      sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
+         &      sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
+         &      sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) ,   &
          !
          &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
@@ -765,4 +782,12 @@
                CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
                CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
+               !
+               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume
+                  CALL iom_get( numrir, jpdom_autoglo, 'sxvl' , sxvl  )
+                  CALL iom_get( numrir, jpdom_autoglo, 'syvl' , syvl  )
+                  CALL iom_get( numrir, jpdom_autoglo, 'sxxvl', sxxvl )
+                  CALL iom_get( numrir, jpdom_autoglo, 'syyvl', syyvl )
+                  CALL iom_get( numrir, jpdom_autoglo, 'sxyvl', sxyvl )
+               ENDIF
             ENDIF
             !
@@ -780,6 +805,9 @@
             sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
             IF( ln_pnd_H12 ) THEN
-               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
-               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
+               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction
+               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume
+               IF ( ln_pnd_lids ) THEN
+                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume
+               ENDIF
             ENDIF
          ENDIF
@@ -862,4 +890,12 @@
             CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
             CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
+            !
+            IF ( ln_pnd_lids ) THEN                                  ! melt pond lid volume
+               CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  )
+               CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  )
+               CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl )
+               CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl )
+               CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl )
+            ENDIF
          ENDIF
          !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_adv_umx.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_adv_umx.F90	(revision 12812)
@@ -60,5 +60,5 @@
 
    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 )
+      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ice_dyn_adv_umx  ***
@@ -85,4 +85,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond concentration
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -334,4 +335,11 @@
             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 )
+            ! lid
+            IF ( ln_pnd_lids ) THEN
+               zamsk = 0._wp
+               zhvar(:,:,:) = pv_il(:,:,:) * z1_aip(:,:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+                  &                                      zhvar, pv_il, zua_ups, zva_ups )
+            ENDIF
          ENDIF
          !
@@ -350,5 +358,5 @@
          ! Remove negative values (conservation is ensured)
          !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
-         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
          !
          ! Make sure ice thickness is not too big
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_rdgrft.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_rdgrft.F90	(revision 12812)
@@ -503,6 +503,6 @@
       REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
       REAL(wp)                  ::   airft1, oirft1, aprft1
-      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges
-      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice
+      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg  ! area etc of new ridges
+      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft  ! area etc of rafted ice
       !
       REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
@@ -581,4 +581,8 @@
                   aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
                   vprft (ji) = v_ip_2d(ji,jl1) * afrft
+                  IF ( ln_pnd_lids ) THEN
+                     vlrdg (ji) = v_il_2d(ji,jl1) * afrdg
+                     vlrft (ji) = v_il_2d(ji,jl1) * afrft
+                  ENDIF
                ENDIF
 
@@ -610,4 +614,7 @@
                   a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
+                  IF ( ln_pnd_lids ) THEN
+                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - vlrdg(ji) - vlrft(ji)
+                  ENDIF
                ENDIF
             ENDIF
@@ -706,4 +713,8 @@
                      a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
                         &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
+                     IF ( ln_pnd_lids ) THEN
+                        v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + (   vlrdg (ji) * rn_fpndrdg * fvol   (ji)   &
+                           &                                   + vlrft (ji) * rn_fpndrft * zswitch(ji)   )
+                     ENDIF
                   ENDIF
                   
@@ -736,5 +747,5 @@
       !----------------
       ! In case ridging/rafting lead to very small negative values (sometimes it happens)
-      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 )
+      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, v_il_2d, ze_s_2d, ze_i_2d )
       !
    END SUBROUTINE rdgrft_shift
@@ -848,4 +859,5 @@
          CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
          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), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
@@ -874,4 +886,5 @@
          CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
          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), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_rhg.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_rhg.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_rhg.F90	(revision 12812)
@@ -110,5 +110,5 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg
       !!-------------------------------------------------------------------
       !
@@ -126,10 +126,11 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist : namdyn_rhg:'
-         WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP = ', ln_rhg_EVP
-         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP    = ', ln_aEVP
-         WRITE(numout,*) '         creep limit                                       rn_creepl  = ', rn_creepl
-         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc     = ', rn_ecc
-         WRITE(numout,*) '         number of iterations for subcycling               nn_nevp    = ', nn_nevp
-         WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast  = ', rn_relast
+         WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP    = ', ln_rhg_EVP
+         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP       = ', ln_aEVP
+         WRITE(numout,*) '         creep limit                                       rn_creepl     = ', rn_creepl
+         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc        = ', rn_ecc
+         WRITE(numout,*) '         number of iterations for subcycling               nn_nevp       = ', nn_nevp
+         WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast     = ', rn_relast
+         WRITE(numout,*) '      check convergence of rheology                        ln_rhg_chkcvg = ', ln_rhg_chkcvg
       ENDIF
       !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_rhg_evp.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icedyn_rhg_evp.F90	(revision 12812)
@@ -41,4 +41,5 @@
    USE prtctl         ! Print control
 
+   USE netcdf         ! NetCDF library for convergence test
    IMPLICIT NONE
    PRIVATE
@@ -49,4 +50,9 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+
+   !! for convergence tests
+   INTEGER ::   ncvgid   ! netcdf file id
+   INTEGER ::   nvarid   ! netcdf variable id
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -125,5 +131,4 @@
       REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
       !
-      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
       REAL(wp) ::   zintb, zintn                                        ! dummy argument
       REAL(wp) ::   zfac_x, zfac_y
@@ -141,5 +146,4 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
-!!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
       REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
       !                                                                 !    ocean surface (ssh_m) if ice is not embedded
@@ -160,6 +164,7 @@
       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
+      !! --- check convergence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
       !! --- diags
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
       !! --- SIMIP diags
@@ -174,5 +179,14 @@
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
       !
-!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
+      ! for diagnostics and convergence tests
+      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) )
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
+            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
+         END DO
+      END DO
+      !
+      !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
       !------------------------------------------------------------------------------!
       ! 0) mask at F points for the ice
@@ -345,10 +359,13 @@
          l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
          !
-!!$         IF(ln_ctl) THEN   ! Convergence test
-!!$            DO jj = 1, jpjm1
-!!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
-!!$               zv_ice(:,jj) = v_ice(:,jj)
-!!$            END DO
-!!$         ENDIF
+         ! convergence test
+         IF(ln_rhg_chkcvg) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
+                  zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
+               END DO
+            END DO
+         ENDIF
 
          ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
@@ -667,11 +684,8 @@
          ENDIF
 
-!!$         IF(ln_ctl) THEN   ! Convergence test
-!!$            DO jj = 2 , jpjm1
-!!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
-!!$            END DO
-!!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
-!!$            CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
-!!$         ENDIF
+         ! convergence test
+         IF(ln_rhg_chkcvg) THEN
+            CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
+         ENDIF
          !
          !                                                ! ==================== !
@@ -734,10 +748,4 @@
       ! 5) diagnostics
       !------------------------------------------------------------------------------!
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
-         END DO
-      END DO
-
       ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
       IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
@@ -796,5 +804,5 @@
          DEALLOCATE( zsig1 , zsig2 , zsig3 )
       ENDIF
-      
+
       ! --- SIMIP --- !
       IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
@@ -852,5 +860,85 @@
       ENDIF
       !
+      ! --- convergence tests --- !
+      IF( ln_rhg_chkcvg ) THEN
+         IF( iom_use('uice_cvg') ) THEN  ! output: u(t=nn_nevp) - u(t=nn_nevp-1)
+            CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
+               &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
+         ENDIF
+      ENDIF      
+      !
+      DEALLOCATE( zmsk00, zmsk15 )
+      !
    END SUBROUTINE ice_dyn_rhg_evp
+
+
+   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE rhg_cvg  ***
+      !!                     
+      !! ** Purpose :   check convergence of oce rheology
+      !!
+      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity
+      !!                during the sub timestepping of rheology so as:
+      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) )
+      !!                This routine is called every sub-iteration, so it is cpu expensive
+      !!
+      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)   
+      !!----------------------------------------------------------------------
+      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
+      !!
+      INTEGER           ::   it, idtime, istatus
+      INTEGER           ::   ji, jj          ! dummy loop indices
+      REAL(wp)          ::   zresm           ! local real 
+      CHARACTER(len=20) ::   clname
+      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence
+      !!----------------------------------------------------------------------
+
+      ! create file
+      IF( kt == nit000 .AND. kiter == 1 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'
+            WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !
+         IF( lwm ) THEN
+            clname = 'ice_cvg.nc'
+            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
+            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
+            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
+            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid )
+            istatus = NF90_ENDDEF(ncvgid)
+         ENDIF
+         !
+      ENDIF
+
+      ! time
+      it = ( kt - 1 ) * kitermax + kiter
+      
+      ! convergence
+      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
+         zresm = 0._wp
+      ELSE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
+                  &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj)
+            END DO
+         END DO
+         zresm = MAXVAL( zres )
+         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
+      ENDIF
+
+      IF( lwm ) THEN
+         ! write variables
+         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
+         ! close file
+         IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid)
+      ENDIF
+      
+   END SUBROUTINE rhg_cvg
 
 
@@ -910,4 +998,5 @@
    END SUBROUTINE rhg_evp_rst
 
+   
 #else
    !!----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/iceistate.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/iceistate.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/iceistate.F90	(revision 12812)
@@ -41,13 +41,16 @@
    !                             !! ** namelist (namini) **
    LOGICAL, PUBLIC  ::   ln_iceini        !: Ice initialization or not
-   LOGICAL, PUBLIC  ::   ln_iceini_file   !: Ice initialization from 2D netcdf file
+   INTEGER, PUBLIC  ::   nn_iceini_file   !: Ice initialization:
+                                  !        0 = Initialise sea ice based on SSTs
+                                  !        1 = Initialise sea ice from single category netcdf file
+                                  !        2 = Initialise sea ice from multi category restart file
    REAL(wp) ::   rn_thres_sst
    REAL(wp) ::   rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n
    REAL(wp) ::   rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s
-   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n
-   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s
+   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n, rn_hld_ini_n
+   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s, rn_hld_ini_s
    !
-   !                              ! if ln_iceini_file = T
-   INTEGER , PARAMETER ::   jpfldi = 9           ! maximum number of files to read
+   !                              ! if nn_iceini_file = 1
+   INTEGER , PARAMETER ::   jpfldi = 10          ! maximum number of files to read
    INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m)
    INTEGER , PARAMETER ::   jp_hts = 2           ! index of snw thickness    (m)
@@ -59,4 +62,5 @@
    INTEGER , PARAMETER ::   jp_apd = 8           ! index of pnd fraction     (-)
    INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m)
+   INTEGER , PARAMETER ::   jp_hld = 10          ! index of pnd lid depth    (m)
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
    !   
@@ -98,8 +102,8 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file
       REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini, zhlid_ini            !data from namelist or nc file
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays
       !!
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d
       !--------------------------------------------------------------------
 
@@ -155,6 +159,8 @@
       a_ip     (:,:,:) = 0._wp
       v_ip     (:,:,:) = 0._wp
-      a_ip_frac(:,:,:) = 0._wp
+      v_il     (:,:,:) = 0._wp
+      a_ip_eff (:,:,:) = 0._wp
       h_ip     (:,:,:) = 0._wp
+      h_il     (:,:,:) = 0._wp
       !
       ! ice velocities
@@ -167,5 +173,5 @@
       IF( ln_iceini ) THEN
          !                             !---------------!
-         IF( ln_iceini_file )THEN      ! Read a file   !
+         IF( nn_iceini_file == 1 )THEN ! Read a file   !
             !                          !---------------!
             WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp
@@ -218,5 +224,11 @@
             IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
                &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
+            !
+            ! pond lid depth
+            IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) &
+               &     si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
+            !
             zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
+            zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1)
             !
             ! change the switch for the following
@@ -243,4 +255,5 @@
                zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
                zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
+               zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:)
             ELSEWHERE
                zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
@@ -253,4 +266,5 @@
                zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
                zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
+               zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:)
             END WHERE
             !
@@ -261,4 +275,9 @@
             zapnd_ini(:,:) = 0._wp
             zhpnd_ini(:,:) = 0._wp
+            zhlid_ini(:,:) = 0._wp
+         ENDIF
+
+         IF ( .NOT.ln_pnd_lids ) THEN
+            zhlid_ini(:,:) = 0._wp
          ENDIF
          
@@ -287,14 +306,18 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini )
          CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti)  , zhlid_ini )
 
          ! allocate temporary arrays
-         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
-            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
+         ALLOCATE( zhi_2d (npti,jpl), zhs_2d (npti,jpl), zai_2d (npti,jpl), &
+            &      zti_2d (npti,jpl), zts_2d (npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), &
+            &      zaip_2d(npti,jpl), zhip_2d(npti,jpl), zhil_2d(npti,jpl) )
          
          ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
-         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   &
-            &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   &
-            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), &
-            &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d )
+         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                  &
+            &              zhi_2d          , zhs_2d          , zai_2d         ,                  &
+            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti),                  &
+            &              s_i_1d(1:npti)  , a_ip_1d(1:npti) , h_ip_1d(1:npti), h_il_1d(1:npti), &
+            &              zti_2d          , zts_2d          , ztsu_2d        ,                  &
+            &              zsi_2d          , zaip_2d         , zhip_2d        , zhil_2d )
 
          ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
@@ -312,8 +335,9 @@
          CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
          CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d  , h_il   )
 
          ! deallocate temporary arrays
          DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
-            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
+            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d )
 
          ! calculate extensive and intensive variables
@@ -357,10 +381,9 @@
 
          ! Melt ponds
-         WHERE( a_i > epsi10 )
-            a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
-         ELSEWHERE
-            a_ip_frac(:,:,:) = 0._wp
+         WHERE( a_i > epsi10 )   ;   a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
+         ELSEWHERE               ;   a_ip_eff(:,:,:) = 0._wp
          END WHERE
          v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
+         v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
           
          ! specific temperatures for coupled runs
@@ -434,5 +457,4 @@
       e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
       sv_i_b (:,:,:)   = sv_i (:,:,:)
-      oa_i_b (:,:,:)   = oa_i (:,:,:)
       u_ice_b(:,:)     = u_ice(:,:)
       v_ice_b(:,:)     = v_ice(:,:)
@@ -463,13 +485,13 @@
       !
       CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
-      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
+      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd, sn_hld
       TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
       !
-      NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
+      NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, &
          &             rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
          &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
          &             rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
-         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
-         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
+         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, rn_hld_ini_n, rn_hld_ini_s, &
+         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, sn_hld, cn_dir
       !!-----------------------------------------------------------------------------
       !
@@ -485,5 +507,5 @@
       slf_i(jp_ati) = sn_ati  ;  slf_i(jp_smi) = sn_smi
       slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_tsu) = sn_tsu   ;   slf_i(jp_tms) = sn_tms
-      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd
+      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd   ;   slf_i(jp_hld) = sn_hld
       !
       IF(lwp) THEN                          ! control print
@@ -493,7 +515,7 @@
          WRITE(numout,*) '   Namelist namini:'
          WRITE(numout,*) '      ice initialization (T) or not (F)                ln_iceini      = ', ln_iceini
-         WRITE(numout,*) '      ice initialization from a netcdf file            ln_iceini_file = ', ln_iceini_file
+         WRITE(numout,*) '      ice initialization from a netcdf file            nn_iceini_file = ', nn_iceini_file
          WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst
-         IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
+         IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN
             WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s 
             WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s
@@ -505,8 +527,9 @@
             WRITE(numout,*) '      initial pnd fraction  in the north-south         rn_apd_ini     = ', rn_apd_ini_n,rn_apd_ini_s
             WRITE(numout,*) '      initial pnd depth     in the north-south         rn_hpd_ini     = ', rn_hpd_ini_n,rn_hpd_ini_s
+            WRITE(numout,*) '      initial pnd lid depth in the north-south         rn_hld_ini     = ', rn_hld_ini_n,rn_hld_ini_s
          ENDIF
       ENDIF
       !
-      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
+      IF( nn_iceini_file == 1 ) THEN                      ! Ice initialization using input file
          !
          ! set si structure
@@ -529,5 +552,10 @@
          rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
          rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
-         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' )
+         rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
+         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 & rn_hld_ini = 0 when no ponds' )
+      ENDIF
+      !
+      IF( .NOT.ln_pnd_lids ) THEN
+         rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
       ENDIF
       !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/iceitd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/iceitd.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/iceitd.F90	(revision 12812)
@@ -411,4 +411,5 @@
       CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
       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), v_il_2d(1:npti,1:jpl), v_il )
       CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
       DO jl = 1, jpl
@@ -483,4 +484,10 @@
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
                   v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
+                  !
+                  IF ( ln_pnd_lids ) THEN                            ! Pond lid volume
+                     ztrans          = v_il_2d(ji,jl1) * zworka(ji)
+                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans
+                     v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans
+                  ENDIF
                ENDIF
                !
@@ -527,5 +534,5 @@
       ! 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 )
+      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, v_il_2d, ze_s_2d, ze_i_2d )
 
       ! at_i must be <= rn_amax
@@ -555,4 +562,5 @@
       CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
       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), v_il_2d(1:npti,1:jpl), v_il )
       CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
       DO jl = 1, jpl
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icerst.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icerst.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icerst.F90	(revision 12812)
@@ -132,4 +132,5 @@
       CALL iom_rstput( iter, nitrst, numriw, 'a_ip' , a_ip  )
       CALL iom_rstput( iter, nitrst, numriw, 'v_ip' , v_ip  )
+      CALL iom_rstput( iter, nitrst, numriw, 'v_il' , v_il  )
       ! Snow enthalpy
       DO jk = 1, nlay_s 
@@ -171,5 +172,5 @@
       INTEGER           ::   jk
       LOGICAL           ::   llok
-      INTEGER           ::   id0, id1, id2, id3, id4   ! local integer
+      INTEGER           ::   id0, id1, id2, id3, id4, id5   ! local integer
       CHARACTER(len=25) ::   znam
       CHARACTER(len=2)  ::   zchar, zchar1
@@ -250,9 +251,17 @@
             v_ip(:,:,:) = 0._wp
          ENDIF
+         ! melt pond lids
+         id3 = iom_varid( numrir, 'v_il' , ldstop = .FALSE. )
+         IF( id3 > 0 ) THEN
+            CALL iom_get( numrir, jpdom_autoglo, 'v_il', v_il)
+         ELSE
+            IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds lids output then set it to zero'
+            v_il(:,:,:) = 0._wp
+         ENDIF
          ! fields needed for Met Office (Jules) coupling
          IF( ln_cpl ) THEN
-            id3 = iom_varid( numrir, 'cnd_ice' , ldstop = .FALSE. )
-            id4 = iom_varid( numrir, 't1_ice'  , ldstop = .FALSE. )
-            IF( id3 > 0 .AND. id4 > 0 ) THEN         ! fields exist
+            id4 = iom_varid( numrir, 'cnd_ice' , ldstop = .FALSE. )
+            id5 = iom_varid( numrir, 't1_ice'  , ldstop = .FALSE. )
+            IF( id4 > 0 .AND. id5 > 0 ) THEN         ! fields exist
                CALL iom_get( numrir, jpdom_autoglo, 'cnd_ice', cnd_ice )
                CALL iom_get( numrir, jpdom_autoglo, 't1_ice' , t1_ice  )
@@ -274,6 +283,6 @@
          CALL ice_istate( nit000 )
          !
-         IF( .NOT.ln_iceini .OR. .NOT.ln_iceini_file ) &
-            &   CALL ctl_stop('STOP', 'ice_rst_read: you need ln_ice_ini=T and ln_iceini_file=T')
+         IF( .NOT.ln_iceini .OR. nn_iceini_file == 0 ) &
+            &   CALL ctl_stop('STOP', 'ice_rst_read: you need ln_ice_ini=T and nn_iceini_file=0')
          !
       ENDIF
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icesbc.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icesbc.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icesbc.F90	(revision 12812)
@@ -116,6 +116,5 @@
       INTEGER  ::   ji, jj, jl      ! dummy loop index
       REAL(wp) ::   zmiss_val       ! missing value retrieved from xios 
-      REAL(wp), DIMENSION(jpi,jpj,jpl)              ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
-      REAL(wp), DIMENSION(:,:)        , ALLOCATABLE ::   zalb, zmsk00      ! 2D workspace
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zalb, zmsk00      ! 2D workspace
       !!--------------------------------------------------------------------
       !
@@ -131,11 +130,7 @@
       CALL iom_miss_val( "icetemp", zmiss_val )
 
-      ! --- cloud-sky and overcast-sky ice albedos --- !
-      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
-
-      ! albedo depends on cloud fraction because of non-linear spectral effects
-!!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument !
-      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-      !
+      ! --- ice albedo --- !
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, cloud_fra, alb_ice )
+
       !
       SELECT CASE( ksbc )   !== fluxes over sea ice ==!
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icestp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icestp.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icestp.F90	(revision 12812)
@@ -252,9 +252,10 @@
       !
       !                                ! Initial sea-ice state
-      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
+
+      IF ( ln_rstart .OR. nn_iceini_file == 2 ) THEN
+         CALL ice_rst_read                      ! start from a restart file
+      ELSE
          CALL ice_istate_init
-         CALL ice_istate( nit000 )
-      ELSE                                    ! start from a restart file
-         CALL ice_rst_read
+         CALL ice_istate( nit000 )              ! start from rest: sea-ice deduced from sst
       ENDIF
       CALL ice_var_glo2eqv
@@ -363,5 +364,4 @@
       v_s_b (:,:,:)   = v_s (:,:,:)     ! snow volume
       sv_i_b(:,:,:)   = sv_i(:,:,:)     ! salt content
-      oa_i_b(:,:,:)   = oa_i(:,:,:)     ! areal age content
       e_s_b (:,:,:,:) = e_s (:,:,:,:)   ! snow thermal energy
       e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy
@@ -372,10 +372,4 @@
          h_i_b(:,:,:) = 0._wp
          h_s_b(:,:,:) = 0._wp
-      END WHERE
-      
-      WHERE( a_ip(:,:,:) >= epsi20 )
-         h_ip_b(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)   ! ice pond thickness
-      ELSEWHERE
-         h_ip_b(:,:,:) = 0._wp
       END WHERE
       !
@@ -421,5 +415,4 @@
       hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
       hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
-      hfx_err_rem(:,:) = 0._wp
       hfx_err_dif(:,:) = 0._wp
       wfx_err_sub(:,:) = 0._wp
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icethd.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icethd.F90	(revision 12812)
@@ -51,4 +51,5 @@
    LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F)
    LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F)
+   LOGICAL ::   ln_leadhfx       !  heat in the leads is used to melt sea-ice before warming the ocean
 
    !! * Substitutions
@@ -164,5 +165,7 @@
             ! 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
+               IF( ln_leadhfx ) 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
+               ELSE                    ;   fhld(ji,jj) = 0._wp
+               ENDIF
                qlead(ji,jj) = 0._wp
             ELSE
@@ -354,5 +357,5 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
          CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            )
@@ -406,5 +409,4 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res       )
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   )
-         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem   )
          CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai     )
          !
@@ -441,4 +443,5 @@
          sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
          v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
+         v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)
          oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti)
          
@@ -460,5 +463,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
-         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
@@ -498,5 +501,4 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     )
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )
          CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai   )
          !
@@ -515,4 +517,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
          !
@@ -536,5 +539,5 @@
       INTEGER  ::   ios   ! Local integer output status for namelist read
       !!
-      NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS
+      NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS, ln_leadhfx
       !!-------------------------------------------------------------------
       !
@@ -552,8 +555,9 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namthd:'
-         WRITE(numout,*) '      activate ice thick change from top/bot (T) or not (F)   ln_icedH  = ', ln_icedH
-         WRITE(numout,*) '      activate lateral melting (T) or not (F)                 ln_icedA  = ', ln_icedA
-         WRITE(numout,*) '      activate ice growth in open-water (T) or not (F)        ln_icedO  = ', ln_icedO
-         WRITE(numout,*) '      activate gravity drainage and flushing (T) or not (F)   ln_icedS  = ', ln_icedS
+         WRITE(numout,*) '      activate ice thick change from top/bot (T) or not (F)                ln_icedH   = ', ln_icedH
+         WRITE(numout,*) '      activate lateral melting (T) or not (F)                              ln_icedA   = ', ln_icedA
+         WRITE(numout,*) '      activate ice growth in open-water (T) or not (F)                     ln_icedO   = ', ln_icedO
+         WRITE(numout,*) '      activate gravity drainage and flushing (T) or not (F)                ln_icedS   = ', ln_icedS
+         WRITE(numout,*) '      heat in the leads is used to melt sea-ice before warming the ocean   ln_leadhfx = ', ln_leadhfx
      ENDIF
       !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icethd_ent.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icethd_ent.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icethd_ent.F90	(revision 12812)
@@ -128,8 +128,8 @@
       ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_do), 
       ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
-      DO ji = 1, npti
-         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice *  &
-            &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) 
-      END DO
+      !DO ji = 1, npti
+      !   hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice *  &
+      !      &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) 
+      !END DO
       
    END SUBROUTINE ice_thd_ent
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icethd_pnd.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icethd_pnd.F90	(revision 12812)
@@ -88,11 +88,11 @@
          !
          IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
-            a_ip_frac_1d(ji) = rn_apnd
             h_ip_1d(ji)      = rn_hpnd    
-            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
+            a_ip_1d(ji)      = rn_apnd * a_i_1d(ji)
+            h_il_1d(ji)      = 0._wp    ! no pond lids whatsoever
          ELSE
-            a_ip_frac_1d(ji) = 0._wp
             h_ip_1d(ji)      = 0._wp    
             a_ip_1d(ji)      = 0._wp
+            h_il_1d(ji)      = 0._wp
          ENDIF
          !
@@ -106,84 +106,202 @@
       !!                ***  ROUTINE pnd_H12  ***
       !!
-      !! ** Purpose    : Compute melt pond evolution
-      !!
-      !! ** Method     : Empirical method. A fraction of meltwater is accumulated in ponds 
-      !!                 and sent to ocean when surface is freezing
-      !!
-      !!                 pond growth:      Vp = Vp + dVmelt
-      !!                    with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
-      !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
-      !!                    with Tp = -2degC
-      !!  
-      !! ** Tunable parameters : (no real expertise yet, ideas?)
+      !! ** Purpose : Compute melt pond evolution
+      !!
+      !! ** Method  : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing
+      !!              We  work with volumes and then redistribute changes into thickness and concentration
+      !!              assuming linear relationship between the two. 
+      !!
+      !! ** Action  : - pond growth:      Vp = Vp + dVmelt                                          --- from Holland et al 2012 ---
+      !!                                     dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
+      !!                                        dh_i  = meltwater from ice surface melt
+      !!                                        dh_s  = meltwater from snow melt
+      !!                                        (1-r) = fraction of melt water that is not flushed
+      !!
+      !!              - limtations:       a_ip must not exceed (1-r)*a_i
+      !!                                  h_ip must not exceed 0.5*h_i
+      !!
+      !!              - pond shrinking:
+      !!                       if lids:   Vp = Vp -dH * a_ip
+      !!                                     dH = lid thickness change. Retrieved from this eq.:    --- from Flocco et al 2010 ---
+      !!
+      !!                                                                   rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H 
+      !!                                                                      H = lid thickness
+      !!                                                                      Lf = latent heat of fusion
+      !!                                                                      Tp = -2C
+      !!
+      !!                                                                And solved implicitely as:
+      !!                                                                   H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0
+      !!
+      !!                    if no lids:   Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)                      --- from Holland et al 2012 ---
+      !!
+      !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi                 --- from Flocco et al 2007 ---
+      !!                                     perm = permability of sea-ice
+      !!                                     visc = water viscosity
+      !!                                     Hp   = height of top of the pond above sea-level
+      !!                                     Hi   = ice thickness thru which there is flushing
+      !!
+      !!              - Corrections:      remove melt ponds when lid thickness is 10 times the pond thickness
+      !!
+      !!              - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip:
+      !!                                  a_ip/a_i = a_ip_frac = h_ip / zaspect
+      !!
+      !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min
       !! 
-      !! ** Note       : Stolen from CICE for quick test of the melt pond
-      !!                 radiation and freshwater interfaces
-      !!                 Coupling can be radiative AND freshwater
-      !!                 Advection, ridging, rafting are called
-      !!
-      !! ** References : Holland, M. M. et al (J Clim 2012)
-      !!-------------------------------------------------------------------
-      REAL(wp), PARAMETER ::   zrmin       = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
-      REAL(wp), PARAMETER ::   zrmax       = 0.70_wp  ! maximum     -           -         -         -            -
-      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio
-      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
-      !
-      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
-      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding
-      REAL(wp) ::   z1_Tp            ! inverse reference temperature
-      REAL(wp) ::   z1_rhow          ! inverse freshwater density
-      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
-      REAL(wp) ::   zfac, zdum
-      !
-      INTEGER  ::   ji   ! loop indices
-      !!-------------------------------------------------------------------
-      z1_rhow        = 1._wp / rhow 
-      z1_zpnd_aspect = 1._wp / zpnd_aspect
-      z1_Tp          = 1._wp / zTp 
+      !! ** Note       :   mostly stolen from CICE
+      !!
+      !! ** References :   Flocco and Feltham (JGR, 2007)
+      !!                   Flocco et al       (JGR, 2010)
+      !!                   Holland et al      (J. Clim, 2012)
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array
+      !!
+      REAL(wp), PARAMETER ::   zaspect =  0.8_wp      ! pond aspect ratio
+      REAL(wp), PARAMETER ::   zTp     = -2._wp       ! reference temperature
+      REAL(wp), PARAMETER ::   zvisc   =  1.79e-3_wp  ! water viscosity
+      !!
+      REAL(wp) ::   zfr_mlt, zdv_mlt                  ! fraction and volume of available meltwater retained for melt ponding
+      REAL(wp) ::   zdv_frz, zdv_flush                ! Amount of melt pond that freezes, flushes
+      REAL(wp) ::   zhp                               ! heigh of top of pond lid wrt ssh
+      REAL(wp) ::   zv_ip_max                         ! max pond volume allowed
+      REAL(wp) ::   zdT                               ! zTp-t_su
+      REAL(wp) ::   zsbr                              ! Brine salinity
+      REAL(wp) ::   zperm                             ! permeability of sea ice
+      REAL(wp) ::   zfac, zdum                        ! temporary arrays
+      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse
+      !!
+      INTEGER  ::   ji, jk                            ! loop indices
+      !!-------------------------------------------------------------------
+      z1_rhow   = 1._wp / rhow 
+      z1_aspect = 1._wp / zaspect
+      z1_Tp     = 1._wp / zTp 
 
       DO ji = 1, npti
-         !                                                        !--------------------------------!
-         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
-            !                                                     !--------------------------------!
-            !--- Remove ponds on thin ice
+         !                                                            !----------------------------------------------------!
+         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction !
+            !                                                         !----------------------------------------------------!
+            !--- Remove ponds on thin ice or tiny ice fractions
             a_ip_1d(ji)      = 0._wp
-            a_ip_frac_1d(ji) = 0._wp
             h_ip_1d(ji)      = 0._wp
-            !                                                     !--------------------------------!
-         ELSE                                                     ! Case ice thickness >= rn_himin !
-            !                                                     !--------------------------------!
-            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
-            !
-            ! available meltwater for melt ponding [m, >0] and fraction
-            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
-            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
-            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper 
-            !
-            !--- Pond gowth ---!
-            ! v_ip should never be negative, otherwise code crashes
-            v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt )
-            !
-            ! melt pond mass flux (<0)
+            h_il_1d(ji)      = 0._wp
+            !
+            ! clem: problem with conservation or not ?
+            !                                                         !--------------------------------!
+         ELSE                                                         ! Case ice thickness >= rn_himin !
+            !                                                         !--------------------------------!
+            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness
+            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
+            !
+            !------------------!
+            ! case ice melting !
+            !------------------!
+            !
+            !--- available meltwater for melt ponding ---!
+            zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
+            zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) in H12 = fraction of melt water that is not flushed
+            zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 
+            !
+            !--- overflow ---!
+            ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
+            !    a_ip_max = zfr_mlt * a_i
+            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
+            zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect
+            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
+
+            ! If pond depth exceeds half the ice thickness then reduce the pond volume
+            !    h_ip_max = 0.5 * h_i
+            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
+            zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji)
+            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
+            
+            !--- Pond growing ---!
+            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
+            !
+            !--- Lid melting ---!
+            IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
+            !
+            !--- mass flux ---!
             IF( zdv_mlt > 0._wp ) THEN
-               zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice
+               zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1]
                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
                !
-               ! adjust ice/snow melting flux to balance melt pond flux (>0)
-               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
+               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux
                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
                wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
             ENDIF
+
+            !-------------------!
+            ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
+            !-------------------!
+            !
+            zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
             !
             !--- Pond contraction (due to refreezing) ---!
-            v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
-            !
-            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
-            !    h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
-            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
-            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
-            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
+            IF( ln_pnd_lids ) THEN
+               !
+               !--- Lid growing and subsequent pond shrinking ---! 
+               zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
+                  &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
+               
+               ! Lid growing
+               v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz )
+               
+               ! Pond shrinking
+               v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
+
+            ELSE
+               ! Pond shrinking
+               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
+            ENDIF
+            !
+            !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
+            ! v_ip     = h_ip * a_ip
+            ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip)
+            a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
+            h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji)
+
+            !---------------!            
+            ! Pond flushing !
+            !---------------!
+            IF( ln_pnd_flush ) THEN
+               ! height of top of the pond above sea-level
+               zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0
+
+               ! Calculate the permeability of the ice (Assur 1958)
+               DO jk = 1, nlay_i
+                  zsbr = - 1.2_wp                                  &
+                     &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    &
+                     &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 &
+                     &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 ! clem: error here the factor was 0.01878 instead of 0.0178 (cf Flocco 2010)
+                  ztmp(jk) = sz_i_1d(ji,jk) / zsbr
+               END DO
+               zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 )
+
+               ! Do the drainage using Darcy's law
+               zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)
+               zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) )
+               v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
+               
+               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
+               a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
+               h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji)
+
+            ENDIF
+
+            !--- Corrections and lid thickness ---!
+            IF( ln_pnd_lids ) THEN
+               !--- retrieve lid thickness from volume ---!
+               IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
+               ELSE                              ;   h_il_1d(ji) = 0._wp
+               ENDIF
+               !--- remove ponds if lids are much larger than ponds ---!
+               IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
+                  a_ip_1d(ji)      = 0._wp
+                  h_ip_1d(ji)      = 0._wp
+                  h_il_1d(ji)      = 0._wp
+               ENDIF
+            ENDIF
             !
          ENDIF
+         
       END DO
       !
@@ -205,5 +323,7 @@
       INTEGER  ::   ios, ioptio   ! Local integer
       !!
-      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
+      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_lids, ln_pnd_flush, rn_apnd_min, rn_apnd_max, &
+         &                          ln_pnd_CST, rn_apnd, rn_hpnd, &
+         &                          ln_pnd_alb
       !!-------------------------------------------------------------------
       !
@@ -221,10 +341,14 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namicethd_pnd:'
-         WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd
-         WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12
-         WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST
-         WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
-         WRITE(numout,*) '            Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
-         WRITE(numout,*) '         Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
+         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd
+         WRITE(numout,*) '         Evolutive  melt pond fraction and depth                  ln_pnd_H12   = ', ln_pnd_H12
+         WRITE(numout,*) '            Melt ponds can have frozen lids                       ln_pnd_lids  = ', ln_pnd_lids
+         WRITE(numout,*) '            Allow ponds to flush thru the ice                     ln_pnd_flush = ', ln_pnd_flush
+         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min
+         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max
+         WRITE(numout,*) '         Prescribed melt pond fraction and depth                  ln_pnd_CST   = ', ln_pnd_CST
+         WRITE(numout,*) '            Prescribed pond fraction                              rn_apnd      = ', rn_apnd
+         WRITE(numout,*) '            Prescribed pond depth                                 rn_hpnd      = ', rn_hpnd
+         WRITE(numout,*) '         Melt ponds affect albedo or not                          ln_pnd_alb   = ', ln_pnd_alb
       ENDIF
       !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/iceupdate.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/iceupdate.F90	(revision 12812)
@@ -94,6 +94,5 @@
       REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2)
       REAL(wp) ::   zqsr             ! New solar flux received by the ocean
-      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                  ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                  ! 2D workspace
       !!---------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('ice_update')
@@ -185,7 +184,6 @@
       ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
       !------------------------------------------------------------------
-      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
-      !
-      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, cloud_fra, alb_ice ) ! ice albedo
+
       !
       IF( lrst_ice ) THEN                       !* write snwice_mass fields in the restart file
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icevar.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icevar.F90	(revision 12812)
@@ -113,4 +113,5 @@
       at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
       vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
+      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 )
       !
       ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction  
@@ -161,6 +162,6 @@
          !
          !                           ! mean melt pond depth
-         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)
-         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp
+         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)   ;   hm_il(:,:) = vt_il(:,:) / at_ip(:,:)
+         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp                     ;   hm_il(:,:) = 0._wp
          END WHERE         
          !
@@ -184,5 +185,7 @@
       REAL(wp) ::   zhmax, z1_zhmax                 !   -      -
       REAL(wp) ::   zlay_i, zlay_s                  !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
+      REAL(wp), PARAMETER ::   zhl_max =  0.015_wp  ! pond lid thickness above which the ponds disappear from the albedo calculation
+      REAL(wp), PARAMETER ::   zhl_min =  0.005_wp  ! pond lid thickness below which the full pond area is used in the albedo calculation
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i, z1_a_ip
       !!-------------------------------------------------------------------
 
@@ -203,4 +206,8 @@
       ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
       END WHERE
+      !
+      WHERE( a_ip(:,:,:) > epsi20 )  ;   z1_a_ip(:,:,:) = 1._wp / a_ip(:,:,:)
+      ELSEWHERE                      ;   z1_a_ip(:,:,:) = 0._wp
+      END WHERE
       !                                           !--- ice thickness
       h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
@@ -217,10 +224,14 @@
       !                                           !--- ice age      
       o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
-      !                                           !--- pond fraction and thickness      
+      !                                           !--- pond and lid thickness      
+      h_ip(:,:,:) = v_ip(:,:,:) * z1_a_ip(:,:,:)
+      h_il(:,:,:) = v_il(:,:,:) * z1_a_ip(:,:,:)
+      !                                           !--- melt pond effective area (used for albedo)
       a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
-      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
-      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
-      END WHERE
-      !
+      WHERE    ( h_il(:,:,:) <= zhl_min )  ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:)       ! lid is very thin.  Expose all the pond
+      ELSEWHERE( h_il(:,:,:) >= zhl_max )  ;   a_ip_eff(:,:,:) = 0._wp                  ! lid is very thick. Cover all the pond up with ice and snow
+      ELSEWHERE                            ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * &   ! lid is in between. Expose part of the pond
+         &                                                       ( h_il(:,:,:) - zhl_min ) / ( zhl_max - zhl_min )
+      END WHERE
       !                                           !---  salinity (with a minimum value imposed everywhere)     
       IF( nn_icesal == 2 ) THEN
@@ -289,4 +300,5 @@
       sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
       v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
+      v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
       !
    END SUBROUTINE ice_var_eqv2glo
@@ -533,4 +545,5 @@
                a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
                v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
+               v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj)
                !
             END DO
@@ -555,5 +568,5 @@
 
 
-   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 )
+   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_zapneg ***
@@ -570,4 +583,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -636,10 +650,10 @@
       WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
       WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
-      !                                                        but it does not change conservation, so keep it this way is ok
+      WHERE( pv_il (:,:,:) < 0._wp )   pv_il (:,:,:) = 0._wp !    but it does not change conservation, so keep it this way is ok
       !
    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 )
+   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_roundoff ***
@@ -654,4 +668,5 @@
       REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -668,4 +683,7 @@
          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
+         IF( ln_pnd_lids ) THEN
+            WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:)   = 0._wp   ! v_il must be >= 0
+         ENDIF
       ENDIF
       !
@@ -786,6 +804,6 @@
    !! ** Purpose :  converting N-cat ice to jpl ice categories
    !!-------------------------------------------------------------------
-   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting 1-cat ice to 1 ice category
@@ -793,6 +811,6 @@
       REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
-      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !!-------------------------------------------------------------------
       ! == thickness and concentration == !
@@ -808,9 +826,10 @@
       pa_ip(:) = patip(:)
       ph_ip(:) = phtip(:)
+      ph_il(:) = phtil(:)
       
    END SUBROUTINE ice_var_itd_1c1c
 
-   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting N-cat ice to 1 ice category
@@ -818,6 +837,6 @@
       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
-      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs
@@ -854,6 +873,10 @@
       ! == ponds == !
       pa_ip(:) = SUM( patip(:,:), dim=2 )
-      WHERE( pa_ip(:) /= 0._wp )   ;   ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
-      ELSEWHERE                    ;   ph_ip(:) = 0._wp
+      WHERE( pa_ip(:) /= 0._wp )
+         ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
+         ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
+      ELSEWHERE
+         ph_ip(:) = 0._wp
+         ph_il(:) = 0._wp
       END WHERE
       !
@@ -862,6 +885,6 @@
    END SUBROUTINE ice_var_itd_Nc1c
    
-   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !!
@@ -885,6 +908,6 @@
       REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti
@@ -997,10 +1020,19 @@
          END WHERE
       END DO
+      ! keep the same v_il/v_i ratio for each category
+      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) )
+      ELSEWHERE                                 ;   zfra(:) = 0._wp
+      END WHERE
+      DO jl = 1, jpl
+         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
+         ELSEWHERE                       ;   ph_il(:,jl) = 0._wp
+         END WHERE
+      END DO
       DEALLOCATE( zfra )
       !
    END SUBROUTINE ice_var_itd_1cMc
 
-   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !!
@@ -1033,6 +1065,6 @@
       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !
       INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2
@@ -1063,4 +1095,5 @@
          pa_ip(:,:) = patip(:,:)
          ph_ip(:,:) = phtip(:,:)
+         ph_il(:,:) = phtil(:,:)
          !                              ! ---------------------- !
       ELSEIF( icat == 1 ) THEN          ! input cat = 1          !
@@ -1068,6 +1101,6 @@
          CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
             &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), &
-            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), &
-            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  )
+            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), &
+            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:)  )
          !                              ! ---------------------- !
       ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
@@ -1075,6 +1108,6 @@
          CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
             &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), &
-            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), &
-            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  )
+            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), &
+            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1)  )
          !                              ! ----------------------- !
       ELSE                              ! input cat /= output cat !
@@ -1218,4 +1251,15 @@
             END WHERE
          END DO
+         ! keep the same v_il/v_i ratio for each category
+         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
+            zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
+         ELSEWHERE
+            zfra(:) = 0._wp
+         END WHERE
+         DO jl = 1, jpl
+            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
+            ELSEWHERE                       ;   ph_il(:,jl) = 0._wp
+            END WHERE
+         END DO
          DEALLOCATE( zfra )
          !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icewri.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/ICE/icewri.F90	(revision 12812)
@@ -116,4 +116,6 @@
       IF( iom_use('icehpnd' ) )   CALL iom_put( 'icehpnd', hm_ip  * zmsk00      )                                           ! melt pond depth
       IF( iom_use('icevpnd' ) )   CALL iom_put( 'icevpnd', vt_ip  * zmsk00      )                                           ! melt pond total volume per unit area
+      IF( iom_use('icehlid' ) )   CALL iom_put( 'icehlid', hm_il  * zmsk00      )                                           ! melt pond lid depth
+      IF( iom_use('icevlid' ) )   CALL iom_put( 'icevlid', vt_il  * zmsk00      )                                           ! melt pond lid total volume per unit area
       ! salt
       IF( iom_use('icesalt' ) )   CALL iom_put( 'icesalt', sm_i                 * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! mean ice salinity
@@ -162,6 +164,8 @@
       IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume
       IF( iom_use('iceapnd_cat' ) )   CALL iom_put( 'iceapnd_cat' ,   a_ip         * zmsk00l                                   ) ! melt pond frac for categories
-      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond frac for categories
+      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories
+      IF( iom_use('icehlid_cat' ) )   CALL iom_put( 'icehlid_cat' ,   h_il         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond lid thickness for categories
       IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories
+      IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories
       IF( iom_use('icealb_cat'  ) )   CALL iom_put( 'icealb_cat'  ,   alb_ice      * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories
 
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/BDY/bdy_oce.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/BDY/bdy_oce.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/BDY/bdy_oce.F90	(revision 12812)
@@ -63,4 +63,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::  aip    !: now ice  pond concentration
       REAL(wp), POINTER, DIMENSION(:,:) ::  hip    !: now ice  pond depth
+      REAL(wp), POINTER, DIMENSION(:,:) ::  hil    !: now ice  pond lid depth
 #if defined key_top
       CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply
@@ -115,4 +116,5 @@
    REAL(wp), DIMENSION(jp_bdy) ::   rice_apnd               !: pond conc.  of incoming sea ice
    REAL(wp), DIMENSION(jp_bdy) ::   rice_hpnd               !: pond thick. of incoming sea ice
+   REAL(wp), DIMENSION(jp_bdy) ::   rice_hlid               !: pond lid thick. of incoming sea ice
    !
    !!----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/BDY/bdydta.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/BDY/bdydta.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/BDY/bdydta.F90	(revision 12812)
@@ -43,5 +43,5 @@
    PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90
 
-   INTEGER , PARAMETER ::   jpbdyfld  = 16    ! maximum number of files to read 
+   INTEGER , PARAMETER ::   jpbdyfld  = 17    ! maximum number of files to read 
    INTEGER , PARAMETER ::   jp_bdyssh = 1     ! 
    INTEGER , PARAMETER ::   jp_bdyu2d = 2     ! 
@@ -60,4 +60,5 @@
    INTEGER , PARAMETER ::   jp_bdyaip = 15    ! 
    INTEGER , PARAMETER ::   jp_bdyhip = 16    ! 
+   INTEGER , PARAMETER ::   jp_bdyhil = 17    ! 
 #if ! defined key_si3
    INTEGER , PARAMETER ::   jpl = 1
@@ -197,4 +198,5 @@
                         dta_bdy(jbdy)%aip(ib,jl) =  a_ip(ii,ij,jl) * tmask(ii,ij,1) 
                         dta_bdy(jbdy)%hip(ib,jl) =  h_ip(ii,ij,jl) * tmask(ii,ij,1) 
+                        dta_bdy(jbdy)%hil(ib,jl) =  h_il(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -302,4 +304,5 @@
                &                                                                         bf_alias(jp_bdya_i)%fnow(:,1,:)     !   ( a_ip = rice_apnd * a_i )
             IF( TRIM(bf_alias(jp_bdyhip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhip)%fnow(:,1,:) = rice_hpnd(jbdy)
+            IF( TRIM(bf_alias(jp_bdyhil)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhil)%fnow(:,1,:) = rice_hlid(jbdy)
             ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
             IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) &
@@ -319,4 +322,8 @@
                bf_alias(jp_bdyaip)%fnow(:,1,:) = 0._wp
                bf_alias(jp_bdyhip)%fnow(:,1,:) = 0._wp
+               bf_alias(jp_bdyhil)%fnow(:,1,:) = 0._wp
+            ENDIF
+            IF ( .NOT.ln_pnd_lids ) THEN
+               bf_alias(jp_bdyhil)%fnow(:,1,:) = 0._wp
             ENDIF
             
@@ -324,12 +331,12 @@
             ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3)            
             IF( ipl /= jpl ) THEN      ! ice: convert N-cat fields (input) into jpl-cat (output)
-               CALL ice_var_itd( bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), &
-                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , &
-                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &
-                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &
-                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), &
-                  &              dta_alias%t_i                  , dta_alias%t_s                  , &
-                  &              dta_alias%tsu                  , dta_alias%s_i                  , &
-                  &              dta_alias%aip                  , dta_alias%hip )
+               CALL ice_var_itd( bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), & ! in
+                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , & ! out
+                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &                                  ! in (optional)
+                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &                                  ! in     -
+                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), bf_alias(jp_bdyhil)%fnow(:,1,:), & ! in     -
+                  &              dta_alias%t_i                  , dta_alias%t_s                  , &                                  ! out    -
+                  &              dta_alias%tsu                  , dta_alias%s_i                  , &                                  ! out    -
+                  &              dta_alias%aip                  , dta_alias%hip                  , dta_alias%hil )                    ! out    -
             ENDIF
          ENDIF
@@ -378,5 +385,5 @@
       !                                                         ! =F => baroclinic velocities in 3D boundary data
       LOGICAL                                ::   ln_zinterp    ! =T => requires a vertical interpolation of the bdydta
-      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd 
+      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid
       INTEGER                                ::   ipk,ipl       !
       INTEGER                                ::   idvar         ! variable ID
@@ -390,11 +397,11 @@
       TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill
       TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
-      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip       
+      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip, bn_hil       
       TYPE(FLD_N), DIMENSION(:), POINTER ::   bn_alias                        ! must be an array to be used with fld_fill
       TYPE(FLD  ), DIMENSION(:), POINTER ::   bf_alias
       !
       NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 
-      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip
-      NAMELIST/nambdy_dta/ rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd
+      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip, bn_hil
+      NAMELIST/nambdy_dta/ rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid
       NAMELIST/nambdy_dta/ ln_full_vel, ln_zinterp
       !!---------------------------------------------------------------------------
@@ -452,6 +459,9 @@
 #if defined key_si3
          IF( .NOT.ln_pnd ) THEN
-            rn_ice_apnd = 0. ; rn_ice_hpnd = 0.
-            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 when no ponds' )
+            rn_ice_apnd = 0. ; rn_ice_hpnd = 0. ; rn_ice_hlid = 0.
+            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 & rn_ice_hlid = 0 when no ponds' )
+         ENDIF
+         IF( .NOT.ln_pnd_lids ) THEN
+            rn_ice_hlid = 0.
          ENDIF
 #endif
@@ -463,5 +473,6 @@
          rice_apnd(jbdy) = rn_ice_apnd
          rice_hpnd(jbdy) = rn_ice_hpnd
-         
+         rice_hlid(jbdy) = rn_ice_hlid
+
          
          DO jfld = 1, jpbdyfld
@@ -562,5 +573,5 @@
             IF(  jfld == jp_bdya_i .OR. jfld == jp_bdyh_i .OR. jfld == jp_bdyh_s .OR. &
                & jfld == jp_bdyt_i .OR. jfld == jp_bdyt_s .OR. jfld == jp_bdytsu .OR. &
-               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip      ) THEN
+               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip .OR. jfld == jp_bdyhil ) THEN
                igrd = 1                                                    ! T point
                ipk = ipl                                                   ! jpl-cat data
@@ -613,4 +624,9 @@
                bf_alias => bf(jp_bdyhip,jbdy:jbdy)                         ! alias for hip structure of bdy number jbdy
                bn_alias => bn_hip                                          ! alias for hip structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyhil ) THEN
+               cl3 = 'hil'
+               bf_alias => bf(jp_bdyhil,jbdy:jbdy)                         ! alias for hil structure of bdy number jbdy
+               bn_alias => bn_hil                                          ! alias for hil structure of nambdy_dta 
             ENDIF
 
@@ -681,4 +697,9 @@
                   ENDIF
                ENDIF
+               IF( jfld == jp_bdyhil ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%hil => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%hil(iszdim,jpl) )
+                  ENDIF
+               ENDIF
             ENDIF
 
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/BDY/bdyice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/BDY/bdyice.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/BDY/bdyice.F90	(revision 12812)
@@ -94,8 +94,8 @@
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
             ! exchange 3d arrays
-            CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. &
-                 &                      , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. &
-                 &                      , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.                &
-                 &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      )
+            CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1.                 &
+               &                        , s_i , 'T', 1., t_su, 'T', 1., v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.  &
+               &                        , a_ip, 'T', 1., v_ip, 'T', 1., v_il, 'T', 1.                                &
+               &                        , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
             ! exchange 4d arrays :   third dimension = 1   and then   third dimension = jpk
             CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
@@ -163,4 +163,5 @@
             a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration
             h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth
+            h_il(ji,jj,  jl) = ( h_il(ji,jj,  jl) * zwgt1 + dta%hil(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond lid depth
             !
             sz_i(ji,jj,:,jl) = s_i(ji,jj,jl)
@@ -170,4 +171,9 @@
                a_ip(ji,jj,jl) = 0._wp
                h_ip(ji,jj,jl) = 0._wp
+               h_il(ji,jj,jl) = 0._wp
+            ENDIF
+
+            IF( .NOT.ln_pnd_lids ) THEN
+               h_il(ji,jj,jl) = 0._wp
             ENDIF
             !
@@ -231,4 +237,5 @@
                a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl)
                h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl)
+               h_il(ji,jj,  jl) = h_il(ib,jb,  jl)
                !
                sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl)
@@ -265,10 +272,6 @@
                !
                ! melt ponds
-               IF( a_i(ji,jj,jl) > epsi10 ) THEN
-                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i (ji,jj,jl)
-               ELSE
-                  a_ip_frac(ji,jj,jl) = 0._wp
-               ENDIF
                v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl)
+               v_il(ji,jj,jl) = h_il(ji,jj,jl) * a_ip(ji,jj,jl)
                !
             ELSE   ! no ice at the boundary
@@ -278,14 +281,11 @@
                h_s (ji,jj,  jl) = 0._wp
                oa_i(ji,jj,  jl) = 0._wp
-               a_ip(ji,jj,  jl) = 0._wp
-               v_ip(ji,jj,  jl) = 0._wp
                t_su(ji,jj,  jl) = rt0
                t_s (ji,jj,:,jl) = rt0
                t_i (ji,jj,:,jl) = rt0 
 
-               a_ip_frac(ji,jj,jl) = 0._wp
-               h_ip     (ji,jj,jl) = 0._wp
-               a_ip     (ji,jj,jl) = 0._wp
-               v_ip     (ji,jj,jl) = 0._wp
+               a_ip(ji,jj,jl) = 0._wp
+               h_ip(ji,jj,jl) = 0._wp
+               h_il(ji,jj,jl) = 0._wp
                
                IF( nn_icesal == 1 ) THEN     ! if constant salinity
@@ -303,4 +303,6 @@
                e_s (ji,jj,:,jl) = 0._wp
                e_i (ji,jj,:,jl) = 0._wp
+               v_ip(ji,jj,  jl) = 0._wp
+               v_il(ji,jj,  jl) = 0._wp
 
             ENDIF
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/LBC/lbc_lnk_multi_generic.h90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/LBC/lbc_lnk_multi_generic.h90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/LBC/lbc_lnk_multi_generic.h90	(revision 12812)
@@ -15,26 +15,30 @@
 #endif
 
-   SUBROUTINE ROUTINE_MULTI( cdname                                                                             &
-      &                    , pt1, cdna1, psgn1, pt2 , cdna2 , psgn2 , pt3 , cdna3 , psgn3 , pt4, cdna4, psgn4   &
-      &                    , pt5, cdna5, psgn5, pt6 , cdna6 , psgn6 , pt7 , cdna7 , psgn7 , pt8, cdna8, psgn8   &
-      &                    , pt9, cdna9, psgn9, pt10, cdna10, psgn10, pt11, cdna11, psgn11                      &
+   SUBROUTINE ROUTINE_MULTI( cdname                                                                               &
+      &                    , pt1 , cdna1 , psgn1 , pt2 , cdna2 , psgn2 , pt3 , cdna3 , psgn3 , pt4 , cdna4 , psgn4   &
+      &                    , pt5 , cdna5 , psgn5 , pt6 , cdna6 , psgn6 , pt7 , cdna7 , psgn7 , pt8 , cdna8 , psgn8   &
+      &                    , pt9 , cdna9 , psgn9 , pt10, cdna10, psgn10, pt11, cdna11, psgn11, pt12, cdna12, psgn12  &
+      &                    , pt13, cdna13, psgn13, pt14, cdna14, psgn14, pt15, cdna15, psgn15, pt16, cdna16, psgn16  &
       &                    , kfillmode, pfillval, lsend, lrecv, ihlcom )
       !!---------------------------------------------------------------------
-      CHARACTER(len=*)   ,                   INTENT(in   ) :: cdname  ! name of the calling subroutine
-      ARRAY_TYPE(:,:,:,:)          , TARGET, INTENT(inout) :: pt1     ! arrays on which the lbc is applied
-      ARRAY_TYPE(:,:,:,:), OPTIONAL, TARGET, INTENT(inout) :: pt2  , pt3  , pt4  , pt5  , pt6  , pt7  , pt8  , pt9  , pt10  , pt11
-      CHARACTER(len=1)                     , INTENT(in   ) :: cdna1   ! nature of pt2D. array grid-points
-      CHARACTER(len=1)   , OPTIONAL        , INTENT(in   ) :: cdna2, cdna3, cdna4, cdna5, cdna6, cdna7, cdna8, cdna9, cdna10, cdna11
-      REAL(wp)                             , INTENT(in   ) :: psgn1   ! sign used across the north fold
-      REAL(wp)           , OPTIONAL        , INTENT(in   ) :: psgn2, psgn3, psgn4, psgn5, psgn6, psgn7, psgn8, psgn9, psgn10, psgn11
-      INTEGER            , OPTIONAL        , INTENT(in   ) :: kfillmode   ! filling method for halo over land (default = constant)
-      REAL(wp)           , OPTIONAL        , INTENT(in   ) :: pfillval    ! background value (used at closed boundaries)
-      LOGICAL, DIMENSION(4), OPTIONAL      , INTENT(in   ) :: lsend, lrecv   ! indicate how communications are to be carried out
-      INTEGER            , OPTIONAL        , INTENT(in   ) :: ihlcom         ! number of ranks and rows to be communicated
+      CHARACTER(len=*)     ,                   INTENT(in   ) ::   cdname  ! name of the calling subroutine
+      ARRAY_TYPE(:,:,:,:)            , TARGET, INTENT(inout) ::   pt1     ! arrays on which the lbc is applied
+      ARRAY_TYPE(:,:,:,:)  , OPTIONAL, TARGET, INTENT(inout) ::   pt2   , pt3   , pt4   , pt5   , pt6   , pt7   , pt8   , pt9  , &
+         &                                                        pt10  , pt11  , pt12  , pt13  , pt14  , pt15  , pt16
+      CHARACTER(len=1)                       , INTENT(in   ) ::   cdna1   ! nature of pt2D. array grid-points
+      CHARACTER(len=1)     , OPTIONAL        , INTENT(in   ) ::   cdna2 , cdna3 , cdna4 , cdna5 , cdna6 , cdna7 , cdna8 , cdna9, &
+         &                                                        cdna10, cdna11, cdna12, cdna13, cdna14, cdna15, cdna16
+      REAL(wp)                               , INTENT(in   ) ::   psgn1   ! sign used across the north fold
+      REAL(wp)             , OPTIONAL        , INTENT(in   ) ::   psgn2 , psgn3 , psgn4 , psgn5 , psgn6 , psgn7 , psgn8 , psgn9, &
+         &                                                        psgn10, psgn11, psgn12, psgn13, psgn14, psgn15, psgn16
+      INTEGER              , OPTIONAL        , INTENT(in   ) ::   kfillmode   ! filling method for halo over land (default = constant)
+      REAL(wp)             , OPTIONAL        , INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries)
+      LOGICAL, DIMENSION(4), OPTIONAL        , INTENT(in   ) ::   lsend, lrecv   ! indicate how communications are to be carried out
+      INTEGER              , OPTIONAL        , INTENT(in   ) ::   ihlcom         ! number of ranks and rows to be communicated
       !!
       INTEGER                          ::   kfld        ! number of elements that will be attributed
-      PTR_TYPE         , DIMENSION(11) ::   ptab_ptr    ! pointer array
-      CHARACTER(len=1) , DIMENSION(11) ::   cdna_ptr    ! nature of ptab_ptr grid-points
-      REAL(wp)         , DIMENSION(11) ::   psgn_ptr    ! sign used across the north fold boundary
+      PTR_TYPE         , DIMENSION(16) ::   ptab_ptr    ! pointer array
+      CHARACTER(len=1) , DIMENSION(16) ::   cdna_ptr    ! nature of ptab_ptr grid-points
+      REAL(wp)         , DIMENSION(16) ::   psgn_ptr    ! sign used across the north fold boundary
       !!---------------------------------------------------------------------
       !
@@ -55,6 +59,11 @@
       IF( PRESENT(psgn10) )   CALL ROUTINE_LOAD( pt10, cdna10, psgn10, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
       IF( PRESENT(psgn11) )   CALL ROUTINE_LOAD( pt11, cdna11, psgn11, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn12) )   CALL ROUTINE_LOAD( pt12, cdna12, psgn12, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn13) )   CALL ROUTINE_LOAD( pt13, cdna13, psgn13, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn14) )   CALL ROUTINE_LOAD( pt14, cdna14, psgn14, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn15) )   CALL ROUTINE_LOAD( pt15, cdna15, psgn15, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn16) )   CALL ROUTINE_LOAD( pt16, cdna16, psgn16, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
       !
-      CALL lbc_lnk_ptr    ( cdname,               ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv, ihlcom )
+      CALL lbc_lnk_ptr( cdname, ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv, ihlcom )
       !
    END SUBROUTINE ROUTINE_MULTI
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/SBC/sbc_ice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/SBC/sbc_ice.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/SBC/sbc_ice.F90	(revision 12812)
@@ -71,4 +71,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz         !: wind speed module at T-point                 [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tsfc_ice       !: sea ice surface skin temperature (on categories)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   cloud_fra      !: cloud cover                                    [-]
 #endif
 
@@ -90,9 +91,10 @@
    ! variables used in the coupled interface
    INTEGER , PUBLIC, PARAMETER ::   jpl = ncat
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice          ! jpi, jpj
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice 
    
    ! already defined in ice.F90 for SI3
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  h_i, h_s
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i_last_couple   !: Sea ice fraction on categories at the last coupling point
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tatm_ice       !: air temperature [K]
@@ -132,5 +134,6 @@
          &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce   (jpi,jpj)     ,   &
          &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce    (jpi,jpj)     ,   &
-         &      emp_ice (jpi,jpj)     , tsfc_ice (jpi,jpj,jpl) , sstfrz     (jpi,jpj)     , STAT= ierr(2) )
+         &      emp_ice (jpi,jpj)     , tsfc_ice (jpi,jpj,jpl) , sstfrz     (jpi,jpj)     ,   &
+         &      cloud_fra(jpi,jpj)    , STAT= ierr(2) )
 #endif
 
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/SBC/sbcblk.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/SBC/sbcblk.F90	(revision 12812)
@@ -80,5 +80,5 @@
    REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
 
-   INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read
+   INTEGER , PARAMETER ::   jpfld   =11           ! maximum number of files to read
    INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
    INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
@@ -90,5 +90,6 @@
    INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
    INTEGER , PARAMETER ::   jp_slp  = 9           ! index of sea level pressure              (Pa)
-   INTEGER , PARAMETER ::   jp_tdif =10           ! index of tau diff associated to HF tau   (N/m2)   at T-point
+   INTEGER , PARAMETER ::   jp_cc   =10           ! index of cloud cover                     (-)      range:0-1
+   INTEGER , PARAMETER ::   jp_tdif =11           ! index of tau diff associated to HF tau   (N/m2)   at T-point
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
@@ -161,5 +162,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument
+      INTEGER  ::   jfpr, jfld            ! dummy loop indice and argument
       INTEGER  ::   ios, ierror, ioptio   ! Local integer
       !!
@@ -168,7 +169,7 @@
       TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
       TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        "
-      TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        "
+      TYPE(FLD_N) ::   sn_slp , sn_tdif, sn_cc                 !       "                        "
       NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
-         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
+         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, sn_cc,         &
          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
@@ -214,5 +215,6 @@
       slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
       slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
-      slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif
+      slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_cc)   = sn_cc
+      slf_i(jp_tdif) = sn_tdif
       !
       lhftau = ln_taudif                     !- add an extra field if HF stress is used
@@ -222,15 +224,24 @@
       ALLOCATE( sf(jfld), STAT=ierror )
       IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
-      DO ifpr= 1, jfld
-         ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
-         IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
-         IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   &
-            &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
-            &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
-
-      END DO
+
       !                                      !- fill the bulk structure with namelist informations
       CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
       !
+      DO jfpr = 1, jfld
+         !
+         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to zero)
+            ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
+            sf(jfpr)%fnow(:,:,1) = 0._wp
+         ELSE                                                  !-- used field --!
+            ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
+            IF( slf_i(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) )
+            IF( slf_i(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(jfpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )                      &
+               &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', &
+               &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
+         ENDIF
+      ENDDO
+      ! fill cloud cover array with constant value if "not used"
+      IF( TRIM(sf(jp_cc)%clrootname) == 'NOT USED' )   sf(jp_cc)%fnow(:,:,1) = cldf_ice
+         
       IF ( ln_wave ) THEN
       !Activated wave module but neither drag nor stokes drift activated
@@ -792,6 +803,5 @@
       REAL(wp) ::   zst3                     ! local variable
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
-      REAL(wp) ::   zztmp, z1_rLsub           !   -      -
-      REAL(wp) ::   zfr1, zfr2               ! local variables
+      REAL(wp) ::   zztmp, z1_rLsub          !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
@@ -801,4 +811,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
       REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
       !!---------------------------------------------------------------------
       !
@@ -902,15 +913,19 @@
       END DO
 
+      ! --- cloud cover --- !
+      cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1)
+      
       ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
-      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
-      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
-      !
-      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
-         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
-      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
-         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
-      ELSEWHERE                                                         ! zero when hs>0
-         qtr_ice_top(:,:,:) = 0._wp 
-      END WHERE
+      ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm
+      !
+      DO jl = 1, jpl
+         WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )     ! linear decrease from hi=0 to 10cm  
+            qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
+         ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )     ! constant (ztri) when hi>10cm
+            qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
+         ELSEWHERE                                                         ! zero when hs>0
+            qtr_ice_top(:,:,jl) = 0._wp 
+         END WHERE
+      ENDDO
       !
       IF(ln_ctl) THEN
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/SBC/sbccpl.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/src/OCE/SBC/sbccpl.F90	(revision 12812)
@@ -48,4 +48,8 @@
    USE lib_mpp        ! distribued memory computing library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+
+#if defined key_oasis3 
+   USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut 
+#endif 
 
    IMPLICIT NONE
@@ -152,5 +156,5 @@
    INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level 
    INTEGER, PARAMETER ::   jps_fice1  = 33   ! first-order ice concentration (for semi-implicit coupling of atmos-ice fluxes)
-   INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area
+   INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area fraction
    INTEGER, PARAMETER ::   jps_ht_p   = 35   ! meltpond thickness
    INTEGER, PARAMETER ::   jps_kice   = 36   ! sea ice effective conductivity
@@ -159,4 +163,12 @@
 
    INTEGER, PARAMETER ::   jpsnd      = 38   ! total number of fields sent 
+
+#if ! defined key_oasis3 
+   ! Dummy variables to enable compilation when oasis3 is not being used 
+   INTEGER                    ::   OASIS_Sent        = -1 
+   INTEGER                    ::   OASIS_SentOut     = -1 
+   INTEGER                    ::   OASIS_ToRest      = -1 
+   INTEGER                    ::   OASIS_ToRestOut   = -1 
+#endif 
 
    !                                  !!** namelist namsbc_cpl **
@@ -184,4 +196,6 @@
    LOGICAL     ::   ln_usecplmask         !  use a coupling mask file to merge data received from several models
                                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   LOGICAL     ::   ln_scale_ice_flux     !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration) 
+
    TYPE ::   DYNARR     
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3   
@@ -248,14 +262,14 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos
       !!
-      NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   & 
+      NAMELIST/namsbc_cpl/  nn_cplmodel  , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux,             &
+         &                  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2   ,  & 
          &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  & 
-         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   & 
-         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr  ,   & 
+         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc ,  & 
+         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr   ,  & 
          &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_tauwoc,  &
-         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal  ,   &
-         &                  sn_rcv_iceflx, sn_rcv_co2   , nn_cplmodel , ln_usecplmask, sn_rcv_mslp ,   &
-         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq , sn_rcv_tauw, nn_cats_cpl  ,   &
+         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal   ,  &
+         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                &
+         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 &
          &                  sn_rcv_ts_ice
-
       !!---------------------------------------------------------------------
       !
@@ -279,4 +293,8 @@
       ENDIF
       IF( lwp .AND. ln_cpl ) THEN                        ! control print
+         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
+         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
+         WRITE(numout,*)'  ln_scale_ice_flux                   = ', ln_scale_ice_flux
+         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
          WRITE(numout,*)'  received fields (mutiple ice categogies)'
          WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
@@ -327,7 +345,4 @@
          WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor 
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd 
-         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
-         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
-         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
       ENDIF
 
@@ -815,4 +830,8 @@
       END SELECT
 
+      ! Initialise ice fractions from last coupling time to zero (needed by Met-Office)
+#if defined key_si3 || defined key_cice
+       a_i_last_couple(:,:,:) = 0._wp
+#endif
       !                                                      ! ------------------------- ! 
       !                                                      !      Ice Meltponds        ! 
@@ -1639,9 +1658,10 @@
       !
       INTEGER  ::   ji, jj, jl   ! dummy loop index
-      REAL(wp) ::   ztri         ! local scalar
       REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
       REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice
       REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice    !!gm , zfrqsr_tr_i
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap_ice_total
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri, zcloud_fra
       !!----------------------------------------------------------------------
       !
@@ -1663,5 +1683,4 @@
          ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
          zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
-         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * picefr(:,:)
       CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
@@ -1675,4 +1694,43 @@
 
 #if defined key_si3
+
+      ! --- evaporation over ice (kg/m2/s) --- !
+      IF (ln_scale_ice_flux) THEN ! typically met-office requirements
+         IF (sn_rcv_emp%clcat == 'yes') THEN
+            WHERE( a_i(:,:,:) > 1.e-10 )  ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+            ELSEWHERE                     ; zevap_ice(:,:,:) = 0._wp
+            END WHERE
+            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
+            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
+            END WHERE
+         ELSE
+            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * SUM( a_i_last_couple, dim=3 ) / picefr(:,:)
+            ELSEWHERE                     ; zevap_ice(:,:,1) = 0._wp
+            END WHERE
+            zevap_ice_total(:,:) = zevap_ice(:,:,1)
+            DO jl = 2, jpl
+               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
+            ENDDO
+         ENDIF
+      ELSE
+         IF (sn_rcv_emp%clcat == 'yes') THEN
+            zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl)
+            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
+            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
+            END WHERE
+         ELSE
+            zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1)
+            zevap_ice_total(:,:) = zevap_ice(:,:,1)
+            DO jl = 2, jpl
+               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
+            ENDDO
+         ENDIF
+      ENDIF
+
+      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
+         ! For conservative case zemp_ice has not been defined yet. Do it now.
+         zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:)
+      ENDIF
+
       ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing)
       zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw )
@@ -1683,11 +1741,5 @@
 
       ! --- evaporation over ocean (used later for qemp) --- !
-      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:)
-
-      ! --- evaporation over ice (kg/m2/s) --- !
-      DO jl=1,jpl
-         IF (sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl)
-         ELSE                                  ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,1 )   ;   ENDIF
-      ENDDO
+      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
 
       ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
@@ -1784,4 +1836,7 @@
       CASE( 'oce only' )         ! the required field is directly provided
          zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
+         ! For Met Office sea ice non-solar fluxes are already delt with by JULES so setting to zero
+         ! here so the only flux is the ocean only one.
+         zqns_ice(:,:,:) = 0._wp 
       CASE( 'conservative' )     ! the required fields are directly provided
          zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
@@ -1914,4 +1969,7 @@
       CASE( 'oce only' )
          zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
+         ! For Met Office sea ice solar fluxes are already delt with by JULES so setting to zero
+         ! here so the only flux is the ocean only one.
+         zqsr_ice(:,:,:) = 0._wp
       CASE( 'conservative' )
          zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
@@ -1992,4 +2050,6 @@
             ENDDO
          ENDIF
+      CASE( 'none' ) 
+         zdqns_ice(:,:,:) = 0._wp
       END SELECT
       
@@ -2007,8 +2067,28 @@
       !                                                      ! ========================= !
       CASE ('coupled')
-         qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
-         qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
+         IF (ln_scale_ice_flux) THEN
+            WHERE( a_i(:,:,:) > 1.e-10_wp )
+               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+            ELSEWHERE
+               qml_ice(:,:,:) = 0.0_wp
+               qcn_ice(:,:,:) = 0.0_wp
+            END WHERE
+         ELSE
+            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
+            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
+         ENDIF
       END SELECT
-      !
+!!$      !                                                      ! ========================= !
+!!$      SELECT CASE( TRIM( sn_rcv_clouds%cldes ) )             !       cloud fraction      !
+!!$      !                                                      ! ========================= !
+!!$         cloud_fra(:,:) = frcv(jpr_clfra)*z3(:,:,1)
+!!$      END SELECT
+      zcloud_fra(:,:) = cldf_ice   ! should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
+      IF( ln_mixcpl ) THEN
+         cloud_fra(:,:) = cloud_fra(:,:) * xcplmask(:,:,0) + zcloud_fra(:,:)* zmsk(:,:)
+      ELSE
+         cloud_fra(:,:) = zcloud_fra(:,:)
+      ENDIF
       !                                                      ! ========================= !
       !                                                      !      Transmitted Qsr      !   [W/m2]
@@ -2017,9 +2097,16 @@
          !
          !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77)
+         !                    !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
+         ztri(:,:) = 0.18 * ( 1.0 - zcloud_fra(:,:) ) + 0.35 * zcloud_fra(:,:)  ! surface transmission when hi>10cm (Grenfell Maykut 77)
          !
-         qtr_ice_top(:,:,:) = ztri * qsr_ice(:,:,:)
-         WHERE( phs(:,:,:) >= 0.0_wp )   qtr_ice_top(:,:,:) = 0._wp            ! snow fully opaque
-         WHERE( phi(:,:,:) <= 0.1_wp )   qtr_ice_top(:,:,:) = qsr_ice(:,:,:)   ! thin ice transmits all solar radiation
+         DO jl = 1, jpl
+            WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )     ! linear decrease from hi=0 to 10cm  
+               zqtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
+            ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )     ! constant (ztri) when hi>10cm
+               zqtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
+            ELSEWHERE                                                         ! zero when hs>0
+               zqtr_ice_top(:,:,jl) = 0._wp
+            END WHERE
+         ENDDO
          !     
       ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
@@ -2193,4 +2280,16 @@
       ENDIF
 
+#if defined key_si3 || defined key_cice
+      ! If this coupling was successful then save ice fraction for use between coupling points. 
+      ! This is needed for some calculations where the ice fraction at the last coupling point 
+      ! is needed. 
+      IF(  info == OASIS_Sent    .OR. info == OASIS_ToRest .OR. & 
+         & info == OASIS_SentOut .OR. info == OASIS_ToRestOut ) THEN 
+         IF ( sn_snd_thick%clcat == 'yes' ) THEN 
+           a_i_last_couple(:,:,1:jpl) = a_i(:,:,1:jpl)
+         ENDIF
+      ENDIF
+#endif
+
       IF( ssnd(jps_fice1)%laction ) THEN
          SELECT CASE( sn_snd_thick1%clcat )
@@ -2256,12 +2355,12 @@
             SELECT CASE( sn_snd_mpnd%clcat )  
             CASE( 'yes' )  
-               ztmp3(:,:,1:jpl) =  a_ip(:,:,1:jpl)
-               ztmp4(:,:,1:jpl) =  v_ip(:,:,1:jpl)  
+               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl)
+               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)  
             CASE( 'no' )  
                ztmp3(:,:,:) = 0.0  
                ztmp4(:,:,:) = 0.0  
                DO jl=1,jpl  
-                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip(:,:,jpl)  
-                 ztmp4(:,:,1) = ztmp4(:,:,1) + v_ip(:,:,jpl) 
+                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)
+                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl)
                ENDDO  
             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )  
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/EXPREF/file_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/EXPREF/file_def_nemo-oce.xml	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/EXPREF/file_def_nemo-oce.xml	(revision 12812)
@@ -15,6 +15,6 @@
 	  <field field_ref="soce" /> 
 	  <field field_ref="ssh"  />
-	  <field field_ref="salgrad"  />
-	  <field field_ref="ke_zint"  />
+	  <field field_ref="socegrad"  />
+	  <field field_ref="eken_int"  />
 	  <field field_ref="relvor"  />
 	  <field field_ref="potvor"  />
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/EXPREF/namelist_cfg	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/EXPREF/namelist_cfg	(revision 12812)
@@ -20,10 +20,10 @@
 &namusr_def    !   User defined :   CANAL configuration: Flat bottom, beta-plane
 !-----------------------------------------------------------------------
-   rn_domszx   =   3600.   !  x horizontal size         [km]
-   rn_domszy   =   1800.   !  y horizontal size         [km]
-   rn_domszz   =   5000.   !  z vertical size            [m]
-   rn_dx       =     30.   !  x horizontal resolution   [km]
-   rn_dy       =     30.   !  y horizontal resolution   [km]
-   rn_dz       =    500.   !  z vertical resolution      [m]
+   rn_domszx   =   2000.   !  x horizontal size         [km]
+   rn_domszy   =   1000.   !  y horizontal size         [km]
+   rn_domszz   =   1000.   !  z vertical size            [m]
+   rn_dx       =     10.   !  x horizontal resolution   [km]
+   rn_dy       =     10.   !  y horizontal resolution   [km]
+   rn_dz       =   1000.   !  z vertical resolution      [m]
    rn_0xratio  =      0.5  !  x-domain ratio of the 0
    rn_0yratio  =      0.5  !  y-domain ratio of the 0
@@ -31,16 +31,23 @@
    rn_ppgphi0  =    38.5   !  Reference latitude      [degrees]
    rn_u10      =      0.   !  10m wind speed              [m/s]
-     rn_windszx =   4000.   !  longitudinal wind extension   [km]
-     rn_windszy =   4000.   !  latitudinal wind extension    [km]
-     rn_uofac  =      0.   !  Uoce multiplicative factor (0.:absolute or 1.:relative winds)
+     rn_windszx =   90.    !  longitudinal wind extension   [km]
+     rn_windszy =   90.    !  latitudinal wind extension    [km]
+!!clem     rn_uofac  =     0.    !  Uoce multiplicative factor (0.:absolute or 1.:relative winds)
    rn_vtxmax   =      1.   !  initial vortex max current  [m/s]
    rn_uzonal   =      1.   !  initial zonal current       [m/s]
-     rn_ujetszx =   4000.   !  longitudinal jet extension   [km]
-     rn_ujetszy =   4000.   !  latitudinal jet extension    [km]
+     rn_ujetszx =   4000.  !  longitudinal jet extension   [km]
+     rn_ujetszy =   400.   !  latitudinal jet extension    [km]
    nn_botcase  =      0    !  bottom definition (0:flat, 1:bump)
-   nn_initcase =      1    !  initial condition case (0:rest, 1:zonal current, 2:current shear, 3: gaussian zonal current,
-      !                    !                          4: geostrophic zonal pulse, 5: vortex)
-   ln_sshnoise =  .false.  !  add random noise on initial ssh
-   rn_lambda   =     50.   ! gaussian lambda
+   nn_initcase =      1    !  initial condition case
+   !                       !          -1 : stratif at rest
+   !                       !           0 : rest
+   !                       !           1 : zonal current
+   !                       !           2 : current shear
+   !                       !           3 : gaussian zonal current
+   !                       !           4 : geostrophic zonal pulse
+   !                       !           5 : baroclinic vortex
+   ln_sshnoise =  .FALSE.  !  add random noise on initial ssh
+   rn_lambda   =     50.   !  gaussian lambda
+   nn_perio    = 1
 /
 !-----------------------------------------------------------------------
@@ -59,6 +66,12 @@
 !-----------------------------------------------------------------------
    ln_linssh   =  .false.  !  =T  linear free surface  ==>>  model level are fixed in time
-   rn_rdt      =   1440.   !  time step for the dynamics (and tracer if nn_acc=0)
-   rn_atfp     =   0.05    !  asselin time filter parameter
+   rn_rdt      =   1200.   !  time step for the dynamics (and tracer if nn_acc=0)
+   rn_atfp     =   0.0     !  asselin time filter parameter
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration                      (default: use namusr_def in namelist_cfg)
+!-----------------------------------------------------------------------
+   ln_write_cfg = .false.   !  (=T) create the domain configuration file
+      cn_domcfg_out = "domain_cfg" ! newly created domain configuration filename
 /
 !!======================================================================
@@ -148,8 +161,8 @@
    ln_traadv_OFF = .false. !  No tracer advection
    ln_traadv_cen = .false. !  2nd order centered scheme
-      nn_cen_h   =  4            !  =2/4, horizontal 2nd order CEN / 4th order CEN
-      nn_cen_v   =  4            !  =2/4, vertical   2nd order CEN / 4th order COMPACT
+      nn_cen_h   =  2            !  =2/4, horizontal 2nd order CEN / 4th order CEN
+      nn_cen_v   =  2            !  =2/4, vertical   2nd order CEN / 4th order COMPACT
    ln_traadv_fct = .false. !  FCT scheme
-      nn_fct_h   =  2            !  =2/4, horizontal 2nd / 4th order
+      nn_fct_h   =  4            !  =2/4, horizontal 2nd / 4th order
       nn_fct_v   =  2            !  =2/4, vertical   2nd / COMPACT 4th order
    ln_traadv_mus = .false. !  MUSCL scheme
@@ -162,5 +175,35 @@
 &namtra_ldf    !   lateral diffusion scheme for tracers                 (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_traldf_OFF   =  .true.  !  No explicit diffusion
+   !                       !  Operator type:
+   ln_traldf_OFF   = .true.    !  No explicit diffusion
+   ln_traldf_lap   = .false.   !    laplacian operator
+   ln_traldf_blp   = .false.   !  bilaplacian operator
+   !
+   !                       !  Direction of action:
+   ln_traldf_lev   = .false.   !  iso-level
+   ln_traldf_hor   = .true.    !  horizontal  (geopotential)
+   ln_traldf_iso   = .false.   !  iso-neutral (standard operator)
+   ln_traldf_triad = .false.   !  iso-neutral (triad    operator)
+   !
+   !                             !  iso-neutral options:
+   ln_traldf_msc   = .false.   !  Method of Stabilizing Correction      (both operators)
+   rn_slpmax       =  0.01     !  slope limit                           (both operators)
+   ln_triad_iso    = .false.   !  pure horizontal mixing in ML              (triad only)
+   rn_sw_triad     = 1         !  =1 switching triad ; =0 all 4 triads used (triad only)
+   ln_botmix_triad = .false.   !  lateral mixing on bottom                  (triad only)
+   !
+   !                       !  Coefficients:
+   nn_aht_ijk_t    = 31         !  space/time variation of eddy coefficient:
+      !                             !   =-20 (=-30)    read in eddy_diffusivity_2D.nc (..._3D.nc) file
+      !                             !   =  0           constant
+      !                             !   = 10 F(k)      =ldf_c1d
+      !                             !   = 20 F(i,j)    =ldf_c2d
+      !                             !   = 21 F(i,j,t)  =Treguier et al. JPO 1997 formulation
+      !                             !   = 30 F(i,j,k)  =ldf_c2d * ldf_c1d
+      !                             !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing)
+      !                        !  time invariant coefficients:  aht0 = 1/2  Ud*Ld   (lap case)
+      !                             !                           or   = 1/12 Ud*Ld^3 (blp case)
+      rn_Ud        = 0.01           !  lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30)
+      rn_Ld        = 200.e+3        !  lateral diffusive length   [m]   (nn_aht_ijk_t= 0, 10)
 /
 !!======================================================================
@@ -183,13 +226,15 @@
       nn_dynkeg     = 0       ! scheme for grad(KE): =0   C2  ;  =1   Hollingsworth correction
    ln_dynadv_cen2 = .false. !  flux form - 2nd order centered scheme
-   ln_dynadv_ubs = .true.  !  flux form - 3rd order UBS      scheme
+   ln_dynadv_ubs  = .true.  !  flux form - 3rd order UBS      scheme
 /
 !-----------------------------------------------------------------------
 &namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_dynvor_ene = .true.  !  energy conserving scheme
-   ln_dynvor_ens = .false. !  enstrophy conserving scheme
-   ln_dynvor_mix = .false. !  mixed scheme
+   ln_dynvor_ene = .false.  !  energy conserving scheme
+   ln_dynvor_ens = .false.  !  enstrophy conserving scheme
+   ln_dynvor_mix = .false.  !  mixed scheme
    ln_dynvor_een = .false.  !  energy & enstrophy scheme
+   ln_dynvor_enT = .false.  !  energy conserving scheme (T-point)
+   ln_dynvor_eeT = .true.   !  energy conserving scheme (een using e3t)
       nn_een_e3f = 0             !  e3f = masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
 /
@@ -210,5 +255,5 @@
          !                          !                     = 1 Boxcar over   nn_baro sub-steps
          !                          !                     = 2 Boxcar over 2*nn_baro  "    "
-      ln_bt_auto    = .false.    ! Number of sub-step defined from:
+      ln_bt_auto    = .true.    ! Number of sub-step defined from:
          nn_baro      =  24         ! =F : the number of sub-step in rn_rdt seconds
 /
@@ -222,8 +267,8 @@
    !                       !  Direction of action  :
    ln_dynldf_lev =  .false.    !  iso-level
-   ln_dynldf_hor =  .true.    !  horizontal (geopotential)
+   ln_dynldf_hor =  .false.    !  horizontal (geopotential)
    ln_dynldf_iso =  .false.    !  iso-neutral
    !                       !  Coefficient
-   nn_ahm_ijk_t  = 20           !  space/time variation of eddy coef
+   nn_ahm_ijk_t  = 31           !  space/time variation of eddy coef
       !                             !  =-30  read in eddy_viscosity_3D.nc file
       !                             !  =-20  read in eddy_viscosity_2D.nc file
@@ -275,7 +320,8 @@
 !!   namdiu       Cool skin and warm layer models                       (default: OFF)
 !!   namdiu       Cool skin and warm layer models                       (default: OFF)
-!!   namflo       float parameters                                      (default: OFF)
-!!   nam_diaharm  Harmonic analysis of tidal constituents               (default: OFF)
-!!   nam_diadct   transports through some sections                      (default: OFF)
+!!   namflo       float parameters                                      ("key_float")
+!!   nam_diaharm  Harmonic analysis of tidal constituents               ("key_diaharm")
+!!   namdct       transports through some sections                      ("key_diadct")
+!!   nam_diatmb   Top Middle Bottom Output                              (default: OFF)
 !!   nam_diatmb   Top Middle Bottom Output                              (default: OFF)
 !!   nam_dia25h   25h Mean Output                                       (default: OFF)
@@ -287,5 +333,5 @@
 !-----------------------------------------------------------------------
    ln_glo_trd  = .false.   ! (T) global domain averaged diag for T, T^2, KE, and PE
-   ln_dyn_trd  = .true.   ! (T) 3D momentum trend output
+   ln_dyn_trd  = .true.    ! (T) 3D momentum trend output
    ln_dyn_mxl  = .false.   ! (T) 2D momentum trends averaged over the mixed layer (not coded yet)
    ln_vor_trd  = .false.   ! (T) 2D barotropic vorticity trends (not coded yet)
@@ -314,8 +360,12 @@
 &nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi")
 !-----------------------------------------------------------------------
+!!   jpni        =   8       !  jpni   number of processors following i (set automatically if < 1)
+!!   jpnj        =   1       !  jpnj   number of processors following j (set automatically if < 1)
 /
 !-----------------------------------------------------------------------
 &namctl        !   Control prints                                       (default: OFF)
 !-----------------------------------------------------------------------
+   ln_timing   = .true.   !  timing by routine write out in timing.output file
+!!   ln_diacfl   = .true.   !  CFL diagnostics write out in cfl_diagnostics.ascii
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/diawri.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/diawri.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/diawri.F90	(revision 12812)
@@ -231,5 +231,5 @@
       IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
 
-      IF ( iom_use("salgrad") .OR. iom_use("salgrad2") ) THEN
+      IF ( iom_use("socegrad") .OR. iom_use("socegrad2") ) THEN
          z3d(:,:,jpk) = 0.
          DO jk = 1, jpkm1
@@ -245,7 +245,7 @@
          END DO
          CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
-         CALL iom_put( "salgrad2",  z3d )          ! square of module of sal gradient
+         CALL iom_put( "socegrad2",  z3d )          ! square of module of sal gradient
          z3d(:,:,:) = SQRT( z3d(:,:,:) )
-         CALL iom_put( "salgrad" ,  z3d )          ! module of sal gradient
+         CALL iom_put( "socegrad" ,  z3d )          ! module of sal gradient
       ENDIF
          
@@ -300,12 +300,12 @@
             END DO
          END DO
-         CALL iom_put( "salt2c", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
-      ENDIF
-      !
-      IF ( iom_use("eken") ) THEN
+         CALL iom_put( "salt2c", rau0 * z2d )          ! vertically integrated squared salt content (PSU*kg/m2)
+      ENDIF
+      !
+      IF ( iom_use("eken") .OR. iom_use("eken_int") ) THEN
          z3d(:,:,jpk) = 0._wp 
          DO jk = 1, jpkm1
-            DO jj = 2, jpj
-               DO ji = 2, jpi
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
                   zztmpx = 0.5 * ( un(ji-1,jj  ,jk) + un(ji,jj,jk) )
                   zztmpy = 0.5 * ( vn(ji  ,jj-1,jk) + vn(ji,jj,jk) )
@@ -316,25 +316,4 @@
          CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
          CALL iom_put( "eken", z3d )                 ! kinetic energy
-      ENDIF
-
-      IF ( iom_use("ke") .or. iom_use("ke_zint") ) THEN
-         !
-         z3d(:,:,jpk) = 0._wp
-         z3d(1,:, : ) = 0._wp
-         z3d(:,1, : ) = 0._wp
-         DO jk = 1, jpkm1
-            DO jj = 2, jpj
-               DO ji = 2, jpi
-                  z3d(ji,jj,jk) = 0.25_wp * ( un(ji  ,jj,jk) * un(ji  ,jj,jk) * e1e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)  &
-                     &                      + un(ji-1,jj,jk) * un(ji-1,jj,jk) * e1e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)  &
-                     &                      + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1e2v(ji,jj  ) * e3v_n(ji,jj  ,jk)  &
-                     &                      + vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1e2v(ji,jj-1) * e3v_n(ji,jj-1,jk)  )  &
-                     &                    * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         
-         CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
-         CALL iom_put( "ke", z3d ) ! kinetic energy
 
          z2d(:,:)  = 0._wp 
@@ -342,10 +321,9 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * z3d(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "ke_zint", z2d )   ! vertically integrated kinetic energy
-
+                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "eken_int", z2d )   ! vertically integrated kinetic energy
       ENDIF
       !
@@ -359,5 +337,5 @@
                DO ji = 1, fs_jpim1   ! vector opt.
                   z3d(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
-                     &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+                     &             - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
                END DO
             END DO
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_istate.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_istate.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_istate.F90	(revision 12812)
@@ -64,9 +64,24 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   '
       !
-      IF (ln_sshnoise) CALL RANDOM_NUMBER(zrandom)
       zjetx = ABS(rn_ujetszx)/2.
       zjety = ABS(rn_ujetszy)/2.
       !
+      zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
+      !
       SELECT CASE(nn_initcase)
+
+      CASE(-1)    ! stratif at rest
+
+         ! sea level:
+         pssh(:,:) = 0.
+         ! temperature:
+         pts(:,:,1,jp_tem) = 25. !!30._wp
+         pts(:,:,2:jpk,jp_tem) = 22. !!24._wp
+         ! salinity:  
+         pts(:,:,:,jp_sal) = 35._wp
+         ! velocities:
+         pu(:,:,:) = 0.
+         pv(:,:,:) = 0.
+
       CASE(0)    ! rest
          
@@ -96,7 +111,7 @@
             zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
             WHERE( ABS(gphit) <= zjety )
-               pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
-            ELSEWHERE
-               pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3   &
+               pssh(:,:) = - rn_uzonal / grav * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
+            ELSEWHERE
+               pssh(:,:) = - rn_uzonal / grav * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3   &
                   &                             + 0.5 * zbeta * zjety * zjety * 1.e6 )
             END WHERE
@@ -107,5 +122,10 @@
          pts(:,:,jpk,jp_sal) = 0.
          DO jk=1, jpkm1
-            pts(:,:,jk,jp_sal) = gphit(:,:)
+            WHERE( ABS(gphit) <= zjety )
+!!$            WHERE( ABS(gphit) <= zjety*0.5 .AND. ABS(glamt) <= zjety*0.5 ) ! for a square of salt
+               pts(:,:,jk,jp_sal) = 35.
+            ELSEWHERE
+               pts(:,:,jk,jp_sal) = 30.
+            END WHERE                    
          END DO
          ! velocities:
@@ -132,8 +152,8 @@
             WHERE( ABS(gphit) <= zjety )
                pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
-                  &        * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
+                  &        * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
             ELSEWHERE
                pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
-                  &        * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 )
+                  &        * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 )
             END WHERE
          END SELECT
@@ -141,7 +161,7 @@
          pts(:,:,:,jp_tem) = 10._wp
          ! salinity:  
-         pts(:,:,:,jp_sal) = 2.
-         DO jk=1, jpkm1
-            WHERE( ABS(gphiv) <= zjety ) pts(:,:,jk,jp_sal) = 2. + SIGN(1.,gphiv(:,:))
+         pts(:,:,:,jp_sal) = 30.
+         DO jk=1, jpkm1
+            WHERE( ABS(gphiv) <= zjety ) pts(:,:,jk,jp_sal) = 30. + SIGN(1.,gphiv(:,:))
          END DO
          ! velocities:
@@ -176,5 +196,5 @@
          ! salinity:  
          DO jk=1, jpkm1
-            pts(:,:,jk,jp_sal) = gphit(:,:)
+            pts(:,:,jk,jp_sal) = pssh(:,:)
          END DO
          ! velocities:
@@ -213,5 +233,5 @@
          zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
          zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
-         zlambda = SQRT(2._wp)*rn_lambda       ! Horizontal scale in meters 
+         zlambda = SQRT(2._wp)*rn_lambda*1.e3       ! Horizontal scale in meters 
          zn2 = 3.e-3**2
          zH = 0.5_wp * 5000._wp
@@ -253,6 +273,6 @@
          ! velocities:
          za = 2._wp * zP0 / zlambda**2
-         DO jj=1, jpj
-            DO ji=1, jpim1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
                zx = glamu(ji,jj) * 1.e3
                zy = gphiu(ji,jj) * 1.e3
@@ -270,6 +290,6 @@
          END DO
          !
-         DO jj=1, jpjm1
-            DO ji=1, jpi
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
                zx = glamv(ji,jj) * 1.e3
                zy = gphiv(ji,jj) * 1.e3
@@ -287,15 +307,14 @@
          END DO
          !            
+         CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
+
       END SELECT
 
       IF (ln_sshnoise) THEN
+         CALL RANDOM_SEED()
          CALL RANDOM_NUMBER(zrandom)
          pssh(:,:) = pssh(:,:) + ( 0.1  * zrandom(:,:) - 0.05 )
       END IF
-      CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. )
-      CALL lbc_lnk(  'usrdef_istate', pts, 'T',  1. )
-      CALL lbc_lnk(   'usrdef_istate', pu, 'U', -1. )
-      CALL lbc_lnk(   'usrdef_istate', pv, 'V', -1. )
-
+  
    END SUBROUTINE usr_def_istate
 
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_nam.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_nam.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_nam.F90	(revision 12812)
@@ -50,4 +50,5 @@
    LOGICAL , PUBLIC ::   ln_sshnoise=.false. ! add random noise on initial ssh
    REAL(wp), PUBLIC ::   rn_lambda  = 50.    ! gaussian lambda
+   INTEGER , PUBLIC ::   nn_perio   =    0   ! periodicity of the channel (0=closed, 1=E-W)
 
    !!----------------------------------------------------------------------
@@ -79,7 +80,7 @@
       !!
       NAMELIST/namusr_def/  rn_domszx, rn_domszy, rn_domszz, rn_dx, rn_dy, rn_dz, rn_0xratio, rn_0yratio   &
-         &                 , nn_fcase, rn_ppgphi0, rn_vtxmax, rn_uzonal, rn_ujetszx, rn_ujetszy   &
-         &                 , rn_u10, rn_windszx, rn_windszy, rn_uofac   &
-         &                 , nn_botcase, nn_initcase, ln_sshnoise, rn_lambda
+         &                 , nn_fcase, rn_ppgphi0, rn_u10, rn_windszx, rn_windszy & !!, rn_uofac   &
+         &                 , rn_vtxmax, rn_uzonal, rn_ujetszx, rn_ujetszy  &
+         &                 , nn_botcase, nn_initcase, ln_sshnoise, rn_lambda, nn_perio
       !!----------------------------------------------------------------------
       !
@@ -151,8 +152,8 @@
          WRITE(numout,*) '      add random noise on initial ssh   ln_sshnoise= ', ln_sshnoise
          WRITE(numout,*) '      Gaussian lambda parameter          rn_lambda = ', rn_lambda
-         WRITE(numout,*) '   '
-         WRITE(numout,*) '   Lateral boundary condition of the global domain'
-         WRITE(numout,*) '      EW_CANAL : closed basin               jperio = ', kperio
+         WRITE(numout,*) '      Periodicity of the basin            nn_perio = ', nn_perio
       ENDIF
+      !                             ! Set the lateral boundary condition of the global domain
+      kperio = nn_perio                    ! EW_CANAL configuration : closed basin
       !
    END SUBROUTINE usr_def_nam
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_sbc.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_sbc.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_sbc.F90	(revision 12812)
@@ -17,5 +17,5 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE phycst          ! physical constants
-   USE usrdef_nam, ONLY : rn_u10, rn_uofac, rn_windszy 
+   USE usrdef_nam, ONLY : rn_u10, rn_uofac, rn_windszy, rn_windszx 
    !
    USE in_out_manager  ! I/O manager
@@ -71,7 +71,4 @@
          !
          utau(:,:) = 0._wp
-         IF( rn_u10 /= 0. .AND. rn_windszy > 0. ) THEN
-            WHERE( ABS(gphit) <= rn_windszy/2. ) utau(:,:) = zrhocd * rn_u10 * rn_u10
-         ENDIF
          vtau(:,:) = 0._wp
          taum(:,:) = 0._wp
@@ -83,4 +80,12 @@
          qsr (:,:) = 0._wp
          !         
+      ENDIF
+
+      IF( rn_u10 /= 0. .AND. rn_windszy > 0. ) THEN
+         IF( nyear == 1 .AND. nmonth == 1 .AND. nday <= 10 ) THEN
+            WHERE( ABS(gphit) <= rn_windszy/2. .AND. ABS(glamt) <= rn_windszx/2. ) utau(:,:) = zrhocd * rn_u10 * rn_u10
+         ELSE
+            utau(:,:) = 0.
+         ENDIF
       ENDIF
 
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_zgr.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_zgr.F90	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/CANAL/MY_SRC/usrdef_zgr.F90	(revision 12812)
@@ -199,7 +199,7 @@
          zmaxlam = MAXVAL(glamt)
          CALL mpp_max( 'usrdef_zgr', zmaxlam )                 ! max over the global domain
-         zscl = rpi / zmaxlam
-         z2d(:,:) = 0.5 * ( 1. - COS( glamt(:,:) * zscl ) )
-         z2d(:,:) = REAL(jpkm1 - NINT( 0.75 * REAL(jpkm1,wp) * z2d(:,:) ), wp)
+         zscl = 0.5 * rpi / zmaxlam
+         z2d(:,:) = COS( glamt(:,:) * zscl )
+         z2d(:,:) = REAL(jpkm1 - NINT( 0.5 * REAL(jpkm1,wp) * z2d(:,:) ), wp)
       END SELECT
       !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg	(revision 12812)
@@ -88,5 +88,5 @@
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file   = .true.          !  netcdf file provided for initialization (T) or not (F)
+   nn_iceini_file   =  1              !  netcdf file provided for initialization
 
    sn_hti = 'initice_60pts'                 , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg_120pts
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg_120pts	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg_120pts	(revision 12812)
@@ -88,5 +88,5 @@
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file   = .true.          !  netcdf file provided for initialization (T) or not (F)
+   nn_iceini_file   =  1              !  netcdf file provided for initialization 
 
    sn_hti = 'initice_120pts'                 , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg_240pts
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg_240pts	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg_240pts	(revision 12812)
@@ -88,5 +88,5 @@
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file   = .true.          !  netcdf file provided for initialization (T) or not (F)
+   nn_iceini_file   =  1              !  netcdf file provided for initialization 
 
    sn_hti = 'initice_240pts'                 , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg_60pts
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg_60pts	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV1D/EXPREF/namelist_ice_cfg_60pts	(revision 12812)
@@ -88,5 +88,5 @@
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file   = .true.          !  netcdf file provided for initialization (T) or not (F)
+   nn_iceini_file   =  1              !  netcdf file provided for initialization 
 
    sn_hti = 'initice_60pts'                 , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV2D/EXPREF/namelist_ice_cfg
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV2D/EXPREF/namelist_ice_cfg	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_ADV2D/EXPREF/namelist_ice_cfg	(revision 12812)
@@ -86,5 +86,5 @@
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file   = .true.          !  netcdf file provided for initialization (T) or not (F)
+   nn_iceini_file   =  1              !  netcdf file provided for initialization 
 
    sn_hti = 'initice'                 , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
Index: /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_AGRIF/EXPREF/namelist_ice_cfg
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_AGRIF/EXPREF/namelist_ice_cfg	(revision 12811)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch/tests/ICE_AGRIF/EXPREF/namelist_ice_cfg	(revision 12812)
@@ -86,5 +86,5 @@
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file   = .true.          !  netcdf file provided for initialization (T) or not (F)
+   nn_iceini_file   =  1              !  netcdf file provided for initialization 
 
    sn_hti = 'initice'                 , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
