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 11916)
+++ /NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynspg_ts.F90	(revision 11917)
@@ -175,5 +175,5 @@
       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)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztautrdu, ztautrdv, zbfrtrdu, zbfrtrdv  ! TAU and BFR trends (if l_trddyn)
       !!----------------------------------------------------------------------
       !
@@ -183,9 +183,12 @@
       !
       IF( l_trddyn ) THEN
-          ALLOCATE( zspgtrdu(jpi,jpj), zspgtrdv(jpi,jpj), zpvotrdu(jpi,jpj), zpvotrdv(jpi,jpj), zbfrtrdu(jpi,jpj), zbfrtrdv(jpi,jpj) )
+          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
@@ -573,4 +576,9 @@
       !
       ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
+      IF( l_trddyn ) THEN
+         zbfrtrdu(:,:) = zu_frc(:,:)
+         zbfrtrdv(:,:) = zv_frc(:,:)
+      ENDIF
+      !
       IF( ln_wd_il ) THEN
          zztmp = -1._wp / rdtbt
@@ -592,4 +600,10 @@
       END IF
       !
+      IF( l_trddyn ) THEN
+         ! bottom friction trend diagnostic: bottom friction due to baroclinic currents
+         zbfrtrdu(:,:) = zu_frc(:,:) - zbfrtrdu(:,:)
+         zbfrtrdv(:,:) = zv_frc(:,:) - zbfrtrdv(:,:) 
+      ENDIF
+      !
       IF( ln_isfcav ) THEN       ! Add TOP stress contribution from baroclinic velocities:      
          IF( ln_bt_fw ) THEN
@@ -622,4 +636,9 @@
       ENDIF
       !       
+      IF( l_trddyn ) THEN
+         ztautrdu(:,:) = zu_frc(:,:)
+         ztautrdv(:,:) = zv_frc(:,:)
+      ENDIF
+      !
       IF( ln_bt_fw ) THEN                        ! Add wind forcing
          DO jj = 2, jpjm1
@@ -638,4 +657,10 @@
          END DO
       ENDIF  
+      !
+      IF( l_trddyn ) THEN
+         ! wind stress trend diagnostic
+         ztautrdu(:,:) = zu_frc(:,:) - ztautrdu(:,:)
+         ztautrdv(:,:) = zv_frc(:,:) - ztautrdv(:,:) 
+      ENDIF
       !
       IF( ln_apr_dyn ) THEN                     ! Add atm pressure forcing
@@ -1357,6 +1382,7 @@
          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, zbfrtrdu, zbfrtrdv )
+         DEALLOCATE( zspgtrdu, zspgtrdv, zpvotrdu, zpvotrdv, ztautrdu, ztautrdv, zbfrtrdu, zbfrtrdv )
       ENDIF
       !
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 11916)
+++ /NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/TRD/trddyn.F90	(revision 11917)
@@ -41,4 +41,5 @@
    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
@@ -87,32 +88,51 @@
          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
-
-         IF( ktrd == jpdyn_bfr ) THEN
-            !
-            ! 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
-            DEALLOCATE( zue, zve, zutrd_bfr, zvtrd_bfr )
-         ENDIF
-
+      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
@@ -140,5 +160,5 @@
          !!gm      IF( ln_dyn_mxl )   CALL trd_mxl_dyn   
          !
-      END SELECT
+      ENDIF
       !
    END SUBROUTINE trd_dyn_3d
@@ -174,5 +194,4 @@
             zvtrd_hpg(:,:,jk) = zvtrd_hpg(:,:,jk) - pvtrd(:,:)
          ENDDO
-         CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_hpg_corr, kt ) 
 
       CASE( jpdyn_pvo_corr )
@@ -183,5 +202,4 @@
             zvtrd_pvo(:,:,jk) = zvtrd_pvo(:,:,jk) - pvtrd(:,:)
          ENDDO
-         CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_pvo_corr, kt ) 
 
       CASE( jpdyn_spg )
@@ -193,5 +211,4 @@
              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 )
@@ -204,7 +221,14 @@
              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( 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 )
@@ -215,17 +239,15 @@
           zvtrd_bfr(:,:,:) = 0.0
           DO jk = 1, jpkm1
-             zutrd_bfr(:,:,jk) = putrd(:,:)
-             zvtrd_bfr(:,:,jk) = pvtrd(:,:)
+             zutrd_bfr(:,:,jk) = putrd(:,:) * umask(:,:,jk)
+             zvtrd_bfr(:,:,jk) = pvtrd(:,:) * vmask(:,:,jk)
           ENDDO
-          CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_bfr, kt ) 
-
-      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 )
+
+      END SELECT
+
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !   2D output of momentum and/or tracers trends using IOM interface
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_dyn_trd )   CALL trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt )
          
-      END SELECT
 
 !!$   CALLS TO THESE ROUTINES FOR 2D DIAGOSTICS NOT CODED YET
