Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 7382)
@@ -41,4 +41,5 @@
          <field id="sstgrad2"     long_name="square of module of sst gradient"                                                        unit="degC2/m2" />
          <field id="sbt"          long_name="sea bottom temperature"                                                                  unit="degC"     />
+         <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="sst_wl"       long_name="Delta SST of warm layer"                                                                 unit="degC"     />
          <field id="sst_cs"       long_name="Delta SST of cool skin"                                                                  unit="degC"     />
@@ -49,5 +50,6 @@
          <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="1e-3" />
+         <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" /> 
@@ -93,5 +95,8 @@
 
          <!-- variables available with key_diaar5 -->
-         <field id="botpres"      long_name="Pressure at sea floor"   standard_name="sea_water_pressure_at_sea_floor"   unit="dbar" />
+         <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 -->
@@ -205,5 +210,9 @@
          <!-- * variable related to ice shelf forcing * -->
          <field id="fwfisf"       long_name="Ice shelf melting"                                            unit="Kg/m2/s"  />
-         <field id="qisf"         long_name="Ice Shelf Heat Flux"                                          unit="W/m2"     />
+         <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"      />
@@ -418,10 +427,16 @@
          <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 eulerian mass transport along i-axis"    standard_name="ocean_mass_x_transport"                          unit="kg/s"        grid_ref="grid_U_3D" />
+         <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="u_salttr"     long_name="ocean eulerian salt transport along i-axis"    standard_name="ocean_salt_x_transport"                          unit="1e-3*kg/s"                        />
+         <field id="uadv_heattr"  long_name="ocean advective heat transport along i-axis"    standard_name="advectice_ocean_heat_x_transport"               unit="W"                                />
+         <field id="uadv_salttr"  long_name="ocean advective salt transport along i-axis"    standard_name="advectice_ocean_salt_x_transport"               unit="1e-3*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="Kg"                                />
+         <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" />
          <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 id="udiff_salttr" long_name="ocean diffusion salt transport along i-axis"   standard_name="ocean_salt_x_transport_due_to_diffusion"         unit="1e-3*kg/s"                                />
       </field_group>
       
@@ -464,6 +479,12 @@
          <field id="v_heattr"     long_name="ocean eulerian heat transport along j-axis"    standard_name="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="1e-3*kg/s"                        />
+         <field id="vadv_heattr"  long_name="ocean advective heat transport along j-axis"   standard_name="advectice_ocean_heat_y_transport"                unit="W"                      />
+         <field id="vadv_salttr"  long_name="ocean advective salt transport along j-axis"   standard_name="advectice_ocean_salt_y_transport"                unit="1e-3*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_x_transport_due_to_bolus_advection"   unit="Kg"                                />
+         <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" />
          <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 id="vdiff_salttr" long_name="ocean diffusion salt transport along j-axis"   standard_name="ocean_salt_y_transport_due_to_diffusion"         unit="1e-3*kg/s"                        />
       </field_group>
       
@@ -625,5 +646,5 @@
 
       <!-- Poleward transport : ptr -->     
-      <field_group id="diaptr" domain_ref="ptr" > 
+      <field_group id="diaptr" domain_ref="ptr" >  
         <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" />
@@ -631,14 +652,14 @@
         <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="degC"     grid_ref="gznl_T_3D" />
-        <field id="zotematl"          long_name="Zonal Mean Temperature : Atlantic"            unit="degC"     grid_ref="gznl_T_3D" />
-        <field id="zotempac"          long_name="Zonal Mean Temperature : Pacific"             unit="degC"     grid_ref="gznl_T_3D" />
-        <field id="zotemind"          long_name="Zonal Mean Temperature : Indian"              unit="degC"     grid_ref="gznl_T_3D" />
-        <field id="zotemipc"          long_name="Zonal Mean Temperature : Pacific+Indian"      unit="degC"     grid_ref="gznl_T_3D" />
-        <field id="zosalglo"          long_name="Zonal Mean Salinity : Global"                 unit="1e-3"     grid_ref="gznl_T_3D" />
-        <field id="zosalatl"          long_name="Zonal Mean Salinity : Atlantic"               unit="1e-3"     grid_ref="gznl_T_3D" />
-        <field id="zosalpac"          long_name="Zonal Mean Salinity : Pacific"                unit="1e-3"     grid_ref="gznl_T_3D" />
-        <field id="zosalind"          long_name="Zonal Mean Salinity : Indian"                 unit="1e-3"     grid_ref="gznl_T_3D" />
-        <field id="zosalipc"          long_name="Zonal Mean Salinity : Pacific+Indian"         unit="1e-3"     grid_ref="gznl_T_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" />
@@ -647,7 +668,53 @@
         <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_group>
 
@@ -668,4 +735,6 @@
       <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" />
@@ -674,6 +743,10 @@
       <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"      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" />
@@ -692,4 +765,47 @@
       <field id="ttrd_atf"      long_name="temperature-trend: asselin time filter"       unit="degC/s" />
       <field id="strd_atf"      long_name="salinity   -trend: asselin time filter"       unit="1e-3/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_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 -->
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 7382)
@@ -6,8 +6,4 @@
    !! History :  3.2  !  2009-11  (S. Masson)  Original code
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
-   !!----------------------------------------------------------------------
-#if defined key_diaar5
-   !!----------------------------------------------------------------------
-   !!   'key_diaar5'  :                           activate ar5 diagnotics
    !!----------------------------------------------------------------------
    !!   dia_ar5       : AR5 diagnostics
@@ -24,4 +20,6 @@
    USE phycst         ! physical constant
    USE in_out_manager  ! I/O manager
+   USE zdfddm
+   USE zdf_oce
 
    IMPLICIT NONE
@@ -29,8 +27,6 @@
 
    PUBLIC   dia_ar5        ! routine called in step.F90 module
-   PUBLIC   dia_ar5_init   ! routine called in opa.F90 module
    PUBLIC   dia_ar5_alloc  ! routine called in nemogcm.F90 module
-
-   LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .TRUE.   ! coupled flag
+   PUBLIC   dia_ar5_hst    ! heat/salt transport
 
    REAL(wp)                         ::   vol0         ! ocean volume (interior domain)
@@ -39,5 +35,10 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity
+
+   LOGICAL  :: l_ar5
       
+   !! * Substitutions
+#  include "zdfddm_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -73,6 +74,8 @@
       INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
       REAL(wp) ::   zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass
+      REAL(wp) ::   zaw, zbw, zrw
       !
       REAL(wp), POINTER, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace 
+      REAL(wp), POINTER, DIMENSION(:,:)     :: zpe                         ! 2D workspace 
       REAL(wp), POINTER, DIMENSION(:,:,:)   :: zrhd , zrhop               ! 3D workspace
       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
@@ -80,123 +83,233 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_ar5')
  
-      CALL wrk_alloc( jpi , jpj              , zarea_ssh , zbotpres )
-      CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    )
-      CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 )
-
-      zarea_ssh(:,:) = area(:,:) * sshn(:,:)
-
-      !                                         ! total volume of liquid seawater
-      zvolssh = SUM( zarea_ssh(:,:) ) 
-      IF( lk_mpp )   CALL mpp_sum( zvolssh )
-      zvol = vol0 + zvolssh
+      IF( kt == nit000 )     CALL dia_ar5_init
+
+      IF( l_ar5 ) THEN
+         CALL wrk_alloc( jpi , jpj              , zarea_ssh , zbotpres )
+         CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    )
+         CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 )
+         zarea_ssh(:,:) = area(:,:) * sshn(:,:)
+      ENDIF
+      !
+      IF( iom_use( 'voltot' ) .OR. iom_use( 'sshtot' )  .OR. iom_use( 'sshdyn' )  ) THEN    
+         !                                         ! total volume of liquid seawater
+         zvolssh = SUM( zarea_ssh(:,:) ) 
+         IF( lk_mpp )   CALL mpp_sum( zvolssh )
+         zvol = vol0 + zvolssh
       
-      CALL iom_put( 'voltot', zvol               )
-      CALL iom_put( 'sshtot', zvolssh / area_tot )
-
-      !                     
-      ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
-      ztsn(:,:,:,jp_sal) = sn0(:,:,:)
-      CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity
-      !
-      zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
-      DO jk = 1, jpkm1
-         zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
-      END DO
-      IF( ln_linssh ) THEN
-         IF( ln_isfcav ) THEN
-            DO ji=1,jpi
-               DO jj=1,jpj
-                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
-               END DO
-            END DO
-         ELSE
-            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
-         END IF
+         CALL iom_put( 'voltot', zvol               )
+         CALL iom_put( 'sshtot', zvolssh / area_tot )
+         CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
+         !
+      ENDIF
+
+      IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) THEN    
+         !                     
+         ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
+         ztsn(:,:,:,jp_sal) = sn0(:,:,:)
+         CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity
+         !
+         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
+         DO jk = 1, jpkm1
+            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+         END DO
+         IF( ln_linssh ) THEN
+            IF( ln_isfcav ) THEN
+               DO ji = 1, jpi
+                  DO jj = 1, jpj
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+                  END DO
+               END DO
+            ELSE
+               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+            END IF
 !!gm
 !!gm   riceload should be added in both ln_linssh=T or F, no?
 !!gm
