Index: NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11904)
+++ NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11934)
@@ -806,10 +806,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"                        />
@@ -819,6 +821,22 @@
      <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_group>
 
@@ -826,12 +844,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"                        />
@@ -839,6 +861,22 @@
      <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_group>
 
Index: NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynhpg.F90	(revision 11904)
+++ NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynhpg.F90	(revision 11934)
@@ -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.1_momentum_trends/src/OCE/DYN/dynspg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynspg.F90	(revision 11904)
+++ NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynspg.F90	(revision 11934)
@@ -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.1_momentum_trends/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynspg_ts.F90	(revision 11904)
+++ NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynspg_ts.F90	(revision 11934)
@@ -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
 
    USE iom   ! to remove
@@ -170,4 +172,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(:,:) :: ztautrdu, ztautrdv, zbfrtrdu, zbfrtrdv  ! TAU and BFR trends (if l_trddyn)
       !!----------------------------------------------------------------------
       !
@@ -175,4 +179,22 @@
       !                                         !* 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), &
+         &          ztautrdu(jpi,jpj), ztautrdv(jpi,jpj), zbfrtrdu(jpi,jpj), zbfrtrdv(jpi,jpj) )
+          zspgtrdu(:,:) = 0._wp
+          zspgtrdv(:,:) = 0._wp
+          zpvotrdu(:,:) = 0._wp
+          zpvotrdv(:,:) = 0._wp
+          ztautrdu(:,:) = 0._wp
+          ztautrdv(:,:) = 0._wp
+          zbfrtrdu(:,:) = 0._wp
+          zbfrtrdv(:,:) = 0._wp
+      ENDIF
+      !
+      zu_trd(:,:) = 0._wp
+      zv_trd(:,:) = 0._wp
+      zu_spg(:,:) = 0._wp
+      zv_spg(:,:) = 0._wp
       !
       zmdi=1.e+20                               !  missing data indicator for masking
@@ -253,5 +275,15 @@
          &                               zu_trd, zv_trd   )   ! ==>> out
       !
+      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
+      !
       IF( .NOT.ln_linssh ) THEN                 !* surface pressure gradient   (variable volume only)
+         !
+         IF( l_trddyn ) THEN
+            zspgtrdu(:,:) = zu_trd(:,:)
+            zspgtrdv(:,:) = zv_trd(:,:)
+         ENDIF
          !
          IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg
@@ -274,4 +306,14 @@
          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
+         !
       ENDIF
       !
@@ -283,7 +325,17 @@
       END DO 
       !
+      IF( l_trddyn ) THEN
+         zbfrtrdu(:,:) = zu_frc(:,:)
+         zbfrtrdv(:,:) = zv_frc(:,:)
+      ENDIF
       !                                   !=  Add bottom stress contribution from baroclinic velocities  =!
       !                                   !  -----------------------------------------------------------  !
       CALL dyn_drg_init( zu_frc, zv_frc,  zCdU_u, zCdU_v )      ! also provide the barotropic drag coefficients
+      !
+      IF( l_trddyn ) THEN
+         ! bottom friction trend diagnostic: bottom friction due to baroclinic currents
+         zbfrtrdu(:,:) = zu_frc(:,:) - zbfrtrdu(:,:)
+         zbfrtrdv(:,:) = zv_frc(:,:) - zbfrtrdv(:,:) 
+      ENDIF
       !
       !                                   !=  Add atmospheric pressure forcing  =!
@@ -312,4 +364,9 @@
       !                                   !=  Add atmospheric pressure forcing  =!
       !                                   !  ----------------------------------  !
+      IF( l_trddyn ) THEN
+         ztautrdu(:,:) = zu_frc(:,:)
+         ztautrdv(:,:) = zv_frc(:,:)
+      ENDIF
+      !
       IF( ln_bt_fw ) THEN                        ! Add wind forcing
          DO jj = 2, jpjm1
@@ -329,4 +386,9 @@
       ENDIF  
       !
+      IF( l_trddyn ) THEN
+         ! wind stress trend diagnostic
+         ztautrdu(:,:) = zu_frc(:,:) - ztautrdu(:,:)
+         ztautrdv(:,:) = zv_frc(:,:) - ztautrdv(:,:) 
+      ENDIF
       !              !----------------!
       !              !==  sssh_frc  ==!   Right-Hand-Side of the barotropic ssh equation   (over the FULL domain)
@@ -591,4 +653,10 @@
          ENDIF
          !
+         IF( l_trddyn ) THEN
+            za2 = wgtbtp2(jn)
+            zspgtrdu(:,:) = zspgtrdu(:,:) + za2 * zu_spg(:,:) * ssumask(:,:)
+            zspgtrdv(:,:) = zspgtrdv(:,:) + za2 * zv_spg(:,:) * ssvmask(:,:)
+         ENDIF
+         !
          ! Add Coriolis trend:
          ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
@@ -596,4 +664,10 @@
          ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
          CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   )
+         !
+         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
@@ -616,4 +690,9 @@
                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
          !
@@ -835,4 +914,12 @@
       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( ztautrdu, ztautrdv, jpdyn_tau, kt )
+         CALL trd_dyn( zbfrtrdu, zbfrtrdv, jpdyn_bfr, kt )
+         DEALLOCATE( zspgtrdu, zspgtrdv, zpvotrdu, zpvotrdv, ztautrdu, ztautrdv, zbfrtrdu, zbfrtrdv )
+      ENDIF
       !
       IF( ln_diatmb ) THEN
