Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 6140)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 6140)
@@ -0,0 +1,326 @@
+MODULE dia25h 
+   !!======================================================================
+   !!                       ***  MODULE  diaharm  ***
+   !! Harmonic analysis of tidal constituents 
+   !!======================================================================
+   !! History :  3.6  !  2014  (E O'Dea)  Original code
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE wrk_nemo        ! working arrays
+#if defined key_zdftke 
+   USE zdf_oce, ONLY: en
+#endif
+   USE zdf_oce, ONLY: avt, avm
+#if defined key_zdfgls
+   USE zdf_oce, ONLY: en
+   USE zdfgls, ONLY: mxln
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   LOGICAL , PUBLIC ::   ln_dia25h     !:  25h mean output
+   PUBLIC   dia_25h_init               ! routine called by nemogcm.F90
+   PUBLIC   dia_25h                    ! routine called by diawri.F90
+
+  !! * variables for calculating 25-hourly means
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h 
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_25h , avm_25h
+#if defined key_zdfgls || key_zdftke
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_25h
+#endif
+#if defined key_zdfgls 
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   rmxln_25h
+#endif
+   INTEGER, SAVE :: cnt_25h     ! Counter for 25 hour means
+
+
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_25h_init 
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_25h_init  ***
+      !!     
+      !! ** Purpose: Initialization of 25h mean namelist 
+      !!        
+      !! ** Method : Read namelist
+      !!   History
+      !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_25h
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER ::   ios                 ! Local integer output status for namelist read
+      INTEGER ::   ierror              ! Local integer for memory allocation
+      !
+      NAMELIST/nam_dia25h/ ln_dia25h
+      !!----------------------------------------------------------------------
+      !
+      REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
+      READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics
+      READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_dia25h )
+
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_25h_init : Output 25 hour mean diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) 'Namelist nam_dia25h : set 25h outputs '
+         WRITE(numout,*) 'Switch for 25h diagnostics (T) or not (F)  ln_dia25h  = ', ln_dia25h
+      ENDIF
+      IF( .NOT. ln_dia25h )   RETURN
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+      ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate tn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate sn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate un_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate vn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate wn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate avt_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate avm_25h' )   ;   RETURN
+      ENDIF
+# if defined key_zdfgls || defined key_zdftke
+      ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN
+      ENDIF
+#endif
+# if defined key_zdfgls 
+      ALLOCATE( rmxln_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate rmxln_25h' )   ;   RETURN
+      ENDIF
+#endif
+      ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate sshn_25h' )   ;   RETURN
+      ENDIF
+      ! ------------------------- !
+      ! 2 - Assign Initial Values !
+      ! ------------------------- !
+      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
+      tn_25h(:,:,:) = tsb(:,:,:,jp_tem)
+      sn_25h(:,:,:) = tsb(:,:,:,jp_sal)
+      sshn_25h(:,:) = sshb(:,:)
+      un_25h(:,:,:) = ub(:,:,:)
+      vn_25h(:,:,:) = vb(:,:,:)
+      wn_25h(:,:,:) = wn(:,:,:)
+      avt_25h(:,:,:) = avt(:,:,:)
+      avm_25h(:,:,:) = avm(:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+         en_25h(:,:,:) = en(:,:,:)
+#endif
+# if defined key_zdfgls
+         rmxln_25h(:,:,:) = mxln(:,:,:)
+#endif
+#if defined key_lim3 || defined key_lim2
+         CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice')
+#endif 
+
+      ! -------------------------- !
+      ! 3 - Return to dia_wri      !
+      ! -------------------------- !
+
+
+   END SUBROUTINE dia_25h_init
+
+
+   SUBROUTINE dia_25h( kt )  
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_25h  ***
+      !!         
+      !!
+      !!--------------------------------------------------------------------
+      !!                   
+      !! ** Purpose :   Write diagnostics with M2/S2 tide removed
+      !!
+      !! ** Method  :   
+      !!      25hr mean outputs for shelf seas
+      !!
+      !! History :
+      !!   ?.0  !  07-04  (A. Hines) New routine, developed from dia_wri_foam
+      !!   3.4  !  02-13  (J. Siddorn) Routine taken from old dia_wri_foam
+      !!   3.6  !  08-14  (E. O'Dea) adapted for VN3.6
+      !!----------------------------------------------------------------------
+      !! * Modules used
+
+      IMPLICIT NONE
+
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+
+
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk
+
+      LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout
+      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi       ! temporary reals
+      INTEGER                          ::   i_steps                               ! no of timesteps per hour
+      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                    ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                                  ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace
+      INTEGER                          ::   iyear0, nimonth0,iday0                ! start year,imonth,day
+
+      !!----------------------------------------------------------------------
+
+      ! 0. Initialisation
+      ! -----------------
+      ! Define frequency of summing to create 25 h mean
+      IF( MOD( 3600,INT(rdt) ) == 0 ) THEN
+         i_steps = 3600/INT(rdt)
+      ELSE
+         CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
+      ENDIF
+
+#if defined key_lim3 || defined key_lim2
+      CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice')
+#endif
+
+      ! local variable for debugging
+      ll_print = ll_print .AND. lwp
+
+      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours
+      ! every day
+      IF( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 ) THEN
+
+         IF (lwp) THEN
+              WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt
+              WRITE(numout,*) '~~~~~~~~~~~~ '
+         ENDIF
+
+         tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem)
+         sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal)
+         sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:)
+         un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:)
+         vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:)
+         wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:)
+         avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:)
+         avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+         en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
+#endif
+# if defined key_zdfgls
+         rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:)
+#endif
+         cnt_25h = cnt_25h + 1
+
+         IF (lwp) THEN
+            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h
+         ENDIF
+
+      ENDIF ! MOD( kt, i_steps ) == 0
+
+         ! Write data for 25 hour mean output streams
+      IF( cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
+
+            IF(lwp) THEN
+               WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
+               WRITE(numout,*) '~~~~~~~~~~~~ '
+            ENDIF
+
+            tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp
+            sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp
+            sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp
+            un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp
+            vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp
+            wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp
+            avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp
+            avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp
+# if defined key_zdfgls || defined key_zdftke
+            en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp
+#endif
+# if defined key_zdfgls
+            rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp
+#endif
+
+            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
+            zmdi=1.e+20 !missing data indicator for masking
+            ! write tracers (instantaneous)
+            zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("temper25h", zw3d)   ! potential temperature
+            zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put( "salin25h", zw3d  )   ! salinity
+            zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+            CALL iom_put( "ssh25h", zw2d )   ! sea surface 
+
+
+            ! Write velocities (instantaneous)
+            zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
+            CALL iom_put("vozocrtx25h", zw3d)    ! i-current
+            zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
+            CALL iom_put("vomecrty25h", zw3d  )   ! j-current
+
+            zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("vomecrtz25h", zw3d )   ! k-current
+            zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("avt25h", zw3d )   ! diffusivity
+            zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("avm25h", zw3d)   ! viscosity
+#if defined key_zdftke || defined key_zdfgls 
+            zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("tke25h", zw3d)   ! tke
+#endif
+#if defined key_zdfgls 
+            zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put( "mxln25h",zw3d)
+#endif
+
+            ! After the write reset the values to cnt=1 and sum values equal current value 
+            tn_25h(:,:,:) = tsn(:,:,:,jp_tem)
+            sn_25h(:,:,:) = tsn(:,:,:,jp_sal)
+            sshn_25h(:,:) = sshn (:,:)
+            un_25h(:,:,:) = un(:,:,:)
+            vn_25h(:,:,:) = vn(:,:,:)
+            wn_25h(:,:,:) = wn(:,:,:)
+            avt_25h(:,:,:) = avt(:,:,:)
+            avm_25h(:,:,:) = avm(:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+            en_25h(:,:,:) = en(:,:,:)
+#endif
+# if defined key_zdfgls
+            rmxln_25h(:,:,:) = mxln(:,:,:)
+#endif
+            cnt_25h = 1
+            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h
+
+      ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000
+
+
+   END SUBROUTINE dia_25h 
+
+   !!======================================================================
+END MODULE dia25h
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 5930)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 6140)
@@ -40,6 +40,4 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity
       
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -99,12 +97,12 @@
       ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
       ztsn(:,:,:,jp_sal) = sn0(:,:,:)
-      CALL eos( ztsn, zrhd, fsdept_n(:,:,:) )                       ! now in situ density using initial salinity
+      CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity
       !
       zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
       DO jk = 1, jpkm1
-         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
-      END DO
-      IF( .NOT.lk_vvl ) THEN
-         IF ( ln_isfcav ) THEN
+         zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+      END DO
+      IF( ln_linssh ) THEN
+         IF( ln_isfcav ) THEN
             DO ji=1,jpi
                DO jj=1,jpj
@@ -115,4 +113,7 @@
             zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
          END IF
+!!gm
+!!gm   riceload should be added in both ln_linssh=T or F, no?
+!!gm
       END IF
       !                                         
@@ -123,5 +124,5 @@
       
       !                                         ! steric sea surface height
-      CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) )                 ! now in situ and potential density
+      CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
       zrhop(:,:,jpk) = 0._wp
       CALL iom_put( 'rhop', zrhop )
@@ -129,7 +130,7 @@
       zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
       DO jk = 1, jpkm1
-         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
-      END DO
-      IF( .NOT.lk_vvl ) THEN
+         zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+      END DO
+      IF( ln_linssh ) THEN
          IF ( ln_isfcav ) THEN
             DO ji=1,jpi
@@ -159,5 +160,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zztmp = area(ji,jj) * fse3t(ji,jj,jk)
+               zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
                ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
                zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal)
@@ -165,6 +166,6 @@
          END DO
       END DO
