Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 2715)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 3294)
@@ -29,4 +29,6 @@
    USE lib_mpp          ! MPP library
    USE lbclnk           ! lateral boundary condition - processor exchanges
+   USE timing           ! preformance summary
+   USE wrk_nemo         ! working arrays
 
    IMPLICIT NONE
@@ -95,5 +97,5 @@
       !!----------------------------------------------------------------------
       INTEGER               ::   dia_ptr_alloc   ! return value
-      INTEGER, DIMENSION(5) ::   ierr
+      INTEGER, DIMENSION(6) ::   ierr
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -121,4 +123,6 @@
          &     ndex_h_ind_30(jpj),   ndex_h_ipc_30(jpj), Stat=ierr(5) )
          !
+     ALLOCATE( btm30(jpi,jpj) , STAT=ierr(6)  )
+         !
       dia_ptr_alloc = MAXVAL( ierr )
       IF(lk_mpp)   CALL mpp_sum( dia_ptr_alloc )
@@ -209,8 +213,4 @@
       !! ** Action  : - p_fval: i-mean poleward flux of pva
       !!----------------------------------------------------------------------
-#if defined key_mpp_mpi
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zwork => wrk_1d_1
-#endif
       !!
       IMPLICIT none
@@ -225,10 +225,12 @@
       INTEGER               ::   ijpjjpk
 #endif
+#if defined key_mpp_mpi
+      REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point
+#endif
       !!--------------------------------------------------------------------
       !
 #if defined key_mpp_mpi
-      IF( wrk_in_use(1, 1) ) THEN
-         CALL ctl_stop('ptr_vjk: ERROR - requested workspace array is unavailable')   ;   RETURN
-      END IF
+      ijpjjpk = jpj*jpk
+      CALL wrk_alloc( jpj*jpk, zwork )
 #endif
 
@@ -257,5 +259,4 @@
       !
 #if defined key_mpp_mpi
-      ijpjjpk = jpj*jpk
       ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk
       zwork(1:ijpjjpk) = RESHAPE( p_fval, ish )
@@ -265,5 +266,5 @@
       !
 #if defined key_mpp_mpi
-      IF( wrk_not_released(1, 1) )   CALL ctl_stop('ptr_vjk: ERROR - failed to release workspace array')
+      CALL wrk_dealloc( jpj*jpk, zwork )
 #endif
       !
@@ -281,8 +282,4 @@
       !! ** Action  : - p_fval: i-sum of e1t*e3t*pta
       !!----------------------------------------------------------------------
-#if defined key_mpp_mpi
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zwork => wrk_1d_1
-#endif
       !!
       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! tracer flux array at T-point
@@ -296,10 +293,12 @@
       INTEGER               ::   ijpjjpk
 #endif
+#if defined key_mpp_mpi
+      REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point
+#endif
       !!-------------------------------------------------------------------- 
       !
 #if defined key_mpp_mpi
-      IF( wrk_in_use(1, 1) ) THEN
-         CALL ctl_stop('ptr_tjk: requested workspace array unavailable')   ;   RETURN
-      ENDIF
+      ijpjjpk = jpj*jpk
+      CALL wrk_alloc( jpj*jpk, zwork )
 #endif
 
@@ -315,5 +314,4 @@
       END DO
 #if defined key_mpp_mpi
-      ijpjjpk = jpj*jpk
       ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk
       zwork(1:ijpjjpk)= RESHAPE( p_fval, ish )
@@ -323,5 +321,5 @@
       !
 #if defined key_mpp_mpi
-      IF( wrk_not_released(1, 1) )   CALL ctl_stop('ptr_tjk: failed to release workspace array')
+      CALL wrk_dealloc( jpj*jpk, zwork )
 #endif
       !    
@@ -343,4 +341,6 @@
       !!----------------------------------------------------------------------
       !
+      IF( nn_timing == 1 )   CALL timing_start('dia_ptr')
+      !
       IF( kt == nit000 .OR. MOD( kt, nn_fptr ) == 0 )   THEN
          !
