Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/context_nemo-BASIC.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/context_nemo-BASIC.xml	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/context_nemo-BASIC.xml	(revision 11110)
@@ -35,4 +35,7 @@
       <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"            />
       <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"            />
+      <axis id="deptht_surface" axis_ref="deptht" >
+         <zoom_axis begin=" 0 " n=" 1 " />
+      </axis>
     </axis_definition>
  
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/context_nemo-CMIP6.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/context_nemo-CMIP6.xml	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/context_nemo-CMIP6.xml	(revision 11110)
@@ -35,4 +35,7 @@
       <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"            />
       <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"            />
+      <axis id="deptht_surface" axis_ref="deptht" >
+         <zoom_axis begin=" 0 " n=" 1 " />
+      </axis>
     </axis_definition>
  
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/context_nemo-CUSTOM.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/context_nemo-CUSTOM.xml	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/context_nemo-CUSTOM.xml	(revision 11110)
@@ -35,4 +35,7 @@
       <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"            />
       <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"            />
+      <axis id="deptht_surface" axis_ref="deptht" >
+         <zoom_axis begin=" 0 " n=" 1 " />
+      </axis>
     </axis_definition>
  
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-oce-CMIP6-10D.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-oce-CMIP6-10D.xml	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-oce-CMIP6-10D.xml	(revision 11110)
@@ -153,8 +153,6 @@
           <field field_ref="uoce"         name="uo"     operation="average" freq_op="10d" > @uoce_e3u / @e3u </field>
 	  <field field_ref="utau"         name="tauuo"   standard_name="surface_downward_x_stress"  />
-<!-- TO FIX
           <field field_ref="ut"           name="uto" standard_name="product_of_xward_sea_water_velocity_and_temperature" long_name="UT"  operation="average" freq_op="10d" > @ut_e3u / @e3u </field>
           <field field_ref="us"           name="uso" standard_name="product_of_xward_sea_water_velocity_and_salinity" long_name="US"    operation="average" freq_op="10d" > @us_e3u / @e3u </field>
--->
           <field field_ref="uoce"         name="u2o" standard_name="square_of_sea_water_x_velocity" long_name="UU"  operation="average" freq_op="10d" > @uoce2_e3u / @e3u </field>
           <!-- available with key_diaar5 -->
@@ -172,8 +170,6 @@
           <field field_ref="voce"         name="vo"     operation="average" freq_op="10d" > @voce_e3v / @e3v </field>
 	  <field field_ref="vtau"         name="tauvo"   standard_name="surface_downward_y_stress" />
-<!-- TO FIX
           <field field_ref="vt"           name="vto" standard_name="product_of_yward_sea_water_velocity_and_temperature" long_name="VT"  operation="average" freq_op="10d" > @vt_e3v / @e3v </field>
           <field field_ref="vs"           name="vso" standard_name="product_of_yward_sea_water_velocity_and_salinity" long_name="VS"  operation="average" freq_op="10d" > @vs_e3v / @e3v </field>
--->
           <field field_ref="voce"         name="v2o" standard_name="square_of_sea_water_y_velocity" long_name="VV"  operation="average" freq_op="10d" > @voce2_e3v / @e3v </field>
           <!-- available with key_diaar5 -->
@@ -195,7 +191,7 @@
 	  <field field_ref="woce"         name="wo" standard_name="upward_sea_water_velocity" long_name="W"                  />
           <field field_ref="woce"         name="w2o" standard_name="square_of_upward_sea_water_velocity" long_name="WW"  operation="average" > woce * woce </field>
-<!-- TO FIX
           <field field_ref="wt"           name="wt"         standard_name="product_of_upward_sea_water_velocity_and_temperature" long_name="WT"         />
           <field field_ref="ws"           name="ws"         standard_name="product_of_upward_sea_water_velocity_and_salinity" long_name="WS"            />
+<!-- TO FIX
           <field field_ref="wpt_dep"      name="zhalf"        />
 -->
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-oce-CMIP6-1M.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-oce-CMIP6-1M.xml	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-oce-CMIP6-1M.xml	(revision 11110)
@@ -153,8 +153,6 @@
           <field field_ref="uoce"         name="uo"     operation="average" freq_op="1mo" > @uoce_e3u / @e3u </field>
 	  <field field_ref="utau"         name="tauuo"   standard_name="surface_downward_x_stress"  />
