Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 4756)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 5260)
@@ -8,5 +8,5 @@
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain
-   USE insitu_tem, ONLY: insitu_t, theta2t
+   USE diainsitutem, ONLY: insitu_t, theta2t
    USE in_out_manager  ! I/O units
    USE iom             ! I/0 library
Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 4756)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 5260)
@@ -21,4 +21,7 @@
    USE timing         ! preformance summary
    USE wrk_nemo       ! working arrays
+   USE fldread        ! type FLD_N
+   USE phycst         ! physical constant
+   USE in_out_manager  ! I/O manager
 
    IMPLICIT NONE
@@ -83,6 +86,4 @@
       CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 )
 
-      CALL iom_put( 'cellthc', fse3t(:,:,:) )
-
       zarea_ssh(:,:) = area(:,:) * sshn(:,:)
 
@@ -104,5 +105,15 @@
          zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
       END DO
-      IF( .NOT.lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+      IF( .NOT.lk_vvl ) THEN
+         IF ( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+               END DO
+            END DO
+         ELSE
+            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+         END IF
+      END IF
       !                                         
       zarho = SUM( area(:,:) * zbotpres(:,:) ) 
@@ -120,5 +131,15 @@
          zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
       END DO
-      IF( .NOT.lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+      IF( .NOT.lk_vvl ) THEN
+         IF ( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+               END DO
+            END DO
+         ELSE
+            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+         END IF
+      END IF
       !    
       zarho = SUM( area(:,:) * zbotpres(:,:) ) 
@@ -145,6 +166,15 @@
       END DO
       IF( .NOT.lk_vvl ) THEN
-         ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
-         zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
+         IF ( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
+                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
+               END DO
+            END DO
+         ELSE
+            ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
+            zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
+         END IF
       ENDIF
       IF( lk_mpp ) THEN  
@@ -181,4 +211,22 @@
       REAL(wp) ::   zztmp  
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
+      ! reading initial file
+      LOGICAL  ::   ln_tsd_init      !: T & S data flag
+      LOGICAL  ::   ln_tsd_tradmp    !: internal damping toward input data flag
+      CHARACTER(len=100)            ::   cn_dir
+      TYPE(FLD_N)                   ::  sn_tem,sn_sal
+      INTEGER  ::   ios=0
+
+      NAMELIST/namtsd/ ln_tsd_init,ln_tsd_tradmp,cn_dir,sn_tem,sn_sal
+      !
+
+      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :
+      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in reference namelist for dia_ar5', lwp )
+      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run
+      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in configuration namelist for dia_ar5', lwp )
+      IF(lwm) WRITE ( numond, namtsd )
+      !
       !!----------------------------------------------------------------------
       !
@@ -200,8 +248,8 @@
       END DO
       IF( lk_mpp )   CALL mpp_sum( vol0 )
-      
-      CALL iom_open ( 'data_1m_salinity_nomask', inum )
-      CALL iom_get  ( inum, jpdom_data, 'vosaline', zsaldta(:,:,:,1), 1  )
-      CALL iom_get  ( inum, jpdom_data, 'vosaline', zsaldta(:,:,:,2), 12 )
+
+      CALL iom_open ( TRIM( cn_dir )//TRIM(sn_sal%clname), inum )
+      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,1), 1  )
+      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 )
       CALL iom_close( inum )
       sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )        
Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 4756)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 5260)
@@ -42,5 +42,4 @@
 #endif
 #if defined key_lim3
-  USE par_ice
   USE ice
 #endif
@@ -113,4 +112,5 @@
   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d  
 
+   !! $Id$
 CONTAINS
 
@@ -1298,4 +1298,5 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
    PUBLIC 
+   !! $Id$
 CONTAINS
 
Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 4756)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 5260)
@@ -7,7 +7,5 @@
    !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
    !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!----------------------------------------------------------------------
-#if ! defined key_coupled
- 
+   !!---------------------------------------------------------------------- 
    !!----------------------------------------------------------------------
    !!   Only for ORCA2 ORCA1 and ORCA025
@@ -29,6 +27,4 @@
 
    PUBLIC dia_fwb    ! routine called by step.F90
-
-   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .TRUE.    !: fresh water budget flag
 
    REAL(wp)               ::   a_fwf ,          &
@@ -453,15 +449,4 @@
    END SUBROUTINE dia_fwb
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                       Dummy Module
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .FALSE.    !: fresh water budget flag
-CONTAINS
-   SUBROUTINE dia_fwb( kt )        ! Empty routine
-      WRITE(*,*) 'dia_fwb: : You should not have seen this print! error?', kt
-   END SUBROUTINE dia_fwb
-#endif
-
    !!======================================================================
 END MODULE diafwb
Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 4756)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 5260)
@@ -18,4 +18,5 @@
    USE daymod
    USE tide_mod
+   !
    USE in_out_manager  ! I/O units
    USE iom             ! I/0 library
@@ -34,9 +35,9 @@
    INTEGER, PARAMETER :: jpdimsparse  = jpincomax*300*24
 
-   !                            !!!namelist variables
+   !                         !!** namelist variables **
    INTEGER ::   nit000_han    ! First time step used for harmonic analysis
    INTEGER ::   nitend_han    ! Last time step used for harmonic analysis
    INTEGER ::   nstep_han     ! Time step frequency for harmonic analysis
-   INTEGER ::   nb_ana           ! Number of harmonics to analyse
+   INTEGER ::   nb_ana        ! Number of harmonics to analyse
 
    INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name
@@ -59,5 +60,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.5 , NEMO Consortium (2013)
-   !! $Id:$
+   !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -119,5 +120,5 @@
             ENDIF
          END DO
-      ENDDO
+      END DO
       !
       IF(lwp) THEN
@@ -158,5 +159,5 @@
       ! ----------------------------
       ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
-      ana_temp(:,:,:,:) = 0.e0
+      ana_temp(:,:,:,:) = 0._wp
 
    END SUBROUTINE dia_harm_init
@@ -179,32 +180,31 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_harm')
 
-      IF ( kt == nit000 ) CALL dia_harm_init
-
-      IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. &
-           (MOD(kt,nstep_han).EQ.0) ) THEN
-
-        ztime = (kt-nit000+1)*rdt 
+      IF( kt == nit000 ) CALL dia_harm_init
+
+      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
+
+         ztime = (kt-nit000+1) * rdt 
        
-        nhc = 0
-        DO jh = 1,nb_ana
-          DO jc = 1,2
-            nhc = nhc+1
-            ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  &
-                    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
-
-            DO jj = 1,jpj
-              DO ji = 1,jpi
-                ! Elevation
-                ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask(ji,jj,1)        
+         nhc = 0
+         DO jh = 1, nb_ana
+            DO jc = 1, 2
+               nhc = nhc+1
+               ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  &
+                  &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
+
+               DO jj = 1,jpj
+                  DO ji = 1,jpi
+                     ! Elevation
+                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask_i(ji,jj)        
 #if defined key_dynspg_ts
-                ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1)
-                ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1)
-#endif
-              END DO
-            END DO
-
-          END DO
-        END DO
-       
+                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask_i(ji,jj)
+                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask_i(ji,jj)
+#endif
+                  END DO
+               END DO
+               !
+            END DO
+         END DO
+         !       
       END IF
 
@@ -249,6 +249,6 @@
          keq = keq + 1
          kun = 0
-         DO jh = 1,nb_ana
-            DO jc = 1,2
+         DO jh = 1, nb_ana
+            DO jc = 1, 2
                kun = kun + 1
                ksp = ksp + 1
@@ -294,9 +294,9 @@
                X1 = ana_amp(ji,jj,jh,1)
                X2 =-ana_amp(ji,jj,jh,2)
-               out_eta(ji,jj,jh       ) = X1 * tmask(ji,jj,1)
-               out_eta(ji,jj,jh+nb_ana) = X2 * tmask(ji,jj,1)
-            ENDDO
-         ENDDO
-      ENDDO
+               out_eta(ji,jj,jh       ) = X1 * tmask_i(ji,jj)
+               out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj)
+            END DO
+         END DO
+      END DO
 
       ! ubar:
@@ -309,6 +309,6 @@
                   kun = kun + 1
                   ztmp4(kun)=ana_temp(ji,jj,kun,2)
-               ENDDO
-            ENDDO
+               END DO
+            END DO
 
             CALL SUR_DETERMINE(jj+1)
@@ -316,6 +316,6 @@
             ! Fill output array
             DO jh = 1, nb_ana
-               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
-               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
+               ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1)
+               ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2)
             END DO
 
@@ -326,8 +326,8 @@
          DO ji = 1, jpi
             DO jh = 1, nb_ana 
-               X1=ana_amp(ji,jj,jh,1)
+               X1= ana_amp(ji,jj,jh,1)
                X2=-ana_amp(ji,jj,jh,2)
-               out_u(ji,jj,jh) = X1 * umask(ji,jj,1)
-               out_u (ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1)
+               out_u(ji,jj,       jh) = X1 * umask_i(ji,jj)
+               out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)
             ENDDO
          ENDDO
@@ -343,6 +343,6 @@
                   kun = kun + 1
                   ztmp4(kun)=ana_temp(ji,jj,kun,3)
-               ENDDO
-            ENDDO
+               END DO
+            END DO
 
             CALL SUR_DETERMINE(jj+1)
@@ -362,9 +362,9 @@
                X1=ana_amp(ji,jj,jh,1)
                X2=-ana_amp(ji,jj,jh,2)
-               out_v(ji,jj,jh)=X1 * vmask(ji,jj,1)
-               out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1)
-            ENDDO
-         ENDDO
-      ENDDO
+               out_v(ji,jj,       jh)=X1 * vmask_i(ji,jj)
+               out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)
+            END DO
+         END DO
+      END DO
 
       CALL dia_wri_harm ! Write results in files
@@ -437,118 +437,118 @@
 #else
       DO jh = 1, nb_ana
-         CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh       ) )
-         CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,jh+nb_ana) )
-      END DO
-#endif
-
+         CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) )
+         CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) )
+      END DO
+#endif
+      !
    END SUBROUTINE dia_wri_harm
 
 
    SUBROUTINE SUR_DETERMINE(init)