-      END IF
-      !                                         
-      zarho = SUM( area(:,:) * zbotpres(:,:) ) 
-      IF( lk_mpp )   CALL mpp_sum( zarho )
-      zssh_steric = - zarho / area_tot
-      CALL iom_put( 'sshthster', zssh_steric )
+         END IF
+         !                                         
+         zarho = SUM( area(:,:) * zbotpres(:,:) ) 
+         IF( lk_mpp )   CALL mpp_sum( zarho )
+         zssh_steric = - zarho / area_tot
+         CALL iom_put( 'sshthster', zssh_steric )
       
-      !                                         ! steric sea surface height
-      CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
-      zrhop(:,:,jpk) = 0._wp
-      CALL iom_put( 'rhop', zrhop )
-      !
-      zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
+         !                                         ! steric sea surface height
+         CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
+         zrhop(:,:,jpk) = 0._wp
+         CALL iom_put( 'rhop', zrhop )
+         !
+         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
+         DO jk = 1, jpkm1
+            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+         END DO
+         IF( ln_linssh ) THEN
+            IF ( ln_isfcav ) THEN
+               DO ji = 1,jpi
+                  DO jj = 1,jpj
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+                  END DO
+               END DO
+            ELSE
+               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+            END IF
+         END IF
+         !    
+         zarho = SUM( area(:,:) * zbotpres(:,:) ) 
+         IF( lk_mpp )   CALL mpp_sum( zarho )
+         zssh_steric = - zarho / area_tot
+         CALL iom_put( 'sshsteric', zssh_steric )
+      
+         !                                         ! ocean bottom pressure
+         zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
+         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
+         CALL iom_put( 'botpres', zbotpres )
+         !
+      ENDIF
+
+      IF( iom_use( 'masstot' ) .OR. iom_use( 'temptot' )  .OR. iom_use( 'saltot' )  ) THEN    
+         !                                         ! Mean density anomalie, temperature and salinity
+         ztemp = 0._wp
+         zsal  = 0._wp
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
+                  ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
+                  zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal)
+               END DO
+            END DO
+         END DO
+         IF( ln_linssh ) THEN
+            IF( ln_isfcav ) THEN
+               DO ji = 1, jpi
+                  DO jj = 1, jpj
+                     ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
+                     zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
+                  END DO
+               END DO
+            ELSE
+               ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
+               zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
+            END IF
+         ENDIF
+         IF( lk_mpp ) THEN  
+            CALL mpp_sum( ztemp )
+            CALL mpp_sum( zsal  )
+         END IF
+         !
+         zmass = rau0 * ( zarho + zvol )                 ! total mass of liquid seawater
+         ztemp = ztemp / zvol                            ! potential temperature in liquid seawater
+         zsal  = zsal  / zvol                            ! Salinity of liquid seawater
+         !
+         CALL iom_put( 'masstot', zmass )
+         CALL iom_put( 'temptot', ztemp )
+         CALL iom_put( 'saltot' , zsal  )
+         !
+      ENDIF
+
+      IF( iom_use( 'tnpeo' )) THEN    
+      ! Work done against stratification by vertical mixing
+      ! Exclude points where rn2 is negative as convection kicks in here and
+      ! work is not being done against stratification
+          CALL wrk_alloc( jpi, jpj, zpe )
+          zpe(:,:) = 0._wp
+          IF( lk_zdfddm ) THEN
+             DO ji=1,jpi
+                DO jj=1,jpj
+                   DO jk=1,jpk
+                      zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
+                         &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
+                      !
+                      zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
+                      zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
+                      !
+                      zpe(ji, jj) = zpe(ji, jj) - MIN(0._wp, rn2(ji,jj,jk)) * &
+                           &       grav * (avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
+                           &       - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
+
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO ji = 1, jpi
+                DO jj = 1, jpj
+                   DO jk = 1, jpk
+                       zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          CALL lbc_lnk( zpe, 'T', 1._wp)         
+          CALL iom_put( 'tnpeo', zpe )
+          CALL wrk_dealloc( jpi, jpj, zpe )
+      ENDIF
+      !
+      IF( l_ar5 ) THEN
+        CALL wrk_dealloc( jpi , jpj              , zarea_ssh , zbotpres )
+        CALL wrk_dealloc( jpi , jpj , jpk        , zrhd      , zrhop    )
+        CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn                 )
+      ENDIF
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5')
+      !
+   END SUBROUTINE dia_ar5
+
+   SUBROUTINE dia_ar5_hst( ktra, cptr, pua, pva ) 
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dia_ar5_htr ***
+      !!----------------------------------------------------------------------
+      !! Wrapper for heat transport calculations
+      !! Called from all advection and/or diffusion routines
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
+      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pua   ! 3D input array of advection/diffusion
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
+      !
+      INTEGER    ::  ji, jj, jk
+      REAL(wp), POINTER, DIMENSION(:,:)  :: z2d
+
+    
+
+      CALL wrk_alloc( jpi, jpj, z2d )
+      z2d(:,:) = pua(:,:,1) 
       DO jk = 1, jpkm1
-         zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
-      END DO
-      IF( ln_linssh ) THEN
-         IF ( ln_isfcav ) THEN
-            DO ji=1,jpi
-               DO jj=1,jpj
-                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
-               END DO
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z2d(ji,jj) = z2d(ji,jj) + pua(ji,jj,jk) 
             END DO
-         ELSE
-            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
-         END IF
-      END IF
-      !    
-      zarho = SUM( area(:,:) * zbotpres(:,:) ) 
-      IF( lk_mpp )   CALL mpp_sum( zarho )
-      zssh_steric = - zarho / area_tot
-      CALL iom_put( 'sshsteric', zssh_steric )
-      
-      !                                         ! ocean bottom pressure
-      zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
-      zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
-      CALL iom_put( 'botpres', zbotpres )
-
-      !                                         ! Mean density anomalie, temperature and salinity
-      ztemp = 0._wp
-      zsal  = 0._wp
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
-               ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
-               zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal)
-            END DO
-         END DO
-      END DO
-      IF( ln_linssh ) THEN
-         IF( ln_isfcav ) THEN
-            DO ji=1,jpi
-               DO jj=1,jpj
-                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
-                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
-               END DO
-            END DO
-         ELSE
-            ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
-            zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
-         END IF
-      ENDIF
-      IF( lk_mpp ) THEN  
-         CALL mpp_sum( ztemp )
-         CALL mpp_sum( zsal  )
-      END IF
-      !
-      zmass = rau0 * ( zarho + zvol )                 ! total mass of liquid seawater
-      ztemp = ztemp / zvol                            ! potential temperature in liquid seawater
-      zsal  = zsal  / zvol                            ! Salinity of liquid seawater
-      !
-      CALL iom_put( 'masstot', zmass )
-      CALL iom_put( 'temptot', ztemp )
-      CALL iom_put( 'saltot' , zsal  )
-      !
-      CALL wrk_dealloc( jpi , jpj              , zarea_ssh , zbotpres )
-      CALL wrk_dealloc( jpi , jpj , jpk        , zrhd      , zrhop    )
-      CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn                 )
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5')
-      !
-   END SUBROUTINE dia_ar5
+         END DO
+       END DO
+       CALL lbc_lnk( z2d, 'U', -1. )
+       IF( cptr == 'adv' ) THEN
+          IF( ktra == jp_tem ) CALL iom_put( "uadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in i-direction
+          IF( ktra == jp_sal ) CALL iom_put( "uadv_salttr" , rau0     * z2d )  ! advective salt transport in i-direction
+       ENDIF
+       IF( cptr == 'ldf' ) THEN
+          IF( ktra == jp_tem ) CALL iom_put( "udiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in i-direction
+          IF( ktra == jp_sal ) CALL iom_put( "udiff_salttr" , rau0     * z2d ) ! diffusive salt transport in i-direction
+       ENDIF
+       !
+       z2d(:,:) = pva(:,:,1) 
+       DO jk = 1, jpkm1
+          DO jj = 2, jpjm1
+             DO ji = fs_2, fs_jpim1   ! vector opt.
+                z2d(ji,jj) = z2d(ji,jj) + pva(ji,jj,jk) 
+             END DO
+          END DO
+       END DO
+       CALL lbc_lnk( z2d, 'V', -1. )
+       IF( cptr == 'adv' ) THEN
+          IF( ktra == jp_tem ) CALL iom_put( "vadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in j-direction
+          IF( ktra == jp_sal ) CALL iom_put( "vadv_salttr" , rau0     * z2d )  ! advective salt transport in j-direction
+       ENDIF
+       IF( cptr == 'ldf' ) THEN
+          IF( ktra == jp_tem ) CALL iom_put( "vdiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in j-direction
+          IF( ktra == jp_sal ) CALL iom_put( "vdiff_salttr" , rau0     * z2d ) ! diffusive salt transport in j-direction
+       ENDIF
+          
+       CALL wrk_dealloc( jpi, jpj, z2d )
+
+   END SUBROUTINE dia_ar5_hst
 
 
@@ -217,59 +330,54 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init')
       !
-      CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta )
-      !                                      ! allocate dia_ar5 arrays
-      IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
-
-      area(:,:) = e1e2t(:,:) * tmask_i(:,:)
-
-      area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot )
-
-      vol0        = 0._wp
-      thick0(:,:) = 0._wp
-      DO jk = 1, jpkm1
-         vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
-         thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( vol0 )
-
-
-      CALL iom_open ( 'sali_ref_clim_monthly', inum )
-      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  )
-      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 )
-      CALL iom_close( inum )
-
-      sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )        
-      sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
-      IF( ln_zps ) THEN               ! z-coord. partial steps
-         DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
-            DO ji = 1, jpi
-               ik = mbkt(ji,jj)
-               IF( ik > 1 ) THEN
-                  zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-                  sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
-               ENDIF
+      l_ar5 = .FALSE.
+      IF(   iom_use( 'voltot'  ) .OR. iom_use( 'sshtot'    )  .OR. iom_use( 'sshdyn' )  .OR.  & 
+         &  iom_use( 'masstot' ) .OR. iom_use( 'temptot'   )  .OR. iom_use( 'saltot' ) .OR.  &    
+         &  iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) L_ar5 = .TRUE.
+  
+      IF( l_ar5 ) THEN
+         !
+         CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta )
+         !                                      ! allocate dia_ar5 arrays
+         IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
+
+         area(:,:) = e1e2t(:,:) * tmask_i(:,:)
+
+         area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot )
+
+         vol0        = 0._wp
+         thick0(:,:) = 0._wp
+         DO jk = 1, jpkm1
+            vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
+            thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( vol0 )
+
+
+         CALL iom_open ( 'sali_ref_clim_monthly', inum )
+         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  )
+         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 )
+         CALL iom_close( inum )
+
+         sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )        
+         sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
+         IF( ln_zps ) THEN               ! z-coord. partial steps
+            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
+               DO ji = 1, jpi
+                  ik = mbkt(ji,jj)
+                  IF( ik > 1 ) THEN
+                     zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                     sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
+                  ENDIF
+               END DO
             END DO
