Index: /trunk/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- /trunk/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 1339)
+++ /trunk/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 1340)
@@ -680,7 +680,9 @@
 &namptr       !   Poleward Transport Diagnostic
 !-----------------------------------------------------------------------
-   ln_diaptr  = .false.    !  Poleward heat and salt transport (T) or not (F)
-   ln_subbas  = .false.    !  Atlantic/Pacific/Indian basins computation (T) or not 
+   ln_diaptr  = .false.     !  Poleward heat and salt transport (T) or not (F)
+   ln_diaznl  = .false.     !  Add zonal means and meridional stream functions
+   ln_subbas  = .false.     !  Atlantic/Pacific/Indian basins computation (T) or not 
                            !  (orca configuration only, need input basins mask file named "subbasins.nc"
-   nf_ptr     =  15        !  Frequency of ptr computation [time step]
-/
+   nf_ptr     =  1         !  Frequency of ptr computation [time step]
+   nf_ptr_wri =  15        !  Frequency of ptr outputs
+/
Index: /trunk/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist
===================================================================
--- /trunk/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 1339)
+++ /trunk/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 1340)
@@ -680,7 +680,9 @@
 &namptr       !   Poleward Transport Diagnostic
 !-----------------------------------------------------------------------
-   ln_diaptr  = .false.    !  Poleward heat and salt transport (T) or not (F)
-   ln_subbas  = .false.    !  Atlantic/Pacific/Indian basins computation (T) or not 
+   ln_diaptr  = .false.     !  Poleward heat and salt transport (T) or not (F)
+   ln_diaznl  = .false.     !  Add zonal means and meridional stream functions
+   ln_subbas  = .false.     !  Atlantic/Pacific/Indian basins computation (T) or not 
                            !  (orca configuration only, need input basins mask file named "subbasins.nc"
-   nf_ptr     =  15        !  Frequency of ptr computation [time step]
-/
+   nf_ptr     =  1         !  Frequency of ptr computation [time step]
+   nf_ptr_wri =  15        !  Frequency of ptr outputs
+/
Index: /trunk/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 1339)
+++ /trunk/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 1340)
@@ -2,8 +2,7 @@
    !!======================================================================
    !!                       ***  MODULE  diaptr  ***
-   !! Ocean physics:  brief description of the purpose of the module
-   !!                 (please no more than 2 lines)
+   !! Ocean physics:  Computes meridonal transports and zonal means
    !!=====================================================================
-   !! History :  9.0  !  03-09  (C. Talandir, G. Madec)  Original code
+   !! History :  9.0  !  03-09  (C. Talandier, G. Madec)  Original code
    !!            9.0  !  06-01  (A. Biastoch)  Allow sub-basins computation
    !!----------------------------------------------------------------------
@@ -14,5 +13,5 @@
    !!   dia_ptr_wri  : Output of poleward fluxes
    !!   ptr_vjk      : "zonal" sum computation of a "meridional" flux array
-   !!   ptr_vtjk     : "zonal" mean computation of a tracer field
+   !!   ptr_tjk      : "zonal" mean computation of a tracer field
    !!   ptr_vj       : "zonal" and vertical sum computation of a "meridional"
    !!                : flux array; Generic interface: ptr_vj_3d, ptr_vj_2d
@@ -20,5 +19,5 @@
    USE oce           ! ocean dynamics and active tracers
    USE dom_oce       ! ocean space and time domain
-   USE ldftra_oce    ! ???
+   USE ldftra_oce    ! ocean active tracers: lateral physics
    USE lib_mpp
    USE in_out_manager
@@ -41,8 +40,10 @@
    PUBLIC   ptr_vjk        ! call by tra_ldf & tra_adv routines
 
-   !!! ** init namelist (namptr)
+!!! ** init namelist (namptr)
    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
    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
@@ -52,18 +53,33 @@
    REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_eiv, pst_eiv      !: heat and salt poleward transport: bolus advection
 #endif
