Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 2715)
@@ -45,7 +45,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -57,33 +56,25 @@
       !!              momentum equations at every time step frequency nn_trd.
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dx             ! Temperature or U trend 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dy             ! Salinity    or V trend
-      INTEGER                     , INTENT(in   ) ::   ktrd                ! tracer trend index
-      CHARACTER(len=3)            , INTENT(in   ) ::   ctype               ! momentum ('DYN') or tracers ('TRA') trends
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dx   ! Temperature or U trend 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dy   ! Salinity    or V trend
+      INTEGER                     , INTENT(in   ) ::   ktrd      ! tracer trend index
+      CHARACTER(len=3)            , INTENT(in   ) ::   ctype     ! momentum ('DYN') or tracers ('TRA') trends
       !!
-      INTEGER  ::   ji, jj                                                 ! loop indices
-      REAL(wp) ::   zmsku, zbtu, zbt                                       ! temporary scalars
-      REAL(wp) ::   zmskv, zbtv                                            !    "         "
-      !!----------------------------------------------------------------------
-
-
-      ! 1. Mask trends
-      ! --------------
-
-      SELECT CASE( ctype )
-      !
-      CASE( 'DYN' )              ! Momentum
+      INTEGER ::   ji, jj   ! loop indices
+      !!----------------------------------------------------------------------
+
+      SELECT CASE( ctype )    !==  Mask trends  ==!
+      !
+      CASE( 'DYN' )                    ! Momentum
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
-               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
-               ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * zmsku
-               ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * zmskv
+               ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
+               ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
             END DO
          END DO
-         ptrd2dx(jpi, : ) = 0.e0      ;      ptrd2dy(jpi, : ) = 0.e0
-         ptrd2dx( : ,jpj) = 0.e0      ;      ptrd2dy( : ,jpj) = 0.e0
-         !
-      CASE( 'TRA' )              ! Tracers
+         ptrd2dx(jpi, : ) = 0._wp      ;      ptrd2dy(jpi, : ) = 0._wp
+         ptrd2dx( : ,jpj) = 0._wp      ;      ptrd2dy( : ,jpj) = 0._wp
+         !
+      CASE( 'TRA' )                    ! Tracers
          ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
          ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
@@ -91,66 +82,30 @@
       END SELECT
       
-      ! 2. Basin averaged tracer/momentum trends
-      ! ----------------------------------------
-
-      SELECT CASE( ctype )
-      !
-      CASE( 'DYN' )              ! Momentum
-         umo(ktrd) = 0.e0
-         vmo(ktrd) = 0.e0
+      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
+      !
+      CASE( 'DYN' )                    ! Momentum
+         umo(ktrd) = 0._wp
+         vmo(ktrd) = 0._wp
          !
          SELECT CASE( ktrd )
-         !
          CASE( jpdyn_trd_swf )         ! surface forcing
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  umo(ktrd) = umo(ktrd) + ptrd2dx(ji,jj) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
-                  vmo(ktrd) = vmo(ktrd) + ptrd2dy(ji,jj) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
-               END DO
-            END DO
-            !
+            umo(ktrd) = SUM( ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1) )
+            vmo(ktrd) = SUM( ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1) )
          END SELECT
          !
       CASE( 'TRA' )              ! Tracers
-         tmo(ktrd) = 0.e0
-         smo(ktrd) = 0.e0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
-               tmo(ktrd) =  tmo(ktrd) + ptrd2dx(ji,jj) * zbt
-               smo(ktrd) =  smo(ktrd) + ptrd2dy(ji,jj) * zbt
-            END DO
-         END DO
-         !
+         tmo(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
+         smo(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
       END SELECT
       
-      ! 3. Basin averaged tracer/momentum square trends
-      ! ----------------------------------------------
-      ! c a u t i o n: field now
-      
-      SELECT CASE( ctype )
+      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
       !
       CASE( 'DYN' )              ! Momentum
-         hke(ktrd) = 0.e0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
-               zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
-               hke(ktrd) = hke(ktrd)   &
-               &   + un(ji,jj,1) * ptrd2dx(ji,jj) * zbtu &
-               &   + vn(ji,jj,1) * ptrd2dy(ji,jj) * zbtv
-            END DO
-         END DO
+         hke(ktrd) = SUM(   un(:,:,1) * ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1)   &
+            &             + vn(:,:,1) * ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1)   )
          !
       CASE( 'TRA' )              ! Tracers
-         t2(ktrd) = 0.e0
-         s2(ktrd) = 0.e0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
-               t2(ktrd) = t2(ktrd) + ptrd2dx(ji,jj) * zbt * tn(ji,jj,1)
-               s2(ktrd) = s2(ktrd) + ptrd2dy(ji,jj) * zbt * sn(ji,jj,1)
-            END DO
-         END DO
+         t2(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) * tn(:,:,1) )
+         s2(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) * sn(:,:,1) )
          !      
       END SELECT
@@ -166,34 +121,28 @@
       !!              momentum equations at every time step frequency nn_trd.
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dx            ! Temperature or U trend 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dy            ! Salinity    or V trend
-      INTEGER,                          INTENT(in   ) ::   ktrd               ! momentum or tracer trend index
-      CHARACTER(len=3),                 INTENT(in   ) ::   ctype              ! momentum ('DYN') or tracers ('TRA') trends
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dx   ! Temperature or U trend 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dy   ! Salinity    or V trend
+      INTEGER,                          INTENT(in   ) ::   ktrd      ! momentum or tracer trend index
+      CHARACTER(len=3),                 INTENT(in   ) ::   ctype     ! momentum ('DYN') or tracers ('TRA') trends
       !!
-      INTEGER ::   ji, jj, jk
-      REAL(wp) ::   zbt, zbtu, zbtv, zmsku, zmskv                             ! temporary scalars
-      !!----------------------------------------------------------------------
-
-      ! 1. Mask the trends
-      ! ------------------
-
-      SELECT CASE( ctype )
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      SELECT CASE( ctype )    !==  Mask the trends  ==!
       !
       CASE( 'DYN' )              ! Momentum        
-         DO jk = 1, jpk
+         DO jk = 1, jpkm1
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-                  zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
-                  zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
-                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * zmsku
-                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * zmskv
+                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
+                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
                END DO
             END DO
          END DO