@@ -349,5 +349,5 @@
             IF( ln_diaznl ) THEN               ! i-mean temperature and salinity
                DO jn = 1, nptr
-                  tn_jk(:,:,jn) = ptr_tjk( tn(:,:,:), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
+                  tn_jk(:,:,jn) = ptr_tjk( tsn(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
                END DO
             ENDIF
@@ -368,5 +368,5 @@
             !
             !                          ! Transports
-            !                                ! local heat & salt transports at T-points  ( tn*mj[vn+v_eiv] )
+            !                                ! local heat & salt transports at T-points  ( tsn*mj[vn+v_eiv] )
             vt(:,:,jpk) = 0._wp   ;   vs(:,:,jpk) = 0._wp
             DO jk= 1, jpkm1
@@ -378,6 +378,6 @@
                      zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) ) * 0.5_wp
 #endif 
-                     vt(:,jj,jk) = zv * tn(:,jj,jk)
-                     vs(:,jj,jk) = zv * sn(:,jj,jk)
+                     vt(:,jj,jk) = zv * tsn(:,jj,jk,jp_tem)
+                     vs(:,jj,jk) = zv * tsn(:,jj,jk,jp_sal)
                   END DO
                END DO
@@ -430,5 +430,11 @@
       ENDIF
       !
-      IF( kt == nitend )   CALL histclo( numptr )      ! Close the file
+#if defined key_mpp_mpi
+      IF( kt == nitend .AND. l_znl_root )   CALL histclo( numptr )      ! Close the file
+#else
+      IF( kt == nitend )                    CALL histclo( numptr )      ! Close the file
+#endif
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr')
       !
    END SUBROUTINE dia_ptr
@@ -449,7 +455,5 @@
       NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, ln_ptrcomp, nn_fptr, nn_fwri
       !!----------------------------------------------------------------------
-
-      !                                      ! allocate dia_ptr arrays
-      IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate arrays' )
+      IF( nn_timing == 1 )   CALL timing_start('dia_ptr_init')
 
       REWIND( numnam )                 ! Read Namelist namptr : poleward transport parameters
@@ -468,55 +472,63 @@
          WRITE(numout,*) '      Frequency of outputs                               nn_fwri    = ', nn_fwri
       ENDIF
-
-      IF( ln_subbas ) THEN   ;   nptr = 5       ! Global, Atlantic, Pacific, Indian, Indo-Pacific
-      ELSE                   ;   nptr = 1       ! Global only
-      ENDIF
-
-      rc_pwatt = rc_pwatt * rau0 * rcp          ! conversion from K.s-1 to PetaWatt
-
-      IF( .NOT. ln_diaptr ) THEN       ! diaptr not used
-        RETURN
-      ENDIF
       
-      IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
-
-      IF( ln_subbas ) THEN                ! load sub-basin mask
-         CALL iom_open( 'subbasins', inum )
-         CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
-         CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
-         CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
-         CALL iom_close( inum )
-         btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
-         WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
-         ELSE WHERE                     ;   btm30(:,:) = tmask(:,:,1)
-         END WHERE
-      ENDIF
-      btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
+      IF( ln_diaptr) THEN  
       
-      DO jn = 1, nptr
-         btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
-      END DO
+         IF( ln_subbas ) THEN   ;   nptr = 5       ! Global, Atlantic, Pacific, Indian, Indo-Pacific
+         ELSE                   ;   nptr = 1       ! Global only
+         ENDIF
+
+         !                                      ! allocate dia_ptr arrays
+         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
+
+         rc_pwatt = rc_pwatt * rau0 * rcp          ! conversion from K.s-1 to PetaWatt
+
+         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
+
+         IF( ln_subbas ) THEN                ! load sub-basin mask
+            CALL iom_open( 'subbasins', inum )
+            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
+            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
+            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
+            CALL iom_close( inum )
+            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
+            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
+            ELSE WHERE                     ;   btm30(:,:) = tmask(:,:,1)
+            END WHERE
+         ENDIF
+         btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
       
-      IF( lk_vvl )   CALL ctl_stop( 'diaptr: error in vvl case as constant i-mean surface is used' )
-
-      !                                   ! i-sum of e1v*e3v surface and its inverse
-      DO jn = 1, nptr
-         sjk(:,:,jn) = ptr_tjk( tmask(:,:,:), btmsk(:,:,jn) )
-         r1_sjk(:,:,jn) = 0._wp
-         WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
-      END DO
+         DO jn = 1, nptr
+            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
+         END DO
+      
+         IF( lk_vvl )   CALL ctl_stop( 'diaptr: error in vvl case as constant i-mean surface is used' )
+
+         !                                   ! i-sum of e1v*e3v surface and its inverse
+         DO jn = 1, nptr
+            sjk(:,:,jn) = ptr_tjk( tmask(:,:,:), btmsk(:,:,jn) )
+            r1_sjk(:,:,jn) = 0._wp
+            WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
+         END DO
+
+      ! Initialise arrays to zero because diatpr is called before they are first calculated
+      ! Note that this means diagnostics will not be exactly correct when model run is restarted.
+      htr_adv(:) = 0._wp ; str_adv(:) =  0._wp ;  htr_ldf(:) = 0._wp ; str_ldf(:) =  0._wp
 
 #if defined key_mpp_mpi 
-      iglo (1) = jpjglo                   ! MPP case using MPI  ('key_mpp_mpi')
-      iloc (1) = nlcj
-      iabsf(1) = njmppt(narea)
-      iabsl(:) = iabsf(:) + iloc(:) - 1
-      ihals(1) = nldj - 1
-      ihale(1) = nlcj - nlej
-      idid (1) = 2
-      CALL flio_dom_set( jpnj, nproc/jpni, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_ptr )
+         iglo (1) = jpjglo                   ! MPP case using MPI  ('key_mpp_mpi')
+         iloc (1) = nlcj
+         iabsf(1) = njmppt(narea)
+         iabsl(:) = iabsf(:) + iloc(:) - 1
+         ihals(1) = nldj - 1
+         ihale(1) = nlcj - nlej
+         idid (1) = 2
+         CALL flio_dom_set( jpnj, nproc/jpni, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_ptr )
 #else
-      nidom_ptr = FLIO_DOM_NONE
-#endif
+         nidom_ptr = FLIO_DOM_NONE
+#endif
+      ENDIF 
+      ! 
+      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr_init')
       ! 
    END SUBROUTINE dia_ptr_init
@@ -531,7 +543,4 @@
       !! ** Method  :   NetCDF file
       !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zphi => wrk_1d_1, zfoo => wrk_1d_2    ! 1D workspace
-      USE wrk_nemo, ONLY:   z_1  => wrk_2d_1                      ! 2D      -
       !!
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
@@ -548,9 +557,11 @@
 #endif
       REAL(wp)           ::   zsto, zout, zdt, zjulian                ! temporary scalars
-      !!----------------------------------------------------------------------
-
-      IF( wrk_in_use(1, 1,2) .OR. wrk_in_use(2, 1) ) THEN
-         CALL ctl_stop('dia_ptr_wri: requested workspace arrays unavailable')   ;   RETURN
-      ENDIF
+      !!
+      REAL(wp), POINTER, DIMENSION(:)   ::   zphi, zfoo    ! 1D workspace
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z_1           ! 2D workspace
+      !!-------------------------------------------------------------------- 
+      !
+      CALL wrk_alloc( jpi      , zphi , zfoo )
+      CALL wrk_alloc( jpi , jpk, z_1 )
 
       ! define time axis
@@ -866,6 +877,6 @@
       ENDIF
       !
-      IF( wrk_not_released(1, 1,2) .OR.    &
-          wrk_not_released(2, 1)    )   CALL ctl_stop('dia_ptr_wri: failed to release workspace arrays')
+      CALL wrk_dealloc( jpi      , zphi , zfoo )
+      CALL wrk_dealloc( jpi , jpk, z_1 )
       !
   END SUBROUTINE dia_ptr_wri