-      IF( .NOT.lk_vvl ) THEN
-         IF ( ln_isfcav ) THEN
+      IF( ln_linssh ) THEN
+         IF( ln_isfcav ) THEN
             DO ji=1,jpi
                DO jj=1,jpj
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 6140)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 6140)
@@ -0,0 +1,206 @@
+MODULE diacfl
+   !!==============================================================================
+   !!                       ***  MODULE  diacfl  ***
+   !! Output CFL diagnostics to ascii file
+   !!==============================================================================
+   !! History :  1.0  !  2010-03  (E. Blockley)  Original code
+   !!                 !  2014-06  (T Graham) Removed CPP key & Updated to vn3.6
+   !! 
+   !!----------------------------------------------------------------------
+   !!   dia_cfl        : Compute and output Courant numbers at each timestep
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers
+   USE dom_oce         ! ocean space and time domain
+   USE lib_mpp         ! distribued memory computing
+   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
+   USE in_out_manager  ! I/O manager
+   USE domvvl     
+   USE timing          ! Performance output
+
+   IMPLICIT NONE
+   PRIVATE
+
+   REAL(wp) :: cu_max, cv_max, cw_max                      ! Run max U Courant number 
+   INTEGER, DIMENSION(3) :: cu_loc, cv_loc, cw_loc         ! Run max locations
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcu_cfl           ! Courant number arrays
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcv_cfl           ! Courant number arrays
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcw_cfl           ! Courant number arrays
+
+   INTEGER  :: numcfl                                       ! outfile unit
+   CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii"      ! ascii filename
+
+   PUBLIC   dia_cfl       ! routine called by step.F90
+   PUBLIC   dia_cfl_init  ! routine called by nemogcm
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+
+CONTAINS
+
+
+   SUBROUTINE dia_cfl ( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_cfl  ***
+      !!
+      !! ** Purpose :  Compute the Courant numbers Cu=u*dt/dx and Cv=v*dt/dy
+      !!               and output to ascii file 'cfl_diagnostics.ascii'
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT(in) ::  kt                            ! ocean time-step index
+
+      REAL(wp) :: zcu_max, zcv_max, zcw_max                 ! max Courant numbers per timestep
+      INTEGER, DIMENSION(3) :: zcu_loc, zcv_loc, zcw_loc    ! max Courant number locations
+
+      REAL(wp) :: dt                                        ! temporary scalars
+      INTEGER, DIMENSION(3) :: zlocu, zlocv, zlocw          ! temporary arrays 
+      INTEGER  :: ji, jj, jk                                ! dummy loop indices
+
+      
+      IF( nn_diacfl == 1) THEN
+         IF( nn_timing == 1 )   CALL timing_start('dia_cfl')
+         ! setup timestep multiplier to account for initial Eulerian timestep
+         IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    dt = rdt
+         ELSE                                        ;    dt = rdt * 2.0
+         ENDIF
+
+             ! calculate Courant numbers
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, fs_jpim1   ! vector opt.
+
+                  ! Courant number for x-direction (zonal current)
+                  zcu_cfl(ji,jj,jk) = ABS(un(ji,jj,jk))*dt/e1u(ji,jj)
+
+                  ! Courant number for y-direction (meridional current)
+                  zcv_cfl(ji,jj,jk) = ABS(vn(ji,jj,jk))*dt/e2v(ji,jj)
+
+                  ! Courant number for z-direction (vertical current)
+                  zcw_cfl(ji,jj,jk) = ABS(wn(ji,jj,jk))*dt/e3w_n(ji,jj,jk)
+               END DO
+            END DO         
+         END DO
+
+         ! calculate maximum values and locations
+         IF( lk_mpp ) THEN
+            CALL mpp_maxloc(zcu_cfl,umask,zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3))
+            CALL mpp_maxloc(zcv_cfl,vmask,zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3))
+            CALL mpp_maxloc(zcw_cfl,tmask,zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3))
+         ELSE
+            zlocu = MAXLOC( ABS( zcu_cfl(:,:,:) ) )
+            zcu_loc(1) = zlocu(1) + nimpp - 1
+            zcu_loc(2) = zlocu(2) + njmpp - 1
+            zcu_loc(3) = zlocu(3)
+            zcu_max = zcu_cfl(zcu_loc(1),zcu_loc(2),zcu_loc(3))
+
+            zlocv = MAXLOC( ABS( zcv_cfl(:,:,:) ) )
+            zcv_loc(1) = zlocv(1) + nimpp - 1
+            zcv_loc(2) = zlocv(2) + njmpp - 1
+            zcv_loc(3) = zlocv(3)
+            zcv_max = zcv_cfl(zcv_loc(1),zcv_loc(2),zcv_loc(3))
+
+            zlocw = MAXLOC( ABS( zcw_cfl(:,:,:) ) )
+            zcw_loc(1) = zlocw(1) + nimpp - 1
+            zcw_loc(2) = zlocw(2) + njmpp - 1
+            zcw_loc(3) = zlocw(3)
+            zcw_max = zcw_cfl(zcw_loc(1),zcw_loc(2),zcw_loc(3))
+         ENDIF
+      
+         ! write out to file
+         IF( lwp ) THEN
+            WRITE(numcfl,FMT='(2x,i4,5x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3)
+            WRITE(numcfl,FMT='(11x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cv', zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3)
+            WRITE(numcfl,FMT='(11x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cw', zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3)
+         ENDIF
+
+         ! update maximum Courant numbers from whole run if applicable
+         IF( zcu_max > cu_max ) THEN
+            cu_max = zcu_max
+            cu_loc = zcu_loc
+         ENDIF
+         IF( zcv_max > cv_max ) THEN
+            cv_max = zcv_max
+            cv_loc = zcv_loc
+         ENDIF
+         IF( zcw_max > cw_max ) THEN
+            cw_max = zcw_max
+            cw_loc = zcw_loc
+         ENDIF
+
+         ! at end of run output max Cu and Cv and close ascii file
+         IF( kt == nitend .AND. lwp ) THEN
+            ! to ascii file
+            WRITE(numcfl,*) '******************************************'
+            WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', cu_max, cu_loc(1), cu_loc(2), cu_loc(3)
+            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
+            WRITE(numcfl,*) '******************************************'
+            WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', cv_max, cv_loc(1), cv_loc(2), cv_loc(3)
+            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
+            WRITE(numcfl,*) '******************************************'
+            WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', cw_max, cw_loc(1), cw_loc(2), cw_loc(3)
+            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
+            CLOSE( numcfl ) 
+
+            ! to ocean output
+            WRITE(numout,*)
+            WRITE(numout,*) 'dia_cfl     : Maximum Courant number information for the run:'
+            WRITE(numout,*) '~~~~~~~~~~~~'
+            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cu', cu_max, 'at (i, j, k) = (', cu_loc(1), cu_loc(2), cu_loc(3), ')'
+            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
+            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cv', cv_max, 'at (i, j, k) = (', cv_loc(1), cv_loc(2), cv_loc(3), ')'
+            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
+            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cw', cw_max, 'at (i, j, k) = (', cw_loc(1), cw_loc(2), cw_loc(3), ')'
+            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
+
+         ENDIF
+
+         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl')
+      ENDIF
+
+   END SUBROUTINE dia_cfl
+
+   SUBROUTINE dia_cfl_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_cfl_init  ***
+      !!                   
+      !! ** Purpose :   create output file, initialise arrays
+      !!----------------------------------------------------------------------
+
+
+      IF( nn_diacfl == 1 ) THEN
+         IF( nn_timing == 1 )   CALL timing_start('dia_cfl_init')
+
+         cu_max=0.0
+         cv_max=0.0
+         cw_max=0.0
+
+         ALLOCATE( zcu_cfl(jpi, jpj, jpk), zcv_cfl(jpi, jpj, jpk), zcw_cfl(jpi, jpj, jpk) )
+
+         zcu_cfl(:,:,:)=0.0
+         zcv_cfl(:,:,:)=0.0
+         zcw_cfl(:,:,:)=0.0
+
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'dia_cfl     : Outputting CFL diagnostics to '//TRIM(clname)
+            WRITE(numout,*) '~~~~~~~~~~~~'
+            WRITE(numout,*)
+
+            ! create output ascii file
+            CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
+            WRITE(numcfl,*) 'Timestep  Direction  Max C     i    j    k'
+            WRITE(numcfl,*) '******************************************'
+         ENDIF
+
+         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl_init')
+
+      ENDIF
+
+   END SUBROUTINE dia_cfl_init
+
+END MODULE diacfl
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 5930)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 6140)
@@ -1,118 +1,110 @@
 MODULE diadct
-  !!=====================================================================
-  !!                       ***  MODULE  diadct  ***
-  !! Ocean diagnostics: Compute the transport trough a sec.
-  !!===============================================================
-  !! History : 
-  !!
-  !!         original  : 02/99 (Y Drillet)
-  !!         addition  : 10/01 (Y Drillet, R Bourdalle Badie)
-  !!                   : 10/05 (M Laborie) F90
-  !!         addition  : 04/07 (G Garric) Ice sections
-  !!         bugfix    : 04/07 (C Bricaud) test on sec%nb_point
-  !!                                      initialisation of ztransp1,ztransp2,...
-  !!         nemo_v_3_4: 09/2011 (C Bricaud)
-  !!
-  !!
-  !!----------------------------------------------------------------------
+   !!======================================================================
+   !!                       ***  MODULE  diadct  ***
+   !! Ocean diagnostics: Compute the transport trough a sec.
+   !!======================================================================
+   !! History :  OPA  ! 02/1999 (Y Drillet)  original code
+   !!                 ! 10/2001 (Y Drillet, R Bourdalle Badie)
+   !!   NEMO     1.0  ! 10/2005 (M Laborie) F90
+   !!            3.0  ! 04/2007 (G Garric) Ice sections
+   !!             -   ! 04/2007 (C Bricaud) test on sec%nb_point, initialisation of ztransp1,ztransp2,...
+   !!            3.4  ! 09/2011 (C Bricaud)
+   !!----------------------------------------------------------------------
 #if defined key_diadct
-  !!----------------------------------------------------------------------
-  !!   'key_diadct' :
-  !!----------------------------------------------------------------------
-  !!----------------------------------------------------------------------
-  !!   dia_dct      :  Compute the transport through a sec.
-  !!   dia_dct_init :  Read namelist.
-  !!   readsec      :  Read sections description and pathway
-  !!   removepoints :  Remove points which are common to 2 procs
-  !!   transport    :  Compute transport for each sections
-  !!   dia_dct_wri  :  Write tranports results in ascii files
-  !!   interp       :  Compute temperature/salinity/density at U-point or V-point
-  !!   
-  !!----------------------------------------------------------------------
-  !! * Modules used
-  USE oce             ! ocean dynamics and tracers
-  USE dom_oce         ! ocean space and time domain
-  USE phycst          ! physical constants
-  USE in_out_manager  ! I/O manager
-  USE daymod          ! calendar
-  USE dianam          ! build name of file
-  USE lib_mpp         ! distributed memory computing library
+   !!----------------------------------------------------------------------
+   !!   'key_diadct' :
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   dia_dct      :  Compute the transport through a sec.
+   !!   dia_dct_init :  Read namelist.
+   !!   readsec      :  Read sections description and pathway
+   !!   removepoints :  Remove points which are common to 2 procs
+   !!   transport    :  Compute transport for each sections
+   !!   dia_dct_wri  :  Write tranports results in ascii files
+   !!   interp       :  Compute temperature/salinity/density at U-point or V-point
+   !!   
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE in_out_manager  ! I/O manager
+   USE daymod          ! calendar
+   USE dianam          ! build name of file
+   USE lib_mpp         ! distributed memory computing library
 #if defined key_lim2
-  USE ice_2
+   USE ice_2
 #endif
 #if defined key_lim3
-  USE ice
+   USE ice
 #endif
-  USE domvvl
-  USE timing          ! preformance summary
-  USE wrk_nemo        ! working arrays
-
-  IMPLICIT NONE
-  PRIVATE
-
-  !! * Routine accessibility
-  PUBLIC   dia_dct      ! routine called by step.F90
-  PUBLIC   dia_dct_init ! routine called by opa.F90
-  PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90 
-  PRIVATE  readsec
-  PRIVATE  removepoints
-  PRIVATE  transport
-  PRIVATE  dia_dct_wri
-
-#include "domzgr_substitute.h90"
-
-  !! * Shared module variables
-  LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
-
-  !! * Module variables
-  INTEGER :: nn_dct        ! Frequency of computation
-  INTEGER :: nn_dctwri     ! Frequency of output
-  INTEGER :: nn_secdebug   ! Number of the section to debug
+   USE domvvl
+   USE timing          ! preformance summary
+   USE wrk_nemo        ! working arrays
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dia_dct      ! routine called by step.F90
+   PUBLIC   dia_dct_init ! routine called by opa.F90
+   PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90 
+   PRIVATE  readsec
+   PRIVATE  removepoints
+   PRIVATE  transport
+   PRIVATE  dia_dct_wri
+
+   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
+
+   INTEGER :: nn_dct        ! Frequency of computation
+   INTEGER :: nn_dctwri     ! Frequency of output
+   INTEGER :: nn_secdebug   ! Number of the section to debug
    
-  INTEGER, PARAMETER :: nb_class_max  = 10
-  INTEGER, PARAMETER :: nb_sec_max    = 150
-  INTEGER, PARAMETER :: nb_point_max  = 2000
-  INTEGER, PARAMETER :: nb_type_class = 10
-  INTEGER, PARAMETER :: nb_3d_vars    = 3 
-  INTEGER, PARAMETER :: nb_2d_vars    = 2 
-  INTEGER            :: nb_sec 
-
-  TYPE POINT_SECTION
-     INTEGER :: I,J
-  END TYPE POINT_SECTION
-
-  TYPE COORD_SECTION
-     REAL(wp) :: lon,lat
-  END TYPE COORD_SECTION
-
-  TYPE SECTION
-     CHARACTER(len=60)                            :: name              ! name of the sec
-     LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
+   INTEGER, PARAMETER :: nb_class_max  = 10
+   INTEGER, PARAMETER :: nb_sec_max    = 150
+   INTEGER, PARAMETER :: nb_point_max  = 2000
+   INTEGER, PARAMETER :: nb_type_class = 10
+   INTEGER, PARAMETER :: nb_3d_vars    = 3 
+   INTEGER, PARAMETER :: nb_2d_vars    = 2 
+   INTEGER            :: nb_sec 
+
+   TYPE POINT_SECTION
+      INTEGER :: I,J
+   END TYPE POINT_SECTION
+
+   TYPE COORD_SECTION
+      REAL(wp) :: lon,lat
+   END TYPE COORD_SECTION
+
+   TYPE SECTION
+      CHARACTER(len=60)                            :: name              ! name of the sec
+      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
                                                                        ! heat transports
