Index: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_momentum_trends/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11612)
+++ NEMO/branches/UKMO/NEMO_4.0_momentum_trends/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11613)
@@ -353,4 +353,5 @@
       <field_group id="grid_U"   grid_ref="grid_U_2D">
         <field id="e3u"          long_name="U-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_U_3D" />
+        <field id="e3u_surf"     long_name="U-cell thickness"   field_ref="e3u"  standard_name="cell_thickness"     unit="m"   grid_ref="grid_U_surface_extract"/>
         <field id="e3u_0"        long_name="Initial U-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_U_3D"/>
         <field id="utau"         long_name="Wind Stress along i-axis"                               standard_name="surface_downward_x_stress"   unit="N/m2"                            />
@@ -407,4 +408,5 @@
       <field_group id="grid_V"   grid_ref="grid_V_2D">
         <field id="e3v"          long_name="V-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_V_3D" />
+        <field id="e3v_surf"     long_name="V-cell thickness"   field_ref="e3v"  standard_name="cell_thickness"     unit="m"   grid_ref="grid_V_surface_extract"/>
         <field id="e3v_0"        long_name="Initial V-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_V_3D" />
         <field id="vtau"         long_name="Wind Stress along j-axis"                               standard_name="surface_downward_y_stress"   unit="N/m2"                            />
@@ -809,10 +811,12 @@
      <!-- variables available with ln_dyn_trd -->
      <field id="utrd_hpg"       long_name="i-trend: hydrostatic pressure gradient"          unit="m/s^2"                        />
-     <field id="utrd_spg"       long_name="i-trend: surface     pressure gradient"          unit="m/s^2"                        />
+     <field id="utrd_hpg_corr"  long_name="i-trend: HPG: correction"       unit="m/s^2"    grid_ref="grid_U_2D" />
+     <field id="utrd_spg2d"     long_name="i-trend: surface pressure gradient: true trend"       unit="m/s^2"    grid_ref="grid_U_2D" />
      <field id="utrd_spgexp"    long_name="i-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        />
-     <field id="utrd_spgflt"    long_name="i-trend: surface pressure gradient (filtered)"   unit="m/s^2"                        />
      <field id="utrd_keg"       long_name="i-trend: KE gradient         or hor. adv."       unit="m/s^2"                        />
      <field id="utrd_rvo"       long_name="i-trend: relative  vorticity or metric term"     unit="m/s^2"                        />
-     <field id="utrd_pvo"       long_name="i-trend: planetary vorticity"                    unit="m/s^2"                        />
+     <field id="utrd_pvo"       long_name="i-trend: planetary vorticity: 3D component"      unit="m/s^2"                        />
+     <field id="utrd_pvo2d"     long_name="i-trend: planetary vorticity: 2D component"      unit="m/s^2"    grid_ref="grid_U_2D" />
+     <field id="utrd_pvo_corr"     long_name="i-trend: planetary vorticity: correction"      unit="m/s^2"    grid_ref="grid_U_2D" />
      <field id="utrd_zad"       long_name="i-trend: vertical  advection"                    unit="m/s^2"                        />
      <field id="utrd_udx"       long_name="i-trend: U.dx[U]"                                unit="m/s^2"                        />
@@ -822,6 +826,23 @@
      <field id="utrd_bfr"       long_name="i-trend: bottom friction (explicit)"             unit="m/s^2"                        />   
      <field id="utrd_bfri"      long_name="i-trend: bottom friction (implicit)"             unit="m/s^2"                        />   
+     <field id="utrd_bfr2d"     long_name="i-trend: bottom friction: 2D component"      unit="m/s^2"    grid_ref="grid_U_2D" />
      <field id="utrd_tot"       long_name="i-trend: total momentum trend before atf"        unit="m/s^2"                        />   
      <field id="utrd_atf"       long_name="i-trend: asselin time filter trend"              unit="m/s^2"                        />   