-<!-- TO FIX
           <field field_ref="ut"           name="uto" standard_name="product_of_xward_sea_water_velocity_and_temperature" long_name="UT"  operation="average" freq_op="1mo" > @ut_e3u / @e3u </field>
           <field field_ref="us"           name="uso" standard_name="product_of_xward_sea_water_velocity_and_salinity" long_name="US"    operation="average" freq_op="1mo" > @us_e3u / @e3u </field>
--->
           <field field_ref="uoce"         name="u2o" standard_name="square_of_sea_water_x_velocity" long_name="UU"  operation="average" freq_op="1mo" > @uoce2_e3u / @e3u </field>
           <!-- available with key_diaar5 -->
@@ -172,8 +170,6 @@
           <field field_ref="voce"         name="vo"     operation="average" freq_op="1mo" > @voce_e3v / @e3v </field>
 	  <field field_ref="vtau"         name="tauvo"   standard_name="surface_downward_y_stress" />
-<!-- TO FIX
           <field field_ref="vt"           name="vto" standard_name="product_of_yward_sea_water_velocity_and_temperature" long_name="VT"  operation="average" freq_op="1mo" > @vt_e3v / @e3v </field>
           <field field_ref="vs"           name="vso" standard_name="product_of_yward_sea_water_velocity_and_salinity" long_name="VS"  operation="average" freq_op="1mo" > @vs_e3v / @e3v </field>
--->
           <field field_ref="voce"         name="v2o" standard_name="square_of_sea_water_y_velocity" long_name="VV"  operation="average" freq_op="1mo" > @voce2_e3v / @e3v </field>
           <!-- available with key_diaar5 -->
@@ -195,7 +191,7 @@
 	  <field field_ref="woce"         name="wo" standard_name="upward_sea_water_velocity" long_name="W"                  />
           <field field_ref="woce"         name="w2o" standard_name="square_of_upward_sea_water_velocity" long_name="WW"  operation="average" > woce * woce </field>
-<!-- TO FIX
           <field field_ref="wt"           name="wt"         standard_name="product_of_upward_sea_water_velocity_and_temperature" long_name="WT"         />
           <field field_ref="ws"           name="ws"         standard_name="product_of_upward_sea_water_velocity_and_salinity" long_name="WS"            />
+<!-- TO FIX
           <field field_ref="wpt_dep"      name="zhalf"        />
 -->
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/domain_def_nemo.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/domain_def_nemo.xml	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/domain_def_nemo.xml	(revision 11110)
@@ -182,4 +182,7 @@
      <!--   zonal mean grid   -->
      <domain id="gznl" long_name="gznl"/>
+     <domain id="ptr" domain_ref="gznl" >
+	 <zoom_domain id="ptr" ibegin="0000" jbegin="0" ni="1" nj="0000" />
+     </domain>  
      <domain id="znl_T" domain_ref="gznl" > <zoom_domain id="znl_T"/> </domain>
      <domain id="znl_W" domain_ref="gznl" > <zoom_domain id="znl_W"/> </domain>
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11110)
@@ -20,5 +20,5 @@
       <field_group id="grid_T" grid_ref="grid_T_2D" >
         <field id="e3t"          long_name="T-cell thickness"                    standard_name="cell_thickness"     unit="m"   grid_ref="grid_T_3D" />
-        <field id="e3t_surf"     long_name="T-cell thickness"   field_ref="e3t"  standard_name="cell_thickness"     unit="m"   grid_ref="grid_T_SFC"/>
+        <field id="e3t_surf"     long_name="T-cell thickness"   field_ref="e3t"  standard_name="cell_thickness"     unit="m"   grid_ref="grid_T_surface_extract"/>
         <field id="e3t_0"        long_name="Initial T-cell thickness"            standard_name="ref_cell_thickness" unit="m"   grid_ref="grid_T_3D" />
 
@@ -511,4 +511,23 @@
       <field id="ahmf_2d"      long_name=" surface f-eddy viscosity coefficient"   unit="m2/s or m4/s" />
       <field id="ahmf_3d"      long_name=" 3D      f-eddy viscosity coefficient"   unit="m2/s or m4/s"  grid_ref="grid_T_3D"/>