-         END DO
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta )
+         ENDIF
+         !
+         CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta )
+         !
+      ENDIF
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')
       !
    END SUBROUTINE dia_ar5_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                         NO diaar5
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .FALSE.   ! coupled flag
-CONTAINS
-   SUBROUTINE dia_ar5_init    ! Dummy routine
-   END SUBROUTINE dia_ar5_init
-   SUBROUTINE dia_ar5( kt )   ! Empty routine
-      INTEGER ::   kt
-      WRITE(*,*) 'dia_ar5: You should not have seen this print! error?', kt
-   END SUBROUTINE dia_ar5
-#endif
 
    !!======================================================================
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 7382)
@@ -9,4 +9,5 @@
    !!            3.3  ! 2010-10  (G. Madec)  dynamical allocation
    !!            3.6  ! 2014-12  (C. Ethe) use of IOM
+   !!            3.6  ! 2016-06  (T. Graham) Addition of diagnostics for CMIP6
    !!----------------------------------------------------------------------
 
@@ -38,13 +39,15 @@
    PUBLIC   dia_ptr_init   ! call in step module
    PUBLIC   dia_ptr        ! call in step module
+   PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines
 
    !                                  !!** namelist  namptr  **
-   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf   !: Heat TRansports (adv, diff, overturn.)
-   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf   !: Salt TRansports (adv, diff, overturn.)
-   
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_adv, htr_ldf, htr_eiv   !: Heat TRansports (adv, diff, Bolus.)
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   str_adv, str_ldf, str_eiv   !: Salt TRansports (adv, diff, Bolus.)
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_ove, str_ove   !: heat Salt TRansports ( overturn.)
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_btr, str_btr   !: heat Salt TRansports ( barotropic )
 
    LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F)
    LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation
-   INTEGER         ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T) 
+   INTEGER, PUBLIC ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T) 
 
    REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
@@ -75,10 +78,20 @@
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::   zv, zsfc               ! local scalar
+      REAL(wp) ::   zsfc,zvfc               ! local scalar
       REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace
-      CHARACTER( len = 10 )  :: cl1
+      REAL(wp), DIMENSION(jpj)     ::  vsum   ! 1D workspace
+      REAL(wp), DIMENSION(jpj,jpts)     ::  tssum   ! 1D workspace
+ 
+      !
+      !overturning calculation
+      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   sjk  , r1_sjk ! i-mean i-k-surface and its inverse
+      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   v_msf, sn_jk  , tn_jk ! i-mean T and S, j-Stream-Function
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace
+
+
+      CHARACTER( len = 12 )  :: cl1
       !!----------------------------------------------------------------------
       !
@@ -109,4 +122,122 @@
             END DO
          ENDIF
