Index: /trunk/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 405)
+++ /trunk/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 406)
@@ -37,22 +37,32 @@
 
    !! * Share Module variables
-   LOGICAL, PUBLIC ::   & !!! ** init namelist (namptr) **
-      ln_diaptr = .FALSE.   !: Poleward transport flag (T) or not (F)
-   INTEGER, PUBLIC ::   &  !!: ** ptr namelist (namptr) **
-      nf_ptr = 15           !: frequency of ptr computation
-   REAL(wp), PUBLIC, DIMENSION(jpj) ::   &   ! poleward transport
-      pht_adv, pst_adv,  &  !: heat and salt: advection
-      pht_ove, pst_ove,  &  !: heat and salt: overturning
-      pht_ldf, pst_ldf,  &  !: heat and salt: lateral diffusion
-      pht_eiv, pst_eiv      !: heat and salt: bolus advection
+   LOGICAL, PUBLIC ::       & !!! ** init namelist (namptr) **
+      ln_diaptr = .FALSE.,  &  !: Poleward transport flag (T) or not (F)
+      ln_subbas = .FALSE.      !: Atlantic/Pacific/Indian basins calculation
+   INTEGER, PUBLIC ::       & !!: ** ptr namelist (namptr) **
+      nf_ptr = 15              !: frequency of ptr computation
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   &   !!: poleward transport
+      pht_adv, pst_adv,     &  !: heat and salt: advection
+      pht_ove, pst_ove,     &  !: heat and salt: overturning
+      pht_ldf, pst_ldf,     &  !: heat and salt: lateral diffusion
+#if defined key_diaeiv
+      pht_eiv, pst_eiv,     &  !: heat and salt: bolus advection
+#endif
+      ht_atl,ht_ind,ht_pac, &  !: heat
+      st_atl,st_ind,st_pac     !: salt
+   REAL(wp),DIMENSION(jpi,jpj) ::   &
+      abasin,pbasin,ibasin     !: return function value
+     
 
    !! Module variables
    REAL(wp), DIMENSION(jpj,jpk) ::   &  
-      tn_jk  , sn_jk  ,  &  ! "zonal" mean temperature and salinity
-      v_msf           ,  &  ! "meridional" Stream-Function
+      tn_jk  , sn_jk  ,  &  !: "zonal" mean temperature and salinity
+      v_msf_atl       ,  &  !: "meridional" Stream-Function
+      v_msf_glo       ,  &  !: "meridional" Stream-Function
+      v_msf_ipc       ,  &  !: "meridional" Stream-Function
 #if defined key_diaeiv
-      v_msf_eiv       ,  &  ! bolus "meridional" Stream-Function
-#endif
-      surf_jk_r             ! inverse of the ocean "zonal" section surface
+      v_msf_eiv       ,  &  !: bolus "meridional" Stream-Function
+#endif
+      surf_jk_r             !: inverse of the ocean "zonal" section surface
 
    !! * Substitutions
@@ -190,16 +200,19 @@
  
       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) * tmask_i(ji,jj+1) * tmask_i(ji,jj)
-            END DO
+           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( lk_mpp)   THEN
-          ish(1) = jpj*jpk ; ish2(1)=jpj ; ish2(2)=jpk
-          zwork(:)= RESHAPE(p_fval, ish )
-          CALL mpp_sum(zwork, jpj*jpk )
-          p_fval(:,:)= RESHAPE(zwork,ish2)
+
+      IF(lk_mpp)   THEN
+         ish(1) = jpj*jpk ; ish2(1)=jpj ; ish2(2)=jpk
+         zwork(:)= RESHAPE(p_fval, ish )
+         CALL mpp_sum(zwork, jpj*jpk )
+         p_fval(:,:)= RESHAPE(zwork,ish2)
       END IF
 
@@ -220,4 +233,5 @@
       !! History :
       !!   9.0  !  03-09  (G. Madec)  Original code
+      !!   9.0  !  06-01  (A. Biastoch)  Allow sub-basins computation
       !!----------------------------------------------------------------------
       !! * arguments
@@ -246,9 +260,9 @@
       END DO
       p_fval(:,:) = p_fval(:,:) * 0.5
-      IF( lk_mpp)   THEN
-          ish(1) = jpj*jpk ; ish2(1)=jpj ; ish2(2)=jpk
-          zwork(:)= RESHAPE(p_fval, ish )
-          CALL mpp_sum(zwork, jpj*jpk )
-          p_fval(:,:)= RESHAPE(zwork,ish2)
+      IF(lk_mpp)   THEN
+         ish(1) = jpj*jpk ; ish2(1)=jpj ; ish2(2)=jpk
+         zwork(:)= RESHAPE(p_fval, ish )
+         CALL mpp_sum(zwork, jpj*jpk )
+         p_fval(:,:)= RESHAPE(zwork,ish2)
       END IF
 
@@ -260,57 +274,211 @@
       !!                  ***  ROUTINE dia_ptr  ***
       !!----------------------------------------------------------------------
+      !! * Moudules used
+      USE ioipsl
+
       !! * Argument
       INTEGER, INTENT(in) ::   kt   ! ocean time step index
 
       !! * Local variables
