Index: /trunk/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 1344)
+++ /trunk/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 1345)
@@ -6,4 +6,5 @@
    !! History :  9.0  !  03-09  (C. Talandier, G. Madec)  Original code
    !!            9.0  !  06-01  (A. Biastoch)  Allow sub-basins computation
+   !!            9.0  !  03-09  (O. Marti) Add fields
    !!----------------------------------------------------------------------
 
@@ -43,18 +44,24 @@
    LOGICAL , PUBLIC                 ::   ln_diaptr = .FALSE.   !: Poleward transport flag (T) or not (F)
    LOGICAL , PUBLIC                 ::   ln_subbas = .FALSE.   !: Atlantic/Pacific/Indian basins calculation
-   LOGICAL , PUBLIC                 ::   ln_diaznl             !: Add zonal means and meridional stream functions
+   LOGICAL , PUBLIC                 ::   ln_diaznl = .FALSE.   !: Add zonal means and meridional stream functions
+   LOGICAL , PUBLIC                 ::   ln_ptrcomp = .FALSE.  !: Add decomposition : overturning (and gyre, soon ...)
    INTEGER , PUBLIC                 ::   nf_ptr = 15           !: frequency of ptr computation
    INTEGER , PUBLIC                 ::   nf_ptr_wri = 15       !: frequency of ptr outputs
 
-   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_adv, pst_adv      !: heat and salt poleward transport: advection
-   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_ove, pst_ove      !: heat and salt poleward transport: overturning
+   REAL(wp), DIMENSION(jpi,jpj), PUBLIC ::   abasin, pbasin, ibasin, dbasin, sbasin !: Sub basin masks
+
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_adv, pst_adv  !: heat and salt poleward transport: advection
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_ove_glo, pst_ove_glo, pht_ove_atl, pst_ove_atl, pht_ove_pac, pst_ove_pac, &
+      &        pht_ove_ind, pst_ove_ind, pht_ove_ipc, pst_ove_ipc  !: heat and salt poleward transport: overturning
    REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_ldf, pst_ldf      !: heat and salt poleward transport: lateral diffusion
 #if defined key_diaeiv
-   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_eiv, pst_eiv      !: heat and salt poleward transport: bolus advection
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_eiv_glo, pst_eiv_glo, pht_eiv_atl, pst_eiv_atl, pht_eiv_pac, pst_eiv_pac, &
+      &        pht_eiv_ind, pst_eiv_ind, pht_eiv_ipc, pst_eiv_ipc !: heat and salt poleward transport: bolus advection
 #endif
    REAL(wp), PUBLIC, DIMENSION(jpj) ::   ht_glo,ht_atl,ht_ind,ht_pac,ht_ipc !: heat
    REAL(wp), PUBLIC, DIMENSION(jpj) ::   st_glo,st_atl,st_ind,st_pac,st_ipc !: salt
 
-   INTEGER :: nidom_diaptr ! domain identifier for IOIPSL
+   INTEGER :: nidom_diaptr = FLIO_DOM_NONE ! domain identifier for IOIPSL
+   INTEGER :: niter
 
    REAL(wp), DIMENSION(jpj,jpk) ::   tn_jk_glo, sn_jk_glo,  &  !: "zonal" mean temperature and salinity
@@ -79,8 +86,7 @@
       &                              surf_jk_r_ipc       
 #if defined key_diaeiv
-   REAL(wp), DIMENSION(jpj,jpk) ::   v_msf_eiv              !: bolus "meridional" Stream-Function
-#endif
-   REAL(wp), DIMENSION(jpi,jpj) ::   abasin, pbasin, ibasin, dbasin, sbasin !: Sub basin masks
-
+   REAL(wp), DIMENSION(jpj,jpk) ::   v_msf_eiv_glo, v_msf_eiv_atl, v_msf_eiv_pac, v_msf_eiv_ind, v_msf_eiv_ipc !: bolus "meridional" Stream-Function
+#endif
+ 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -123,5 +129,5 @@
       END DO
       !
-      IF( lk_mpp )   CALL mpp_sum( p_fval, ijpj )     !!bug  I presume
+      IF( lk_mpp )   CALL mpp_sum( p_fval, ijpj, ncomm_znl)     !!bug  I presume
       !
    END FUNCTION ptr_vj_3d
@@ -150,10 +156,10 @@
       p_fval(:) = 0.e0
       DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! Vector opt.
+         DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
             p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj+1) * tmask_i(ji,jj)
          END DO
       END DO
       !
-      IF( lk_mpp )   CALL mpp_sum( p_fval, ijpj )     !!bug  I presume
+      IF( lk_mpp )   CALL mpp_sum( p_fval, ijpj, ncomm_znl )     !!bug  I presume
       ! 
    END FUNCTION ptr_vj_2d
@@ -186,5 +192,5 @@
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
+               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
                   p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk)   &
                      &                                          * tmask_i(ji,jj+1) * tmask_i(ji,jj) &
@@ -196,5 +202,5 @@
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
+               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
                   p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk)   &
                      &                                          * tmask_i(ji,jj+1) * tmask_i(ji,jj)
@@ -207,5 +213,5 @@
          ish(1) = jpj*jpk  ;  ish2(1) = jpj  ;  ish2(2) = jpk
          zwork(:)= RESHAPE( p_fval, ish )
-         CALL mpp_sum( zwork, jpj*jpk )
+         CALL mpp_sum( zwork, jpj*jpk, ncomm_znl )
          p_fval(:,:)= RESHAPE( zwork, ish2 )
       END IF
@@ -239,5 +245,5 @@
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! Vector opt.
+               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
                   p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk)                  &
                      &                          * e1t(ji,jj) * fse3t(ji,jj,jk)   &
@@ -250,5 +256,5 @@
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! Vector opt.
+               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
                   p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk)                  &
                      &                          * e1t(ji,jj) * fse3t(ji,jj,jk)   &
@@ -262,5 +268,5 @@
          ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk
          zwork(:)= RESHAPE( p_fval, ish )
-         CALL mpp_sum( zwork, jpj*jpk )
+         CALL mpp_sum( zwork, jpj*jpk, ncomm_znl )
          p_fval(:,:)= RESHAPE(zwork,ish2)
       END IF
@@ -284,140 +290,183 @@
       IF( kt == nit000 .OR. MOD( kt, nf_ptr ) == 0 )   THEN
 
-         zsverdrup = 1.e-6
-         zpwatt    = 1.e-15
-         zggram    = 1.e-6
-         
-         IF ( ln_diaznl ) THEN
-            ! "zonal" mean temperature and salinity at V-points
-            tn_jk_glo(:,:) = ptr_tjk( tn(:,:,:) ) * surf_jk_r_glo(:,:)
-            sn_jk_glo(:,:) = ptr_tjk( sn(:,:,:) ) * surf_jk_r_glo(:,:)
-            
-            IF (ln_subbas) THEN 
-               tn_jk_atl(:,:) = ptr_tjk( tn(:,:,:), abasin(:,:) ) * surf_jk_r_atl(:,:)
-               sn_jk_atl(:,:) = ptr_tjk( sn(:,:,:), abasin(:,:) ) * surf_jk_r_atl(:,:)
-               tn_jk_pac(:,:) = ptr_tjk( tn(:,:,:), pbasin(:,:) ) * surf_jk_r_pac(:,:)
-               sn_jk_pac(:,:) = ptr_tjk( sn(:,:,:), pbasin(:,:) ) * surf_jk_r_pac(:,:)
-               tn_jk_ind(:,:) = ptr_tjk( tn(:,:,:), ibasin(:,:) ) * surf_jk_r_ind(:,:)
-               sn_jk_ind(:,:) = ptr_tjk( sn(:,:,:), ibasin(:,:) ) * surf_jk_r_ind(:,:)
-               tn_jk_ipc(:,:) = ptr_tjk( tn(:,:,:), dbasin(:,:) ) * surf_jk_r_ipc(:,:)
-               sn_jk_ipc(:,:) = ptr_tjk( sn(:,:,:), dbasin(:,:) ) * surf_jk_r_ipc(:,:)
-            ENDIF
-         ENDIF
-         
-         !--------------------------------------------------------
-         ! overturning calculation:
-         
-         ! horizontal integral and vertical dz 
-
-#if defined key_diaeiv
-         v_msf_glo(:,:) = ptr_vjk( vn(:,:,:)+v_eiv(:,:,:) ) 
-         IF( ln_subbas .AND. ln_diaznl ) THEN
-            v_msf_atl(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), abasin(:,:)*sbasin(:,:) ) 
-            v_msf_pac(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), pbasin(:,:)*sbasin(:,:) ) 
-            v_msf_ind(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), ibasin(:,:)*sbasin(:,:) ) 
-            v_msf_ipc(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), dbasin(:,:)*sbasin(:,:) ) 
-         ENDIF
+         IF ( MOD( kt, nf_ptr ) == 0 ) THEN 
+
+            zsverdrup = 1.e-6
+            zpwatt    = 1.e-15
+            zggram    = 1.e-6
+
+            IF ( ln_diaznl ) THEN
+               ! "zonal" mean temperature and salinity at V-points
+               tn_jk_glo(:,:) = ptr_tjk( tn(:,:,:) ) * surf_jk_r_glo(:,:)
+               sn_jk_glo(:,:) = ptr_tjk( sn(:,:,:) ) * surf_jk_r_glo(:,:)
+
+               IF (ln_subbas) THEN 
+                  tn_jk_atl(:,:) = ptr_tjk( tn(:,:,:), abasin(:,:) ) * surf_jk_r_atl(:,:)
+                  sn_jk_atl(:,:) = ptr_tjk( sn(:,:,:), abasin(:,:) ) * surf_jk_r_atl(:,:)
+                  tn_jk_pac(:,:) = ptr_tjk( tn(:,:,:), pbasin(:,:) ) * surf_jk_r_pac(:,:)
+                  sn_jk_pac(:,:) = ptr_tjk( sn(:,:,:), pbasin(:,:) ) * surf_jk_r_pac(:,:)
+                  tn_jk_ind(:,:) = ptr_tjk( tn(:,:,:), ibasin(:,:) ) * surf_jk_r_ind(:,:)
+                  sn_jk_ind(:,:) = ptr_tjk( sn(:,:,:), ibasin(:,:) ) * surf_jk_r_ind(:,:)
+                  tn_jk_ipc(:,:) = ptr_tjk( tn(:,:,:), dbasin(:,:) ) * surf_jk_r_ipc(:,:)
+                  sn_jk_ipc(:,:) = ptr_tjk( sn(:,:,:), dbasin(:,:) ) * surf_jk_r_ipc(:,:)
+               ENDIF
+            ENDIF
+
+            !--------------------------------------------------------
+            ! overturning calculation:
+
+            ! horizontal integral and vertical dz 
+
+#if defined key_diaeiv
+            v_msf_glo(:,:) = ptr_vjk( vn(:,:,:)+v_eiv(:,:,:) ) 
+            IF( ln_subbas .AND. ln_diaznl ) THEN
+               v_msf_atl(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), abasin(:,:)*sbasin(:,:) ) 
+               v_msf_pac(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), pbasin(:,:)*sbasin(:,:) ) 
+               v_msf_ind(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), ibasin(:,:)*sbasin(:,:) ) 
+               v_msf_ipc(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), dbasin(:,:)*sbasin(:,:) ) 
+            ENDIF
 #else
