Index: NEMO/trunk/src/OCE/DIA/diaptr.F90
===================================================================
--- NEMO/trunk/src/OCE/DIA/diaptr.F90	(revision 12276)
+++ NEMO/trunk/src/OCE/DIA/diaptr.F90	(revision 12377)
@@ -46,6 +46,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional)
 
-   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F)
-   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation
+   LOGICAL , PUBLIC ::   l_diaptr        !: tracers  trend flag (set from namelist in trdini)
    INTEGER, PARAMETER, PUBLIC ::   nptr = 5  ! (glo, atl, pac, ind, ipc)
 
@@ -60,6 +59,7 @@
    REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:) :: p_fval2d
 
+   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini)
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -69,8 +69,10 @@
 CONTAINS
 
-   SUBROUTINE dia_ptr( pvtr )
+   SUBROUTINE dia_ptr( kt, Kmm, pvtr )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dia_ptr  ***
       !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index     
+      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
       !
@@ -92,5 +94,9 @@
       !
       IF( ln_timing )   CALL timing_start('dia_ptr')
-      !
+
+      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init
+      !
+      IF( .NOT. l_diaptr )   RETURN
+
       IF( PRESENT( pvtr ) ) THEN
          IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF
@@ -111,14 +117,10 @@
             zmask(:,:,:) = 0._wp
             zts(:,:,:,:) = 0._wp
-            DO jk = 1, jpkm1
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpi
-                     zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk)
-                     zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
-                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
-                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
-                  ENDDO
-               ENDDO
-             ENDDO
+            DO_3D_10_11( 1, jpkm1 )
+               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
+               zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
+               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid
+               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
+            END_3D
          ENDIF
          IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
@@ -186,14 +188,10 @@
          zts(:,:,:,:) = 0._wp
          IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface 
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zsfc = e1t(ji,jj) * e3t_n(ji,jj,jk)
-                     zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
-                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc
-                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc
-                  END DO
-               END DO
-            END DO
+            DO_3D_11_11( 1, jpkm1 )
+               zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm)
+               zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
+               zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc
+               zts(ji,jj,jk,jp_sal) = ts(ji,jj,jk,jp_sal,Kmm) * zsfc
+            END_3D
             !
             DO jn = 1, nptr
@@ -280,13 +278,9 @@
          IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
             zts(:,:,:,:) = 0._wp
-            DO jk = 1, jpkm1
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpi
-                     zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk)
-                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
-                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
-                  ENDDO
-               ENDDO
-             ENDDO
+            DO_3D_10_11( 1, jpkm1 )
+               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
+               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid
+               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
+            END_3D
              CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) )
              CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) )
@@ -326,20 +320,18 @@
       !! ** Purpose :   Initialization, namelist read
       !!----------------------------------------------------------------------
-      INTEGER ::  inum, jn, ios, ierr           ! local integers
-      !!
-      NAMELIST/namptr/ ln_diaptr, ln_subbas
+      INTEGER ::  inum, jn           ! local integers
+      !!
       REAL(wp), DIMENSION(jpi,jpj) :: zmsk
       !!----------------------------------------------------------------------
 
-
-      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
-      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist' )
-
-      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport
-      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist' )
-      IF(lwm) WRITE ( numond, namptr )
-
+      l_diaptr = .FALSE.
+      IF(   iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  &
+         &  iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  &
+         &  iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  &
+         &  iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  & 
+         &  iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  &
+         &  iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) )  l_diaptr  = .TRUE.
+
+ 
       IF(lwp) THEN                     ! Control print
          WRITE(numout,*)
@@ -347,8 +339,8 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namptr : set ptr parameters'
-         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr
+         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr
       ENDIF
 
-      IF( ln_diaptr ) THEN  
+      IF( l_diaptr ) THEN  
          !
          IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
@@ -389,4 +381,6 @@
          hstr_vtr(:,:,:) = 0._wp           !
          !
+         ll_init = .FALSE.
+         !
       ENDIF 
       ! 
@@ -394,5 +388,5 @@
 
 
-   SUBROUTINE dia_ptr_hst( ktra, cptr, pva ) 
+   SUBROUTINE dia_ptr_hst( ktra, cptr, pvflx ) 
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_ptr_hst ***
@@ -403,5 +397,5 @@
       INTEGER                         , INTENT(in )  :: ktra  ! tracer index
       CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pvflx   ! 3D input array of advection/diffusion
       INTEGER                                        :: jn    !
 