+         IF( iom_use("sopstove") .OR. iom_use("sophtove") .OR. iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN
+            ! define fields multiplied by scalar
+            zmask(:,:,:) = 0._wp
+            zts(:,:,:,:) = 0._wp
+            zvn(:,:,:) = 0._wp
+            DO jk = 1, jpkm1
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpi
+                     zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk)
+                     zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
+                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
+                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
+                     zvn(ji,jj,jk)        = vn(ji,jj,jk)         * zvfc
+                  ENDDO
+               ENDDO
+             ENDDO
+         ENDIF
+         IF( iom_use("sopstove") .OR. iom_use("sophtove") ) THEN
+             sjk(:,:,1) = ptr_sjk( zmask(:,:,:), btmsk(:,:,1) )
+             r1_sjk(:,:,1) = 0._wp
+             WHERE( sjk(:,:,1) /= 0._wp )   r1_sjk(:,:,1) = 1._wp / sjk(:,:,1)
+
+             ! i-mean T and S, j-Stream-Function, global
+             tn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_tem) ) * r1_sjk(:,:,1)
+             sn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_sal) ) * r1_sjk(:,:,1)
+             v_msf(:,:,1) = ptr_sjk( zvn(:,:,:) )
+
+             htr_ove(:,1) = SUM( v_msf(:,:,1)*tn_jk(:,:,1) ,2 )
+             str_ove(:,1) = SUM( v_msf(:,:,1)*sn_jk(:,:,1) ,2 )
+
+             z2d(1,:) = htr_ove(:,1) * rc_pwatt        !  (conversion in PW)
+             DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+             ENDDO
+             cl1 = 'sophtove'
+             CALL iom_put( TRIM(cl1), z2d )
+             z2d(1,:) = str_ove(:,1) * rc_ggram        !  (conversion in Gg)
+             DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+             ENDDO
+             cl1 = 'sopstove'
+             CALL iom_put( TRIM(cl1), z2d )
+             IF( ln_subbas ) THEN
+                DO jn = 2, nptr
+                    sjk(:,:,jn) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
+                    r1_sjk(:,:,jn) = 0._wp
+                    WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
+
+                    ! i-mean T and S, j-Stream-Function, basin
+                    tn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
+                    sn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
+                    v_msf(:,:,jn) = ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) ) 
+                    htr_ove(:,jn) = SUM( v_msf(:,:,jn)*tn_jk(:,:,jn) ,2 )
+                    str_ove(:,jn) = SUM( v_msf(:,:,jn)*sn_jk(:,:,jn) ,2 )
+
+                    z2d(1,:) = htr_ove(:,jn) * rc_pwatt !  (conversion in PW)
+                    DO ji = 1, jpi
+                        z2d(ji,:) = z2d(1,:)
+                    ENDDO
+                    cl1 = TRIM('sophtove_'//clsubb(jn))
+                    CALL iom_put( cl1, z2d )
+                    z2d(1,:) = str_ove(:,jn) * rc_ggram        ! (conversion in Gg)
+                    DO ji = 1, jpi
+                        z2d(ji,:) = z2d(1,:)
+                    ENDDO
+                    cl1 = TRIM('sopstove_'//clsubb(jn))
+                    CALL iom_put( cl1, z2d )
+                END DO
+             ENDIF
+         ENDIF
+         IF( iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN
+         ! Calculate barotropic heat and salt transport here 
+             sjk(:,1,1) = ptr_sj( zmask(:,:,:), btmsk(:,:,1) )
+             r1_sjk(:,1,1) = 0._wp
+             WHERE( sjk(:,1,1) /= 0._wp )   r1_sjk(:,1,1) = 1._wp / sjk(:,1,1)
+            
+            vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,1))
+            tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,1) )
+            tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,1) )
+            htr_btr(:,1) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,1)
+            str_btr(:,1) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,1)
+            z2d(1,:) = htr_btr(:,1) * rc_pwatt        !  (conversion in PW)
+            DO ji = 2, jpi
+               z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sophtbtr'
+            CALL iom_put( TRIM(cl1), z2d )
+            z2d(1,:) = str_btr(:,1) * rc_ggram        !  (conversion in Gg)
+            DO ji = 2, jpi
+              z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sopstbtr'
+            CALL iom_put( TRIM(cl1), z2d )
+            IF( ln_subbas ) THEN
+                DO jn = 2, nptr
+                    sjk(:,1,jn) = ptr_sj( zmask(:,:,:), btmsk(:,:,jn) )
+                    r1_sjk(:,1,jn) = 0._wp
+                    WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
+                    vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,jn))
+                    tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
+                    tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
+                    htr_btr(:,jn) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,jn)
+                    str_btr(:,jn) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,jn)
+                    z2d(1,:) = htr_btr(:,jn) * rc_pwatt !  (conversion in PW)
+                    DO ji = 1, jpi
+                        z2d(ji,:) = z2d(1,:)
+                    ENDDO
+                    cl1 = TRIM('sophtbtr_'//clsubb(jn))
+                    CALL iom_put( cl1, z2d )
+                    z2d(1,:) = str_btr(:,jn) * rc_ggram        ! (conversion in Gg)
+                    DO ji = 1, jpi
+                        z2d(ji,:) = z2d(1,:)
+                    ENDDO
+                    cl1 = TRIM('sopstbtr_'//clsubb(jn))
+                    CALL iom_put( cl1, z2d )
+               ENDDO
+            ENDIF !ln_subbas
+         ENDIF !iom_use("sopstbtr....)
          !
       ELSE
@@ -148,5 +279,5 @@
          !                                ! Advective and diffusive heat and salt transport
          IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN   
-            z2d(1,:) = htr_adv(:) * rc_pwatt        !  (conversion in PW)
+            z2d(1,:) = htr_adv(:,1) * rc_pwatt        !  (conversion in PW)
             DO ji = 1, jpi
                z2d(ji,:) = z2d(1,:)
@@ -154,5 +285,5 @@
             cl1 = 'sophtadv'                 
             CALL iom_put( TRIM(cl1), z2d )
-            z2d(1,:) = str_adv(:) * rc_ggram        ! (conversion in Gg)
+            z2d(1,:) = str_adv(:,1) * rc_ggram        ! (conversion in Gg)
             DO ji = 1, jpi
                z2d(ji,:) = z2d(1,:)
@@ -160,8 +291,24 @@
             cl1 = 'sopstadv'
             CALL iom_put( TRIM(cl1), z2d )
+            IF( ln_subbas ) THEN
+              DO jn=2,nptr
+               z2d(1,:) = htr_adv(:,jn) * rc_pwatt        !  (conversion in PW)
+               DO ji = 1, jpi
+                 z2d(ji,:) = z2d(1,:)
+               ENDDO
+               cl1 = TRIM('sophtadv_'//clsubb(jn))                 
+               CALL iom_put( cl1, z2d )
+               z2d(1,:) = str_adv(:,jn) * rc_ggram        ! (conversion in Gg)
+               DO ji = 1, jpi
+                  z2d(ji,:) = z2d(1,:)
+               ENDDO
+               cl1 = TRIM('sopstadv_'//clsubb(jn))                 
+               CALL iom_put( cl1, z2d )              
+              ENDDO
+            ENDIF
          ENDIF
          !
          IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN   
-            z2d(1,:) = htr_ldf(:) * rc_pwatt        !  (conversion in PW) 
+            z2d(1,:) = htr_ldf(:,1) * rc_pwatt        !  (conversion in PW) 
             DO ji = 1, jpi
                z2d(ji,:) = z2d(1,:)
@@ -169,5 +316,5 @@
             cl1 = 'sophtldf'
             CALL iom_put( TRIM(cl1), z2d )
-            z2d(1,:) = str_ldf(:) * rc_ggram        !  (conversion in Gg)
+            z2d(1,:) = str_ldf(:,1) * rc_ggram        !  (conversion in Gg)
             DO ji = 1, jpi
                z2d(ji,:) = z2d(1,:)
@@ -175,4 +322,51 @@
             cl1 = 'sopstldf'
             CALL iom_put( TRIM(cl1), z2d )
+            IF( ln_subbas ) THEN
+              DO jn=2,nptr
+               z2d(1,:) = htr_ldf(:,jn) * rc_pwatt        !  (conversion in PW)
+               DO ji = 1, jpi
+                 z2d(ji,:) = z2d(1,:)
+               ENDDO
+               cl1 = TRIM('sophtldf_'//clsubb(jn))                 
+               CALL iom_put( cl1, z2d )
+               z2d(1,:) = str_ldf(:,jn) * rc_ggram        ! (conversion in Gg)
+               DO ji = 1, jpi
+                  z2d(ji,:) = z2d(1,:)
+               ENDDO
+               cl1 = TRIM('sopstldf_'//clsubb(jn))                 
+               CALL iom_put( cl1, z2d )              
+              ENDDO
+            ENDIF
+         ENDIF
+
+         IF( iom_use("sophteiv") .OR. iom_use("sopsteiv") ) THEN 
+            z2d(1,:) = htr_eiv(:,1) * rc_pwatt        !  (conversion in PW) 
+            DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sophteiv'
+            CALL iom_put( TRIM(cl1), z2d )
+            z2d(1,:) = str_eiv(:,1) * rc_ggram        !  (conversion in Gg)
+            DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sopsteiv'
+            CALL iom_put( TRIM(cl1), z2d )
+            IF( ln_subbas ) THEN
+               DO jn=2,nptr
+                  z2d(1,:) = htr_eiv(:,jn) * rc_pwatt        !  (conversion in PW)
+                  DO ji = 1, jpi
+                     z2d(ji,:) = z2d(1,:)
+                  ENDDO
+                  cl1 = TRIM('sophteiv_'//clsubb(jn))                 
+                  CALL iom_put( cl1, z2d )
+                  z2d(1,:) = str_eiv(:,jn) * rc_ggram        ! (conversion in Gg)
+                  DO ji = 1, jpi
+                     z2d(ji,:) = z2d(1,:)
+                  ENDDO
+                  cl1 = TRIM('sopsteiv_'//clsubb(jn)) 
+                  CALL iom_put( cl1, z2d )              
+               ENDDO
+            ENDIF
          ENDIF
          !
@@ -254,6 +448,9 @@
          ! Initialise arrays to zero because diatpr is called before they are first calculated
          ! Note that this means diagnostics will not be exactly correct when model run is restarted.
-         htr_adv(:) = 0._wp  ;  str_adv(:) =  0._wp  
-         htr_ldf(:) = 0._wp  ;  str_ldf(:) =  0._wp 
+         htr_adv(:,:) = 0._wp  ;  str_adv(:,:) =  0._wp 
+         htr_ldf(:,:) = 0._wp  ;  str_ldf(:,:) =  0._wp 
+         htr_eiv(:,:) = 0._wp  ;  str_eiv(:,:) =  0._wp 
+         htr_ove(:,:) = 0._wp  ;   str_ove(:,:) =  0._wp
+         htr_btr(:,:) = 0._wp  ;   str_btr(:,:) =  0._wp
          !
       ENDIF 
@@ -261,4 +458,71 @@
    END SUBROUTINE dia_ptr_init
 
+   SUBROUTINE dia_ptr_hst( ktra, cptr, pva ) 
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dia_ptr_hst ***
+      !!----------------------------------------------------------------------
+      !! Wrapper for heat and salt transport calculations to calculate them for each basin
+      !! Called from all advection and/or diffusion routines
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
+      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
+      INTEGER                                        :: jn    !
+
+      IF( cptr == 'adv' ) THEN
+         IF( ktra == jp_tem )  htr_adv(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_sal )  str_adv(:,1) = ptr_sj( pva(:,:,:) )
+      ENDIF
+      IF( cptr == 'ldf' ) THEN
+         IF( ktra == jp_tem )  htr_ldf(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_sal )  str_ldf(:,1) = ptr_sj( pva(:,:,:) )
+      ENDIF
+      IF( cptr == 'eiv' ) THEN
+         IF( ktra == jp_tem )  htr_eiv(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_sal )  str_eiv(:,1) = ptr_sj( pva(:,:,:) )
+      ENDIF
+      !
+      IF( ln_subbas ) THEN
+         !
+         IF( cptr == 'adv' ) THEN
+             IF( ktra == jp_tem ) THEN 
+                DO jn = 2, nptr
+                   htr_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                END DO
+             ENDIF
+             IF( ktra == jp_sal ) THEN 
+                DO jn = 2, nptr
+                   str_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                END DO
+             ENDIF
+         ENDIF
+         IF( cptr == 'ldf' ) THEN
+             IF( ktra == jp_tem ) THEN 
+                DO jn = 2, nptr
+                    htr_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                 END DO
+             ENDIF
+             IF( ktra == jp_sal ) THEN 
+                DO jn = 2, nptr
+                   str_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                END DO
+             ENDIF
+         ENDIF
+         IF( cptr == 'eiv' ) THEN
+             IF( ktra == jp_tem ) THEN 
+                DO jn = 2, nptr
+                    htr_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                 END DO
+             ENDIF
+             IF( ktra == jp_sal ) THEN 
+                DO jn = 2, nptr
+                   str_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                END DO
+             ENDIF
+         ENDIF
+         !
+      ENDIF
+   END SUBROUTINE dia_ptr_hst
+
 
    FUNCTION dia_ptr_alloc()
@@ -271,7 +535,10 @@
       ierr(:) = 0
       !
-      ALLOCATE( btmsk(jpi,jpj,nptr) ,           &
-         &      htr_adv(jpj) , str_adv(jpj) ,   &
-         &      htr_ldf(jpj) , str_ldf(jpj) , STAT=ierr(1)  )
+      ALLOCATE( btmsk(jpi,jpj,nptr) ,              &
+         &      htr_adv(jpj,nptr) , str_adv(jpj,nptr) ,   &
+         &      htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) ,   &
+         &      htr_ove(jpj,nptr) , str_ove(jpj,nptr) ,   &
+         &      htr_btr(jpj,nptr) , str_btr(jpj,nptr) ,   &
+         &      htr_ldf(jpj,nptr) , str_ldf(jpj,nptr) , STAT=ierr(1)  )
          !
       ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2))
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7382)
@@ -302,10 +302,13 @@
       CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
       !
-      IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
+      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
          z3d(:,:,jpk) = 0.e0
+         z2d(:,:) = 0.e0
          DO jk = 1, jpkm1
             z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
+            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
          END DO
          CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
+         CALL iom_put( "u_masstr_vint", z2d )             ! mass transport in i-direction vertical sum
       ENDIF
       
@@ -370,4 +373,34 @@
          CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
       ENDIF
+
+      ! Vertical integral of temperature
+      IF( iom_use("tosmint") ) THEN
+         z2d(:,:)=0._wp
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( z2d, 'T', -1. )
+         CALL iom_put( "tosmint", z2d ) 
+      ENDIF
+
+      ! Vertical integral of salinity
+      IF( iom_use("somint") ) THEN
+         z2d(:,:)=0._wp
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( z2d, 'T', -1. )
+         CALL iom_put( "somint", z2d ) 
+      ENDIF
+
+      CALL iom_put( "bn2", rn2 )  !Brunt-Vaisala buoyancy frequency (N^2)
       !
       CALL wrk_dealloc( jpi , jpj      , z2d )
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 7382)
@@ -24,5 +24,5 @@
    USE ldfslp          ! lateral diffusion: slope of iso-neutral surfaces
    USE ldfc1d_c2d      ! lateral diffusion: 1D & 2D cases 
-   USE diaar5, ONLY:   lk_diaar5
+   USE diaptr
    !
    USE trc_oce, ONLY: lk_offline ! offline flag
@@ -730,36 +730,80 @@
       CALL iom_put( "woce_eiv", zw3d )
       !
+      !      
+      !
+      CALL wrk_alloc( jpi,jpj,   zw2d )
+      !
+      zztmp = 0.5_wp * rau0 * rcp 
+      IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
+        zw2d(:,:)   = 0._wp 
+        zw3d(:,:,:) = 0._wp 
+        DO jk = 1, jpkm1
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
+                    &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) ) 
+                 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
+              END DO
+           END DO
+        END DO
+        CALL lbc_lnk( zw2d, 'U', -1. )
+        CALL lbc_lnk( zw3d, 'U', -1. )
+        CALL iom_put( "ueiv_heattr"  , zztmp * zw2d )                  ! heat transport in i-direction
+        CALL iom_put( "ueiv_heattr3d", zztmp * zw3d )                  ! heat transport in i-direction
+      ENDIF
+      zw2d(:,:)   = 0._wp 
+      zw3d(:,:,:) = 0._wp 
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
+                  &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) ) 
+               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zw2d, 'V', -1. )
+      CALL iom_put( "veiv_heattr", zztmp * zw2d )                  !  heat transport in j-direction
+      CALL iom_put( "veiv_heattr", zztmp * zw3d )                  !  heat transport in j-direction
+      !
+      IF( ln_diaptr )  CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
+      !
+      zztmp = 0.5_wp * 0.5
+      IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
+        zw2d(:,:) = 0._wp 
+        zw3d(:,:,:) = 0._wp 
+        DO jk = 1, jpkm1
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
+                    &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji+1,jj,jk,jp_sal) ) 
+                 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
+              END DO
+           END DO
+        END DO
+        CALL lbc_lnk( zw2d, 'U', -1. )
+        CALL lbc_lnk( zw3d, 'U', -1. )
+        CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction
+        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction
+      ENDIF
+      zw2d(:,:) = 0._wp 
+      zw3d(:,:,:) = 0._wp 
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
+                  &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji,jj+1,jk,jp_sal) ) 
+               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zw2d, 'V', -1. )
+      CALL iom_put( "veiv_salttr", zztmp * zw2d )                  !  salt transport in j-direction
+      CALL iom_put( "veiv_salttr", zztmp * zw3d )                  !  salt transport in j-direction
+      !
+      IF( ln_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
+      !
+      CALL wrk_dealloc( jpi,jpj,   zw2d )
       CALL wrk_dealloc( jpi,jpj,jpk,   zw3d )
-      !      
-      !
-      IF( lk_diaar5 ) THEN                               !==  eiv heat transport: calculate and output  ==!
-         CALL wrk_alloc( jpi,jpj,   zw2d )
-         !
-         zztmp = 0.5_wp * rau0 * rcp 
-         zw2d(:,:) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zw2d(ji,jj) = zw2d(ji,jj) + zztmp * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
-                     &                              * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) ) 
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( zw2d, 'U', -1. )
-         CALL iom_put( "ueiv_heattr", zw2d )                  ! heat transport in i-direction
-         zw2d(:,:) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zw2d(ji,jj) = zw2d(ji,jj) + zztmp * ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
-                     &                              * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) ) 
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( zw2d, 'V', -1. )
-         CALL iom_put( "veiv_heattr", zw2d )                  !  heat transport in i-direction
-         !
-         CALL wrk_dealloc( jpi,jpj,   zw2d )
-      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop( 'ldf_eiv_dia')      
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 7382)
@@ -90,6 +90,9 @@
       INTEGER, INTENT( in ) :: kt                   ! ocean time step
       !