-   !!---------------------------------------------------------------------------------
-   !!                      *** ROUTINE SUR_DETERMINE ***
-   !!    
-   !!    
-   !!       
-   !!---------------------------------------------------------------------------------
-   INTEGER, INTENT(in) ::   init 
-   !
-   INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
-   REAL(wp)                        :: zval1, zval2, zx1
-   REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
-   INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
-   !---------------------------------------------------------------------------------
-   CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
-   CALL wrk_alloc( jpincomax , ipos2 , ipivot        )
+      !!---------------------------------------------------------------------------------
+      !!                      *** ROUTINE SUR_DETERMINE ***
+      !!    
+      !!    
+      !!       
+      !!---------------------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   init 
+      !
+      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
+      REAL(wp)                        :: zval1, zval2, zx1
+      REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
+      INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
+      !---------------------------------------------------------------------------------
+      CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
+      CALL wrk_alloc( jpincomax , ipos2 , ipivot        )
             
-   IF( init == 1 ) THEN
-      IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
-      IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
-      !
-      ztmp3(:,:) = 0._wp
-      !
-      DO jk1_sd = 1, nsparse
-         DO jk2_sd = 1, nsparse
-            nisparse(jk2_sd) = nisparse(jk2_sd)
-            njsparse(jk2_sd) = njsparse(jk2_sd)
-            IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
-               ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  &
-                                                        + valuesparse(jk1_sd)*valuesparse(jk2_sd)
-            ENDIF
-         END DO
-      END DO
-
-      DO jj_sd = 1 ,ninco
-          ipos1(jj_sd) = jj_sd
-          ipos2(jj_sd) = jj_sd
-      ENDDO
-
-      DO ji_sd = 1 , ninco
-
-         !find greatest non-zero pivot:
-         zval1 = ABS(ztmp3(ji_sd,ji_sd))
-
-         ipivot(ji_sd) = ji_sd
-         DO jj_sd = ji_sd, ninco
-            zval2 = ABS(ztmp3(ji_sd,jj_sd))
-            IF( zval2.GE.zval1 )THEN
-               ipivot(ji_sd) = jj_sd
-               zval1         = zval2
-            ENDIF
-         ENDDO
-
-         DO ji1_sd = 1, ninco
-            zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd)
-            zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd))
-            ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd)
-            ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
-         ENDDO
-
-         ipos2(ji_sd)         = ipos1(ipivot(ji_sd))
-         ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
-         ipos1(ji_sd)         = ipos2(ji_sd)
-         ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
-         zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd)
-         DO jj_sd = 1, ninco
-            ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
-         ENDDO
-
+      IF( init == 1 ) THEN
+         IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
+         IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
+         !
+         ztmp3(:,:) = 0._wp
+         !
+         DO jk1_sd = 1, nsparse
+            DO jk2_sd = 1, nsparse
+               nisparse(jk2_sd) = nisparse(jk2_sd)
+               njsparse(jk2_sd) = njsparse(jk2_sd)
+               IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
+                  ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  &
+                     &                                     + valuesparse(jk1_sd)*valuesparse(jk2_sd)
+               ENDIF
+            END DO
+         END DO
+         !
+         DO jj_sd = 1 ,ninco
+            ipos1(jj_sd) = jj_sd
+            ipos2(jj_sd) = jj_sd
+         END DO
+         !
+         DO ji_sd = 1 , ninco
+            !
+            !find greatest non-zero pivot:
+            zval1 = ABS(ztmp3(ji_sd,ji_sd))
+            !
+            ipivot(ji_sd) = ji_sd
+            DO jj_sd = ji_sd, ninco
+               zval2 = ABS(ztmp3(ji_sd,jj_sd))
+               IF( zval2.GE.zval1 )THEN
+                  ipivot(ji_sd) = jj_sd
+                  zval1         = zval2
+               ENDIF
+            END DO
+            !
+            DO ji1_sd = 1, ninco
+               zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd)
+               zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd))
+               ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd)
+               ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
+            END DO
+            !
+            ipos2(ji_sd)         = ipos1(ipivot(ji_sd))
+            ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
+            ipos1(ji_sd)         = ipos2(ji_sd)
+            ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
+            zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd)
+            DO jj_sd = 1, ninco
+               ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
+            END DO
+            !
+            DO ji2_sd = ji_sd+1, ninco
+               zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
+               DO jj_sd=1,ninco
+                  ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
+               END DO
+            END DO
+            !
+         END DO
+         !
+      ENDIF ! End init==1
+
+      DO ji_sd = 1, ninco
+         ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
          DO ji2_sd = ji_sd+1, ninco
-            zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
-            DO jj_sd=1,ninco
-               ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
-            ENDDO
-         ENDDO
-
-      ENDDO
-
-   ENDIF ! End init==1
-
-   DO ji_sd = 1, ninco
-      ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
-      DO ji2_sd = ji_sd+1, ninco
-         ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
-      ENDDO
-   ENDDO
-
-   !system solving: 
-   ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
-   ji_sd = ninco
-   DO ji_sd = ninco-1, 1, -1
-      zx1=0.
-      DO jj_sd = ji_sd+1, ninco
-         zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
-      ENDDO
-      ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
-   ENDDO
-
-   DO jj_sd =1, ninco
-      ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
-   ENDDO
-
-   CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
-   CALL wrk_dealloc( jpincomax , ipos2 , ipivot        )
-
-  END SUBROUTINE SUR_DETERMINE
+            ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
+         END DO
+      END DO
+
+      !system solving: 
+      ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
+      ji_sd = ninco
+      DO ji_sd = ninco-1, 1, -1
+         zx1 = 0._wp
+         DO jj_sd = ji_sd+1, ninco
+            zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
+         END DO
+         ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
+      END DO
+
+      DO jj_sd =1, ninco
+         ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
+      END DO
+
+      CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
+      CALL wrk_dealloc( jpincomax , ipos2 , ipivot        )
+      !
+   END SUBROUTINE SUR_DETERMINE
 
 #else
Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4756)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 5260)
@@ -9,21 +9,27 @@
 
    !!----------------------------------------------------------------------
+   !!   dia_hsb       : Diagnose the conservation of ocean heat and salt contents, and volume
+   !!   dia_hsb_rst   : Read or write DIA file in restart file
+   !!   dia_hsb_init  : Initialization of the conservation diagnostic
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
    USE sbc_oce         ! surface thermohaline fluxes
-   USE in_out_manager  ! I/O manager
+   USE sbcrnf          ! river runoff
+   USE sbcisf          ! ice shelves
    USE domvvl          ! vertical scale factors
    USE traqsr          ! penetrative solar radiation
    USE trabbc          ! bottom boundary condition 
-   USE lib_mpp         ! distributed memory computing library
    USE trabbc          ! bottom boundary condition
    USE bdy_par         ! (for lk_bdy)
+   USE restart         ! ocean restart
+   !
+   USE iom             ! I/O manager
+   USE in_out_manager  ! I/O manager
+   USE lib_fortran     ! glob_sum
+   USE lib_mpp         ! distributed memory computing library
    USE timing          ! preformance summary
-   USE iom             ! I/O manager
-   USE lib_fortran     ! glob_sum
-   USE restart         ! ocean restart
-   USE wrk_nemo         ! work arrays
-   USE sbcrnf         ! river runoffd
+   USE wrk_nemo        ! work arrays
 
    IMPLICIT NONE
@@ -36,15 +42,15 @@
    LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
 
-   REAL(dp)                                ::   surf_tot                !
-   REAL(dp)                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends
-   REAL(dp)                                ::   frc_wn_t      , frc_wn_s ! global forcing trends
-   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   surf      , ssh_ini              !
-   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
-   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini
+   REAL(wp) ::   surf_tot              ! ocean surface
+   REAL(wp) ::   frc_t, frc_s, frc_v   ! global forcing trends
+   REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends
+   !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf          , ssh_ini          !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
-
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -52,5 +58,4 @@
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -67,27 +72,37 @@
       !!---------------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !!
-      INTEGER    ::   jk                          ! dummy loop indice
-      REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
-      REAL(dp)   ::   zdiff_hc1   , zdiff_sc1     ! -   -   -   -   -   -   -   - 
-      REAL(dp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
-      REAL(dp)   ::   zerr_hc1    , zerr_sc1       ! heat and salt content misfit
-      REAL(dp)   ::   zvol_tot                    ! volume
-      REAL(dp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
-      REAL(dp)   ::   z_frc_trd_v                 !    -     -
-      REAL(dp)   ::   z_wn_trd_t , z_wn_trd_s   !    -     -
-      REAL(dp)   ::   z_ssh_hc , z_ssh_sc   !    -     -
+      !
+      INTEGER    ::   ji, jj, jk                  ! dummy loop indice
+      REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
+      REAL(wp)   ::   zdiff_hc1   , zdiff_sc1     !  -         -     -        - 
+      REAL(wp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
+      REAL(wp)   ::   zerr_hc1    , zerr_sc1      ! heat and salt content misfit
+      REAL(wp)   ::   zvol_tot                    ! volume
+      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
+      REAL(wp)   ::   z_frc_trd_v                 !    -     -
+      REAL(wp)   ::   z_wn_trd_t , z_wn_trd_s     !    -     -
+      REAL(wp)   ::   z_ssh_hc , z_ssh_sc         !    -     -
+      REAL(wp), DIMENSION(:,:), POINTER ::   z2d0, z2d1
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('dia_hsb')      
-
+      CALL wrk_alloc( jpi,jpj,   z2d0, z2d1 )
+      !
+      tsn(:,:,:,1) = tsn(:,:,:,1) * tmask(:,:,:) ; tsb(:,:,:,1) = tsb(:,:,:,1) * tmask(:,:,:) ;
+      tsn(:,:,:,2) = tsn(:,:,:,2) * tmask(:,:,:) ; tsb(:,:,:,2) = tsb(:,:,:,2) * tmask(:,:,:) ;
       ! ------------------------- !
       ! 1 - Trends due to forcing !
       ! ------------------------- !
-      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) ) ! volume fluxes
-      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )       ! heat fluxes
-      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )       ! salt fluxes
-      ! Add runoff heat & salt input
+      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) ) * surf(:,:) ) ! volume fluxes
+      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )                               ! heat fluxes
+      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )                               ! salt fluxes
+      ! Add runoff    heat & salt input
       IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) )
       IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) )