-         ptrd3dx(jpi, : ,:) = 0.e0      ;      ptrd3dy(jpi, : ,:) = 0.e0
-         ptrd3dx( : ,jpj,:) = 0.e0      ;      ptrd3dy( : ,jpj,:) = 0.e0
+         ptrd3dx(jpi, : ,:) = 0._wp      ;      ptrd3dy(jpi, : ,:) = 0._wp
+         ptrd3dx( : ,jpj,:) = 0._wp      ;      ptrd3dy( : ,jpj,:) = 0._wp
          !
       CASE( 'TRA' )              ! Tracers
-         DO jk = 1, jpk
+         DO jk = 1, jpkm1
             ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
             ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
@@ -202,69 +151,39 @@
       END SELECT   
 
-      ! 2. Basin averaged tracer/momentum trends
-      ! ----------------------------------------
-      
-      SELECT CASE( ctype )
+      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
       !
       CASE( 'DYN' )              ! Momentum
-         umo(ktrd) = 0.e0
-         vmo(ktrd) = 0.e0
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-                  umo(ktrd) = umo(ktrd) + ptrd3dx(ji,jj,jk) * zbtu
-                  vmo(ktrd) = vmo(ktrd) + ptrd3dy(ji,jj,jk) * zbtv
-               END DO
-            END DO
+         umo(ktrd) = 0._wp
+         vmo(ktrd) = 0._wp
+         DO jk = 1, jpkm1
+            umo(ktrd) = umo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk) )
+            vmo(ktrd) = vmo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk) )
          END DO
          !
       CASE( 'TRA' )              ! Tracers
-         tmo(ktrd) = 0.e0
-         smo(ktrd) = 0.e0
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) 
-                  tmo(ktrd) =  tmo(ktrd) + ptrd3dx(ji,jj,jk) * zbt
-                  smo(ktrd) =  smo(ktrd) + ptrd3dy(ji,jj,jk) * zbt
-               END DO
-            END DO
+         tmo(ktrd) = 0._wp
+         smo(ktrd) = 0._wp
+         DO jk = 1, jpkm1
+            tmo(ktrd) = tmo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
+            smo(ktrd) = smo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
          END DO
          !
       END SELECT
 
-      ! 3. Basin averaged tracer/momentum square trends
-      ! -----------------------------------------------
-      ! c a u t i o n: field now
-      
-      SELECT CASE( ctype )
+      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
       !
       CASE( 'DYN' )              ! Momentum
-         hke(ktrd) = 0.e0
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-                  hke(ktrd) = hke(ktrd)   &
-                  &   + un(ji,jj,jk) * ptrd3dx(ji,jj,jk) * zbtu &
-                  &   + vn(ji,jj,jk) * ptrd3dy(ji,jj,jk) * zbtv
-               END DO
-            END DO
+         hke(ktrd) = 0._wp
+         DO jk = 1, jpkm1
+            hke(ktrd) = hke(ktrd) + SUM(   un(:,:,jk) * ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk)   &
+               &                         + vn(:,:,jk) * ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk)   )
          END DO
          !
       CASE( 'TRA' )              ! Tracers