-         v_msf_glo(:,:) = ptr_vjk( vn(:,:,:) ) 
-         IF( ln_subbas .AND. ln_diaznl ) THEN
-            v_msf_atl(:,:) = ptr_vjk( vn (:,:,:), abasin(:,:)*sbasin(:,:) ) 
-            v_msf_pac(:,:) = ptr_vjk( vn (:,:,:), pbasin(:,:)*sbasin(:,:) ) 
-            v_msf_ind(:,:) = ptr_vjk( vn (:,:,:), ibasin(:,:)*sbasin(:,:) ) 
-            v_msf_ipc(:,:) = ptr_vjk( vn (:,:,:), dbasin(:,:)*sbasin(:,:) ) 
-         ENDIF
-#endif
-         
-#if defined key_diaeiv
-         v_msf_eiv(:,:) = ptr_vjk( v_eiv(:,:,:) ) 
-#endif
-         
-         ! "Meridional" Stream-Function
-         DO jk = 2,jpk 
-            v_msf_glo(:,jk) = v_msf_glo(:,jk-1) + v_msf_glo(:,jk)
-         END DO
-         v_msf_glo(:,:) = v_msf_glo(:,:) * zsverdrup
-#if defined key_diaeiv
-         ! Bolus "Meridional" Stream-Function
-         DO jk = 2,jpk
-            v_msf_eiv(:,jk) = v_msf_eiv(:,jk-1) + v_msf_eiv(:,jk)
-         END DO
-         v_msf_eiv(:,:) = v_msf_eiv(:,:) * zsverdrup
-#endif
-         !
-         IF( ln_subbas .AND. ln_diaznl ) THEN
+            v_msf_glo(:,:) = ptr_vjk( vn(:,:,:) ) 
+            IF( ln_subbas .AND. ln_diaznl ) THEN
+               v_msf_atl(:,:) = ptr_vjk( vn (:,:,:), abasin(:,:)*sbasin(:,:) ) 
+               v_msf_pac(:,:) = ptr_vjk( vn (:,:,:), pbasin(:,:)*sbasin(:,:) ) 
+               v_msf_ind(:,:) = ptr_vjk( vn (:,:,:), ibasin(:,:)*sbasin(:,:) ) 
+               v_msf_ipc(:,:) = ptr_vjk( vn (:,:,:), dbasin(:,:)*sbasin(:,:) ) 
+            ENDIF
+#endif
+
+#if defined key_diaeiv
+            v_msf_eiv_glo(:,:) = ptr_vjk( v_eiv(:,:,:) )
+            IF (ln_subbas ) THEN 
+               v_msf_eiv_atl(:,:) = ptr_vjk( v_eiv(:,:,:), abasin(:,:)*sbasin(:,:) )
+               v_msf_eiv_pac(:,:) = ptr_vjk( v_eiv(:,:,:), pbasin(:,:)*sbasin(:,:) )
+               v_msf_eiv_ind(:,:) = ptr_vjk( v_eiv(:,:,:), ibasin(:,:)*sbasin(:,:) )
+               v_msf_eiv_ipc(:,:) = ptr_vjk( v_eiv(:,:,:), dbasin(:,:)*sbasin(:,:) )
+            END IF
+#endif
+
+            ! "Meridional" Stream-Function
             DO jk = 2,jpk 
-               v_msf_atl(:,jk) = v_msf_atl(:,jk-1) + v_msf_atl(:,jk)
-               v_msf_pac(:,jk) = v_msf_pac(:,jk-1) + v_msf_pac(:,jk)
-               v_msf_ind(:,jk) = v_msf_ind(:,jk-1) + v_msf_ind(:,jk)
-               v_msf_ipc(:,jk) = v_msf_ipc(:,jk-1) + v_msf_ipc(:,jk)
+               v_msf_glo(:,jk) = v_msf_glo(:,jk-1) + v_msf_glo(:,jk)
             END DO
-            v_msf_atl(:,:) = v_msf_atl(:,:) * zsverdrup
-            v_msf_pac(:,:) = v_msf_pac(:,:) * zsverdrup
-            v_msf_ind(:,:) = v_msf_ind(:,:) * zsverdrup
-            v_msf_ipc(:,:) = v_msf_ipc(:,:) * zsverdrup
-         ENDIF
-            
-         ! Transports
-         ! T times V on T points (include bolus velocities)
+            v_msf_glo(:,:) = v_msf_glo(:,:) * zsverdrup
+#if defined key_diaeiv
+            ! Bolus "Meridional" Stream-Function
+            DO jk = 2,jpk
+               v_msf_eiv_glo(:,jk) = v_msf_eiv_glo(:,jk-1) + v_msf_eiv_glo(:,jk)
+            END DO
+            v_msf_eiv_glo(:,:) = v_msf_eiv_glo(:,:) * zsverdrup
+            IF ( ln_subbas ) THEN 
+               DO jk = 2,jpk
+                  v_msf_eiv_atl(:,jk) = v_msf_eiv_atl(:,jk-1) + v_msf_eiv_atl(:,jk)
+                  v_msf_eiv_pac(:,jk) = v_msf_eiv_pac(:,jk-1) + v_msf_eiv_pac(:,jk)
+                  v_msf_eiv_ind(:,jk) = v_msf_eiv_ind(:,jk-1) + v_msf_eiv_ind(:,jk)
+                  v_msf_eiv_ipc(:,jk) = v_msf_eiv_ipc(:,jk-1) + v_msf_eiv_ipc(:,jk)
+               END DO
+            ENDIF
+#endif
+            !
+            IF( ln_subbas .AND. ln_diaznl ) THEN
+               DO jk = 2,jpk 
+                  v_msf_atl(:,jk) = v_msf_atl(:,jk-1) + v_msf_atl(:,jk)
+                  v_msf_pac(:,jk) = v_msf_pac(:,jk-1) + v_msf_pac(:,jk)
+                  v_msf_ind(:,jk) = v_msf_ind(:,jk-1) + v_msf_ind(:,jk)
+                  v_msf_ipc(:,jk) = v_msf_ipc(:,jk-1) + v_msf_ipc(:,jk)
+               END DO
+               v_msf_atl(:,:) = v_msf_atl(:,:) * zsverdrup
+               v_msf_pac(:,:) = v_msf_pac(:,:) * zsverdrup
+               v_msf_ind(:,:) = v_msf_ind(:,:) * zsverdrup
+               v_msf_ipc(:,:) = v_msf_ipc(:,:) * zsverdrup
+            ENDIF
+
+            ! Transports
+            ! T times V on T points (include bolus velocities)
 #if defined key_diaeiv 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + u_eiv(ji,jj,:) + u_eiv(ji,jj-1,:) )*0.5