+      ! Add ice shelf heat & salt input
+      IF( nn_isf .GE. 1 )  THEN
+          z_frc_trd_t = z_frc_trd_t &
+              &   + glob_sum( ( risf_tsc(:,:,jp_tem) - rdivisf * fwfisf(:,:) * (-1.9) * r1_rau0 ) * surf(:,:) )
+          z_frc_trd_s = z_frc_trd_s + (1.0_wp - rdivisf) * glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) )
+      ENDIF
 
       ! Add penetrative solar radiation
@@ -97,6 +112,17 @@
       !
       IF( .NOT. lk_vvl ) THEN
-         z_wn_trd_t = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) )
-         z_wn_trd_s = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) )
+         IF ( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem)
+                  z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal)
+               ENDDO
+            ENDDO
+         ELSE
+            z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem)
+            z2d1(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal)
+         END IF
+         z_wn_trd_t = - glob_sum( z2d0 ) 
+         z_wn_trd_s = - glob_sum( z2d1 )
       ENDIF
 
@@ -113,7 +139,7 @@
       ! 2 -  Content variations !
       ! ------------------------ !
-      zdiff_v2 = 0.d0
-      zdiff_hc = 0.d0
-      zdiff_sc = 0.d0
+      zdiff_v2 = 0._wp
+      zdiff_hc = 0._wp
+      zdiff_sc = 0._wp
 
       ! volume variation (calculated with ssh)
@@ -122,6 +148,17 @@
       ! heat & salt content variation (associated with ssh)
       IF( .NOT. lk_vvl ) THEN
-         z_ssh_hc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) )
-         z_ssh_sc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) )
+         IF ( ln_isfcav ) THEN
+            DO ji = 1, jpi
+               DO jj = 1, jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) ) 
+                  z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) ) 
+               END DO
+            END DO
+         ELSE
+            z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) 
+            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 
+         END IF
+         z_ssh_hc = glob_sum( z2d0 ) 
+         z_ssh_sc = glob_sum( z2d1 ) 
       ENDIF
 
@@ -153,8 +190,15 @@
       ! 3 - Diagnostics writing !
       ! ----------------------- !
-      zvol_tot   = 0.d0                                                   ! total ocean volume
+      zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors)
       DO jk = 1, jpkm1
          zvol_tot  = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
       END DO
+
+!!gm to be added ?
+!      IF( .NOT. lk_vvl ) THEN            ! fixed volume, add the ssh contribution
+!        zvol_tot = zvol_tot + glob_sum( surf(:,:) * sshn(:,:) )
+!      ENDIF
+!!gm end
+
 
       IF( lk_vvl ) THEN
@@ -183,90 +227,10 @@
       IF( lrst_oce )   CALL dia_hsb_rst( kt, 'WRITE' )
 
+      CALL wrk_dealloc( jpi,jpj,   z2d0, z2d1 )
+
       IF( nn_timing == 1 )   CALL timing_stop('dia_hsb')
-!
+      !
    END SUBROUTINE dia_hsb
 
-
-   SUBROUTINE dia_hsb_init
-      !!---------------------------------------------------------------------------
-      !!                  ***  ROUTINE dia_hsb  ***
-      !!     
-      !! ** Purpose: Initialization for the heat salt volume budgets
-      !!	
-      !! ** Method : Compute initial heat content, salt content and volume
-      !!
-      !! ** Action : - Compute initial heat content, salt content and volume
-      !!             - Initialize forcing trends
-      !!             - Compute coefficients for conversion
-      !!---------------------------------------------------------------------------
-      INTEGER            ::   jk       ! dummy loop indice
-      INTEGER            ::   ierror   ! local integer
-      !!
-      NAMELIST/namhsb/ ln_diahsb
-      !
-      INTEGER  ::   ios
-      !!----------------------------------------------------------------------
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
-         WRITE(numout,*) '~~~~~~~~ '
-      ENDIF
-
-      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
-      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist
-      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namhsb )
-
-      !
-      IF(lwp) THEN                   ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
-         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
-         WRITE(numout,*)
-      ENDIF
-
-      IF( .NOT. ln_diahsb )   RETURN
-!      IF( .NOT. lk_mpp_rep ) &
-!        CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', &
-!             &         ' whereas the global sum to be precise must be done in double precision ',&
-!             &         ' please add key_mpp_rep')
-
-      ! ------------------- !
-      ! 1 - Allocate memory !
-      ! ------------------- !
-      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &
-         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
-      ENDIF
-
-      IF(.NOT. lk_vvl ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
-      ENDIF
-
-      ! ----------------------------------------------- !
-      ! 2 - Time independant variables and file opening !
-      ! ----------------------------------------------- !
-      IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
-      IF(lwp) WRITE(numout,*) '~~~~~~~'
-      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
-      surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area
-
-      IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
-      !
-      ! ---------------------------------- !
-      ! 4 - initial conservation variables !
-      ! ---------------------------------- !
-      CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
-      !
-   END SUBROUTINE dia_hsb_init
 
    SUBROUTINE dia_hsb_rst( kt, cdrw )
@@ -281,6 +245,6 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
      !
-     INTEGER ::   jk   ! 
-     INTEGER ::   id1   ! local integers
+     INTEGER ::   ji, jj, jk   ! dummy loop indices
+     INTEGER ::   id1          ! local integers
      !!----------------------------------------------------------------------
      !
@@ -317,12 +281,21 @@
              sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content
           END DO
-          frc_v = 0.d0                                           ! volume       trend due to forcing
-          frc_t = 0.d0                                           ! heat content   -    -   -    -   
-          frc_s = 0.d0                                           ! salt content   -    -   -    -        
+          frc_v = 0._wp                                           ! volume       trend due to forcing
+          frc_t = 0._wp                                           ! heat content   -    -   -    -   
+          frc_s = 0._wp                                           ! salt content   -    -   -    -        
           IF( .NOT. lk_vvl ) THEN
-             ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
-             ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
-             frc_wn_t = 0.d0                                       ! initial heat content misfit due to free surface
-             frc_wn_s = 0.d0                                       ! initial salt content misfit due to free surface
+             IF ( ln_isfcav ) THEN
+                DO ji=1,jpi
+                   DO jj=1,jpj
+                      ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh
+                      ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh
+                   ENDDO
+                ENDDO
+             ELSE
+                ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
+                ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
+             END IF
+             frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
+             frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface
           ENDIF
        ENDIF
@@ -354,4 +327,86 @@
    END SUBROUTINE dia_hsb_rst
 
+
+   SUBROUTINE dia_hsb_init
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!     
+      !! ** Purpose: Initialization for the heat salt volume budgets
+      !!	
+      !! ** Method : Compute initial heat content, salt content and volume
+      !!
+      !! ** Action : - Compute initial heat content, salt content and volume
+      !!             - Initialize forcing trends
+      !!             - Compute coefficients for conversion
+      !!---------------------------------------------------------------------------
+      INTEGER ::   jk       ! dummy loop indice
+      INTEGER ::   ierror   ! local integer
+      INTEGER ::   ios
+      !
+      NAMELIST/namhsb/ ln_diahsb
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
+
+      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
+      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist
+      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namhsb )
+
+      !
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
+         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
+         WRITE(numout,*)
+      ENDIF
+
+      IF( .NOT. ln_diahsb )   RETURN
+         !      IF( .NOT. lk_mpp_rep ) &
+         !        CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', &
+         !             &         ' whereas the global sum to be precise must be done in double precision ',&
+         !             &         ' please add key_mpp_rep')
+
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &
+         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      IF(.NOT. lk_vvl ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      ! ----------------------------------------------- !
+      ! 2 - Time independant variables and file opening !
+      ! ----------------------------------------------- !
+      IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
+      IF(lwp) WRITE(numout,*) '~~~~~~~'
+      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)      ! masked surface grid cell area
+      surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area
+
+      IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
+      !
+      ! ---------------------------------- !
+      ! 4 - initial conservation variables !
+      ! ---------------------------------- !
+      CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
+      !
+   END SUBROUTINE dia_hsb_init
+
    !!======================================================================
 END MODULE diahsb
Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90	(revision 4756)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90	(revision 5260)
@@ -1,3 +1,3 @@
-MODULE insitu_tem
+MODULE diainsitutem
    !!======================================================================
    !!                       ***  MODULE  diaharm  ***
@@ -122,3 +122,3 @@
    END SUBROUTINE ATG
 
-END MODULE insitu_tem
+END MODULE diainsitutem
Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 4756)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 5260)
@@ -8,4 +8,5 @@
    !!            3.2  ! 2010-03  (O. Marti, S. Flavoni) Add fields
    !!            3.3  ! 2010-10  (G. Madec)  dynamical allocation
+   !!            3.6  ! 2014-12  (C. Ethe) use of IOM
    !!----------------------------------------------------------------------
 
@@ -13,59 +14,37 @@
    !!   dia_ptr      : Poleward Transport Diagnostics module
    !!   dia_ptr_init : Initialization, namelist read
-   !!   dia_ptr_wri  : Output of poleward fluxes
-   !!   ptr_vjk      : "zonal" sum computation of a "meridional" flux array
-   !!   ptr_tjk      : "zonal" mean computation of a tracer field
-   !!   ptr_vj       : "zonal" and vertical sum computation of a "meridional" flux array
-   !!                   (Generic interface to ptr_vj_3d, ptr_vj_2d)
+   !!   ptr_sjk      : "zonal" mean computation of a field - tracer or flux array
+   !!   ptr_sj       : "zonal" and vertical sum computation of a "meridional" flux array
+   !!                   (Generic interface to ptr_sj_3d, ptr_sj_2d)
    !!----------------------------------------------------------------------
    USE oce              ! ocean dynamics and active tracers
    USE dom_oce          ! ocean space and time domain
    USE phycst           ! physical constants
-   USE ldftra_oce       ! ocean active tracers: lateral physics
-   USE dianam           !
+   !
    USE iom              ! IOM library
-   USE ioipsl           ! IO-IPSL library
    USE in_out_manager   ! I/O manager
    USE lib_mpp          ! MPP library
-   USE lbclnk           ! lateral boundary condition - processor exchanges
    USE timing           ! preformance summary