-         t2(ktrd) = 0.e0
-         s2(ktrd) = 0.e0
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-                  t2(ktrd) = t2(ktrd) + ptrd3dx(ji,jj,jk) * zbt * tn(ji,jj,jk)
-                  s2(ktrd) = s2(ktrd) + ptrd3dy(ji,jj,jk) * zbt * sn(ji,jj,jk)
-               END DO
-            END DO
+         t2(ktrd) = 0._wp
+         s2(ktrd) = 0._wp
+         DO jk = 1, jpkm1
+            t2(ktrd) = t2(ktrd) + SUM( ptrd3dx(ji,jj,jk) * tn(ji,jj,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
+            s2(ktrd) = s2(ktrd) + SUM( ptrd3dy(ji,jj,jk) * sn(ji,jj,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
          END DO
          !
@@ -272,5 +191,4 @@
       !
    END SUBROUTINE trd_3d
-
 
 
@@ -281,9 +199,5 @@
       !! ** Purpose :   Read the namtrd namelist
       !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zmskt
-#if  defined key_trddyn
-      REAL(wp) ::   zmsku, zmskv
-#endif
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
       !!----------------------------------------------------------------------
 
@@ -295,12 +209,7 @@
 
       ! Total volume at t-points:
-      tvolt = 0.e0
+      tvolt = 0._wp
       DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zmskt = tmask(ji,jj,jk) * tmask_i(ji,jj)
-               tvolt = tvolt + zmskt * e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk)
-            END DO
-         END DO
+         tvolt = SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
       END DO
       IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
@@ -310,17 +219,15 @@
 #if  defined key_trddyn
       ! Initialization of potential to kinetic energy conversion
-      rpktrd = 0.e0
+      rpktrd = 0._wp
 
       ! Total volume at u-, v- points:
-      tvolu = 0.e0
-      tvolv = 0.e0
+      tvolu = 0._wp
+      tvolv = 0._wp
 
       DO jk = 1, jpk
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
-               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
-               tvolu = tvolu + zmsku * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-               tvolv = tvolv + zmskv * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
+               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
+               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
             END DO
          END DO
@@ -344,10 +251,17 @@
       !! ** Purpose :  write dynamic trends in ocean.output 
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt                                  ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   ze1e2w, zcof, zbe1ru, zbe2rv, zbtr, ztz, zth   !    "      scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zkepe, zkx, zky, zkz   ! temporary arrays
-      !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   zkepe => wrk_3d_1 , zkx => wrk_3d_2   ! 3D workspace
+      USE wrk_nemo, ONLY:   zky   => wrk_3d_3 , zkz => wrk_3d_4   !  -      -
+      !
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zcof         ! local scalar
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(3, 1,2,3,4) ) THEN
+         CALL ctl_stop('trd_dwr: requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
 
       ! I. Momentum trends
@@ -359,63 +273,45 @@
          ! --------------------------------------------------
          ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
-
-         zkx(:,:,:)   = 0.e0
-         zky(:,:,:)   = 0.e0
-         zkz(:,:,:)   = 0.e0
-         zkepe(:,:,:) = 0.e0
+         zkx  (:,:,:) = 0._wp
+         zky  (:,:,:) = 0._wp
+         zkz  (:,:,:) = 0._wp
+         zkepe(:,:,:) = 0._wp
    
          CALL eos( tsn, rhd, rhop )       ! now potential and in situ densities
 
-         ! Density flux at w-point
+         zcof = 0.5_wp / rau0             ! Density flux at w-point
+         zkz(:,:,1) = 0._wp
          DO jk = 2, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ze1e2w = 0.5 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) * tmask_i(ji,jj)
-                  zkz(ji,jj,jk) = ze1e2w / rau0 * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) )
+            zkz(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
+         END DO
+         
+         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zkx(ji,jj,jk) = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
+                  zky(ji,jj,jk) = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
                END DO
             END DO
          END DO
-         zkz(:,:,1) = 0.e0
          
-         ! Density flux at u and v-points
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zcof   = 0.5 / rau0
-                  zbe1ru = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
-                  zbe2rv = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
-                  zkx(ji,jj,jk) = zbe1ru * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
-                  zky(ji,jj,jk) = zbe2rv * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
+         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
+                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
+                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
+                     &              / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
          END DO
-         
-         ! Density flux divergence at t-point
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
-                  ztz = - zbtr * (    zkz(ji,jj,jk) - zkz(ji,jj,jk+1) )
-                  zth = - zbtr * (  ( zkx(ji,jj,jk) - zkx(ji-1,jj,jk) )   &
-                    &             + ( zky(ji,jj,jk) - zky(ji,jj-1,jk) )  )
-                  zkepe(ji,jj,jk) = (zth + ztz) * tmask(ji,jj,jk) * tmask_i(ji,jj)
-               END DO
-            END DO
-         END DO
-         zkepe( : , : ,jpk) = 0.e0
-         zkepe( : ,jpj, : ) = 0.e0
-         zkepe(jpi, : , : ) = 0.e0
 
          ! I.2 Basin averaged kinetic energy trend
          ! ----------------------------------------
-         peke = 0.e0
-         DO jk = 1,jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  peke = peke + zkepe(ji,jj,jk) * grav * fsdept(ji,jj,jk)   &
-                     &                     * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         peke = 0._wp
+         DO jk = 1, jpkm1
+            peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
+         END DO
+         peke = grav * peke
 
          ! I.3 Sums over the global domain
@@ -542,4 +438,6 @@
          !
       ENDIF
+      !
+      IF( wrk_not_released(3, 1,2,3,4) )   CALL ctl_stop('trd_dwr: failed to release workspace arrays')
       !
    END SUBROUTINE trd_dwr
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90	(revision 2715)
@@ -15,5 +15,4 @@
    !!   'key_trdmld'                          mixed layer trend diagnostics
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
    !!   trd_mld          : T and S cumulated trends averaged over the mixed layer
    !!   trd_mld_zint     : T and S trends vertical integration
@@ -23,4 +22,5 @@
    USE dom_oce         ! ocean space and time domain variables
    USE trdmod_oce      ! ocean variables trends
+   USE trdmld_oce      ! ocean variables trends
    USE ldftra_oce      ! ocean active tracers lateral physics
    USE zdf_oce         ! ocean vertical physics
@@ -37,4 +37,5 @@
    USE prtctl          ! Print control
    USE restart         ! for lrst_oce
+   USE lib_mpp         ! MPP library
 
    IMPLICIT NONE
@@ -47,5 +48,6 @@
    CHARACTER (LEN=40) ::  clhstnam         ! name of the trends NetCDF file
    INTEGER ::   nh_t, nmoymltrd
-   INTEGER ::   nidtrd, ndextrd1(jpi*jpj)
+   INTEGER ::   nidtrd
+   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::   ndextrd1
    INTEGER ::   ndimtrd1                        
    INTEGER ::   ionce, icount                   
@@ -58,8 +60,18 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
+
+   INTEGER FUNCTION trd_mld_alloc()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mld_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( ndextrd1(jpi*jpj) , STAT=trd_mld_alloc )
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( trd_mld_alloc )
+      IF( trd_mld_alloc /= 0 )   CALL ctl_warn('trd_mld_alloc: failed to allocate array ndextrd1')
+   END FUNCTION trd_mld_alloc
+
 
    SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )
@@ -81,12 +93,18 @@
       !!            surface and the control surface is called "mixed-layer"
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   ktrd                             ! ocean trend index
-      CHARACTER(len=2), INTENT( in ) :: ctype                     ! surface/bottom (2D arrays) or
-      !                                                           ! interior (3D arrays) physics
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::  pttrdmld ! temperature trend 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::  pstrdmld ! salinity trend 
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   zvlmsk => wrk_2d_10     ! 2D workspace
+      !
+      INTEGER                         , INTENT( in ) ::   ktrd       ! ocean trend index
+      CHARACTER(len=2)                , INTENT( in ) ::   ctype      ! 2D surface/bottom or 3D interior physics
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   pttrdmld   ! temperature trend 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   pstrdmld   ! salinity trend 
+      !
       INTEGER ::   ji, jj, jk, isum
-      REAL(wp), DIMENSION(jpi,jpj) ::  zvlmsk
-      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 10) ) THEN
+         CALL ctl_stop('trd_mld_zint : requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
 
       ! I. Definition of control surface and associated fields
@@ -176,4 +194,6 @@
          smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,1) * wkx(:,:,1)            
       END SELECT
+      !
+      IF( wrk_not_released(2, 10) )   CALL ctl_stop('trd_mld_zint: failed to release workspace arrays')
       !
    END SUBROUTINE trd_mld_zint
@@ -227,18 +247,19 @@
       !!       - See NEMO documentation (in preparation)
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY: ztmltot => wrk_2d_1,  zsmltot => wrk_2d_2 ! dT/dt over the anlysis window (including Asselin)
+      USE wrk_nemo, ONLY: ztmlres => wrk_2d_3,  zsmlres => wrk_2d_4 ! residual = dh/dt entrainment term
+      USE wrk_nemo, ONLY: ztmlatf => wrk_2d_5,  zsmlatf => wrk_2d_6 ! needed for storage only
+      USE wrk_nemo, ONLY: ztmltot2 => wrk_2d_7, ztmlres2 => wrk_2d_8, ztmltrdm2 => wrk_2d_9    ! \  working arrays to diagnose the trends
+      USE wrk_nemo, ONLY: zsmltot2 => wrk_2d_10, zsmlres2 => wrk_2d_11, zsmltrdm2 => wrk_2d_12 !  > associated with the time meaned ML T & S
+      USE wrk_nemo, ONLY: ztmlatf2 => wrk_2d_13, zsmlatf2 => wrk_2d_14    
+      USE wrk_nemo, ONLY: wrk_3d_1, wrk_3d_2                     ! /
+      !
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      !!
+      !
       INTEGER :: ji, jj, jk, jl, ik, it, itmod
       LOGICAL :: lldebug = .TRUE.
       REAL(wp) :: zavt, zfn, zfn2