@@ -410,10 +404,10 @@
          IF( ktra == jp_tem )  THEN
              DO jn = 1, nptr
-                hstr_adv(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_adv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
          IF( ktra == jp_sal )  THEN
              DO jn = 1, nptr
-                hstr_adv(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_adv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
@@ -423,10 +417,10 @@
          IF( ktra == jp_tem )  THEN
              DO jn = 1, nptr
-                hstr_ldf(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_ldf(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
          IF( ktra == jp_sal )  THEN
              DO jn = 1, nptr
-                hstr_ldf(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_ldf(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
@@ -436,10 +430,10 @@
          IF( ktra == jp_tem )  THEN
              DO jn = 1, nptr
-                hstr_eiv(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_eiv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
          IF( ktra == jp_sal )  THEN
              DO jn = 1, nptr
-                hstr_eiv(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_eiv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
@@ -449,10 +443,10 @@
          IF( ktra == jp_tem )  THEN
              DO jn = 1, nptr
-                hstr_vtr(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_vtr(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
          IF( ktra == jp_sal )  THEN
              DO jn = 1, nptr
-                hstr_vtr(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_vtr(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
@@ -486,5 +480,5 @@
 
 
-   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval )
+   FUNCTION ptr_sj_3d( pvflx, pmsk )   RESULT ( p_fval )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ptr_sj_3d  ***
@@ -492,10 +486,10 @@
       !! ** Purpose :   i-k sum computation of a j-flux array
       !!
-      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
-      !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
-      !!
-      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)  ::   pva   ! mask flux array at V-point
+      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
+      !!              pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
+      !!
+      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)  ::   pvflx  ! mask flux array at V-point
       REAL(wp), INTENT(in), DIMENSION(jpi,jpj)      ::   pmsk   ! Optional 2D basin mask
       !
@@ -509,11 +503,7 @@
       ijpj = jpj
       p_fval(:) = 0._wp
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! Vector opt.
-               p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
-            END DO
-         END DO
-      END DO
+      DO_3D_00_00( 1, jpkm1 )
+         p_fval(jj) = p_fval(jj) + pvflx(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
+      END_3D
 #if defined key_mpp_mpi
       CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl)
@@ -523,16 +513,16 @@
 
 
-   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval )
+   FUNCTION ptr_sj_2d( pvflx, pmsk )   RESULT ( p_fval )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ptr_sj_2d  ***
       !!
-      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array
-      !!
-      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
-      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
-      !!
-      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!----------------------------------------------------------------------
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point
+      !! ** Purpose :   "zonal" and vertical sum computation of a j-flux array
+      !!
+      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
+      !!      pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
+      !!
+      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
+      !!----------------------------------------------------------------------
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pvflx  ! mask flux array at V-point
       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
       !
@@ -546,9 +536,7 @@
       ijpj = jpj
       p_fval(:) = 0._wp
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! Vector opt.
-            p_fval(jj) = p_fval(jj) + pva(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj)
-         END DO
-      END DO
+      DO_2D_00_00
+         p_fval(jj) = p_fval(jj) + pvflx(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj)
+      END_2D
 #if defined key_mpp_mpi
       CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl )
@@ -577,9 +565,7 @@
       p_fval(:,:) = 0._wp
       DO jc = 1, jpnj ! looping over all processors in j axis
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! Vector opt.
-               p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj)
-            END DO
-         END DO
+         DO_2D_00_00
+            p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj)
+         END_2D
          CALL lbc_lnk( 'diaptr', p_fval, 'U', -1. )
       END DO
@@ -595,7 +581,7 @@
       !! ** Purpose :   i-sum computation of an array
       !!
-      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
-      !!
-      !! ** Action  : - p_fval: i-mean poleward flux of pva
+      !! ** Method  : - i-sum of field using the interior 2D vmask (pmsk).
+      !!
+      !! ** Action  : - p_fval: i-sum of masked field
       !!----------------------------------------------------------------------
       !!
@@ -618,11 +604,7 @@
       p_fval(:,:) = 0._wp
       !
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! Vector opt.
-               p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
-            END DO
-         END DO
-      END DO
+      DO_3D_00_00( 1, jpkm1 )
+         p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
+      END_3D
       !
 #if defined key_mpp_mpi