-   USE wrk_nemo         ! working arrays
 
    IMPLICIT NONE
    PRIVATE
 
-   INTERFACE ptr_vj
-      MODULE PROCEDURE ptr_vj_3d, ptr_vj_2d
+   INTERFACE ptr_sj
+      MODULE PROCEDURE ptr_sj_3d, ptr_sj_2d
    END INTERFACE
 
-   PUBLIC   dia_ptr_init   ! call in opa module
+   PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines
+   PUBLIC   ptr_sjk        ! 
+   PUBLIC   dia_ptr_init   ! call in step module
    PUBLIC   dia_ptr        ! call in step module
-   PUBLIC   ptr_vj         ! call by tra_ldf & tra_adv routines
-   PUBLIC   ptr_vjk        ! call by tra_ldf & tra_adv routines
 
    !                                  !!** namelist  namptr  **
-   LOGICAL , PUBLIC ::   ln_diaptr     !: Poleward transport flag (T) or not (F)
-   LOGICAL , PUBLIC ::   ln_subbas     !: Atlantic/Pacific/Indian basins calculation
-   LOGICAL , PUBLIC ::   ln_diaznl     !: Add zonal means and meridional stream functions
-   LOGICAL , PUBLIC ::   ln_ptrcomp    !: Add decomposition : overturning (and gyre, soon ...)
-   INTEGER , PUBLIC ::   nn_fptr       !: frequency of ptr computation  [time step]
-   INTEGER , PUBLIC ::   nn_fwri       !: frequency of ptr outputs      [time step]
-
-   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf, htr_ove   !: Heat TRansports (adv, diff, overturn.)
-   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf, str_ove   !: Salt TRansports (adv, diff, overturn.)
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf   !: Heat TRansports (adv, diff, overturn.)
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf   !: Salt TRansports (adv, diff, overturn.)
    
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   btmsk                  ! T-point basin interior masks
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   btm30                  ! mask out Southern Ocean (=0 south of 30°S)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htr  , str             ! adv heat and salt transports (approx)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_jk, sn_jk , v_msf   ! i-mean T and S, j-Stream-Function
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sjk  , r1_sjk          ! i-mean i-k-surface and its inverse        
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htr_eiv, str_eiv       ! bolus adv heat ans salt transports ('key_diaeiv')
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_msf_eiv              ! bolus j-streamfuction              ('key_diaeiv')
-
-
-   INTEGER ::   niter       !
-   INTEGER ::   nidom_ptr   !
-   INTEGER ::   numptr      ! logical unit for Poleward TRansports
-   INTEGER ::   nptr        ! = 1 (ln_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (ln_subbas=T) 
+
+   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F)
+   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation
+   INTEGER         ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T) 
 
    REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
@@ -73,12 +52,11 @@
    REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg
 
-   REAL(wp), TARGET, DIMENSION(:),   ALLOCATABLE, SAVE :: p_fval1d
-   REAL(wp), TARGET, DIMENSION(:,:), ALLOCATABLE, SAVE :: p_fval2d
-
-   !! Integer, 1D workspace arrays. Not common enough to be implemented in 
-   !! wrk_nemo module.
-   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex  , ndex_atl     , ndex_pac     , ndex_ind     , ndex_ipc
-   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::         ndex_atl_30  , ndex_pac_30  , ndex_ind_30  , ndex_ipc_30
-   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex_h, ndex_h_atl_30, ndex_h_pac_30, ndex_h_ind_30, ndex_h_ipc_30
+   CHARACTER(len=3), ALLOCATABLE, SAVE, DIMENSION(:)     :: clsubb
+   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
+   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   :: btm30   ! mask out Southern Ocean (=0 south of 30°S)
+
+   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)     :: p_fval1d
+   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: p_fval2d
+
 
    !! * Substitutions
@@ -92,4 +70,198 @@
 CONTAINS
 
+   SUBROUTINE dia_ptr( pvtr )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_ptr  ***
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
+      !
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::   zv, zsfc               ! local scalar
+      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace
+      CHARACTER( len = 10 )  :: cl1
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('dia_ptr')
+
+      !
+      IF( PRESENT( pvtr ) ) THEN
+         IF( iom_use("zomsfglo") ) THEN    ! effective MSF
+            z3d(1,:,:) = ptr_sjk( pvtr(:,:,:) )  ! zonal cumulative effective transport
+            DO jk = 2, jpkm1 
+              z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)   ! effective j-Stream-Function (MSF)
+            END DO
+            DO ji = 1, jpi
+               z3d(ji,:,:) = z3d(1,:,:)
+            ENDDO
+            cl1 = TRIM('zomsf'//clsubb(1) )
+            CALL iom_put( cl1, z3d * rc_sv )
+            DO jn = 2, nptr                                    ! by sub-basins
+               z3d(1,:,:) =  ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) 
+               DO jk = 2, jpkm1 
+                  z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)    ! effective j-Stream-Function (MSF)
+               END DO
+               DO ji = 1, jpi
+                  z3d(ji,:,:) = z3d(1,:,:)
+               ENDDO
+               cl1 = TRIM('zomsf'//clsubb(jn) )
+               CALL iom_put( cl1, z3d * rc_sv )
+            END DO
+         ENDIF
+         !
+      ELSE
+         !
+         IF( iom_use("zotemglo") ) THEN    ! i-mean i-k-surface 
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zsfc = e1t(ji,jj) * fse3t(ji,jj,jk)
+                     zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
+                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc
+                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc
+                  ENDDO
+               ENDDO
+            ENDDO
+            DO jn = 1, nptr
+               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
+               cl1 = TRIM('zosrf'//clsubb(jn) )
+               CALL iom_put( cl1, zmask )
+               !
+               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) &
+                  &            / MAX( zmask(1,:,:), 10.e-15 )
+               DO ji = 1, jpi
+                  z3d(ji,:,:) = z3d(1,:,:)
+               ENDDO
+               cl1 = TRIM('zotem'//clsubb(jn) )
+               CALL iom_put( cl1, z3d )
+               !
+               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) &
+                  &            / MAX( zmask(1,:,:), 10.e-15 )
+               DO ji = 1, jpi
+                  z3d(ji,:,:) = z3d(1,:,:)
+               ENDDO
+               cl1 = TRIM('zosal'//clsubb(jn) )
+               CALL iom_put( cl1, z3d )
+            END DO
+         ENDIF
+         !
+         !                                ! Advective and diffusive heat and salt transport
+         IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN   
+            z2d(1,:) = htr_adv(:) * rc_pwatt        !  (conversion in PW)
+            DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sophtadv'                 
+            CALL iom_put( TRIM(cl1), z2d )
+            z2d(1,:) = str_adv(:) * rc_ggram        ! (conversion in Gg)
+            DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sopstadv'
+            CALL iom_put( TRIM(cl1), z2d )
+         ENDIF
+         !
+         IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN   
+            z2d(1,:) = htr_ldf(:) * rc_pwatt        !  (conversion in PW) 
+            DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sophtldf'
+            CALL iom_put( TRIM(cl1), z2d )
+            z2d(1,:) = str_ldf(:) * rc_ggram        !  (conversion in Gg)
+            DO ji = 1, jpi
+               z2d(ji,:) = z2d(1,:)
+            ENDDO
+            cl1 = 'sopstldf'
+            CALL iom_put( TRIM(cl1), z2d )
+         ENDIF
+         !
+      ENDIF
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr')
+      !
+   END SUBROUTINE dia_ptr
+
+
+   SUBROUTINE dia_ptr_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_ptr_init  ***
+      !!                   
+      !! ** Purpose :   Initialization, namelist read
+      !!----------------------------------------------------------------------
+      INTEGER ::  jn           ! local integers
+      INTEGER ::  inum, ierr   ! local integers
+      INTEGER ::  ios          ! Local integer output status for namelist read
+      !!
+      NAMELIST/namptr/ ln_diaptr, ln_subbas
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
+      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport
+      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namptr )
+
+      IF(lwp) THEN                     ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namptr : set ptr parameters'
+         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr
+         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas
+      ENDIF
+
+      IF( ln_diaptr ) THEN  
+         !
+         IF( ln_subbas ) THEN 
+            nptr = 5            ! Global, Atlantic, Pacific, Indian, Indo-Pacific
+            ALLOCATE( clsubb(nptr) )
+            clsubb(1) = 'glo' ;  clsubb(2) = 'atl'  ;  clsubb(3) = 'pac'  ;  clsubb(4) = 'ind'  ;  clsubb(5) = 'ipc'
+         ELSE               
+            nptr = 1       ! Global only
+            ALLOCATE( clsubb(nptr) )
+            clsubb(1) = 'glo' 
+         ENDIF
+
+         !                                      ! allocate dia_ptr arrays
+         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
+
+         rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt
+
+         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
+
+         IF( ln_subbas ) THEN                ! load sub-basin mask
+            CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  )
+            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
+            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
+            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
+            CALL iom_close( inum )
+            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
+            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
+            ELSE WHERE                     ;   btm30(:,:) = ssmask(:,:)
+            END WHERE
+         ENDIF
+   
+         btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
+      
+         DO jn = 1, nptr
+            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
+         END DO
+
+         ! Initialise arrays to zero because diatpr is called before they are first calculated
+         ! Note that this means diagnostics will not be exactly correct when model run is restarted.
+         htr_adv(:) = 0._wp  ;  str_adv(:) =  0._wp  
+         htr_ldf(:) = 0._wp  ;  str_ldf(:) =  0._wp 
+         !
+      ENDIF 
+      ! 
+   END SUBROUTINE dia_ptr_init
+
+
    FUNCTION dia_ptr_alloc()
       !!----------------------------------------------------------------------
@@ -97,5 +269,5 @@
       !!----------------------------------------------------------------------
       INTEGER               ::   dia_ptr_alloc   ! return value