+     <!-- thickness weighted versions -->
+     <field id="utrd_hpg_e3u"       unit="m2/s^2" > utrd_hpg * e3u </field>
+     <field id="utrd_spg_e3u"       unit="m2/s^2" > utrd_spg * e3u </field>
+     <field id="utrd_spgexp_e3u"    unit="m2/s^2" > utrd_spgexp * e3u </field>
+     <field id="utrd_keg_e3u"       unit="m2/s^2" > utrd_keg * e3u </field>
+     <field id="utrd_rvo_e3u"       unit="m2/s^2" > utrd_rvo * e3u </field>
+     <field id="utrd_pvo_e3u"       unit="m2/s^2" > utrd_pvo * e3u </field>
+     <field id="utrd_zad_e3u"       unit="m2/s^2" > utrd_zad * e3u </field>
+     <field id="utrd_ldf_e3u"       unit="m2/s^2" > utrd_ldf * e3u </field>
+     <field id="utrd_zdf_e3u"       unit="m2/s^2" > utrd_zdf * e3u </field>
+     <field id="utrd_tot_e3u"       unit="m2/s^2" > utrd_tot * e3u </field>
+     <field id="utrd_atf_e3u"       unit="m2/s^2" > utrd_atf * e3u </field>
+     <field id="utrd_bta_e3u"       unit="m2/s^2" > utrd_bta * e3u </field>  
+     <field id="utrd_bfr_e3u"       unit="m2/s^2" > utrd_bfr * e3u </field>  
+     <field id="utrd_bfri_e3u"      unit="m2/s^2" > utrd_bfri * e3u </field>  
+     <field id="utrd_tau_e3u"       unit="m2/s^2" > utrd_tau * e3u_surf </field>  
    </field_group>
 
@@ -829,12 +850,16 @@
      <!-- variables available with ln_dyn_trd -->
      <field id="vtrd_hpg"       long_name="j-trend: hydrostatic pressure gradient"          unit="m/s^2"                        />
+     <field id="vtrd_hpg_corr"  long_name="j-trend: HPG: correction"       unit="m/s^2"    grid_ref="grid_V_2D" />
      <field id="vtrd_spg"       long_name="j-trend: surface     pressure gradient"          unit="m/s^2"                        />
+     <field id="vtrd_spgbt"     long_name="j-trend: surface pressure gradient: barotropic part"  unit="m/s^2"    grid_ref="grid_V_2D" />
+     <field id="vtrd_spg2d"     long_name="j-trend: surface pressure gradient: true trend"       unit="m/s^2"    grid_ref="grid_V_2D" />
      <field id="vtrd_spgexp"    long_name="j-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        />
-     <field id="vtrd_spgflt"    long_name="j-trend: surface pressure gradient (filtered)"   unit="m/s^2"                        />
      <field id="vtrd_keg"       long_name="j-trend: KE gradient         or hor. adv."       unit="m/s^2"                        />
      <field id="vtrd_rvo"       long_name="j-trend: relative  vorticity or metric term"     unit="m/s^2"                        />
-     <field id="vtrd_pvo"       long_name="j-trend: planetary vorticity"                    unit="m/s^2"                        />
+     <field id="vtrd_pvo"       long_name="j-trend: planetary vorticity: 3D component"      unit="m/s^2"                        />
+     <field id="vtrd_pvo2d"     long_name="j-trend: planetary vorticity: 2D component"      unit="m/s^2"    grid_ref="grid_V_2D" />
+     <field id="vtrd_pvo_corr"     long_name="j-trend: planetary vorticity: correction"      unit="m/s^2"    grid_ref="grid_V_2D" />
      <field id="vtrd_zad"       long_name="j-trend: vertical  advection"                    unit="m/s^2"                        />
-     <field id="vtrd_vdy"       long_name="i-trend: V.dx[V]"                                unit="m/s^2"                        />
+     <field id="vtrd_vdy"       long_name="j-trend: V.dx[V]"                                unit="m/s^2"                        />
      <field id="vtrd_ldf"       long_name="j-trend: lateral   diffusion"                    unit="m/s^2"                        />
      <field id="vtrd_zdf"       long_name="j-trend: vertical  diffusion"                    unit="m/s^2"                        />
@@ -842,6 +867,23 @@
      <field id="vtrd_bfr"       long_name="j-trend: bottom friction (explicit)"             unit="m/s^2"                        />   
      <field id="vtrd_bfri"      long_name="j-trend: bottom friction (implicit)"             unit="m/s^2"                        />   
+     <field id="vtrd_bfr2d"     long_name="j-trend: bottom friction: 2D component"      unit="m/s^2"    grid_ref="grid_V_2D" />
      <field id="vtrd_tot"       long_name="j-trend: total momentum trend before atf"        unit="m/s^2"                        />   
      <field id="vtrd_atf"       long_name="j-trend: asselin time filter trend"              unit="m/s^2"                        />   