-      INTEGER               :: ji, jj               ! loop index
+      INTEGER               :: ji, jj, jk           ! loop index
+      INTEGER               :: ikt, ikb             ! loop index
       REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 
+      REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d
+      REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d
       !!---------------------------------------------------------------------
       !
@@ -161,6 +164,45 @@
          CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
          CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
-         CALL lbc_lnk(fwfisf(:,:)         ,'T',1.)
-         CALL lbc_lnk(qisf(:,:)           ,'T',1.)
+         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
+         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
+
+!=============================================================================================================================================
+         IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
+            CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
+            CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        )
+
+            zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s)
+            zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2)
+            zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2)
+            zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2)
+
+            DO jj = 1,jpj
+               DO ji = 1,jpi
+                  ikt = misfkt(ji,jj)
+                  ikb = misfkb(ji,jj)
+                  DO jk = ikt, ikb - 1
+                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
+                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
+                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
+                  END DO
+                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * e3t_n(ji,jj,jk)
+                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * e3t_n(ji,jj,jk)
+                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * e3t_n(ji,jj,jk)
+               END DO
+            END DO
+
+            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
+            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
+            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
+            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  ))
+
+            CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
+            CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        )
+         END IF
+
+         ! output
+         CALL iom_put('qlatisf'  , qisf)
+         CALL iom_put('fwfisf', fwfisf)
+!=============================================================================================================================================
 
          IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    !
