Index: /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6671)
+++ /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6672)
@@ -404,4 +404,5 @@
          <field id="u_masstr_vint" long_name="vertical integral of ocean eulerian mass transport along i-axis"    standard_name="vertical_integral_of_ocean_mass_x_transport"  unit="kg/s" />
          <field id="u_heattr"     long_name="ocean eulerian heat transport along i-axis"    standard_name="ocean_heat_x_transport"                          unit="W"                                />
+         <field id="uadv_heattr"  long_name="ocean advective heat transport along i-axis"    standard_name="advectice_ocean_heat_x_transport"                          unit="W"                                />
          <field id="u_salttr"     long_name="ocean eulerian salt transport along i-axis"    standard_name="ocean_salt_x_transport"                          unit="0.001*kg/s"                        />
          <field id="ueiv_heattr"  long_name="ocean bolus heat transport along i-axis"       standard_name="ocean_heat_x_transport_due_to_bolus_advection"   unit="W"                                />
@@ -441,4 +442,5 @@
          <field id="v_masstr"     long_name="Ocean Mass Y Transport"    standard_name="ocean_mass_y_transport"                          unit="kg/s"        grid_ref="grid_V_3D" />
          <field id="v_heattr"     long_name="Ocean Heat X Transport"    standard_name="ocean_heat_y_transport"                          unit="W"                                />
+         <field id="vadv_heattr"  long_name="ocean advective heat transport along j-axis"    standard_name="advectice_ocean_heat_y_transport"                          unit="W"                                />
          <field id="v_salttr"     long_name="ocean eulerian salt transport along i-axis"    standard_name="ocean_salt_y_transport"                          unit="0.001*kg/s"                        />
          <field id="veiv_heattr"  long_name="ocean bolus heat transport along j-axis"       standard_name="ocean_heat_y_transport_due_to_bolus_advection"   unit="W"                                />
@@ -491,4 +493,23 @@
       </field_group>
           
+      <!-- product fields -->
+      <field_group id="diaprod">
+	<field id="ut"           long_name="product_of_sea_water_x_velocity_and_potential_temperature"      unit="degree_C m/s"      grid_ref="grid_U_3D"   />
+        <field id="ut_e3u"       long_name="product_of_sea_water_x_velocity_and_potential_temperature * e3u"  unit="degree_C m2/s"   grid_ref="grid_U_3D" > ut * e3u </field >
+	<field id="us"           long_name="product_of_sea_water_x_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_U_3D"   />
+        <field id="us_e3u"       long_name="product_of_sea_water_x_velocity_and_salinity * e3u"             unit="PSU m2/s"      grid_ref="grid_U_3D" > us * e3u </field >
+	<field id="urhop"        long_name="product_of_sea_water_x_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_U_3D"   />
+        <field id="urhop_e3u"    long_name="product_of_sea_water_x_velocity_and_potential_density * e3u"    unit="(kg/m3).(m2/s)"   grid_ref="grid_U_3D" > urhop * e3u </field >
+	<field id="vt"           long_name="product_of_sea_water_y_velocity_and_potential_temperature"      unit="degree_C m/s"      grid_ref="grid_V_3D"   />
+        <field id="vt_e3v"       long_name="product_of_sea_water_y_velocity_and_potential_temperature * e3v"  unit="degree_C m2/s"   grid_ref="grid_V_3D" > vt * e3v </field >
+	<field id="vs"           long_name="product_of_sea_water_y_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_V_3D"   />
+        <field id="vs_e3v"       long_name="product_of_sea_water_y_velocity_and_salinity * e3t"             unit="PSU m2/s"      grid_ref="grid_V_3D" > vs * e3v </field >
+	<field id="vrhop"        long_name="product_of_sea_water_y_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_V_3D"   />
+        <field id="vrhop_e3v"    long_name="product_of_sea_water_y_velocity_and_potential_density * e3t"    unit="(kg/m3).(m2/s)"  grid_ref="grid_V_3D" > vrhop * e3v </field >
+	<field id="wt"           long_name="product_of_upward_sea_water_velocity_and_potential_temperature" unit="degree_C m/s"      grid_ref="grid_W_3D"   />
+	<field id="ws"           long_name="product_of_upward_sea_water_velocity_and_salinity"              unit="PSU m/s"       grid_ref="grid_W_3D"   />
+	<field id="wrhop"        long_name="product_of_upward_sea_water_velocity_and_potential_density"     unit="(kg/m3).(m/s)" grid_ref="grid_W_3D"   />
+      </field_group>
+
       <!-- scalar variables available with key_diaar5 -->
 
@@ -659,4 +680,14 @@
         <field id="sopsteiv_ind"      long_name="Salt Transport from mesoscale eddy advection: Indian"             unit="Giga g/s"       grid_ref="gznl_T_2D" />
         <field id="sopsteiv_ipc"      long_name="Salt Transport from mesoscale eddy advection: Pacific+Indian"     unit="Giga g/s"       grid_ref="gznl_T_2D" />       
+        <field id="sopht_vt"          long_name="Heat Transport"                     unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopht_vt_atl"      long_name="Heat Transport: Atlantic"           unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopht_vt_ind"      long_name="Heat Transport: Indian"             unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopht_vt_pac"      long_name="Heat Transport: Pacific"            unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopht_vt_ipc"      long_name="Heat Transport: Indo-Pacific"       unit="PW"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vs"          long_name="Salt Transport"                     unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vt_atl"      long_name="Salt Transport: Atlantic"           unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vt_ind"      long_name="Salt Transport: Indian"             unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vt_pac"      long_name="Salt Transport: Pacific"            unit="Giga g/s"       grid_ref="gznl_T_2D" />
+        <field id="sopst_vt_ipc"      long_name="Salt Transport: Indo-Pacific"       unit="Giga g/s"       grid_ref="gznl_T_2D" />
       </field_group>
 