-     LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation
-     LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
-     TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec
-     INTEGER                                      :: nb_class          ! number of boundaries for density classes
-     INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
-     CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
-     REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
-                                                     zsigp           ,&! potential density classes    (99 if you don't want)
-                                                     zsal            ,&! salinity classes   (99 if you don't want)
-                                                     ztem            ,&! temperature classes(99 if you don't want)
-                                                     zlay              ! level classes      (99 if you don't want)
-     REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
-     REAL(wp)                                         :: slopeSection  ! slope of the section
-     INTEGER                                          :: nb_point      ! number of points in the section
-     TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections
-  END TYPE SECTION
-
-  TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
- 
-  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
-  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d  
-
+      LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation
+      LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
+      TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec
+      INTEGER                                      :: nb_class          ! number of boundaries for density classes
+      INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
+      CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
+      REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
+                                                      zsigp           ,&! potential density classes    (99 if you don't want)
+                                                      zsal            ,&! salinity classes   (99 if you don't want)
+                                                      ztem            ,&! temperature classes(99 if you don't want)
+                                                      zlay              ! level classes      (99 if you don't want)
+      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
+      REAL(wp)                                         :: slopeSection  ! slope of the section
+      INTEGER                                          :: nb_point      ! number of points in the section
+      TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections
+   END TYPE SECTION
+
+   TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
+ 
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d  
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
-
  
   INTEGER FUNCTION diadct_alloc() 
@@ -130,4 +122,5 @@
  
   END FUNCTION diadct_alloc 
+
 
   SUBROUTINE dia_dct_init
@@ -208,8 +201,6 @@
      !!               Reinitialise all relevant arrays to zero 
      !!---------------------------------------------------------------------
-     !! * Arguments
-     INTEGER,INTENT(IN)        ::kt
-
-     !! * Local variables
+     INTEGER,INTENT(in)        ::kt
+     !
      INTEGER             :: jsec,            &! loop on sections
                             itotal            ! nb_sec_max*nb_type_class*nb_class_max
@@ -220,6 +211,6 @@
      REAL(wp), POINTER, DIMENSION(:)    :: zwork !   "  
      REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   "
-
      !!---------------------------------------------------------------------    
+     !
      IF( nn_timing == 1 )   CALL timing_start('dia_dct')
 
@@ -619,142 +610,140 @@
                             zumid_ice, zvmid_ice,                &!U/V ice velocity 
                             zTnorm                                !transport of velocity through one cell's sides 
-     REAL(wp)            :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep !temperature/salinity/potential density/ssh/depth at u/v point
+     REAL(wp)            :: ztn, zsn, zrhoi, zrhop, zsshn, zdep !temperature/salinity/potential density/ssh/depth at u/v point
 
      TYPE(POINT_SECTION) :: k
-     !!--------------------------------------------------------
-
-     IF( ld_debug )WRITE(numout,*)'      Compute transport'
-
-     !---------------------------!
-     !  COMPUTE TRANSPORT        !
-     !---------------------------!
-     IF(sec%nb_point .NE. 0)THEN   
-
-        !----------------------------------------------------------------------------------------------------
-        !Compute sign for velocities:
-        !
-        !convention:
-        !   non horizontal section: direction + is toward left hand of section
-        !       horizontal section: direction + is toward north of section
-        !
-        !
-        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
-        !       ----------------      -----------------     ---------------             --------------
-        !
-        !   isgnv=1         direction +      
-        !  ______         _____             ______                                                   
-        !        |           //|            |                  |                         direction +   
-        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
-        !        |_______  //         ______|    \\            | ---\                        |
-        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
-        !               |             |          __\\|         |                    
-        !               |             |     direction +        |                      isgnv=1                                 
-        !                                                      
-        !----------------------------------------------------------------------------------------------------
-        isgnu = 1
-        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
-        ELSE                                ; isgnv =  1
-        ENDIF
-        IF( sec%slopeSection .GE. 9999. )     isgnv =  1
-
-        IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
-
-        !--------------------------------------!
-        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
-        !--------------------------------------!
-        DO jseg=1,MAX(sec%nb_point-1,0)
+      !!--------------------------------------------------------
+      !
+      IF( ld_debug )WRITE(numout,*)'      Compute transport'
+
+      !---------------------------!
+      !  COMPUTE TRANSPORT        !
+      !---------------------------!
+      IF(sec%nb_point .NE. 0)THEN   
+
+         !----------------------------------------------------------------------------------------------------
+         !Compute sign for velocities:
+         !
+         !convention:
+         !   non horizontal section: direction + is toward left hand of section
+         !       horizontal section: direction + is toward north of section
+         !
+         !
+         !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
+         !       ----------------      -----------------     ---------------             --------------
+         !
+         !   isgnv=1         direction +      
+         !  ______         _____             ______                                                   
+         !        |           //|            |                  |                         direction +   
+         !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
+         !        |_______  //         ______|    \\            | ---\                        |
+         !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
+         !               |             |          __\\|         |                    
+         !               |             |     direction +        |                      isgnv=1                                 
+         !                                                      
+         !----------------------------------------------------------------------------------------------------
+         isgnu = 1
+         IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
+         ELSE                                ; isgnv =  1
+         ENDIF
+         IF( sec%slopeSection .GE. 9999. )     isgnv =  1
+
+         IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
+
+         !--------------------------------------!
+         ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
+         !--------------------------------------!
+         DO jseg=1,MAX(sec%nb_point-1,0)
               
-           !-------------------------------------------------------------------------------------------
-           ! Select the appropriate coordinate for computing the velocity of the segment
-           !
-           !                      CASE(0)                                    Case (2)
-           !                      -------                                    --------
-           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)      
-           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
-           !                                                                            |
-           !                                                                            |
-           !                                                                            |
-           !                      Case (3)                                            U(i,j)
-           !                      --------                                              |
-           !                                                                            |
-           !  listPoint(jseg+1) F(i,j+1)                                                |
-           !                        |                                                   |
-           !                        |                                                   |
-           !                        |                                 listPoint(jseg+1) F(i,j-1)
-           !                        |                                            
-           !                        |                                            
-           !                     U(i,j+1)                                            
-           !                        |                                       Case(1)     
-           !                        |                                       ------      
-           !                        |                                            
-           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
-           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
-           ! listPoint(jseg)     F(i,j)
-           ! 
-           !-------------------------------------------------------------------------------------------
-
-           SELECT CASE( sec%direction(jseg) )
-           CASE(0)  ;   k = sec%listPoint(jseg)
-           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
-           CASE(2)  ;   k = sec%listPoint(jseg)
-           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
-           END SELECT
-
-           !---------------------------| 
-           !     LOOP ON THE LEVEL     | 
-           !---------------------------| 
-           !Sum of the transport on the vertical  
-           DO jk=1,mbathy(k%I,k%J) 
- 
-              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 
-              SELECT CASE( sec%direction(jseg) ) 
-              CASE(0,1) 
-                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
-                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
-                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
-                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
-                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
-              CASE(2,3) 
-                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
-                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
-                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
-                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
-                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
-              END SELECT 
- 
-              zfsdep= fsdept(k%I,k%J,jk) 
+            !-------------------------------------------------------------------------------------------
+            ! Select the appropriate coordinate for computing the velocity of the segment
+            !
+            !                      CASE(0)                                    Case (2)
+            !                      -------                                    --------
+            !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)      
+            !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
+            !                                                                            |
+            !                                                                            |
+            !                                                                            |
+            !                      Case (3)                                            U(i,j)
+            !                      --------                                              |
+            !                                                                            |
+            !  listPoint(jseg+1) F(i,j+1)                                                |
+            !                        |                                                   |
+            !                        |                                                   |
+            !                        |                                 listPoint(jseg+1) F(i,j-1)
+            !                        |                                            
+            !                        |                                            
+            !                     U(i,j+1)                                            
+            !                        |                                       Case(1)     
+            !                        |                                       ------      
+            !                        |                                            
+            !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
+            !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
+            ! listPoint(jseg)     F(i,j)
+            ! 
+            !-------------------------------------------------------------------------------------------
+
+            SELECT CASE( sec%direction(jseg) )
+            CASE(0)   ;    k = sec%listPoint(jseg)
+            CASE(1)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
+            CASE(2)   ;    k = sec%listPoint(jseg)
+            CASE(3)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
+            END SELECT
+
+            !---------------------------| 
+            !     LOOP ON THE LEVEL     | 
+            !---------------------------| 
+            DO jk = 1, mbathy(k%I,k%J)            !Sum of the transport on the vertical
+            !           ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 
+            SELECT CASE( sec%direction(jseg) )
+               CASE(0,1) 
+                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
+                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
+                  zrhop = interp(k%I,k%J,jk,'V',rhop) 
+                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
+                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
+               CASE(2,3) 
+                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
+                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
+                  zrhop = interp(k%I,k%J,jk,'U',rhop) 
+                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
+                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
+               END SELECT 
+               !
+               zdep= gdept_n(k%I,k%J,jk) 
   
-              !compute velocity with the correct direction 
-              SELECT CASE( sec%direction(jseg) ) 
-              CASE(0,1)   
-                 zumid=0. 
-                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
-              CASE(2,3) 
-                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
-                 zvmid=0. 
-              END SELECT 
- 
-              !zTnorm=transport through one cell; 
-              !velocity* cell's length * cell's thickness 
-              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
-                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
-
-#if ! defined key_vvl
-              !add transport due to free surface 
-              IF( jk==1 )THEN 
-                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
-                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
-              ENDIF 
-#endif
+               SELECT CASE( sec%direction(jseg) )                !compute velocity with the correct direction 
+               CASE(0,1)   
+                  zumid=0._wp
+                  zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
+               CASE(2,3) 
+                  zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
+                  zvmid=0._wp
+               END SELECT 
+ 
+               !zTnorm=transport through one cell; 
+               !velocity* cell's length * cell's thickness 
+               zTnorm = zumid*e2u(k%I,k%J) * e3u_n(k%I,k%J,jk)     & 
+                  &   + zvmid*e1v(k%I,k%J) * e3v_n(k%I,k%J,jk) 
+
+!!gm  THIS is WRONG  no transport due to ssh in linear free surface case !!!!!
+               IF( ln_linssh ) THEN              !add transport due to free surface 
+                  IF( jk==1 ) THEN 
+                     zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk)   & 
+                        &            + zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
+                  ENDIF 
+               ENDIF
+!!gm end
               !COMPUTE TRANSPORT  
  
               transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 
   
-              IF ( sec%llstrpond ) THEN 
+              IF( sec%llstrpond ) THEN 
                  transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp
                  transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zsn * zrhop * 0.001
               ENDIF
    
-           ENDDO !end of loop on the level
+           END DO !end of loop on the level
 
 #if defined key_lim2 || defined key_lim3
@@ -797,5 +786,5 @@
                  transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
-              ENDDO
+              END DO
 #endif
    
@@ -803,10 +792,11 @@
 #endif
  
-        ENDDO !end of loop on the segment
+        END DO !end of loop on the segment
 
      ENDIF !end of sec%nb_point =0 case
      !
   END SUBROUTINE transport
-  
+
+
   SUBROUTINE dia_dct_sum(sec,jsec) 
      !!------------------------------------------------------------- 
@@ -828,5 +818,4 @@
      !! 
      !!------------------------------------------------------------- 
-     !! * arguments 
      TYPE(SECTION),INTENT(INOUT) :: sec 
      INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
@@ -834,5 +823,5 @@
      TYPE(POINT_SECTION) :: k 
      INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes  
-     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
+     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/ssh/potential density /depth at u/v point 
      !!------------------------------------------------------------- 
  
@@ -903,5 +892,5 @@
               END SELECT 
  
-              zfsdep= fsdept(k%I,k%J,jk) 
+              zdep= gdept_n(k%I,k%J,jk) 
   
               !------------------------------- 
@@ -932,6 +921,6 @@
                     ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
  
-                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
-                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
+                    ((( zdep .GE. sec%zlay(jclass)) .AND.                & 
+                    (   zdep .LE. sec%zlay(jclass+1))) .OR.              & 
                     ( sec%zlay(jclass) .EQ. 99. ))                         & 
                                                                    ))   THEN 
@@ -991,5 +980,5 @@
 #endif 
   
-        ENDDO !end of loop on the segment 
+        END DO !end of loop on the segment 
  
      ELSE  !if sec%nb_point =0 
@@ -1000,5 +989,6 @@
  
   END SUBROUTINE dia_dct_sum 
-  
+
+
   SUBROUTINE dia_dct_wri(kt,ksec,sec)
      !!-------------------------------------------------------------
@@ -1138,13 +1128,15 @@
                               sec%transport(9,1),sec%transport(10,1), &
                               sec%transport(9,1)+sec%transport(10,1) 
-     ENDIF
+      ENDIF
                                               
-118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
-119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
-
-     CALL wrk_dealloc(nb_type_class , zsumclasses )  
-  END SUBROUTINE dia_dct_wri
-
-  FUNCTION interp(ki, kj, kk, cd_point, ptab)
+118   FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
+119   FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
+
+      CALL wrk_dealloc(nb_type_class , zsumclasses )  
+      !
+   END SUBROUTINE dia_dct_wri
+
+
+   FUNCTION interp(ki, kj, kk, cd_point, ptab)
   !!----------------------------------------------------------------------
   !!
@@ -1214,5 +1206,5 @@
   !*local declations
   INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
-  REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real
+  REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu            ! local real
   REAL(wp):: zet1, zet2                                        ! weight for interpolation 
   REAL(wp):: zdep1,zdep2                                       ! differences of depth
@@ -1241,7 +1233,7 @@
   IF( ln_sco )THEN   ! s-coordinate case
 
-     zdepu = ( fsdept(ii1,ij1,kk) +  fsdept(ii2,ij2,kk) ) /2 
-     zdep1 = fsdept(ii1,ij1,kk) - zdepu
-     zdep2 = fsdept(ii2,ij2,kk) - zdepu
+     zdepu = ( gdept_n(ii1,ij1,kk) +  gdept_n(ii2,ij2,kk) ) * 0.5_wp 
+     zdep1 = gdept_n(ii1,ij1,kk) - zdepu
+     zdep2 = gdept_n(ii2,ij2,kk) - zdepu
 
      ! weights
@@ -1255,17 +1247,7 @@
   ELSE       ! full step or partial step case 
 
-#if defined key_vvl
-
-     ze3t  = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk) 
-     zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk)
-     zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk)
-
-#else
-
-     ze3t  = fse3t(ii2,ij2,kk)   - fse3t(ii1,ij1,kk) 
-     zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk)
-     zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk)
-
-#endif
+     ze3t  = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk) 
+     zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
+     zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
 
      IF(kk .NE. 1)THEN