-      REAL(wp) ,DIMENSION(jpi,jpj) ::     &
-           ztmltot,  zsmltot,             & ! dT/dt over the anlysis window (including Asselin)
-           ztmlres,  zsmlres,             & ! residual = dh/dt entrainment term
-           ztmlatf,  zsmlatf,             & ! needed for storage only
-           ztmltot2, ztmlres2, ztmltrdm2, & ! \  working arrays to diagnose the trends
-           zsmltot2, zsmlres2, zsmltrdm2, & !  > associated with the time meaned ML T & S
-           ztmlatf2, zsmlatf2               ! /
-      REAL(wp), DIMENSION(jpi,jpj,jpltrd) ::  &
-           ztmltrd2, zsmltrd2               ! only needed for mean diagnostics
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztmltrd2, zsmltrd2   ! only needed for mean diagnostics
 #if defined key_dimgout
       INTEGER ::  iyear,imon,iday
@@ -247,4 +268,16 @@
       !!----------------------------------------------------------------------
       
+      ! Check that the workspace arrays are all OK to be used
+      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14)  .OR. &
+          wrk_in_use(3, 1,2)                                 ) THEN
+         CALL ctl_stop('trd_mld : requested workspace arrays unavailable')   ;   RETURN
+      ELSE IF(jpltrd > jpk) THEN
+         ! ARPDBG, is this reasonable or will this cause trouble in the future?
+         CALL ctl_stop('trd_mld : no. of mixed-layer trends (jpltrd) exceeds no. of model levels so cannot use 3D workspaces.')
+         RETURN         
+      END IF
+      ! Set-up pointers into sub-arrays of 3d-workspaces
+      ztmltrd2 => wrk_3d_1(1:,:,1:jpltrd)
+      zsmltrd2 => wrk_3d_2(1:,:,1:jpltrd)
 
       ! ======================================================================
@@ -707,4 +740,8 @@
       IF( lrst_oce )   CALL trd_mld_rst_write( kt ) 
 
+      IF( wrk_not_released(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14)  .OR. &
+          wrk_not_released(3, 1,2)                                )   &
+          CALL ctl_stop('trd_mld : failed to release workspace arrays.')
+      !
    END SUBROUTINE trd_mld
 
@@ -716,15 +753,10 @@
       !! ** Purpose :   computation of vertically integrated T and S budgets
       !!      from ocean surface down to control surface (NetCDF output)
-      !!
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       INTEGER :: jl
       INTEGER :: inum   ! logical unit
-
       REAL(wp) ::   zjulian, zsto, zout
-
       CHARACTER (LEN=40) ::   clop
       CHARACTER (LEN=12) ::   clmxl, cltu, clsu
-
       !!----------------------------------------------------------------------
 
@@ -763,4 +795,8 @@
          nwarn = nwarn + 1
       END IF
+
+      !                                   ! allocate trdmld arrays
+      IF( trd_mld_alloc()    /= 0 )   CALL ctl_stop( 'STOP', 'trd_mld_init : unable to allocate trdmld     arrays' )
+      IF( trdmld_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_mld_init : unable to allocate trdmld_oce arrays' )
 
       ! I.2 Initialize arrays to zero or read a restart file
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld_oce.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld_oce.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld_oce.F90	(revision 2715)
@@ -4,10 +4,12 @@
    !! Ocean trends :   set tracer and momentum trend variables
    !!======================================================================
-   !! History :  9.0  !  04-08  (C. Talandier)  New trends organization
+   !! History :  1.0  ! 2004-08  (C. Talandier)  New trends organization
    !!----------------------------------------------------------------------
-   USE par_oce         ! ocean parameters
+   USE par_oce        ! ocean parameters
 
    IMPLICIT NONE
    PRIVATE
+
+   PUBLIC   trdmld_oce_alloc    ! Called in trdmld.F90
 
 #if defined key_trdmld
@@ -17,10 +19,6 @@
 #endif
    !!* mixed layer trends indices
-   INTEGER, PARAMETER, PUBLIC ::   jpltrd = 11    !: number of mixed-layer trends arrays
-   INTEGER, PUBLIC   & 
-#if !defined key_agrif
-      , PARAMETER  &
-#endif
-::   jpktrd = jpk   !: max level for mixed-layer trends diag.
+   INTEGER, PARAMETER, PUBLIC ::   jpltrd = 11      !: number of mixed-layer trends arrays
+   INTEGER, PUBLIC            ::   jpktrd           !: max level for mixed-layer trends diag.
    !
    INTEGER, PUBLIC, PARAMETER ::   jpmld_xad =  1   !:  zonal      
@@ -46,10 +44,10 @@
    CHARACTER(LEN=80) , PUBLIC :: clname, ctrd(jpltrd+1,2)
 
-   INTEGER , PUBLIC, DIMENSION(jpi,jpj)     ::   nmld   !: mixed layer depth indexes 
-   INTEGER , PUBLIC, DIMENSION(jpi,jpj)     ::   nbol   !: mixed-layer depth indexes when read from file
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nmld   !: mixed layer depth indexes 
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nbol   !: mixed-layer depth indexes when read from file
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   wkx    !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wkx    !:
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  &
       rmld   ,                      & !: mld depth (m) corresponding to nmld
       tml    , sml  ,               & !: \ "now" mixed layer temperature/salinity
@@ -66,5 +64,5 @@
       rmld_sum, rmldbn                !: needed to compute the leap-frog time mean of the ML depth
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  &
       tmlatfb, tmlatfn ,            & !: "before" Asselin contribution at begining of the averaging
       smlatfb, smlatfn,             & !: period (i.e. last contrib. from previous such period) and 
@@ -72,5 +70,5 @@
       tmlatfm, smlatfm                !: accumulator for Asselin trends (needed for storage only)
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpltrd) ::  &
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::  &
       tmltrd,                       & !: \ physical contributions to the total trend (for T/S),
       smltrd,                       & !: / cumulated over the current analysis window