+
+      <!-- product fields -->
+      <field_group id="diaprod">
+	<field id="ut"           long_name="product_of_sea_water_x_velocity_and_potential_temperature"      unit="degree_C m/s"      grid_ref="grid_U_3D"   />
+        <field id="ut_e3u"       long_name="product_of_sea_water_x_velocity_and_potential_temperature * e3u"  unit="degree_C m2/s"   grid_ref="grid_U_3D" > ut * e3u </field >
+	<field id="us"           long_name="product_of_sea_water_x_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_U_3D"   />
+        <field id="us_e3u"       long_name="product_of_sea_water_x_velocity_and_salinity * e3u"             unit="PSU m2/s"      grid_ref="grid_U_3D" > us * e3u </field >
+	<field id="urhop"        long_name="product_of_sea_water_x_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_U_3D"   />
+        <field id="urhop_e3u"    long_name="product_of_sea_water_x_velocity_and_potential_density * e3u"    unit="(kg/m3).(m2/s)"   grid_ref="grid_U_3D" > urhop * e3u </field >
+	<field id="vt"           long_name="product_of_sea_water_y_velocity_and_potential_temperature"      unit="degree_C m/s"      grid_ref="grid_V_3D"   />
+        <field id="vt_e3v"       long_name="product_of_sea_water_y_velocity_and_potential_temperature * e3v"  unit="degree_C m2/s"   grid_ref="grid_V_3D" > vt * e3v </field >
+	<field id="vs"           long_name="product_of_sea_water_y_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_V_3D"   />
+        <field id="vs_e3v"       long_name="product_of_sea_water_y_velocity_and_salinity * e3t"             unit="PSU m2/s"      grid_ref="grid_V_3D" > vs * e3v </field >
+	<field id="vrhop"        long_name="product_of_sea_water_y_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_V_3D"   />
+        <field id="vrhop_e3v"    long_name="product_of_sea_water_y_velocity_and_potential_density * e3t"    unit="(kg/m3).(m2/s)"  grid_ref="grid_V_3D" > vrhop * e3v </field >
+	<field id="wt"           long_name="product_of_upward_sea_water_velocity_and_potential_temperature" unit="degree_C m/s"      grid_ref="grid_W_3D"   />
+	<field id="ws"           long_name="product_of_upward_sea_water_velocity_and_salinity"              unit="PSU m/s"       grid_ref="grid_W_3D"   />
+	<field id="wrhop"        long_name="product_of_upward_sea_water_velocity_and_potential_density"     unit="(kg/m3).(m/s)" grid_ref="grid_W_3D"   />
+      </field_group>
 
       <field_group id="scalar"  grid_ref="grid_T_2D"  >
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/grid_def_nemo.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/grid_def_nemo.xml	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/cfgs/SHARED/grid_def_nemo.xml	(revision 11110)
@@ -57,4 +57,46 @@
          <axis id="nfloat" />
        </grid>