Index: NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynvor.F90	(revision 11904)
+++ NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynvor.F90	(revision 11934)
@@ -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.1_momentum_trends/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynzdf.F90	(revision 11904)
+++ NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/DYN/dynzdf.F90	(revision 11934)
@@ -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.1_momentum_trends/src/OCE/TRD/trd_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/TRD/trd_oce.F90	(revision 11904)
+++ NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/TRD/trd_oce.F90	(revision 11934)
@@ -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.1_momentum_trends/src/OCE/TRD/trddyn.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/TRD/trddyn.F90	(revision 11904)
+++ NEMO/branches/UKMO/NEMO_4.0.1_momentum_trends/src/OCE/TRD/trddyn.F90	(revision 11934)
@@ -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,13 @@
    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
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_bfr, zvtrd_bfr
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:)  , SAVE :: zutrd_tau, zvtrd_tau
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -44,7 +52,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, 
@@ -55,4 +63,6 @@
       INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
       INTEGER                   , INTENT(in   ) ::   kt             ! time step
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zue, zve       ! temporary 2D arrays
+      INTEGER                                   ::   jk
       !!----------------------------------------------------------------------
       !
@@ -63,23 +73,197 @@
 !!gm NB : here a lbc_lnk should probably be added
 
+      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( jpdyn_spg ) 
+         ! 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 )
+
+      CASE( jpdyn_bfr )
+         !
+         ! Add 3D part of BFR trend minus its depth-mean part to depth-mean trend already saved.
+         ALLOCATE( zue(jpi,jpj), zve(jpi,jpj) )
+         zue(:,:) = e3u_a(:,:,1) * putrd(:,:,1) * umask(:,:,1)
+         zve(:,:) = e3v_a(:,:,1) * pvtrd(:,:,1) * vmask(:,:,1)
+         DO jk = 2, jpkm1
+            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * putrd(:,:,jk) * umask(:,:,jk)
+            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * pvtrd(:,:,jk) * vmask(:,:,jk)
+         END DO
+         DO jk = 1, jpkm1
+            putrd(:,:,jk) = zutrd_bfr(:,:,jk) + putrd(:,:,jk) - zue(:,:) * r1_hu_a(:,:)
+            pvtrd(:,:,jk) = zvtrd_bfr(:,:,jk) + pvtrd(:,:,jk) - zve(:,:) * r1_hv_a(:,:)
+         END DO
+         ! Update locally saved BFR trends to add to ZDF trend.
+         zutrd_bfr(:,:,:) = putrd(:,:,:) 
+         zvtrd_bfr(:,:,:) = pvtrd(:,:,:)
+
+      CASE( jpdyn_zdf ) 
+         ! ZDF trend: Remove barotropic component and add wind stress and bottom friction
+         !            trends from dynspg_ts. Also adding on the bottom stress for the 
+         !            baroclinic solution in the case of explicit bottom friction. 
+         ALLOCATE( zue(jpi,jpj), zve(jpi,jpj) )
+         zue(:,:) = e3u_a(:,:,1) * putrd(:,:,1) * umask(:,:,1)
+         zve(:,:) = e3v_a(:,:,1) * pvtrd(:,:,1) * vmask(:,:,1)
+         DO jk = 2, jpkm1
+            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * putrd(:,:,jk) * umask(:,:,jk)
+            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * pvtrd(:,:,jk) * vmask(:,:,jk)
+         END DO
+         DO jk = 1, jpkm1
+            putrd(:,:,jk) = zutrd_tau(:,:) + zutrd_bfr(:,:,jk) + putrd(:,:,jk) - zue(:,:) * r1_hu_a(:,:)
+            pvtrd(:,:,jk) = zutrd_tau(:,:) + zvtrd_bfr(:,:,jk) + pvtrd(:,:,jk) - zve(:,:) * r1_hv_a(:,:)
+         END DO
+         DEALLOCATE( zue, zve, zutrd_tau, zvtrd_tau, zutrd_bfr, zvtrd_bfr )
+
+      END SELECT
+
+      IF ( ktrd /= jpdyn_hpg_save .AND. ktrd /= jpdyn_pvo_save ) THEN
+         !
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+         !   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   
+         !
+      ENDIF
+      !
+   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
+
+      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
+
+      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_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_3d( zutrd_pvo, zvtrd_pvo, jpdyn_pvo, kt )
+          DEALLOCATE( zutrd_pvo, zvtrd_pvo )
+
+      CASE( jpdyn_tau )
+          !
+          ! Save 2D wind forcing trend locally (to be added to ZDF trend)
+          ! and output as a trend in its own right.
+          ALLOCATE( zutrd_tau(jpi,jpj), zvtrd_tau(jpi,jpj) )
+          zutrd_tau(:,:) = putrd(:,:)
+          zvtrd_tau(:,:) = pvtrd(:,:)
+
+      CASE( jpdyn_bfr )
+          !
+          ! Create 3D BFR trend from 2D field and also output 2D field as diagnostic in own right.
+          ALLOCATE( zutrd_bfr(jpi,jpj,jpk), zvtrd_bfr(jpi,jpj,jpk) )
+          zutrd_bfr(:,:,:) = 0.0
+          zvtrd_bfr(:,:,:) = 0.0
+          DO jk = 1, jpkm1
+             zutrd_bfr(:,:,jk) = putrd(:,:) * umask(:,:,jk)
+             zvtrd_bfr(:,:,jk) = pvtrd(:,:) * vmask(:,:,jk)
+          ENDDO
+
+      END SELECT
+
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !   3D output of momentum and/or tracers trends using IOM interface
+      !   2D output of momentum and/or tracers trends using IOM interface
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt )
+      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 )
+
+!!$   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 +272,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 +294,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 +329,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
 
    !!======================================================================