@@ -177,8 +219,4 @@
          END IF
          ! 
-         ! output
-         CALL iom_put('qisf'  , qisf)
-         CALL iom_put('fwfisf', fwfisf)
-
          ! deallocation
          CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  )
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7382)
@@ -26,4 +26,6 @@
    USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff.
    USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! trends manager: tracers 
    !
    USE in_out_manager ! I/O manager
@@ -86,4 +88,5 @@
       INTEGER ::   jk   ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
       !!----------------------------------------------------------------------
       !
@@ -127,4 +130,10 @@
       IF( ln_diaptr )   CALL dia_ptr( zvn )                                    ! diagnose the effective MSF 
 !!gm ???
+      !
+      IF( l_trdtra )   THEN                    !* Save ta and sa trends
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+      ENDIF
       !
       SELECT CASE ( nadv )                      !==  compute advection trend and add it to general trend  ==!
@@ -145,5 +154,14 @@
       END SELECT
       !
-      !                                         ! print mean trends (used for debugging)
+      IF( l_trdtra )   THEN                      ! save the advective trends for further diagnostics
+         DO jk = 1, jpkm1
+            ztrdt(:,:,jk) = tsa(:,:,jk,jp_tem) - ztrdt(:,:,jk)
+            ztrds(:,:,jk) = tsa(:,:,jk,jp_sal) - ztrds(:,:,jk)
+         END DO
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_totad, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_totad, ztrds )
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
+      ENDIF
+      !                                              ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv  - Ta: ', mask1=tmask,               &
          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90	(revision 7382)
@@ -18,4 +18,5 @@
    USE trdtra         ! trends manager: tracers 
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
    !
    USE in_out_manager ! I/O manager
@@ -32,4 +33,8 @@
    
    REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
+
+   LOGICAL :: l_trd   ! flag to compute trends
+   LOGICAL :: l_ptr   ! flag to compute poleward transport
+   LOGICAL :: l_hst   ! flag to compute heat/salt transport
 
    !! * Substitutions
@@ -88,4 +93,12 @@
       ENDIF
       !
+      l_trd = .FALSE.
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )        l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
+      !
       !                    
       zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers
@@ -184,14 +197,13 @@
          END DO
          !                             ! trend diagnostics
-         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN
+         IF( l_trd ) THEN
             CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
             CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
             CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
          END IF
-         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-           IF( jn == jp_tem )   htr_adv(:) = ptr_sj( zwy(:,:,:) )
-           IF( jn == jp_sal )   str_adv(:) = ptr_sj( zwy(:,:,:) )
-         ENDIF
+         !                                 ! "Poleward" heat and salt transports 
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         !                                 !  heat and salt transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
          !
       END DO
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7382)
@@ -20,6 +20,9 @@
    USE trdtra         ! tracers trends
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
+   USE phycst, ONLY: rau0_rcp
    !
    USE in_out_manager ! I/O manager
+   USE iom
    USE lib_mpp        ! MPP library
    USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
@@ -36,4 +39,6 @@
 
    LOGICAL  ::   l_trd   ! flag to compute trends
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
    REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
 
@@ -80,5 +85,6 @@
       REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrdz
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrdz, zptry
+      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d
       !!----------------------------------------------------------------------
       !
@@ -94,11 +100,20 @@
       !
       l_trd = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
-      !
-      IF( l_trd )  THEN
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA'   .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )     l_trd = .TRUE.
+      IF(   cdtype == 'TRA'   .AND. ln_diaptr )                                              l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
+      !
+      IF( l_trd .OR. l_hst )  THEN
          CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
          ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
       ENDIF
       !
+      IF( l_ptr ) THEN  
+         CALL wrk_alloc( jpi, jpj, jpk, zptry )
+         zptry(:,:,:) = 0._wp
+      ENDIF
       !                          ! surface & bottom value : flux set to zero one for all
       zwz(:,:, 1 ) = 0._wp            
@@ -161,12 +176,9 @@
          CALL lbc_lnk( zwi, 'T', 1. )  ! Lateral boundary conditions on zwi  (unchanged sign)
          !                
-         IF( l_trd )  THEN             ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
             ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
          END IF
          !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
-           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
-         ENDIF
+         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:) 
          !
          !        !==  anti-diffusive flux : high order minus low order  ==!
@@ -292,24 +304,30 @@
          END DO
          !
-         IF( l_trd ) THEN     ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes)
             ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
             ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
             ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
-            !
+         ENDIF
+            !
+         IF( l_trd ) THEN 
             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
             CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
             !
-            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
          END IF
-         !                    ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-           IF( jn == jp_tem )   htr_adv(:) = htr_adv(:) + ptr_sj( zwy(:,:,:) )
-           IF( jn == jp_sal )   str_adv(:) = str_adv(:) + ptr_sj( zwy(:,:,:) )
+         !                                !  heat/salt transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
+
+         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( l_ptr ) THEN  
+            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
          ENDIF
          !
       END DO                     ! end of tracer loop
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
+                              CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
+      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct')
@@ -357,4 +375,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdx, ztrdy, ztrdz
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zptry
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrs
       !!----------------------------------------------------------------------
@@ -373,11 +392,20 @@
       !
       l_trd = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
-      !
-      IF( l_trd )  THEN
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
+      !
+      IF( l_trd .OR. l_hst )  THEN
          CALL wrk_alloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
          ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp
       ENDIF
       !
+      IF( l_ptr ) THEN  
+         CALL wrk_alloc( jpi, jpj,jpk, zptry )
+         zptry(:,:,:) = 0._wp
+      ENDIF
       zwi(:,:,:) = 0._wp
       z_rzts = 1._wp / REAL( kn_fct_zts, wp )
@@ -445,12 +473,9 @@
          CALL lbc_lnk( zwi, 'T', 1. )     ! Lateral boundary conditions on zwi  (unchanged sign)
          !                
-         IF( l_trd )  THEN                ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd .OR. l_hst )  THEN                ! trend diagnostics (contribution of upstream fluxes)
             ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
          END IF
          !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
-           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
-         ENDIF
+         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:)
 
          ! 3. anti-diffusive flux : high order minus low order
@@ -568,27 +593,33 @@
          END DO
 
-         !                                 ! trend diagnostics (contribution of upstream fluxes)
-         IF( l_trd )  THEN 
+        !
+         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes)
             ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
             ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
             ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
-            !
-            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )   
-            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )  
-            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
-            !
-            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
+         ENDIF
+            !
+         IF( l_trd ) THEN 
+            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
+            !
          END IF
-         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:)
-           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:)
+         !                                             ! heat/salt transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
+
+         !                                            ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( l_ptr ) THEN  
+            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
          ENDIF
          !
       END DO
       !
-      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
-      CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
-      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
+                              CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
+                              CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
+                              CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
+      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts')
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7382)
@@ -23,5 +23,8 @@
    USE sbcrnf         ! river runoffs
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
+
    !
+   USE iom
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
@@ -40,4 +43,8 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
    
+   LOGICAL  ::   l_trd   ! flag to compute trends
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -116,4 +123,12 @@
       ENDIF 
       !      
+      l_trd = .FALSE.
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
+      !
       DO jn = 1, kjpt            !==  loop over the tracers  ==!
          !
@@ -192,14 +207,12 @@
          END DO        
          !                                ! trend diagnostics
-         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   &
-            &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN
+         IF( l_trd )  THEN
             CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
             CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
          END IF
-         !                                ! "Poleward" heat and salt transports
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
-            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
-         ENDIF
+         !                                 ! "Poleward" heat and salt transports 
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         !                                 !  heat transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
          !
          !                          !* Vertical advective fluxes
@@ -262,7 +275,5 @@
          END DO
          !                                ! send trends for diagnostic
-         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     &
-            &( cdtype == 'TRC' .AND. l_trdtrc )      )   &
-            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
+         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
          !
       END DO                     ! end of tracer loop
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 7382)
@@ -34,6 +34,9 @@
    PUBLIC   tra_adv_qck   ! routine called by step.F90
 
-   LOGICAL  :: l_trd           ! flag to compute trends
    REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
+
+   LOGICAL  ::   l_trd   ! flag to compute trends
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+
 
    !! * Substitutions
@@ -103,5 +106,8 @@
       !
       l_trd = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
+      !
       !
       !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
@@ -224,5 +230,5 @@
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
+         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
          !
       END DO
@@ -347,10 +353,7 @@
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
+         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
-           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
-         ENDIF
+         IF( l_ptr )                     CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
          !
       END DO
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 7382)
@@ -19,5 +19,8 @@
    USE trdtra         ! trends manager: tracers 
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
+
    !
+   USE iom
    USE lib_mpp        ! I/O library
    USE lbclnk         ! ocean lateral boundary condition (or mpp link)