@@ -1288,7 +1270,6 @@
 
   ENDIF
-
-
-  END FUNCTION interp
+      !
+   END FUNCTION interp
 
 #else
@@ -1311,3 +1292,4 @@
 #endif
 
+   !!======================================================================
 END MODULE diadct
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90	(revision 5930)
+++ 	(revision )
@@ -1,229 +1,0 @@
-MODULE diadimg
-   !!======================================================================
-   !!                     ***  MODULE  diadimg  ***
-   !! Ocean diagnostics :  write ocean output files in dimg direct access format (mpp)
-   !!=====================================================================
-# if defined key_dimgout
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers 
-   USE dom_oce         ! ocean space and time domain
-   USE in_out_manager  ! I/O manager
-   USE daymod          ! calendar
-   USE lib_mpp
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC dia_wri_dimg            ! called by trd_mld (eg)
-   PUBLIC dia_wri_dimg_alloc      ! called by nemo_alloc in nemogcm.F90
-
-
-   !! These workspace arrays are inside the module so that we can make them
-   !! allocatable in a clean way. Not done in wrk_nemo because these are of KIND(sp).
-   REAL(sp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: z42d    ! 2d temporary workspace (sp)
-   REAL(sp), ALLOCATABLE, SAVE, DIMENSION(:)   :: z4dep   ! vertical level (sp)
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   FUNCTION dia_wri_dimg_alloc()
-      !!---------------------------------------------------------------------
-      !!        *** ROUTINE dia_wri_dimg_alloc ***
-      !!
-      !!---------------------------------------------------------------------
-      INTEGER :: dia_wri_dimg_alloc   ! return value
-      !!---------------------------------------------------------------------
-      !
-      IF( .NOT. ALLOCATED( z42d ) )THEN
-
-         ALLOCATE( z42d(jpi,jpj), z4dep(jpk), STAT=dia_wri_dimg_alloc )
-
-         IF( lk_mpp                  )   CALL mpp_sum ( dia_wri_dimg_alloc )
-         IF( dia_wri_dimg_alloc /= 0 )   CALL ctl_warn('dia_wri_dimg_alloc: allocation of array failed.')
-
-      ELSE
-
-         dia_wri_dimg_alloc = 0
-
-      ENDIF
-      !
-  END FUNCTION dia_wri_dimg_alloc
-
-
-  SUBROUTINE dia_wri_dimg( cd_name, cd_text, ptab, klev, cd_type , ksubi )
-    !!-------------------------------------------------------------------------
-    !!        *** ROUTINE dia_wri_dimg ***
-    !!
-    !! ** Purpose :   write ptab in the dimg file cd_name, with comment cd_text.
-    !!              ptab has klev x 2D fields
-    !!
-    !! ** Action :   Define header variables from the config parameters
-    !!               Open the dimg file on unit inum = 14 ( IEEE I4R4 file )
-    !!               Write header on record 1
-    !!               Write ptab on the following klev records
-    !!
-    !! History :  2003-12 (J.M. Molines ) : Original. Replace ctl_opn, writn2d
-    !!---------------------------------------------------------------------------
-    CHARACTER(len=*),INTENT(in) ::   &
-         &                            cd_name,  &  ! dimg file name
-         &                            cd_text      ! comment to write on record #1
-    INTEGER, INTENT(in) ::            klev         ! number of level in ptab to write
-    REAL(wp),INTENT(in), DIMENSION(:,:,:) :: ptab  ! 3D array to write 
-    CHARACTER(LEN=1),INTENT(in) ::    cd_type      ! either 'T', 'W' or '2' , depending on the vertical
-    !                                              ! grid for ptab. 2 stands for 2D file
-    INTEGER, INTENT(in), OPTIONAL, DIMENSION(klev) :: ksubi 
-
-    !! * Local declarations
-    INTEGER :: jk, jn           ! dummy loop indices
-    INTEGER :: irecl4,             &    ! record length in bytes
-         &       inum,             &    ! logical unit (set to 14)
-         &       irec,             &    ! current record to be written
-         &       irecend                ! record number where nclit... are stored
-    REAL(sp)                    :: zdx,zdy,zspval,zwest,ztimm
-    REAL(sp)                    :: zsouth
-
-    CHARACTER(LEN=80) :: clname                ! name of file in case of dimgnnn
-    CHARACTER(LEN=4) :: clver='@!01'           ! dimg string identifier
-    !!---------------------------------------------------------------------------
-
-    !                                      ! allocate dia_wri_dimg array
-    IF( dia_wri_dimg_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_wri_dimg : unable to allocate arrays' )
-
-    !! * Initialisations
-
-    irecl4 = MAX(jpi*jpj*sp , 84+(18+1+jpk)*sp )
-
-    zspval=0.0_sp    ! special values on land
-    !  the 'numerical' grid is described. The geographical one is in a grid file
-    zdx=1._sp
-    zdy=1._sp
-    zsouth=njmpp * 1._sp
-    zwest=nimpp * 1._sp
-    !  time in days since the historical begining of the run (nit000 = 0 ) 
-    ztimm=adatrj
-
-    SELECT CASE ( cd_type )
-
-    CASE ( 'T')
-       z4dep(:)=gdept_1d(:)
-
-    CASE ( 'W' )
-       z4dep(:)=gdepw_1d(:)
-
-    CASE ( '2' )
-       z4dep(1:klev) =(/(jk, jk=1,klev)/)
-
-    CASE ( 'I' )
-       z4dep(1:klev) = ksubi(1:klev)
-
-    CASE DEFAULT
-       IF(lwp) WRITE(numout,*) ' E R R O R : bad cd_type in dia_wri_dimg '
-       STOP 'dia_wri_dimg'
-
-    END SELECT
-
-    IF ( ln_dimgnnn  ) THEN
-     irecl4 = MAX(jpi*jpj*sp , 84+(18+jpk)*sp + 8*jpnij*sp  )
-       WRITE(clname,'(a,a,i3.3)') TRIM(cd_name),'.',narea
-       CALL ctl_opn(inum, clname,'UNKNOWN','UNFORMATTED','DIRECT',irecl4,numout,lwp)
-       WRITE(inum,REC=1 ) clver, cd_text, irecl4, &
-            &     jpi,jpj, klev, 1 , 1 ,            &
-            &     zwest, zsouth, zdx, zdy, zspval,  &
-            &     z4dep(1:klev),                    &
-            &     ztimm,                            &
-            &     narea, jpnij,jpiglo,jpjglo,jpizoom, jpjzoom,    &    ! extension to dimg for mpp output
-            &     nlcit,nlcjt, nldit, nldjt, nleit, nlejt, nimppt, njmppt  !
-
-       !! * Write klev levels
-       IF ( cd_type == 'I' ) THEN
-
-          DO jk = 1, klev
-             irec =1 + jk
-             z42d(:,:) = ptab(:,:,ksubi(jk))
-             WRITE(inum,REC=irec)  z42d(:,:)
-          END DO
-       ELSE
-          DO jk = 1, klev
-             irec =1 + jk
-             z42d(:,:) = ptab(:,:,jk)
-             WRITE(inum,REC=irec)  z42d(:,:)
-          END DO
-       ENDIF
-    ELSE
-       clver='@!03'           ! dimg string identifier
-       ! note that version @!02 is optimized with respect to record length.
-       ! The vertical dep variable is reduced to klev instead of klev*jpnij :
-       !   this is OK for jpnij < 181 (jpk=46)
-       ! for more processors, irecl4 get huge and that's why we switch to '@!03':
-       !  In this case we just add an extra integer to the standard dimg structure,
-       !  which is a record number where the arrays nlci etc... starts (1 per record)
-       
-       !! Standard dimgproc (1 file per variable, all procs. write to this file )
-       !! * Open file
-       CALL ctl_opn(inum, cd_name,'UNKNOWN','UNFORMATTED','DIRECT',irecl4,numout,lwp)
-
-       !! * Write header on record #1
-       irecend=1 + klev*jpnij 
-       IF(lwp) WRITE(inum,REC=1 ) clver, cd_text, irecl4, &
-            &     jpi,jpj, klev, 1 , 1 ,            &
-            &     zwest, zsouth, zdx, zdy, zspval,  &
-            &     z4dep(1:klev),       &
-            &     ztimm,                            &
-            &     narea, jpnij,jpiglo,jpjglo,jpizoom, jpjzoom, irecend 
-       IF (lwp ) THEN
-         WRITE(inum,REC=irecend + 1 ) nlcit
-         WRITE(inum,REC=irecend + 2 ) nlcjt
-         WRITE(inum,REC=irecend + 3 ) nldit
-         WRITE(inum,REC=irecend + 4 ) nldjt
-         WRITE(inum,REC=irecend + 5 ) nleit
-         WRITE(inum,REC=irecend + 6 ) nlejt
-         WRITE(inum,REC=irecend + 7 ) nimppt
-         WRITE(inum,REC=irecend + 8 ) njmppt
-       ENDIF
-      !   &    ! extension to dimg for mpp output
-      !   &     nlcit,nlcjt, nldit, nldjt, nleit, nlejt, nimppt, njmppt  !
-
-       !! * Write klev levels
-       IF ( cd_type == 'I' ) THEN
-
-          DO jk = 1, klev
-             irec =1 + klev * (narea -1) + jk
-             z42d(:,:) = ptab(:,:,ksubi(jk))
-             WRITE(inum,REC=irec)  z42d(:,:)
-          END DO
-       ELSE
-          DO jk = 1, klev
-             irec =1 + klev * (narea -1) + jk
-             z42d(:,:) = ptab(:,:,jk)
-             WRITE(inum,REC=irec)  z42d(:,:)
-          END DO
-       ENDIF
-    ENDIF
-
-    !! * Close the file
-    CLOSE(inum)
-
-  END SUBROUTINE dia_wri_dimg
-
-#  else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                       Empty module
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dia_wri_dimg( cd_name, cd_exper, ptab, klev, cd_type )
-      REAL, DIMENSION(:,:,:) :: ptab
-      INTEGER :: klev
-      CHARACTER(LEN=80) :: cd_name, cd_exper,cd_type
-      WRITE(*,*) ' This print must never occur ', cd_name, cd_exper,ptab, klev, cd_type
-      WRITE(*,*) ' this routine is here just for compilation '
-   END SUBROUTINE dia_wri_dimg
-# endif
-   !!======================================================================
-END MODULE diadimg
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 5930)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 6140)
@@ -33,5 +33,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -40,5 +39,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -80,5 +78,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zwei  = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+                  zwei  = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
                   a_salb = a_salb + ( tsb(ji,jj,jk,jp_sal) - zsm0 ) * zwei
                END DO
@@ -106,5 +104,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zwei  = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+                  zwei  = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
                   a_saln = a_saln + ( tsn(ji,jj,jk,jp_sal) - zsm0 ) * zwei
                   zvol  = zvol  + zwei
@@ -116,5 +114,5 @@
          
          ! Conversion in m3
-         a_fwf    = a_fwf * rdttra(1) * 1.e-3 
+         a_fwf    = a_fwf * rdt * 1.e-3 
          
          ! fwf correction to bring back the mean ssh to zero
@@ -186,5 +184,5 @@
                   zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
                   zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
-                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
+                  zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
 
                   IF( un(ji,jj,jk) > 0.e0 ) THEN 
@@ -238,5 +236,5 @@
                   zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
                   zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
-                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
+                  zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
                   
                   IF( un(ji,jj,jk) > 0.e0 ) THEN 
@@ -290,5 +288,5 @@
                   zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
                   zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
-                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
+                  zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
                   
                   IF( un(ji,jj,jk) > 0.e0 ) THEN 
@@ -342,5 +340,5 @@
                   zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
                   zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
-                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
+                  zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
                   
                   IF( un(ji,jj,jk) > 0.e0 ) THEN 
@@ -404,5 +402,5 @@
          WRITE(inum,*)
          WRITE(inum,*)    'Net freshwater budget '