-   REAL(wp), PUBLIC, DIMENSION(jpj) ::   ht_atl,ht_ind,ht_pac  !: heat
-   REAL(wp), PUBLIC, DIMENSION(jpj) ::   st_atl,st_ind,st_pac  !: salt
-
-     
-
-   REAL(wp), DIMENSION(jpj,jpk) ::   tn_jk  , sn_jk  ,  &  !: "zonal" mean temperature and salinity
-      &                              v_msf_atl       ,  &  !: "meridional" Stream-Function
+   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
+
+   REAL(wp), DIMENSION(jpj,jpk) ::   tn_jk_glo, sn_jk_glo,  &  !: "zonal" mean temperature and salinity
+      &                              tn_jk_atl, sn_jk_atl,  &
+      &                              tn_jk_pac, sn_jk_pac,  &
+      &                              tn_jk_ind, sn_jk_ind,  &
+      &                              tn_jk_ipc, sn_jk_ipc,  &
       &                              v_msf_glo       ,  &  !: "meridional" Stream-Function
-      &                              v_msf_ipc       ,  &  !: "meridional" Stream-Function
-      &                              surf_jk_r             !: inverse of the ocean "zonal" section surface
-#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     !: return function value
+      &                              v_msf_atl       ,  &  
+      &                              v_msf_pac       ,  &  
+      &                              v_msf_ind       ,  &  
+      &                              v_msf_ipc       ,  &  
+      &                              surf_jk_glo     ,  &  !: Ocean "zonal" section surface
+      &                              surf_jk_atl     ,  &  
+      &                              surf_jk_pac     ,  &  
+      &                              surf_jk_ind     ,  &         
+      &                              surf_jk_ipc     ,  & 
+      &                              surf_jk_r_glo   ,  &  !: inverse of the ocean "zonal" section surface
+      &                              surf_jk_r_atl   ,  &  
+      &                              surf_jk_r_pac   ,  &  
+      &                              surf_jk_r_ind   ,  &         
+      &                              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
 
    !! * Substitutions
@@ -144,5 +160,5 @@
 
 
-   FUNCTION ptr_vjk( pva )   RESULT ( p_fval )
+   FUNCTION ptr_vjk( pva, bmask )   RESULT ( p_fval )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ptr_vjk  ***
@@ -156,4 +172,5 @@
       !!----------------------------------------------------------------------
       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pva   ! mask flux array at V-point
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL :: bmask ! Optional 2D basin mask
       !!
       INTEGER                      ::   ji, jj, jk   ! dummy loop arguments
@@ -166,12 +183,24 @@
       p_fval(:,:) = 0.e0
       !
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-           DO ji = fs_2, fs_jpim1
-            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)
-           END DO
-         END DO
-      END DO
+      IF (PRESENT (bmask)) THEN 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  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) &
+                     &                                          * bmask(ji,jj)
+               END DO
+            END DO
+         END DO
+      ELSE 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  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)
+               END DO
+            END DO
+         END DO
+      END IF
       !
       IF(lk_mpp) THEN
@@ -184,18 +213,18 @@
    END FUNCTION ptr_vjk
 
-
-   FUNCTION ptr_vtjk( pva )   RESULT ( p_fval )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE ptr_vtjk  ***
+   FUNCTION ptr_tjk( pta, bmask )   RESULT ( p_fval )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE ptr_tjk  ***
       !!
       !! ** Purpose :   "zonal" mean computation of a tracer field
       !!
-      !! ** Method  : - i-sum of mj(pva) using the interior 2D vmask (vmask_i)
+      !! ** Method  : - i-sum of mj(pta) using tmask
       !!      multiplied by the inverse of the surface of the "zonal" ocean
       !!      section
       !!
-      !! ** 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
+      !! ** Action  : - p_fval: i-k-mean poleward flux of pta
+      !!----------------------------------------------------------------------
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! tracer flux array at T-point
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL :: bmask ! Optional 2D basin mask
       !!
       INTEGER                     ::   ji, jj, jk   ! dummy loop arguments
@@ -207,13 +236,26 @@
       !
       p_fval(:,:) = 0.e0
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! Vector opt.
-               p_fval(jj,jk) = p_fval(jj,jk) + ( pva(ji,jj,jk) + pva(ji,jj+1,jk) )              &
-                  &                          * e1v(ji,jj) * fse3v(ji,jj,jk) * vmask(ji,jj,jk)   &
-                  &                          * tmask_i(ji,jj+1) * tmask_i(ji,jj)
-            END DO
-         END DO
-      END DO
+      IF (PRESENT (bmask)) THEN 
+         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)                  &
+                     &                          * e1t(ji,jj) * fse3t(ji,jj,jk)   &
+                     &                          * tmask_i(ji,jj)                 &
+                     &                          * bmask(ji,jj)
+               END DO
+            END DO
+         END DO
+      ELSE 
+         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)                  &
+                     &                          * e1t(ji,jj) * fse3t(ji,jj,jk)   &
+                     &                          * tmask_i(ji,jj)
+               END DO
+            END DO
+         END DO
+      END IF
       p_fval(:,:) = p_fval(:,:) * 0.5
       IF(lk_mpp) THEN
