Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/field_def_nemo-opa.xml
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/field_def_nemo-opa.xml	(revision 3293)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/field_def_nemo-opa.xml	(revision 3293)
@@ -0,0 +1,1194 @@
+<!--
+<?xml version="1.0"?>
+       -->
+ 
+    <!-- $id$ -->
+    
+    <!-- 
+============================================================================================================
+=                                  definition of all existing variables                                    =
+=                                            DO NOT CHANGE                                                 =
+============================================================================================================
+    -->
+   <field_definition level="1" prec="4" operation="average" enabled=".TRUE." default_value="1.e20" > <!-- time step automaticaly defined -->
+
+    <!-- 
+============================================================================================================
+                                  Physical ocean and ice model variables
+============================================================================================================
+    -->
+
+      <!-- T grid -->
+      
+      <field_group id="grid_T" grid_ref="grid_T_2D" >
+         <field id="e3t"          long_name="T-cell thickness"   standard_name="cell_thickness"   unit="m"   grid_ref="grid_T_3D"/>
+         <field id="e3t_0"        long_name="Initial T-cell thickness"   standard_name="ref_cell_thickness"   unit="m"   grid_ref="grid_T_3D"/>
+
+         <field id="toce"         long_name="Sea Water Potential Temperature"         standard_name="sea_water_potential_temperature"   unit="degree_C"     grid_ref="grid_T_3D"/>
+         <field id="toce_e3t"     long_name="temperature * e3t"                                                     unit="degree_C*m"   grid_ref="grid_T_3D" > toce * e3t </field >
+         <field id="soce"         long_name="Sea Water Salinity"            standard_name="sea_water_salinity"      unit="0.001"     grid_ref="grid_T_3D"/>
+         <field id="soce_e3t"     long_name="salinity * e3t"                                                        unit="0.001*m"   grid_ref="grid_T_3D" > soce * e3t </field >
+
+         <field id="sst"          long_name="Sea Surface Temperature"             standard_name="sea_surface_temperature"             unit="degree_C"     />
+         <field id="sst2"         long_name="Square of Sea Surface Temperature"   standard_name="square_of_sea_surface_temperature"   unit="degree_C2"     > sst * sst </field >
+         <field id="sstmax"       long_name="max of sea surface temperature"   field_ref="sst"   operation="maximum"                                  />
+         <field id="sstmin"       long_name="min of sea surface temperature"   field_ref="sst"   operation="minimum"                                  />
+         <field id="sstgrad"      long_name="module of sst gradient"                                                                  unit="degree_C/m"   />
+         <field id="sstgrad2"     long_name="square of module of sst gradient"                                                        unit="degree_C2/m2" />
+         <field id="sbt"          long_name="sea bottom temperature"                                                                  unit="degree_C"     />
+         <field id="tosmint"      long_name="vertical integral of temperature times density"   standard_name="integral_wrt_depth_of_product_of_density_and_potential_temperature"  unit="(kg m2) degree_C" />
+	 
+         <field id="sss"          long_name="Sea Surface Salinity"             standard_name="sea_surface_salinity"   unit="0.001" />
+         <field id="sss2"         long_name="square of sea surface salinity"                                          unit="1e-6"  > sss * sss </field >
+         <field id="sssmax"       long_name="max of sea surface salinity"   field_ref="sss"   operation="maximum"                 />
+         <field id="sssmin"       long_name="min of sea surface salinity"   field_ref="sss"   operation="minimum"                 />
+         <field id="sbs"          long_name="sea bottom salinity"                                                     unit="0.001" />
+         <field id="somint"       long_name="vertical integral of salinity times density"   standard_name="integral_wrt_depth_of_product_of_density_and_salinity"  unit="(kg m2) x (1e-3)" /> 
+
+         <field id="taubot"       long_name="bottom stress module"                                                    unit="N/m2" /> 
+
+         <field id="ssh"          long_name="sea surface height"             standard_name="sea_surface_height_above_geoid"             unit="m" />
+         <field id="ssh2"         long_name="square of sea surface height"   standard_name="square_of_sea_surface_height_above_geoid"   unit="m2" > ssh * ssh </field >
+         <field id="sshmax"       long_name="max of sea surface height"   field_ref="ssh"   operation="maximum"                                  />
+
+         <field id="mldkz5"       long_name="Turbocline depth (Kz = 5e-4)"                       standard_name="ocean_mixed_layer_thickness_defined_by_vertical_tracer_diffusivity"                unit="m"          />
+         <field id="mldr10_1"     long_name="Mixed Layer Depth (dsigma = 0.01 wrt 10m)"          standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"          />
+         <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="maximum"                                                                          />
+         <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="minimum"                                                                          />
+         <field id="heatc"        long_name="Heat content vertically integrated"                 standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
+         <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                                   unit="0.001*kg/m2" />
+         <field id="mixlength"    long_name="Mixing Length at surface"         standard_name="mixing_length"   unit="m" />
+         <field id="mixlength2"   long_name="Mixing Length at surface"                                          unit="m2"  > mixlength * mixlength </field >
+
+         <!-- EOS -->
+         <field id="alpha"        long_name="thermal expansion"                                                         unit="degree_C-1" grid_ref="grid_T_3D" />
+         <field id="beta"         long_name="haline contraction"                                                        unit="1e3"    grid_ref="grid_T_3D" />
+         <field id="rhop"         long_name="potential density (sigma0)"        standard_name="sea_water_sigma_theta"   unit="kg/m3"  grid_ref="grid_T_3D" />
+
+         <!-- Energy - horizontal divergence -->
+         <field id="eken"         long_name="kinetic energy"          standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2"  grid_ref="grid_T_3D" />
+         <field id="hdiv"         long_name="horizontal divergence"                                                          unit="s-1"    grid_ref="grid_T_3D" />
+
+         <!-- variables available with MLE -->
+         <field id="Lf_NHpf"      long_name="MLE: Lf = N H / f"   unit="m" />
+
+         <!-- next variables available with key_diahth -->
+         <field id="mlddzt"       long_name="Thermocline Depth (depth of max dT/dz)"         standard_name="depth_at_maximum_upward_derivative_of_sea_water_potential_temperature"             unit="m"                         />
+         <field id="mldr10_3"     long_name="Mixed Layer Depth (dsigma = 0.03 wrt 10m)"      standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"                         />
+         <field id="mldr0_1"      long_name="Mixed Layer Depth (dsigma = 0.01 wrt sfc)"      standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"                         />
+         <field id="mldr0_3"      long_name="Mixed Layer Depth (dsigma = 0.03 wrt sfc)"      standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"                         />
+         <field id="mld_dt02"     long_name="Mixed Layer Depth (|dT| = 0.2 wrt 10m)"         standard_name="ocean_mixed_layer_thickness_defined_by_temperature"                                unit="m"                         />
+         <field id="topthdep"     long_name="Top of Thermocline Depth (dT = -0.2 wrt 10m)"   standard_name="ocean_mixed_layer_thickness_defined_by_temperature"                                unit="m"                         />
+         <field id="pycndep"      long_name="Pycnocline Depth (dsigma[dT=-0.2] wrt 10m)"     standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"                         />
+         <field id="BLT"          long_name="Barrier Layer Thickness"                                                                                                                          unit="m"                          > topthdep - pycndep </field>
+         <field id="tinv"         long_name="Max of vertical invertion of temperature"                                                                                                         unit="degree_C"                      />
+         <field id="depti"        long_name="Depth of max. vert. inv. of temperature"                                                                                                          unit="m"                         />
+         <field id="20d"          long_name="Depth of 20C isotherm"                          standard_name="depth_of_isosurface_of_sea_water_potential_temperature"                            unit="m"   />
+         <field id="28d"          long_name="Depth of 28C isotherm"                          standard_name="depth_of_isosurface_of_sea_water_potential_temperature"                            unit="m"   />
+         <field id="hc300"        long_name="Heat content 0-300m"                            standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"                      />
+
+         <!-- variables available with key_diaar5 -->
+         <field id="botpres"      long_name="Sea Water Pressure at Sea Floor"   standard_name="sea_water_pressure_at_sea_floor"   unit="dbar" />
+         <field id="sshdyn"       long_name="dynamic sea surface height"     standard_name="dynamic_sea_surface_height_above_geoid"     unit="m" />
+         <field id="sshdyn2"      long_name="square of dynamic sea surface height"     standard_name="dynamic_sea_surface_height_above_geoid_squared"     unit="m2" > sshdyn * sshdyn </field>
+         <field id="tnpeo"      long_name="Tendency of ocean potential energy content"          unit="W/m2"                           />
+
+         <!-- variables available with key_vvl -->
+         <field id="tpt_dep"      long_name="T-point depth"                  standard_name="depth_below_geoid"   unit="m"   grid_ref="grid_T_3D" />
+         <field id="e3tdef"       long_name="T-cell thickness deformation"                                       unit="%"   grid_ref="grid_T_3D" />
+      </field_group>
+
+      <!-- Tides -->
+
+      <field_group id="Tides_T" grid_ref="grid_T_2D" operation="once" >
+         <!-- tidal composante -->
+         <field id="M2x"          long_name="M2 Elevation harmonic real part "                             unit="m"        />
+         <field id="M2y"          long_name="M2 Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="S2x"          long_name="S2 Elevation harmonic real part "                             unit="m"        />
+         <field id="S2y"          long_name="S2 Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="N2x"          long_name="N2 Elevation harmonic real part "                             unit="m"        />
+         <field id="N2y"          long_name="N2 Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="K1x"          long_name="K1 Elevation harmonic real part "                             unit="m"        />
+         <field id="K1y"          long_name="K1 Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="O1x"          long_name="O1 Elevation harmonic real part "                             unit="m"        />
+         <field id="O1y"          long_name="O1 Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="Q1x"          long_name="Q1 Elevation harmonic real part "                             unit="m"        />
+         <field id="Q1y"          long_name="Q1 Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="M4x"          long_name="M4 Elevation harmonic real part "                             unit="m"        />
+         <field id="M4y"          long_name="M4 Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="K2x"          long_name="K2 Elevation harmonic real part "                             unit="m"        />
+         <field id="K2y"          long_name="K2 Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="P1x"          long_name="P1 Elevation harmonic real part "                             unit="m"        />
+         <field id="P1y"          long_name="P1 Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="Mfx"          long_name="Mf Elevation harmonic real part "                             unit="m"        />
+         <field id="Mfy"          long_name="Mf Elevation harmonic imaginary part"                         unit="m"        />
+         <field id="Mmx"          long_name="Mm Elevation harmonic real part "                             unit="m"        />
+         <field id="Mmy"          long_name="Mm Elevation harmonic imaginary part"                         unit="m"        />
+      </field_group>
+	 
+      <field_group id="Tides_U" grid_ref="grid_U_2D" operation="once" >
+         <field id="M2x_u"        long_name="M2 current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="M2y_u"        long_name="M2 current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="S2x_u"        long_name="S2 current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="S2y_u"        long_name="S2 current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="N2x_u"        long_name="N2 current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="N2y_u"        long_name="N2 current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="K1x_u"        long_name="K1 current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="K1y_u"        long_name="K1 current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="O1x_u"        long_name="O1 current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="O1y_u"        long_name="O1 current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="Q1x_u"        long_name="Q1 current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="Q1y_u"        long_name="Q1 current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="M4x_u"        long_name="M4 current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="M4y_u"        long_name="M4 current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="K2x_u"        long_name="K2 current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="K2y_u"        long_name="K2 current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="P1x_u"        long_name="P1 current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="P1y_u"        long_name="P1 current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="Mfx_u"        long_name="Mf current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="Mfy_u"        long_name="Mf current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="Mmx_u"        long_name="Mm current barotrope along i-axis harmonic real part "        unit="m/s"      />
+         <field id="Mmy_u"        long_name="Mm current barotrope along i-axis harmonic imaginary part "   unit="m/s"      />
+      </field_group>
+	 
+      <field_group id="Tides_V" grid_ref="grid_V_2D" operation="once" >
+         <field id="M2x_v"        long_name="M2 current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="M2y_v"        long_name="M2 current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="S2x_v"        long_name="S2 current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="S2y_v"        long_name="S2 current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="N2x_v"        long_name="N2 current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="N2y_v"        long_name="N2 current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="K1x_v"        long_name="K1 current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="K1y_v"        long_name="K1 current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="O1x_v"        long_name="O1 current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="O1y_v"        long_name="O1 current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="Q1x_v"        long_name="Q1 current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="Q1y_v"        long_name="Q1 current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="M4x_v"        long_name="M4 current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="M4y_v"        long_name="M4 current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="K2x_v"        long_name="K2 current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="K2y_v"        long_name="K2 current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="P1x_v"        long_name="P1 current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="P1y_v"        long_name="P1 current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="Mfx_v"        long_name="Mf current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="Mfy_v"        long_name="Mf current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />
+         <field id="Mmx_v"        long_name="Mm current barotrope along j-axis harmonic real part "        unit="m/s"      />
+         <field id="Mmy_v"        long_name="Mm current barotrope along j-axis harmonic imaginary part "   unit="m/s"      />	 
+      </field_group>
+
+      <!-- SBC -->
+      
+      <field_group id="SBC" grid_ref="grid_T_2D" > <!-- time step automaticaly defined based on nn_fsbc -->
+         <field id="empmr"        long_name="Net Upward Water Flux"                standard_name="water_flux_out_of_sea_ice_and_sea_water"                              unit="kg/m2/s"   />
+         <field id="empbmr"       long_name="Net Upward Water Flux at pre. tstep"  standard_name="water_flux_out_of_sea_ice_and_sea_water"                              unit="kg/m2/s"   />
+         <field id="emp_oce"      long_name="Evap minus Precip over ocean"         standard_name="evap_minus_precip_over_sea_water"                                     unit="kg/m2/s"   />
+         <field id="emp_ice"      long_name="Evap minus Precip over ice"           standard_name="evap_minus_precip_over_sea_ice"                                       unit="kg/m2/s"   />
+         <field id="saltflx"      long_name="Downward salt flux"                                                                                                        unit="0.001/m2/s" />
+         <field id="fmmflx"       long_name="Water flux due to freezing/melting"                                                                                        unit="kg/m2/s"   />
+         <field id="snowpre"      long_name="Snow precipitation"                   standard_name="snowfall_flux"                                                        unit="kg/m2/s"   />
+         <field id="runoffs"      long_name="Water Flux into Sea Water From Rivers" standard_name="water_flux_into_sea_water_from_rivers"                                unit="kg/m2/s"   />
+         <field id="precip"       long_name="Total precipitation"                  standard_name="precipitation_flux"                                                   unit="kg/m2/s"   />
+ 
+         <field id="qt"           long_name="Net Downward Heat Flux"                standard_name="surface_downward_heat_flux_in_sea_water"                              unit="W/m2"                           />
+         <field id="qns"          long_name="non solar Downward Heat Flux"                                                                                               unit="W/m2"                           />
+         <field id="qsr"          long_name="Shortwave Radiation"                   standard_name="net_downward_shortwave_flux_at_sea_water_surface"                     unit="W/m2"                           />
+         <field id="qsr3d"        long_name="Shortwave Radiation 3D distribution"   standard_name="downwelling_shortwave_flux_in_sea_water"                              unit="W/m2"      grid_ref="grid_T_3D" />
+         <field id="qrp"          long_name="Surface Heat Flux: Damping"            standard_name="heat_flux_into_sea_water_due_to_newtonian_relaxation"                 unit="W/m2"                           />
+         <field id="erp"          long_name="Surface Water Flux: Damping"           standard_name="water_flux_out_of_sea_water_due_to_newtonian_relaxation"              unit="kg/m2/s"                        />
+         <field id="taum"         long_name="wind stress module"                    standard_name="magnitude_of_surface_downward_stress"                                 unit="N/m2"                           />
+         <field id="wspd"         long_name="wind speed module"                     standard_name="wind_speed"                                                           unit="m/s"                            />
+         
+         <!-- * variable relative to atmospheric pressure forcing : available with ln_apr_dyn -->
+         <field id="ssh_ib"       long_name="Inverse barometer sea surface height"  standard_name="sea_surface_height_correction_due_to_air_pressure_at_low_frequency"   unit="m"        />
+
+         <!-- * variable related to ice shelf forcing * -->
+         <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  />
+         <field id="fwfisf3d"     long_name="Ice shelf melting"                             unit="kg/m2/s"  grid_ref="grid_T_3D" />
+         <field id="qlatisf"      long_name="Ice shelf latent heat flux"                    unit="W/m2"     />
+         <field id="qlatisf3d"    long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
+         <field id="qhcisf"       long_name="Ice shelf heat content flux"                   unit="W/m2"     />
+         <field id="qhcisf3d"     long_name="Ice shelf heat content flux"                   unit="W/m2"     grid_ref="grid_T_3D" />
+         <field id="isfgammat"    long_name="transfert coefficient for isf (temperature)"   unit="m/s"      />
+         <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)"      unit="m/s"      />
+         <field id="stbl"         long_name="salinity in the Losh tbl"                      unit="0.001"     />
+         <field id="ttbl"         long_name="temperature in the Losh tbl"                   unit="degree_C"     />
+
+         <!-- *_oce variables available with ln_blk_clio or ln_blk_core -->
+         <field id="qlw_oce"      long_name="Longwave Downward Heat Flux over open ocean"  standard_name="surface_net_downward_longwave_flux"   unit="W/m2"  />
+         <field id="qsb_oce"      long_name="Sensible Downward Heat Flux over open ocean"  standard_name="surface_downward_sensible_heat_flux"  unit="W/m2"  />
+         <field id="qla_oce"      long_name="Latent Downward Heat Flux over open ocean"    standard_name="surface_downward_latent_heat_flux"    unit="W/m2"  />
+         <field id="qemp_oce"     long_name="Downward Heat Flux from E-P over open ocean"                                                       unit="W/m2"  />
+         <field id="taum_oce"     long_name="wind stress module over open ocean"           standard_name="magnitude_of_surface_downward_stress" unit="N/m2"  />
+
+         <!-- available key_oasis3 -->
+         <field id="snow_ao_cea"  long_name="Snow over ice-free ocean (cell average)"   standard_name="snowfall_flux"                             unit="kg/m2/s"  />
+         <field id="snow_ai_cea"  long_name="Snow over sea-ice (cell average)"          standard_name="snowfall_flux"                             unit="kg/m2/s"  />
+         <field id="subl_ai_cea"  long_name="Sublimation over sea-ice (cell average)"   standard_name="surface_snow_and_ice_sublimation_flux"     unit="kg/m2/s"  />
+         <field id="icealb_cea"   long_name="Ice albedo (cell average)"                 standard_name="sea_ice_albedo"                            unit="1"        />
+         <field id="calving_cea"  long_name="Calving"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s"  />
+         <field id="iceberg_cea"  long_name="Iceberg"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s"  />
+         <field id="iceshelf_cea" long_name="Iceshelf"                                  standard_name="water_flux_into_sea_water_from_iceshelf"   unit="kg/m2/s"  />
+
+         <!-- available if key_oasis3 + conservative method -->
+         <field id="rain"          long_name="Liquid precipitation"                                     standard_name="rainfall_flux"                                                                 unit="kg/m2/s"  />
+         <field id="rain_ao_cea"   long_name="Liquid precipitation over ice-free ocean (cell average)"  standard_name="rainfall_flux"                                                                 unit="kg/m2/s"  />
+         <field id="evap_ao_cea"   long_name="Evaporation over ice-free ocean (cell average)"           standard_name="water_evaporation_flux"                                                        unit="kg/m2/s"  />
+         <field id="isnwmlt_cea"   long_name="Snow over Ice melting (cell average)"                     standard_name="surface_snow_melt_flux"                                                        unit="kg/m2/s"  />
+         <field id="fsal_virt_cea" long_name="Virtual salt flux due to ice formation (cell average)"    standard_name="virtual_salt_flux_into_sea_water_due_to_sea_ice_thermodynamics"                unit="kg/m2/s"  />
+         <field id="fsal_real_cea" long_name="Real salt flux due to ice formation (cell average)"       standard_name="downward_sea_ice_basal_salt_flux"                                              unit="kg/m2/s"  />
+         <field id="hflx_rain_cea" long_name="heat flux due to rainfall"                                standard_name="temperature_flux_due_to_rainfall_expressed_as_heat_flux_into_sea_water"        unit="W/m2"     />
+         <field id="hflx_prec_cea" long_name="heat flux due to precipitation"                           standard_name="temperature_flux_due_to_precipitation_expressed_as_heat_flux"                  unit="W/m2"     />
+         <field id="hflx_evap_cea" long_name="heat flux due to evaporation"                             standard_name="temperature_flux_due_to_evaporation_expressed_as_heat_flux_out_of_sea_water"   unit="W/m2"     />
+         <field id="hflx_snow_cea" long_name="heat flux due to snow falling"                            standard_name="heat_flux_onto_ocean_and_ice_due_to_snow_thermodynamics"                       unit="W/m2"     />
+         <field id="hflx_snow_ai_cea" long_name="heat flux due to snow falling over ice"                standard_name="heat_flux_onto_ice_due_to_snow_thermodynamics"                                 unit="W/m2"     />
+         <field id="hflx_snow_ao_cea" long_name="heat flux due to snow falling over ice-free ocean"     standard_name="heat_flux_onto_sea_water_due_to_snow_thermodynamics"                           unit="W/m2"     />
+         <field id="hflx_ice_cea"  long_name="heat flux due to ice thermodynamics"                      standard_name="heat_flux_into_sea_water_due_to_sea_ice_thermodynamics"                        unit="W/m2"     />
+         <field id="hflx_rnf_cea"  long_name="heat flux due to runoffs"                                 standard_name="temperature_flux_due_to_runoff_expressed_as_heat_flux_into_sea_water"          unit="W/m2"     />
+         <field id="hflx_cal_cea"  long_name="heat flux due to calving"                                 standard_name="heat_flux_into_sea_water_due_to_calving"                                       unit="W/m2"     />
+         <field id="hflx_icb_cea"  long_name="heat flux due to iceberg"                                 standard_name="heat_flux_into_sea_water_due_to_icebergs"                                      unit="W/m2"     />
+         <field id="hflx_isf_cea"  long_name="heat flux due to iceshelf"                                standard_name="heat_flux_into_sea_water_due_to_iceshelf"                                      unit="W/m2"     />
+         <field id="bicemel_cea"   long_name="Rate of Melt at Sea Ice Base (cell average)"              standard_name="tendency_of_sea_ice_amount_due_to_basal_melting"                               unit="kg/m2/s"  />
+         <field id="licepro_cea"   long_name="Lateral Sea Ice Growth Rate (cell average)"               standard_name="tendency_of_sea_ice_amount_due_to_lateral_growth_of_ice_floes"                 unit="kg/m2/s"  />
+         <field id="snowmel_cea"   long_name="Snow Melt Rate (cell average)"                            standard_name="surface_snow_melt_flux"                                                        unit="kg/m2/s"  />
+         <field id="sntoice_cea"   long_name="Snow-Ice Formation Rate (cell average)"                   standard_name="tendency_of_sea_ice_amount_due_to_snow_conversion"                             unit="kg/m2/s"  />
+         <field id="ticemel_cea"   long_name="Rate of Melt at Upper Surface of Sea Ice (cell average)"  standard_name="tendency_of_sea_ice_amount_due_to_surface_melting"                             unit="kg/m2/s"  />
+
+         <!-- ice fields -->
+         <field id="ice_cover"    long_name="Ice fraction"                                                 standard_name="sea_ice_area_fraction"                              unit="1"            />
+
+         <field id="qsr_ai_cea"   long_name="Air-Ice downward solar heat flux (cell average)"              standard_name="surface_downwelling_shortwave_flux_in_air"          unit="W/m2"         />
+         <field id="qns_ai_cea"   long_name="Air-Ice downward non-solar heat flux (cell average)"                                                                             unit="W/m2"         />
+         <field id="qla_ai_cea"   long_name="Air-Ice downward Latent heat flux (cell average)"             standard_name="surface_downward_latent_heat_flux"                  unit="W/m2"         />
+         
+         <field id="qsr_io_cea"   long_name="Ice-Oce downward solar heat flux (cell average)"              standard_name="net_downward_shortwave_flux_at_sea_water_surface"   unit="W/m2"         />
+         <field id="qns_io_cea"   long_name="Ice-Oce downward non-solar heat flux (cell average)"                                                                             unit="W/m2"         />
+         
+         <field id="snowthic_cea" long_name="Snow thickness (cell average)"                                standard_name="surface_snow_thickness"                             unit="m"            />
+         <field id="icethic_cea"  long_name="Ice thickness (cell average)"                                 standard_name="sea_ice_thickness"                                  unit="m"            />
+         <field id="iceprod_cea"  long_name="Ice production (cell average)"                                                                                                   unit="m/s"          />
+         <field id="iiceconc"     long_name="Ice concentration"                                            standard_name="sea_ice_area_fraction"                              unit="1"            />
+         
+         <field id="ice_pres"     long_name="Ice presence"                                                                                                                    unit=""             />
+         <field id="ist_cea"      long_name="Ice surface temperature (cell average)"                       standard_name="surface_temperature"                                unit="degree_C"         />
+         <field id="ist_ipa"      long_name="Ice surface temperature (ice presence average)"               standard_name="surface_temperature"                                unit="degree_C"         />      
+         <field id="uice_ipa"     long_name="Ice velocity along i-axis at I-point (ice presence average)"  standard_name="sea_ice_x_velocity"                                 unit="m/s"          />      
+         <field id="vice_ipa"     long_name="Ice velocity along j-axis at I-point (ice presence average)"  standard_name="sea_ice_y_velocity"                                 unit="m/s"          />      
+         
+         <field id="utau_ice"     long_name="Wind stress along i-axis over the ice at i-point"             standard_name="surface_downward_x_stress"                          unit="N/m2"         />
+         <field id="vtau_ice"     long_name="Wind stress along j-axis over the ice at i-point"             standard_name="surface_downward_y_stress"                          unit="N/m2"         />
+         
+         <field id="u_imasstr"    long_name="Sea-ice mass transport along i-axis"                          standard_name="sea_ice_x_transport"                                unit="kg/s"         />
+         <field id="v_imasstr"    long_name="Sea-ice mass transport along j-axis"                          standard_name="sea_ice_y_transport"                                unit="kg/s"         />
+         <field id="emp_x_sst"    long_name="Concentration/Dilution term on SST"                                                                                              unit="kg*degree_C/m2/s" />
+         <field id="emp_x_sss"    long_name="Concentration/Dilution term on SSS"                                                                                              unit="kg*0.001/m2/s" />        
+         <field id="rnf_x_sst"    long_name="Runoff term on SST"                                                                                                              unit="kg*degC/m2/s" />
+         <field id="rnf_x_sss"    long_name="Runoff term on SSS"                                                                                                              unit="kg*1e-3/m2/s" />
+       
+         <field id="iceconc"      long_name="ice concentration"                                            standard_name="sea_ice_area_fraction"                              unit="%"            />
+       	 <field id="isst"         long_name="sea surface temperature"                                      standard_name="sea_surface_temperature"                            unit="degree_C"         />
+         <field id="isss"         long_name="sea surface salinity"                                         standard_name="sea_surface_salinity"                               unit="0.001"         /> 
+         <field id="qt_oce"       long_name="total flux at ocean surface"                                  standard_name="surface_downward_heat_flux_in_sea_water"            unit="W/m2"         />
+         <field id="qsr_oce"      long_name="solar heat flux at ocean surface"                             standard_name="net_downward_shortwave_flux_at_sea_water_surface"   unit="W/m2"         />
+         <field id="qns_oce"      long_name="non-solar heat flux at ocean surface"                                                                                            unit="W/m2"         />
+         <field id="qt_ice"       long_name="total heat flux at ice surface: sum over categories"          standard_name="surface_downward_heat_flux_in_air"                  unit="W/m2"         />
+         <field id="qsr_ice"      long_name="solar heat flux at ice surface: sum over categories"          standard_name="surface_downwelling_shortwave_flux_in_air"          unit="W/m2"         />
+         <field id="qns_ice"      long_name="non-solar heat flux at ice surface: sum over categories"                                                                         unit="W/m2"         />
+         <field id="qtr_ice"      long_name="solar heat flux transmitted through ice: sum over categories"                                                                    unit="W/m2"         />
+         <field id="qemp_ice"     long_name="Downward Heat Flux from E-P over ice"                                                                                            unit="W/m2"         />
+         <field id="micesalt"     long_name="Mean ice salinity"                                                                                                               unit="0.001"         />
+         <field id="miceage"      long_name="Mean ice age"                                                                                                                    unit="years"        />
+         <field id="alb_ice"      long_name="Mean albedo over sea ice"                                                                                                        unit=""             />
+         <field id="albedo"       long_name="Mean albedo over sea ice and ocean"                                                                                              unit=""             />
+
+         <field id="iceage_cat"   long_name="Ice age for categories"                                       unit="days"   axis_ref="ncatice" />
+         <field id="iceconc_cat"  long_name="Ice concentration for categories"                             unit="%"      axis_ref="ncatice" />
+         <field id="icethic_cat"  long_name="Ice thickness for categories"                                 unit="m"      axis_ref="ncatice" />
+         <field id="snowthic_cat" long_name="Snow thicknessi for categories"                               unit="m"      axis_ref="ncatice" />
+         <field id="salinity_cat" long_name="Sea-Ice Bulk salinity for categories"                         unit="g/kg"   axis_ref="ncatice" />
+         <field id="brinevol_cat" long_name="Brine volume for categories"                                  unit="%"      axis_ref="ncatice" />
+         <field id="icetemp_cat"  long_name="Ice temperature for categories"                               unit="degree_C"   axis_ref="ncatice" />
+         <field id="snwtemp_cat"  long_name="Snow temperature for categories"                              unit="degree_C"   axis_ref="ncatice" />
+
+         <field id="micet"        long_name="Mean ice temperature"                                         unit="degree_C"     />
+         <field id="icehc"        long_name="ice total heat content"                                       unit="10^9J"    /> 
+         <field id="isnowhc"      long_name="snow total heat content"                                      unit="10^9J"    />
+         <field id="icest"        long_name="ice surface temperature"                                      unit="degree_C"     />
+         <field id="ibrinv"       long_name="brine volume"                                                 unit="%"        />
+         <field id="icecolf"      long_name="frazil ice collection thickness"                              unit="m"        />
+         <field id="icestr"       long_name="ice strength"                                                 unit="N/m"      />
+         <field id="icevel"       long_name="ice velocity"                                                 unit="m/s"      />
+         <field id="idive"        long_name="divergence"                                                   unit="1e-8s-1"  />
+         <field id="ishear"       long_name="shear"                                                        unit="1e-8s-1"  />
+         <field id="icevolu"      long_name="ice volume"                                                   unit="m"        />
+         <field id="snowvol"      long_name="snow volume"                                                  unit="m"        />
+
+         <field id="icetrp"       long_name="ice volume transport"                                         unit="m/day"          />
+         <field id="snwtrp"       long_name="snw volume transport"                                         unit="m/day"          />
+         <field id="saltrp"       long_name="salt content transport"                                       unit="0.001*kg/m2/day" />
+         <field id="deitrp"       long_name="advected ice enthalpy"                                        unit="W/m2"           />
+         <field id="destrp"       long_name="advected snw enthalpy"                                        unit="W/m2"           />
+
+         <field id="sfxbri"       long_name="brine salt flux"                                              unit="0.001*kg/m2/day" />
+         <field id="sfxdyn"       long_name="salt flux from ridging rafting"                               unit="0.001*kg/m2/day" />
+         <field id="sfxres"       long_name="salt flux from lipupdate (resultant)"                         unit="0.001*kg/m2/day" />
+         <field id="sfxbog"       long_name="salt flux from bot growth"                                    unit="0.001*kg/m2/day" />
+         <field id="sfxbom"       long_name="salt flux from bot melt"                                      unit="0.001*kg/m2/day" />
+         <field id="sfxsum"       long_name="salt flux from surf melt"                                     unit="0.001*kg/m2/day" />
+         <field id="sfxsni"       long_name="salt flux from snow-ice formation"                            unit="0.001*kg/m2/day" />
+         <field id="sfxopw"       long_name="salt flux from open water ice formation"                      unit="0.001*kg/m2/day" />
+         <field id="sfxsub"       long_name="salt flux from sublimation"                                   unit="0.001*kg/m2/day" />
+         <field id="sfx"          long_name="salt flux total"                                              unit="0.001*kg/m2/day" />
+
+         <field id="vfxbog"       long_name="daily bottom thermo ice prod."                                unit="m/day"   />
+         <field id="vfxdyn"       long_name="daily  dynamic ice prod."                                     unit="m/day"   />
+         <field id="vfxopw"       long_name="daily lateral thermo ice prod."                               unit="m/day"   />
+         <field id="vfxsni"       long_name="daily snowice ice prod."                                      unit="m/day"   />
+         <field id="vfxsum"       long_name="surface melt"                                                 unit="m/day"   />
+         <field id="vfxbom"       long_name="bottom melt"                                                  unit="m/day"   />
+         <field id="vfxres"       long_name="daily resultant ice prod./melting from limupdate"             unit="m/day"   />
+         <field id="vfxice"       long_name="ice melt/growth"                                              unit="m/day"   />
+         <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="m/day"   />
+         <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   />
+         <field id="vfxsub_err"   long_name="excess of snw sublimation sent to ocean"                      unit="m/day"   />
+         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   />
+         <field id="vfxthin"      long_name="daily thermo ice prod. for thin ice(<20cm) + open water"      unit="m/day"   />
+
+         <field id="afxtot"       long_name="area tendency (total)"                                        unit="day-1"   />
+         <field id="afxdyn"       long_name="area tendency (dynamics)"                                     unit="day-1"   />
+         <field id="afxthd"       long_name="area tendency (thermo)"                                       unit="day-1"   />
+
+         <field id="hfxsum"       long_name="heat fluxes causing surface ice melt"                         unit="W/m2"  />
+         <field id="hfxbom"       long_name="heat fluxes causing bottom ice melt"                          unit="W/m2"  />
+         <field id="hfxbog"       long_name="heat fluxes causing bottom ice growth"                        unit="W/m2"  />
+         <field id="hfxdif"       long_name="heat fluxes causing ice temperature change"                   unit="W/m2"  />
+         <field id="hfxopw"       long_name="heat fluxes causing open water ice formation"                 unit="W/m2"  />
+         <field id="hfxsnw"       long_name="heat fluxes causing snow melt"                                unit="W/m2"  />
+         <field id="hfxerr"       long_name="heat fluxes error after heat diffusion"                       unit="W/m2"  />
+         <field id="hfxerr_rem"   long_name="heat fluxes error after remapping"                            unit="W/m2"  />
+         <field id="hfxout"       long_name="total heat fluxes received by the ocean"                      unit="W/m2"  />
+         <field id="hfxin"        long_name="total heat fluxes at the ice/ocean surface"                   unit="W/m2"  />
+
+         <!-- heat flux associated with mass exchange -->
+         <field id="hfxthd"       long_name="heat fluxes from ice-ocean mass exchange during thermo"       unit="W/m2"  />
+         <field id="hfxdyn"       long_name="heat fluxes from ice-ocean mass exchange during dynamic"      unit="W/m2"  />
+         <field id="hfxres"       long_name="heat fluxes from ice-ocean mass exchange during resultant"    unit="W/m2"  />
+         <field id="hfxsub"       long_name="heat fluxes from ice-atm. mass exchange during sublimation"   unit="W/m2"  />
+         <field id="hfxspr"       long_name="heat fluxes from ice-atm. mass exchange during snow precip"   unit="W/m2"  />
+
+         <!-- diags -->
+         <field id="hfxdhc"       long_name="Heat content variation in snow and ice"                       unit="W/m2" />
+         <field id="hfxtur"       long_name="turbulent heat flux at the ice base"                          unit="W/m2" />
+	 <!-- sbcssm variables -->
+         <field id="sst_m"    unit="degree_C" />
+         <field id="sss_m"    unit="psu"  />
+         <field id="ssu_m"    unit="m/s"  />
+         <field id="ssv_m"    unit="m/s"  />
+         <field id="ssh_m"    unit="m"    />
+         <field id="e3t_m"    unit="m"    />
+         <field id="frq_m"    unit="-"    />
+
+      </field_group>
+
+      <!-- U grid -->
+      
+      <field_group id="grid_U"   grid_ref="grid_U_2D">
+         <field id="e3u"          long_name="U-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_U_3D" />
+         <field id="e3u_0"        long_name="Initial U-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_U_3D"/>
+         <field id="utau"         long_name="Surface Downward X Stress"                               standard_name="surface_downward_x_stress"   unit="N/m2"                            />
+         <field id="uoce"         long_name="Sea Water X Velocity"                             standard_name="sea_water_x_velocity"        unit="m/s"        grid_ref="grid_U_3D" />
+         <field id="uoce_e3u"     long_name="ocean current along i-axis * e3u"                                                                   unit="m2/s"       grid_ref="grid_U_3D"  > uoce * e3u </field>
+         <field id="ssu"          long_name="ocean surface current along i-axis"                                                                 unit="m/s"                             />
+         <field id="sbu"          long_name="ocean bottom current along i-axis"                                                                  unit="m/s"                             />
+         <field id="ubar"         long_name="ocean barotropic current along i-axis"                                                              unit="m/s"                             />
+         <field id="uocetr_eff"   long_name="Effective ocean transport along i-axis"                 standard_name="ocean_volume_x_transport"    unit="m3/s"       grid_ref="grid_U_3D" />
+         <field id="uocet"        long_name="ocean transport along i-axis times temperature (CRS)"                                               unit="degree_C*m/s"   grid_ref="grid_U_3D" />
+         <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_U_3D" />
+
+         <!-- variables available with MLE -->
+         <field id="psiu_mle"     long_name="MLE streamfunction along i-axis"   unit="m3/s"   grid_ref="grid_U_3D" />
+
+         <!-- uoce_eiv: available with key_traldf_eiv and key_diaeiv -->
+         <field id="uoce_eiv"      long_name="EIV ocean current along i-axis"   standard_name="bolus_sea_water_x_velocity"   unit="m/s"   grid_ref="grid_U_3D" />
+         <field id="ueiv_masstr"   long_name="EIV Ocean Mass X Transport"    standard_name="bolus_ocean_mass_x_transport"                          unit="kg/s"        grid_ref="grid_U_3D" />
+         <field id="ueiv_heattr3d" long_name="ocean bolus heat transport along i-axis"    standard_name="ocean_heat_x_transport_due_to_bolus_advection"   unit="W"    grid_ref="grid_U_3D" />
+         <field id="ueiv_salttr3d" long_name="ocean bolus salt transport along i-axis"    standard_name="ocean_salt_x_transport_due_to_bolus_advection"   unit="kg"   grid_ref="grid_U_3D" />
+
+         <!-- uoce_eiv: available with key_trabbl -->
+         <field id="uoce_bbl"     long_name="BBL ocean current along i-axis"    unit="m/s"  />
+         <field id="ahu_bbl"      long_name="BBL diffusive flux along i-axis"   unit="m3/s" />
+
+         <!-- variable for ice shelves -->
+         <field id="utbl"         long_name="zonal current in the Losh tbl"     unit="m/s" />
+
+         <!-- variables available with key_diaar5 -->
+         <field id="u_masstr"     long_name="Ocean Mass X Transport"    standard_name="ocean_mass_x_transport"                          unit="kg/s"        grid_ref="grid_U_3D" />
+         <field id="u_masstr_vint" long_name="vertical integral of ocean eulerian mass transport along i-axis"    standard_name="vertical_integral_of_ocean_mass_x_transport"  unit="kg/s" />
+         <field id="u_heattr"     long_name="ocean eulerian heat transport along i-axis"    standard_name="ocean_heat_x_transport"                          unit="W"                                />
+         <field id="uadv_heattr"  long_name="ocean advective heat transport along i-axis"    standard_name="advectice_ocean_heat_x_transport"                          unit="W"                                />
+         <field id="u_salttr"     long_name="ocean eulerian salt transport along i-axis"    standard_name="ocean_salt_x_transport"                          unit="0.001*kg/s"                        />
+         <field id="ueiv_heattr"  long_name="ocean bolus heat transport along i-axis"       standard_name="ocean_heat_x_transport_due_to_bolus_advection"   unit="W"                                />
+         <field id="ueiv_salttr"  long_name="ocean bolus salt transport along i-axis"       standard_name="ocean_salt_x_transport_due_to_bolus_advection"   unit="W"                                />
+         <field id="udiff_heattr" long_name="ocean diffusion heat transport along i-axis"   standard_name="ocean_heat_x_transport_due_to_diffusion"         unit="W"                                />
+      </field_group>
+      
+      <!-- V grid -->
+      
+      <field_group id="grid_V"   grid_ref="grid_V_2D">
+         <field id="e3v"          long_name="V-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_V_3D" />
+         <field id="e3v_0"        long_name="Initial V-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_V_3D"/>
+         <field id="vtau"         long_name="Surface Downward Y Stress"                               standard_name="surface_downward_y_stress"   unit="N/m2"                            />
+         <field id="voce"         long_name="Sea Water Y Velocity"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" />
+         <field id="voce_e3v"     long_name="ocean current along j-axis * e3v"                                                                   unit="m2/s"       grid_ref="grid_V_3D"  > voce * e3v </field>
+         <field id="ssv"          long_name="ocean surface current along j-axis"                                                                 unit="m/s"                             />
+         <field id="sbv"          long_name="ocean bottom current along j-axis"                                                                  unit="m/s"                             />
+         <field id="vbar"         long_name="ocean barotropic current along j-axis"                                                              unit="m/s"                             />
+         <field id="vocetr_eff"   long_name="Effective ocean transport along j-axis"                 standard_name="ocean_volume_y_transport"    unit="m3/s"       grid_ref="grid_V_3D" />
+         <field id="vocet"        long_name="ocean transport along j-axis times temperature (CRS)"                                               unit="degree_C*m/s"   grid_ref="grid_V_3D" />
+         <field id="voces"        long_name="ocean transport along j-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_V_3D" />
+
+         <!-- variables available with MLE -->
+         <field id="psiv_mle"     long_name="MLE streamfunction along j-axis"   unit="m3/s"   grid_ref="grid_V_3D" />
+
+         <!-- voce_eiv: available with key_traldf_eiv and key_diaeiv -->
+         <field id="voce_eiv"     long_name="EIV ocean current along j-axis"   standard_name="bolus_sea_water_y_velocity"   unit="m/s"   grid_ref="grid_V_3D" />
+         <field id="veiv_masstr"   long_name="EIV Ocean Mass Y Transport"    standard_name="bolus_ocean_mass_y_transport"                          unit="kg/s"        grid_ref="grid_V_3D" />
+         <field id="veiv_heattr3d" long_name="ocean bolus heat transport along j-axis"    standard_name="ocean_heat_y_transport_due_to_bolus_advection"   unit="W"    grid_ref="grid_V_3D" />
+         <field id="veiv_salttr3d" long_name="ocean bolus salt transport along j-axis"    standard_name="ocean_salt_y_transport_due_to_bolus_advection"   unit="kg"   grid_ref="grid_V_3D" />
+
+         <!-- voce_eiv: available with key_trabbl -->
+         <field id="voce_bbl"     long_name="BBL ocean current along j-axis"    unit="m/s"  />
+         <field id="ahv_bbl"      long_name="BBL diffusive flux along j-axis"   unit="m3/s" />
+
+         <!-- variable for ice shelves -->
+         <field id="vtbl"         long_name="meridional current in the Losh tbl"   unit="m/s" />
+
+         <!-- variables available with key_diaar5 -->
+         <field id="v_masstr"     long_name="Ocean Mass Y Transport"    standard_name="ocean_mass_y_transport"                          unit="kg/s"        grid_ref="grid_V_3D" />
+         <field id="v_heattr"     long_name="Ocean Heat X Transport"    standard_name="ocean_heat_y_transport"                          unit="W"                                />
+         <field id="vadv_heattr"  long_name="ocean advective heat transport along j-axis"    standard_name="advectice_ocean_heat_y_transport"                          unit="W"                                />
+         <field id="v_salttr"     long_name="ocean eulerian salt transport along i-axis"    standard_name="ocean_salt_y_transport"                          unit="0.001*kg/s"                        />
+         <field id="veiv_heattr"  long_name="ocean bolus heat transport along j-axis"       standard_name="ocean_heat_y_transport_due_to_bolus_advection"   unit="W"                                />
+         <field id="veiv_salttr"  long_name="ocean bolus salt transport along j-axis"       standard_name="ocean_salt_y_transport_due_to_bolus_advection"   unit="W"                                />
+         <field id="vdiff_heattr" long_name="ocean diffusion heat transport along j-axis"   standard_name="ocean_heat_y_transport_due_to_diffusion"         unit="W"                                />
+      </field_group>
+      
+      <!-- W grid -->
+      
+      <field_group id="grid_W" grid_ref="grid_W_3D">
+        <field id="e3w"          long_name="W-cell thickness"                     standard_name="cell_thickness"              unit="m"    />
+        <field id="woce"         long_name="ocean vertical velocity"              standard_name="upward_sea_water_velocity"   unit="m/s"  />
+        <field id="wocetr_eff"   long_name="effective ocean vertical transport"                                               unit="m3/s" />
+
+        <!-- woce_eiv: available with key_traldf_eiv and key_diaeiv -->
+        <field id="woce_eiv"     long_name="EIV ocean vertical velocity"   standard_name="bolus_upward_sea_water_velocity"   unit="m/s" />
+        <field id="weiv_masstr"  long_name="EIV Upward Ocean Mass Transport"  standard_name="bolus_upward_ocean_mass_transport"             unit="kg/s"   />
+        <field id="weiv_heattr3d" long_name="ocean bolus heat transport"    standard_name="ocean_heat_z_transport_due_to_bolus_advection"   unit="W"    grid_ref="grid_W_3D" />
+        <field id="weiv_salttr3d" long_name="ocean bolus salt transport"    standard_name="ocean_salt_z_transport_due_to_bolus_advection"   unit="kg"   grid_ref="grid_W_3D" />
+
+        <field id="avt"          long_name="vertical eddy diffusivity"   standard_name="ocean_vertical_heat_diffusivity"       unit="m2/s" />
+        <field id="logavt"       long_name="logarithm of vertical eddy diffusivity"   standard_name="ocean_vertical_heat_diffusivity"       unit="m2/s" />
+        <field id="avm"          long_name="vertical eddy viscosity"     standard_name="ocean_vertical_momentum_diffusivity"   unit="m2/s" />
+
+        <!-- avs: available with key_zdfddm -->
+        <field id="avs"          long_name="salt vertical eddy diffusivity"   standard_name="ocean_vertical_salt_diffusivity"   unit="m2/s" />
+        <field id="logavs"       long_name="logarithm of salt vertical eddy diffusivity"   standard_name="ocean_vertical_heat_diffusivity"       unit="m2/s" />
+
+        <!-- avt_evd and avm_evd: available with ln_zdfevd -->
+        <field id="avt_evd"      long_name="convective enhancement to vertical diffusivity"   standard_name="ocean_vertical_tracer_diffusivity_due_to_convection"     unit="m2/s" />
+        <field id="avm_evd"      long_name="convective enhancement to vertical viscosity"     standard_name="ocean_vertical_momentum_diffusivity_due_to_convection"   unit="m2/s" />
+
+        <!-- avt_tide: available with key_zdftmx -->
+        <field id="av_tide"      long_name="tidal vertical diffusivity"   standard_name="ocean_vertical_tracer_diffusivity_due_to_tides"   unit="m2/s" />
+
+       <!-- variables available with key_zdftmx_new -->
+        <field id="av_ratio"     long_name="S over T diffusivity ratio"            standard_name="salinity_over_temperature_diffusivity_ratio"                     unit="1"    />
+        <field id="av_wave"      long_name="wave-induced vertical diffusivity"     standard_name="ocean_vertical_tracer_diffusivity_due_to_internal_waves"         unit="m2/s" />
+        <field id="bflx_tmx"     long_name="wave-induced buoyancy flux"            standard_name="buoyancy_flux_due_to_internal_waves"                             unit="W/kg" />
+        <field id="pcmap_tmx"    long_name="power consumed by wave-driven mixing"  standard_name="vertically_integrated_power_consumption_by_wave_driven_mixing"   unit="W/m2"      grid_ref="grid_W_2D" />
+        <field id="emix_tmx"     long_name="power density available for mixing"    standard_name="power_available_for_mixing_from_breaking_internal_waves"         unit="W/kg" />
+
+        <!-- variables available with key_diaar5 -->   
+        <field id="w_masstr"     long_name="Upward Ocean Mass Transport"             standard_name="upward_ocean_mass_transport"             unit="kg/s"   />
+        <field id="w_masstr2"    long_name="square of vertical mass transport"   standard_name="square_of_upward_ocean_mass_transport"   unit="kg2/s2" />
+
+        <!-- aht2d and  aht2d_eiv: available with key_traldf_eiv and key_traldf_c2d -->
+        <field id="aht2d"        long_name="lateral eddy diffusivity"       standard_name="ocean_tracer_xy_laplacian_diffusivity"      unit="m2/s"   grid_ref="grid_W_2D" />
+        <field id="aht2d_eiv"    long_name="EIV lateral eddy diffusivity"   standard_name="ocean_tracer_bolus_laplacian_diffusivity"   unit="m2/s"   grid_ref="grid_W_2D" />
+
+         <!-- Variable from eosbn2 -->
+         <field id="bn2"          long_name="Squared buoyancy frequency measuring ocean vertical stratification"        standard_name="square_of_brunt_vaisala_frequency_in_sea_water"   unit="1/s^2"  />
+      </field_group>
+          
+      <!-- product fields -->
+      <field_group id="diaprod">
+	<field id="ut"           long_name="product_of_sea_water_x_velocity_and_potential_temperature"      unit="degree_C m/s"      grid_ref="grid_U_3D"   />
+        <field id="ut_e3u"       long_name="product_of_sea_water_x_velocity_and_potential_temperature * e3u"  unit="degree_C m2/s"   grid_ref="grid_U_3D" > ut * e3u </field >
+	<field id="us"           long_name="product_of_sea_water_x_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_U_3D"   />
+        <field id="us_e3u"       long_name="product_of_sea_water_x_velocity_and_salinity * e3u"             unit="PSU m2/s"      grid_ref="grid_U_3D" > us * e3u </field >
+	<field id="urhop"        long_name="product_of_sea_water_x_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_U_3D"   />
+        <field id="urhop_e3u"    long_name="product_of_sea_water_x_velocity_and_potential_density * e3u"    unit="(kg/m3).(m2/s)"   grid_ref="grid_U_3D" > urhop * e3u </field >
+	<field id="vt"           long_name="product_of_sea_water_y_velocity_and_potential_temperature"      unit="degree_C m/s"      grid_ref="grid_V_3D"   />
+        <field id="vt_e3v"       long_name="product_of_sea_water_y_velocity_and_potential_temperature * e3v"  unit="degree_C m2/s"   grid_ref="grid_V_3D" > vt * e3v </field >
+	<field id="vs"           long_name="product_of_sea_water_y_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_V_3D"   />
+        <field id="vs_e3v"       long_name="product_of_sea_water_y_velocity_and_salinity * e3t"             unit="PSU m2/s"      grid_ref="grid_V_3D" > vs * e3v </field >
+	<field id="vrhop"        long_name="product_of_sea_water_y_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_V_3D"   />
+        <field id="vrhop_e3v"    long_name="product_of_sea_water_y_velocity_and_potential_density * e3t"    unit="(kg/m3).(m2/s)"  grid_ref="grid_V_3D" > vrhop * e3v </field >
+	<field id="wt"           long_name="product_of_upward_sea_water_velocity_and_potential_temperature" unit="degree_C m/s"      grid_ref="grid_W_3D"   />
+	<field id="ws"           long_name="product_of_upward_sea_water_velocity_and_salinity"              unit="PSU m/s"       grid_ref="grid_W_3D"   />
+	<field id="wrhop"        long_name="product_of_upward_sea_water_velocity_and_potential_density"     unit="(kg/m3).(m/s)" grid_ref="grid_W_3D"   />
+      </field_group>
+
+      <!-- scalar variables available with key_diaar5 -->
+
+      <field_group id="scalar"  domain_ref="1point" >
+         <field id="voltot"     long_name="Sea Water Volume"                             standard_name="sea_water_volume"                               unit="m3"   />
+         <field id="sshtot"     long_name="global mean ssh"                              standard_name="global_average_sea_level_change"                unit="m"    />
+         <field id="sshsteric"  long_name="global mean ssh steric"                       standard_name="global_average_steric_sea_level_change"         unit="m"    />
+         <field id="sshthster"  long_name="Global Average Thermosteric Sea Level Change" standard_name="global_average_thermosteric_sea_level_change"   unit="m"    />
+         <field id="masstot"    long_name="global total mass"                            standard_name="sea_water_mass"                                 unit="kg"   />
+         <field id="temptot"    long_name="Global Average Sea Water Potential Temperature"      standard_name="sea_water_potential_temperature"         unit="degree_C" />
+         <field id="saltot"     long_name="Global Average Sea Water Salinity"                   standard_name="sea_water_salinity"                      unit="0.001" />
+         <field id="fram_trans" long_name="Sea Ice Mass Transport Through Fram Strait"   standard_name="sea_ice_transport_across_line"                  unit="kg/s" />
+
+      	 <!-- available with ln_diahsb -->
+    	 <field id="bgtemper"     long_name="drift in global mean temperature wrt timestep 1"                 standard_name="change_over_time_in_sea_water_potential_temperature"   unit="degC"     />
+    	 <field id="bgsaline"     long_name="drift in global mean salinity wrt timestep 1"                    standard_name="change_over_time_in_sea_water_practical_salinity"      unit="0.001"    />
+    	 <field id="bgheatco"     long_name="drift in global mean heat content wrt timestep 1"                                                                                      unit="1.e20J"   />
+     	 <field id="bgheatfx"     long_name="drift in global mean heat flux    wrt timestep 1"                                                                                      unit="W/m2"     />
+   	 <field id="bgsaltco"     long_name="drift in global mean salt content wrt timestep 1"                                                                                      unit="1e-3*km3" />
+    	 <field id="bgvolssh"     long_name="drift in global mean ssh volume wrt timestep 1"                                                                                        unit="km3"      />
+         <field id="bgvole3t"     long_name="drift in global mean volume variation (e3t) wrt timestep 1"                                                                            unit="km3"      />
+    	 <field id="bgfrcvol"     long_name="global mean volume from forcing"                                                                                                       unit="km3"      />
+    	 <field id="bgfrctem"     long_name="global mean heat content from forcing"                                                                                                 unit="1.e20J"   />
+         <field id="bgfrchfx"     long_name="global mean heat flux from forcing"                                                                                                    unit="W/m2"     />
+    	 <field id="bgfrcsal"     long_name="global mean salt content from forcing"                                                                                                 unit="1e-3*km3" />
+    	 <field id="bgmistem"     long_name="global mean temperature error due to free surface (no vvl)"                                                                            unit="degC"     />
+    	 <field id="bgmissal"     long_name="global mean salinity error due to free surface (no vvl)"                                                                               unit="1e-3"     />
+      </field_group>
+
+      <!-- LIM3 scalar variables -->
+
+      <field_group id="SBC_scalar"  domain_ref="1point" >
+         <!-- available with ln_limdiaout -->
+         <field id="ibgfrcvoltop"    long_name="global mean ice/snow forcing at interface ice/snow-atm (volume equivalent ocean volume)"   unit="km3"       />
+         <field id="ibgfrcvolbot"    long_name="global mean ice/snow forcing at interface ice/snow-ocean (volume equivalent ocean volume)" unit="km3"       />
+         <field id="ibgfrctemtop"    long_name="global mean heat on top of ice/snw/ocean-atm "                                             unit="1e20J"     />
+         <field id="ibgfrctembot"    long_name="global mean heat below ice (on top of ocean) "                                             unit="1e20J"     />
+         <field id="ibgfrcsal"       long_name="global mean ice/snow forcing (salt equivalent ocean volume)"                               unit="pss*km3"   />
+         <field id="ibgfrchfxtop"    long_name="global mean heat flux on top of ice/snw/ocean-atm "                                        unit="W/m2"      />
+         <field id="ibgfrchfxbot"    long_name="global mean heat flux below ice (on top of ocean) "                                        unit="W/m2"      />
+ 
+         <field id="ibgvolume"    long_name="drift in ice/snow volume (equivalent ocean volume)"            unit="km3"        />
+         <field id="ibgsaltco"    long_name="drift in ice salt content (equivalent ocean volume)"           unit="pss*km3"    />
+         <field id="ibgheatco"    long_name="drift in ice/snow heat content"                                unit="1e20J"      />
+         <field id="ibgheatfx"    long_name="drift in ice/snow heat flux"                                   unit="W/m2"       />
+
+         <field id="ibgvol_tot"    long_name="global mean ice volume"                                       unit="km3"        />
+         <field id="sbgvol_tot"    long_name="global mean snow volume"                                      unit="km3"        />
+         <field id="ibgarea_tot"   long_name="global mean ice area"                                         unit="km2"        />
+         <field id="ibgsalt_tot"   long_name="global mean ice salt content"                                 unit="1e-3*km3"   />
+         <field id="ibgheat_tot"   long_name="global mean ice heat content"                                 unit="1e20J"      />
+         <field id="sbgheat_tot"   long_name="global mean snow heat content"                                unit="1e20J"      />
+      </field_group>
+  
+      <!-- variables available with key_float -->
+
+      <field_group id="floatvar" domain_ref="1point" axis_ref="nfloat" operation="instant" >
+        <field id="traj_lon"      long_name="floats longitude"                                                           unit="degrees_east"  />
+        <field id="traj_lat"      long_name="floats latitude"                                                            unit="degrees_north" />
+        <field id="traj_dep"      long_name="floats depth"                                                               unit="m"             />
+        <field id="traj_temp"     long_name="floats temperature"       standard_name="sea_water_potential_temperature"   unit="degree_C"          />
+        <field id="traj_salt"     long_name="floats salinity"          standard_name="sea_water_practical_salinity"      unit="0.001"          />
+        <field id="traj_dens"     long_name="floats in-situ density"   standard_name="sea_water_density"                 unit="kg/m3"         />
+        <field id="traj_group"    long_name="floats group"                                                               unit="1"             />
+      </field_group>
+
+      <!-- variables available with iceberg trajectories -->
+
+      <field_group id="icbvar" domain_ref="grid_T"  > 
+        <field id="berg_melt"          long_name="icb melt rate of icebergs"                       unit="kg/m2/s"                    />
+        <field id="berg_buoy_melt"     long_name="icb buoyancy component of iceberg melt rate"     unit="kg/m2/s"                    />
+        <field id="berg_eros_melt"     long_name="icb erosion component of iceberg melt rate"      unit="kg/m2/s"                    />
+        <field id="berg_conv_melt"     long_name="icb convective component of iceberg melt rate"   unit="kg/m2/s"                    />
+        <field id="berg_virtual_area"  long_name="icb virtual coverage by icebergs"                unit="m2"                         />
+        <field id="bits_src"           long_name="icb mass source of bergy bits"                   unit="kg/m2/s"                    />
+        <field id="bits_melt"          long_name="icb melt rate of bergy bits"                     unit="kg/m2/s"                    />
+        <field id="bits_mass"          long_name="icb bergy bit density field"                     unit="kg/m2"                      />
+        <field id="berg_mass"          long_name="icb iceberg density field"                       unit="kg/m2"                      />
+        <field id="calving"            long_name="icb calving mass input"                          unit="kg/s"                       />
+        <field id="berg_floating_melt" long_name="icb melt rate of icebergs + bits"                unit="kg/m2/s"                    />
+        <field id="berg_real_calving"  long_name="icb calving into iceberg class"                  unit="kg/s"     axis_ref="icbcla" />
+        <field id="berg_stored_ice"    long_name="icb accumulated ice mass by class"               unit="kg"       axis_ref="icbcla" />
+      </field_group>
+
+      <!-- Poleward transport : ptr -->     
+      <field_group id="diaptr"> 
+        <field id="zomsfglo"          long_name="Meridional Stream-Function: Global"           unit="Sv"       grid_ref="gznl_W_3D" />
+        <field id="zomsfatl"          long_name="Meridional Stream-Function: Atlantic"         unit="Sv"       grid_ref="gznl_W_3D" />
+        <field id="zomsfpac"          long_name="Meridional Stream-Function: Pacific"          unit="Sv"       grid_ref="gznl_W_3D" />
+        <field id="zomsfind"          long_name="Meridional Stream-Function: Indian"           unit="Sv"       grid_ref="gznl_W_3D" />
+        <field id="zomsfipc"          long_name="Meridional Stream-Function: Pacific+Indian"   unit="Sv"       grid_ref="gznl_W_3D" />
+        <field id="zotemglo"          long_name="Zonal Mean Temperature : Global"              unit="degree_C"     grid_ref="gznl_T_3D" />
+        <field id="zotematl"          long_name="Zonal Mean Temperature : Atlantic"            unit="degree_C"     grid_ref="gznl_T_3D" />
+        <field id="zotempac"          long_name="Zonal Mean Temperature : Pacific"             unit="degree_C"     grid_ref="gznl_T_3D" />
+        <field id="zotemind"          long_name="Zonal Mean Temperature : Indian"              unit="degree_C"     grid_ref="gznl_T_3D" />
+        <field id="zotemipc"          long_name="Zonal Mean Temperature : Pacific+Indian"      unit="degree_C"     grid_ref="gznl_T_3D" />
+        <field id="zosalglo"          long_name="Zonal Mean Salinity : Global"                 unit="0.001"     grid_ref="gznl_T_3D" />
+        <field id="zosalatl"          long_name="Zonal Mean Salinity : Atlantic"               unit="0.001"     grid_ref="gznl_T_3D" />
+        <field id="zosalpac"          long_name="Zonal Mean Salinity : Pacific"                unit="0.001"     grid_ref="gznl_T_3D" />
+        <field id="zosalind"          long_name="Zonal Mean Salinity : Indian"                 unit="0.001"     grid_ref="gznl_T_3D" />
+        <field id="zosalipc"          long_name="Zonal Mean Salinity : Pacific+Indian"         unit="0.001"     grid_ref="gznl_T_3D" />
+        <field id="zosrfglo"          long_name="Zonal Mean Surface"                           unit="m2"       grid_ref="gznl_T_3D" />
+        <field id="zosrfatl"          long_name="Zonal Mean Surface : Atlantic"                unit="m2"       grid_ref="gznl_T_3D" />
+        <field id="zosrfpac"          long_name="Zonal Mean Surface : Pacific"                 unit="m2"       grid_ref="gznl_T_3D" />
+        <field id="zosrfind"          long_name="Zonal Mean Surface : Indian"                  unit="m2"       grid_ref="gznl_T_3D" />
+        <field id="zosrfipc"          long_name="Zonal Mean Surface : Pacific+Indian"          unit="m2"       grid_ref="gznl_T_3D" />
+        <field id="sophtadv"          long_name="Advective Heat Transport"                     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtadv_atl"      long_name="Advective Heat Transport: Atlantic"           unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtadv_pac"      long_name="Advective Heat Transport: Pacific"            unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtadv_ind"      long_name="Advective Heat Transport: Indian"             unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtadv_ipc"      long_name="Advective Heat Transport: Pacific+Indian"     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtldf"          long_name="Diffusive Heat Transport"                     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtldf_atl"      long_name="Diffusive Heat Transport: Atlantic"           unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtldf_pac"      long_name="Diffusive Heat Transport: Pacific"            unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtldf_ind"      long_name="Diffusive Heat Transport: Indian"             unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtldf_ipc"      long_name="Diffusive Heat Transport: Pacific+Indian"     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtove"          long_name="Overturning Heat Transport"                     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtove_atl"      long_name="Overturning Heat Transport: Atlantic"           unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtove_pac"      long_name="Overturning Heat Transport: Pacific"            unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtove_ind"      long_name="Overturning Heat Transport: Indian"             unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtove_ipc"      long_name="Overturning Heat Transport: Pacific+Indian"     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtbtr"          long_name="Barotropic Heat Transport"                     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtbtr_atl"      long_name="Barotropic Heat Transport: Atlantic"           unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtbtr_pac"      long_name="Barotropic Heat Transport: Pacific"            unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtbtr_ind"      long_name="Barotropic Heat Transport: Indian"             unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophtbtr_ipc"      long_name="Barotropic Heat Transport: Pacific+Indian"     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophteiv"          long_name="Heat Transport from mesoscale eddy advection"                     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophteiv_atl"      long_name="Heat Transport from mesoscale eddy advection: Atlantic"           unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophteiv_pac"      long_name="Heat Transport from mesoscale eddy advection: Pacific"            unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophteiv_ind"      long_name="Heat Transport from mesoscale eddy advection: Indian"             unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sophteiv_ipc"      long_name="Heat Transport from mesoscale eddy advection: Pacific+Indian"     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopstadv"          long_name="Advective Salt Transport"                     unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstadv_atl"      long_name="Advective Salt Transport: Atlantic"           unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstadv_pac"      long_name="Advective Salt Transport: Pacific"            unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstadv_ind"      long_name="Advective Salt Transport: Indian"             unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstadv_ipc"      long_name="Advective Salt Transport: Pacific+Indian"     unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstove"          long_name="Overturning Salt Transport"                     unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstove_atl"      long_name="Overturning Salt Transport: Atlantic"           unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstove_pac"      long_name="Overturning Salt Transport: Pacific"            unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstove_ind"      long_name="Overturning Salt Transport: Indian"             unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstove_ipc"      long_name="Overturning Salt Transport: Pacific+Indian"     unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstbtr"          long_name="Barotropic Salt Transport"                     unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstbtr_atl"      long_name="Barotropic Salt Transport: Atlantic"           unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstbtr_pac"      long_name="Barotropic Salt Transport: Pacific"            unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstbtr_ind"      long_name="Barotropic Salt Transport: Indian"             unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstbtr_ipc"      long_name="Barotropic Salt Transport: Pacific+Indian"     unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstldf"          long_name="Diffusive Salt Transport"                     unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstldf_atl"      long_name="Diffusive Salt Transport: Atlantic"           unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstldf_pac"      long_name="Diffusive Salt Transport: Pacific"            unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstldf_ind"      long_name="Diffusive Salt Transport: Indian"             unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopstldf_ipc"      long_name="Diffusive Salt Transport: Pacific+Indian"     unit="Giga g/s" grid_ref="gznl_T_2D" />
+        <field id="sopsteiv"          long_name="Salt Transport from mesoscale eddy advection"                     unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopsteiv_atl"      long_name="Salt Transport from mesoscale eddy advection: Atlantic"           unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopsteiv_pac"      long_name="Salt Transport from mesoscale eddy advection: Pacific"            unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopsteiv_ind"      long_name="Salt Transport from mesoscale eddy advection: Indian"             unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopsteiv_ipc"      long_name="Salt Transport from mesoscale eddy advection: Pacific+Indian"     unit="Giga g/s"       grid_ref="gznl_T_2D" />       
+        <field id="sopht_vt"          long_name="Heat Transport"                     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopht_vt_atl"      long_name="Heat Transport: Atlantic"           unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopht_vt_ind"      long_name="Heat Transport: Indian"             unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopht_vt_pac"      long_name="Heat Transport: Pacific"            unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopht_vt_ipc"      long_name="Heat Transport: Indo-Pacific"       unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vs"          long_name="Salt Transport"                     unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vs_atl"      long_name="Salt Transport: Atlantic"           unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vs_ind"      long_name="Salt Transport: Indian"             unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vs_pac"      long_name="Salt Transport: Pacific"            unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vs_ipc"      long_name="Salt Transport: Indo-Pacific"       unit="Giga g/s"       grid_ref="gznl_T_2D" />
+      </field_group>
+
+    <!-- 
+============================================================================================================
+                  Physical ocean model trend diagnostics : temperature, KE, PE, momentum
+============================================================================================================
+    -->
+
+    <field_group id="trendT" grid_ref="grid_T_3D">
+      <!-- variables available with ln_tra_trd -->
+      <field id="ttrd_xad"      long_name="temperature-trend: i-advection"                                                                                          unit="degC/s"                        />
+      <field id="strd_xad"      long_name="salinity   -trend: i-advection"                                                                                          unit="1e-3/s"                        />
+      <field id="ttrd_yad"      long_name="temperature-trend: j-advection"                                                                                          unit="degC/s"                        />
+      <field id="strd_yad"      long_name="salinity   -trend: j-advection"                                                                                          unit="1e-3/s"                        />
+      <field id="ttrd_zad"      long_name="temperature-trend: k-advection"                                                                                          unit="degC/s"                        />
+      <field id="strd_zad"      long_name="salinity   -trend: k-advection"                                                                                          unit="1e-3/s"                        />
+      <field id="ttrd_ad"       long_name="temperature-trend: advection"               standard_name="tendency_of_sea_water_temperature_due_to_advection"           unit="degC/s"                         > sqrt( ttrd_xad^2 + ttrd_yad^2 + ttrd_zad^2 ) </field>
+      <field id="strd_ad"       long_name="salinity   -trend: advection"               standard_name="tendency_of_sea_water_salinity_due_to_advection"              unit="1e-3/s"                         > sqrt( strd_xad^2 + strd_yad^2 + strd_zad^2 ) </field>
+      <field id="ttrd_totad"    long_name="temperature-trend: total advection"         standard_name="tendency_of_sea_water_salinity_due_to_advection"              unit="degC/s"                        />
+      <field id="strd_totad"    long_name="salinity   -trend: total advection"         standard_name="tendency_of_sea_water_salinity_due_to_advection"              unit="1e-3/s"                        />
+      <field id="ttrd_sad"      long_name="temperature-trend: surface adv. (no-vvl)"                                                                                unit="degC/s"   grid_ref="grid_T_2D" />
+      <field id="strd_sad"      long_name="salinity   -trend: surface adv. (no-vvl)"                                                                                unit="1e-3/s"   grid_ref="grid_T_2D" />
+      <field id="ttrd_ldf"      long_name="temperature-trend: lateral  diffusion"      standard_name="tendency_of_sea_water_temperature_due_to_horizontal_mixing"   unit="degC/s"                        />
+      <field id="strd_ldf"      long_name="salinity   -trend: lateral  diffusion"      standard_name="tendency_of_sea_water_salinity_due_to_horizontal_mixing"      unit="1e-3/s"                        />
+      <field id="ttrd_zdf"      long_name="temperature-trend: vertical diffusion"      standard_name="tendency_of_sea_water_temperature_due_to_vertical_mixing"     unit="degC/s"                        />
+      <field id="strd_zdf"      long_name="salinity   -trend: vertical diffusion"      standard_name="tendency_of_sea_water_salinity_due_to_vertical_mixing"        unit="1e-3/s"                        />
+      <field id="ttrd_evd"      long_name="temperature-trend: EVD convection"                                                                                       unit="degC/s"                        />
+      <field id="strd_evd"      long_name="salinity   -trend: EVD convection"                                                                                       unit="1e-3/s"                        />
+
+      <!-- ln_traldf_iso=T only (iso-neutral diffusion) -->
+      <field id="ttrd_iso_x"    long_name="temperature-trend: isopycnal diffusion x-component"                                                                      unit="degC/s"                        />
+      <field id="strd_iso_x"    long_name="salinity   -trend: isopycnal diffusion x-component"                                                                      unit="1e-3/s"                        />
+      <field id="ttrd_iso_y"    long_name="temperature-trend: isopycnal diffusion y-component"                                                                      unit="degC/s"                        />
+      <field id="strd_iso_y"    long_name="salinity   -trend: isopycnal diffusion y-component"                                                                      unit="1e-3/s"                        />
+      <field id="ttrd_iso_z1"   long_name="temperature-trend: isopycnal diffusion z-component part 1"                                                               unit="degC/s"                        />
+      <field id="strd_iso_z1"   long_name="salinity   -trend: isopycnal diffusion z-component part 1"                                                               unit="1e-3/s"                        />
+      <field id="ttrd_iso_z"    long_name="temperature-trend: isopycnal diffusion z-component"                 unit="degC/s" > ttrd_iso_z1 + ttrd_zdf - ttrd_zdfp </field>
+      <field id="strd_iso_z"    long_name="salinity   -trend: isopycnal diffusion z-component"                 unit="1e-3/s" > strd_iso_z1 + strd_zdf - strd_zdfp </field>
+      <field id="ttrd_iso"      long_name="temperature-trend: isopycnal diffusion"                             unit="degC/s" > ttrd_ldf + ttrd_zdf - ttrd_zdfp </field>
+      <field id="strd_iso"      long_name="salinity   -trend: isopycnal diffusion"                             unit="1e-3/s" > strd_ldf + strd_zdf - strd_zdfp </field>
+      <field id="ttrd_zdfp"     long_name="temperature-trend: pure vert. diffusion"   unit="degC/s" />
+      <field id="strd_zdfp"     long_name="salinity   -trend: pure vert. diffusion"   unit="1e-3/s" />
+
+      <!-- -->
+      <field id="ttrd_dmp"      long_name="temperature-trend: interior restoring"        unit="degree_C/s" />
+      <field id="strd_dmp"      long_name="salinity   -trend: interior restoring"        unit="0.001/s" />
+      <field id="ttrd_bbl"      long_name="temperature-trend: bottom boundary layer"     unit="degree_C/s" />
+      <field id="strd_bbl"      long_name="salinity   -trend: bottom boundary layer"     unit="0.001/s" />
+      <field id="ttrd_npc"      long_name="temperature-trend: non-penetrative conv."     unit="degree_C/s" />
+      <field id="strd_npc"      long_name="salinity   -trend: non-penetrative conv."     unit="0.001/s" />
+      <field id="ttrd_qns"      long_name="temperature-trend: non-solar flux + runoff"   unit="degC/s" grid_ref="grid_T_2D"/>
+      <field id="strd_cdt"      long_name="salinity   -trend: C/D term       + runoff"   unit="degC/s" grid_ref="grid_T_2D"/>
+      <field id="ttrd_qsr"      long_name="temperature-trend: solar penetr. heating"     unit="degree_C/s" />
+      <field id="ttrd_bbc"      long_name="temperature-trend: geothermal heating"        unit="degree_C/s" />
+      <field id="ttrd_atf"      long_name="temperature-trend: asselin time filter"       unit="degree_C/s" />
+      <field id="strd_atf"      long_name="salinity   -trend: asselin time filter"       unit="0.001/s" />
+      <field id="ttrd_tot"      long_name="temperature-trend: total model trend"         unit="degC/s" />
+      <field id="strd_tot"      long_name="salinity   -trend: total model trend"         unit="1e-3/s" />
+
+      <!-- Thickness weighted versions: -->
+      <field id="ttrd_xad_e3t"      unit="degC/s * m" >  ttrd_xad * e3t </field>
+      <field id="strd_xad_e3t"      unit="1e-3/s * m" >  strd_xad * e3t </field>
+      <field id="ttrd_yad_e3t"      unit="degC/s * m" >  ttrd_yad * e3t </field>
+      <field id="strd_yad_e3t"      unit="1e-3/s * m" >  strd_yad * e3t </field>
+      <field id="ttrd_zad_e3t"      unit="degC/s * m" >  ttrd_zad * e3t </field>
+      <field id="strd_zad_e3t"      unit="1e-3/s * m" >  strd_zad * e3t </field>
+      <field id="ttrd_ad_e3t"       unit="degC/s * m" >  ttrd_ad  * e3t </field>
+      <field id="strd_ad_e3t"       unit="1e-3/s * m" >  strd_ad  * e3t </field>
+      <field id="ttrd_totad_e3t"    unit="degC/s * m" >  ttrd_totad  * e3t </field>
+      <field id="strd_totad_e3t"    unit="1e-3/s * m" >  strd_totad  * e3t </field>
+      <field id="ttrd_ldf_e3t"      unit="degC/s * m" >  ttrd_ldf * e3t </field>
+      <field id="strd_ldf_e3t"      unit="1e-3/s * m" >  strd_ldf * e3t </field>
+      <field id="ttrd_zdf_e3t"      unit="degC/s * m" >  ttrd_zdf * e3t </field>
+      <field id="strd_zdf_e3t"      unit="1e-3/s * m" >  strd_zdf * e3t </field>
+      <field id="ttrd_evd_e3t"      unit="degC/s * m" >  ttrd_evd * e3t </field>
+      <field id="strd_evd_e3t"      unit="1e-3/s * m" >  strd_evd * e3t </field>
+
+      <!-- ln_traldf_iso=T only (iso-neutral diffusion) -->
+      <field id="ttrd_iso_x_e3t"    unit="degC/s * m"  >  ttrd_iso_x * e3t </field>
+      <field id="strd_iso_x_e3t"    unit="1e-3/s * m"  >  strd_iso_x * e3t </field>
+      <field id="ttrd_iso_y_e3t"    unit="degC/s * m"  >  ttrd_iso_y * e3t </field>
+      <field id="strd_iso_y_e3t"    unit="1e-3/s * m"  >  strd_iso_y * e3t </field>
+      <field id="ttrd_iso_z_e3t"    unit="degC/s * m"  >  ttrd_iso_z * e3t </field>
+      <field id="strd_iso_z_e3t"    unit="1e-3/s * m"  >  strd_iso_z * e3t </field>
+      <field id="ttrd_iso_e3t"      unit="degC/s * m"  >  ttrd_iso * e3t </field>
+      <field id="strd_iso_e3t"      unit="1e-3/s * m"  >  strd_iso * e3t </field>
+      <field id="ttrd_zdfp_e3t"     unit="degC/s * m"  >  ttrd_zdfp * e3t </field>
+      <field id="strd_zdfp_e3t"     unit="1e-3/s * m"  >  strd_zdfp * e3t </field>
+
+      <!-- -->
+      <field id="ttrd_dmp_e3t"      unit="degC/s * m"  >  ttrd_dmp * e3t </field>
+      <field id="strd_dmp_e3t"      unit="1e-3/s * m"  >  strd_dmp * e3t </field>
+      <field id="ttrd_bbl_e3t"      unit="degC/s * m"  >  ttrd_bbl * e3t </field>
+      <field id="strd_bbl_e3t"      unit="1e-3/s * m"  >  strd_bbl * e3t </field>
+      <field id="ttrd_npc_e3t"      unit="degC/s * m"  >  ttrd_npc * e3t </field>
+      <field id="strd_npc_e3t"      unit="1e-3/s * m"  >  strd_npc * e3t </field>
+      <field id="ttrd_qns_e3t"      unit="degC/s * m"  >  ttrd_qns * e3t_surf </field>
+      <field id="strd_cdt_e3t"      unit="degC/s * m"  >  strd_cdt * e3t_surf </field>
+      <field id="ttrd_qsr_e3t"      unit="degC/s * m"  >  ttrd_qsr * e3t </field>
+      <field id="ttrd_bbc_e3t"      unit="degC/s * m"  >  ttrd_bbc * e3t </field>
+      <field id="ttrd_atf_e3t"      unit="degC/s * m"  >  ttrd_atf * e3t </field>
+      <field id="strd_atf_e3t"      unit="1e-3/s * m"  >  strd_atf * e3t </field>
+      <field id="ttrd_tot_e3t"      unit="degC/s * m"  >  ttrd_tot * e3t </field>
+      <field id="strd_tot_e3t"      unit="1e-3/s * m"  >  strd_tot * e3t </field>
+
+      <!-- variables available with ln_KE_trd -->
+      <field id="ketrd_hpg"     long_name="ke-trend: hydrostatic pressure gradient"          unit="W/s^3"                        />
+      <field id="ketrd_spg"     long_name="ke-trend: surface     pressure gradient"          unit="W/s^3"                        />
+      <field id="ketrd_spgexp"  long_name="ke-trend: surface pressure gradient (explicit)"   unit="W/s^3"                        />
+      <field id="ketrd_spgflt"  long_name="ke-trend: surface pressure gradient (filter)"     unit="W/s^3"                        />
+      <field id="ssh_flt"       long_name="filtered contribution to ssh (dynspg_flt)"        unit="m"       grid_ref="grid_T_2D" />
+      <field id="w0"            long_name="surface vertical velocity"                        unit="m/s"     grid_ref="grid_T_2D" />
+      <field id="pw0_exp"       long_name="surface pressure flux due to ssh"                 unit="W/s^2"   grid_ref="grid_T_2D" />
+      <field id="pw0_flt"       long_name="surface pressure flux due to filtered ssh"        unit="W/s^2"   grid_ref="grid_T_2D" />
+      <field id="ketrd_keg"     long_name="ke-trend: KE gradient         or hor. adv."       unit="W/s^3"                        />
+      <field id="ketrd_rvo"     long_name="ke-trend: relative  vorticity or metric term"     unit="W/s^3"                        />
+      <field id="ketrd_pvo"     long_name="ke-trend: planetary vorticity"                    unit="W/s^3"                        />
+      <field id="ketrd_zad"     long_name="ke-trend: vertical  advection"                    unit="W/s^3"                        />
+      <field id="ketrd_udx"     long_name="ke-trend: U.dx[U]"                                unit="W/s^3"                        />
+      <field id="ketrd_ldf"     long_name="ke-trend: lateral   diffusion"                    unit="W/m^3"                        />
+      <field id="ketrd_zdf"     long_name="ke-trend: vertical  diffusion"                    unit="W/s^3"                        />
+      <field id="ketrd_tau"     long_name="ke-trend: wind stress "                           unit="W/s^3"   grid_ref="grid_T_2D" />
+      <field id="ketrd_bfr"     long_name="ke-trend: bottom friction (explicit)"             unit="W/s^3"                        />   
+      <field id="ketrd_bfri"    long_name="ke-trend: bottom friction (implicit)"             unit="W/s^3"                        />   
+      <field id="ketrd_atf"     long_name="ke-trend: asselin time filter trend"              unit="W/s^3"                        />  
+      <field id="ketrd_convP2K" long_name="ke-trend: conversion (potential to kinetic)"      unit="W/s^3"                        />
+      <field id="KE"            long_name="kinetic energy: u(n)*u(n+1)/2"                    unit="W/s^2"                        />   
+      <field id="ketrd_eiv"     long_name="Kinetic energy trend from mesoscale eddy advection"   unit="W/s^2" grid_ref="grid_T_2D" />   
+
+      <!-- variables available with ln_PE_trd -->
+      <field id="petrd_xad"     long_name="pe-trend: i-advection"                unit="W/m^3"                        />
+      <field id="petrd_yad"     long_name="pe-trend: j-advection"                unit="W/m^3"                        />
+      <field id="petrd_zad"     long_name="pe-trend: k-advection"                unit="W/m^3"                        />
+      <field id="petrd_sad"     long_name="pe-trend: surface adv. (no-vvl)"      unit="W/m^3"   grid_ref="grid_T_2D" />
+      <field id="petrd_ldf"     long_name="pe-trend: lateral  diffusion"         unit="W/m^3"                        />
+      <field id="petrd_zdf"     long_name="pe-trend: vertical diffusion"         unit="W/m^3"                        />
+      <field id="petrd_zdfp"    long_name="pe-trend: pure vert. diffusion"       unit="W/m^3"                        />
+      <field id="petrd_dmp"     long_name="pe-trend: interior restoring"         unit="W/m^3"                        />
+      <field id="petrd_bbl"     long_name="pe-trend: bottom boundary layer"      unit="W/m^3"                        />
+      <field id="petrd_npc"     long_name="pe-trend: non-penetrative conv."      unit="W/m^3"                        />
+      <field id="petrd_nsr"     long_name="pe-trend: surface forcing + runoff"   unit="W/m^3"                        />
+      <field id="petrd_qsr"     long_name="pe-trend: solar penetr. heating"      unit="W/m^3"                        />
+      <field id="petrd_bbc"     long_name="pe-trend: geothermal heating"         unit="W/m^3"                        />
+      <field id="petrd_atf"     long_name="pe-trend: asselin time filter"        unit="W/m^3"                        />
+      <field id="PEanom"        long_name="potential energy anomaly"             unit="1"                            />   
+      <field id="alphaPE"       long_name="partial deriv. of PEanom wrt T"       unit="degree_C-1"                       />   
+      <field id="betaPE"        long_name="partial deriv. of PEanom wrt S"       unit="1e3"                          />   
+      <field id="petrd_zdfp_e3t"    unit="W/m^2" grid_ref="grid_T_3D"  > petrd_zdfp * e3t </field>
+    </field_group>
+
+    <field_group id="trendU" grid_ref="grid_U_3D">
+     <!-- variables available with ln_dyn_trd -->
+     <field id="utrd_hpg"       long_name="i-trend: hydrostatic pressure gradient"          unit="m/s^2"                        />
+     <field id="utrd_spg"       long_name="i-trend: surface     pressure gradient"          unit="m/s^2"                        />
+     <field id="utrd_spgexp"    long_name="i-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        />
+     <field id="utrd_spgflt"    long_name="i-trend: surface pressure gradient (filtered)"   unit="m/s^2"                        />
+     <field id="utrd_keg"       long_name="i-trend: KE gradient         or hor. adv."       unit="m/s^2"                        />
+     <field id="utrd_rvo"       long_name="i-trend: relative  vorticity or metric term"     unit="m/s^2"                        />
+     <field id="utrd_pvo"       long_name="i-trend: planetary vorticity"                    unit="m/s^2"                        />
+     <field id="utrd_zad"       long_name="i-trend: vertical  advection"                    unit="m/s^2"                        />
+     <field id="utrd_udx"       long_name="i-trend: U.dx[U]"                                unit="m/s^2"                        />
+     <field id="utrd_ldf"       long_name="i-trend: lateral   diffusion"                    unit="m/s^2"                        />
+     <field id="utrd_zdf"       long_name="i-trend: vertical  diffusion"                    unit="m/s^2"                        />
+     <field id="utrd_tau"       long_name="i-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_U_2D" />
+     <field id="utrd_bfr"       long_name="i-trend: bottom friction (explicit)"             unit="m/s^2"                        />   
+     <field id="utrd_bfri"      long_name="i-trend: bottom friction (implicit)"             unit="m/s^2"                        />   
+     <field id="utrd_tot"       long_name="i-trend: total momentum trend before atf"        unit="m/s^2"                        />   
+     <field id="utrd_atf"       long_name="i-trend: asselin time filter trend"              unit="m/s^2"                        />   
+    </field_group>
+
+    <field_group id="trendV" grid_ref="grid_V_3D">
+     <!-- variables available with ln_dyn_trd -->
+     <field id="vtrd_hpg"       long_name="j-trend: hydrostatic pressure gradient"          unit="m/s^2"                        />
+     <field id="vtrd_spg"       long_name="j-trend: surface     pressure gradient"          unit="m/s^2"                        />
+     <field id="vtrd_spgexp"    long_name="j-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        />
+     <field id="vtrd_spgflt"    long_name="j-trend: surface pressure gradient (filtered)"   unit="m/s^2"                        />
+     <field id="vtrd_keg"       long_name="j-trend: KE gradient         or hor. adv."       unit="m/s^2"                        />
+     <field id="vtrd_rvo"       long_name="j-trend: relative  vorticity or metric term"     unit="m/s^2"                        />
+     <field id="vtrd_pvo"       long_name="j-trend: planetary vorticity"                    unit="m/s^2"                        />
+     <field id="vtrd_zad"       long_name="j-trend: vertical  advection"                    unit="m/s^2"                        />
+     <field id="vtrd_vdy"       long_name="i-trend: V.dx[V]"                                unit="m/s^2"                        />
+     <field id="vtrd_ldf"       long_name="j-trend: lateral   diffusion"                    unit="m/s^2"                        />
+     <field id="vtrd_zdf"       long_name="j-trend: vertical  diffusion"                    unit="m/s^2"                        />
+     <field id="vtrd_tau"       long_name="j-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_V_2D" />
+     <field id="vtrd_bfr"       long_name="j-trend: bottom friction (explicit)"             unit="m/s^2"                        />   
+     <field id="vtrd_bfri"      long_name="j-trend: bottom friction (implicit)"             unit="m/s^2"                        />   
+     <field id="vtrd_tot"       long_name="j-trend: total momentum trend before atf"        unit="m/s^2"                        />   
+     <field id="vtrd_atf"       long_name="j-trend: asselin time filter trend"              unit="m/s^2"                        />   
+    </field_group>
+
+    <!-- 
+============================================================================================================
+                                     Biogeochemistry model variables
+============================================================================================================
+    -->
+
+      <!-- ptrc on T grid -->
+
+     <field_group id="ptrc_T" grid_ref="grid_T_3D">
+       <field id="DIC"          long_name="Dissolved inorganic Concentration"        unit="mmol/m3" />
+       <field id="DIC_E3T"      long_name="DIC * E3T"                                unit="mmol/m2" > DIC * e3t </field >
+       <field id="Alkalini"     long_name="Total Alkalinity Concentration"           unit="mmol/m3" />
+       <field id="Alkalini_E3T" long_name="Alkalini * E3T"                           unit="mmol/m2"  > Alkalini * e3t </field >
+       <field id="O2"           long_name="Oxygen Concentration"                     unit="mmol/m3" />
+       <field id="O2_E3T"       long_name="O2 * E3T"                                 unit="mmol/m2"  > O2 * e3t </field >
+       <field id="CaCO3"        long_name="Calcite Concentration"                    unit="mmol/m3" />
+       <field id="CaCO3_E3T"    long_name="CaCO3 * E3T"                              unit="mmol/m2"  > CaCO3 * e3t </field >
+       <field id="PO4"          long_name="Phosphate Concentration"                  unit="mmol/m3" />
+       <field id="PO4_E3T"      long_name="PO4 * E3T"                                unit="mmol/m2"  > PO4 * e3t </field >
+       <field id="POC"          long_name="Small organic carbon Concentration"       unit="mmol/m3" />
+       <field id="POC_E3T"      long_name="POC * E3T"                                unit="mmol/m2"  > POC * e3t </field >
+       <field id="Si"           long_name="Silicate Concentration"                   unit="mmol/m3" />
+       <field id="Si_E3T"       long_name="Si * E3T"                                 unit="mmol/m2"  > Si * e3t </field >
+       <field id="PHY"          long_name="(Nano)Phytoplankton Concentration"        unit="mmol/m3" />
+       <field id="PHY_E3T"      long_name="PHY * E3T"                                unit="mmol/m2"  > PHY * e3t </field >
+       <field id="ZOO"          long_name="(Micro)Zooplankton Concentration"         unit="mmol/m3" />
+       <field id="ZOO_E3T"      long_name="ZOO2 * E3T"                               unit="mmol/m2"  > ZOO * e3t </field >
+       <field id="DOC"          long_name="Dissolved organic Concentration"          unit="mmol/m3" />
+       <field id="DOC_E3T"      long_name="DOC * E3T"                                unit="mmol/m2"  > DOC * e3t </field >
+       <field id="PHY2"         long_name="Diatoms Concentration"                    unit="mmol/m3" />
+       <field id="PHY2_E3T"     long_name="PHY2 * E3T"                               unit="mmol/m2"  > PHY2 * e3t </field >
+       <field id="ZOO2"         long_name="Mesozooplankton Concentration"            unit="mmol/m3" />
+       <field id="ZOO2_E3T"     long_name="ZOO2 * E3T"                               unit="mmol/m2"  > ZOO2 * e3t </field >
+       <field id="DSi"          long_name="Diatoms Silicate Concentration"           unit="mmol/m3" />
+       <field id="DSi_E3T"      long_name="Dsi * E3T"                                unit="mmol/m2"  > DSi * e3t </field >
+       <field id="Fer"          long_name="Dissolved Iron Concentration"             unit="mmol/m3" />
+       <field id="Fer_E3T"      long_name="Fer * E3T"                                unit="mmol/m2"  > Fer * e3t </field >
+       <field id="BFe"          long_name="Big iron particles Concentration"         unit="mmol/m3" />
+       <field id="BFe_E3T"      long_name="BFe * E3T"                                unit="mmol/m2"  > BFe * e3t </field >
+       <field id="GOC"          long_name="Big organic carbon Concentration"         unit="mmol/m3" />
+       <field id="GOC_E3T"      long_name="GOC * E3T"                                unit="mmol/m2"  > GOC * e3t </field >
+       <field id="SFe"          long_name="Small iron particles Concentration"       unit="mmol/m3" />
+       <field id="SFe_E3T"      long_name="SFe * E3T"                                unit="mmol/m2"  > SFe * e3t </field >
+       <field id="DFe"          long_name="Diatoms iron  Concentration"              unit="mmol/m3" />
+       <field id="DFe_E3T"      long_name="DFe * E3T"                                unit="mmol/m2"  > DFe * e3t </field >
+       <field id="GSi"          long_name="Sinking biogenic Silicate Concentration"  unit="mmol/m3" />
+       <field id="GSi_E3T"      long_name="GSi * E3T"                                unit="mmol/m2"  > GSi * e3t </field >
+       <field id="NFe"          long_name="Nano iron Concentration"                  unit="mmol/m3" />
+       <field id="NFe_E3T"      long_name="NFe * E3T"                                unit="mmol/m2"  > NFe * e3t </field >
+       <field id="NCHL"         long_name="Nano chlorophyl Concentration"            unit="mg/m3"   />
+       <field id="NCHL_E3T"     long_name="NCHL * E3T"                               unit="mmol/m2"  > NCHL * e3t </field >
+       <field id="DCHL"         long_name="Diatoms chlorophyl Concentration"         unit="mg/m3"   />
+       <field id="DCHL_E3T"     long_name="DCHL * E3T"                               unit="mmol/m2"  > DCHL * e3t </field >
+       <field id="NO3"          long_name="Nitrate Concentration"                    unit="mmol/m3" />
+       <field id="NO3_E3T"      long_name="NO3 * E3T"                                unit="mmol/m2"  > NO3 * e3t </field >
+       <field id="NH4"          long_name="Ammonium Concentration"                   unit="mmol/m3" />
+       <field id="NH4_E3T"      long_name="NH4 * E3T"                                unit="mmol/m2"  > NH4 * e3t </field >
+
+       <!-- PISCES with Kriest parametisation : variables available with key_kriest -->
+       <field id="Num"         long_name="Number of organic particles"              unit="1" />
+       <field id="Num_E3T"     long_name="Num * E3T"                                unit="m"  > Num * e3t </field >
+
+       <!-- PISCES light : variables available with key_pisces_reduced -->
+       <field id="DET"         long_name="Detritus"                                 unit="mmol-N/m3" />
+       <field id="DET_E3T"     long_name="DET * E3T"                                unit="mmol-N/m2"  > DET * e3t </field >
+       <field id="DOM"         long_name="Dissolved Organic Matter"                 unit="mmol-N/m3" />
+       <field id="DOM_E3T"     long_name="DOM * E3T"                                unit="mmol-N/m2"  > DOM * e3t </field >
+
+       <!-- CFC11 : variables available with key_cfc -->
+       <field id="CFC11"       long_name="CFC-11 Concentration"                     unit="umol/m3" />
+       <field id="CFC11_E3T"   long_name="CFC11 * E3T"                              unit="umol/m2"  > CFC11 * e3t </field >
+
+       <!-- Bomb C14 : variables available with key_c14b -->
+       <field id="C14B"       long_name="Bomb C14 Concentration"                    unit="1"         />
+       <field id="C14B_E3T"   long_name="C14B * E3T"                                unit="m"  > C14B * e3t </field >
+
+       <!-- AGE : variables available with key_age -->
+       <field id="Age"        long_name="Sea Water Age Since Surface contact" standard_name="sea_water_age_since_surface_contact"  unit="yr"         />
+       <field id="Age_e3t"    long_name="Age * e3t"                                 unit="yr * m"  > Age * e3t </field >
+
+     </field_group>
+
+     <!-- PISCES additional diagnostics on T grid  -->
+
+     <field_group id="diad_T" grid_ref="grid_T_2D">
+       <field id="PH"          long_name="PH"                                      unit="1"          grid_ref="grid_T_3D" />
+       <field id="CO3"         long_name="Bicarbonates"                            unit="mol/m3"     grid_ref="grid_T_3D" />
+       <field id="CO3sat"      long_name="CO3 saturation"                          unit="mol/m3"     grid_ref="grid_T_3D" />
+       <field id="PAR"         long_name="Photosynthetically Available Radiation"  unit="W/m2"       grid_ref="grid_T_3D" />
+       <field id="PARDM"       long_name="Daily mean PAR"                          unit="W/m2"       grid_ref="grid_T_3D" />
+       <field id="PPPHY"       long_name="Primary production of nanophyto"         unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="PPPHY2"      long_name="Primary production of diatoms"           unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="PPNEWN"      long_name="New Primary production of nanophyto"     unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="PPNEWD"      long_name="New Primary production of diatoms"       unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="PBSi"        long_name="Primary production of Si diatoms"        unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="PFeN"        long_name="Primary production of nano iron"         unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="PFeD"        long_name="Primary production of diatoms iron"      unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="xfracal"     long_name="Calcifying fraction"                     unit="1"          grid_ref="grid_T_3D" />
+       <field id="PCAL"        long_name="Calcite production"                      unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="DCAL"        long_name="Calcite dissolution"                     unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="GRAZ1"       long_name="Grazing by microzooplankton"             unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="GRAZ2"       long_name="Grazing by mesozooplankton"              unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="REMIN"       long_name="Oxic remineralization of OM"             unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="DENIT"       long_name="Anoxic remineralization of OM"           unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="Nfix"        long_name="Nitrogen fixation"                       unit="mol/m3/s"   grid_ref="grid_T_3D" />
+       <field id="Mumax"       long_name="Maximum growth rate"                     unit="s-1"        grid_ref="grid_T_3D" />
+       <field id="MuN"         long_name="Realized growth rate for nanophyto"      unit="s-1"        grid_ref="grid_T_3D" />
+       <field id="MuD"         long_name="Realized growth rate for diatomes"       unit="s-1"        grid_ref="grid_T_3D" />
+       <field id="LNnut"       long_name="Nutrient limitation term in Nanophyto"   unit=""           grid_ref="grid_T_3D" />
+       <field id="LDnut"       long_name="Nutrient limitation term in Diatoms"     unit=""           grid_ref="grid_T_3D" />
+       <field id="LNFe"        long_name="Iron limitation term in Nanophyto"       unit=""           grid_ref="grid_T_3D" />
+       <field id="LDFe"        long_name="Iron limitation term in Diatoms"         unit=""           grid_ref="grid_T_3D" />
+       <field id="LNlight"     long_name="Light limitation term in Nanophyto"      unit=""           grid_ref="grid_T_3D" />
+       <field id="LDlight"     long_name="Light limitation term in Diatoms"        unit=""           grid_ref="grid_T_3D" />
+       <field id="Fe2"         long_name="Iron II concentration"                   unit="nmol/m3"    grid_ref="grid_T_3D" />
+       <field id="Fe3"         long_name="Iron III concentration"                  unit="nmol/m3"    grid_ref="grid_T_3D" />
+       <field id="FeL1"        long_name="Complexed Iron concentration with L1"    unit="nmol/m3"    grid_ref="grid_T_3D" />
+       <field id="FeL2"        long_name="Complexed Iron concentration with L2"    unit="nmol/m3"    grid_ref="grid_T_3D" />
+       <field id="FeP"         long_name="Precipitated Iron III"                   unit="nmol/m3"    grid_ref="grid_T_3D" />
+       <field id="TL1"         long_name="Total L1 concentration"                  unit="nmol/m3"    grid_ref="grid_T_3D" />
+       <field id="TL2"         long_name="Total L2 concentration"                  unit="nmol/m3"    grid_ref="grid_T_3D" />
+       <field id="pdust"       long_name="dust concentration"                      unit="g/m3"                            />
+       <field id="Totlig"      long_name="Total ligand concentation"               unit="nmol/m3"    grid_ref="grid_T_3D" />
+       <field id="Biron"       long_name="Bioavailable iron"                       unit="nmol/m3"    grid_ref="grid_T_3D" />
+       <field id="Sdenit"      long_name="Nitrate reduction in the sediments"      unit="molN/m2/s"                       />
+       <field id="SedCal"      long_name="Calcite burial in the sediments"         unit="molC/m2/s"                       />
+       <field id="SedSi"       long_name="Silicon burial in the sediments"         unit="molSi/m2/s"                      />
+       <field id="SedC"        long_name="Organic C burial in the sediments"       unit="molC/m2/s"                       />
+       <field id="Ironice"     long_name="Iron input/uptake due to sea ice"        unit="mol/m2/s"                        />
+       <field id="HYDR"        long_name="Iron input from hydrothemal vents"       unit="mol/m2/s"   grid_ref="grid_T_3D" />
+       <field id="EPC100"      long_name="Export of carbon particles at 100 m"     unit="mol/m2/s"                        />
+       <field id="EPFE100"     long_name="Export of biogenic iron at 100 m"        unit="mol/m2/s"                        />
+       <field id="EPSI100"     long_name="Export of Silicate at 100 m"             unit="mol/m2/s"                        />
+       <field id="EPCAL100"    long_name="Export of Calcite at 100 m"              unit="mol/m2/s"                        />
+       <field id="EXPC"        long_name="Export of carbon"                        unit="mol/m2/s"   grid_ref="grid_T_3D" />
+       <field id="EXPFE"       long_name="Export of biogenic iron"                 unit="mol/m2/s"   grid_ref="grid_T_3D" />
+       <field id="EXPSI"       long_name="Export of Silicate"                      unit="mol/m2/s"   grid_ref="grid_T_3D" />
+       <field id="EXPCAL"      long_name="Export of Calcite"                       unit="mol/m2/s"   grid_ref="grid_T_3D" />
+       <field id="Cflx"        long_name="DIC flux"                                unit="mol/m2/s"                        />
+       <field id="Oflx"        long_name="Oxygen flux"                             unit="mol/m2/s"                        />
+       <field id="Kg"          long_name="Gas transfer"                            unit="mol/m2/s/uatm"                   />
+       <field id="Dpco2"       long_name="Delta CO2"                               unit="uatm"                            />
+       <field id="Dpo2"        long_name="Delta O2"                                unit="uatm"                            />
+       <field id="Heup"        long_name="Euphotic layer depth"                    unit="m"                               />
+       <field id="Irondep"     long_name="Iron deposition from dust"               unit="mol/m2/s"                        />
+       <field id="Ironsed"     long_name="Iron deposition from sediment"           unit="mol/m2/s"   grid_ref="grid_T_3D" />
+
+
+       <!-- PISCES with Kriest parametisation : variables available with key_kriest -->
+       <field id="EPN100"      long_name="Particulate number flux at 100 m"        unit="mol/m2/s"                        />
+       <field id="EXPN"        long_name="Particulate number flux"                 unit="mol/m2/s"   grid_ref="grid_T_3D" />
+       <field id="XNUM"        long_name="Number of particles in aggregats"        unit="1"          grid_ref="grid_T_3D" />
+       <field id="WSC"         long_name="sinking speed of mass flux"              unit="m2/s"       grid_ref="grid_T_3D" />
+       <field id="WSN"         long_name="sinking speed of number flux"            unit="m2/s"       grid_ref="grid_T_3D" />
+
+       <!-- dbio_T on T grid : variables available with key_diaar5 -->
+       <field id="TPP"         long_name="Total Primary production of phyto"                   unit="mol/m3/s"  grid_ref="grid_T_3D" />
+       <field id="TPNEW"       long_name="New Primary production of phyto"                     unit="mol/m3/s"  grid_ref="grid_T_3D" />
+       <field id="TPBFE"       long_name="Total biogenic iron production"                      unit="mol/m3/s"  grid_ref="grid_T_3D" />
+       <field id="INTDIC"      long_name="DIC content"                                         unit="kg/m2"                          />
+       <field id="O2MIN"       long_name="Oxygen minimum concentration"                        unit="mol/m3"                          />
+       <field id="ZO2MIN"      long_name="Depth of oxygen minimum concentration"               unit="m"                              />
+       <field id="INTNFIX"     long_name="Nitrogen fixation rate : vert. integrated"           unit="mol/m2/s"                       />
+       <field id="INTPPPHY"    long_name="Vertically integrated primary production by nanophy" unit="mol/m2/s"                       />
+       <field id="INTPPPHY2"   long_name="Vertically integrated primary production by diatom"  unit="mol/m2/s"                       />
+       <field id="INTPP"       long_name="Vertically integrated primary production by phyto"   unit="mol/m2/s"                       />
+       <field id="INTPNEW"     long_name="Vertically integrated new primary production"        unit="mol/m2/s"                       />
+       <field id="INTPBFE"     long_name="Vertically integrated of biogenic iron production"   unit="mol/m2/s"                       />
+       <field id="INTPBSI"     long_name="Vertically integrated of biogenic Si production"     unit="mol/m2/s"                       />
+       <field id="INTPCAL"     long_name="Vertically integrated of calcite production"         unit="mol/m2/s"                       />
+
+       <!-- PISCES light : variables available with key_pisces_reduced -->
+       <field id="FNO3PHY"     long_name="FNO3PHY"                                 unit=""          grid_ref="grid_T_3D" /> 
+       <field id="FNH4PHY"     long_name="FNH4PHY"                                 unit=""          grid_ref="grid_T_3D" /> 
+       <field id="FNH4NO3"     long_name="FNH4NO3"                                 unit=""          grid_ref="grid_T_3D" /> 
+       <field id="TNO3PHY"     long_name="TNO3PHY"                                 unit=""  /> 
+       <field id="TNH4PHY"     long_name="TNH4PHY"                                 unit=""  /> 
+       <field id="TPHYDOM"     long_name="TPHYDOM"                                 unit=""  /> 
+       <field id="TPHYNH4"     long_name="TPHYNH4"                                 unit=""  /> 
+       <field id="TPHYZOO"     long_name="TPHYZOO"                                 unit=""  /> 
+       <field id="TPHYDET"     long_name="TPHYDET"                                 unit=""  /> 
+       <field id="TDETZOO"     long_name="TDETZOO"                                 unit=""  /> 
+       <field id="TZOODET"     long_name="TZOODET"                                 unit=""  /> 
+       <field id="TZOOBOD"     long_name="TZOOBOD"                                 unit=""  /> 
+       <field id="TZOONH4"     long_name="TZOONH4"                                 unit=""  /> 
+       <field id="TZOODOM"     long_name="TZOODOM"                                 unit=""  /> 
+       <field id="TNH4NO3"     long_name="TNH4NO3"                                 unit=""  /> 
+       <field id="TDOMNH4"     long_name="TDOMNH4"                                 unit=""  /> 
+       <field id="TDETNH4"     long_name="TDETNH4"                                 unit=""  /> 
+       <field id="TPHYTOT"     long_name="TPHYTOT"                                 unit=""  /> 
+       <field id="TZOOTOT"     long_name="TZOOTOT"                                 unit=""  /> 
+       <field id="SEDPOC"      long_name="SEDPOC"                                  unit=""  /> 
+       <field id="TDETSED"     long_name="TDETSED"                                 unit=""  /> 
+
+       <!-- CFC11 : variables available with key_cfc -->
+       <field id="qtrCFC11"    long_name="Air-sea flux of CFC-11"                  unit="mol/m2/s"   />
+       <field id="qintCFC11"   long_name="Cumulative air-sea flux of CFC-11"       unit="mol/m2"     />
+
+       <!-- Bomb C14 : variables available with key_c14b -->
+       <field id="qtrC14b"     long_name="Air-sea flux of Bomb C14"                unit="mol/m2/s"   />
+       <field id="qintC14b"    long_name="Cumulative air-sea flux of Bomb C14"     unit="mol/m2"     />
+       <field id="fdecay"      long_name="Radiactive decay of Bomb C14"            unit="mol/m3"  grid_ref="grid_T_3D"  />
+     </field_group>
+
+     <field_group id="PISCES_scalar"  domain_ref="1point" >
+       <field id="pno3tot"         long_name="global mean nitrate concentration"                  unit="mol/m3"   />
+       <field id="ppo4tot"         long_name="global mean phosphorus concentration"               unit="mol/m3"   />
+       <field id="psiltot"         long_name="global mean silicate concentration"                 unit="mol/m3"   />
+       <field id="palktot"         long_name="global mean alkalinity concentration"               unit="mol/m3"   />
+       <field id="pfertot"         long_name="global mean iron concentration"                     unit="mol/m3"   />
+       <field id="tcflx"           long_name="total Flux of Carbon out of the ocean"              unit="mol/s"   />
+       <field id="tcflxcum"        long_name="cumulative total Flux of Carbon out of the ocean"   unit="mol/s"   />
+       <field id="tcexp"           long_name="total Carbon export at 100m"                        unit="mol/s"   />
+       <field id="tintpp"          long_name="global total integrated primary production"         unit="mol/s"   />
+       <field id="tnfix"           long_name="global total nitrogen fixation"                     unit="mol/s"   />
+       <field id="tdenit"          long_name="Total denitrification"                              unit="mol/s"   />
+     </field_group>
+
+    <!-- 
+============================================================================================================
+                                        Definitions for iodef_demo.xml
+============================================================================================================
+    -->
+
+    <field_group id="mooring" >
+      <field field_ref="toce"         name="thetao"   long_name="sea_water_potential_temperature"      />
+      <field field_ref="soce"         name="so"       long_name="sea_water_salinity"                   />
+      <field field_ref="uoce"         name="uo"       long_name="sea_water_x_velocity"                 />
+      <field field_ref="voce"         name="vo"       long_name="sea_water_y_velocity"                 />
+      <field field_ref="woce"         name="wo"       long_name="sea_water_z_velocity"                 />
+      <field field_ref="avt"          name="difvho"   long_name="ocean_vertical_heat_diffusivity"      />
+      <field field_ref="avm"          name="difvmo"   long_name="ocean_vertical_momentum_diffusivity"  />
+   
+      <field field_ref="sst"          name="tos"      long_name="sea_surface_temperature"                       />
+      <field field_ref="sst2"         name="tossq"    long_name="square_of_sea_surface_temperature"             />
+      <field field_ref="sstgrad"      name="tosgrad"  long_name="module_of_sea_surface_temperature_gradient"    />
+      <field field_ref="sss"          name="sos"      long_name="sea_surface_salinity"                          />
+      <field field_ref="ssh"          name="zos"      long_name="sea_surface_height_above_geoid"                />
+      <field field_ref="empmr"        name="wfo"      long_name="water_flux_into_sea_water"                     />
+      <field field_ref="qsr"          name="rsntds"   long_name="surface_net_downward_shortwave_flux"           />
+      <field field_ref="qt"           name="tohfls"   long_name="surface_net_downward_total_heat_flux"          />
+      <field field_ref="taum"                                                                                   />
+      <field field_ref="20d"                                                                                    />
+      <field field_ref="mldkz5"                                                                                 />
+      <field field_ref="mldr10_1"                                                                               />
+      <field field_ref="mldr10_3"                                                                               />
+      <field field_ref="mldr0_1"                                                                                />
+      <field field_ref="mldr0_3"                                                                                />
+      <field field_ref="mld_dt02"                                                                               />
+      <field field_ref="topthdep"                                                                               />
+      <field field_ref="pycndep"                                                                                />
+      <field field_ref="tinv"                                                                                   />
+      <field field_ref="depti"                                                                                  />
+      <field field_ref="BLT"          name="blt"      long_name="barrier_layer_thickness"                       />
+      <field field_ref="utau"         name="tauuo"    long_name="surface_downward_x_stress"                     />
+      <field field_ref="vtau"         name="tauvo"    long_name="surface_downward_y_stress"                     />
+    </field_group>
+
+    <field_group id="groupT" >
+      <field field_ref="toce"         name="thetao"   long_name="sea_water_potential_temperature"               />
+      <field field_ref="soce"         name="so"       long_name="sea_water_salinity"                            />
+      <field field_ref="sst"          name="tos"      long_name="sea_surface_temperature"                       />
+      <field field_ref="sst2"         name="tossq"    long_name="square_of_sea_surface_temperature"             />
+      <field field_ref="sss"          name="sos"      long_name="sea_surface_salinity"                          />
+      <field field_ref="ssh"          name="zos"      long_name="sea_surface_height_above_geoid"                />
+      <field field_ref="empmr"        name="wfo"      long_name="water_flux_into_sea_water"                     />
+      <field field_ref="qsr"          name="rsntds"   long_name="surface_net_downward_shortwave_flux"           />
+      <field field_ref="qt"           name="tohfls"   long_name="surface_net_downward_total_heat_flux"          />
+      <field field_ref="taum"                                                                                   />
+      <field field_ref="20d"                                                                                    />
+      <field field_ref="mldkz5"                                                                                 />
+      <field field_ref="mldr10_1"                                                                               />
+      <field field_ref="mldr10_3"                                                                               />
+      <field field_ref="mld_dt02"                                                                               />
+      <field field_ref="topthdep"                                                                               />
+      <field field_ref="pycndep"                                                                                />
+      <field field_ref="tinv"                                                                                   />
+      <field field_ref="depti"                                                                                  />
+      <field field_ref="BLT"          name="blt"      long_name="Barrier Layer Thickness"                       />
+    </field_group>
+    
+    <field_group id="groupU" >
+      <field field_ref="uoce"         name="uo"      long_name="sea_water_x_velocity"      />
+      <field field_ref="ssu"          name="uos"     long_name="sea_surface_x_velocity"    />
+      <field field_ref="utau"         name="tauuo"   long_name="surface_downward_x_stress" />
+    </field_group>
+    
+    <field_group id="groupV" >
+      <field field_ref="voce"         name="vo"      long_name="sea_water_y_velocity"      />
+      <field field_ref="ssv"          name="vos"     long_name="sea_surface_y_velocity"    />
+      <field field_ref="vtau"         name="tauvo"   long_name="surface_downward_y_stress" />
+    </field_group>
+    
+    <field_group id="groupW" >
+      <field field_ref="woce"         name="wo"       long_name="ocean vertical velocity"  />
+    </field_group>
+    
+
+   </field_definition>
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/stpctl.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/stpctl.F90	(revision 3293)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/stpctl.F90	(revision 3293)
@@ -0,0 +1,223 @@
+MODULE stpctl
+   !!======================================================================
+   !!                       ***  MODULE  stpctl  ***
+   !! Ocean run control :  gross check of the ocean time stepping
+   !!======================================================================
+   !! History :  OPA  ! 1991-03  (G. Madec) Original code
+   !!            6.0  ! 1992-06  (M. Imbard)
+   !!            8.0  ! 1997-06  (A.M. Treguier)
+   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
+   !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp_ctl      : Control the run
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE sol_oce         ! ocean space and time domain variables 
+   USE sbc_oce         ! surface boundary conditions variables
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! distributed memory computing
+   USE lib_fortran     ! Fortran routines library 
+   USE dynspg_oce      ! pressure gradient schemes 
+   USE c1d             ! 1D vertical configuration
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC stp_ctl           ! routine called by step.F90
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: stpctl.F90 7847 2017-03-30 13:30:29Z cetlod $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE stp_ctl( kt, kindic )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE stp_ctl  ***
+      !!                     
+      !! ** Purpose :   Control the run
+      !!
+      !! ** Method  : - Save the time step in numstp
+      !!              - Print it each 50 time steps
+      !!              - Print solver statistics in numsol 
+      !!              - Stop the run IF problem for the solver ( indec < 0 )
+      !!
+      !! ** Actions :   'time.step' file containing the last ocean time-step
+      !!                
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
+      INTEGER, INTENT( inout ) ::   kindic  ! indicator of solver convergence
+      !!
+      CHARACTER(len = 32) ::        clfname ! time stepping output file name
+      INTEGER  ::   ji, jj, jk              ! dummy loop indices
+      INTEGER  ::   ii, ij, ik              ! temporary integers
+      REAL(wp) ::   zumax, zsmin, zssh2, zsshmax     ! temporary scalars
+      INTEGER, DIMENSION(3) ::   ilocu      ! 
+      INTEGER, DIMENSION(2) ::   ilocs      ! 
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'stp_ctl : time-stepping control'
+         WRITE(numout,*) '~~~~~~~'
+         ! open time.step file with special treatment for SAS
+         IF ( nn_components == jp_iam_sas ) THEN
+            clfname = 'time.step.sas'
+         ELSE
+            clfname = 'time.step'
+         ENDIF
+         CALL ctl_opn( numstp, TRIM(clfname), 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+      ENDIF
+
+      IF(lwp) WRITE ( numstp, '(1x, i8)' )   kt      !* save the current time step in numstp
+      IF(lwp) REWIND( numstp )                       !  --------------------------
+
+      !                                              !* Test maximum of velocity (zonal only)
+      !                                              !  ------------------------
+      !! zumax = MAXVAL( ABS( un(:,:,:) ) )                ! slower than the following loop on NEC SX5
+      zumax = 0.e0
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zumax = MAX(zumax,ABS(un(ji,jj,jk)))
+          END DO 
+        END DO 
+      END DO        
+      IF( lk_mpp )   CALL mpp_max( zumax )                 ! max over the global domain
+      !
+      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' abs(U) max: ', zumax
+      !
+      IF( zumax > 20.e0 ) THEN
+         IF( lk_mpp ) THEN
+            CALL mpp_maxloc(ABS(un),umask,zumax,ii,ij,ik)
+         ELSE
+            ilocu = MAXLOC( ABS( un(:,:,:) ) )
+            ii = ilocu(1) + nimpp - 1
+            ij = ilocu(2) + njmpp - 1
+            ik = ilocu(3)
+         ENDIF
+         IF(lwp) THEN
+            WRITE(numout,cform_err)
+            WRITE(numout,*) ' stpctl: the zonal velocity is larger than 20 m/s'
+            WRITE(numout,*) ' ====== '
+            WRITE(numout,9400) kt, zumax, ii, ij, ik
+            WRITE(numout,*)
+            WRITE(numout,*) '          output of last fields in numwso'
+         ENDIF
+         kindic = -3
+      ENDIF
+9400  FORMAT (' kt=',i6,' max abs(U): ',1pg11.4,', i j k: ',3i5)
+
+      !                                              !* Test minimum of salinity
+      !                                              !  ------------------------
+      !! zsmin = MINVAL( tsn(:,:,1,jp_sal), mask = tmask(:,:,1) == 1.e0 )  slower than the following loop on NEC SX5
+      zsmin = 100.e0
+      DO jj = 2, jpjm1
+         DO ji = 1, jpi
+            IF( tmask(ji,jj,1) == 1) zsmin = MIN(zsmin,tsn(ji,jj,1,jp_sal))
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_min( zsmin )                ! min over the global domain
+      !
+      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' SSS min:', zsmin
+      !
+      IF( zsmin < 0.) THEN 
+         IF (lk_mpp) THEN
+            CALL mpp_minloc ( tsn(:,:,1,jp_sal),tmask(:,:,1), zsmin, ii,ij )
+         ELSE
+            ilocs = MINLOC( tsn(:,:,1,jp_sal), mask = tmask(:,:,1) == 1.e0 )
+            ii = ilocs(1) + nimpp - 1
+            ij = ilocs(2) + njmpp - 1
+         ENDIF
+         !
+         IF(lwp) THEN
+            WRITE(numout,cform_err)
+            WRITE(numout,*) 'stp_ctl : NEGATIVE sea surface salinity'
+            WRITE(numout,*) '======= '
+            WRITE(numout,9500) kt, zsmin, ii, ij
+            WRITE(numout,*)
+            WRITE(numout,*) '          output of last fields in numwso'
+         ENDIF
+         kindic = -3
+      ENDIF
+9500  FORMAT (' kt=',i6,' min SSS: ',1pg11.4,', i j: ',2i5)
+
+      
+      IF( lk_c1d )  RETURN          ! No log file in case of 1D vertical configuration
+
+      ! log file (solver or ssh statistics)
+      ! --------
+      IF( lk_dynspg_flt ) THEN      ! elliptic solver statistics (if required)
+         !
+         IF(lwp) WRITE(numsol,9200) kt, niter, res, SQRT(epsr)/eps       ! Solver
+         !
+         IF( kindic < 0 .AND. zsmin > 0.e0 .AND. zumax <= 20.e0 ) THEN   ! create a abort file if problem found 
+            IF(lwp) THEN
+               WRITE(numout,*) ' stpctl: the elliptic solver DO not converge or explode'
+               WRITE(numout,*) ' ====== '
+               WRITE(numout,9200) kt, niter, res, sqrt(epsr)/eps
+               WRITE(numout,*)
+               WRITE(numout,*) ' stpctl: output of last fields'
+               WRITE(numout,*) ' ======  '
+            ENDIF
+         ENDIF
+         !
+      ELSE                                   !* ssh statistics (and others...)
+         IF( kt == nit000 .AND. lwp ) THEN   ! open ssh statistics file (put in solver.stat file)
+            CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+         ENDIF
+         !
+         zsshmax = 0.e0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( tmask(ji,jj,1) == 1) zsshmax = MAX( zsshmax, ABS(sshn(ji,jj)) )
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_max( zsshmax )                ! min over the global domain
+         !
+         IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' ssh max:', zsshmax
+         !
+         IF( zsshmax > 10.e0 ) THEN 
+            IF (lk_mpp) THEN
+              CALL mpp_maxloc( ABS(sshn(:,:)),tmask(:,:,1),zsshmax,ii,ij)
+            ELSE
+              ilocs = MAXLOC( ABS(sshn(:,:)) )
+              ii = ilocs(1) + nimpp - 1
+              ij = ilocs(2) + njmpp - 1
+            ENDIF
+            !
+            IF(lwp) THEN
+               WRITE(numout,cform_err)
+               WRITE(numout,*) 'stp_ctl : the ssh is larger than 10m'
+               WRITE(numout,*) '======= '
+               WRITE(numout,9600) kt, zsshmax, ii, ij
+               WRITE(numout,*)
+               WRITE(numout,*) '          output of last fields in numwso'
+            ENDIF
+           kindic = -3
+       ENDIF
+9600   FORMAT (' kt=',i6,' max ssh: ',1pg11.4,', i j: ',2i5)
+        !
+        zssh2 = glob_sum( sshn(:,:) * sshn(:,:) )
+        !
+        IF(lwp) WRITE(numsol,9300) kt, zssh2, zumax, zsmin      ! ssh statistics
+        !
+        IF( isnan( zssh2 ) .OR. isnan( zsmin ) .OR. isnan( zumax ) ) THEN           
+            kindic = -3                                       
+            IF(lwp) WRITE(numout,*) " ssh too big. output.abort "
+        ENDIF   
+        !
+      ENDIF
+9200  FORMAT('it:', i8, ' iter:', i4, ' r: ',d23.16, ' b: ',d23.16 )
+9300  FORMAT(' it :', i8, ' ssh2: ', d23.16, ' Umax: ',d23.16,' Smin: ',d23.16)
+      !
+   END SUBROUTINE stp_ctl
+
+   !!======================================================================
+END MODULE stpctl
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/zdftke.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/zdftke.F90	(revision 3293)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/zdftke.F90	(revision 3293)
@@ -0,0 +1,970 @@
+MODULE zdftke
+   !!======================================================================
+   !!                       ***  MODULE  zdftke  ***
+   !! Ocean physics:  vertical mixing coefficient computed from the tke 
+   !!                 turbulent closure parameterization
+   !!=====================================================================
+   !! History :  OPA  !  1991-03  (b. blanke)  Original code
+   !!            7.0  !  1991-11  (G. Madec)   bug fix
+   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
+   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
+   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
+   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
+   !!            8.0  !  1997-07  (G. Madec)   lbc
+   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
+   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
+   !!             -   !  2004-10  (C. Ethe )  1D configuration
+   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
+   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
+   !!                 !           - tke penetration (wind steering)
+   !!                 !           - suface condition for tke & mixing length
+   !!                 !           - Langmuir cells
+   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
+   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
+   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term 
+   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl 
+   !!                 !                                + cleaning of the parameters + bugs correction
+   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
+   !!----------------------------------------------------------------------
+#if defined key_zdftke   ||   defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   'key_zdftke'                                   TKE vertical physics
+   !!----------------------------------------------------------------------
+   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
+   !!   tke_tke       : tke time stepping: update tke at now time step (en)
+   !!   tke_avn       : compute mixing length scale and deduce avm and avt
+   !!   zdf_tke_init  : initialization, namelist read, and parameters control
+   !!   tke_rst       : read/write tke restart in ocean restart file
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean: dynamics and active tracers variables
+   USE phycst         ! physical constants
+   USE dom_oce        ! domain: ocean
+   USE domvvl         ! domain: variable volume layer
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE zdf_oce        ! vertical physics: ocean variables
+   USE zdfmxl         ! vertical physics: mixed layer
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+#if defined key_agrif
+   USE agrif_opa_interp
+   USE agrif_opa_update
+#endif
+!CE
+#if defined key_lim3
+    USE ice, ONLY: htm_i, ht_i
+#endif
+#if defined key_lim2 || defined key_cice
+    USE sbc_ice, ONLY: ht_i
+#endif
+!CE
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_tke        ! routine called in step module
+   PUBLIC   zdf_tke_init   ! routine called in opa module
+   PUBLIC   tke_rst        ! routine called in step module
+
+   LOGICAL , PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
+
+   !                      !!** Namelist  namzdf_tke  **
+   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
+!CE
+   INTEGER  ::   nn_mxl0   ! type of scaling under sea-ice (=0/1/2/3)
+   REAL(wp) ::   rn_hice   ! ice thickness value when scaling under sea-ice
+!CE
+   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
+   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
+   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
+   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
+   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation 
+   REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
+   REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
+   REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
+   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
+   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
+   INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1)
+   REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean
+   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
+   REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells
+
+   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
+   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
+   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
+   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
+
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
+#if defined key_c1d
+   !                                                                        !!** 1D cfg only  **   ('key_c1d')
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_dis, e_mix   !: dissipation and mixing turbulent lengh scales
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers
+#endif
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION zdf_tke_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_tke_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE(                                                                    &
+#if defined key_c1d
+         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
+         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
+#endif
+         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc      )
+         !
+      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
+      IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
+      !
+   END FUNCTION zdf_tke_alloc
+
+
+   SUBROUTINE zdf_tke( kt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_tke  ***
+      !!
+      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
+      !!              coefficients using a turbulent closure scheme (TKE).
+      !!
+      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
+      !!              is computed from a prognostic equation :
+      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
+      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
+      !!                  + avt N^2                      ! stratif. destruc.
+      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
+      !!      with the boundary conditions:
+      !!         surface: en = max( rn_emin0, rn_ebb * taum )
+      !!         bottom : en = rn_emin
+      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 
+      !!
+      !!        The now Turbulent kinetic energy is computed using the following 
+      !!      time stepping: implicit for vertical diffusion term, linearized semi
+      !!      implicit for kolmogoroff dissipation term, and explicit forward for 
+      !!      both buoyancy and shear production terms. Therefore a tridiagonal 
+      !!      linear system is solved. Note that buoyancy and shear terms are
+      !!      discretized in a energy conserving form (Bruchard 2002).
+      !!
+      !!        The dissipative and mixing length scale are computed from en and
+      !!      the stratification (see tke_avn)
+      !!
+      !!        The now vertical eddy vicosity and diffusivity coefficients are
+      !!      given by: 
+      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
+      !!              avt = max( avmb, pdl * avm                 )  
+      !!              eav = max( avmb, avm )
+      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
+      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1 
+      !!
+      !! ** Action  :   compute en (now turbulent kinetic energy)
+      !!                update avt, avmu, avmv (before vertical eddy coef.)
+      !!
+      !! References : Gaspar et al., JGR, 1990,
+      !!              Blanke and Delecluse, JPO, 1991
+      !!              Mellor and Blumberg, JPO 2004
+      !!              Axell, JGR, 2002
+      !!              Bruchard OM 2002
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!----------------------------------------------------------------------
+      !
+      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
+         avt (:,:,:) = avt_k (:,:,:) 
+         avm (:,:,:) = avm_k (:,:,:) 
+         avmu(:,:,:) = avmu_k(:,:,:) 
+         avmv(:,:,:) = avmv_k(:,:,:) 
+      ENDIF 
+      !
+      CALL tke_tke      ! now tke (en)
+      !
+      CALL tke_avn      ! now avt, avm, avmu, avmv
+      !
+      avt_k (:,:,:) = avt (:,:,:) 
+      avm_k (:,:,:) = avm (:,:,:) 
+      avmu_k(:,:,:) = avmu(:,:,:) 
+      avmv_k(:,:,:) = avmv(:,:,:) 
+      !
+#if defined key_agrif
+      ! Update child grid f => parent grid 
+      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
+#endif      
+     ! 
+   END SUBROUTINE zdf_tke
+
+
+   SUBROUTINE tke_tke
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tke_tke  ***
+      !!
+      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
+      !!
+      !! ** Method  : - TKE surface boundary condition
+      !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
+      !!              - source term due to shear (saved in avmu, avmv arrays)
+      !!              - Now TKE : resolution of the TKE equation by inverting
+      !!                a tridiagonal linear system by a "methode de chasse"
+      !!              - increase TKE due to surface and internal wave breaking
+      !!
+      !! ** Action  : - en : now turbulent kinetic energy)
+      !!              - avmu, avmv : production of TKE by shear at u and v-points
+      !!                (= Kz dz[Ub] * dz[Un] )
+      !! ---------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
+!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar
+!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar
+      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
+      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
+      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
+      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
+      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
+      REAL(wp) ::   ztau  , zdif                    !    -         -
+      REAL(wp) ::   zus   , zwlc  , zind            !    -         -
+      REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
+!!bfr      REAL(wp) ::   zebot                           !    -         -
+      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw
+      !!--------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tke_tke')
+      !
+      CALL wrk_alloc( jpi,jpj, imlc )    ! integer
+      CALL wrk_alloc( jpi,jpj, zhlc ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      !
+      zbbrau = rn_ebb / rau0       ! Local constant initialisation
+      zfact1 = -.5_wp * rdt 
+      zfact2 = 1.5_wp * rdt * rn_ediss
+      zfact3 = 0.5_wp       * rn_ediss
+      !
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Surface boundary condition on tke
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF ( ln_isfcav ) THEN
+         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1)
+            END DO
+         END DO
+      END IF
+      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+         END DO
+      END DO
+      
+!!bfr   - start commented area
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Bottom boundary condition on tke
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      ! Tests to date have found the bottom boundary condition on tke to have very little effect.
+      ! The condition is coded here for completion but commented out until there is proof that the
+      ! computational cost is justified
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
+!CDIR NOVERRCHK
+!!    DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+!!       DO ji = fs_2, fs_jpim1   ! vector opt.
+!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
+!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
+!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + &
+!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
+!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
+!!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
+!!       END DO
+!!    END DO
+!!bfr   - end commented area
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002)
+         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+         !
+         !                        !* total energy produce by LC : cumulative sum over jk
+         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)
+         DO jk = 2, jpk
+            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
+         END DO
+         !                        !* finite Langmuir Circulation depth
+         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
+         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
+         DO jk = jpkm1, 2, -1
+            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
+               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
+                  zus  = zcof * taum(ji,jj)
+                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
+               END DO
+            END DO
+         END DO
+         !                               ! finite LC depth
+         DO jj = 1, jpj 
+            DO ji = 1, jpi
+               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
+            END DO
+         END DO
+         zcof = 0.016 / SQRT( zrhoa * zcdrag )
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
+                  !                                           ! vertical velocity due to LC
+                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
+                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
+                  !                                           ! TKE Langmuir circulation source term
+                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   &
+                     &   zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Now Turbulent kinetic energy (output in en)
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
+      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
+      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
+      !
+      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
+         DO jj = 1, jpj                 ! here avmu, avmv used as workspace
+            DO ji = 1, jpi
+               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
+                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
+                  &                            / (  fse3uw_n(ji,jj,jk)               &
+                  &                              *  fse3uw_b(ji,jj,jk)  )
+               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
+                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
+                  &                            / (  fse3vw_n(ji,jj,jk)               &
+                  &                              *  fse3vw_b(ji,jj,jk)  )
+            END DO
+         END DO
+      END DO
+      !
+      DO jk = 2, jpkm1           !* Matrix and right hand side in en
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zcof   = zfact1 * tmask(ji,jj,jk)
+# if defined key_zdftmx_new
+               ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability)
+               zzd_up = zcof * ( MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp ) )   &  ! upper diagonal
+                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
+               zzd_lw = zcof * ( MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp ) )   &  ! lower diagonal
+                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
+# else
+               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
+                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
+               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
+                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
+# endif
+                  !                                                           ! shear prod. at w-point weightened by mask
+               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+                  !
+               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
+               zd_lw(ji,jj,jk) = zzd_lw
+               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
+               !
+               !                                   ! right hand side in en
+               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
+                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
+                  &                                 * wmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
+      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
+         END DO
+      END DO
+      DO jk = 3, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 2, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      DO jk = 2, jpkm1                             ! set the minimum value of tke
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  TKE due to surface and internal wave breaking
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                     &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+               END DO
+            END DO
+         END DO
+      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               jk = nmln(ji,jj)
+               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                  &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+            END DO
+         END DO
+      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
+                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
+                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress 
+                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
+                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
+                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                     &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
+      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
+      CALL wrk_dealloc( jpi,jpj, zhlc ) 
+      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
+      !
+   END SUBROUTINE tke_tke
+
+
+   SUBROUTINE tke_avn
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tke_avn  ***
+      !!
+      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
+      !!
+      !! ** Method  :   At this stage, en, the now TKE, is known (computed in 
+      !!              the tke_tke routine). First, the now mixing lenth is 
+      !!      computed from en and the strafification (N^2), then the mixings
+      !!      coefficients are computed.
+      !!              - Mixing length : a first evaluation of the mixing lengh
+      !!      scales is:
+      !!                      mxl = sqrt(2*en) / N  
+      !!      where N is the brunt-vaisala frequency, with a minimum value set
+      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
+      !!        The mixing and dissipative length scale are bound as follow : 
+      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
+      !!                        zmxld = zmxlm = mxl
+      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
+      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 
+      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
+      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
+      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to 
+      !!                    the surface to obtain ldown. the resulting length 
+      !!                    scales are:
+      !!                        zmxld = sqrt( lup * ldown ) 
+      !!                        zmxlm = min ( lup , ldown )
+      !!              - Vertical eddy viscosity and diffusivity:
+      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
+      !!                      avt = max( avmb, pdlr * avm )  
+      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
+      !!
+      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
+      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
+      REAL(wp) ::   zdku, zpdlr, zri, zsqen     !   -      -
+      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
+      !!--------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tke_avn')
+
+      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Mixing length
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
+      !
+      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
+      zmxlm(:,:,:)  = rmxl_min    
+      zmxld(:,:,:)  = rmxl_min
+      !
+      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
+        !
+         zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
+         !
+         SELECT CASE( nn_mxl0 )             ! Type of scaling under sea-ice
+         !
+         CASE( 0 )                      ! No scaling under sea-ice
+            zmxlm(:,:,1) = zraug * taum(:,:) * tmask(:,:,1) 
+         !
+         CASE( 1 )                           ! scaling with constant sea-ice thickness
+            zmxlm(:,:,1) =  rn_hice 
+         !
+         CASE( 2 )                                 ! scaling with mean sea-ice thickness
+#if defined key_lim3
+            zmxlm(:,:,1) = ( 1. - fr_i(:,:) ) * zraug * taum(:,:) * tmask(:,:,1) + fr_i(:,:) * htm_i(:,:) * 2.
+#elif ( defined key_lim2 || defined key_cice )
+            zmxlm(:,:,1) = ( 1. - fr_i(:,:) ) * zraug * taum(:,:) * tmask(:,:,1) + fr_i(:,:) * MAXVAL( ht_i, 3 )
+#else
+            zmxlm(:,:,1) = zraug * taum(:,:) * tmask(:,:,1) 
+#endif
+            !
+         CASE( 3 )                                 ! scaling with max sea-ice thickness
+#if defined key_lim3 || defined key_lim2 || defined key_cice 
+            zmxlm(:,:,1) = ( 1. - fr_i(:,:) ) * zraug * taum(:,:) * tmask(:,:,1) + fr_i(:,:) * MAXVAL( ht_i, 3 )
+#else
+            zmxlm(:,:,1) = zraug * taum(:,:) * tmask(:,:,1) 
+#endif
+            !
+         END SELECT
+         ! 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
+            END DO
+         END DO
+         !
+      ELSE 
+         zmxlm(:,:,1) = rn_mxl0
+      ENDIF
+      !
+      IF( iom_use('mixlength') )   THEN   
+          CALL lbc_lnk( zmxlm(:,:,1) , 'T', 1. )
+          CALL iom_put( "mixlength" , zmxlm(:,:,1) ) 
+      ENDIF
+      !
+!CDIR NOVERRCHK
+      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
+!CDIR NOVERRCHK
+         DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
+               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
+            END DO
+         END DO
+      END DO
+      !
+      !                     !* Physical limits for the mixing length
+      !
+      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value 
+      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
+      !
+      SELECT CASE ( nn_mxl )
+      !
+      ! where wmask = 0 set zmxlm == fse3w
+      CASE ( 0 )           ! bounded by the distance to surface and bottom
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
+                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
+                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
+                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
+                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
+               END DO
+            END DO
+         END DO
+         !
+      CASE ( 1 )           ! bounded by the vertical scale factor
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = zemxl
+                  zmxld(ji,jj,jk) = zemxl
+               END DO
+            END DO
+         END DO
+         !
+      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
+         DO jk = 2, jpkm1         ! from the surface to the bottom :
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = zemxl
+                  zmxld(ji,jj,jk) = zemxl
+               END DO
+            END DO
+         END DO
+         !
+      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
+         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
+                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = zemlm
+                  zmxld(ji,jj,jk) = zemlp
+               END DO
+            END DO
+         END DO
+         !
+      END SELECT
+      !
+# if defined key_c1d
+      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
+      e_mix(:,:,:) = zmxlm(:,:,:)
+# endif
+
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!CDIR NOVERRCHK
+      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
+!CDIR NOVERRCHK
+         DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zsqen = SQRT( en(ji,jj,jk) )
+               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
+               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
+               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
+               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
+      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
+               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
+      !
+      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
+                  !                                          ! shear
+                  zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) )   &
+                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
+                  zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) )   &
+                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
+                  !                                          ! local Richardson number
+                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
+                  zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  )
+!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
+!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  )
+                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
+# if defined key_c1d
+                  e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number
+                  e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri
+# endif
+              END DO
+            END DO
+         END DO
+      ENDIF
+      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
+         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
+            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
+      !
+   END SUBROUTINE tke_avn
+
+
+   SUBROUTINE zdf_tke_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_tke_init  ***
+      !!                     
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!              viscosity when using a tke turbulent closure scheme
+      !!
+      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
+      !!              called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namlist namzdf_tke
+      !!
+      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      INTEGER ::   ios, ierr
+      !!
+      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
+         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
+!CE
+         &                 nn_mxl0 , rn_hice,  &
+!CE
+         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
+         &                 nn_etau , nn_htau  , rn_efr   
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
+      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
+      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namzdf_tke )
+      !
+      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
+# if defined key_zdftmx_new
+      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
+      rn_emin  = 1.e-10_wp
+      rmxl_min = 1.e-03_wp
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+# else
+      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
+# endif
+      !
+      IF(lwp) THEN                    !* Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
+         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
+         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
+         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
+         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
+         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
+         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
+         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
+         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
+         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
+!CE
+         IF( ln_mxl0 ) THEN
+            WRITE(numout,*) '      type of scaling under sea-ice                  nn_mxl0   = ', nn_mxl0
+            IF( nn_mxl0 == 1 ) &
+            WRITE(numout,*) '      ice thickness value when scaling under sea-ice rn_hice   = ', rn_hice
+         ENDIF
+!CE
+         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
+         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
+         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
+         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
+         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
+         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
+         WRITE(numout,*)
+         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
+      ENDIF
+      !
+      !                              ! allocate tke arrays
+      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
+      !
+      !                               !* Check of some namelist values
+      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
+      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
+      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
+      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
+
+      IF( ln_mxl0 ) THEN
+         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
+         rn_mxl0 = rmxl_min
+      ENDIF
+      
+      IF( nn_etau == 2  ) THEN
+          ierr = zdf_mxl_alloc()
+          nmln(:,:) = nlb10           ! Initialization of nmln
+      ENDIF
+
+      !                               !* depth of penetration of surface tke
+      IF( nn_etau /= 0 ) THEN      
+         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
+         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
+            htau(:,:) = 10._wp
+         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
+            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )            
+         END SELECT
+      ENDIF
+      !                               !* set vertical eddy coef. to the background value
+      DO jk = 1, jpk
+         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
+         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
+         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
+         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
+      END DO
+      dissl(:,:,:) = 1.e-12_wp
+      !                              
+      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
+      !
+   END SUBROUTINE zdf_tke_init
+
+
+   SUBROUTINE tke_rst( kt, cdrw )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE tke_rst  ***
+     !!                     
+     !! ** Purpose :   Read or write TKE file (en) in restart file
+     !!
+     !! ** Method  :   use of IOM library
+     !!                if the restart does not contain TKE, en is either 
+     !!                set to rn_emin or recomputed 
+     !!----------------------------------------------------------------------
+     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+     !
+     INTEGER ::   jit, jk   ! dummy loop indices
+     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
+     !!----------------------------------------------------------------------
+     !
+     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+        !                                   ! ---------------
+        IF( ln_rstart ) THEN                   !* Read the restart file
+           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
+           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
+           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
+           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
+           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
+           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
+           !
+           IF( id1 > 0 ) THEN                       ! 'en' exists
+              CALL iom_get( numror, jpdom_autoglo, 'en', en )
+              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
+                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
+                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
+                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
+                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
+                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
+              ELSE                                                 ! one at least array is missing
+                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
+              ENDIF
+           ELSE                                     ! No TKE array found: initialisation
+              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
+              en (:,:,:) = rn_emin * tmask(:,:,:)
+              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
+              !
+              avt_k (:,:,:) = avt (:,:,:)
+              avm_k (:,:,:) = avm (:,:,:)
+              avmu_k(:,:,:) = avmu(:,:,:)
+              avmv_k(:,:,:) = avmv(:,:,:)
+              !
+              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
+           ENDIF
+        ELSE                                   !* Start from rest
+           en(:,:,:) = rn_emin * tmask(:,:,:)
+           DO jk = 1, jpk                           ! set the Kz to the background value
+              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
+              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
+              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
+              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
+           END DO
+        ENDIF
+        !
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+        !                                   ! -------------------
+        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
+        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
+        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
+        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
+        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
+        !
+     ENDIF
+     !
+   END SUBROUTINE tke_rst
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :                                        NO TKE scheme
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
+CONTAINS
+   SUBROUTINE zdf_tke_init           ! Dummy routine
+   END SUBROUTINE zdf_tke_init
+   SUBROUTINE zdf_tke( kt )          ! Dummy routine
+      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
+   END SUBROUTINE zdf_tke
+   SUBROUTINE tke_rst( kt, cdrw )
+     CHARACTER(len=*) ::   cdrw
+     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
+   END SUBROUTINE tke_rst
+#endif
+
+   !!======================================================================
+END MODULE zdftke