+     <!-- thickness weighted versions -->
+     <field id="vtrd_hpg_e3v"       unit="m2/s^2" > vtrd_hpg * e3v </field>
+     <field id="vtrd_spg_e3v"       unit="m2/s^2" > vtrd_spg * e3v </field>
+     <field id="vtrd_spgexp_e3v"    unit="m2/s^2" > vtrd_spgexp * e3v </field>
+     <field id="vtrd_keg_e3v"       unit="m2/s^2" > vtrd_keg * e3v </field>
+     <field id="vtrd_rvo_e3v"       unit="m2/s^2" > vtrd_rvo * e3v </field>
+     <field id="vtrd_pvo_e3v"       unit="m2/s^2" > vtrd_pvo * e3v </field>
+     <field id="vtrd_zad_e3v"       unit="m2/s^2" > vtrd_zad * e3v </field>
+     <field id="vtrd_ldf_e3v"       unit="m2/s^2" > vtrd_ldf * e3v </field>
+     <field id="vtrd_zdf_e3v"       unit="m2/s^2" > vtrd_zdf * e3v </field>
+     <field id="vtrd_tot_e3v"       unit="m2/s^2" > vtrd_tot * e3v </field>
+     <field id="vtrd_atf_e3v"       unit="m2/s^2" > vtrd_atf * e3v </field>
+     <field id="vtrd_bta_e3v"       unit="m2/s^2" > vtrd_bta * e3v </field>  
+     <field id="vtrd_bfr_e3v"       unit="m2/s^2" > vtrd_bfr * e3v </field>  
+     <field id="vtrd_bfri_e3v"      unit="m2/s^2" > vtrd_bfri * e3v </field>  
+     <field id="vtrd_tau_e3v"       unit="m2/s^2" > vtrd_tau * e3v_surf </field>  
    </field_group>
 
Index: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynhpg.F90	(revision 11612)
+++ NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynhpg.F90	(revision 11613)
@@ -115,5 +115,5 @@
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt )
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg_save, kt )
          DEALLOCATE( ztrdu , ztrdv )
       ENDIF
Index: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynspg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynspg.F90	(revision 11612)
+++ NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynspg.F90	(revision 11613)
@@ -81,5 +81,5 @@
       IF( ln_timing )   CALL timing_start('dyn_spg')
       !
-      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
+      IF( l_trddyn .AND. nspg == np_EXP )   THEN                      ! temporary save of ta and sa trends
          ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
          ztrdu(:,:,:) = ua(:,:,:)
@@ -165,8 +165,12 @@
       END SELECT
       !                    
-      IF( l_trddyn )   THEN                  ! save the surface pressure gradient trends for further diagnostics
+      IF( l_trddyn .AND. nspg == np_EXP )   THEN     ! save the surface pressure gradient trends for further diagnostics
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt )
+         ! In this case we also need to finalise the write-out of the PVO 3D diagnostic with zero correction
+         ztrdu(:,:,1) = 0._wp
+         ztrdv(:,:,1) = 0._wp
+         CALL trd_dyn( ztrdu(:,:,1), ztrdv(:,:,1), jpdyn_pvo, kt )
          DEALLOCATE( ztrdu , ztrdv ) 
       ENDIF
Index: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynspg_ts.F90	(revision 11612)
+++ NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynspg_ts.F90	(revision 11613)
@@ -63,4 +63,6 @@
    USE restart         ! only for lrst_oce
    USE diatmb          ! Top,middle,bottom output
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
 
    IMPLICIT NONE
@@ -172,4 +174,6 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zspgtrdu, zspgtrdv, zpvotrdu, zpvotrdv  ! SPG and PVO trends (if l_trddyn)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zbfrtrdu, zbfrtrdv                      ! BFR trends (if l_trddyn)
       !!----------------------------------------------------------------------
       !
@@ -177,4 +181,14 @@
       !                                         !* Allocate temporary arrays
       IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
+      !
+      IF( l_trddyn ) THEN
+          ALLOCATE( zspgtrdu(jpi,jpj), zspgtrdv(jpi,jpj), zpvotrdu(jpi,jpj), zpvotrdv(jpi,jpj), zbfrtrdu(jpi,jpj), zbfrtrdv(jpi,jpj) )
+          zspgtrdu(:,:) = 0._wp
+          zspgtrdv(:,:) = 0._wp
+          zpvotrdu(:,:) = 0._wp
+          zpvotrdv(:,:) = 0._wp
+          zbfrtrdu(:,:) = 0._wp
+          zbfrtrdv(:,:) = 0._wp
+      ENDIF
       !
       zmdi=1.e+20                               !  missing data indicator for masking