@@ -83,7 +81,56 @@
 #endif
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+  INTEGER FUNCTION trdmld_oce_alloc()
+     !!----------------------------------------------------------------------
+     !!                 ***  FUNCTION trdmld_oce_alloc   ***
+     !!----------------------------------------------------------------------
+     USE lib_mpp
+     INTEGER :: ierr(5)
+     !!----------------------------------------------------------------------
+
+     ! Initialise jpktrd here as can no longer do it in MODULE body since
+     ! jpk is now a variable.
+     jpktrd = jpk   !: max level for mixed-layer trends diag.
+
+     ierr(:) = 0
+
+#if   defined  key_trdmld   ||   defined key_esopa
+     ALLOCATE( nmld(jpi,jpj), nbol(jpi,jpj),       &
+        &      wkx(jpi,jpj,jpk), rmld(jpi,jpj),    & 
+        &      tml(jpi,jpj)    , sml(jpi,jpj),     & 
+        &      tmlb(jpi,jpj)   , smlb(jpi,jpj) ,   &
+        &      tmlbb(jpi,jpj)  , smlbb(jpi,jpj), STAT = ierr(1) )
+
+     ALLOCATE( tmlbn(jpi,jpj)  , smlbn(jpi,jpj),   &
+        &      tmltrdm(jpi,jpj), smltrdm(jpi,jpj), &
+        &      tml_sum(jpi,jpj), tml_sumb(jpi,jpj),&
+        &      tmltrd_atf_sumb(jpi,jpj)           , STAT=ierr(2) )
+
+     ALLOCATE( sml_sum(jpi,jpj), sml_sumb(jpi,jpj), &
+        &      smltrd_atf_sumb(jpi,jpj),            &
+        &      rmld_sum(jpi,jpj), rmldbn(jpi,jpj),  &
+        &      tmlatfb(jpi,jpj), tmlatfn(jpi,jpj), STAT = ierr(3) )
+
+     ALLOCATE( smlatfb(jpi,jpj), smlatfn(jpi,jpj), & 
+        &      tmlatfm(jpi,jpj), smlatfm(jpi,jpj), &
+        &      tmltrd(jpi,jpj,jpltrd),   smltrd(jpi,jpj,jpltrd), STAT=ierr(4))
+
+     ALLOCATE( tmltrd_sum(jpi,jpj,jpltrd),tmltrd_csum_ln(jpi,jpj,jpltrd),      &
+        &      tmltrd_csum_ub(jpi,jpj,jpltrd), smltrd_sum(jpi,jpj,jpltrd),     &
+        &      smltrd_csum_ln(jpi,jpj,jpltrd), smltrd_csum_ub(jpi,jpj,jpltrd), STAT=ierr(5) )
+#endif
+      !
+      trdmld_oce_alloc = MAXVAL( ierr )
+      IF( lk_mpp                )   CALL mpp_sum ( trdmld_oce_alloc )
+      IF( trdmld_oce_alloc /= 0 )   CALL ctl_warn('trdmld_oce_alloc: failed to allocate arrays')
+      !
+   END FUNCTION trdmld_oce_alloc
+
    !!======================================================================
 END MODULE trdmld_oce
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 2715)
@@ -24,4 +24,5 @@
    USE trdmld                  ! ocean active mixed layer tracers trends 
    USE in_out_manager          ! I/O manager
+   USE lib_mpp         ! MPP library
 
    IMPLICIT NONE
@@ -39,5 +40,5 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -51,20 +52,29 @@
       !!              integral constraints
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt                                ! time step
-      INTEGER, INTENT( in ) ::   ktrd                              ! tracer trend index
-      CHARACTER(len=3), INTENT( in ) ::   ctype                    ! momentum or tracers trends type 'DYN'/'TRA'
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   ptrdx ! Temperature or U trend 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   ptrdy ! Salinity    or V trend
+      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY: ztswu => wrk_2d_1,  &
+                          ztswv => wrk_2d_2,  &
+                          ztbfu => wrk_2d_3,  &
+                          ztbfv => wrk_2d_4,  &
+                          z2dx  => wrk_2d_5,  &
+                          z2dy  => wrk_2d_6
+      !
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
+      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type 'DYN'/'TRA'
+      INTEGER                   , INTENT(in   ) ::   kt      ! time step
+      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
       !!
-      INTEGER ::   ji, jj
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztswu, ztswv               ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztbfu, ztbfv               ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2dx, z2dy                 ! workspace arrays
-      !!----------------------------------------------------------------------
-
-      z2dx(:,:)   = 0.e0   ;   z2dy(:,:)   = 0.e0                  ! initialization of workspace arrays
-
-      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restarting with Euler time stepping)
-      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
+      INTEGER ::   ji, jj   ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF(wrk_in_use(2, 1,2,3,4,5,6))THEN
+         CALL ctl_warn('trd_mod: Requested workspace arrays already in use.')   ;   RETURN
+      END IF
+
+      z2dx(:,:) = 0._wp   ;   z2dy(:,:) = 0._wp                            ! initialization of workspace arrays
+
+      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt   ! = rdtra (restart with Euler time stepping)
+      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt   ! = 2 rdttra (leapfrog)
       ENDIF
 
@@ -84,15 +94,15 @@
             CASE ( jptra_trd_dmp )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_dmp, ctype )   ! damping
             CASE ( jptra_trd_qsr )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_qsr, ctype )   ! penetrative solar radiat.
-            CASE ( jptra_trd_nsr )   
-               z2dx(:,:) = ptrdx(:,:,1)   ;   z2dy(:,:) = ptrdy(:,:,1)
-               CALL trd_icp( z2dx, z2dy, jpicpt_nsr, ctype )                               ! non solar radiation
+            CASE ( jptra_trd_nsr )   ;   z2dx(:,:) = ptrdx(:,:,1)   
+                                         z2dy(:,:) = ptrdy(:,:,1)
+                                         CALL trd_icp( z2dx , z2dy , jpicpt_nsr, ctype )   ! non solar radiation
             CASE ( jptra_trd_xad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_xad, ctype )   ! x- horiz adv
             CASE ( jptra_trd_yad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_yad, ctype )   ! y- horiz adv