-      INTEGER, DIMENSION(6) ::   ierr
+      INTEGER, DIMENSION(3) ::   ierr
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -103,25 +275,10 @@
       ALLOCATE( btmsk(jpi,jpj,nptr) ,           &
          &      htr_adv(jpj) , str_adv(jpj) ,   &
-         &      htr_ldf(jpj) , str_ldf(jpj) ,   &
-         &      htr_ove(jpj) , str_ove(jpj),    &
-         &      htr(jpj,nptr) , str(jpj,nptr) , &
-         &      tn_jk(jpj,jpk,nptr) , sn_jk (jpj,jpk,nptr) , v_msf(jpj,jpk,nptr) , &
-         &      sjk  (jpj,jpk,nptr) , r1_sjk(jpj,jpk,nptr) , STAT=ierr(1)  )
-         !
-#if defined key_diaeiv
-      ALLOCATE( htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) , &
-         &      v_msf_eiv(jpj,jpk,nptr) , STAT=ierr(2) )
-#endif
-      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(3))
-      !
-      ALLOCATE(ndex(jpj*jpk),        ndex_atl(jpj*jpk), ndex_pac(jpj*jpk), &
-         &     ndex_ind(jpj*jpk),    ndex_ipc(jpj*jpk),                    &
-         &     ndex_atl_30(jpj*jpk), ndex_pac_30(jpj*jpk), Stat=ierr(4))
-
-      ALLOCATE(ndex_ind_30(jpj*jpk), ndex_ipc_30(jpj*jpk),                   &
-         &     ndex_h(jpj),          ndex_h_atl_30(jpj), ndex_h_pac_30(jpj), &
-         &     ndex_h_ind_30(jpj),   ndex_h_ipc_30(jpj), Stat=ierr(5) )
-         !
-     ALLOCATE( btm30(jpi,jpj) , STAT=ierr(6)  )
+         &      htr_ldf(jpj) , str_ldf(jpj) , STAT=ierr(1)  )
+         !
+      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2))
+      !
+      ALLOCATE( btm30(jpi,jpj), STAT=ierr(3)  )
+
          !
       dia_ptr_alloc = MAXVAL( ierr )
@@ -131,7 +288,7 @@
 
 
-   FUNCTION ptr_vj_3d( pva )   RESULT ( p_fval )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE ptr_vj_3d  ***
+   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE ptr_sj_3d  ***
       !!
       !! ** Purpose :   i-k sum computation of a j-flux array
@@ -142,6 +299,7 @@
       !! ** 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
-      !!
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pva   ! mask flux array at V-point
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
+      !
       INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
       INTEGER                  ::   ijpj         ! ???
@@ -153,21 +311,31 @@
       ijpj = jpj
       p_fval(:) = 0._wp
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! Vector opt.
-               p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) 
-            END DO
-         END DO
-      END DO
+      IF( PRESENT( pmsk ) ) THEN 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! Vector opt.
+                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) * pmsk(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) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) 
+               END DO
+            END DO
+         END DO
+      ENDIF
 #if defined key_mpp_mpi
       IF(lk_mpp)   CALL mpp_sum( p_fval, ijpj, ncomm_znl)
 #endif
       !
-   END FUNCTION ptr_vj_3d
-
-
-   FUNCTION ptr_vj_2d( pva )   RESULT ( p_fval )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE ptr_vj_2d  ***
+   END FUNCTION ptr_sj_3d
+
+
+   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE ptr_sj_2d  ***
       !!
       !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array
@@ -178,7 +346,7 @@
       !! ** Action  : - p_fval: i-k-mean poleward flux of pva
       !!----------------------------------------------------------------------
-      IMPLICIT none
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point
-      !!
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pva   ! mask flux array at V-point
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
+      !
       INTEGER                  ::   ji,jj       ! dummy loop arguments
       INTEGER                  ::   ijpj        ! ???
@@ -190,24 +358,31 @@
       ijpj = jpj
       p_fval(:) = 0._wp
-      DO jj = 2, jpjm1
-         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)
-         END DO
-      END DO
+      IF( PRESENT( pmsk ) ) THEN 
+         DO jj = 2, jpjm1
+            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) * pmsk(ji,jj)
+            END DO
+         END DO
+      ELSE
+         DO jj = 2, jpjm1
+            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)
+            END DO
+         END DO
+      ENDIF
 #if defined key_mpp_mpi
       CALL mpp_sum( p_fval, ijpj, ncomm_znl )
 #endif
       ! 
-   END FUNCTION ptr_vj_2d
-
-
-   FUNCTION ptr_vjk( pva, pmsk )   RESULT ( p_fval )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE ptr_vjk  ***
-      !!
-      !! ** Purpose :   i-sum computation of a j-velocity array
+   END FUNCTION ptr_sj_2d
+
+
+   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE ptr_sjk  ***
+      !!
+      !! ** Purpose :   i-sum computation of an array
       !!
       !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
-      !!              pva is supposed to be a masked flux (i.e. * vmask)
       !!
       !! ** Action  : - p_fval: i-mean poleward flux of pva
@@ -215,5 +390,5 @@
       !!
       IMPLICIT none
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pva    ! mask flux array at V-point
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pta    ! mask flux array at V-point
       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask
       !!
@@ -224,15 +399,8 @@
       INTEGER, DIMENSION(2) ::   ish2
       INTEGER               ::   ijpjjpk
-#endif
-#if defined key_mpp_mpi
-      REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point
+      REAL(wp), DIMENSION(jpj*jpk) ::   zwork    ! mask flux array at V-point
 #endif
       !!--------------------------------------------------------------------
       !
-#if defined key_mpp_mpi
-      ijpjjpk = jpj*jpk
-      CALL wrk_alloc( jpj*jpk, zwork )
-#endif
-
       p_fval => p_fval2d
 
@@ -244,5 +412,5 @@
 !!gm here, use of tmask_i  ==> no need of loop over nldi, nlei....
                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) * pmsk(ji,jj)
+                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj)
                END DO
             END DO
@@ -252,5 +420,5 @@
             DO jj = 2, jpjm1
                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)
+                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
@@ -266,630 +434,6 @@
 #endif
       !