@@ -381,5 +395,5 @@
 !!gm  Is it correct to do so ?   I think so...
       
-      
+
       !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
       !                                   ! --------------------------------------------------------
@@ -387,4 +401,7 @@
       zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes 
       zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
+      !
+      zu_trd(:,:) = 0._wp
+      zv_trd(:,:) = 0._wp
       !
       SELECT CASE( nvor_scheme )
@@ -393,12 +410,12 @@
             DO ji = 2, jpim1   ! vector opt.
                zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    &
-                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   &
-                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   )
-                  !
+                    &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   &
+                    &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   )
+               !
                zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    &
-                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   & 
-                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   ) 
-            END DO  
-         END DO  
+                    &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   & 
+                    &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   ) 
+            END DO
+         END DO
          !         
       CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
@@ -443,7 +460,15 @@
       END SELECT
       !
+      IF( l_trddyn ) THEN
+         ! send correction to baroclinic planetary vorticity trend to trd_dyn
+         CALL trd_dyn( zu_trd, zv_trd, jpdyn_pvo_corr, kt )
+      ENDIF
       !                                   !* Right-Hand-Side of the barotropic momentum equation
       !                                   ! ----------------------------------------------------
       IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient
+         IF( l_trddyn ) THEN
+            zspgtrdu(:,:) = zu_trd(:,:)
+            zspgtrdv(:,:) = zv_trd(:,:)
+         ENDIF
          IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters
             DO jj = 2, jpjm1
@@ -505,4 +530,14 @@
                END DO
             END DO
+         ENDIF
+         !
+         IF( l_trddyn ) THEN
+            zspgtrdu(:,:) = zu_trd(:,:) - zspgtrdu(:,:) 
+            zspgtrdv(:,:) = zv_trd(:,:) - zspgtrdv(:,:) 
+            ! send correction to HPG trend to trd_dyn
+            CALL trd_dyn( zspgtrdu, zspgtrdv, jpdyn_hpg_corr, kt )
+            ! reset temporary arrays for use later
+            zspgtrdu(:,:) = 0._wp
+            zspgtrdv(:,:) = 0._wp
          ENDIF
          !
@@ -1056,4 +1091,10 @@
          END SELECT
          !
+         IF( l_trddyn ) THEN
+            za2 = wgtbtp2(jn)
+            zpvotrdu(:,:) = zpvotrdu(:,:) + za2 * zu_trd(:,:) * ssumask(:,:)
+            zpvotrdv(:,:) = zpvotrdv(:,:) + za2 * zv_trd(:,:) * ssvmask(:,:)
+         ENDIF
+         !
          ! Add tidal astronomical forcing if defined
          IF ( ln_tide .AND. ln_tide_pot ) THEN
@@ -1077,4 +1118,10 @@
                END DO
             END DO
+            !
+            IF( l_trddyn ) THEN
+               za2 = wgtbtp2(jn)
+               zbfrtrdu(:,:) = zbfrtrdu(:,:) + za2 * zCdU_u(:,:) * un_e(:,:) * hur_e(:,:)
+               zbfrtrdv(:,:) = zbfrtrdv(:,:) + za2 * zCdU_v(:,:) * vn_e(:,:) * hvr_e(:,:)
+            ENDIF
          ENDIF 
          !
@@ -1101,5 +1148,10 @@
            END DO
          END IF
-
+         !
+         IF( l_trddyn ) THEN
+            za2 = wgtbtp2(jn)
+            zspgtrdu(:,:) = zspgtrdu(:,:) + za2 * zwx(:,:) * ssumask(:,:)
+            zspgtrdv(:,:) = zspgtrdv(:,:) + za2 * zwy(:,:) * ssvmask(:,:)
+         ENDIF
          !
          ! Set next velocities:
@@ -1302,4 +1354,10 @@
       IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
       IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
+      IF( l_trddyn ) THEN
+         CALL trd_dyn( zspgtrdu, zspgtrdv, jpdyn_spg, kt )
+         CALL trd_dyn( zpvotrdu, zpvotrdv, jpdyn_pvo, kt )
+         CALL trd_dyn( zbfrtrdu, zbfrtrdv, jpdyn_bfr, kt )
+         DEALLOCATE( zspgtrdu, zspgtrdv, zpvotrdu, zpvotrdv, zbfrtrdu, zbfrtrdv )
+      ENDIF
       !
       IF( ln_diatmb ) THEN