+        <!--  -->
+       <grid id="grid_EqT" >
+         <domain id="EqT" />
+       </grid>
+        <!--  -->
+       <grid id="gznl_T_2D">
+         <domain id="ptr" />
+       </grid>
+        <!--  -->
+       <grid id="gznl_T_3D">
+         <domain id="ptr" />
+         <axis axis_ref="deptht" />
+       </grid>
+        <!--  -->
+       <grid id="gznl_W_2D">
+         <domain id="ptr" />
+       </grid>
+        <!--  -->
+       <grid id="gznl_W_3D">
+         <domain id="ptr" />
+         <axis axis_ref="depthw" />
+       </grid>
+       <grid id="vert_sum">
+         <domain id="grid_T"/>
+         <scalar>
+            <reduce_axis operation="sum" />
+         </scalar>
+       </grid>
+       <grid id="zoom_300">
+         <domain id="grid_T" />
+         <axis axis_ref="deptht300"/>
+       </grid>
+       <grid id="zoom_300_sum">
+         <domain id="grid_T" />
+         <scalar>
+            <reduce_axis operation="sum" />
+         </scalar>
+       </grid>
+       <grid id="grid_T_surface_extract">
+         <domain id="grid_T" />
+         <axis   axis_ref="deptht_surface" />
+       </grid>
 
     </grid_definition>
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DIA/diaprod.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DIA/diaprod.F90	(revision 11110)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/DIA/diaprod.F90	(revision 11110)
@@ -0,0 +1,224 @@
+MODULE diaprod
+! Requires key_iom_put
+# if defined key_iomput
+   !!======================================================================
+   !!                     ***  MODULE  diaprod  ***
+   !! Ocean diagnostics :  write ocean product diagnostics
+   !!=====================================================================
+   !! History :  3.4  ! 2012  (D. Storkey)  Original code
+   !!            4.0  ! 2019  (D. Storkey)
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dia_prod      : calculate and write out product diagnostics
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE domvvl          ! for thickness weighted diagnostics if key_vvl
+   USE eosbn2          ! equation of state  (eos call)
+   USE phycst          ! physical constants
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager  ! I/O manager
+   USE iom
+   USE ioipsl
+   USE lib_mpp         ! MPP library
+   USE timing          ! preformance summary
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dia_prod                 ! routines called by step.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.4 , NEMO Consortium (2012)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_prod( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_prod  ***
+      !!                   
+      !! ** Purpose :   Write out product diagnostics (uT, vS etc.)
+      !!
+      !! ** Method  :  use iom_put
+      !!               Product diagnostics are not thickness-weighted in this routine.
+      !!               They should be thickness-weighted using XIOS if key_vvl is set. 
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !!
+      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
+      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   ! 
+      !!
+      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhop    ! potential density
+      !!----------------------------------------------------------------------
+      ! 
+      IF( ln_timing )   CALL timing_start('dia_prod')
+      ! 
+      ALLOCATE( z2d(jpi,jpj), z3d(jpi,jpj,jpk), zrhop(jpi,jpj,jpk) )
+      !
+
+      IF( iom_use("urhop") .OR. iom_use("vrhop") .OR. iom_use("wrhop") &
+#if ! defined key_diaar5
+     &  .OR. iom_use("rhop") &
+#endif
+     & ) THEN
+         CALL eos( tsn, z3d, zrhop )                 ! now in situ and potential density
+         zrhop(:,:,:) = zrhop(:,:,:)-1000.e0         ! reference potential density to 1000 to avoid precision issues in rhop2 calculation
+         zrhop(:,:,jpk) = 0._wp
+#if ! defined key_diaar5
+         CALL iom_put( 'rhop', zrhop )
+#else
+         ! If key_diaar5 set then there is already an iom_put call to output rhop.
+         ! Really should be a standard diagnostics option?
+#endif
+      ENDIF
+
+      IF( iom_use("ut") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "ut", z3d )                  ! product of temperature and zonal velocity at U points
+      ENDIF
+
+      IF( iom_use("vt") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "vt", z3d )                  ! product of temperature and meridional velocity at V points
+      ENDIF
+
+      IF( iom_use("wt") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,1) = wn(ji,jj,1) * tsn(ji,jj,1,jp_tem)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk-1,jp_tem) + tsn(ji,jj,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "wt", z3d )                  ! product of temperature and vertical velocity at W points
+      ENDIF
+
+      IF( iom_use("us") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "us", z3d )                  ! product of salinity and zonal velocity at U points
+      ENDIF
+
+      IF( iom_use("vs") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "vs", z3d )                  ! product of salinity and meridional velocity at V points
+      ENDIF
+
+      IF( iom_use("ws") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,1) = wn(ji,jj,1) * tsn(ji,jj,1,jp_sal)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk-1,jp_sal) + tsn(ji,jj,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "ws", z3d )                  ! product of salinity and vertical velocity at W points
+      ENDIF
+
+      IF( iom_use("urhop") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk) + zrhop(ji+1,jj,jk) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "urhop", z3d )                  ! product of density and zonal velocity at U points
+      ENDIF
+
+      IF( iom_use("vrhop") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk) + zrhop(ji,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "vrhop", z3d )                  ! product of density and meridional velocity at V points
+      ENDIF
+
+      IF( iom_use("wrhop") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,1) = wn(ji,jj,1) * zrhop(ji,jj,1)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk-1) + zrhop(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "wrhop", z3d )                  ! product of density and vertical velocity at W points
+      ENDIF
+
+      !
+      DEALLOCATE( z2d, z3d, zrhop )
+      !
+      IF( ln_timing )   CALL timing_stop('dia_prod')
+      !
+   END SUBROUTINE dia_prod
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :                                         NO diaprod
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER :: lk_diaprod = .FALSE.   ! coupled flag
+CONTAINS
+   SUBROUTINE dia_prod( kt )   ! Empty routine
+      INTEGER ::   kt
+      WRITE(*,*) 'dia_prod: You should not have seen this print! error?', kt
+   END SUBROUTINE dia_prod
+#endif
+   !!======================================================================
+END MODULE diaprod
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/step.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/step.F90	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/step.F90	(revision 11110)
@@ -218,4 +218,5 @@
                          CALL dia_ar5 ( kstp )        ! ar5 diag
       IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
+                         CALL dia_prod( kstp )        ! ocean model: product diagnostics
                          CALL dia_wri ( kstp )        ! ocean model: outputs
       !
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/step_oce.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/step_oce.F90	(revision 11109)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/OCE/step_oce.F90	(revision 11110)
@@ -80,4 +80,5 @@
    USE diahsb          ! heat, salt and volume budgets    (dia_hsb routine)
    USE diaharm
+   USE diaprod
    USE diacfl
    USE diaobs          ! Observation operator