-         WRITE(inum,9010) '  fwf    = ',a_fwf,   ' m3 =', a_fwf   /(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
+         WRITE(inum,9010) '  fwf    = ',a_fwf,   ' m3 =', a_fwf   /(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
          WRITE(inum,*)
          WRITE(inum,9010) '  zarea =',zarea
@@ -460,5 +458,5 @@
       ENDIF
 
-      IF( nn_timing == 1 )   CALL timing_start('dia_fwb')
+      IF( nn_timing == 1 )   CALL timing_stop('dia_fwb')
 
  9005 FORMAT(1X,A,ES24.16)
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 5930)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 6140)
@@ -21,5 +21,4 @@
    USE ioipsl          ! NetCDF IPSL library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE diadimg         ! To write dimg
    USE timing          ! preformance summary
    USE wrk_nemo        ! working arrays
@@ -135,5 +134,5 @@
       DO jk=1,nb_ana
        DO ji=1,jpmax_harmo
-          IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN
+          IF (TRIM(tname(jk)) ==  Wave(ji)%cname_tide) THEN
              name(jk) = ji
              EXIT
@@ -194,7 +193,7 @@
                   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)        
-                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*umask_i(ji,jj)
-                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*vmask_i(ji,jj)
+                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)        
+                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)
+                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)
                   END DO
                END DO
@@ -324,6 +323,6 @@
                X1= ana_amp(ji,jj,jh,1)
                X2=-ana_amp(ji,jj,jh,2)
-               out_u(ji,jj,       jh) = X1 * umask_i(ji,jj)
-               out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)
+               out_u(ji,jj,       jh) = X1 * ssumask(ji,jj)
+               out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj)
             ENDDO
          ENDDO
@@ -358,6 +357,6 @@
                X1=ana_amp(ji,jj,jh,1)
                X2=-ana_amp(ji,jj,jh,2)
-               out_v(ji,jj,       jh)=X1 * vmask_i(ji,jj)
-               out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)
+               out_v(ji,jj,       jh)=X1 * ssvmask(ji,jj)
+               out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj)
             END DO
          END DO
@@ -384,17 +383,6 @@
       !!----------------------------------------------------------------------
 
-#if defined key_dimgout
-      cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'
-      cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'
-      cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'
-#endif
-
       IF(lwp) WRITE(numout,*) '  '
       IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
-#if defined key_dimgout
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  Output files: ', TRIM(cdfile_name_T)
-      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_U)
-      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_V)
-#endif
       IF(lwp) WRITE(numout,*) '  '
 
@@ -402,39 +390,24 @@
       !/////////////
       !
-#if defined key_dimgout
-      cltext='Elevation amplitude and phase'
-      CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')
-#else
       DO jh = 1, nb_ana
       CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
       CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
       END DO
-#endif
 
       ! B) ubar
       !/////////
       !
-#if defined key_dimgout
-      cltext='ubar amplitude and phase'
-      CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')
-#else
       DO jh = 1, nb_ana
       CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
       CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) )
       END DO
-#endif
 
       ! C) vbar
       !/////////
       !
-#if defined key_dimgout
-      cltext='vbar amplitude and phase'
-      CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')
-#else
       DO jh = 1, nb_ana
          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
@@ -488,5 +461,5 @@
             DO jj_sd = ji_sd, ninco
                zval2 = ABS(ztmp3(ji_sd,jj_sd))
-               IF( zval2.GE.zval1 )THEN
+               IF( zval2 >= zval1 )THEN
                   ipivot(ji_sd) = jj_sd
                   zval1         = zval2
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 5930)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 6140)
@@ -46,10 +46,10 @@
    REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends
    !
-   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf          , ssh_ini          !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf 
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf_ini      , 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"
    !!----------------------------------------------------------------------
@@ -100,9 +100,5 @@
       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) * surf(:,:) )
-          z_frc_trd_s = z_frc_trd_s + glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) )
-      ENDIF
-
+      IF( ln_isf    ) z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) )
       ! Add penetrative solar radiation
       IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) )
@@ -110,12 +106,12 @@
       IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * surf(:,:) )
       !
-      IF( .NOT. lk_vvl ) THEN
-         IF ( ln_isfcav ) THEN
+      IF( ln_linssh ) THEN
+         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
+               END DO
+            END DO
          ELSE
             z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem)
@@ -130,5 +126,5 @@
       frc_s = frc_s + z_frc_trd_s * rdt
       !                                          ! Advection flux through fixed surface (z=0)
-      IF( .NOT. lk_vvl ) THEN
+      IF( ln_linssh ) THEN
          frc_wn_t = frc_wn_t + z_wn_trd_t * rdt
          frc_wn_s = frc_wn_s + z_wn_trd_s * rdt
@@ -138,4 +134,5 @@
       ! 2 -  Content variations !
       ! ------------------------ !
+      ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
       zdiff_v2 = 0._wp
       zdiff_hc = 0._wp
@@ -143,9 +140,9 @@
 
       ! volume variation (calculated with ssh)
-      zdiff_v1 = glob_sum( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
+      zdiff_v1 = glob_sum_full( surf(:,:) * sshn(:,:) - surf_ini(:,:) * ssh_ini(:,:) )
 
       ! heat & salt content variation (associated with ssh)
-      IF( .NOT. lk_vvl ) THEN
-         IF ( ln_isfcav ) THEN
+      IF( ln_linssh ) THEN
+         IF( ln_isfcav ) THEN
             DO ji = 1, jpi
                DO jj = 1, jpj
@@ -158,26 +155,26 @@
             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 ) 
+         z_ssh_hc = glob_sum_full( z2d0 ) 
+         z_ssh_sc = glob_sum_full( z2d1 ) 
       ENDIF
 
       DO jk = 1, jpkm1
          ! volume variation (calculated with scale factors)
-         zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk) &
-            &                           * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
+         zdiff_v2 = zdiff_v2 + glob_sum_full( surf(:,:) * tmask(:,:,jk)            &
+            &                           * e3t_n(:,:,jk) - surf_ini(:,:) * e3t_ini(:,:,jk) )
          ! heat content variation
-         zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * tmask(:,:,jk) & 
-            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )
+         zdiff_hc = zdiff_hc + glob_sum_full(  surf(:,:) * tmask(:,:,jk)                                   & 
+            &                           * e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - surf_ini(:,:) * hc_loc_ini(:,:,jk) ) 
          ! salt content variation
-         zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * tmask(:,:,jk)   &
-            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) )
+         zdiff_sc = zdiff_sc + glob_sum_full( surf    (:,:) * tmask(:,:,jk)                           &
+                                        * e3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - surf_ini(:,:) * sc_loc_ini(:,:,jk) )
       ENDDO
 
       ! Substract forcing from heat content, salt content and volume variations
       zdiff_v1 = zdiff_v1 - frc_v
-      IF( lk_vvl )   zdiff_v2 = zdiff_v2 - frc_v
+      IF( .NOT.ln_linssh )   zdiff_v2 = zdiff_v2 - frc_v
       zdiff_hc = zdiff_hc - frc_t
       zdiff_sc = zdiff_sc - frc_s
-      IF( .NOT. lk_vvl ) THEN
+      IF( ln_linssh ) THEN
          zdiff_hc1 = zdiff_hc + z_ssh_hc 
          zdiff_sc1 = zdiff_sc + z_ssh_sc
@@ -191,14 +188,25 @@
       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) )
+         zvol_tot  = zvol_tot + glob_sum_full( surf(:,:) * tmask(:,:,jk) * e3t_n(:,:,jk) )
       END DO
 
 !!gm to be added ?
-!      IF( .NOT. lk_vvl ) THEN            ! fixed volume, add the ssh contribution
+!      IF( ln_linssh ) THEN            ! fixed volume, add the ssh contribution
 !        zvol_tot = zvol_tot + glob_sum( surf(:,:) * sshn(:,:) )
 !      ENDIF
 !!gm end
 
-      IF( lk_vvl ) THEN
+      IF( ln_linssh ) THEN
+        CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content variation (C) 
+        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu)
+        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J) 
+        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9   )              ! Salt content variation (psu*km3)
+        CALL iom_put( 'bgvolssh' , zdiff_v1  * 1.e-9   )              ! volume ssh variation (km3)  
+        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
+        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C) 
+        CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu) 
+        CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot )              ! hc  - error due to free surface (C)
+        CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot )              ! sc  - error due to free surface (psu)
+      ELSE
         CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature variation (C) 
         CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    variation (psu)
@@ -210,21 +218,10 @@
         CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C) 
         CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu) 
-      ELSE
-        CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content variation (C) 
-        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu)
-        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J) 
-        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content variation (psu*km3)
-        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)  
-        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
-        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C) 
-        CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu) 
-        CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot )              ! hc  - error due to free surface (C)
-        CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot )              ! sc  - error due to free surface (psu)
       ENDIF
       !
       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')
       !
@@ -257,13 +254,14 @@
            CALL iom_get( numror, 'frc_t', frc_t )
            CALL iom_get( numror, 'frc_s', frc_s )
-           IF( .NOT. lk_vvl ) THEN
+           IF( ln_linssh ) THEN
               CALL iom_get( numror, 'frc_wn_t', frc_wn_t )
               CALL iom_get( numror, 'frc_wn_s', frc_wn_s )
            ENDIF
+           CALL iom_get( numror, jpdom_autoglo, 'surf_ini', surf_ini ) ! ice sheet coupling
            CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini )
            CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini )
            CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini )
            CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini )
-           IF( .NOT. lk_vvl ) THEN
+           IF( ln_linssh ) THEN
               CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
               CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
@@ -273,14 +271,16 @@
           IF(lwp) WRITE(numout,*) ' dia_hsb at initial state '
           IF(lwp) WRITE(numout,*) '~~~~~~~'
-          ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh
+          surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
+          ssh_ini(:,:) = sshn(:,:)                          ! initial ssh
           DO jk = 1, jpk
-             e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                        ! initial vertical scale factors
-             hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk)   ! initial heat content
-             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content
+             ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
+             e3t_ini   (:,:,jk) = e3t_n(:,:,jk)                      * tmask(:,:,jk)  ! initial vertical scale factors
+             hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial heat content
+             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial salt content
           END DO
           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
+          IF( ln_linssh ) THEN
              IF ( ln_isfcav ) THEN
                 DO ji=1,jpi
@@ -308,13 +308,14 @@
         CALL iom_rstput( kt, nitrst, numrow, 'frc_t'   , frc_t     )
         CALL iom_rstput( kt, nitrst, numrow, 'frc_s'   , frc_s     )
-        IF( .NOT. lk_vvl ) THEN
+        IF( ln_linssh ) THEN
            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t )
            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s )
         ENDIF
+        CALL iom_rstput( kt, nitrst, numrow, 'surf_ini', surf_ini )      ! ice sheet coupling
         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini )
         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini )
         CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini )
         CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
-        IF( .NOT. lk_vvl ) THEN
+        IF( ln_linssh ) THEN
            CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
            CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
@@ -379,11 +380,11 @@
       ! 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 )
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
+         &      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( ln_linssh )   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
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90	(revision 5930)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90	(revision 6140)
@@ -20,4 +20,5 @@
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
+   !
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
@@ -31,5 +32,6 @@
    PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90
 
-   LOGICAL , PUBLIC, PARAMETER          ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag
+   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag
+   
    ! note: following variables should move to local variables once iom_put is always used 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m]
@@ -38,6 +40,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [W]
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -52,5 +52,5 @@
       !!---------------------------------------------------------------------
       !
-      ALLOCATE(hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc)
+      ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc )
       !
       IF( lk_mpp           )   CALL mpp_sum ( dia_hth_alloc )
@@ -108,7 +108,7 @@
       IF( kt == nit000 ) THEN
          !                                      ! allocate dia_hth array
-         IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
-
-         IF(.not. ALLOCATED(ik20))THEN
+         IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' )
+
+         IF(.NOT. ALLOCATED(ik20) ) THEN
             ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), &
                &      zabs2(jpi,jpj),   &
@@ -187,5 +187,5 @@
             DO ji = 1, jpi
                !
-               zzdep = fsdepw(ji,jj,jk)
+               zzdep = gdepw_n(ji,jj,jk)
                zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
                zzdep = zzdep * tmask(ji,jj,1)
@@ -223,5 +223,5 @@
             DO ji = 1, jpi
                !
-               zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)
+               zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1)
                !
                zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
@@ -270,10 +270,10 @@
          DO ji = 1, jpi
             !
-            zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
+            zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
             !
             iid = ik20(ji,jj)
             IF( iid /= 1 ) THEN 
-               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
-                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
+               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
+                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
                   &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
                   &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
@@ -285,6 +285,6 @@
             iid = ik28(ji,jj)
             IF( iid /= 1 ) THEN 
-               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
-                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
+               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
+                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
                   &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   &
                   &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
@@ -311,11 +311,11 @@
       END DO
       ! surface boundary condition
-      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
-      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
+      IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)  
+      ELSE                   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
       ENDIF
       ! integration down to ilevel
       DO jk = 1, ilevel