Index: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynvor.F90	(revision 11612)
+++ NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynvor.F90	(revision 11613)
@@ -133,5 +133,5 @@
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo_save, kt )
          !
          IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
Index: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynzdf.F90	(revision 11612)
+++ NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynzdf.F90	(revision 11613)
@@ -25,4 +25,5 @@
    !
    USE in_out_manager ! I/O manager
+   USE iom             ! IOM library
    USE lib_mpp        ! MPP library
    USE prtctl         ! Print control
@@ -69,4 +70,5 @@
       INTEGER  ::   ji, jj, jk         ! dummy loop indices
       INTEGER  ::   iku, ikv           ! local integers
+      INTEGER  ::   ikbu, ikbv         ! local integers
       REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
       REAL(wp) ::   zzws, ze3va        !   -      -
@@ -99,10 +101,4 @@
       !
       !
-      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
-         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
-         ztrdu(:,:,:) = ua(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:)
-      ENDIF
-      !
       !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va)
       !
@@ -131,4 +127,9 @@
             va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
          END DO
+         IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
+            ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
+            ztrdu(:,:,:) = ua(:,:,:)
+            ztrdv(:,:,:) = va(:,:,:)
+         ENDIF
          DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -141,4 +142,9 @@
             END DO
          END DO
+         IF( l_trddyn )   THEN                      ! save explicit part of bottom friction trends
+            ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / r2dt 
+            ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / r2dt 
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfr, kt )
+         ENDIF
          IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
             DO jj = 2, jpjm1        
@@ -153,4 +159,10 @@
             END DO
          END IF
+      ENDIF
+      !
+      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
+         IF( .NOT. ALLOCATED(ztrdu) ) ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
+         ztrdu(:,:,:) = ua(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:)
       ENDIF
       !
@@ -486,7 +498,26 @@
       !
       IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
-         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
-         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
+
+         ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / r2dt 
+         ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / r2dt 
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
+         !
+         IF( ln_drgimp ) THEN
+            ztrdu(:,:,:) = 0._wp    ;   ztrdv(:,:,:) = 0._wp
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = 2, jpim1
+                     ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
+                     ikbv = mbkv(ji,jj)
+                     ztrdu(ji,jj,jk) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) )  & 
+     &                               * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu)
+                     ztrdv(ji,jj,jk) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) )  &
+     &                               * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv)
+                  END DO
+               END DO
+            END DO
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfri, kt )
+         ENDIF
+         !
          DEALLOCATE( ztrdu, ztrdv ) 
       ENDIF
Index: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/TRD/trd_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/TRD/trd_oce.F90	(revision 11612)
+++ NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/TRD/trd_oce.F90	(revision 11613)
@@ -59,19 +59,23 @@
    !
    !                                                  !!!* Momentum trends indices
-   INTEGER, PUBLIC, PARAMETER ::   jptot_dyn  = 13     !: Total trend nb: change it when adding/removing one indice below
+   INTEGER, PUBLIC, PARAMETER ::   jptot_dyn  = 17     !: Total trend nb: change it when adding/removing one indice below
    !                               ===============     !  
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg  =  1     !: hydrostatic pressure gradient 
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_spg  =  2     !: surface     pressure gradient
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_keg  =  3     !: kinetic energy gradient  or horizontal advection
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_rvo  =  4     !: relative  vorticity      or metric term
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_pvo  =  5     !: planetary vorticity
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_zad  =  6     !: vertical advection
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_ldf  =  7     !: horizontal diffusion   
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_zdf  =  8     !: vertical   diffusion
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfr  =  9     !: bottom  stress 
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_atf  = 10     !: Asselin time filter
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_tau  = 11     !: surface stress
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfri = 12     !: implicit bottom friction (ln_drgimp=.TRUE.)
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_ken  = 13     !: use for calculation of KE
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg   =  1     !: hydrostatic pressure gradient 
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg_save =  2  !: hydrostatic pressure gradient (saved value)
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg_corr =  3  !: hydrostatic pressure gradient (initial correction)
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_spg   =  4     !: surface     pressure gradient
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_keg   =  5     !: kinetic energy gradient  or horizontal advection
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_rvo   =  6     !: relative  vorticity      or metric term
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_pvo   =  7     !: planetary vorticity
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_pvo_save =  8  !: planetary vorticity (saved value)
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_pvo_corr =  9  !: planetary vorticity (initial correction)
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_zad   = 10     !: vertical advection
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_ldf   = 11     !: horizontal diffusion   
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_zdf   = 12     !: vertical   diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfr   = 13     !: bottom  stress 
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_atf   = 14     !: Asselin time filter
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_tau   = 15     !: surface stress
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfri  = 16     !: implicit bottom friction (ln_drgimp=.TRUE.)
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_ken   = 17     !: use for calculation of KE
    !
    !!----------------------------------------------------------------------