-#if defined key_mpp_mpi
-      CALL wrk_dealloc( jpj*jpk, zwork )
-#endif
-      !
-   END FUNCTION ptr_vjk
-
-
-   FUNCTION ptr_tjk( pta, pmsk )   RESULT ( p_fval )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE ptr_tjk  ***
-      !!
-      !! ** Purpose :   i-sum computation of e1t*e3t * a tracer field
-      !!
-      !! ** Method  : - i-sum of mj(pta) using tmask
-      !!
-      !! ** Action  : - p_fval: i-sum of e1t*e3t*pta
-      !!----------------------------------------------------------------------
-      !!
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! tracer flux array at T-point
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)     ::   pmsk   ! Optional 2D basin mask
-      !!
-      INTEGER                           :: ji, jj, jk   ! dummy loop arguments
-      REAL(wp), POINTER, DIMENSION(:,:) :: p_fval       ! return function value
-#if defined key_mpp_mpi
-      INTEGER, DIMENSION(1) ::   ish
-      INTEGER, DIMENSION(2) ::   ish2
-      INTEGER               ::   ijpjjpk
-#endif
-#if defined key_mpp_mpi
-      REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point
-#endif
-      !!-------------------------------------------------------------------- 
-      !
-#if defined key_mpp_mpi
-      ijpjjpk = jpj*jpk
-      CALL wrk_alloc( jpj*jpk, zwork )
-#endif
-
-      p_fval => p_fval2d
-
-      p_fval(:,:) = 0._wp
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            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) * pmsk(ji,jj)
-            END DO
-         END DO
-      END DO
-#if defined key_mpp_mpi
-      ijpjjpk = jpj*jpk
-      ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk
-      zwork(1:ijpjjpk)= RESHAPE( p_fval, ish )
-      CALL mpp_sum( zwork, ijpjjpk, ncomm_znl )
-      p_fval(:,:)= RESHAPE( zwork, ish2 )
-#endif
-      !
-#if defined key_mpp_mpi
-      CALL wrk_dealloc( jpj*jpk, zwork )
-#endif
-      !    
-   END FUNCTION ptr_tjk
-
-
-   SUBROUTINE dia_ptr( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dia_ptr  ***
-      !!----------------------------------------------------------------------
-      USE oce,     vt  =>   ua   ! use ua as workspace
-      USE oce,     vs  =>   va   ! use va as workspace
-      IMPLICIT none
-      !!
-      INTEGER, INTENT(in) ::   kt   ! ocean time step index
-      !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::   zv               ! local scalar
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('dia_ptr')
-      !
-      IF( kt == nit000 .OR. MOD( kt, nn_fptr ) == 0 )   THEN
-         !
-         IF( MOD( kt, nn_fptr ) == 0 ) THEN 
-            !
-            IF( ln_diaznl ) THEN               ! i-mean temperature and salinity
-               DO jn = 1, nptr
-                  tn_jk(:,:,jn) = ptr_tjk( tsn(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
-                  sn_jk(:,:,jn) = ptr_tjk( tsn(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
-               END DO
-            ENDIF
-            !
-            !                          ! horizontal integral and vertical dz 
-            !                                ! eulerian velocity
-            v_msf(:,:,1) = ptr_vjk( vn(:,:,:) ) 
-            DO jn = 2, nptr
-               v_msf(:,:,jn) = ptr_vjk( vn(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) 
-            END DO
-#if defined key_diaeiv
-            DO jn = 1, nptr                  ! bolus velocity
-               v_msf_eiv(:,:,jn) = ptr_vjk( v_eiv(:,:,:), btmsk(:,:,jn) )   ! here no btm30 for MSFeiv
-            END DO
-            !                                ! add bolus stream-function to the eulerian one
-            v_msf(:,:,:) = v_msf(:,:,:) + v_msf_eiv(:,:,:)
-#endif
-            !
-            !                          ! Transports
-            !                                ! local heat & salt transports at T-points  ( tsn*mj[vn+v_eiv] )
-            vt(:,:,jpk) = 0._wp   ;   vs(:,:,jpk) = 0._wp
-            DO jk= 1, jpkm1
-               DO jj = 2, jpj
-                  DO ji = 1, jpi
-#if defined key_diaeiv 
-                     zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) + v_eiv(ji,jj,jk) + v_eiv(ji,jj-1,jk) ) * 0.5_wp
-#else
-                     zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) ) * 0.5_wp
-#endif 
-                     vt(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_tem)
-                     vs(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_sal)
-                  END DO
-               END DO
-            END DO
-!!gm useless as overlap areas are not used in ptr_vjk
-            CALL lbc_lnk( vs, 'V', -1. )   ;   CALL lbc_lnk( vt, 'V', -1. )
-!!gm
-            !                                ! heat & salt advective transports (approximation)
-            htr(:,1) = SUM( ptr_vjk( vt(:,:,:) ) , 2 ) * rc_pwatt   ! SUM over jk + conversion
-            str(:,1) = SUM( ptr_vjk( vs(:,:,:) ) , 2 ) * rc_ggram
-            DO jn = 2, nptr 
-               htr(:,jn) = SUM( ptr_vjk( vt(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) , 2 ) * rc_pwatt   ! mask Southern Ocean
-               str(:,jn) = SUM( ptr_vjk( vs(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) , 2 ) * rc_ggram   ! mask Southern Ocean
-            END DO
-
-            IF( ln_ptrcomp ) THEN            ! overturning transport
-               htr_ove(:) = SUM( v_msf(:,:,1) * tn_jk(:,:,1), 2 ) * rc_pwatt   ! SUM over jk + conversion
-               str_ove(:) = SUM( v_msf(:,:,1) * sn_jk(:,:,1), 2 ) * rc_ggram
-            END IF
-            !                                ! Advective and diffusive transport
-            htr_adv(:) = htr_adv(:) * rc_pwatt        ! these are computed in tra_adv... and tra_ldf... routines 
-            htr_ldf(:) = htr_ldf(:) * rc_pwatt        ! here just the conversion in PW and Gg
-            str_adv(:) = str_adv(:) * rc_ggram
-            str_ldf(:) = str_ldf(:) * rc_ggram
-
-#if defined key_diaeiv
-            DO jn = 1, nptr                  ! Bolus component
-               htr_eiv(:,jn) = SUM( v_msf_eiv(:,:,jn) * tn_jk(:,:,jn), 2 ) * rc_pwatt   ! SUM over jk
-               str_eiv(:,jn) = SUM( v_msf_eiv(:,:,jn) * sn_jk(:,:,jn), 2 ) * rc_ggram   ! SUM over jk
-            END DO
-#endif
-            !                                ! "Meridional" Stream-Function
-            DO jn = 1, nptr
-               DO jk = 2, jpk 
-                  v_msf    (:,jk,jn) = v_msf    (:,jk-1,jn) + v_msf    (:,jk,jn)       ! Eulerian j-Stream-Function
-#if defined key_diaeiv
-                  v_msf_eiv(:,jk,jn) = v_msf_eiv(:,jk-1,jn) + v_msf_eiv(:,jk,jn)       ! Bolus    j-Stream-Function
-
-#endif
-               END DO
-            END DO
-            v_msf    (:,:,:) = v_msf    (:,:,:) * rc_sv       ! converte in Sverdrups
-#if defined key_diaeiv
-            v_msf_eiv(:,:,:) = v_msf_eiv(:,:,:) * rc_sv
-#endif
-         ENDIF
-         !
-         CALL dia_ptr_wri( kt )                        ! outputs
-         !
-      ENDIF
-      !
-#if defined key_mpp_mpi
-      IF( kt == nitend .AND. l_znl_root )   CALL histclo( numptr )      ! Close the file
-#else
-      IF( kt == nitend )                    CALL histclo( numptr )      ! Close the file
-#endif
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr')
-      !
-   END SUBROUTINE dia_ptr
-
-
-   SUBROUTINE dia_ptr_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dia_ptr_init  ***
-      !!                   
-      !! ** Purpose :   Initialization, namelist read
-      !!----------------------------------------------------------------------
-      INTEGER ::   jn           ! dummy loop indices 
-      INTEGER ::   inum, ierr   ! local integers
-      INTEGER ::   ios          ! Local integer output status for namelist read
-#if defined key_mpp_mpi
-      INTEGER, DIMENSION(1) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
-#endif
-      !!
-      NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, ln_ptrcomp, nn_fptr, nn_fwri
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
-      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport
-      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namptr )
-
-      IF(lwp) THEN                     ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namptr : set ptr parameters'
-         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr
-         WRITE(numout,*) '      Overturning heat & salt transport                  ln_ptrcomp = ', ln_ptrcomp
-         WRITE(numout,*) '      T & S zonal mean and meridional stream function    ln_diaznl  = ', ln_diaznl 
-         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas
-         WRITE(numout,*) '      Frequency of computation                           nn_fptr    = ', nn_fptr
-         WRITE(numout,*) '      Frequency of outputs                               nn_fwri    = ', nn_fwri
-      ENDIF
-      
-      IF( ln_diaptr) THEN  
-     
-         IF( nn_timing == 1 )   CALL timing_start('dia_ptr_init')
-      
-         IF( ln_subbas ) THEN   ;   nptr = 5       ! Global, Atlantic, Pacific, Indian, Indo-Pacific
-         ELSE                   ;   nptr = 1       ! Global only
-         ENDIF
-
-         !                                      ! allocate dia_ptr arrays
-         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
-
-         rc_pwatt = rc_pwatt * rau0 * rcp          ! conversion from K.s-1 to PetaWatt
-
-         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
-
-         IF( ln_subbas ) THEN                ! load sub-basin mask
-            CALL iom_open( 'subbasins', inum )
-            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
-            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
-            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
-            CALL iom_close( inum )
-            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
-            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
-            ELSE WHERE                     ;   btm30(:,:) = tmask(:,:,1)
-            END WHERE
-         ENDIF
-         btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
-      
-         DO jn = 1, nptr
-            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
-         END DO
-      
-         IF( lk_vvl )   CALL ctl_stop( 'diaptr: error in vvl case as constant i-mean surface is used' )
-
-         !                                   ! i-sum of e1v*e3v surface and its inverse
-         DO jn = 1, nptr
-            sjk(:,:,jn) = ptr_tjk( tmask(:,:,:), btmsk(:,:,jn) )
-            r1_sjk(:,:,jn) = 0._wp
-            WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
-         END DO
-
-      ! Initialise arrays to zero because diatpr is called before they are first calculated
-      ! Note that this means diagnostics will not be exactly correct when model run is restarted.
-      htr_adv(:) = 0._wp ; str_adv(:) =  0._wp ;  htr_ldf(:) = 0._wp ; str_ldf(:) =  0._wp
-
-#if defined key_mpp_mpi 
-         iglo (1) = jpjglo                   ! MPP case using MPI  ('key_mpp_mpi')
-         iloc (1) = nlcj
-         iabsf(1) = njmppt(narea)
-         iabsl(:) = iabsf(:) + iloc(:) - 1
-         ihals(1) = nldj - 1
-         ihale(1) = nlcj - nlej
-         idid (1) = 2
-         CALL flio_dom_set( jpnj, nproc/jpni, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_ptr )
-#else
-         nidom_ptr = FLIO_DOM_NONE
-#endif
-      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr_init')
-      !
-      ENDIF 
-      ! 
-   END SUBROUTINE dia_ptr_init
-
-
-   SUBROUTINE dia_ptr_wri( kt )
-      !!---------------------------------------------------------------------
-      !!                ***  ROUTINE dia_ptr_wri  ***
-      !!
-      !! ** Purpose :   output of poleward fluxes
-      !!
-      !! ** Method  :   NetCDF file
-      !!----------------------------------------------------------------------
-      !!
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !!
-      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
-      !!
-      CHARACTER (len=40) ::   clhstnam, clop, clop_once, cl_comment   ! temporary names
-      INTEGER            ::   iline, it, itmod, ji, jj, jk            !
-#if defined key_iomput
-      INTEGER            ::   inum                                    ! temporary logical unit
-#endif
-      REAL(wp)           ::   zsto, zout, zdt, zjulian                ! temporary scalars
-      !!
-      REAL(wp), POINTER, DIMENSION(:)   ::   zphi, zfoo    ! 1D workspace
-      REAL(wp), POINTER, DIMENSION(:,:) ::   z_1           ! 2D workspace
-      !!-------------------------------------------------------------------- 
-      !
-      CALL wrk_alloc( jpj      , zphi , zfoo )
-      CALL wrk_alloc( jpj , jpk, z_1 )
-
-      ! define time axis
-      it    = kt / nn_fptr
-      itmod = kt - nit000 + 1
-      
-      ! Initialization
-      ! --------------
-      IF( kt == nit000 ) THEN
-         niter = ( nit000 - 1 ) / nn_fptr
-         zdt = rdt
-         IF( nacc == 1 )   zdt = rdtmin
-         !
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'dia_ptr_wri : poleward transport and msf writing: initialization , niter = ', niter
-            WRITE(numout,*) '~~~~~~~~~~~~'
-         ENDIF
-
-         ! Reference latitude (used in plots)
-         ! ------------------
-         !                                           ! =======================
-         IF( cp_cfg == "orca" ) THEN                 !   ORCA configurations
-            !                                        ! =======================
-            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(1:jpj) = 0._wp
-            DO ji = mi0(iline), mi1(iline) 
-               zphi(1:jpj) = 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(mj0(jpjdta-1)) + ( zphi(mj0(jpjdta-1))-zphi(mj0(jpjdta-2)) ) * 0.5_wp
-                     zphi( jj ) = MIN( zphi(jj), 90._wp )
-                  END DO
-               END IF
-               IF( jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) THEN
-                  DO jj = mj0(jpjdta-1), mj1(jpjdta-1) 
-                     zphi( jj ) = 88.5_wp
-                  END DO
-                  DO jj = mj0(jpjdta  ), mj1(jpjdta  ) 
-                     zphi( jj ) = 89.5_wp
-                  END DO
-               END IF
-            END DO
-            ! provide the correct zphi to all local domains
-#if defined key_mpp_mpi
-            CALL mpp_sum( zphi, jpj, ncomm_znl )        
-#endif
-            !                                        ! =======================
-         ELSE                                        !   OTHER configurations 
-            !                                        ! =======================
-            zphi(1:jpj) = gphiv(1,:)             ! assume lat/lon coordinate, select the first i-line
-            !
-         ENDIF
-         !
-         ! Work only on westmost processor (will not work if mppini2 is used)
-#if defined key_mpp_mpi
-         IF( l_znl_root ) THEN 
-#endif
-            !
-            ! OPEN netcdf file 
-            ! ----------------
-            ! Define frequency of output and means
-            zsto = nn_fptr * 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 = nn_fwri * zdt
-            zfoo(1:jpj) = 0._wp
-
-            CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )  ! Compute julian date from starting date of the run
-            zjulian = zjulian - adatrj                         ! set calendar origin to the beginning of the experiment
-
-#if defined key_iomput
-            ! Requested by IPSL people, use by their postpro...
-            IF(lwp) THEN
-               CALL dia_nam( clhstnam, nn_fwri,' ' )
-               CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
-               WRITE(inum,*) clhstnam
-               CLOSE(inum)
-            ENDIF
-#endif
-
-            CALL dia_nam( clhstnam, nn_fwri, '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*nn_fptr, nhoridz, numptr, domain_id=nidom_ptr)
-            ! Vertical grids : gdept_1d, gdepw_1d
-            CALL histvert( numptr, "deptht", "Vertical T levels",   &
-               &                   "m", jpk, gdept_1d, ndepidzt, "down" )
-            CALL histvert( numptr, "depthw", "Vertical W levels",   &
-               &                   "m", jpk, gdepw_1d, ndepidzw, "down" )
-            !
-            CALL wheneq ( jpj*jpk, MIN(sjk(:,:,1), 1._wp), 1, 1., ndex  , ndim  )      ! Lat-Depth
-            CALL wheneq ( jpj    , MIN(sjk(:,1,1), 1._wp), 1, 1., ndex_h, ndim_h )     ! Lat
-
-            IF( ln_subbas ) THEN
-               z_1(:,1) = 1._wp
-               WHERE ( gphit(jpi/2,:) < -30._wp )   z_1(:,1) = 0._wp
-               DO jk = 2, jpk
-                  z_1(:,jk) = z_1(:,1)
-               END DO
-               !                       ! Atlantic (jn=2)
-               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,2)         , 1._wp), 1, 1., ndex_atl     , ndim_atl      ) ! Lat-Depth
-               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,2)*z_1(:,:), 1._wp), 1, 1., ndex_atl_30  , ndim_atl_30   ) ! Lat-Depth
-               CALL wheneq ( jpj    , MIN(sjk(:,1,2)*z_1(:,1), 1._wp), 1, 1., ndex_h_atl_30, ndim_h_atl_30 ) ! Lat
-               !                       ! Pacific (jn=3)
-               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,3)         , 1._wp), 1, 1., ndex_pac     , ndim_pac      ) ! Lat-Depth
-               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,3)*z_1(:,:), 1._wp), 1, 1., ndex_pac_30  , ndim_pac_30   ) ! Lat-Depth
-               CALL wheneq ( jpj    , MIN(sjk(:,1,3)*z_1(:,1), 1._wp), 1, 1., ndex_h_pac_30, ndim_h_pac_30 ) ! Lat
-               !                       ! Indian (jn=4)
-               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,4)         , 1._wp), 1, 1., ndex_ind     , ndim_ind      ) ! Lat-Depth
-               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,4)*z_1(:,:), 1._wp), 1, 1., ndex_ind_30  , ndim_ind_30   ) ! Lat-Depth
-               CALL wheneq ( jpj    , MIN(sjk(:,1,4)*z_1(:,1), 1._wp), 1, 1., ndex_h_ind_30, ndim_h_ind_30 ) ! Lat
-               !                       ! Indo-Pacific (jn=5)
-               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,5)         , 1._wp), 1, 1., ndex_ipc     , ndim_ipc      ) ! Lat-Depth
-               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,5)*z_1(:,:), 1._wp), 1, 1., ndex_ipc_30  , ndim_ipc_30   ) ! Lat-Depth
-               CALL wheneq ( jpj    , MIN(sjk(:,1,5)*z_1(:,1), 1._wp), 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
-            IF( ln_diaznl ) THEN             !  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 )
-
-               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
-#if defined key_mpp_mpi
-      END IF
-#endif
-
-#if defined key_mpp_mpi
-      IF( MOD( itmod, nn_fptr ) == 0 .AND. l_znl_root ) THEN
-#else
-      IF( MOD( itmod, nn_fptr ) == 0  ) THEN
-#endif
-         niter = niter + 1
-
-         IF( ln_diaznl ) THEN 
-            CALL histwrite( numptr, "zosrfglo", niter, sjk  (:,:,1) , ndim, ndex )
-            CALL histwrite( numptr, "zotemglo", niter, tn_jk(:,:,1)  , ndim, ndex )
-            CALL histwrite( numptr, "zosalglo", niter, sn_jk(:,:,1)  , ndim, ndex )
-
-            IF (ln_subbas) THEN 
-               CALL histwrite( numptr, "zosrfatl", niter, sjk(:,:,2), ndim_atl, ndex_atl )
-               CALL histwrite( numptr, "zosrfpac", niter, sjk(:,:,3), ndim_pac, ndex_pac )
-               CALL histwrite( numptr, "zosrfind", niter, sjk(:,:,4), ndim_ind, ndex_ind )
-               CALL histwrite( numptr, "zosrfipc", niter, sjk(:,:,5), ndim_ipc, ndex_ipc )
-
-               CALL histwrite( numptr, "zotematl", niter, tn_jk(:,:,2)  , ndim_atl, ndex_atl )
-               CALL histwrite( numptr, "zosalatl", niter, sn_jk(:,:,2)  , ndim_atl, ndex_atl )
-               CALL histwrite( numptr, "zotempac", niter, tn_jk(:,:,3)  , ndim_pac, ndex_pac )
-               CALL histwrite( numptr, "zosalpac", niter, sn_jk(:,:,3)  , ndim_pac, ndex_pac )
-               CALL histwrite( numptr, "zotemind", niter, tn_jk(:,:,4)  , ndim_ind, ndex_ind )
-               CALL histwrite( numptr, "zosalind", niter, sn_jk(:,:,4)  , ndim_ind, ndex_ind )
-               CALL histwrite( numptr, "zotemipc", niter, tn_jk(:,:,5)  , ndim_ipc, ndex_ipc )
-               CALL histwrite( numptr, "zosalipc", niter, sn_jk(:,:,5)  , ndim_ipc, ndex_ipc )
-            END IF
-         ENDIF
-
-         ! overturning outputs:
-         CALL histwrite( numptr, "zomsfglo", niter, v_msf(:,:,1), ndim, ndex )
-         IF( ln_subbas .AND. ln_diaznl ) THEN
-            CALL histwrite( numptr, "zomsfatl", niter, v_msf(:,:,2) , ndim_atl_30, ndex_atl_30 )
-            CALL histwrite( numptr, "zomsfpac", niter, v_msf(:,:,3) , ndim_pac_30, ndex_pac_30 )
-            CALL histwrite( numptr, "zomsfind", niter, v_msf(:,:,4) , ndim_ind_30, ndex_ind_30 )
-            CALL histwrite( numptr, "zomsfipc", niter, v_msf(:,:,5) , ndim_ipc_30, ndex_ipc_30 )
-         ENDIF
-#if defined key_diaeiv
-         CALL histwrite( numptr, "zomsfeiv", niter, v_msf_eiv(:,:,1), ndim  , ndex   )
-#endif
-
-         ! heat transport outputs:
-         IF( ln_subbas ) THEN
-            CALL histwrite( numptr, "sohtatl", niter, htr(:,2)  , ndim_h_atl_30, ndex_h_atl_30 )
-            CALL histwrite( numptr, "sohtpac", niter, htr(:,3)  , ndim_h_pac_30, ndex_h_pac_30 )
-            CALL histwrite( numptr, "sohtind", niter, htr(:,4)  , ndim_h_ind_30, ndex_h_ind_30 )
-            CALL histwrite( numptr, "sohtipc", niter, htr(:,5)  , ndim_h_ipc_30, ndex_h_ipc_30 )
-            CALL histwrite( numptr, "sostatl", niter, str(:,2)  , ndim_h_atl_30, ndex_h_atl_30 )
-            CALL histwrite( numptr, "sostpac", niter, str(:,3)  , ndim_h_pac_30, ndex_h_pac_30 )
-            CALL histwrite( numptr, "sostind", niter, str(:,4)  , ndim_h_ind_30, ndex_h_ind_30 )
-            CALL histwrite( numptr, "sostipc", niter, str(:,5)  , ndim_h_ipc_30, ndex_h_ipc_30 )
-         ENDIF
-
-         CALL histwrite( numptr, "sophtadv", niter, htr_adv     , ndim_h, ndex_h )
-         CALL histwrite( numptr, "sophtldf", niter, htr_ldf     , ndim_h, ndex_h )
-         CALL histwrite( numptr, "sopstadv", niter, str_adv     , ndim_h, ndex_h )
-         CALL histwrite( numptr, "sopstldf", niter, str_ldf     , ndim_h, ndex_h )
-         IF( ln_ptrcomp ) THEN 
-            CALL histwrite( numptr, "sopstove", niter, str_ove(:) , ndim_h, ndex_h )
-            CALL histwrite( numptr, "sophtove", niter, htr_ove(:) , ndim_h, ndex_h )
-         ENDIF
-#if defined key_diaeiv
-         CALL histwrite( numptr, "sophteiv", niter, htr_eiv(:,1)  , ndim_h, ndex_h )
-         CALL histwrite( numptr, "sopsteiv", niter, str_eiv(:,1)  , ndim_h, ndex_h )
-#endif
-         !
-      ENDIF
-      !
-      CALL wrk_dealloc( jpj      , zphi , zfoo )
-      CALL wrk_dealloc( jpj , jpk, z_1 )
-      !
-  END SUBROUTINE dia_ptr_wri
+   END FUNCTION ptr_sjk
+
 
    !!======================================================================
Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 4756)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 5260)
@@ -44,6 +44,4 @@
    USE in_out_manager  ! I/O manager
    USE diadimg         ! dimg direct access file format output