-            CASE ( jptra_trd_zad )                                                         ! z- vertical adv 
-               CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype )   
-               ! compute the surface flux condition wn(:,:,1)*tn(:,:,1)
-               z2dx(:,:) = wn(:,:,1)*tn(:,:,1)/fse3t(:,:,1)
-               z2dy(:,:) = wn(:,:,1)*sn(:,:,1)/fse3t(:,:,1)
-               CALL trd_icp( z2dx , z2dy , jpicpt_zl1, ctype )                             ! 1st z- vertical adv 
+            CASE ( jptra_trd_zad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype )   ! z- vertical adv 
+                                         CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype )   
+                                         ! compute the surface flux condition wn(:,:,1)*tn(:,:,1)
+                                         z2dx(:,:) = wn(:,:,1)*tn(:,:,1)/fse3t(:,:,1)
+                                         z2dy(:,:) = wn(:,:,1)*sn(:,:,1)/fse3t(:,:,1)
+                                         CALL trd_icp( z2dx , z2dy , jpicpt_zl1, ctype )   ! 1st z- vertical adv 
             END SELECT
          END IF
@@ -113,6 +123,6 @@
                ! subtract surface forcing/bottom friction trends 
                ! from vertical diffusive momentum trends
-               ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
-               ztbfu(:,:) = 0.e0   ;   ztbfv(:,:) = 0.e0  
+               ztswu(:,:) = 0._wp   ;   ztswv(:,:) = 0._wp
+               ztbfu(:,:) = 0._wp   ;   ztbfv(:,:) = 0._wp 
                DO jj = 2, jpjm1   
                   DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -121,7 +131,6 @@
                      ztswv(ji,jj) = vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
                      ! bottom friction contribution now handled explicitly
-                     !
-                     ptrdx(ji,jj,1     ) = ptrdx(ji,jj,1     ) - ztswu(ji,jj)
-                     ptrdy(ji,jj,1     ) = ptrdy(ji,jj,1     ) - ztswv(ji,jj)
+                     ptrdx(ji,jj,1) = ptrdx(ji,jj,1) - ztswu(ji,jj)
+                     ptrdy(ji,jj,1) = ptrdy(ji,jj,1) - ztswv(ji,jj)
                   END DO
                END DO
@@ -218,4 +227,6 @@
       ENDIF
       !
+      IF( wrk_not_released(2, 1,2,3,4,5,6) )   CALL ctl_warn('trd_mod: Failed to release workspace arrays.')
+      !
    END SUBROUTINE trd_mod
 
@@ -228,9 +239,9 @@
    USE trdicp          ! ocean bassin integral constraints properties
    USE trdmld          ! ocean active mixed layer tracers trends 
-
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE trd_mod(ptrd3dx, ptrd3dy, ktrd , ctype, kt )   ! Empty routine
-      REAL    ::   ptrd3dx(:,:,:), ptrd3dy(:,:,:)
-      INTEGER ::   ktrd, kt                            
+      REAL(wp) ::   ptrd3dx(:,:,:), ptrd3dy(:,:,:)
+      INTEGER  ::   ktrd, kt                            
       CHARACTER(len=3) ::  ctype                  
       WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd3dx(1,1,1), ptrd3dy(1,1,1)
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod_oce.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod_oce.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod_oce.F90	(revision 2715)
@@ -75,5 +75,5 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE trdmod_oce
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod_trc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod_trc.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod_trc.F90	(revision 2715)
@@ -4,9 +4,4 @@
    !!  Dummy module
    !!======================================================================
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   Dummy module                                             NO TOP use
@@ -21,4 +16,8 @@
    END SUBROUTINE trd_mod_trc
 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE trdmod_trc
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 2715)
@@ -4,7 +4,7 @@
    !! Ocean diagnostics:  ocean tracers trends
    !!=====================================================================
-   !! History :  9.0  !  2004-08  (C. Talandier) Original code
-   !!                 !  2005-04  (C. Deltel)    Add Asselin trend in the ML budget
-   !!      NEMO  3.3  !  2010-06  (C. Ethe) merge TRA-TRC 
+   !! History :  1.0  !  2004-08  (C. Talandier) Original code
+   !!            2.0  !  2005-04  (C. Deltel)    Add Asselin trend in the ML budget
+   !!            3.3  !  2010-06  (C. Ethe) merge TRA-TRC 
    !!----------------------------------------------------------------------
 #if  defined key_trdtra || defined key_trdmld || defined key_trdmld_trc 
@@ -12,16 +12,17 @@
    !!   trd_tra      : Call the trend to be computed
    !!----------------------------------------------------------------------
-   USE dom_oce            ! ocean domain 
-   USE trdmod_oce         ! ocean active mixed layer tracers trends 
-   USE trdmod             ! ocean active mixed layer tracers trends 
-   USE trdmod_trc         ! ocean passive mixed layer tracers trends 
+   USE dom_oce          ! ocean domain 
+   USE trdmod_oce       ! ocean active mixed layer tracers trends 
+   USE trdmod           ! ocean active mixed layer tracers trends 
+   USE trdmod_trc       ! ocean passive mixed layer tracers trends 
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp          ! MPP library
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC trd_tra          ! called by all  traXX modules
+   PUBLIC   trd_tra          ! called by all  traXX modules
  
-   !! * Module declaration
-   REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: trdtx, trdty, trdt  !:
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: trdtx, trdty, trdt  !:
 
    !! * Substitutions
@@ -29,10 +30,20 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
+
+   INTEGER FUNCTION trd_tra_alloc()
+      !!----------------------------------------------------------------------------
+      !!                  ***  FUNCTION trd_tra_alloc  ***
+      !!----------------------------------------------------------------------------
+      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , STAT= trd_tra_alloc )
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( trd_tra_alloc )
+      IF( trd_tra_alloc /= 0 )   CALL ctl_warn('trd_tra_alloc: failed to allocate arrays')
+   END FUNCTION trd_tra_alloc
+
 
    SUBROUTINE trd_tra( kt, ctype, ktra, ktrd, ptrd, pun, ptra )
@@ -50,4 +61,7 @@
       !!        nn_ctls > 1  : use fixed level surface jk = nn_ctls
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   ztrds => wrk_3d_10   ! 3D workspace
+      !
       INTEGER                         , INTENT(in)           ::  kt      ! time step
       CHARACTER(len=3)                , INTENT(in)           ::  ctype   ! tracers trends type 'TRA'/'TRC'