Index: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/TRD/trddyn.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/TRD/trddyn.F90	(revision 11612)
+++ NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/TRD/trddyn.F90	(revision 11613)
@@ -18,5 +18,4 @@
    USE sbc_oce        ! surface boundary condition: ocean
    USE zdf_oce        ! ocean vertical physics: variables
-!!gm   USE zdfdrg         ! ocean vertical physics: bottom friction
    USE trd_oce        ! trends: ocean variables
    USE trdken         ! trends: Kinetic ENergy 
@@ -35,4 +34,11 @@
    PUBLIC trd_dyn        ! called by all dynXXX modules
 
+   INTERFACE trd_dyn
+      module procedure trd_dyn_3d, trd_dyn_2d
+   END INTERFACE
+
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_hpg, zvtrd_hpg
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_pvo, zvtrd_pvo
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -44,7 +50,7 @@
 CONTAINS
 
-   SUBROUTINE trd_dyn( putrd, pvtrd, ktrd, kt )
+   SUBROUTINE trd_dyn_3d( putrd, pvtrd, ktrd, kt )
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_mod  ***
+      !!                  ***  ROUTINE trd_dyn_3d  ***
       !! 
       !! ** Purpose :   Dispatch momentum trend computation, e.g. 3D output, 
@@ -63,23 +69,147 @@
 !!gm NB : here a lbc_lnk should probably be added
 
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !   3D output of momentum and/or tracers trends using IOM interface
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt )
+      SELECT CASE( ktrd )
+      CASE( jpdyn_hpg_save ) 
+         !
+         ! save 3D HPG trends to possibly have barotropic part corrected later before writing out
+         ALLOCATE( zutrd_hpg(jpi,jpj,jpk), zvtrd_hpg(jpi,jpj,jpk) )
+         zutrd_hpg(:,:,:) = putrd(:,:,:)
+         zvtrd_hpg(:,:,:) = pvtrd(:,:,:)
+
+      CASE( jpdyn_pvo_save ) 
+         !
+         ! save 3D coriolis trends to possibly have barotropic part corrected later before writing out
+         ALLOCATE( zutrd_pvo(jpi,jpj,jpk), zvtrd_pvo(jpi,jpj,jpk) )
+         zutrd_pvo(:,:,:) = putrd(:,:,:)
+         zvtrd_pvo(:,:,:) = pvtrd(:,:,:)
+
+      CASE DEFAULT
+
+         IF( ktrd == jpdyn_spg ) THEN
+            ! For explicit scheme SPG trends come here as 3D fields
+            ! Add SPG trend to 3D HPG trend and also output as 2D diagnostic in own right.
+            CALL trd_dyn_iom_2d( putrd(:,:,1), pvtrd(:,:,1), jpdyn_spg, kt ) 
+            putrd(:,:,:) = putrd(:,:,:) + zutrd_hpg(:,:,:) 
+            pvtrd(:,:,:) = pvtrd(:,:,:) + zvtrd_hpg(:,:,:) 
+            DEALLOCATE( zutrd_hpg, zvtrd_hpg )
+         ENDIF
+
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         !   3D output of momentum and/or tracers trends using IOM interface
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         IF( ln_dyn_trd )   CALL trd_dyn_iom_3d( putrd, pvtrd, ktrd, kt )
+
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         !  Integral Constraints Properties for momentum and/or tracers trends
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt )
+
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         !  Kinetic Energy trends
+         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+         IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt )
+
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         !  Vorticity trends
+         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+         IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt )
+
+         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+         !  Mixed layer trends for active tracers
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         !!gm      IF( ln_dyn_mxl )   CALL trd_mxl_dyn   
+         !
+      END SELECT
+      !
+   END SUBROUTINE trd_dyn_3d
+
+
+   SUBROUTINE trd_dyn_2d( putrd, pvtrd, ktrd, kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mod  ***
+      !! 
+      !! ** Purpose :   Dispatch momentum trend computation, e.g. 2D output, 
+      !!              integral constraints, barotropic vorticity, kinetic enrgy, 
+      !!              and/or mixed layer budget.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends 
+      INTEGER                 , INTENT(in   ) ::   ktrd           ! trend index
+      INTEGER                 , INTENT(in   ) ::   kt             ! time step
+      INTEGER                                 ::   jk
+      !!----------------------------------------------------------------------
+      !
+      putrd(:,:) = putrd(:,:) * umask(:,:,1)                       ! mask the trends
+      pvtrd(:,:) = pvtrd(:,:) * vmask(:,:,1)
+      !
+
+!!gm NB : here a lbc_lnk should probably be added
+
+      SELECT CASE(ktrd)
+
+      CASE ( jpdyn_hpg_corr )
+         !
+         ! Remove "first-guess" SPG trend from 3D HPG trend. 
+         DO jk = 1, jpkm1
+            zutrd_hpg(:,:,jk) = zutrd_hpg(:,:,jk) - putrd(:,:)
+            zvtrd_hpg(:,:,jk) = zvtrd_hpg(:,:,jk) - pvtrd(:,:)
+         ENDDO
+         CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_hpg_corr, kt ) 
+
+      CASE( jpdyn_pvo_corr )
+         !
+         ! Remove "first-guess" barotropic coriolis trend from 3D PVO trend. 
+         DO jk = 1, jpkm1
+            zutrd_pvo(:,:,jk) = zutrd_pvo(:,:,jk) - putrd(:,:)
+            zvtrd_pvo(:,:,jk) = zvtrd_pvo(:,:,jk) - pvtrd(:,:)
+         ENDDO
+         CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_pvo_corr, kt ) 
+
+      CASE( jpdyn_spg )
+          !
+          ! For split-explicit scheme SPG trends come here as 2D fields
+          ! Add SPG trend to 3D HPG trend and also output as 2D diagnostic in own right.
+          DO jk = 1, jpkm1
+             zutrd_hpg(:,:,jk) = zutrd_hpg(:,:,jk) + putrd(:,:)
+             zvtrd_hpg(:,:,jk) = zvtrd_hpg(:,:,jk) + pvtrd(:,:)
+          ENDDO
+          CALL trd_dyn_iom_2d( putrd(:,:), pvtrd(:,:), jpdyn_spg, kt ) 
+          CALL trd_dyn_3d( zutrd_hpg, zvtrd_hpg, jpdyn_hpg, kt )
+          DEALLOCATE( zutrd_hpg, zvtrd_hpg )
+
+      CASE( jpdyn_pvo )
+          !
+          ! Add 2D PVO trend to 3D PVO trend and also output as diagnostic in own right.
+          DO jk = 1, jpkm1
+             zutrd_pvo(:,:,jk) = zutrd_pvo(:,:,jk) + putrd(:,:)
+             zvtrd_pvo(:,:,jk) = zvtrd_pvo(:,:,jk) + pvtrd(:,:)
+          ENDDO
+          CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_pvo, kt ) 
+          CALL trd_dyn_3d( zutrd_pvo, zvtrd_pvo, jpdyn_pvo, kt )
+          DEALLOCATE( zutrd_pvo, zvtrd_pvo )
+
+      CASE DEFAULT 
+
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         !   2D output of momentum and/or tracers trends using IOM interface
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         IF( ln_dyn_trd )   CALL trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt )
          
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !  Integral Constraints Properties for momentum and/or tracers trends
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt )
-
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !  Kinetic Energy trends
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt )
-
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !  Vorticity trends
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt )
+      END SELECT
+
+!!$   CALLS TO THESE ROUTINES FOR 2D DIAGOSTICS NOT CODED YET
+!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!$      !  Integral Constraints Properties for momentum and/or tracers trends
+!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!$      IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt )
+!!$
+!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!$      !  Kinetic Energy trends
+!!$      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!!$      IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt )
+!!$
+!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!$      !  Vorticity trends
+!!$      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!!$      IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -88,8 +218,8 @@
 !!gm      IF( ln_dyn_mxl )   CALL trd_mxl_dyn   
       !