-         zthick(:,:) = zthick(:,:) + fse3t(:,:,jk)
-         htc3  (:,:) = htc3  (:,:) + fse3t(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
+         zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk)
+         htc3  (:,:) = htc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
       END DO
       ! deepest layer
@@ -323,6 +323,6 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) )  &
-                                                                   * tmask(ji,jj,ilevel+1)
+            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  &
+               &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
          END DO
       END DO
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/dianam.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/dianam.F90	(revision 5930)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/dianam.F90	(revision 6140)
@@ -72,5 +72,5 @@
 
       IF( llfsec .OR. kfreq < 0 ) THEN   ;   inbsec = kfreq                       ! output frequency already in seconds
-      ELSE                               ;   inbsec = kfreq * NINT( rdttra(1) )   ! from time-step to seconds
+      ELSE                               ;   inbsec = kfreq * NINT( rdt )   ! from time-step to seconds
       ENDIF
       iddss = NINT( rday          )                                         ! number of seconds in 1 day
@@ -116,6 +116,6 @@
       ! date of the beginning and the end of the run
 
-      zdrun = rdttra(1) / rday * REAL( nitend - nit000, wp )                ! length of the run in days
-      zjul  = fjulday - rdttra(1) / rday
+      zdrun = rdt / rday * REAL( nitend - nit000, wp )                ! length of the run in days
+      zjul  = fjulday - rdt / rday
       CALL ju2ymds( zjul        , iyear1, imonth1, iday1, zsec1 )           ! year/month/day of the beginning of run
       CALL ju2ymds( zjul + zdrun, iyear2, imonth2, iday2, zsec2 )           ! year/month/day of the end       of run
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 5930)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 6140)
@@ -59,7 +59,5 @@
    REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: p_fval2d
 
-
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -118,11 +116,11 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zsfc = e1t(ji,jj) * fse3t(ji,jj,jk)
+                     zsfc = e1t(ji,jj) * e3t_n(ji,jj,jk)
                      zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
                      zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc
                      zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc
-                  ENDDO
-               ENDDO
-            ENDDO
+                  END DO
+               END DO
+            END DO
             DO jn = 1, nptr
                zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 6140)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 6140)
@@ -0,0 +1,169 @@
+MODULE diatmb 
+   !!======================================================================
+   !!                       ***  MODULE  diaharm  ***
+   !! Harmonic analysis of tidal constituents 
+   !!======================================================================
+   !! History :  3.6  !  2014  (E O'Dea)  Original code
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE wrk_nemo        ! working arrays
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   LOGICAL , PUBLIC ::   ln_diatmb     !: Top Middle and Bottom output 
+   PUBLIC   dia_tmb_init            ! routine called by nemogcm.F90
+   PUBLIC   dia_tmb                 ! routine called by diawri.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_tmb_init 
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_tmb_init  ***
+      !!     
+      !! ** Purpose: Initialization of tmb namelist 
+      !!        
+      !! ** Method : Read namelist
+      !!   History
+      !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_tmb
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER ::   ios                 ! Local integer output status for namelist read
+      !
+      NAMELIST/nam_diatmb/ ln_diatmb
+      !!----------------------------------------------------------------------
+      !
+      REWIND ( numnam_ref )              ! Read Namelist nam_diatmb in reference namelist : TMB diagnostics
+      READ   ( numnam_ref, nam_diatmb, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in reference namelist', lwp )
+ 
+      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics
+      READ  ( numnam_cfg, nam_diatmb, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_diatmb )
+
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_tmb_init : Output Top, Middle, Bottom Diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) 'Namelist nam_diatmb : set tmb outputs '
+         WRITE(numout,*) 'Switch for TMB diagnostics (T) or not (F)  ln_diatmb  = ', ln_diatmb
+      ENDIF
+
+   END SUBROUTINE dia_tmb_init
+
+   SUBROUTINE dia_calctmb( pinfield,pouttmb )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_tmb  ***
+      !!                   
+      !! ** Purpose :    Find the Top, Mid and Bottom fields of water Column
+      !!
+      !! ** Method  :   
+      !!      use mbathy to find surface, mid and bottom of model levels
+      !!
+      !! History :
+      !!   3.6  !  08-14  (E. O'Dea) Routine based on dia_wri_foam
+      !!----------------------------------------------------------------------
+      !! * Modules used
+
+      ! Routine to map 3d field to top, middle, bottom
+      IMPLICIT NONE
+
+
+      ! Routine arguments
+      REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN   ) :: pinfield    ! Input 3d field and mask
+      REAL(wp), DIMENSION(jpi, jpj, 3  ), INTENT(  OUT) :: pouttmb     ! Output top, middle, bottom
+
+
+
+      ! Local variables
+      INTEGER :: ji,jj,jk  ! Dummy loop indices
+
+      ! Local Real
+      REAL(wp)                         ::   zmdi  !  set masked values
+
+      zmdi=1.e+20 !missing data indicator for masking
+
+      ! Calculate top
+      pouttmb(:,:,1) = pinfield(:,:,1)*tmask(:,:,1)  + zmdi*(1.0-tmask(:,:,1))
+
+      ! Calculate middle
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            jk              = max(1,mbathy(ji,jj)/2)
+            pouttmb(ji,jj,2) = pinfield(ji,jj,jk)*tmask(ji,jj,jk)  + zmdi*(1.0-tmask(ji,jj,jk))
+         END DO
+      END DO
+
+      ! Calculate bottom
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            jk              = max(1,mbathy(ji,jj) - 1)
+            pouttmb(ji,jj,3) = pinfield(ji,jj,jk)*tmask(ji,jj,jk)  + zmdi*(1.0-tmask(ji,jj,jk))
+         END DO
+      END DO
+
+   END SUBROUTINE dia_calctmb
+
+
+
+   SUBROUTINE dia_tmb
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_tmb  ***
+      !! ** Purpose :   Write diagnostics for Top, Mid and Bottom of water Column
+      !!
+      !! ** Method  :   
+      !!      use mbathy to find surface, mid and bottom of model levels
+      !!      calls calctmb to retrieve TMB values before sending to iom_put
+      !!
+      !! History :
+      !!   3.6  !  08-14  (E. O'Dea) 
+      !!         
+      !!--------------------------------------------------------------------
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb    ! temporary workspace 
+      REAL(wp)                         ::   zmdi      ! set masked values
+
+      zmdi=1.e+20 !missing data indicator for maskin
+
+      IF (ln_diatmb) THEN
+         CALL wrk_alloc( jpi , jpj, 3 , zwtmb )
+         CALL dia_calctmb(  tsn(:,:,:,jp_tem),zwtmb )
+         !ssh already output but here we output it masked
+         CALL iom_put( "sshnmasked" , sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )   ! tmb Temperature
+         CALL iom_put( "top_temp" , zwtmb(:,:,1) )    ! tmb Temperature
+         CALL iom_put( "mid_temp" , zwtmb(:,:,2) )    ! tmb Temperature
+         CALL iom_put( "bot_temp" , zwtmb(:,:,3) )    ! tmb Temperature
+!         CALL iom_put( "sotrefml" , hmld_tref(:,:) )    ! "T criterion Mixed Layer Depth
+
+         CALL dia_calctmb(  tsn(:,:,:,jp_sal),zwtmb )
+         CALL iom_put( "top_sal" , zwtmb(:,:,1) )    ! tmb Salinity 
+         CALL iom_put( "mid_sal" , zwtmb(:,:,2) )    ! tmb Salinity
+         CALL iom_put( "bot_sal" , zwtmb(:,:,3) )    ! tmb Salinity
+
+         CALL dia_calctmb(  un(:,:,:),zwtmb )
+         CALL iom_put( "top_u" , zwtmb(:,:,1) )    ! tmb  U Velocity
+         CALL iom_put( "mid_u" , zwtmb(:,:,2) )    ! tmb  U Velocity
+         CALL iom_put( "bot_u" , zwtmb(:,:,3) )    ! tmb  U Velocity
+!Called in  dynspg_ts.F90        CALL iom_put( "baro_u" , un_b )    ! Barotropic  U Velocity
+
+         CALL dia_calctmb(  vn(:,:,:),zwtmb )
+         CALL iom_put( "top_v" , zwtmb(:,:,1) )    ! tmb  V Velocity
+         CALL iom_put( "mid_v" , zwtmb(:,:,2) )    ! tmb  V Velocity
+         CALL iom_put( "bot_v" , zwtmb(:,:,3) )    ! tmb  V Velocity
+!Called in  dynspg_ts.F90       CALL iom_put( "baro_v" , vn_b )    ! Barotropic  V Velocity
+      ELSE
+         CALL ctl_warn('dia_tmb: tmb diagnostic is set to false you should not have seen this')
+      ENDIF
+
+   END SUBROUTINE dia_tmb
+   !!======================================================================
+END MODULE diatmb
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 5930)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6140)
@@ -30,4 +30,5 @@
    USE zdf_oce         ! ocean vertical physics
    USE ldftra          ! lateral physics: eddy diffusivity coef.
+   USE ldfdyn          ! lateral physics: eddy viscosity   coef.
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
@@ -40,7 +41,9 @@
    USE zdfddm          ! vertical  physics: double diffusion
    USE diahth          ! thermocline diagnostics
+   !
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
-   USE diadimg         ! dimg direct access file format output
+   USE diatmb          ! Top,middle,bottom output
+   USE dia25h          ! 25h Mean output
    USE iom
    USE ioipsl
@@ -53,4 +56,6 @@
    USE lib_mpp         ! MPP library
    USE timing          ! preformance summary
+   USE diurnal_bulk    ! diurnal warm layer
+   USE cool_skin       ! Cool skin
    USE wrk_nemo        ! working array
 
@@ -74,5 +79,4 @@
    !! * Substitutions
 #  include "zdfddm_substitute.h90"
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -97,15 +101,8 @@
   END FUNCTION dia_wri_alloc
 
-#if defined key_dimgout
-   !!----------------------------------------------------------------------
-   !!   'key_dimgout'                                      DIMG output file
-   !!----------------------------------------------------------------------
-#   include "diawri_dimg.h90"
-
-#else
    !!----------------------------------------------------------------------
    !!   Default option                                   NetCDF output file
    !!----------------------------------------------------------------------
-# if defined key_iomput
+#if defined key_iomput
    !!----------------------------------------------------------------------
    !!   'key_iomput'                                        use IOM library
@@ -143,9 +140,9 @@
       ENDIF
 
-      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(:,:,:) )
+      IF( ln_linssh ) THEN
+         CALL iom_put( "e3t" , e3t_n(:,:,:) )
+         CALL iom_put( "e3u" , e3u_n(:,:,:) )
+         CALL iom_put( "e3v" , e3v_n(:,:,:) )
+         CALL iom_put( "e3w" , e3w_n(:,:,:) )
       ENDIF
 
@@ -204,10 +201,4 @@
          CALL iom_put( "sbu", z2d )                ! bottom i-current
       ENDIF
-
-      IF ( ln_dynspg_ts ) THEN
-         CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
-      ELSE
-         CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current
-      ENDIF
       
       CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
@@ -223,10 +214,4 @@
       ENDIF
 
-      IF ( ln_dynspg_ts ) THEN
-         CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current
-      ELSE
-         CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current
-      ENDIF
-
       CALL iom_put( "woce", wn )                   ! vertical velocity
       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
@@ -266,5 +251,5 @@
             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)
+                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -278,5 +263,5 @@
             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)
+                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -290,11 +275,11 @@
             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   = 1._wp / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
+                  zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)    &
+                     &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * e3u_n(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) )  &
+                  zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)    &
+                     &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) )  &
                      &          *  zztmp 
                   !
@@ -311,5 +296,5 @@
          z3d(:,:,jpk) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
+            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
          END DO
          CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
@@ -346,5 +331,5 @@
          z3d(:,:,jpk) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
+            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
          END DO
          CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
@@ -380,4 +365,13 @@
       CALL wrk_dealloc( jpi , jpj, jpk , z3d )
       !
+      ! If we want tmb values 
+
+      IF (ln_diatmb) THEN
+         CALL dia_tmb 
+      ENDIF 
+      IF (ln_dia25h) THEN
+         CALL dia_25h( kt )
+      ENDIF 
+
       IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
       !
@@ -410,5 +404,5 @@
       INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
       INTEGER  ::   jn, ierror                               ! local integers
-      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
+      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
       !
       REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
@@ -418,6 +412,6 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_wri')
       !
-                     CALL wrk_alloc( jpi,jpj      , zw2d )
-      IF( lk_vvl )   CALL wrk_alloc( jpi,jpj,jpk  , zw3d )
+                             CALL wrk_alloc( jpi,jpj      , zw2d )
+      IF( .NOT.ln_linssh )   CALL wrk_alloc( jpi,jpj,jpk  , zw3d )
       !
       ! Output the initial state and forcings