@@ -32,5 +35,8 @@
    PUBLIC   tra_adv_ubs   ! routine called by traadv module
 
-   LOGICAL :: l_trd  ! flag to compute trends or not
+   LOGICAL :: l_trd   ! flag to compute trends
+   LOGICAL :: l_ptr   ! flag to compute poleward transport
+   LOGICAL :: l_hst   ! flag to compute heat transport
+
 
    !! * Substitutions
@@ -109,5 +115,10 @@
       !
       l_trd = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
+      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
       !
       ztw (:,:, 1 ) = 0._wp      ! surface & bottom value : set to zero for all tracers
@@ -176,9 +187,10 @@
              CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) )
          END IF
-         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( ztv(:,:,:) )
-            IF( jn == jp_sal )  str_adv(:) = ptr_sj( ztv(:,:,:) )
-         ENDIF
+         !     
+         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
+         !                                !  heati/salt transport
+         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
+         !
          !
          !                       !== vertical advective trend  ==!
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7382)
@@ -24,4 +24,5 @@
    USE ldfslp         ! iso-neutral slopes
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
    !
    USE in_out_manager ! I/O manager
@@ -36,4 +37,7 @@
 
    PUBLIC   tra_ldf_iso   ! routine called by step.F90
+
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat transport
 
    !! * Substitutions
@@ -107,7 +111,4 @@
       REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
       REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
-#if defined key_diaar5
-      REAL(wp) ::   zztmp   ! local scalar
-#endif
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zdkt, zdk1t, z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdit, zdjt, zftu, zftv, ztfw 
@@ -127,5 +128,12 @@
          ah_wslp2(:,:,:) = 0._wp
       ENDIF
-      !                                               ! set time step size (Euler/Leapfrog)
+      !   
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
+      !
+      !                                            ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
       ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
@@ -369,41 +377,8 @@
             !
             !                             ! "Poleward" diffusive heat or salt transports (T-S case only)
-            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
-               ! note sign is reversed to give down-gradient diffusive transports (#1043)
-               IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) )
-               IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) )
-            ENDIF
-            !
-            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
-              !
-              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
-                  z2d(:,:) = zftu(ji,jj,1) 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
-                        END DO
-                     END DO
-                  END DO
-!!gm CAUTION I think there is an error of sign when using BLP operator....
-!!gm         a multiplication by zsign is required (to be checked twice !)
-                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
-                  CALL lbc_lnk( z2d, 'U', -1. )
-                  CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
-                  !
-                  z2d(:,:) = zftv(ji,jj,1) 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
-                        END DO
-                     END DO
-                  END DO
-                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
-                  CALL lbc_lnk( z2d, 'V', -1. )
-                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
-               END IF
-               !
-            ENDIF
+               ! note sign is reversed to give down-gradient diffusive transports )
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -zftv(:,:,:)  )
+            !                          ! Diffusive heat transports
+            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) )
             !
          ENDIF                                                    !== end pass selection  ==!
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90	(revision 7382)
@@ -17,4 +17,5 @@
    USE traldf_triad   ! iso-neutral lateral diffusion (triad    operator)     (tra_ldf_triad routine)
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
    USE trc_oce        ! share passive tracers/Ocean variables
    USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
@@ -25,4 +26,5 @@
    USE timing         ! Timing
    USE wrk_nemo       ! Memory allocation
+   USE iom
 
    IMPLICIT NONE
@@ -39,4 +41,7 @@
    INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11   ,   np_blp_i  = 21  ! standard iso-neutral or geopotential operator
    INTEGER, PARAMETER, PUBLIC ::   np_lap_it = 12   ,   np_blp_it = 22  ! triad    iso-neutral or geopotential operator
+
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat transport
 
    !! * Substitutions
@@ -95,4 +100,10 @@
       CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev ) 
       !
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
+      !
       !                                !==  Initialization of metric arrays used for all tracers  ==!
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
@@ -150,8 +161,7 @@
          IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==!
              ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==!
-            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
-               IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( -ztv(:,:,:) )
-               IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( -ztv(:,:,:) )
-            ENDIF
+
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:)  )
+            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) )
          ENDIF
          !                          ! ==================
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90	(revision 7382)
@@ -20,4 +20,5 @@
    USE traldf_iso     ! lateral diffusion (Madec operator)         (tra_ldf_iso routine)
    USE diaptr         ! poleward transport diagnostics
+   USE diaar5         ! AR5 diagnostics
    USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
    !
@@ -35,4 +36,8 @@
 
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels
+
+   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL  ::   l_hst   ! flag to compute heat transport
+
 
    !! * Substitutions
@@ -89,7 +94,4 @@
       REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt
       REAL(wp) ::   zah, zah_slp, zaei_slp
-#if defined key_diaar5
-      REAL(wp) ::   zztmp              ! local scalar
-#endif
       REAL(wp), POINTER, DIMENSION(:,:  ) :: z2d                                            ! 2D workspace
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     -
@@ -112,5 +114,12 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
       ENDIF
-      !                                               ! set time step size (Euler/Leapfrog)
+      !   
+      l_hst = .FALSE.
+      l_ptr = .FALSE.
+      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
+      !
+      !                                                        ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
       ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
@@ -416,38 +425,7 @@
             !
             !                          ! "Poleward" diffusive heat or salt transports (T-S case only)
-            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
-               IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( zftv(:,:,:) )        ! 3.3  names
-               IF( jn == jp_sal)   str_ldf(:) = ptr_sj( zftv(:,:,:) )
-            ENDIF
-            !
-            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
-              !
-              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
-                  z2d(:,:) = zftu(ji,jj,1) 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
-                        END DO
-                     END DO
-                  END DO
-                  z2d(:,:) = rau0_rcp * z2d(:,:) 
-                  CALL lbc_lnk( z2d, 'U', -1. )
-                  CALL iom_put( "udiff_heattr", z2d )                  ! heat i-transport
-                  !
-                  z2d(:,:) = zftv(ji,jj,1) 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
-                        END DO
-                     END DO
-                  END DO
-                  z2d(:,:) = rau0_rcp * z2d(:,:)     
-                  CALL lbc_lnk( z2d, 'V', -1. )
-                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat j-transport
-               ENDIF
-               !
-            ENDIF
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:)  )
+            !                          ! Diffusive heat transports
+            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) )
             !
          ENDIF                                                    !== end pass selection  ==!
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90	(revision 7382)
@@ -33,5 +33,5 @@
 # endif
    !                                                  !!!* Active tracers trends indexes
-   INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 14     !: Total trend nb: change it when adding/removing one indice below
+   INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 20     !: Total trend nb: change it when adding/removing one indice below
    !                               ===============     !  
    INTEGER, PUBLIC, PARAMETER ::   jptra_xad  =  1     !: x- horizontal advection
@@ -39,22 +39,25 @@
    INTEGER, PUBLIC, PARAMETER ::   jptra_zad  =  3     !: z- vertical   advection
    INTEGER, PUBLIC, PARAMETER ::   jptra_sad  =  4     !: z- vertical   advection
-   INTEGER, PUBLIC, PARAMETER ::   jptra_ldf  =  5     !: lateral       diffusion
-   INTEGER, PUBLIC, PARAMETER ::   jptra_zdf  =  6     !: vertical      diffusion
-   INTEGER, PUBLIC, PARAMETER ::   jptra_zdfp =  7     !: "PURE" vert.  diffusion (ln_traldf_iso=T)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_bbc  =  8     !: Bottom Boundary Condition (geoth. heating) 
-   INTEGER, PUBLIC, PARAMETER ::   jptra_bbl  =  9     !: Bottom Boundary Layer (diffusive and/or advective)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_npc  = 10     !: non-penetrative convection treatment
-   INTEGER, PUBLIC, PARAMETER ::   jptra_dmp  = 11     !: internal restoring (damping)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_qsr  = 12     !: penetrative solar radiation
-   INTEGER, PUBLIC, PARAMETER ::   jptra_nsr  = 13     !: non solar radiation / C/D on salinity  (+runoff if ln_rnf=T)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_atf  = 14     !: Asselin time filter
+   INTEGER, PUBLIC, PARAMETER ::   jptra_totad  =  5   !: total         advection
+   INTEGER, PUBLIC, PARAMETER ::   jptra_ldf  =  6     !: lateral       diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jptra_zdf  =  7     !: vertical      diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jptra_zdfp =  8     !: "PURE" vert.  diffusion (ln_traldf_iso=T)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_evd  =  9     !: EVD term (convection)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_bbc  = 10     !: Bottom Boundary Condition (geoth. heating) 
+   INTEGER, PUBLIC, PARAMETER ::   jptra_bbl  = 11     !: Bottom Boundary Layer (diffusive and/or advective)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_npc  = 12     !: non-penetrative convection treatment
+   INTEGER, PUBLIC, PARAMETER ::   jptra_dmp  = 13     !: internal restoring (damping)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_qsr  = 14     !: penetrative solar radiation
+   INTEGER, PUBLIC, PARAMETER ::   jptra_nsr  = 15     !: non solar radiation / C/D on salinity  (+runoff if ln_rnf=T)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_atf  = 16     !: Asselin time filter
+   INTEGER, PUBLIC, PARAMETER ::   jptra_tot  = 17     !: Model total trend
    !
    !                                                  !!!* Passive tracers trends indices (use if "key_top" defined)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_sms  = 15     !: sources m. sinks