Index: /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90
===================================================================
--- /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90	(revision 6672)
+++ /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90	(revision 6672)
@@ -0,0 +1,240 @@
+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
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   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 diadimg         ! dimg direct access file format output
+   USE iom
+   USE ioipsl
+   USE lib_mpp         ! MPP library
+   USE timing          ! preformance summary
+   USE wrk_nemo        ! working array
+   USE diaptr
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dia_prod                 ! routines called by step.F90
+
+   !! * Substitutions
+#  include "zdfddm_substitute.h90"
+#  include "domzgr_substitute.h90"
+#  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( nn_timing == 1 )   CALL timing_start('dia_prod')
+      ! 
+      CALL wrk_alloc( jpi , jpj      , z2d )
+      CALL wrk_alloc( jpi , jpj, jpk , z3d )
+      CALL wrk_alloc( jpi , jpj, jpk , zrhop )
+      !
+
+      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") .OR. iom_use("sopht_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
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk) * e1v(ji,jj)
+               END DO
+            END DO
+         END DO
+         CALL dia_ptr_ohst_components( jp_tem, 'vts', z3d)
+      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
+
+      !
+      CALL wrk_dealloc( jpi , jpj      , z2d )
+      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
+      CALL wrk_dealloc( jpi , jpj, jpk , zrhop )
+      !
+      IF( nn_timing == 1 )   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: /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 6671)
+++ /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 6672)
@@ -42,6 +42,6 @@
 
    !                                  !!** namelist  namptr  **
-   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_adv, htr_ldf, htr_eiv, htr_vt   !: Heat TRansports (adv, diff, Bolus.)
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   str_adv, str_ldf, str_eiv, str_vs   !: Salt TRansports (adv, diff, Bolus.)
 
    LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F)
@@ -207,4 +207,35 @@
                ENDDO
                cl1 = TRIM('sopstldf_'//clsubb(jn))                 
+               CALL iom_put( cl1, z2d )              
+              ENDDO
+            ENDIF
+         ENDIF
+
+         IF( iom_use("sopht_vt") .OR. iom_use("sopst_vs") ) THEN   
+            z2d(1,:) = htr_vt(:,1) * rc_pwatt        !  (conversion in PW) 
+            DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sopht_vt'
+            CALL iom_put( TRIM(cl1), z2d )
+            z2d(1,:) = str_vs(:,1) * rc_ggram        !  (conversion in Gg)
+            DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sopst_vs'
+            CALL iom_put( TRIM(cl1), z2d )
+            IF( ln_subbas ) THEN
+              DO jn=2,nptr
+               z2d(1,:) = htr_vt(:,jn) * rc_pwatt        !  (conversion in PW)
+               DO ji = 1, jpi
+                 z2d(ji,:) = z2d(1,:)
+               ENDDO
+               cl1 = TRIM('sopht_vt_'//clsubb(jn))                 
+               CALL iom_put( cl1, z2d )
+               z2d(1,:) = str_vs(:,jn) * rc_ggram        ! (conversion in Gg)
+               DO ji = 1, jpi
+                  z2d(ji,:) = z2d(1,:)
+               ENDDO
+               cl1 = TRIM('sopst_vs_'//clsubb(jn))                 
                CALL iom_put( cl1, z2d )              
               ENDDO
@@ -327,4 +358,5 @@
          htr_ldf(:,:) = 0._wp  ;  str_ldf(:,:) =  0._wp 
          htr_eiv(:,:) = 0._wp  ;  str_eiv(:,:) =  0._wp 
+         htr_vt(:,:) = 0._wp  ;   str_vs(:,:) =  0._wp 
          !
       ENDIF 
@@ -357,4 +389,8 @@
          IF( ktra == jp_sal )  str_eiv(:,1) = ptr_sj( pva(:,:,:) )
       ENDIF
+      IF( cptr == 'vts' ) THEN
+         IF( ktra == jp_tem )  htr_vt(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_sal )  str_vs(:,1) = ptr_sj( pva(:,:,:) )
+      ENDIF
       !
       IF( ln_subbas ) THEN
@@ -393,4 +429,16 @@
                 DO jn = 2, nptr
                    str_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                END DO
+             ENDIF
+         ENDIF
+         IF( cptr == 'vts' ) THEN
+             IF( ktra == jp_tem ) THEN 
+                DO jn = 2, nptr
+                    htr_vt(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                 END DO
+             ENDIF
+             IF( ktra == jp_sal ) THEN 
+                DO jn = 2, nptr
+                   str_vs(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
                 END DO
              ENDIF
@@ -414,4 +462,5 @@
          &      htr_adv(jpj,nptr) , str_adv(jpj,nptr) ,   &
          &      htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) ,   &
+         &      htr_vt(jpj,nptr)  , str_vs(jpj,nptr)  ,   &
          &      htr_ldf(jpj,nptr) , str_ldf(jpj,nptr) , STAT=ierr(1)  )
          !
Index: /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6671)
+++ /branches/UKMO/v3_6_extra_CMIP6_diagnostics/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6672)
@@ -234,4 +234,5 @@
       IF( lk_diaar5  )      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
       !