-   END SUBROUTINE trd_dyn
-
-
-   SUBROUTINE trd_dyn_iom( putrd, pvtrd, ktrd, kt )
+   END SUBROUTINE trd_dyn_2d
+
+
+   SUBROUTINE trd_dyn_iom_3d( putrd, pvtrd, ktrd, kt )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_dyn_iom  ***
@@ -110,6 +240,4 @@
       CASE( jpdyn_hpg )   ;   CALL iom_put( "utrd_hpg", putrd )    ! hydrostatic pressure gradient
                               CALL iom_put( "vtrd_hpg", pvtrd )
-      CASE( jpdyn_spg )   ;   CALL iom_put( "utrd_spg", putrd )    ! surface pressure gradient
-                              CALL iom_put( "vtrd_spg", pvtrd )
       CASE( jpdyn_pvo )   ;   CALL iom_put( "utrd_pvo", putrd )    ! planetary vorticity
                               CALL iom_put( "vtrd_pvo", pvtrd )
@@ -147,36 +275,45 @@
                               CALL iom_put( "vtrd_tau", z2dy )
                               DEALLOCATE( z2dx , z2dy )
-!!gm  to be changed : computation should be done in dynzdf.F90
-!!gm                + missing the top friction 
-!                              !                                    ! bottom stress tends (implicit case)
-!                              IF( ln_drgimp ) THEN
-!                                 ALLOCATE( z3dx(jpi,jpj,jpk) , z3dy(jpi,jpj,jpk) )
-!                             z3dx(:,:,:) = 0._wp   ;   z3dy(:,:,:) = 0._wp  ! after velocity known (now filed at this stage)
-!                            DO jk = 1, jpkm1
-!                                    DO jj = 2, jpjm1
-!                                       DO ji = 2, jpim1
-!                                      ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
-!                                          ikbv = mbkv(ji,jj)
-!                                          z3dx(ji,jj,jk) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) ) & 
-!                                               &         * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu)
-!                                          z3dy(ji,jj,jk) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) ) &
-!                                               &         * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv)
-!                                    END DO
-!                                 END DO
-!                              END DO
-!                              CALL lbc_lnk_multi( 'trddyn', z3dx, 'U', -1., z3dy, 'V', -1. )
-!                              CALL iom_put( "utrd_bfr", z3dx )
-!                              CALL iom_put( "vtrd_bfr", z3dy )
-!                                 DEALLOCATE( z3dx , z3dy )
-!                              ENDIF
-!!gm end
-      CASE( jpdyn_bfr )       ! called if ln_drgimp=F
-                              CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case)
+      CASE( jpdyn_bfr )   ;   CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case)
                               CALL iom_put( "vtrd_bfr", pvtrd )