-   USE diaar5, ONLY :   lk_diaar5
-   USE dynadv, ONLY :   ln_dynadv_vec
    USE diatmb          ! Top,middle,bottom output
    USE dia25h          ! 25h Mean output
@@ -82,5 +80,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id $
+   !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -91,7 +89,5 @@
       INTEGER, DIMENSION(2) :: ierr
       !!----------------------------------------------------------------------
-      !
       ierr = 0
-      !
       ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
          &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
@@ -133,5 +129,5 @@
       REAL(wp)                     ::   zztmp, zztmpx, zztmpy   ! 
       !!
-      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d       ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
       !!----------------------------------------------------------------------
@@ -148,82 +144,190 @@
       ENDIF
 
-      IF( lk_vvl ) THEN
-         z3d(:,:,:) = tsn(:,:,:,jp_tem) * fse3t_n(:,:,:)
-         CALL iom_put( "toce" , z3d                        )   ! heat content
-         CALL iom_put( "sst"  , z3d(:,:,1)                 )   ! sea surface heat content
-         z3d(:,:,1) = tsn(:,:,1,jp_tem) * z3d(:,:,1)
-         CALL iom_put( "sst2" , z3d(:,:,1)                 )   ! sea surface content of squared temperature
-         z3d(:,:,:) = tsn(:,:,:,jp_sal) * fse3t_n(:,:,:)            
-         CALL iom_put( "soce" , z3d                        )   ! salinity content
-         CALL iom_put( "sss"  , z3d(:,:,1)                 )   ! sea surface salinity content
-         z3d(:,:,1) = tsn(:,:,1,jp_sal) * z3d(:,:,1)
-         CALL iom_put( "sss2" , z3d(:,:,1)                 )   ! sea surface content of squared salinity
-      ELSE
-         CALL iom_put( "toce" , tsn(:,:,:,jp_tem)          )   ! temperature
-         CALL iom_put( "sst"  , tsn(:,:,1,jp_tem)          )   ! sea surface temperature
-         CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature
-         CALL iom_put( "soce" , tsn(:,:,:,jp_sal)          )   ! salinity
-         CALL iom_put( "sss"  , tsn(:,:,1,jp_sal)          )   ! sea surface salinity
-         CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity
-      END IF
-      IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN
-         CALL iom_put( "uoce" , un(:,:,:) * fse3u_n(:,:,:) )    ! i-transport
-         CALL iom_put( "voce" , vn(:,:,:) * fse3v_n(:,:,:) )    ! j-transport
-      ELSE
-         CALL iom_put( "uoce" , un                         )    ! i-current
-         CALL iom_put( "voce" , vn                         )    ! j-current
-      END IF
-      CALL iom_put(    "avt"  , avt                        )    ! T vert. eddy diff. coef.
-      CALL iom_put(    "avm"  , avmu                       )    ! T vert. eddy visc. coef.
-      IF( lk_zdfddm ) THEN
-         CALL iom_put( "avs" , fsavs(:,:,:)                          )    ! S vert. eddy diff. coef.
-      ENDIF
-
-      DO jj = 2, jpjm1                                    ! sst gradient
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zztmp      = tsn(ji,jj,1,jp_tem)
-            zztmpx     = ( tsn(ji+1,jj  ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) / e1u(ji-1,jj  )
-            zztmpy     = ( tsn(ji  ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) / e2v(ji  ,jj-1)
-            z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
-               &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
-         END DO
-      END DO
-      CALL lbc_lnk( z2d, 'T', 1. )
-      CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
-!CDIR NOVERRCHK
-      z2d(:,:) = SQRT( z2d(:,:) )
-      CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
-
-      IF( lk_diaar5 ) THEN
+      IF( .NOT.lk_vvl ) THEN
+         CALL iom_put( "e3t" , fse3t_n(:,:,:) )
+         CALL iom_put( "e3u" , fse3u_n(:,:,:) )
+         CALL iom_put( "e3v" , fse3v_n(:,:,:) )
+         CALL iom_put( "e3w" , fse3w_n(:,:,:) )
+      ENDIF
+      
+      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
+      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
+      IF ( iom_use("sbt") ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_tem)
+            END DO
+         END DO
+         CALL iom_put( "sbt", z2d )                ! bottom temperature
+      ENDIF
+      
+      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
+      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
+      IF ( iom_use("sbs") ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_sal)
+            END DO
+         END DO
+         CALL iom_put( "sbs", z2d )                ! bottom salinity
+      ENDIF
+         
+      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current
+      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current
+      IF ( iom_use("sbu") ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               z2d(ji,jj) = un(ji,jj,MAX(mbathy(ji,jj),1))
+            END DO
+         END DO
+         CALL iom_put( "sbu", z2d )                ! bottom i-current
+      ENDIF
+      
+      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
+      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
+      IF ( iom_use("sbv") ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               z2d(ji,jj) = vn(ji,jj,MAX(mbathy(ji,jj),1))
+            END DO
+         END DO
+         CALL iom_put( "sbv", z2d )                ! bottom j-current
+      ENDIF
+
+      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
+      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
+      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
+
+      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
+         DO jj = 2, jpjm1                                    ! sst gradient
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zztmp      = tsn(ji,jj,1,jp_tem)
+               zztmpx     = ( tsn(ji+1,jj  ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) / e1u(ji-1,jj  )
+               zztmpy     = ( tsn(ji  ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) / e2v(ji  ,jj-1)
+               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
+                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
+            END DO
+         END DO
+         CALL lbc_lnk( z2d, 'T', 1. )
+         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
+         z2d(:,:) = SQRT( z2d(:,:) )
+         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
+      ENDIF
+         
+      ! clem: heat and salt content
+      IF( iom_use("heatc") ) THEN
+         z2d(:,:)  = 0._wp 
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2)
+      ENDIF
+
+      IF( iom_use("saltc") ) THEN
+         z2d(:,:)  = 0._wp 
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2)
+      ENDIF
+      !
+      IF ( iom_use("eken") ) THEN
+         rke(:,:,jk) = 0._wp                               !      kinetic energy 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zztmp   = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk)    &
+                     &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * fse3u(ji  ,jj,jk) )  &
+                     &          *  zztmp 
+                  !
+                  zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk)    &
+                     &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) )  &
+                     &          *  zztmp 
+                  !
+                  rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
+                  !
+               ENDDO
+            ENDDO
+         ENDDO
+         CALL lbc_lnk( rke, 'T', 1. )
+         CALL iom_put( "eken", rke )           
+      ENDIF
+         
+      IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
          z3d(:,:,jpk) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk)