@@ -435,16 +429,14 @@
 
       ! Define frequency of output and means
-      zdt = rdt
-      IF( nacc == 1 ) zdt = rdtmin
       clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
 #if defined key_diainstant
-      zsto = nwrite * zdt
+      zsto = nwrite * rdt
       clop = "inst("//TRIM(clop)//")"
 #else
-      zsto=zdt
+      zsto=rdt
       clop = "ave("//TRIM(clop)//")"
 #endif
-      zout = nwrite * zdt
-      zmax = ( nitend - nit000 + 1 ) * zdt
+      zout = nwrite * rdt
+      zmax = ( nitend - nit000 + 1 ) * rdt
 
       ! Define indices of the horizontal output zoom and vertical limit storage
@@ -488,5 +480,5 @@
          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
             &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
-            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
+            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
          CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
             &           "m", ipk, gdept_1d, nz_T, "down" )
@@ -524,5 +516,5 @@
          CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
             &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
-            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
+            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
          CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
             &           "m", ipk, gdept_1d, nz_U, "down" )
@@ -537,5 +529,5 @@
          CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
             &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
-            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
+            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
          CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
             &          "m", ipk, gdept_1d, nz_V, "down" )
@@ -550,5 +542,5 @@
          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
             &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
-            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
+            &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
          CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
             &          "m", ipk, gdepw_1d, nz_W, "down" )
@@ -562,5 +554,5 @@
          CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
             &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-         IF(  lk_vvl  ) THEN
+         IF(  .NOT.ln_linssh  ) THEN
             CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
@@ -583,5 +575,5 @@
          CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         IF(  .NOT. lk_vvl  ) THEN
+         IF(  ln_linssh  ) THEN
             CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
             &                                                                  , "KgC/m2/s",  &  ! sosst_cd
@@ -729,9 +721,9 @@
       ENDIF
 
-      IF( lk_vvl ) THEN
-         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
-         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
-         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
-         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
+      IF( .NOT.ln_linssh ) THEN
+         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
+         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
+         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
+         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
       ELSE
          CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
@@ -740,8 +732,8 @@
          CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
       ENDIF
-      IF( lk_vvl ) THEN
-         zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
-         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
-         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
+      IF( .NOT.ln_linssh ) THEN
+         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
+         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
+         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
          CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
       ENDIF
@@ -752,5 +744,5 @@
                                                                                   ! (includes virtual salt flux beneath ice 
                                                                                   ! in linear free surface case)
-      IF( .NOT. lk_vvl ) THEN
+      IF( ln_linssh ) THEN
          zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
          CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
@@ -837,12 +829,10 @@
       ENDIF
       !
-                     CALL wrk_dealloc( jpi , jpj        , zw2d )
-      IF( lk_vvl )   CALL wrk_dealloc( jpi , jpj , jpk  , zw3d )
+                             CALL wrk_dealloc( jpi , jpj        , zw2d )
+      IF( .NOT.ln_linssh )   CALL wrk_dealloc( jpi , jpj , jpk  , zw3d )
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
       !
    END SUBROUTINE dia_wri
-# endif
-
 #endif
 
@@ -867,5 +857,5 @@
       INTEGER  ::   id_i , nz_i, nh_i       
       INTEGER, DIMENSION(1) ::   idex             ! local workspace
-      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
+      REAL(wp) ::   zsto, zout, zmax, zjulian
       !!----------------------------------------------------------------------
       ! 
@@ -876,9 +866,8 @@
       clname = cdfile_name
       IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
-      zdt  = rdt
       zsto = rdt
       clop = "inst(x)"           ! no use of the mask value (require less cpu time)
       zout = rdt
-      zmax = ( nitend - nit000 + 1 ) * zdt
+      zmax = ( nitend - nit000 + 1 ) * rdt
 
       IF(lwp) WRITE(numout,*)
@@ -895,5 +884,5 @@
       zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
       CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
-          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
+          1, jpi, 1, jpj, nit000-1, zjulian, rdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
       CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
           "m", jpk, gdept_1d, nz_i, "down")
@@ -913,4 +902,14 @@
       CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
          &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