+      CASE( jpdyn_bfri)   ;   CALL iom_put( "utrd_bfri", putrd )    ! bottom friction (implicit case)
+                              CALL iom_put( "vtrd_bfri", pvtrd )
       CASE( jpdyn_atf )   ;   CALL iom_put( "utrd_atf", putrd )        ! asselin filter trends 
                               CALL iom_put( "vtrd_atf", pvtrd )
       END SELECT
       !
-   END SUBROUTINE trd_dyn_iom
+   END SUBROUTINE trd_dyn_iom_3d
+
+
+   SUBROUTINE trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_dyn_iom  ***
+      !! 
+      !! ** Purpose :   output 2D trends using IOM
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends
+      INTEGER                 , INTENT(in   ) ::   ktrd           ! trend index
+      INTEGER                 , INTENT(in   ) ::   kt             ! time step
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      INTEGER ::   ikbu, ikbv   ! local integers
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3dx, z3dy   ! 3D workspace 
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE( ktrd )
+      CASE( jpdyn_spg )      ;   CALL iom_put( "utrd_spg2d", putrd )      ! surface pressure gradient
+                                 CALL iom_put( "vtrd_spg2d", pvtrd )
+      CASE( jpdyn_pvo )      ;   CALL iom_put( "utrd_pvo2d", putrd )      ! planetary vorticity (barotropic part)
+                                 CALL iom_put( "vtrd_pvo2d", pvtrd )
+      CASE( jpdyn_hpg_corr ) ;   CALL iom_put( "utrd_hpg_corr", putrd )   ! horizontal pressure gradient correction
+                                 CALL iom_put( "vtrd_hpg_corr", pvtrd )
+      CASE( jpdyn_pvo_corr ) ;   CALL iom_put( "utrd_pvo_corr", putrd )   ! planetary vorticity correction
+                                 CALL iom_put( "vtrd_pvo_corr", pvtrd )
+      CASE( jpdyn_bfr )      ;   CALL iom_put( "utrd_bfr2d", putrd )      ! bottom friction due to barotropic currents
+                                 CALL iom_put( "vtrd_bfr2d", pvtrd )
+      END SELECT
+      !
+   END SUBROUTINE trd_dyn_iom_2d
 
    !!======================================================================