-      INTEGER ::   jk               ! dummy loop
+      INTEGER ::   jk,jj,ji               ! dummy loop
       REAL(wp) ::    &
          zsverdrup,  &              ! conversion from m3/s to Sverdrup
          zpwatt,     &              ! conversion from W    to PW
          zggram                     ! conversion from g    to Pg
-      !!----------------------------------------------------------------------
-      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(:,:)
-
-      ! "zonal" mean mass flux at V-points
-      v_msf(:,:) = ptr_vjk( vn(:,:,:) ) 
+
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &
+         v_atl , v_ipc,                    &
+         vt_atl, vt_pac, vt_ind,           &
+         vs_atl, vs_pac, vs_ind,           &
+         zv_eiv
+      CHARACTER (len=32) ::   &
+         clnam = 'subbasins.nc'                
+      INTEGER ::  itime,inum,ipi,ipj,ipk       ! temporary integer
+      INTEGER, DIMENSION (1) ::   istep
+      REAL(wp) ::    zdate0,zsecond,zdt        ! temporary scalars
+      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
+         zlamt, zphit, zdta             ! temporary workspace (NetCDF read)
+      REAL(wp), DIMENSION(jpk) ::   &
+         zdept                          ! temporary workspace (NetCDF read)
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 .OR. MOD( kt, nf_ptr ) == 0 )   THEN
+
+         zsverdrup = 1.e-6
+         zpwatt    = 1.e-15
+         zggram    = 1.e-6
+         ipi       = jpidta
+         ipj       = jpjdta
+         ipk       = 1
+         itime     = 1
+         zsecond   = 0.e0
+         zdate0    = 0.e0
+   
+# if defined key_diaeiv
+         zv_eiv(:,:,:) = v_eiv(:,:,:)
+# else
+         zv_eiv(:,:,:) = 0.e0
+# endif
+
+         ! "zonal" mean temperature and salinity at V-points
+         tn_jk(:,:) = ptr_vtjk( tn(:,:,:) ) * surf_jk_r(:,:)
+         sn_jk(:,:) = ptr_vtjk( sn(:,:,:) ) * surf_jk_r(:,:)
+
+         !--------------------------------------------------------
+         ! overturning calculation:
+ 
+         IF( ln_subbas ) THEN              ! Basins computation
+
+            IF( kt == nit000 ) THEN                ! load basin mask
+               itime = 1
+               ipi   = jpidta
+               ipj   = jpjdta
+               ipk   = 1
+               zdt   = 0.e0
+               istep = 0
+               clnam = 'subbasins.nc'
+
+               CALL flinopen(clnam,1,jpidta,1,jpjdta,.FALSE.,ipi,ipj, &
+                  &          ipk,zlamt,zphit,zdept,itime,istep,zdate0,zdt,inum)
+
+               ! get basins:
+               abasin (:,:) = 0.e0
+               pbasin (:,:) = 0.e0
+               ibasin (:,:) = 0.e0
+
+               ! Atlantic basin
+               CALL flinget(inum,'atlmsk',jpidta,jpjdta,1,itime,1,   &
+                  &         0,1,jpidta,1,jpjdta,zdta(:,:))
+               DO jj = 1, nlcj                                 ! interior values
+                  DO ji = 1, nlci
+                     abasin (ji,jj) = zdta( mig(ji), mjg(jj) )
+                  END DO
+               END DO
+
+               ! Pacific basin
+               CALL flinget(inum,'pacmsk',jpidta,jpjdta,1,itime,1,   &
+                  &         0,1,jpidta,1,jpjdta,zdta(:,:))
+               DO jj = 1, nlcj                                 ! interior values
+                  DO ji = 1, nlci
+                     pbasin (ji,jj) = zdta( mig(ji), mjg(jj) )
+                  END DO
+               END DO
+
+               ! Indian basin
+               CALL flinget(inum,'indmsk',jpidta,jpjdta,1,itime,1,   &
+                  &         0,1,jpidta,1,jpjdta,zdta(:,:))
+               DO jj = 1, nlcj                                 ! interior values
+                  DO ji = 1, nlci
+                     ibasin (ji,jj) = zdta( mig(ji), mjg(jj) )
+                  END DO
+               END DO
+
+               CALL flinclo(inum)
+
+            ENDIF
+
+            ! basin separation:
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! basin separated velocity
+                  v_atl(ji,jj,:) = (vn(ji,jj,:)+zv_eiv(ji,jj,:))*abasin(ji,jj)
+                  v_ipc(ji,jj,:) = (vn(ji,jj,:)+zv_eiv(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   &
+                     &              + (zv_eiv(ji,jj,:) + zv_eiv(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   &
+                     &              + (zv_eiv(ji,jj,:) + zv_eiv(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 
+         v_msf_glo(:,:) = ptr_vjk( vn(:,:,:) ) 
 #if defined key_diaeiv
-      ! "zonal" mean bolus mass flux at V-points
-      v_msf_eiv(:,:) = ptr_vjk( v_eiv(:,:,:) ) 
-      ! Bolus "Meridional" Stream-Function
-      DO jk = jpkm1, 1 , -1
-         v_msf_eiv(:,jk) = v_msf_eiv(:,jk-1) + v_msf_eiv(:,jk)
-      END DO
-      v_msf_eiv(:,:) = v_msf_eiv(:,:) * zsverdrup
-#endif
-
-      ! poleward transport: overturning component
-      pht_ove(:) = SUM( v_msf(:,:) * tn_jk(:,:), 2 )   ! SUM over jk
-      pst_ove(:) = SUM( v_msf(:,:) * sn_jk(:,:), 2 )   ! SUM over jk
-
-      ! 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
-      pht_eiv(:) = pht_eiv(:) * zpwatt
-      pst_adv(:) = pst_adv(:) * zggram
-      pst_ove(:) = pst_ove(:) * zggram
-      pst_ldf(:) = pst_ldf(:) * zggram
-      pst_eiv(:) = pst_eiv(:) * zggram
-
-      ! "Meridional" Stream-Function
-      DO jk = jpkm1, 1, -1
-         v_msf(:,jk) = v_msf(:,jk-1) + v_msf(:,jk)
-      END DO
-      v_msf(:,:) = v_msf(:,:) * zsverdrup
-
-      ! output
-      CALL dia_ptr_wri( kt )
+         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
+
+         ! 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
+      
+         ! 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
+            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
+         ENDIF
+
+         ! outputs
+         CALL dia_ptr_wri( kt )
+
+      ENDIF
+
+      ! Close the file
+      IF( kt == nitend ) CALL histclo( numptr )
 
    END SUBROUTINE dia_ptr
@@ -333,8 +501,7 @@
       !!----------------------------------------------------------------------
       !! * local declarations
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         z_1             ! temporary workspace
-
-      NAMELIST/namptr/ ln_diaptr, nf_ptr
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_1         ! temporary workspace
+
+      NAMELIST/namptr/ ln_diaptr, ln_subbas, nf_ptr
       !!----------------------------------------------------------------------
 
@@ -351,4 +518,5 @@
          WRITE(numout,*) '          Namelist namptr : set ptr parameters'
          WRITE(numout,*) '             Switch for ptr diagnostic (T) or not (F) ln_diaptr = ', ln_diaptr
+         WRITE(numout,*) '             Atla/Paci/Ind basins computation         ln_subbas = ', ln_subbas
          WRITE(numout,*) '             Frequency of computation                    nf_ptr = ', nf_ptr
       ENDIF
@@ -383,6 +551,5 @@
 
       !! * Save variables   
-      INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw   &
-          , ndex(1)
+      INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw, ndex(1)
 
       !! * Local variables
@@ -435,4 +602,8 @@
             DO ji = mi0(iline), mi1(iline) 
                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.)
+
             END DO
             ! provide the correct zphi to all local domains
@@ -458,5 +629,5 @@
          CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian )
 
-         CALL dia_nam( clhstnam, nwrite, 'diaptr' )
+         CALL dia_nam( clhstnam, nf_ptr, 'diaptr' )
          IF(lwp)WRITE( numout,*)" Name of diaptr NETCDF file ",clhstnam
 
@@ -479,6 +650,12 @@
          !  Meridional Stream-Function (eulerian and bolus)
          
-         CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: global","Sv" ,   &
+         CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: Global","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" ,   &
+               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
+            CALL histdef( numptr, "zomsfipc", "Meridional Stream-Function: Indo-Pacific","Sv" ,&
+               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
+         ENDIF
 
          !  Heat transport 
@@ -490,4 +667,13 @@
          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"      ,  &
+               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            CALL histdef( numptr, "sohtpac", "Heat Transport Pacific"      ,   &
+               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            CALL histdef( numptr, "sohtind", "Heat Transport Indic"      ,     &
+               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+         ENDIF
+
 
          !  Salt transport 
@@ -509,4 +695,13 @@
             "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"      ,    &
+               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            CALL histdef( numptr, "sostpac", "Salt Transport Pacific"      ,     &
+               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+            CALL histdef( numptr, "sostind", "Salt Transport Indic"      ,       &
+               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
+         ENDIF
+         
 
          CALL histend( numptr )
@@ -515,17 +710,39 @@
 
       IF( MOD( kt, nf_ptr ) == 0 ) THEN
+
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'dia_ptr : write Poleward Transports at time-step : ', kt
+            WRITE(numout,*) '~~~~~~~~'
+            WRITE(numout,*)
+         ENDIF
 
          ! define time axis
          it= kt - nit000 + 1
-         ndex(1) = 1
+         ndex(1) = 0
          CALL histwrite( numptr, "zotemglo", it, tn_jk    , jpj*jpk, ndex )
          CALL histwrite( numptr, "zosalglo", it, sn_jk    , jpj*jpk, ndex )
-         CALL histwrite( numptr, "zomsfglo", it, v_msf    , jpj*jpk, ndex )
-         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 )
+         ! 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 )
@@ -536,7 +753,4 @@
       ENDIF
 
-      ! Close the file
-      IF( kt == nitend )   CALL histclo( numptr )           ! Netcdf write
-
    END SUBROUTINE dia_ptr_wri
 