+            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
          END DO
          CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
-         zztmp = 0.5 * rcp
+      ENDIF
+      
+      IF( iom_use("u_heattr") ) THEN
          z2d(:,:) = 0.e0 
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
                END DO
             END DO
          END DO
          CALL lbc_lnk( z2d, 'U', -1. )
-         CALL iom_put( "u_heattr", z2d )                  ! heat transport in i-direction
-         DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk)
-         END DO
-         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
+         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
+      ENDIF
+
+      IF( iom_use("u_salttr") ) THEN
          z2d(:,:) = 0.e0 
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
                END DO
             END DO
          END DO
+         CALL lbc_lnk( z2d, 'U', -1. )
+         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
+      ENDIF
+
+      
+      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
+         z3d(:,:,jpk) = 0.e0
+         DO jk = 1, jpkm1
+            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
+         END DO
+         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
+      ENDIF
+      
+      IF( iom_use("v_heattr") ) THEN
+         z2d(:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
          CALL lbc_lnk( z2d, 'V', -1. )
-         CALL iom_put( "v_heattr", z2d )                  !  heat transport in i-direction
+         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction
+      ENDIF
+
+      IF( iom_use("v_salttr") ) THEN
+         z2d(:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( z2d, 'V', -1. )
+         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
       ENDIF
       !
@@ -500,23 +604,22 @@
          ENDIF
 
-#if ! defined key_coupled 
-         CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-#endif
-
-
-
-#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 ) 
-         CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-#endif
+         IF( .NOT. lk_cpl ) THEN
+            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         ENDIF
+
+         IF( lk_cpl .AND. nn_ice <= 1 ) THEN
+            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         ENDIF
+         
          clmx ="l_max(only(x))"    ! max index on a period
          CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX 
@@ -533,15 +636,10 @@
 #endif
 
-#if defined key_coupled 
-# if defined key_lim3
-         Must be adapted to LIM3
-# endif 
-# if defined key_lim2
-         CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-# endif 
-#endif 
+         IF( lk_cpl .AND. nn_ice == 2 ) THEN
+            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         ENDIF
 
          CALL histend( nid_T, snc4chunks=snc4set )
@@ -634,5 +732,4 @@
       ENDIF
 
-      ! Write fields on T grid
       IF( lk_vvl ) THEN
          CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
@@ -645,5 +742,4 @@
          CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
          CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
-
       ENDIF
       IF( lk_vvl ) THEN
@@ -695,18 +791,18 @@
       ENDIF
 
-#if ! defined key_coupled
-      CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
-      CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-      IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
-      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
-#endif
-#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 ) 
-      CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
-      CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
+      IF( .NOT. lk_cpl ) THEN
+         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
+         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
          IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
-      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
-#endif
-      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
-      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
+         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
+      ENDIF
+      IF( lk_cpl .AND. nn_ice <= 1 ) THEN
+         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
+         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
+         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
+         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
+      ENDIF
+!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
+!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
 
 #if defined key_diahth
@@ -717,16 +813,9 @@
 #endif
 
-#if defined key_coupled 
-# if defined key_lim3
-      Must be adapted for LIM3
-      CALL histwrite( nid_T, "soicetem", it, tn_ice        , ndim_hT, ndex_hT )   ! surf. ice temperature
-      CALL histwrite( nid_T, "soicealb", it, alb_ice       , ndim_hT, ndex_hT )   ! ice albedo
-# endif
-# if defined key_lim2
-      CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
-      CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
-# endif
-#endif
-         ! Write fields on U grid
+      IF( lk_cpl .AND. nn_ice == 2 ) THEN
+         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
+         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
+      ENDIF
+
       CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
       IF( ln_traldf_gdia ) THEN
@@ -750,5 +839,4 @@
       CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
 
-         ! Write fields on V grid
       CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
       IF( ln_traldf_gdia ) THEN
@@ -765,5 +853,4 @@
       CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
 
-         ! Write fields on W grid
       CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
       IF( ln_traldf_gdia ) THEN