-   INTEGER, PUBLIC, PARAMETER ::   jptra_radn = 16     !: corr. trn<0 in trcrad
-   INTEGER, PUBLIC, PARAMETER ::   jptra_radb = 17     !: corr. trb<0 in trcrad (like atf)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_sms  = 18     !: sources m. sinks
+   INTEGER, PUBLIC, PARAMETER ::   jptra_radn = 19     !: corr. trn<0 in trcrad
+   INTEGER, PUBLIC, PARAMETER ::   jptra_radb = 20     !: corr. trb<0 in trcrad (like atf)
    !
    !                                                  !!!* Momentum trends indices
-   INTEGER, PUBLIC, PARAMETER ::   jptot_dyn  = 15     !: Total trend nb: change it when adding/removing one indice below
+   INTEGER, PUBLIC, PARAMETER ::   jptot_dyn  = 13     !: Total trend nb: change it when adding/removing one indice below
    !                               ===============     !  
    INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg  =  1     !: hydrostatic pressure gradient 
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90	(revision 7382)
@@ -90,5 +90,5 @@
 !!gm end
       !
-      IF( .NOT.ln_linssh .AND. ( l_trdtra .OR. l_trddyn ) )  CALL ctl_stop( 'trend diagnostics with variable volume not validated' )
+!      IF( .NOT.ln_linssh .AND. ( l_trdtra .OR. l_trddyn ) )  CALL ctl_stop( 'trend diagnostics with variable volume not validated' )
       
 !!gm  : Potential BUG : 3D output only for vector invariant form!  add a ctl_stop or code the flux form case
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90	(revision 7382)
@@ -28,4 +28,5 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! Memory allocation
+   USE ldfslp         ! Isopycnal slopes
 
    IMPLICIT NONE
@@ -180,11 +181,12 @@
 !                              CALL iom_put( "ketrd_bfri", zke2d )
 !         ENDIF
-         CASE( jpdyn_ken )   ;                                          ! kinetic energy
-                                 ! called in dynnxt.F90 before asselin time filter with putrd=ua and pvtrd=va
-                                 zke(:,:,:) = 0.5_wp * zke(:,:,:)
-                                 CALL iom_put( "KE", zke )
-                                 !
-                                 CALL ken_p2k( kt , zke )
-                                 CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w
+        CASE( jpdyn_ken )   ;   ! kinetic energy
+                    ! called in dynnxt.F90 before asselin time filter
+                    ! with putrd=ua and pvtrd=va
+                    zke(:,:,:) = 0.5_wp * zke(:,:,:)
+                    CALL iom_put( "KE", zke )
+                    !
+                    CALL ken_p2k( kt , zke )
+                      CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w
          !
       END SELECT
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 7382)
@@ -39,4 +39,5 @@
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   trdtx, trdty, trdt   ! use to store the temperature trends
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_evd  ! store avt_evd to calculate EVD trend
 
    !! * Substitutions
@@ -54,5 +55,5 @@
       !!                  ***  FUNCTION trd_tra_alloc  ***
       !!---------------------------------------------------------------------
-      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , STAT= trd_tra_alloc )
+      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , avt_evd(jpi,jpj,jpk), STAT= trd_tra_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( trd_tra_alloc )
@@ -127,5 +128,5 @@
             zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
             DO jk = 2, jpk
-               zwt(:,:,jk) =   avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
+               zwt(:,:,jk) = avt_k(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
                zws(:,:,jk) = fsavs(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
             END DO
@@ -137,4 +138,18 @@
             END DO
             CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt )  
+            !
+            !                         ! Also calculate EVD trend at this point. 
+            zwt(:,:,:) = 0._wp   ;   zws(:,:,:) = 0._wp            ! vertical diffusive fluxes
+            DO jk = 2, jpk
+               zwt(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
+               zws(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
+            END DO
+            !
+            ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
+            DO jk = 1, jpkm1
+               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t_n(:,:,jk)
+               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t_n(:,:,jk) 
+            END DO
+            CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt )  
             !
             CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt )
@@ -311,4 +326,6 @@
                                   CALL wrk_dealloc( jpi, jpj, z2dx, z2dy )
                                ENDIF
+      CASE( jptra_totad  ) ;   CALL iom_put( "ttrd_totad" , ptrdx )      ! total   advection
+                               CALL iom_put( "strd_totad" , ptrdy )
       CASE( jptra_ldf  )   ;   CALL iom_put( "ttrd_ldf" , ptrdx )        ! lateral diffusion
                                CALL iom_put( "strd_ldf" , ptrdy )
@@ -317,4 +334,6 @@
       CASE( jptra_zdfp )   ;   CALL iom_put( "ttrd_zdfp", ptrdx )        ! PURE vertical diffusion (no isoneutral contribution)
                                CALL iom_put( "strd_zdfp", ptrdy )
+      CASE( jptra_evd )    ;   CALL iom_put( "ttrd_evd", ptrdx )         ! EVD trend (convection)
+                               CALL iom_put( "strd_evd", ptrdy )
       CASE( jptra_dmp  )   ;   CALL iom_put( "ttrd_dmp" , ptrdx )        ! internal restoring (damping)
                                CALL iom_put( "strd_dmp" , ptrdy )
@@ -323,10 +342,12 @@
       CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc" , ptrdx )        ! static instability mixing
                                CALL iom_put( "strd_npc" , ptrdy )
-      CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx )        ! surface forcing + runoff (ln_rnf=T)
-                               CALL iom_put( "strd_cdt" , ptrdy )
+      CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx(:,:,1) )        ! surface forcing + runoff (ln_rnf=T)
+                               CALL iom_put( "strd_cdt" , ptrdy(:,:,1) )        ! output as 2D surface fields
       CASE( jptra_qsr  )   ;   CALL iom_put( "ttrd_qsr" , ptrdx )        ! penetrative solar radiat. (only on temperature)
       CASE( jptra_bbc  )   ;   CALL iom_put( "ttrd_bbc" , ptrdx )        ! geothermal heating   (only on temperature)
       CASE( jptra_atf  )   ;   CALL iom_put( "ttrd_atf" , ptrdx )        ! asselin time Filter
                                CALL iom_put( "strd_atf" , ptrdy )
+      CASE( jptra_tot  )   ;   CALL iom_put( "ttrd_tot" , ptrdx )        ! model total trend
+                               CALL iom_put( "strd_tot" , ptrdy )
       END SELECT
       !
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 7382)
@@ -17,4 +17,6 @@
    USE dom_oce         ! ocean space and time domain variables
    USE zdf_oce         ! ocean vertical physics variables
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! trends manager: tracers 
    !
    USE in_out_manager  ! I/O manager
@@ -111,4 +113,5 @@
       zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd
       CALL iom_put( "avt_evd", zavt_evd )              ! output this change
+      IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd )
       !
       CALL wrk_dealloc( jpi,jpj,jpk,   zavt_evd, zavm_evd ) 
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7382)
@@ -490,5 +490,4 @@
       IF( lk_floats     )   CALL     flo_init   ! drifting Floats
                             CALL dia_cfl_init   ! Initialise CFL diagnostics
-      IF( lk_diaar5     )   CALL dia_ar5_init   ! ar5 diag
                             CALL dia_ptr_init   ! Poleward TRansports initialization
       IF( lk_diadct     )   CALL dia_dct_init   ! Sections tranports
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7382)
@@ -234,5 +234,5 @@
       IF(.NOT.ln_cpl )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
       IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
-      IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
+                         CALL dia_ar5( kstp )         ! ar5 diag
       IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
                          CALL dia_wri( kstp )         ! ocean model: outputs
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/SAS_SRC/diawri.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/SAS_SRC/diawri.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/SAS_SRC/diawri.F90	(revision 7382)
@@ -36,5 +36,4 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
-   USE diaar5, ONLY :   lk_diaar5
    USE iom
    USE ioipsl
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 7382)
@@ -114,5 +114,4 @@
    USE zdfmxl , ONLY :   hmlpt       =>   hmlpt       !: mixed layer depth at t-points (m)
 
-   USE diaar5 , ONLY :   lk_diaar5  =>   lk_diaar5
 #else
    !!----------------------------------------------------------------------