@@ -224,5 +266,5 @@
       END IF
       !
-   END FUNCTION ptr_vtjk
+   END FUNCTION ptr_tjk
 
 
@@ -237,83 +279,124 @@
          &          zpwatt,     &              ! conversion from W    to PW
          &          zggram                     ! conversion from g    to Pg
-
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &
-         v_atl , v_ipc,                    &
-         vt_atl, vt_pac, vt_ind,           &
-         vs_atl, vs_pac, vs_ind             
-      INTEGER ::  inum       ! temporary logical unit
-      !!----------------------------------------------------------------------
-
-      IF( kt == nit000 .OR. MOD( kt - nit000 + 1, nf_ptr ) == 0 )   THEN
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: vt, vs
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 .OR. MOD( kt, nf_ptr ) == 0 )   THEN
 
          zsverdrup = 1.e-6
          zpwatt    = 1.e-15
          zggram    = 1.e-6
-   
-         ! "zonal" mean temperature and salinity at V-points
-         tn_jk(:,:) = ptr_vtjk( tn(:,:,:) ) * surf_jk_r(:,:)
-         sn_jk(:,:) = ptr_vtjk( sn(:,:,:) ) * surf_jk_r(:,:)
-
+         
+         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:
- 
-         IF( ln_subbas ) THEN              ! Basins computation
-
-            IF( kt == nit000 ) THEN                ! load sub-basin mask
-               CALL iom_open( 'subbasins', inum )
-               CALL iom_get( inum, jpdom_data, 'atlmsk', abasin )      ! Atlantic basin
-               CALL iom_get( inum, jpdom_data, 'pacmsk', pbasin )      ! Pacific basin
-               CALL iom_get( inum, jpdom_data, 'indmsk', ibasin )      ! Indian basin
-               CALL iom_close( inum )
-            ENDIF
-
-            ! basin separation:
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! basin separated velocity
-                  v_atl(ji,jj,:) = vn(ji,jj,:)*abasin(ji,jj)
-                  v_ipc(ji,jj,:) = vn(ji,jj,:)*(pbasin(ji,jj)+ibasin(ji,jj))
-
-                  ! basin separated T times V on T points
-                  vt_ind(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) )*0.5
-                  vt_atl(ji,jj,:) = vt_ind(ji,jj,:) * abasin(ji,jj)
-                  vt_pac(ji,jj,:) = vt_ind(ji,jj,:) * pbasin(ji,jj)
-                  vt_ind(ji,jj,:) = vt_ind(ji,jj,:) * ibasin(ji,jj)
-
-                  ! basin separated S times V on T points
-                  vs_ind(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) )*0.5
-                  vs_atl(ji,jj,:) = vs_ind(ji,jj,:) * abasin(ji,jj)
-                  vs_pac(ji,jj,:) = vs_ind(ji,jj,:) * pbasin(ji,jj)
-                  vs_ind(ji,jj,:) = vs_ind(ji,jj,:) * ibasin(ji,jj)
-               END DO
-            END DO
-
-         ENDIF
-
+         
          ! 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
-         IF( ln_subbas ) THEN
-            v_msf_atl(:,:) = ptr_vjk( v_atl (:,:,:) ) 
-            v_msf_ipc(:,:) = ptr_vjk( v_ipc (:,:,:) ) 
-            ht_atl(:) = SUM( ptr_vjk( vt_atl(:,:,:)), 2 )
-            ht_pac(:) = SUM( ptr_vjk( vt_pac(:,:,:)), 2 )
-            ht_ind(:) = SUM( ptr_vjk( vt_ind(:,:,:)), 2 )
-            st_atl(:) = SUM( ptr_vjk( vs_atl(:,:,:)), 2 )
-            st_pac(:) = SUM( ptr_vjk( vs_pac(:,:,:)), 2 )
-            st_ind(:) = SUM( ptr_vjk( vs_ind(:,:,:)), 2 )
-         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
+            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
+            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
+            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(:,:), 2 )   ! SUM over jk
-         pst_ove(:) = SUM( v_msf_glo(:,:) * sn_jk(:,:), 2 )   ! SUM over jk
-#if defined key_diaeiv
-         pht_eiv(:) = SUM( v_msf_eiv(:,:) * tn_jk(:,:), 2 )   ! SUM over jk
-         pst_eiv(:) = SUM( v_msf_eiv(:,:) * sn_jk(:,:), 2 )   ! SUM over jk
-#endif
-      
+         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
@@ -332,30 +415,9 @@
             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