+         !
+      CALL histdef( id_i, "ahtu"    , "u-eddy diffusivity"    , "m2/s"    ,   &   ! zonal current
+         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
+      CALL histdef( id_i, "ahtv"    , "v-eddy diffusivity"    , "m2/s"    ,   &   ! meridonal current
+         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
+      CALL histdef( id_i, "ahmt"    , "t-eddy viscosity"      , "m2/s"    ,   &   ! zonal current
+         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
+      CALL histdef( id_i, "ahmf"    , "f-eddy viscosity"      , "m2/s"    ,   &   ! meridonal current
+         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
+         !
       CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater 
          &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -925,6 +924,8 @@
       CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
          &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-      IF( lk_vvl ) THEN
-         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      ,   &   ! t-point depth
+      IF( .NOT.ln_linssh ) THEN
+         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      , &   ! t-point depth
+            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
+         CALL histdef( id_i, "vovvle3t", "T point thickness"     , "m"      , &   ! t-point depth
             &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
       ENDIF
@@ -952,4 +953,10 @@
       CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
       CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
+      !
+      CALL histwrite( id_i, "ahtu"    , kt, ahtu             , jpi*jpj*jpk, idex )    ! aht at u-point
+      CALL histwrite( id_i, "ahtv"    , kt, ahtv             , jpi*jpj*jpk, idex )    !  -  at v-point
+      CALL histwrite( id_i, "ahmt"    , kt, ahmt             , jpi*jpj*jpk, idex )    ! ahm at t-point
+      CALL histwrite( id_i, "ahmf"    , kt, ahmf             , jpi*jpj*jpk, idex )    !  -  at f-point
+      !
       CALL histwrite( id_i, "sowaflup", kt, emp-rnf          , jpi*jpj    , idex )    ! freshwater budget
       CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
@@ -959,8 +966,12 @@
       CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
 
+      IF(  .NOT.ln_linssh  ) THEN             
+         CALL histwrite( id_i, "vovvldep", kt, gdept_n(:,:,:), jpi*jpj*jpk, idex )!  T-cell depth 
+         CALL histwrite( id_i, "vovvle3t", kt, e3t_n (:,:,:) , jpi*jpj*jpk, idex )!  T-cell thickness  
+      END IF 
       ! 3. Close the file
       ! -----------------
       CALL histclo( id_i )
-#if ! defined key_iomput && ! defined key_dimgout
+#if ! defined key_iomput
       IF( ninist /= 1  ) THEN
          CALL histclo( nid_T )
@@ -972,4 +983,5 @@
       ! 
    END SUBROUTINE dia_wri_state
+
    !!======================================================================
 END MODULE diawri
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri_dimg.h90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri_dimg.h90	(revision 5930)
+++ 	(revision )
@@ -1,372 +1,0 @@
-  !!----------------------------------------------------------------------
-  !!                        ***  diawri_dimg.h90  ***
-  !!----------------------------------------------------------------------
-  !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-  !! $Id $
-  !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-  !!----------------------------------------------------------------------
-
-  SUBROUTINE dia_wri( kt )
-    !!----------------------------------------------------------------------
-    !!           *** routine dia_wri ***
-    !!
-    !! ** Purpose : output dynamics and tracer fields on direct access file
-    !!              suitable for MPP computing
-    !!
-    !! ** define key : 'key_dimgout'
-    !!
-    !! **  Method : Default is to cumulate the values over the interval between
-    !!      2 output, and each nwrite time-steps  the mean value is  computed 
-    !!      and written to the direct access file.
-    !!     If 'key_diainstant' is defined, no mean values are computed and the 
-    !!     instantaneous fields are dump.
-    !!       Each processor creates its own file with its local data
-    !!     Merging all the files is performed off line by a dedicated program
-    !!
-    !! ** Arguments :
-    !!     kt      : time-step number
-    !!     kindinc :  error condition indicator : >=0 :  OK, < 0 : error.
-    !!
-    !! ** Naming convention for files
-    !!
-    !! {cexper}_{var}_y----m--d--.dimg
-    !!   cexper is the name of the experience, given in the namelist
-    !!   var can be either U, V, T, S, KZ, SSH, ...
-    !!   var can also be 2D, which means that each level of the file is a 2D field as described below
-    !!    y----m--d--  is the date at the time of the dump
-    !!    For mpp output, each processor dumps its own memory, on appropriate record range 
-    !!    (direct access : for level jk of a klev field on proc narea irec = 1+ klev*(narea -1) + jk )
-    !!    To be tested with a lot of procs !!!!
-    !!
-    !!  level 1:  utau(:,:) * umask(:,:,1) zonal stress in N.m-2
-    !!  level 2:  vtau(:,:) * vmask(:,:,1) meridional stress in N. m-2
-    !!  level 3:  qsr + qns                total heat flux (W/m2)
-    !!  level 4:  ( emp (:,:)-rnf(:,:) )   E-P flux (mm/day)
-    !!  level 5:  tb  (:,:,1)-sst          model SST -forcing sst (degree C) ! deprecated
-    !!  level 6:  bsfb(:,:)		   streamfunction (m**3/s)
-    !!  level 7:  qsr (:,:)		   solar flux (W/m2)
-    !!  level 8:  qrp (:,:)                relax component of T flux.
-    !!  level 9:  erp (:,:)                relax component of S flux
-    !!  level 10: hmld(:,:)                turbocline depth
-    !!  level 11: hmlp(:,:)                mixed layer depth
-    !!  level 12: fr_i(:,:)                ice fraction (between 0 and 1)
-    !!  level 13: sst(:,:)                 the observed SST we relax to. ! deprecated
-    !!  level 14: qct(:,:)                 equivalent flux due to treshold SST
-    !!  level 15: fbt(:,:)                 feedback term .
-    !!  level 16: ( emp * sss )            concentration/dilution term on salinity
-    !!  level 17: ( emp * sst )            concentration/dilution term on temperature
-    !!  level 17: fsalt(:,:)               Ice=>ocean net freshwater
-    !!  level 18: gps(:,:)                 the surface pressure (m).
-    !!  level 19: spgu(:,:)                the surface pressure gradient in X direction.
-    !!  level 20: spgv(:,:)                the surface pressure gradient in Y direction.
-    !! 
-    !! History:  OPA  ! 1997-02 ( Clipper Group ) dimg files
-    !!            -   ! 2003-12 ( J.M. Molines) f90, mpp output for OPA9.0
-    !!   NEMO    1.0  ! 2005-05  (S. Theetten) add emps fsalt move gps spgu spgv 2 lines below
-    !!            -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
-    !!----------------------------------------------------------------------
-    USE lib_mpp
-    !!
-    INTEGER ,INTENT(in) :: kt
-    !!
-#if defined key_diainstant
-    LOGICAL, PARAMETER :: ll_dia_inst=.TRUE.  !: for instantaneous output
-#else
-    LOGICAL, PARAMETER :: ll_dia_inst=.FALSE. !: for average output
-#endif
-    INTEGER              , SAVE                    ::  nmoyct 
-    REAL(wp), ALLOCATABLE, SAVE, DIMENSION (:,:,:) ::  um , vm, wm   ! mean u, v, w fields
-    REAL(wp), ALLOCATABLE, SAVE, DIMENSION (:,:,:) ::  avtm          ! mean kz fields
-    REAL(wp), ALLOCATABLE, SAVE, DIMENSION (:,:,:) ::  tm , sm       ! mean t, s fields
-    REAL(wp), ALLOCATABLE, SAVE, DIMENSION (:,:,:) ::  fsel          ! mean 2d fields
-    
-    INTEGER :: inbsel, jk
-    INTEGER :: iyear,imon,iday
-    INTEGER :: ialloc
-    REAL(wp) :: zdtj
-    CHARACTER(LEN=80) :: clname
-    CHARACTER(LEN=80) :: cltext
-    CHARACTER(LEN=80) :: clmode
-    CHARACTER(LEN= 4) :: clver
-    !!----------------------------------------------------------------------
-    IF( nn_timing == 1 )   CALL timing_start('dia_wri')
-    !
-    !  Initialization
-    !  ---------------
-    !
-    IF( .NOT. ALLOCATED(um) )THEN
-       ALLOCATE(um(jpi,jpj,jpk), vm(jpi,jpj,jpk), &
-                wm(jpi,jpj,jpk),                  &
-                avtm(jpi,jpj,jpk),                &
-                tm(jpi,jpj,jpk), sm(jpi,jpj,jpk), &
-                fsel(jpi,jpj,jpk),                &
-                STAT=ialloc )
-       !
-       IF( lk_mpp      )   CALL mpp_sum ( ialloc  )
-       IF( ialloc /= 0 )   CALL ctl_warn('dia_wri( diawri_dimg.h90) : failed to allocate arrays')
-    ENDIF
-
-
-    inbsel = 18
-
-    IF( inbsel >  jpk ) THEN
-       IF(lwp) WRITE(numout,*)  ' STOP inbsel =',inbsel,' is larger than jpk=',jpk
-       STOP
-    ENDIF
-
-    iyear = ndastp/10000
-    imon = (ndastp-iyear*10000)/100
-    iday = ndastp - imon*100 - iyear*10000
-
-    !     
-    !! dimg format V1.0 should start with the 4 char. string '@!01'
-    !!
-    clver='@!01'
-    !
-    IF( .NOT. ll_dia_inst ) THEN
-       !
-       !! * Mean output section
-       !! ----------------------
-       !
-       IF( kt == nit000 .AND. lwp ) WRITE(numout,*) &
-            'THE OUTPUT FILES CONTAINS THE AVERAGE OF EACH FIELD'
-       !
-       IF( kt == nit000  ) THEN
-          ! reset arrays for average computation
-          nmoyct = 0
-          !
-          um(:,:,:) = 0._wp
-          vm(:,:,:) = 0._wp
-          wm(:,:,:) = 0._wp
-          avtm(:,:,:) = 0._wp
-          tm(:,:,:) = 0._wp
-          sm(:,:,:) = 0._wp
-          fsel(:,:,:) = 0._wp
-          !
-       ENDIF
-
-       !  cumulate values
-       !  ---------------
-
-       nmoyct = nmoyct+1
-       !
-       um(:,:,:)=um(:,:,:) + un (:,:,:)
-       vm(:,:,:)=vm(:,:,:) + vn (:,:,:)
-       wm(:,:,:)=wm(:,:,:) + wn (:,:,:)
-       avtm(:,:,:)=avtm(:,:,:) + avt (:,:,:)
-       tm(:,:,:)=tm(:,:,:) + tsn(:,:,:,jp_tem)
-       sm(:,:,:)=sm(:,:,:) + tsn(:,:,:,jp_sal)
-       !
-       fsel(:,:,1 ) = fsel(:,:,1 ) + utau(:,:) * umask(:,:,1)
-       fsel(:,:,2 ) = fsel(:,:,2 ) + vtau(:,:) * vmask(:,:,1)
-       fsel(:,:,3 ) = fsel(:,:,3 ) + qsr (:,:) + qns  (:,:) 
-       fsel(:,:,4 ) = fsel(:,:,4 ) + ( emp(:,:)-rnf(:,:) ) 
-       !        fsel(:,:,5 ) = fsel(:,:,5 ) + tsb(:,:,1,jp_tem)  !RB not used
-       fsel(:,:,6 ) = fsel(:,:,6 ) + sshn(:,:) 
-       fsel(:,:,7 ) = fsel(:,:,7 ) + qsr(:,:)
-       IF( ln_ssr ) THEN
-          IF( nn_sstr /= 0 )   fsel(:,:,8 ) = fsel(:,:,8 ) + qrp (:,:)
-          IF( nn_sssr /= 0 )   fsel(:,:,9 ) = fsel(:,:,9 ) + erp (:,:)
-       ENDIF
-       fsel(:,:,10) = fsel(:,:,10) + hmld(:,:)
-       fsel(:,:,11) = fsel(:,:,11) + hmlp(:,:)
-       fsel(:,:,12) = fsel(:,:,12) + fr_i(:,:)
-       !        fsel(:,:,13) = fsel(:,:,13)   !RB not used
-       !        fsel(:,:,14) = fsel(:,:,14) + qct(:,:)
-       !        fsel(:,:,15) = fsel(:,:,15) + fbt(:,:)
-       fsel(:,:,16) = fsel(:,:,16) + ( emp(:,:)*tsn(:,:,1,jp_sal) ) 
-       fsel(:,:,17) = fsel(:,:,17) + ( emp(:,:)*tsn(:,:,1,jp_tem) ) 
-       !
-       ! Output of dynamics and tracer fields and selected fields
-       ! --------------------------------------------------------
-       !
-       !
-       zdtj=rdt/86400.   ! time step in days
-       WRITE(clmode,'(f5.1,a)' ) nwrite*zdtj,' days average'
-
-       !       iwrite=NINT(adatrj/rwrite)
-       !      IF (abs(adatrj-iwrite*rwrite) < zdtj/2.      &
-
-       IF(  ( MOD (kt-nit000+1,nwrite) ==  0 )          &
-            &   .OR. ( kt == 1 .AND. ninist ==1 )  ) THEN
-          ! it is time to make a dump on file 
-          ! compute average
-          um(:,:,:) = um(:,:,:) / nmoyct
-          vm(:,:,:) = vm(:,:,:) / nmoyct
-          wm(:,:,:) = wm(:,:,:) / nmoyct
-          avtm(:,:,:) = avtm(:,:,:) / nmoyct
-          tm(:,:,:) = tm(:,:,:) / nmoyct
-          sm(:,:,:) = sm(:,:,:) / nmoyct
-          !
-          fsel(:,:,:) = fsel(:,:,:) / nmoyct
-          !
-          ! note : the surface pressure is not averaged, but rather 
-          ! computed from the averaged gradients.
-          !
-          ! mask mean field with tmask except utau vtau (1,2)
-          DO jk=3,inbsel
-            fsel(:,:,jk)=fsel(:,:,jk)*tmask(:,:,1)
-          END DO
-       ENDIF
-       !
-    ELSE   ! ll_dia_inst true
-       !
-       !! * Instantaneous output section
-       !! ------------------------------
-       !
-       IF( kt == nit000 .AND. lwp ) WRITE(numout,*) &
-            'THE OUTPUT FILES CONTAINS INSTANTANEOUS VALUES OF EACH FIELD'
-       !
-       zdtj=rdt/86400.   ! time step in days
-       !  iwrite=NINT(adatrj/rwrite)
-       clmode='instantaneous'
-       !     IF (abs(adatrj-iwrite*rwrite) < zdtj/2.  &
-       IF (  ( MOD (kt-nit000+1,nwrite) ==  0 )          &
-            &   .OR. ( kt == 1 .AND. ninist == 1 )  ) THEN
-          !
-          ! transfer wp arrays to sp arrays for dimg files
-          fsel(:,:,:) = 0._wp
-          !
-          fsel(:,:,1 ) = utau(:,:) * umask(:,:,1)
-          fsel(:,:,2 ) = vtau(:,:) * vmask(:,:,1)
-          fsel(:,:,3 ) = (qsr (:,:) + qns (:,:)) * tmask(:,:,1)
-          fsel(:,:,4 ) = ( emp(:,:)-rnf(:,:) ) * tmask(:,:,1) 
-          !         fsel(:,:,5 ) = (tsb(:,:,1,jp_tem) - sf_sst(1)%fnow(:,:) ) *tmask(:,:,1) !RB not used
-
-          fsel(:,:,6 ) = sshn(:,:)
-          fsel(:,:,7 ) = qsr (:,:) * tmask(:,:,1)
-          IF( ln_ssr ) THEN
-             IF( nn_sstr /= 0 )   fsel(:,:,8 ) = qrp (:,:) * tmask(:,:,1)
-             IF( nn_sssr /= 0 )   fsel(:,:,9 ) = erp (:,:) * tmask(:,:,1)
-          ENDIF
-          fsel(:,:,10) = hmld(:,:) * tmask(:,:,1)
-          fsel(:,:,11) = hmlp(:,:) * tmask(:,:,1)
-          fsel(:,:,12) = fr_i(:,:) * tmask(:,:,1)
-          !         fsel(:,:,13) = sf_sst(1)%fnow(:,:) !RB not used
-          !         fsel(:,:,14) =  qct(:,:)
-          !         fsel(:,:,15) =  fbt(:,:)
-          fsel(:,:,16) = ( emp(:,:)-tsn(:,:,1,jp_sal) ) * tmask(:,:,1) 
-          fsel(:,:,17) = ( emp(:,:)-tsn(:,:,1,jp_tem) ) * tmask(:,:,1) 
-          !
-          !         qct(:,:) = 0._wp
-       ENDIF
-    ENDIF
-    !
-    ! Opening of the datrj.out file with the absolute time in day of each dump
-    ! this file gives a record of the dump date for post processing ( ASCII file )
-    !
-    IF(  ( MOD (kt-nit000+1,nwrite) ==  0 )          &
-         &   .OR. ( kt == 1 .AND. ninist == 1 )  ) THEN
-
-       IF( lwp) WRITE(numout,*)'Days since the begining of the run :',adatrj
-
-       !! * U section
-
-       WRITE(clname,9000) TRIM(cexper),'U',iyear,imon,iday
-       cltext=TRIM(cexper)//' U(m/s) '//TRIM(clmode)
-       !
-       IF( ll_dia_inst) THEN   ;   CALL dia_wri_dimg(clname, cltext, un, jpk, 'T')
-       ELSE                    ;   CALL dia_wri_dimg(clname, cltext, um, jpk, 'T')
-       ENDIF
-
-       !! * V section
-
-       WRITE(clname,9000) TRIM(cexper),'V',iyear,imon,iday
-       cltext=TRIM(cexper)//' V(m/s) '//TRIM(clmode)
-       !
-       IF( ll_dia_inst) THEN
-          CALL dia_wri_dimg(clname, cltext, vn, jpk, 'T')
-       ELSE
-          CALL dia_wri_dimg(clname, cltext, vm, jpk, 'T')
-       ENDIF
-       !
-
-       !! * KZ section
-
-       WRITE(clname,9000) TRIM(cexper),'KZ',iyear,imon,iday
-       cltext=TRIM(cexper)//' KZ(m2/s) '//TRIM(clmode)
-
-       IF( ll_dia_inst) THEN
-          CALL dia_wri_dimg(clname, cltext, avt, jpk, 'W')
-       ELSE
-          CALL dia_wri_dimg(clname, cltext, avtm, jpk, 'W')
-       ENDIF
-       !
-
-       !! * W section
-
-       WRITE(clname,9000) TRIM(cexper),'W',iyear,imon,iday
-       cltext=TRIM(cexper)//' W(m/s) '//TRIM(clmode)
-
-       IF( ll_dia_inst) THEN
-          CALL dia_wri_dimg(clname, cltext, wn, jpk, 'W')
-       ELSE
-          CALL dia_wri_dimg(clname, cltext, wm, jpk, 'W')
-       ENDIF
-
-       !! * T section
-
-       WRITE(clname,9000) TRIM(cexper),'T',iyear,imon,iday
-       cltext=TRIM(cexper)//' T (DegC) '//TRIM(clmode)
-
-       IF( ll_dia_inst) THEN
-          CALL dia_wri_dimg(clname, cltext, tsn(:,:,:,jp_tem), jpk, 'T')
-       ELSE
-          CALL dia_wri_dimg(clname, cltext, tm               , jpk, 'T')
-       ENDIF
-       !
-
-       !! * S section
-
-       WRITE(clname,9000) TRIM(cexper),'S',iyear,imon,iday
-       cltext=TRIM(cexper)//' S (PSU) '//TRIM(clmode)
-
-       IF( ll_dia_inst) THEN
-          CALL dia_wri_dimg(clname, cltext, tsn(:,:,:,jp_sal), jpk, 'T')
-       ELSE
-          CALL dia_wri_dimg(clname, cltext, sm               , jpk, 'T')
-       ENDIF
-       !
-
-       !! * 2D section
-
-       WRITE(clname,9000) TRIM(cexper),'2D',iyear,imon,iday
-       cltext='2D fields '//TRIM(clmode)
-
-       IF( ll_dia_inst) THEN
-          CALL dia_wri_dimg(clname, cltext, fsel, inbsel, '2')
-       ELSE
-          CALL dia_wri_dimg(clname, cltext, fsel, inbsel, '2')
-       ENDIF
-
-       IF( lk_mpp )   CALL mppsync   ! synchronization in mpp
-
-       !! * Log message in numout 
-
-       IF( lwp)WRITE(numout,*) ' '
-       IF( lwp)WRITE(numout,*) ' **** WRITE in dimg file ',kt
-
-       IF( lwp .AND.        ll_dia_inst) WRITE(numout,*) '    instantaneous fields'
-       IF( lwp .AND. .NOT.  ll_dia_inst) WRITE(numout,*) '    average fields with ',nmoyct,'pdt'
-       !
-       !
-       !! * Reset cumulating arrays  and counter to 0 after writing
-       !
-       IF( .NOT. ll_dia_inst ) THEN
-          nmoyct = 0
-          !
-          um(:,:,:) = 0._wp
-          vm(:,:,:) = 0._wp
-          wm(:,:,:) = 0._wp
-          tm(:,:,:) = 0._wp
-          sm(:,:,:) = 0._wp
-          fsel(:,:,:) = 0._wp
-          avtm(:,:,:) = 0._wp
-       ENDIF
-    ENDIF
-    !
-    IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
-    !
-9000 FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc")
-    !
-  END SUBROUTINE dia_wri