@@ -57,8 +71,14 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::  pun     ! velocity 
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::  ptra    ! Tracer variable 
-      !!
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ztrds    !  
-      !!----------------------------------------------------------------------
-
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(3, 10) ) THEN
+         CALL ctl_stop('trd_tra: requested workspace array unavailable')   ;   RETURN
+      ENDIF
+
+      IF( .NOT. ALLOCATED( trdtx ) ) THEN       ! allocate trdtra arrays
+         IF( trd_tra_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' )
+      ENDIF
+      
       ! Control of optional arguments
       IF( ctype == 'TRA' .AND. ktra == jp_tem ) THEN 
@@ -118,5 +138,8 @@
       ENDIF
       !
+      IF( wrk_not_released(3, 10) )   CALL ctl_stop('trd_tra: failed to release workspace array')
+      !
    END SUBROUTINE trd_tra
+
 
    SUBROUTINE trd_tra_adv( pf, pun, ptn, cdir, ptrd )
@@ -130,13 +153,13 @@
       !!                k-advective trends = -un. di+1[T] = -( dk+1[fi] - tn dk+1[un] )
       !!----------------------------------------------------------------------
-      REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk)           ::   pf      ! advective flux in one direction
-      REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk)           ::   pun     ! now velocity  in one direction
-      REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk)           ::   ptn     ! now or before tracer 
-      CHARACTER(len=1), INTENT(in )                                   ::   cdir    ! X/Y/Z direction
-      REAL(wp)        , INTENT(out), DIMENSION(jpi,jpj,jpk)           ::   ptrd    ! advective trend in one direction
-      !!
-      INTEGER                          ::   ji, jj, jk   ! dummy loop indices
-      INTEGER                          ::   ii, ij, ik   ! index shift function of the direction
-      REAL(wp)                         ::   zbtr         ! temporary scalar
+      REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk) ::   pf      ! advective flux in one direction
+      REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk) ::   pun     ! now velocity  in one direction
+      REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk) ::   ptn     ! now or before tracer 
+      CHARACTER(len=1), INTENT(in )                         ::   cdir    ! X/Y/Z direction
+      REAL(wp)        , INTENT(out), DIMENSION(jpi,jpj,jpk) ::   ptrd    ! advective trend in one direction
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ii, ij, ik   ! index shift function of the direction
+      REAL(wp) ::   zbtr         ! local scalar
       !!----------------------------------------------------------------------
 
@@ -167,10 +190,8 @@
 #   else
    !!----------------------------------------------------------------------
-   !!   Default case :                                         Empty module
+   !!   Default case :          Dummy module           No trend diagnostics
    !!----------------------------------------------------------------------
    USE par_oce      ! ocean variables trends
-
 CONTAINS
-
    SUBROUTINE trd_tra( kt, ctype, ktra, ktrd, ptrd, pu, ptra )
       !!----------------------------------------------------------------------
@@ -182,11 +203,6 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::  pu      ! velocity 
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::  ptra    ! Tracer variable 
-      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd(1,1,1)
-      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ptra(1,1,1)
-      WRITE(*,*) ' "   ": You should not have seen this print! error ?', pu(1,1,1)
-      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ktrd
-      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ktra
-      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ctype
-      WRITE(*,*) ' "   ": You should not have seen this print! error ?', kt
+      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd(1,1,1), ptra(1,1,1), pu(1,1,1),   &
+         &                                                               ktrd, ktra, ctype, kt
    END SUBROUTINE trd_tra
 #   endif
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90	(revision 2715)
@@ -4,6 +4,6 @@
    !! Ocean diagnostics:  momentum trends
    !!=====================================================================
-   !! History :  9.0  !  04-06  (L. Brunier, A-M. Treguier) Original code 
-   !!                 !  04-08  (C. Talandier) New trends organization
+   !! History :  1.0  !  04-2006  (L. Brunier, A-M. Treguier) Original code 
+   !!            2.0  !  04-2008  (C. Talandier) New trends organization
    !!----------------------------------------------------------------------
 #if defined key_trdvor   ||   defined key_esopa
@@ -26,4 +26,5 @@
    USE ioipsl          ! NetCDF library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! MPP library
 
    IMPLICIT NONE
@@ -37,17 +38,18 @@
    PUBLIC   trd_vor_zint   ! routine called by dynamics routines
    PUBLIC   trd_vor_init   ! routine called by opa.F90
-
-   INTEGER ::   nh_t, nmoydpvor, nidvor, nhoridvor, ndexvor1(jpi*jpj), ndimvor1, icount   ! needs for IOIPSL output
+   PUBLIC   trd_vor_alloc  ! routine called by nemogcm.F90
+
+   INTEGER ::   nh_t, nmoydpvor, nidvor, nhoridvor, ndimvor1, icount   ! needs for IOIPSL output
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) ::   ndexvor1   ! needed for IOIPSL output
    INTEGER ::   ndebug     ! (0/1) set it to 1 in case of problem to have more print
 
-   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avr      ! average
-   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrb     ! before vorticity (kt-1)
-   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrbb    ! vorticity at begining of the nwrite-1 timestep averaging period
-   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrbn    ! after vorticity at time step after the
-   REAL(wp), DIMENSION(jpi,jpj) ::   rotot        ! begining of the NWRITE-1 timesteps
-   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrtot   !
-   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrres   !
-
-   REAL(wp), DIMENSION(jpi,jpj,jpltot_vor) ::   vortrd  ! curl of trends
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   vor_avr      ! average
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   vor_avrb     ! before vorticity (kt-1)
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   vor_avrbb    ! vorticity at begining of the nwrite-1 timestep averaging period
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   vor_avrbn    ! after vorticity at time step after the
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   rotot        ! begining of the NWRITE-1 timesteps
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   vor_avrtot   !
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   vor_avrres   !
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   vortrd       ! curl of trends
          
    CHARACTER(len=12) ::   cvort
@@ -63,4 +65,18 @@
    !!----------------------------------------------------------------------
 CONTAINS