-         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 ) THEN
-            DO jk = 2,jpk 
-               v_msf_atl(:,jk) = v_msf_atl(:,jk-1) + v_msf_atl(:,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_ipc(:,:) = v_msf_ipc(:,:) * zsverdrup
+            st_ipc(:) = st_ipc(:) * zggram
          ENDIF
 
@@ -377,7 +439,7 @@
       !! ** Purpose :   Initialization, namelist read
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_1         ! temporary workspace
-
-      NAMELIST/namptr/ ln_diaptr, ln_subbas, nf_ptr
+      NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, nf_ptr, nf_ptr_wri
+      INTEGER ::  inum       ! temporary logical unit
+      INTEGER, DIMENSION (1) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
       !!----------------------------------------------------------------------
 
@@ -395,11 +457,62 @@
          WRITE(numout,*) '             Atla/Paci/Ind basins computation         ln_subbas = ', ln_subbas
          WRITE(numout,*) '             Frequency of computation                    nf_ptr = ', nf_ptr
+         WRITE(numout,*) '             Frequency of outputs                    nf_ptr_wri = ', nf_ptr_wri
       ENDIF
 
+      IF( ln_subbas ) THEN                ! load sub-basin mask
+         CALL iom_open( 'subbasins', inum )
+         CALL iom_get( inum, jpdom_data, 'atlmsk', abasin )      ! Atlantic basin
+         CALL iom_get( inum, jpdom_data, 'pacmsk', pbasin )      ! Pacific basin
+         CALL iom_get( inum, jpdom_data, 'indmsk', ibasin )      ! Indian basin
+         CALL iom_close( inum )
+         dbasin(:,:) = MAX ( pbasin(:,:), ibasin(:,:) )
+         sbasin(:,:) = tmask (:,:,1)
+         WHERE ( gphit (:,:) < -30.e0) sbasin(:,:) = 0.e0
+      ENDIF
+      
       ! inverse of the ocean "zonal" v-point section
-      z_1(:,:,:) = 1.e0
-      surf_jk_r(:,:) = ptr_vtjk( z_1(:,:,:) )
-      WHERE( surf_jk_r(:,:) /= 0.e0 )   surf_jk_r(:,:) = 1.e0 / surf_jk_r(:,:)
-
+      surf_jk_glo(:,:) = ptr_tjk( tmask(:,:,:) )
+      surf_jk_r_glo(:,:) = 0.e0
+      WHERE( surf_jk_glo(:,:) /= 0.e0 )   surf_jk_r_glo(:,:) = 1.e0 / surf_jk_glo(:,:)
+      
+      IF (ln_subbas) THEN
+         surf_jk_atl(:,:) = ptr_tjk( tmask (:,:,:), abasin(:,:) )
+         surf_jk_r_atl(:,:) = 0.e0
+         WHERE( surf_jk_atl(:,:) /= 0.e0 )   surf_jk_r_atl(:,:) = 1.e0 / surf_jk_atl(:,:)
+         !
+         surf_jk_pac(:,:) = ptr_tjk( tmask (:,:,:), pbasin(:,:) )
+         surf_jk_r_pac(:,:) = 0.e0
+         WHERE( surf_jk_pac(:,:) /= 0.e0 )   surf_jk_r_pac(:,:) = 1.e0 / surf_jk_pac(:,:)
+         ! 
+         surf_jk_ind(:,:) = ptr_tjk( tmask (:,:,:), ibasin(:,:) )
+         surf_jk_r_ind(:,:) = 0.e0
+         WHERE( surf_jk_ind(:,:) /= 0.e0 )   surf_jk_r_ind(:,:) = 1.e0 / surf_jk_ind(:,:)
+         !
+         surf_jk_ipc(:,:) = ptr_tjk( tmask (:,:,:), dbasin(:,:) )
+         surf_jk_r_ipc(:,:) = 0.e0
+         WHERE( surf_jk_ipc(:,:) /= 0.e0 )   surf_jk_r_ipc(:,:) = 1.e0 / surf_jk_ipc(:,:)
+      END IF
+
+      
+      !!----------------------------------------------------------------------
+
+      iglo (1) = jpjglo
+      iloc (1) = nlcj
+      iabsf(1) = njmppt(narea)
+      iabsl(:) = iabsf(:) + iloc(:) - 1
+      ihals(1) = nldj - 1
+      ihale(1) = nlcj - nlej
+      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)
+      
    END SUBROUTINE dia_ptr_init
 