-               vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + v_eiv(ji,jj,:) + v_eiv(ji,jj-1,:) )*0.5
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + u_eiv(ji,jj,:) + u_eiv(ji,jj-1,:) )*0.5
+                  vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + v_eiv(ji,jj,:) + v_eiv(ji,jj-1,:) )*0.5
+               END DO
             END DO
-         END DO
 #else
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) )*0.5
-               vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:)  )*0.5
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) )*0.5
+                  vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:)  )*0.5
+               END DO
             END DO
-         END DO
 #endif 
-         
-         ht_glo(:) = SUM( ptr_vjk( vt(:,:,:)), 2 )
-         st_glo(:) = SUM( ptr_vjk( vs(:,:,:)), 2 )
-         
-         IF ( ln_subbas ) THEN 
-            ht_atl(:) = SUM( ptr_vjk( vt (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
-            ht_pac(:) = SUM( ptr_vjk( vt (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
-            ht_ind(:) = SUM( ptr_vjk( vt (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
-            ht_ipc(:) = SUM( ptr_vjk( vt (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
-            st_atl(:) = SUM( ptr_vjk( vs (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
-            st_pac(:) = SUM( ptr_vjk( vs (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
-            st_ind(:) = SUM( ptr_vjk( vs (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
-            st_ipc(:) = SUM( ptr_vjk( vs (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
-         ENDIF
-         
-         ! poleward tracer transports: 
-         ! overturning components:
-         pht_ove(:) = SUM( v_msf_glo(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
-         pst_ove(:) = SUM( v_msf_glo(:,:) * sn_jk_glo(:,:), 2 )   ! SUM over jk
-#if defined key_diaeiv
-         pht_eiv(:) = SUM( v_msf_eiv(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
-         pst_eiv(:) = SUM( v_msf_eiv(:,:) * sn_jk_glo(:,:), 2 )   ! SUM over jk
-#endif
-
-         ! conversion in PW and G g
-         zpwatt = zpwatt * rau0 * rcp
-         pht_adv(:) = pht_adv(:) * zpwatt  
-         pht_ove(:) = pht_ove(:) * zpwatt
-         pht_ldf(:) = pht_ldf(:) * zpwatt
-         pst_adv(:) = pst_adv(:) * zggram
-         pst_ove(:) = pst_ove(:) * zggram
-         pst_ldf(:) = pst_ldf(:) * zggram
-#if defined key_diaeiv
-         pht_eiv(:) = pht_eiv(:) * zpwatt
-         pst_eiv(:) = pst_eiv(:) * zggram
-#endif
-         IF( ln_subbas ) THEN
-            ht_atl(:) = ht_atl(:) * zpwatt
-            ht_pac(:) = ht_pac(:) * zpwatt
-            ht_ind(:) = ht_ind(:) * zpwatt
-            ht_ipc(:) = ht_ipc(:) * zpwatt
-            st_atl(:) = st_atl(:) * zggram 
-            st_pac(:) = st_pac(:) * zggram
-            st_ind(:) = st_ind(:) * zggram
-            st_ipc(:) = st_ipc(:) * zggram
+
+            ht_glo(:) = SUM( ptr_vjk( vt(:,:,:)), 2 )
+            st_glo(:) = SUM( ptr_vjk( vs(:,:,:)), 2 )
+
+            IF ( ln_subbas ) THEN 
+               ht_atl(:) = SUM( ptr_vjk( vt (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
+               ht_pac(:) = SUM( ptr_vjk( vt (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
+               ht_ind(:) = SUM( ptr_vjk( vt (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
+               ht_ipc(:) = SUM( ptr_vjk( vt (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
+               st_atl(:) = SUM( ptr_vjk( vs (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
+               st_pac(:) = SUM( ptr_vjk( vs (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
+               st_ind(:) = SUM( ptr_vjk( vs (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
+               st_ipc(:) = SUM( ptr_vjk( vs (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
+            ENDIF
+
+            ! poleward tracer transports: 
+            ! overturning components:
+            IF ( ln_ptrcomp ) THEN 
+               pht_ove_glo(:) = SUM( v_msf_glo(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
+               pst_ove_glo(:) = SUM( v_msf_glo(:,:) * sn_jk_glo(:,:), 2 )  
+               IF ( ln_subbas ) THEN 
+                  pht_ove_atl(:) = SUM( v_msf_atl(:,:) * tn_jk_atl(:,:), 2 )   ! SUM over jk
+                  pst_ove_atl(:) = SUM( v_msf_atl(:,:) * sn_jk_atl(:,:), 2 )  
+                  pht_ove_pac(:) = SUM( v_msf_pac(:,:) * tn_jk_pac(:,:), 2 )   ! SUM over jk
+                  pst_ove_pac(:) = SUM( v_msf_pac(:,:) * sn_jk_pac(:,:), 2 )  
+                  pht_ove_ind(:) = SUM( v_msf_ind(:,:) * tn_jk_ind(:,:), 2 )   ! SUM over jk
+                  pst_ove_ind(:) = SUM( v_msf_ind(:,:) * sn_jk_ind(:,:), 2 )  
+                  pht_ove_ipc(:) = SUM( v_msf_ipc(:,:) * tn_jk_ipc(:,:), 2 )   ! SUM over jk
+                  pst_ove_ipc(:) = SUM( v_msf_ipc(:,:) * sn_jk_ipc(:,:), 2 )  
+               END IF
+            END IF
+
+            ! Bolus component
+#if defined key_diaeiv
+            pht_eiv_glo(:) = SUM( v_msf_eiv_glo(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
+            pst_eiv_glo(:) = SUM( v_msf_eiv_glo(:,:) * sn_jk_glo(:,:), 2 )   ! SUM over jk
+            IF ( ln_subbas ) THEN 
+               pht_eiv_atl(:) = SUM( v_msf_eiv_glo(:,:) * tn_jk_atl(:,:), 2 )   ! SUM over jk
+               pst_eiv_atl(:) = SUM( v_msf_eiv_glo(:,:) * sn_jk_atl(:,:), 2 )   ! SUM over jk
+               pht_eiv_pac(:) = SUM( v_msf_eiv_pac(:,:) * tn_jk_pac(:,:), 2 )   ! SUM over jk
+               pst_eiv_pac(:) = SUM( v_msf_eiv_pac(:,:) * sn_jk_pac(:,:), 2 )   ! SUM over jk
+               pht_eiv_ind(:) = SUM( v_msf_eiv_ind(:,:) * tn_jk_ind(:,:), 2 )   ! SUM over jk
+               pst_eiv_ind(:) = SUM( v_msf_eiv_ind(:,:) * sn_jk_ind(:,:), 2 )   ! SUM over jk
+               pht_eiv_ipc(:) = SUM( v_msf_eiv_ipc(:,:) * tn_jk_ipc(:,:), 2 )   ! SUM over jk
+               pst_eiv_ipc(:) = SUM( v_msf_eiv_ipc(:,:) * sn_jk_ipc(:,:), 2 )   ! SUM over jk
+            ENDIF
+#endif
+
+            ! conversion in PW and G g
+            zpwatt = zpwatt * rau0 * rcp
+            pht_adv(:) = pht_adv(:) * zpwatt  
+            pht_ldf(:) = pht_ldf(:) * zpwatt
+            pst_adv(:) = pst_adv(:) * zggram
+            pst_ldf(:) = pst_ldf(:) * zggram
+            IF ( ln_ptrcomp ) THEN 
+               pht_ove_glo(:) = pht_ove_glo(:) * zpwatt
+               pst_ove_glo(:) = pst_ove_glo(:) * zggram
+            END IF
+#if defined key_diaeiv
+            pht_eiv_glo(:) = pht_eiv_glo(:) * zpwatt
+            pst_eiv_glo(:) = pst_eiv_glo(:) * zggram
+#endif
+            IF( ln_subbas ) THEN
+               ht_atl(:) = ht_atl(:) * zpwatt
+               ht_pac(:) = ht_pac(:) * zpwatt
+               ht_ind(:) = ht_ind(:) * zpwatt
+               ht_ipc(:) = ht_ipc(:) * zpwatt
+               st_atl(:) = st_atl(:) * zggram 
+               st_pac(:) = st_pac(:) * zggram
+               st_ind(:) = st_ind(:) * zggram
+               st_ipc(:) = st_ipc(:) * zggram
+            ENDIF
          ENDIF
 
@@ -439,6 +488,6 @@
       !! ** Purpose :   Initialization, namelist read
       !!----------------------------------------------------------------------
-      NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, nf_ptr, nf_ptr_wri
-      INTEGER ::  inum       ! temporary logical unit
+      NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, ln_ptrcomp, nf_ptr, nf_ptr_wri
+      INTEGER :: inum       ! temporary logical unit
       INTEGER, DIMENSION (1) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
       !!----------------------------------------------------------------------
@@ -458,4 +507,11 @@
          WRITE(numout,*) '             Frequency of computation                    nf_ptr = ', nf_ptr
          WRITE(numout,*) '             Frequency of outputs                    nf_ptr_wri = ', nf_ptr_wri
+      ENDIF
+
+      !
+      ! Define MPI communicator for zonal sum
+      !
+      IF( lk_mpp )  THEN
+         CALL mpp_ini_znl
       ENDIF
 
@@ -505,13 +561,13 @@
       idid (1) = 2
 
-      IF(lwp) THEN
-          WRITE(numout,*)
-          WRITE(numout,*) 'diaptr_init :   iloc  = ', iloc 
-          WRITE(numout,*) '~~~~~~~~~~~     iabsf = ', iabsf
-          WRITE(numout,*) '                ihals = ', ihals
-          WRITE(numout,*) '                ihale = ', ihale
-      ENDIF
-
-      CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_diaptr)
+!-$$      IF(lwp) THEN
+!-$$          WRITE(numout,*)
+!-$$          WRITE(numout,*) 'dia_ptr_init :   iloc  = ', iloc 
+!-$$          WRITE(numout,*) '~~~~~~~~~~~~     iabsf = ', iabsf
+!-$$          WRITE(numout,*) '                 ihals = ', ihals
+!-$$          WRITE(numout,*) '                 ihale = ', ihale
+!-$$      ENDIF 
+
+      CALL flio_dom_set ( jpnj, nproc/jpni, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_diaptr)
       
    END SUBROUTINE dia_ptr_init
@@ -544,10 +600,24 @@
 
       ! define time axis
-      it    = kt
+      it    = kt / nf_ptr
       itmod = kt - nit000 + 1
+      
+!-$$      IF(lwp) THEN
+!-$$         WRITE(numout,*)
+!-$$         WRITE(numout,*) 'dia_ptr_wri : kt = ', kt, 'it = ', it, ' itmod = ', itmod, ' niter = ', niter
+!-$$         WRITE(numout,*) '~~~~~~~~~~~~'
+!-$$      ENDIF
 
       ! Initialization
       ! --------------
       IF( kt == nit000 ) THEN
+
+         niter = (nit000 - 1) / nf_ptr
+
+!-$$         IF(lwp) THEN
+!-$$            WRITE(numout,*)
+!-$$            WRITE(numout,*) 'dia_ptr_wri : poleward transport and msf writing: initialization , niter = ', niter
+!-$$            WRITE(numout,*) '~~~~~~~~~~~~'
+!-$$         ENDIF
 
          zdt = rdt
@@ -562,20 +632,18 @@
             IF( jp_cfg == 05  )   iline = 192   ! i-line that passes near the North Pole
             IF( jp_cfg == 025 )   iline = 384   ! i-line that passes near the North Pole
+            IF( jp_cfg == 1   )   iline =  96   ! i-line that passes near the North Pole
             IF( jp_cfg == 2   )   iline =  48   ! i-line that passes near the North Pole
             IF( jp_cfg == 4   )   iline =  24   ! i-line that passes near the North Pole
-            zphi(:) = -9999.9999e0
+            zphi(:) = 0.e0
             DO ji = mi0(iline), mi1(iline) 
-               WRITE(numout,*) 'diaptr : ', nproc, narea, iline, ji, mi0(iline), mi1(iline), &
-                  &  mj0(jpjdta-1), mj1(jpjdta-1), mj0(jpjdta), mj1(jpjdta), '--'
-               CALL flush (numout)
                zphi(:) = gphiv(ji,:)         ! if iline is in the local domain
                ! Correct highest latitude for some configurations - will work if domain is parallelized in J ?
                IF( jp_cfg == 05 ) THEN
                   DO jj = mj0(jpjdta), mj1(jpjdta) 
-                     zphi( jj ) = zphi(jpjdta-1) + (zphi(jpjdta-1)-zphi(jpjdta-2))/2.
+                     zphi( jj ) = zphi(mj0(jpjdta-1)) + (zphi(mj0(jpjdta-1))-zphi(mj0(jpjdta-2)))/2.
                      zphi( jj ) = MIN( zphi(jj), 90.)
                   END DO
                END IF
-               IF( jp_cfg == 2 .OR. jp_cfg == 4 ) THEN
+               IF( jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) THEN
                   DO jj = mj0(jpjdta-1), mj1(jpjdta-1) 
                      zphi( jj ) = 88.5e0
@@ -587,268 +655,272 @@
             END DO
             ! provide the correct zphi to all local domains
-            DO jj = 1, jpj
-               WRITE(numout,*) 'diaptr(1) ', nproc, jj, mjg(jj), zphi(jj), '--'
-               CALL flush (numout)
-            ENDDO
-            IF( lk_mpp )   CALL mpp_sum( zphi, jpj )        
+            IF( lk_mpp )   CALL mpp_sum( zphi, jpj, ncomm_znl )        
 
             !                                        ! =======================
          ELSE                                        !   OTHER configurations  zjulian = zjulian - adatrj 
-                                                     !   set calendar origin to the beginning of the experiment
+            !   set calendar origin to the beginning of the experiment
             !                                        ! =======================
             zphi(:) = gphiv(1,:)             ! assume lat/lon coordinate, select the first i-line
             !
          ENDIF
-         DO jj = 1, jpj
-            WRITE(numout,*) 'diaptr(2) ', nproc, jj, mjg(jj), zphi(jj), '--'
-            CALL flush (numout)
-         ENDDO
-
-         ! OPEN netcdf file 
-         ! ----------------
-         ! Define frequency of output and means
-         zsto = nf_ptr * zdt
-         IF( ln_mskland )   THEN    ! put 1.e+20 on land (very expensive!!)
-            clop      = "ave(only(x))"
-            clop_once = "once(only(x))"
-         ELSE                       ! no use of the mask value (require less cpu time)
-            clop      = "ave(x)"       
-            clop_once = "once"
-         ENDIF
-
-         zout = nf_ptr_wri * zdt
-         zfoo(:) = 0.e0
-
-         ! Compute julian date from starting date of the run
-
-         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
-         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
-
-         CALL dia_nam( clhstnam, nf_ptr_wri, 'diaptr' )
-         IF(lwp)WRITE( numout,*)" Name of diaptr NETCDF file ",clhstnam
-
-         ! Horizontal grid : zphi()
-         CALL histbeg(clhstnam, 1, zfoo, jpj, zphi,   &
-            1, 1, 1, jpj, nit000-1, zjulian, zdt, nhoridz, numptr, domain_id=nidom_diaptr )
-         ! Vertical grids : gdept_0, gdepw_0
-         CALL histvert( numptr, "deptht", "Vertical T levels",   &
-            "m", jpk, gdept_0, ndepidzt, "down" )
-         CALL histvert( numptr, "depthw", "Vertical W levels",   &
-            "m", jpk, gdepw_0, ndepidzw, "down" )
-
          !
-         CALL wheneq ( jpj*jpk, MIN(surf_jk_glo(:,:), 1.e0), 1, 1., ndex  , ndim  )      ! Lat-Depth
-         CALL wheneq ( jpj    , MIN(surf_jk_glo(:,1), 1.e0), 1, 1., ndex_h, ndim_h )     ! Lat
-
-         IF (ln_subbas) THEN
-            z_1 (:,1) = 1.0e0
-            WHERE ( gphit (jpi/2,:) .LT. -30 ) z_1 (:,1) = 0.e0
-            DO jk = 2, jpk
-               z_1 (:,jk) = z_1 (:,1)
-            END DO
-
-            CALL wheneq ( jpj*jpk, MIN(surf_jk_atl(:,:)         , 1.e0), 1, 1., ndex_atl     , ndim_atl      ) ! Lat-Depth
-            CALL wheneq ( jpj*jpk, MIN(surf_jk_atl(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_atl_30  , ndim_atl_30   ) ! Lat-Depth
-            CALL wheneq ( jpj    , MIN(surf_jk_atl(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_atl_30, ndim_h_atl_30 ) ! Lat
-
-            CALL wheneq ( jpj*jpk, MIN(surf_jk_pac(:,:)         , 1.e0), 1, 1., ndex_pac     , ndim_pac      ) ! Lat-Depth
-            CALL wheneq ( jpj*jpk, MIN(surf_jk_pac(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_pac_30  , ndim_pac_30   ) ! Lat-Depth
-            CALL wheneq ( jpj    , MIN(surf_jk_pac(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_pac_30, ndim_h_pac_30 ) ! Lat
-
-            CALL wheneq ( jpj*jpk, MIN(surf_jk_ind(:,:)         , 1.e0), 1, 1., ndex_ind     , ndim_ind      ) ! Lat-Depth
-            CALL wheneq ( jpj*jpk, MIN(surf_jk_ind(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_ind_30  , ndim_ind_30   ) ! Lat-Depth
-            CALL wheneq ( jpj    , MIN(surf_jk_ind(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_ind_30, ndim_h_ind_30 ) ! Lat
-
-            CALL wheneq ( jpj*jpk, MIN(surf_jk_ipc(:,:)         , 1.e0), 1, 1., ndex_ipc     , ndim_ipc      ) ! Lat-Depth
-            CALL wheneq ( jpj*jpk, MIN(surf_jk_ipc(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_ipc_30  , ndim_ipc_30   ) ! Lat-Depth
-            CALL wheneq ( jpj    , MIN(surf_jk_ipc(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_ipc_30, ndim_h_ipc_30 ) ! Lat
-
-         ENDIF
-
-         ! 
-#if defined key_diaeiv
-         cl_comment = ' (Bolus part included)'
+         ! Work only on westmost processor (will not work if mppini2 is used)
+         IF ( l_znl_root ) THEN 
+            !
+            ! OPEN netcdf file 
+            ! ----------------
+            ! Define frequency of output and means
+            zsto = nf_ptr * zdt
+            IF( ln_mskland )   THEN    ! put 1.e+20 on land (very expensive!!)
+               clop      = "ave(only(x))"
+               clop_once = "once(only(x))"
+            ELSE                       ! no use of the mask value (require less cpu time)
+               clop      = "ave(x)"       
+               clop_once = "once"
+            ENDIF
+
+            zout = nf_ptr_wri * zdt
+            zfoo(:) = 0.e0
+
+            ! Compute julian date from starting date of the run
+
+            CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
+            zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
+
+            CALL dia_nam( clhstnam, nf_ptr_wri, 'diaptr' )
+            IF(lwp)WRITE( numout,*)" Name of diaptr NETCDF file : ", clhstnam
+
+            ! Horizontal grid : zphi()
+            CALL histbeg(clhstnam, 1, zfoo, jpj, zphi,   &
+               1, 1, 1, jpj, niter, zjulian, zdt*nf_ptr, nhoridz, numptr, domain_id=nidom_diaptr )
+            ! Vertical grids : gdept_0, gdepw_0
+            CALL histvert( numptr, "deptht", "Vertical T levels",   &
+               "m", jpk, gdept_0, ndepidzt, "down" )
+            CALL histvert( numptr, "depthw", "Vertical W levels",   &
+               "m", jpk, gdepw_0, ndepidzw, "down" )
+
+            !
+            CALL wheneq ( jpj*jpk, MIN(surf_jk_glo(:,:), 1.e0), 1, 1., ndex  , ndim  )      ! Lat-Depth
+            CALL wheneq ( jpj    , MIN(surf_jk_glo(:,1), 1.e0), 1, 1., ndex_h, ndim_h )     ! Lat
+
+            IF (ln_subbas) THEN
+               z_1 (:,1) = 1.0e0
+               WHERE ( gphit (jpi/2,:) .LT. -30 ) z_1 (:,1) = 0.e0
+               DO jk = 2, jpk
+                  z_1 (:,jk) = z_1 (:,1)
+               END DO
+
+               CALL wheneq ( jpj*jpk, MIN(surf_jk_atl(:,:)         , 1.e0), 1, 1., ndex_atl     , ndim_atl      ) ! Lat-Depth
+               CALL wheneq ( jpj*jpk, MIN(surf_jk_atl(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_atl_30  , ndim_atl_30   ) ! Lat-Depth
+               CALL wheneq ( jpj    , MIN(surf_jk_atl(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_atl_30, ndim_h_atl_30 ) ! Lat
+
+               CALL wheneq ( jpj*jpk, MIN(surf_jk_pac(:,:)         , 1.e0), 1, 1., ndex_pac     , ndim_pac      ) ! Lat-Depth
+               CALL wheneq ( jpj*jpk, MIN(surf_jk_pac(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_pac_30  , ndim_pac_30   ) ! Lat-Depth
+               CALL wheneq ( jpj    , MIN(surf_jk_pac(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_pac_30, ndim_h_pac_30 ) ! Lat
+
+               CALL wheneq ( jpj*jpk, MIN(surf_jk_ind(:,:)         , 1.e0), 1, 1., ndex_ind     , ndim_ind      ) ! Lat-Depth
+               CALL wheneq ( jpj*jpk, MIN(surf_jk_ind(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_ind_30  , ndim_ind_30   ) ! Lat-Depth
+               CALL wheneq ( jpj    , MIN(surf_jk_ind(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_ind_30, ndim_h_ind_30 ) ! Lat
+
+               CALL wheneq ( jpj*jpk, MIN(surf_jk_ipc(:,:)         , 1.e0), 1, 1., ndex_ipc     , ndim_ipc      ) ! Lat-Depth
+               CALL wheneq ( jpj*jpk, MIN(surf_jk_ipc(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_ipc_30  , ndim_ipc_30   ) ! Lat-Depth
+               CALL wheneq ( jpj    , MIN(surf_jk_ipc(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_ipc_30, ndim_h_ipc_30 ) ! Lat
+
+            ENDIF
+
+            ! 
+#if defined key_diaeiv
+            cl_comment = ' (Bolus part included)'
 #else
-         cl_comment = ' '
-#endif
-         !  Zonal mean T and S
-
-         IF ( ln_diaznl ) THEN 
-            CALL histdef( numptr, "zotemglo", "Zonal Mean Temperature","C" ,   &
-               1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-            CALL histdef( numptr, "zosalglo", "Zonal Mean Salinity","PSU"  ,   &
-               1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-
-            CALL histdef( numptr, "zosrfglo", "Zonal Mean Surface","m^2"   ,   &
-               1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
+            cl_comment = '                      '
+#endif
+            !  Zonal mean T and S
+
+            IF ( ln_diaznl ) THEN 
+               CALL histdef( numptr, "zotemglo", "Zonal Mean Temperature","C" ,   &
+                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+               CALL histdef( numptr, "zosalglo", "Zonal Mean Salinity","PSU"  ,   &
+                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+
+               CALL histdef( numptr, "zosrfglo", "Zonal Mean Surface","m^2"   ,   &
+                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
+
+               IF (ln_subbas) THEN 
+                  CALL histdef( numptr, "zotematl", "Zonal Mean Temperature: Atlantic","C" ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+                  CALL histdef( numptr, "zosalatl", "Zonal Mean Salinity: Atlantic","PSU"  ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+                  CALL histdef( numptr, "zosrfatl", "Zonal Mean Surface: Atlantic","m^2"   ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
+
+                  CALL histdef( numptr, "zotempac", "Zonal Mean Temperature: Pacific","C"  ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+                  CALL histdef( numptr, "zosalpac", "Zonal Mean Salinity: Pacific","PSU"   ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+                  CALL histdef( numptr, "zosrfpac", "Zonal Mean Surface: Pacific","m^2"    ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
+
+                  CALL histdef( numptr, "zotemind", "Zonal Mean Temperature: Indian","C"   ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+                  CALL histdef( numptr, "zosalind", "Zonal Mean Salinity: Indian","PSU"    ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+                  CALL histdef( numptr, "zosrfind", "Zonal Mean Surface: Indian","m^2"     ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
+
+                  CALL histdef( numptr, "zotemipc", "Zonal Mean Temperature: Pacific+Indian","C" ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+                  CALL histdef( numptr, "zosalipc", "Zonal Mean Salinity: Pacific+Indian","PSU"  ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
+                  CALL histdef( numptr, "zosrfipc", "Zonal Mean Surface: Pacific+Indian","m^2"   ,   &
+                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
+               ENDIF
+
+            ENDIF
+
+            !  Meridional Stream-Function (Eulerian and Bolus)
+
+            CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: Global"//TRIM(cl_comment),"Sv" ,   &
+               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
+            IF( ln_subbas .AND. ln_diaznl ) THEN
+               CALL histdef( numptr, "zomsfatl", "Meridional Stream-Function: Atlantic"//TRIM(cl_comment),"Sv" ,   &
+                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
+               CALL histdef( numptr, "zomsfpac", "Meridional Stream-Function: Pacific"//TRIM(cl_comment),"Sv"  ,   &
+                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
+               CALL histdef( numptr, "zomsfind", "Meridional Stream-Function: Indian"//TRIM(cl_comment),"Sv"   ,   &
+                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
+               CALL histdef( numptr, "zomsfipc", "Meridional Stream-Function: Indo-Pacific"//TRIM(cl_comment),"Sv" ,&
+                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
+            ENDIF
+
+            !  Heat transport 
+
+            CALL histdef( numptr, "sophtadv", "Advective Heat Transport"      ,   &
+               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            CALL histdef( numptr, "sophtldf", "Diffusive Heat Transport"      ,   &
+               "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            IF ( ln_ptrcomp ) THEN 
+               CALL histdef( numptr, "sophtove", "Overturning Heat Transport"    ,   &
+                  "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            END IF
+            IF( ln_subbas ) THEN
+               CALL histdef( numptr, "sohtatl", "Heat Transport Atlantic"//TRIM(cl_comment),  &
+                  "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+               CALL histdef( numptr, "sohtpac", "Heat Transport Pacific"//TRIM(cl_comment) ,  &
+                  "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+               CALL histdef( numptr, "sohtind", "Heat Transport Indian"//TRIM(cl_comment)  ,  &
+                  "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+               CALL histdef( numptr, "sohtipc", "Heat Transport Pacific+Indian"//TRIM(cl_comment), &
+                  "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            ENDIF
+
+
+            !  Salt transport 
+
+            CALL histdef( numptr, "sopstadv", "Advective Salt Transport"      ,   &
+               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            CALL histdef( numptr, "sopstldf", "Diffusive Salt Transport"      ,   &
+               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            IF ( ln_ptrcomp ) THEN 
+               CALL histdef( numptr, "sopstove", "Overturning Salt Transport"    ,   &
+                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            END IF
+#if defined key_diaeiv
+            ! Eddy induced velocity
+            CALL histdef( numptr, "zomsfeiv", "Bolus Meridional Stream-Function: global",   &
+               "Sv"      , 1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
+            CALL histdef( numptr, "sophteiv", "Bolus Advective Heat Transport",   &
+               "PW"      , 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            CALL histdef( numptr, "sopsteiv", "Bolus Advective Salt Transport",   &
+               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+#endif
+            IF( ln_subbas ) THEN
+               CALL histdef( numptr, "sostatl", "Salt Transport Atlantic"//TRIM(cl_comment)      ,  &
+                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+               CALL histdef( numptr, "sostpac", "Salt Transport Pacific"//TRIM(cl_comment)      ,   &
+                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+               CALL histdef( numptr, "sostind", "Salt Transport Indian"//TRIM(cl_comment)      ,    &
+                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+               CALL histdef( numptr, "sostipc", "Salt Transport Pacific+Indian"//TRIM(cl_comment),  &
+                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            ENDIF
+
+            CALL histend( numptr )
+
+         END IF
+      END IF
+
+      IF( MOD( itmod, nf_ptr ) == 0 .AND. l_znl_root ) THEN
+
+         niter = niter + 1
+
+!-$$         IF(lwp) THEN
+!-$$            WRITE(numout,*)
+!-$$            WRITE(numout,*) 'dia_ptr_wri : write Poleward Transports at time-step : kt = ', kt, &
+!-$$               & 'it = ', it, ' itmod = ', itmod, ' niter = ', niter
+!-$$            WRITE(numout,*) '~~~~~~~~~~'
+!-$$            WRITE(numout,*)
+!-$$         ENDIF
+
+         IF (ln_diaznl ) THEN 
+            CALL histwrite( numptr, "zosrfglo", niter, surf_jk_glo , ndim, ndex )
+            CALL histwrite( numptr, "zotemglo", niter, tn_jk_glo  , ndim, ndex )
+            CALL histwrite( numptr, "zosalglo", niter, sn_jk_glo  , ndim, ndex )
 
             IF (ln_subbas) THEN 
-               CALL histdef( numptr, "zotematl", "Zonal Mean Temperature: Atlantic","C" ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-               CALL histdef( numptr, "zosalatl", "Zonal Mean Salinity: Atlantic","PSU"  ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-               CALL histdef( numptr, "zosrfatl", "Zonal Mean Surface: Atlantic","m^2"   ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
-
-               CALL histdef( numptr, "zotempac", "Zonal Mean Temperature: Pacific","C"  ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-               CALL histdef( numptr, "zosalpac", "Zonal Mean Salinity: Pacific","PSU"   ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-               CALL histdef( numptr, "zosrfpac", "Zonal Mean Surface: Pacific","m^2"    ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
-
-               CALL histdef( numptr, "zotemind", "Zonal Mean Temperature: Indian","C"   ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-               CALL histdef( numptr, "zosalind", "Zonal Mean Salinity: Indian","PSU"    ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-               CALL histdef( numptr, "zosrfind", "Zonal Mean Surface: Indian","m^2"     ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
-
-               CALL histdef( numptr, "zotemipc", "Zonal Mean Temperature: Pacific+Indian","C" ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-               CALL histdef( numptr, "zosalipc", "Zonal Mean Salinity: Pacific+Indian","PSU"  ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
-               CALL histdef( numptr, "zosrfipc", "Zonal Mean Surface: Pacific+Indian","m^2"   ,   &
-                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
-            ENDIF
-
-         ENDIF
-
-         !  Meridional Stream-Function (Eulerian and Bolus)
-
-         CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: Global"//TRIM(cl_comment),"Sv" ,   &
-            1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
-         IF( ln_subbas .AND. ln_diaznl ) THEN
-            CALL histdef( numptr, "zomsfatl", "Meridional Stream-Function: Atlantic"//TRIM(cl_comment),"Sv" ,   &
-               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
-            CALL histdef( numptr, "zomsfpac", "Meridional Stream-Function: Pacific"//TRIM(cl_comment),"Sv"  ,   &
-               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
-            CALL histdef( numptr, "zomsfind", "Meridional Stream-Function: Indian"//TRIM(cl_comment),"Sv"   ,   &
-               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
-            CALL histdef( numptr, "zomsfipc", "Meridional Stream-Function: Indo-Pacific"//TRIM(cl_comment),"Sv" ,&
-               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
-         ENDIF
-
-         !  Heat transport 
-
-         CALL histdef( numptr, "sophtadv", "Advective Heat Transport"      ,   &
-            "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-         CALL histdef( numptr, "sophtldf", "Diffusive Heat Transport"      ,   &
-            "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-         CALL histdef( numptr, "sophtove", "Overturning Heat Transport"    ,   &
-            "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-         IF( ln_subbas ) THEN
-            CALL histdef( numptr, "sohtatl", "Heat Transport Atlantic"//TRIM(cl_comment),  &
-               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-            CALL histdef( numptr, "sohtpac", "Heat Transport Pacific"//TRIM(cl_comment) ,  &
-               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-            CALL histdef( numptr, "sohtind", "Heat Transport Indian"//TRIM(cl_comment)  ,  &
-               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-            CALL histdef( numptr, "sohtipc", "Heat Transport Pacific+Indian"//TRIM(cl_comment), &
-               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-         ENDIF
-
-
-         !  Salt transport 
-
-         CALL histdef( numptr, "sopstadv", "Advective Salt Transport"      ,   &
-            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-         CALL histdef( numptr, "sopstldf", "Diffusive Salt Transport"      ,   &
-            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-         CALL histdef( numptr, "sopstove", "Overturning Salt Transport"    ,   &
-            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-
-#if defined key_diaeiv
-         ! Eddy induced velocity
-         CALL histdef( numptr, "zomsfeiv", "Bolus Meridional Stream-Function: global",   &
-            "Sv"      , 1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
-         CALL histdef( numptr, "sophteiv", "Bolus Advective Heat Transport",   &
-            "PW"      , 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-         CALL histdef( numptr, "sopsteiv", "Bolus Advective Salt Transport",   &
-            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-#endif
-         IF( ln_subbas ) THEN
-            CALL histdef( numptr, "sostatl", "Salt Transport Atlantic"//TRIM(cl_comment)      ,  &
-               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-            CALL histdef( numptr, "sostpac", "Salt Transport Pacific"//TRIM(cl_comment)      ,   &
-               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-            CALL histdef( numptr, "sostind", "Salt Transport Indian"//TRIM(cl_comment)      ,    &
-               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-            CALL histdef( numptr, "sostipc", "Salt Transport Pacific+Indian"//TRIM(cl_comment),  &
-               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-         ENDIF
-
-         CALL histend( numptr )
-
-      ENDIF
-
-      IF( MOD( itmod, nf_ptr ) == 0 ) THEN
-
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'dia_ptr : compute Poleward Transports at time-step : ', kt
-            WRITE(numout,*) '~~~~~~~~'
-            WRITE(numout,*)
-         ENDIF
-
-         IF (ln_diaznl ) THEN 
-            CALL histwrite( numptr, "zosrfglo", it, surf_jk_glo , ndim, ndex )
-            CALL histwrite( numptr, "zotemglo", it, tn_jk_glo  , ndim, ndex )
-            CALL histwrite( numptr, "zosalglo", it, sn_jk_glo  , ndim, ndex )
-
-            IF (ln_subbas) THEN 
-               CALL histwrite( numptr, "zosrfatl", it, surf_jk_atl, ndim_atl, ndex_atl )
-               CALL histwrite( numptr, "zosrfpac", it, surf_jk_pac, ndim_pac, ndex_pac )
-               CALL histwrite( numptr, "zosrfind", it, surf_jk_ind, ndim_ind, ndex_ind )
-               CALL histwrite( numptr, "zosrfipc", it, surf_jk_ipc, ndim_ipc, ndex_ipc )
-
-               CALL histwrite( numptr, "zotematl", it, tn_jk_atl  , ndim_atl, ndex_atl )
-               CALL histwrite( numptr, "zosalatl", it, sn_jk_atl  , ndim_atl, ndex_atl )
-               CALL histwrite( numptr, "zotempac", it, tn_jk_pac  , ndim_pac, ndex_pac )
-               CALL histwrite( numptr, "zosalpac", it, sn_jk_pac  , ndim_pac, ndex_pac )
-               CALL histwrite( numptr, "zotemind", it, tn_jk_ind  , ndim_ind, ndex_ind )
-               CALL histwrite( numptr, "zosalind", it, sn_jk_ind  , ndim_ind, ndex_ind )
-               CALL histwrite( numptr, "zotemipc", it, tn_jk_ipc  , ndim_ipc, ndex_ipc )
-               CALL histwrite( numptr, "zosalipc", it, sn_jk_ipc  , ndim_ipc, ndex_ipc )
+               CALL histwrite( numptr, "zosrfatl", niter, surf_jk_atl, ndim_atl, ndex_atl )
+               CALL histwrite( numptr, "zosrfpac", niter, surf_jk_pac, ndim_pac, ndex_pac )
+               CALL histwrite( numptr, "zosrfind", niter, surf_jk_ind, ndim_ind, ndex_ind )
+               CALL histwrite( numptr, "zosrfipc", niter, surf_jk_ipc, ndim_ipc, ndex_ipc )
+
+               CALL histwrite( numptr, "zotematl", niter, tn_jk_atl  , ndim_atl, ndex_atl )
+               CALL histwrite( numptr, "zosalatl", niter, sn_jk_atl  , ndim_atl, ndex_atl )
+               CALL histwrite( numptr, "zotempac", niter, tn_jk_pac  , ndim_pac, ndex_pac )
+               CALL histwrite( numptr, "zosalpac", niter, sn_jk_pac  , ndim_pac, ndex_pac )
+               CALL histwrite( numptr, "zotemind", niter, tn_jk_ind  , ndim_ind, ndex_ind )
+               CALL histwrite( numptr, "zosalind", niter, sn_jk_ind  , ndim_ind, ndex_ind )
+               CALL histwrite( numptr, "zotemipc", niter, tn_jk_ipc  , ndim_ipc, ndex_ipc )
+               CALL histwrite( numptr, "zosalipc", niter, sn_jk_ipc  , ndim_ipc, ndex_ipc )
             END IF
          ENDIF
 
          ! overturning outputs:
-         CALL histwrite( numptr, "zomsfglo", it, v_msf_glo, ndim, ndex )
+         CALL histwrite( numptr, "zomsfglo", niter, v_msf_glo, ndim, ndex )
          IF( ln_subbas .AND. ln_diaznl ) THEN
-            CALL histwrite( numptr, "zomsfatl", it, v_msf_atl , ndim_atl_30, ndex_atl_30 )
-            CALL histwrite( numptr, "zomsfpac", it, v_msf_pac , ndim_pac_30, ndex_pac_30 )
-            CALL histwrite( numptr, "zomsfind", it, v_msf_ind , ndim_ind_30, ndex_ind_30 )
-            CALL histwrite( numptr, "zomsfipc", it, v_msf_ipc , ndim_ipc_30, ndex_ipc_30 )
+            CALL histwrite( numptr, "zomsfatl", niter, v_msf_atl , ndim_atl_30, ndex_atl_30 )
+            CALL histwrite( numptr, "zomsfpac", niter, v_msf_pac , ndim_pac_30, ndex_pac_30 )
+            CALL histwrite( numptr, "zomsfind", niter, v_msf_ind , ndim_ind_30, ndex_ind_30 )
+            CALL histwrite( numptr, "zomsfipc", niter, v_msf_ipc , ndim_ipc_30, ndex_ipc_30 )
          ENDIF
 #if defined key_diaeiv
-         CALL histwrite( numptr, "zomsfeiv", it, v_msf_eiv, ndim  , ndex   )
-#endif
-
+         CALL histwrite( numptr, "zomsfeiv", niter, v_msf_eiv_glo, ndim  , ndex   )
+#endif
+
+
+         ! heat transport outputs:
+         IF( ln_subbas ) THEN
+            CALL histwrite( numptr, "sohtatl", niter, ht_atl  , ndim_h_atl_30, ndex_h_atl_30 )
+            CALL histwrite( numptr, "sohtpac", niter, ht_pac  , ndim_h_pac_30, ndex_h_pac_30 )
+            CALL histwrite( numptr, "sohtind", niter, ht_ind  , ndim_h_ind_30, ndex_h_ind_30 )
+            CALL histwrite( numptr, "sohtipc", niter, ht_ipc  , ndim_h_ipc_30, ndex_h_ipc_30 )
+            CALL histwrite( numptr, "sostatl", niter, st_atl  , ndim_h_atl_30, ndex_h_atl_30 )
+            CALL histwrite( numptr, "sostpac", niter, st_pac  , ndim_h_pac_30, ndex_h_pac_30 )
+            CALL histwrite( numptr, "sostind", niter, st_ind  , ndim_h_ind_30, ndex_h_ind_30 )
+            CALL histwrite( numptr, "sostipc", niter, st_ipc  , ndim_h_ipc_30, ndex_h_ipc_30 )
+         ENDIF
+
+         CALL histwrite( numptr, "sophtadv", niter, pht_adv     , ndim_h, ndex_h )
+         CALL histwrite( numptr, "sophtldf", niter, pht_ldf     , ndim_h, ndex_h )
+         CALL histwrite( numptr, "sopstadv", niter, pst_adv     , ndim_h, ndex_h )
+         CALL histwrite( numptr, "sopstldf", niter, pst_ldf     , ndim_h, ndex_h )
+         IF ( ln_ptrcomp ) THEN 
+            CALL histwrite( numptr, "sopstove", niter, pst_ove_glo , ndim_h, ndex_h )
+            CALL histwrite( numptr, "sophtove", niter, pht_ove_glo , ndim_h, ndex_h )
+         ENDIF
+#if defined key_diaeiv
+         CALL histwrite( numptr, "sophteiv", niter, pht_eiv_glo  , ndim_h, ndex_h )
+         CALL histwrite( numptr, "sopsteiv", niter, pst_eiv_glo  , ndim_h, ndex_h )
+#endif
+         !
       ENDIF
-
-      ! heat transport outputs:
-      IF( ln_subbas ) THEN
-         CALL histwrite( numptr, "sohtatl", it, ht_atl  , ndim_h_atl_30, ndex_h_atl_30 )
-         CALL histwrite( numptr, "sohtpac", it, ht_pac  , ndim_h_pac_30, ndex_h_pac_30 )
-         CALL histwrite( numptr, "sohtind", it, ht_ind  , ndim_h_ind_30, ndex_h_ind_30 )
-         CALL histwrite( numptr, "sohtipc", it, ht_ipc  , ndim_h_ipc_30, ndex_h_ipc_30 )
-         CALL histwrite( numptr, "sostatl", it, st_atl  , ndim_h_atl_30, ndex_h_atl_30 )
-         CALL histwrite( numptr, "sostpac", it, st_pac  , ndim_h_pac_30, ndex_h_pac_30 )
-         CALL histwrite( numptr, "sostind", it, st_ind  , ndim_h_ind_30, ndex_h_ind_30 )
-         CALL histwrite( numptr, "sostipc", it, st_ipc  , ndim_h_ipc_30, ndex_h_ipc_30 )
-      ENDIF
-
-      CALL histwrite( numptr, "sophtadv", it, pht_adv  , ndim_h, ndex_h )
-      CALL histwrite( numptr, "sophtldf", it, pht_ldf  , ndim_h, ndex_h )
-      CALL histwrite( numptr, "sophtove", it, pht_ove  , ndim_h, ndex_h )
-      CALL histwrite( numptr, "sopstadv", it, pst_adv  , ndim_h, ndex_h )
-      CALL histwrite( numptr, "sopstldf", it, pst_ldf  , ndim_h, ndex_h )
-      CALL histwrite( numptr, "sopstove", it, pst_ove  , ndim_h, ndex_h )
-#if defined key_diaeiv
-      CALL histwrite( numptr, "sophteiv", it, pht_eiv  , ndim_h, ndex_h )
-      CALL histwrite( numptr, "sopsteiv", it, pst_eiv  , ndim_h, ndex_h )
-#endif
-
       !
    END SUBROUTINE dia_ptr_wri
Index: /trunk/NEMO/OPA_SRC/lib_mpp.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/lib_mpp.F90	(revision 1344)
+++ /trunk/NEMO/OPA_SRC/lib_mpp.F90	(revision 1345)
@@ -17,4 +17,5 @@
    !!             -   !  2008  (R. Benshila) add mpp_ini_ice
    !!            3.2  !  2009  (R. Benshila) SHMEM suppression, north fold in lbc_nfd
+   !!            3.2  !  2009  (O. Marti)    add mpp_ini_znl 
    !!----------------------------------------------------------------------
 #if   defined key_mpp_mpi  
@@ -27,5 +28,5 @@
    !!   mpp_lnk_e   : interface (defined in lbclnk) for message passing of 2d array with extra halo (mpp_lnk_2d_e)
    !!   mpprecv     :
-   !!   mppsend     :
+   !!   mppsend     :   SUBROUTINE mpp_ini_znl
    !!   mppscatter  :
    !!   mppgather   :
@@ -70,5 +71,5 @@
    PUBLIC   mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e
    PUBLIC   mpprecv, mppsend, mppscatter, mppgather
-   PUBLIC   mppobc, mpp_ini_ice, mpp_isl
+   PUBLIC   mppobc, mpp_ini_ice, mpp_isl, mpp_ini_znl
 #if defined key_oasis3 || defined key_oasis4
    PUBLIC   mppsize, mpprank
@@ -119,12 +120,20 @@
 !!gm question : Pourquoi toutes les variables ice sont public???
    ! variables used in case of sea-ice
-   INTEGER, PUBLIC ::   ngrp_ice        !: group ID for the ice processors (for rheology)
    INTEGER, PUBLIC ::   ncomm_ice       !: communicator made by the processors with sea-ice
-   INTEGER, PUBLIC ::   ndim_rank_ice   !: number of 'ice' processors
-   INTEGER, PUBLIC ::   n_ice_root      !: number (in the comm_ice) of proc 0 in the ice comm
+   INTEGER ::   ngrp_ice        !  group ID for the ice processors (for rheology)
+   INTEGER ::   ndim_rank_ice   !  number of 'ice' processors
+   INTEGER ::   n_ice_root      !  number (in the comm_ice) of proc 0 in the ice comm
    INTEGER, DIMENSION(:), ALLOCATABLE ::   nrank_ice     ! dimension ndim_rank_ice
+
+   ! variables used for zonal integration
+   INTEGER, PUBLIC ::   ncomm_znl       !: communicator made by the processors on the same zonal average
+   LOGICAL, PUBLIC ::   l_znl_root      ! True on the 'left'most processor on the same row
+   INTEGER ::   ngrp_znl        ! group ID for the znl processors
+   INTEGER ::   ndim_rank_znl   ! number of processors on the same zonal average
+   INTEGER, DIMENSION(:), ALLOCATABLE ::   nrank_znl  ! dimension ndim_rank_znl, number of the procs into the same znl domain
    
    ! North fold condition in mpp_mpi with jpni > 1
    INTEGER ::   ngrp_world        ! group ID for the world processors
+   INTEGER ::   ngrp_opa          ! group ID for the opa processors
    INTEGER ::   ngrp_north        ! group ID for the northern processors (to be fold)
    INTEGER ::   ncomm_north       ! communicator made by the processors belonging to ngrp_north
@@ -355,5 +364,5 @@
             t3we(:,jl,:,1) = ptab(iihom +jl,:,:)
          END DO
-      END SELECT
+      END SELECT  
       !
       !                           ! Migrations
@@ -1244,5 +1253,5 @@
 
 
-   SUBROUTINE mppmin_int( ktab )
+   SUBROUTINE mppmin_int( ktab, kcom )
       !!----------------------------------------------------------------------
       !!                  ***  routine mppmin_int  ***
@@ -1252,9 +1261,13 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(inout) ::   ktab      ! ???
-      !!
-      INTEGER ::  ierror, iwork
-      !!----------------------------------------------------------------------
-      !
-      CALL mpi_allreduce( ktab, iwork, 1, mpi_integer, mpi_min, mpi_comm_opa, ierror )
+      INTEGER , INTENT( in  ), OPTIONAL        ::   kcom        ! input array
+      !!
+      INTEGER ::  ierror, iwork, localcomm
+      !!----------------------------------------------------------------------
+      !
+      localcomm = mpi_comm_opa
+      IF( PRESENT(kcom) )   localcomm = kcom
+      !
+     CALL mpi_allreduce( ktab, iwork, 1, mpi_integer, mpi_min, localcomm, ierror )
       !
       ktab = iwork
@@ -1983,4 +1996,98 @@
 
 
+   SUBROUTINE mpp_ini_znl
+      !!----------------------------------------------------------------------
+      !!               ***  routine mpp_ini_znl  ***
+      !!
+      !! ** Purpose :   Initialize special communicator for computing zonal sum
+      !!
+      !! ** Method  : - Look for processors in the same row
+      !!              - Put their number in nrank_znl
+      !!              - Create group for the znl processors
+      !!              - Create a communicator for znl processors
+      !!              - Determine if processor should write znl files
+      !!
+      !! ** output
+      !!      ndim_rank_znl = number of processors on the same row
+      !!      ngrp_znl = group ID for the znl processors
+      !!      ncomm_znl = communicator for the ice procs.
+      !!      n_znl_root = number (in the world) of proc 0 in the ice comm.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr
+      INTEGER :: jproc
+      INTEGER :: ii
+      INTEGER, DIMENSION(jpnij) :: kwork
+      !
+      !-$$     WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ngrp_world     : ', ngrp_world
+      !-$$     WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - mpi_comm_world : ', mpi_comm_world
+      !-$$     WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - mpi_comm_opa   : ', mpi_comm_opa
+      !
+      IF ( jpnj == 1 ) THEN
+         ngrp_znl  = ngrp_world
+         ncomm_znl = mpi_comm_opa
+      ELSE
+         !
+         CALL MPI_ALLGATHER ( njmpp, 1, mpi_integer, kwork, 1, mpi_integer, mpi_comm_opa, ierr )
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - kwork pour njmpp : ', kwork
+         !-$$        CALL flush(numout)
+         !
+         ! Count number of processors on the same row
+         ndim_rank_znl = 0
+         DO jproc=1,jpnij
+            IF ( kwork(jproc) == njmpp ) THEN
+               ndim_rank_znl = ndim_rank_znl + 1
+            ENDIF
+         END DO
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ndim_rank_znl : ', ndim_rank_znl
+         !-$$        CALL flush(numout)
+         ! Allocate the right size to nrank_znl
+#if ! defined key_agrif
+         IF (ALLOCATED(nrank_znl)) DEALLOCATE(nrank_znl)
+#else
+         DEALLOCATE(nrank_znl)
+#endif
+         ALLOCATE(nrank_znl(ndim_rank_znl))
+         ii = 0     
+         nrank_znl (:) = 0
+         DO jproc=1,jpnij
+            IF ( kwork(jproc) == njmpp) THEN
+               ii = ii + 1
+               nrank_znl(ii) = jproc -1 
+            ENDIF
+         END DO
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - nrank_znl : ', nrank_znl
+         !-$$        CALL flush(numout)
+
+         ! Create the opa group
+         CALL MPI_COMM_GROUP(mpi_comm_opa,ngrp_opa,ierr)
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ngrp_opa : ', ngrp_opa
+         !-$$        CALL flush(numout)
+
+         ! Create the znl group from the opa group
+         CALL MPI_GROUP_INCL  ( ngrp_opa, ndim_rank_znl, nrank_znl, ngrp_znl, ierr )
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ngrp_znl ', ngrp_znl
+         !-$$        CALL flush(numout)
+
+         ! Create the znl communicator from the opa communicator, ie the pool of procs in the same row
+         CALL MPI_COMM_CREATE ( mpi_comm_opa, ngrp_znl, ncomm_znl, ierr )
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ncomm_znl ', ncomm_znl
+         !-$$        CALL flush(numout)
+         !
+      END IF
+
+      ! Determines if processor if the first (starting from i=1) on the row
+      IF ( jpni == 1 ) THEN 
+         l_znl_root = .TRUE.
+      ELSE
+         l_znl_root = .FALSE.
+         kwork (1) = nimpp
+         CALL mpp_min ( kwork(1), kcom = ncomm_znl)
+         IF ( nimpp == kwork(1)) l_znl_root = .TRUE.
+      END IF
+
+   END SUBROUTINE mpp_ini_znl
+
+
    SUBROUTINE mpp_ini_north
       !!----------------------------------------------------------------------
@@ -2493,4 +2600,9 @@
    END SUBROUTINE mpp_ini_ice
 
+   SUBROUTINE mpp_ini_znl
+      INTEGER :: kcom
+      WRITE(*,*) 'mpp_ini_znl: You should not have seen this print! error?'
+   END SUBROUTINE mpp_ini_znl
+
    SUBROUTINE mpp_comm_free( kcom )
       INTEGER :: kcom