+
+   INTEGER FUNCTION trd_vor_alloc()
+      !!----------------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_vor_alloc  ***
+      !!----------------------------------------------------------------------------
+      ALLOCATE( vor_avr   (jpi,jpj) , vor_avrb(jpi,jpj) , vor_avrbb (jpi,jpj) ,   &
+         &      vor_avrbn (jpi,jpj) , rotot   (jpi,jpj) , vor_avrtot(jpi,jpj) ,   &
+         &      vor_avrres(jpi,jpj) , vortrd  (jpi,jpj,jpltot_vor) ,              &
+         &      ndexvor1  (jpi*jpj)                                ,   STAT= trd_vor_alloc )
+         !
+      IF( lk_mpp             )   CALL mpp_sum ( trd_vor_alloc )
+      IF( trd_vor_alloc /= 0 )   CALL ctl_warn('trd_vor_alloc: failed to allocate arrays')
+   END FUNCTION trd_vor_alloc
+
 
    SUBROUTINE trd_vor_zint_2d( putrdvor, pvtrdvor, ktrd )
@@ -91,19 +107,23 @@
       !!      trends output in netCDF format using ioipsl
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   zudpvor => wrk_2d_1 , zvdpvor => wrk_2d_2   ! total cmulative trends
+      !
       INTEGER                     , INTENT(in   ) ::   ktrd       ! ocean trend index
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   putrdvor   ! u vorticity trend 
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pvtrdvor   ! v vorticity trend
-      !!
+      !
       INTEGER ::   ji, jj       ! dummy loop indices
       INTEGER ::   ikbu, ikbv   ! local integers
-      REAL(wp), DIMENSION(jpi,jpj) ::   zudpvor, zvdpvor   ! total cmulative trends
-      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 1,2) ) THEN
+         CALL ctl_stop('trd_vor_zint_2d: requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
 
       ! Initialization
-      zudpvor(:,:) = 0._wp
-      zvdpvor(:,:) = 0._wp
-      !
-      CALL lbc_lnk( putrdvor,  'U' , -1. )         ! lateral boundary condition on input momentum trends
-      CALL lbc_lnk( pvtrdvor,  'V' , -1. )
+      zudpvor(:,:) = 0._wp                 ;   zvdpvor(:,:) = 0._wp
+      CALL lbc_lnk( putrdvor, 'U', -1. )   ;   CALL lbc_lnk( pvtrdvor, 'V', -1. )      ! lateral boundary condition
+      
 
       !  =====================================
@@ -147,4 +167,6 @@
       ENDIF
       !
+      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('trd_vor_zint_2d : failed to release workspace arrays.')
+      !
    END SUBROUTINE trd_vor_zint_2d
 
@@ -177,13 +199,19 @@
       !!      trends output in netCDF format using ioipsl
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   zubet   => wrk_2d_1,   zvbet => wrk_2d_2   ! Beta.V 
+      USE wrk_nemo, ONLY:   zudpvor => wrk_2d_3, zvdpvor => wrk_2d_4   ! total cmulative trends
+      !
       INTEGER                         , INTENT(in   ) ::   ktrd       ! ocean trend index
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   putrdvor   ! u vorticity trend 
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvtrdvor   ! v vorticity trend
-      !!
-      INTEGER ::   ji, jj, jk
-      REAL(wp), DIMENSION(jpi,jpj) ::   zubet  , zvbet     ! Beta.V 
-      REAL(wp), DIMENSION(jpi,jpj) ::   zudpvor, zvdpvor   ! total cmulative trends
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
       !!----------------------------------------------------------------------
      
+      IF( wrk_in_use(2, 1,2,3,4) ) THEN
+         CALL ctl_stop('trd_vor_zint_3d: requested workspace arrays unavailable.')   ;   RETURN
+      ENDIF
+
       ! Initialization
       zubet  (:,:) = 0._wp
@@ -192,6 +220,6 @@
       zvdpvor(:,:) = 0._wp
       !
-      CALL lbc_lnk( putrdvor, 'U' , -1. )         ! lateral boundary condition on input momentum trends
-      CALL lbc_lnk( pvtrdvor, 'V' , -1. )
+      CALL lbc_lnk( putrdvor, 'U', -1. )         ! lateral boundary condition on input momentum trends
+      CALL lbc_lnk( pvtrdvor, 'V', -1. )
 
       !  =====================================
@@ -248,4 +276,6 @@
       ENDIF
       !
+      IF( wrk_not_released(2, 1,2,3,4) )   CALL ctl_stop('trd_vor_zint_3d: failed to release workspace arrays')
+      !
    END SUBROUTINE trd_vor_zint_3d
 
@@ -258,11 +288,17 @@
       !!               and make outputs (NetCDF or DIMG format)
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   zun => wrk_2d_1 , zvn => wrk_2d_2 ! 2D workspace
+      !
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !!
+      !
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       INTEGER  ::   it, itmod        ! local integers
       REAL(wp) ::   zmean            ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zun, zvn   ! 2D workspace
-      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 1,2) ) THEN
+         CALL ctl_stop('trd_vor: requested workspace arrays unavailable.')   ;   RETURN
+      ENDIF
 
       !  =================
@@ -431,4 +467,6 @@
       IF( kt == nitend )   CALL histclo( nidvor )
       !
+      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('trd_vor: failed to release workspace arrays')
+      !
    END SUBROUTINE trd_vor
 
@@ -466,4 +504,7 @@
          WRITE(numout,*) ' '
       ENDIF
+
+      IF( trd_vor_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_vor_init : unable to allocate trdvor arrays' )
+
 
       ! cumulated trends array init
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor_oce.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor_oce.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor_oce.F90	(revision 2715)
@@ -4,5 +4,5 @@
    !! Ocean trends :   set vorticity trend variables
    !!======================================================================
-   !! History :  9.0  ! ???
+   !! History :  9.0  !  04-2006  (L. Brunier, A-M. Treguier) Original code 
    !!----------------------------------------------------------------------
 
@@ -14,9 +14,9 @@
 
 #if defined key_trdvor
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trdvor = .TRUE.     !: momentum trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdvor = .TRUE.    !: momentum trend flag
 #else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trdvor = .FALSE.    !: momentum trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdvor = .FALSE.   !: momentum trend flag
 #endif
-   !!* vorticity trends index
+   !                                               !!* vorticity trends index
    INTEGER, PUBLIC, PARAMETER ::   jpltot_vor = 11  !: Number of vorticity trend terms
    !
@@ -36,5 +36,5 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE trdvor_oce