@@ -415,14 +528,21 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!
-      INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw, ndex(1)
-
-      CHARACTER (len=40)       ::   clhstnam, clop                   ! temporary names
-      INTEGER                  ::   iline, it, ji, itmod             !
+      INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw
+      INTEGER, SAVE :: ndim  , ndim_atl     , ndim_pac     , ndim_ind     , ndim_ipc
+      INTEGER, SAVE ::         ndim_atl_30  , ndim_pac_30  , ndim_ind_30  , ndim_ipc_30
+      INTEGER, SAVE :: ndim_h, ndim_h_atl_30, ndim_h_pac_30, ndim_h_ind_30, ndim_h_ipc_30
+      INTEGER, SAVE, DIMENSION (jpj*jpk) :: ndex  , ndex_atl     , ndex_pac     , ndex_ind     , ndex_ipc
+      INTEGER, SAVE, DIMENSION (jpj*jpk) ::         ndex_atl_30  , ndex_pac_30  , ndex_ind_30  , ndex_ipc_30
+      INTEGER, SAVE, DIMENSION (jpj)     :: ndex_h, ndex_h_atl_30, ndex_h_pac_30, ndex_h_ind_30, ndex_h_ipc_30
+
+      CHARACTER (len=40)       ::   clhstnam, clop, clop_once, cl_comment   ! temporary names
+      INTEGER                  ::   iline, it, itmod, ji, jj, jk      !
       REAL(wp)                 ::   zsto, zout, zdt, zjulian   ! temporary scalars
       REAL(wp), DIMENSION(jpj) ::   zphi, zfoo
-      !!----------------------------------------------------------------------
-      
+      REAL(wp), DIMENSION(jpj,jpk) :: z_1
+      !!----------------------------------------------------------------------
+
       ! define time axis
-      it = kt
+      it    = kt
       itmod = kt - nit000 + 1
 
@@ -430,5 +550,5 @@
       ! --------------
       IF( kt == nit000 ) THEN
-         
+
          zdt = rdt
          IF( nacc == 1 ) zdt = rdtmin
@@ -444,21 +564,44 @@
             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(:) = 0.e0
+            zphi(:) = -9999.9999e0
             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 ORCA05
-               IF( jp_cfg == 05  ) zphi(jpj) = zphi(jpjm1) + (zphi(jpjm1)-zphi(jpj-2))/2.
-               IF( jp_cfg == 05  ) zphi(jpj) = MIN( zphi(jpj), 90.)
-
+               ! 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 ) = MIN( zphi(jj), 90.)
+                  END DO
+               END IF
+               IF( jp_cfg == 2 .OR. jp_cfg == 4 ) THEN
+                  DO jj = mj0(jpjdta-1), mj1(jpjdta-1) 
+                     zphi( jj ) = 88.5e0
+                  END DO
+                  DO jj = mj0(jpjdta  ), mj1(jpjdta  ) 
+                     zphi( jj ) = 89.5e0
+                  END DO
+               END IF
             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 )        
 
             !                                        ! =======================
-         ELSE                                        !   OTHER configurations
+         ELSE                                        !   OTHER configurations  zjulian = zjulian - adatrj 
+                                                     !   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 
@@ -466,8 +609,13 @@
          ! Define frequency of output and means
          zsto = nf_ptr * zdt
-         IF( ln_mskland )   THEN   ;   clop = "ave(only(x))"   ! put 1.e+20 on land (very expensive!!)
-         ELSE                      ;   clop = "ave(x)"         ! no use of the mask value (require less cpu time)
-         ENDIF
-         zout = 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
 
@@ -477,10 +625,10 @@
          zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
 
-         CALL dia_nam( clhstnam, nf_ptr, 'diaptr' )
+         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 )
+            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",   &
@@ -488,20 +636,95 @@
          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
-         
-         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 )
-
-         !  Meridional Stream-Function (eulerian and bolus)
-         
-         CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: Global","Sv" ,   &
+
+         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 ) THEN
-            CALL histdef( numptr, "zomsfatl", "Meridional Stream-Function: Atlantic","Sv" ,   &
+         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, "zomsfipc", "Meridional Stream-Function: Indo-Pacific","Sv" ,&
+            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
@@ -516,9 +739,11 @@
             "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
          IF( ln_subbas ) THEN
-            CALL histdef( numptr, "sohtatl", "Heat Transport Atlantic"      ,  &
+            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"      ,   &
+            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 Indic"      ,     &
+            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
@@ -544,12 +769,14 @@
 #endif
          IF( ln_subbas ) THEN
-            CALL histdef( numptr, "sostatl", "Salt Transport Atlantic"      ,    &
+            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"      ,     &
+            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 Indic"      ,       &
+            CALL histdef( numptr, "sostind", "Salt Transport Indian"//TRIM(cl_comment)      ,    &
                "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
-         ENDIF
-         
+            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 )
 
@@ -560,41 +787,68 @@
          IF(lwp) THEN
             WRITE(numout,*)
-            WRITE(numout,*) 'dia_ptr : write Poleward Transports at time-step : ', kt
+            WRITE(numout,*) 'dia_ptr : compute Poleward Transports at time-step : ', kt
             WRITE(numout,*) '~~~~~~~~'
             WRITE(numout,*)
          ENDIF
 
-         ndex(1) = 0
-         CALL histwrite( numptr, "zotemglo", it, tn_jk    , jpj*jpk, ndex )
-         CALL histwrite( numptr, "zosalglo", it, sn_jk    , jpj*jpk, ndex )
+         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 )
+            END IF
+         ENDIF
+
          ! overturning outputs:
-         CALL histwrite( numptr, "zomsfglo", it, v_msf_glo , jpj*jpk, ndex )
-         IF( ln_subbas ) THEN
-            CALL histwrite( numptr, "zomsfatl", it, v_msf_atl , jpj*jpk, ndex )
-            CALL histwrite( numptr, "zomsfipc", it, v_msf_ipc , jpj*jpk, ndex )
-         ENDIF
-         ! heat transport outputs:
-         IF( ln_subbas ) THEN
-            CALL histwrite( numptr, "sohtatl", it, ht_atl  , jpj, ndex )
-            CALL histwrite( numptr, "sohtpac", it, ht_pac  , jpj, ndex )
-            CALL histwrite( numptr, "sohtind", it, ht_ind  , jpj, ndex )
-            CALL histwrite( numptr, "sostatl", it, st_atl  , jpj, ndex )
-            CALL histwrite( numptr, "sostpac", it, st_pac  , jpj, ndex )
-            CALL histwrite( numptr, "sostind", it, st_ind  , jpj, ndex )
-         ENDIF
-
-         CALL histwrite( numptr, "sophtadv", it, pht_adv  , jpj, ndex )
-         CALL histwrite( numptr, "sophtldf", it, pht_ldf  , jpj, ndex )
-         CALL histwrite( numptr, "sophtove", it, pht_ove  , jpj, ndex )
-         CALL histwrite( numptr, "sopstadv", it, pst_adv  , jpj, ndex )
-         CALL histwrite( numptr, "sopstldf", it, pst_ldf  , jpj, ndex )
-         CALL histwrite( numptr, "sopstove", it, pst_ove  , jpj, ndex )
-#if defined key_diaeiv
-         CALL histwrite( numptr, "zomsfeiv", it, v_msf_eiv, jpj*jpk, ndex )
-         CALL histwrite( numptr, "sophteiv", it, pht_eiv  , jpj    , ndex )
-         CALL histwrite( numptr, "sopsteiv", it, pst_eiv  , jpj    , ndex )
-#endif
- 
+         CALL histwrite( numptr, "zomsfglo", it, 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 )
+         ENDIF
+#if defined key_diaeiv
+         CALL histwrite( numptr, "zomsfeiv", it, v_msf_eiv, ndim  , ndex   )
+#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
